Electric dipole moment constraints on CP-violating heavy-quark Yukawas at next-to-leading order

Electric dipole moments are sensitive probes of new phases in the Higgs Yukawa couplings. We calculate the complete two-loop QCD anomalous dimension matrix for the mixing of CP-odd scalar and tensor operators and apply our results for a phenomenological study of CP violation in the bottom and charm Yukawa couplings. We find large shifts of the induced Wilson coefficients at next-to-leading-logarithmic order. Using the experimental bound on the electric dipole moment of the neutron, we update the constraints on CP-violating phases in the bottom and charm quark Yukawas.


Introduction
With the discovery of the Higgs boson in 2012, the precise determination of its couplings to all other standard model (SM) particles became a primary goal of particle physics. Of special interest are quark-Yukawa couplings. In the SM, they are real and proportional to the quark masses; any deviation from this relation would indicate physics beyond the SM. Such new contributions to Yukawa couplings are o en unavoidable in extensions of the SM that predict new particles at the LHC.
Non-SM Yukawa interactions may in fact be related to the dynamics underlying baryogenesis. For instance, new sources of CP violation are required for electroweak baryogenesis (see Ref. [1] for a review). Various models of electroweak baryogenesis require a sizeable phase in the top Yukawa. Naively, such a phase is ruled out by its contribution to the electric dipole moment (EDM) of the neutron. However, contributions of the other Yukawas can cancel the contribution to EDMs without spoiling baryogenesis [2]. is motivates a detailed study of CP-violating contributions to all Yukawa couplings.
In Ref. [3], EDM constraints on the third generation (top, bo om, tau) Yukawa couplings were presented and compared to bounds from Higgs production and decay at the LHC (see Ref. [4] for a more targeted collider study for the tau Yukawa). Presently, anomalous bo om Yukawas are more stringently constrained by h → bb than by hadronic EDMs [3]. EDM and collider constraints on the electron Yukawa were subsequently studied in Ref. [5]. A more generic approach was taken in a series of works [6][7][8], in which EDM constraints on Higgs-quark and Higgs-gluon couplings were studied in the SM e ective eld theory (EFT) approach. However, these analyses neglect logarithmically enhanced e ects that were shown to be large in Ref. [3], and two-loop Barr-Zee contributions to the light (up, down, and strange) Yukawas [9].
In the present work, we address CP violation in the bo om-and charm-quark Yukawas. e leading-logarithmic (LL) analysis for the case of the bo om quark was performed in Ref. [3]. Although not discussed there, the residual perturbative uncertainty is signi cant, exceeding the one in the non-perturbative hadronic matrix elements. Recent progress in la ice determinations of these matrix elements [10][11][12][13] motivates a next-to-leading-logarithmic (NLL) analysis in order to reduce the perturbative error of EDM predictions for the case of beyond-the-SM CP violation in the bo om and/or charm Yukawa. We calculate the complete anomalous dimension matrix for the mixing of CP-odd scalar and tensor operators up to next-to-leading (two-loop) order, and apply our results for a phenomenological study of CP violation in the bo om and charm Yukawa couplings.
is paper is organized as follows. In Section 2 we de ne the e ective theory needed for our calculation, and present the calculation of the renormalisation-group (RG) evolution of the Wilson coe cients in Section 3. We illustrate the impact of our calculation on the constraints on the CP phases in Section 4, and conclude in Section 5. e appendices contain further details on our work. In Section A we collect the requisite unphysical operators, in Section B we present all relevant renormalisation constants, and in Section C we present the full anomalous-dimension matrix. We discuss e ect on our results of a change of the renormalisation scheme in Section D. Section E contains the expansion of the resummed results.

E ective theory below the weak scale
Our starting point is a Lagrangian in which the Higgs particle couples to bo om or charm quarks di erently than in the SM. Such a modi cation can originate from TeV-scale new physics that can be parameterised by higher-dimensional operators, e.g., dimension-six operators of the form Here, H denotes the Higgs doublet in the unbroken phase of electroweak gauge symmetry, while Q L, 3 and d R, 3 represent the le -handed quark doublet and the right-handed down-quark eld of the third generation, respectively. e presence of such operators induces anomalous couplings of the Higgs particle to quarks in the electroweak broken phase. At leading order in the electroweak interactions the above dimension-six operators lead to a modi cation of the SM Yukawa interactions which, for the bo om quark, can be conveniently parameterised as (2) should be thought of as the dimension-four part of the so-called Higgs E ective Field eory [14] or the electroweak chiral Lagrangian [15] in unitarity gauge for the electroweak sector. In this sense, the SM extended by a pseudoscalar Higgs coupling is not a renormalizable theory. However, as we consider in this work only loops induced by the strong coupling, we do not encounter additional divergences since the operators of the form (1) mix only into themselves under QCD [16]. is picture would change if one considered higher-order electroweak corrections. e basic idea underlying this work is to calculate the e ect of CP-odd phases in the Higgs Yukawa couplings of the bo om and charm quark on hadronic EDMs. EDMs receive contributions from partonic CP-violating electric and chromoelectric dipole operators, with coe cients d q andd q . ey are traditionally de ned via the e ective Lagrangian valid at hadronic energies µ 2 GeV [17], with σ µν = i 2 [γ µ , γ ν ], as well as T a the fundamental generators of SU(n c ) with Tr[T a , T b ] = δ ab /2 and n c = 3 the number of colors. is Lagrangian also includes the purely gluonic Weinberg operator [18]. Its contributions are subdominant because of its small nuclear matrix elements [17,19], but are kept for completeness.
CP-violating Yukawa couplings contribute to d q andd q via Barr-Zee-type diagrams [20] with heavyquark loops, see Fig. 1. e Weinberg operator is induced via a nite threshold correction [21,22]. In the case of a CP-violating bo om Yukawa, expanding the loop functions for small bo om mass and matching directly onto the Lagrangian of Eq. (3), we nd up to higher orders in x b ≡ m 2 b /M 2 h . However, as already noted in Ref. [3], such a naive evaluation of the gluonic diagram leads to an uncertainty of a factor of order ve. e uncertainty is related to the ambiguity in choosing the proper value of the strong coupling α s (µ); namely, at which dynamical scale should it be evaluated -the weak scale, the bo om-quark mass, or the hadronic scale?
is scale dependence is related to logarithms of the large scale ratios and can be reduced by resummation of the large logarithms, which is easiest performed in an e ective theory (EFT) framework. e LL series then reproduces the quadratic logarithm in Eq. (4), while also resumming all higher-order terms. e uncertainties a er the LL resummation are still large, at the order of two at the Wilson coe cient level. Hence, in this work we extend the LL analysis of Ref. [3] to NLL and discuss the remaining theory uncertainty via the residual scale and scheme dependence in detail. In addition, we consider also modi cations of the charm-quark Yukawa.
To construct the e ective Lagrangian originating from anomalous avor-conserving, CP-violating Higgs couplings to quarks we integrate out the heavy degrees of freedom of the SM (the top quark, the weak gauge bosons, and the Higgs) and match onto an e ective ve-avor Lagrangian. EDMs are then induced by non-renormalizable operators that are CP odd. e corresponding Lagrangian reads where the sums run over all quarks with masses below the weak scale (q, q = u, d, s, c, b). e linearly independent operators are and e basis of all CP-odd operators is closed under the renormalisation group ow of both QCD and QED as they both conserve CP. Figure 2: Sample tree-level Feynman diagrams contributing to the calculation of the initial conditions of the RG evolution at the electroweak scale for the case of modi ed bo om-quark Yukawas (indicated by the red square). Light quarks are denoted by the label q = u, d, s.
In these equations, the γ 5 matrix is de ned by where µνρσ is the totally antisymmetric Levi-Civita tensor in four space-time dimensions with 0123 = − 0123 = 1, and we use the notation G a,µν = 1 2 µνρσ G a ρσ . We have included the factor 1/2 in the contributions of the O qq 3(4) operators to the e ective Lagrangian to account for the double counting implied by the sums in Eq. (5). e non-standard sign convention for O q 3 is related to our de nition of the covariant derivative, Eq. (59).
Note that whenever possible we de ned the operators directly in terms of the Levi-Civita tensors instead of γ 5 matrices.
is reduces the number of γ matrices entering the computation of the anomalous dimensions. In d = 4 dimensions our de nition is equivalent to the usual one because of the relation valid in d = 4. While such evanescent di erences in the de nition of the operators do not a ect the leading-order anomalous dimensions, they do a ect the next-to-leading-order results that we compute in this work. For further details regarding our treatment of γ 5 see section 3.2.

Renormalisation group evolution
Our goal is the summation of all leading and next-to-leading logarithms via RG improved perturbation theory in the ve-, four-, and three-avor e ective theory. e calculation proceeds in several steps. First, we integrate out the Higgs and weak gauge bosons together with the top quark at the electroweak scale, µ ∼ M h = 125.10 GeV. A er the RG ow, the heavy (bo om and charm) quarks are integrated out at their respective masses, m b (m b ) = 4.18 GeV and m c (m c ) = 1.27 GeV (all numerical input values are taken from Ref. [23]). We then match to an e ective three-avor theory where only the light-quark operators are present. Finally, we evaluate the Wilson coe cients in the three-avor theory at µ str = 2 GeV where the hadronic matrix elements of the electric dipole operators are given by the la ice calculations. e RG evolution between the di erent matching scales is computed using the appropriate anomalous dimensions, following the general formalism outlined in Ref. [24]. e actual calculation was performed with self-wri en FORM [25] routines, implementing the two-loop recursion presented in Refs. [26,27]. e amplitudes were generated using QGRAF [28]. Figure 3: Sample one-loop Feynman diagrams contributing to the calculation of the initial conditions of the RG evolution at the electroweak scale for the case of modi ed bo om-quark Yukawas (indicated by the red square). Light quarks are denoted by the label q = u, d, s.

Initial conditions at the weak scale
We augment the SM with avor-conserving anomalous Higgs Yukawas parameterised as in Eq. (2) and at a scale µ ew ≈ M h integrate out the heavy degrees of freedom of the SM. Up to quadratic order in the strong coupling constant, we nd the following non-zero initial conditions for the Wilson coe cients at scale µ ew : C q 1 (µ ew ) = − 1 + α s 4π C q 2 (µ ew ) = α s 4π We see that if a single CP-violating Yukawa is switched on, e.g., the q-Yukawa with q = b, c, then at tree-level only the operators O q 1 and O q q 1 are induced by the anomalous Higgs coupling to q, i.e., to bo om or charm quarks, see Fig. 2. At one loop, also the operators O q 2... 4 and O qq 4 receive non-zero initial conditions (Fig. 3). e quark masses in the expressions above are understood to be evaluated at the matching scale µ = µ ew . Note that the one-loop expressions depend on the de nition of evanescent operators; our choice is given in App. A. Our predictions of physical observables does not depend on the renormalisation scheme to O(α s ); we discuss the residual scheme dependence in Section 4. Figure 4: Sample two-loop Feynman diagrams whose divergent parts contribute to the calculation of the two-loop anomalous dimensions involving bo om quarks. e empty squares symbolize the insertion of an e ective operator. e diagrams of the kind shown in the middle panel, involving three di erent quark avors, vanish due to the odd number of Dirac matrices in the trace. Diagrams of the kind shown in the right panel contain traces of Dirac matrices with γ 5 , which are treated in the Larin scheme.

Calculation of the anomalous dimensions
e RG evolution below the weak scale in the presence of the heavy quarks can be calculated from the mixing of the operators in Eq. (5). e one-loop results are fully known (see, e.g., Ref. [29]). e two-loop mixing of the dipole operators has been presented in Refs. [30,31]. As a cross check we have reproduced all these results in the literature (see App. D for details). e two-loop mixing of the four-fermion operators among themselves and into the dipole operators is presented here for the rst time.
We calculate the anomalous dimensions by extracting the divergent parts of the insertions of all operators into appropriate Greens functions (see Fig. 4 for sample diagrams involving the bo om quark). We use dimensional regularisation, working in d = 4 − 2 space-time dimensions, and employ the method of infrared rearrangement (IRA) described in Ref. [32] to disentangle infrared from ultraviolet poles. For the whole computation we work in generic R ξg gauge for SU(3) c , which provides further checks for the correctness of the computation. e appearance of γ 5 in closed fermion loops requires special care, since the use of a naively anticommuting γ 5 , with {γ µ , γ 5 } = 0 for all µ (NDR scheme) is algebraically inconsistent if traces with γ 5 have to be evaluated [33]. Since such traces appear in our calculation we use Larin's prescription [34] throughout: e Levi-Civita tensor appearing either in the de nition of γ 5 or in the e ective operators is taken outside the momentum integrals and acts as a projector on the four-dimensional physical subspace a er the momentum integration has been performed. e momentum integration, including the Dirac algebra, is then performed in d space-time dimensions. e contraction of the Levi-Civita tensor with the results of the momentum integrals generates several evanescent operators that a ect the two-loop anomalous dimensions in the usual way; they are listed in App. A. No more than one Levi-Civita tensor appears in any diagram of the computation, so there is no need to contract multiple tensors.
e RG ow of the Wilson coe cients is governed by the RG equation where the superscript t denotes transposition and γ (0) , γ (1) are the one-and two-loop anomalousdimension matrices, respectively. For our NLL analysis we require the two-loop mixing of the fourquark operators among themselves and into the dipole operators. Below we collectively use the subscripts qq and q in the γ's to indicate subblocks of the full anomalous-dimension matrix. e ordering in each case is With the de nition of evanescent operators given in App. A, we nd the mixing among four-fermion operators with two di erent avors to be We nd for their mixing into four-quark operators with one quark avor and into dipoles For the anomalous dimensions in the q → q sector we nd e two-loop mixing among the dipole operators O q 3 and O q 4 has been calculated before; the relation between our results and the literature is discussed in App. D. All remaining results are new. e two-loop ADM that involve the Weinberg operator and do not necessarily vanish are γ W →qq , and γ (1) W →W . ey are not needed in our analysis, since the Weinberg operator contributes only via a nite threshold correction, see Section 3.3 (the two-loop and three-loop self mixings of the Weinberg operator for the case of N f = 0 have recently been published for pure Yang-Mills theory in Ref. [35]). All remaining subblocks are zero, i.e., e appearance of quark-mass ratios in these matrices is related to the explicit factors of quark masses in Eq. (9). ese mass ratios are scale-and scheme independent to the order we are working. ey could, in principle, be avoided altogether by de ning several dipole operators with the same eld content, but di erent quark mass factors.
For completeness, we collect the one-loop anomalous dimensions in App. C.

Matching at the heavy-quark thresholds
If m q m q , the dipole operators with external q-quark lines receive matching corrections at the q -quark threshold. (In practice, q = b or q = c.) We write the e ective Lagrangian of the theory in which all heavy q quarks have been integrated out as Here, q = u, d, s denotes one of the light quarks. CP-violating four-fermion operators involving only light quarks are present in principle and denoted by the ellipsis; however, their Wilson coe cients are suppressed by a power of m q /m q . Let us call the theory with and without the q quark the f -avor and (f − 1)-avor theories, respectively. At the threshold scale µ f the amplitudes in both theories must match each other; we write the equality between two general amplitudes in these theories as Here, angle brackets denote appropriate matrix elements. We expand the Wilson coe cients and matrix elements in the strong coupling of the f -avor theory: where the ellipses denote higher orders in the strong coupling constant. We expressed the higher-order matrix elements in terms of tree-level matrix elements, denoted by the superscript "(0)", via the coe cients r ij . ere are various subtleties to keep in mind when calculating the threshold corrections. e strong coupling constant itself receives a non-vanishing threshold correction at one loop, Figure 5: Sample Feynman diagrams whose nite parts contribute to the calculation of the matching corrections at the respective heavy-quark thresholds (shown here for the case of the bo om quark).
Similarly, the gluon eld renormalisation receives a non-zero threshold correction. Finally, the anomalous dimensions of the dipole operators depend explicitly on the number of active quark avors. e quark masses, on the other hand, are a ected only at the two-loop level (see, e.g., Ref. [36]). By explicit calculation of various one-loop diagrams (see Fig. 5), we nd, evaluating Eq. (26), the following non-zero threshold corrections as a power series in the strong coupling in the (f − 1)-avor theory.

Numerics
In the last three sections we presented all the ingredients needed to consistently perform the resummation of large logarithms appearing in hadronic Barr-Zee-type diagrams from avor-conserving CP-violating Higgs Yukawas at NLL: the next-to-leading-order (NLO) (one-loop) matching at the electroweak scale, the NLO (two-loop) anomalous-dimension matrix, and the (one-loop) threshold corrections for the Wilson coe cients at the heavy-quark thresholds. In this section, we implement the NLL evolution numerically and discuss its impact on a set of hadronic EDMs. We rst present values for the partonic Wilson coe cients in dependence of the CP-violating phase and discuss the theoretical uncertainty in detail. en we give bounds on the phases using experimental input.

Wilson coe icients
To compute the e ect of modi ed Yukawa couplings on hadronic EDMs, we need the values of the induced Wilson coe cients of the electric dipole, chromoelectric dipole, and Weinberg operators at the scale µ str = 2 GeV where the matrix elements of these operators are evaluated. We consider two cases: rst, we only turn on a CP-violating bo om Yukawa and second, only a CP-violating charm Yukawa. e dependence on the matching scales µ ew and µ b(c) of the dipole Wilson coe cients evaluated 80 120 160  at scale µ str = 2 GeV cancels at NLO in the expansion in α s . 1 However, the RG evolution induces a residual dependence on these scales. is dependence is formally of higher order in α s and can be used to assess the remaining theoretical uncertainty of the prediction. In Fig. 6 we show the value of the dipole Wilson coe cients, evaluated at 2 GeV, as a function of the electroweak matching scale, For the dipole operators, the la er terms are subleading. For purposes of illustration we thus choose to plot the case φ b(c) = π/2, se ing κ b(c) = 1 and factoring out the global Higgs-mass dependence. Focusing rst on the Wilson coe cients C q 3 4 5 6 7 8 9 10 Bo om: u-quark dipoles reduced by going from LL to NLL. For the coe cients C q 4 of the chromoelectric dipole operator (blue lines) the scale dependence δC q 4 /C q 4 is reduced by factor of approximately ve. In all cases, however, the shi from the LL to the NLL results is substantial and larger than indicated by the residual scale dependence of both the LL and NLL result.
It is interesting to note that the RG evolution, in the case of a modi ed charm Yukawa, leads to a small Wilson coe cient of the electric dipole operator, such that the contribution to the hadronic EDMs at low scales is completely dominated by the chromoelectric dipole operator.
Next, we consider the residual dependence on the matching scale at which the heavy quark is integrated out. Similarly to Fig. 6, in Fig. 7 we show the Wilson coe cients at 2 GeV, this time varying the heavy-quark matching scales within 2.5 GeV ≤ µ b ≤ 10 GeV and 1 GeV ≤ µ c ≤ 2.5 GeV. While the scale dependence of the electric dipole coe cients C q 3 is still mild (red lines), the chromoelectric coe cients C q 4 (blue lines) show a large residual scale dependence that is barely reduced in going from LL to NLL. Again, all Wilson coe cients show a large shi when including the NLL corrections.
While the large scale dependence at NLL can be partly understood by the appearance of many new,
One may then wonder whether these large entries in the anomalous-dimension matrix are an artefact of a badly chosen renormalisation scheme. In fact, as discussed for instance in Refs. [37,38], these entries depend on the de nition of evanescent operators. Needless to say, this scheme dependence cancels up to the order to which the calculation is performed. Nevertheless, the residual scheme dependence can, in principle, be large. We have tested this by converting our anomalous dimension to di erent schemes chosen such that many of the large entries vanish. While the residual scale dependence is indeed somewhat smaller in these schemes, the central values of the Wilson coe cients strongly depend on the choice of scheme, taking values in approximately the same range as indicated by the scale dependence of the results in our original calculation.
All these observations hint at a slow convergence of the perturbation series. We therefore adopt the following prescription for the numerical values of the Wilson coe cients, including our best estimate of the associated remaining theory uncertainty: we obtain the "central values" of the Wilson coe cients as the value for µ ew = M h and µ b(c) = m b(c) (m b(c) ) and assign as the theory uncertainty either half the range of the NLL scale variations, or half the shi between LL and NLL, whichever is larger.
We then nd the three-avour Wilson coe cients, evaluated at 2 GeV, of electric dipole, chromoelectric dipole, and Weinberg operators to be: Here, the subscripts q = b, c refer to the cases of the bo om and charm quark, respectively. e values for the a and b coe cients and their respective uncertainties are given, for the two cases, in the second and third columns of table 1. We see that in both cases the residual uncertainty on the dipole Wilson coe cients is of the order of 30%, while the contributions induced by the Weinberg operator have larger uncertainties. ey could potentially be reduced by a next-to-next-to-leading-logarithmic calculation.  In this section we derive constraints on new CP-violating phases in the bo om or charm Yukawa under the simplifying assumption that only one such phase is present. A global t for the general case including constraints from the LHC will be presented in a future publication [39].

Bounds on CP phases from hadronic EDMs
We start with a generic parameterisation of a nuclear dipole moment, X, as Here, the f X,i with i = 1, . . . , n are the hadronic input parameters entering the prediction of the given d X . We denote their uncertainties by ∆f X,i . e a i , with i = 3u, 3d, 4, w, and b i , with i = 3, 4, w, are de ned in Section 4.1 and given in table 1. ey parameterise the Wilson coe cients that contribute to d X if either the bo om or the charm quark Yukawa is modi ed. By standard quadratic error propagation we compute the total theory uncertainty To derive the allowed con dence level (CL) intervals from the measurements of dipole moments and to combine them we construct a χ 2 function of two parameters, κ q , φ q or equivalently of κ q sin φ q , κ q cos φ q : where we have neglected the tiny SM contribution to any EDM. e allowed 68.27% CL region for the two-parameters are then given by the region χ 2 (κ q , φ q ) − χ 2 min ≤ 2.30. e relation between the coe cients d q ,d q , and w in Eq. (3) and the Wilson coe cients in the three-avour EFT is In the expression for d q above we have included the electroweak contribution of Eq. (4). From this we obtain the neutron EDM as where the matrix elements of the electric dipole operator are parameterised by g u T = −0.24(3), g d T = 0.85(8), g s T = −0.012 (18). ese values are calculated using la ice QCD with 2 + 1 active avors and are converted to the MS scheme at 2 GeV [40,41] (see Ref. [42] for a N f = 2 + 1 + 1 result). By using the la ice values of the three-avour theory we include the e ect of threshold corrections associated to the charm quark. To capture this e ect in the four-avour theory the matrix elements of four-quark operators with charm quarks would also need to be evaluated via la ice methods. e hadronic matrix elements of the chromoelectric dipole and the Weinberg operators are estimated using QCD sum rules and chiral techniques [17,19]. Notice that the sign of the hadronic matrix element of the Weinberg is not known, and thus the allowed CL intervals will depend on it. For prospects on obtaining the la er via la ice calculations see Refs. [12,13]. e experimental 90% CL exclusion bound obtained in Ref. [43] is |d n | < 2.9×10 −26 e cm. Using the central values of the Wilson coe cients in table 1 and the two-parameter χ 2 we compute the allowed 68.27% CL intervals for the bo om-and charm-quark cases. e label "sign w " indicates whether the sign of the Weinberg-operator contribution in Eq. (38) is taken to be positive or negative. We show the CL intervals for the cases in which: i) no theory uncertainty is included (no theory error label), ii) only the short-distance theory uncertainty is included (with short-distance theory error label) iii) the short-distance theory uncertainty is added in quadrature with the present theory uncertainties of the hadronic input (with theory error label). For brevity we introduce the short-hand notation For the bo om case we nd the allowed 68.27% CL regions to be: For the charm case we nd: Due to the large theory uncertainties there is no 68.27% CL interval for the case sign w = − when the full theory uncertainties are included.
Other hadronic EDMs give, in principle, complementary bounds. For instance, the contribution to the mercury EDM is given by [19] d Hg Using the current upper experimental 95% CL bound |d Hg | < 3.1 × 10 −29 e cm [44] we compute the allowed 68.27% CL intervals from the two-parameter χ 2 . e presently large hadronic uncertainty in Eq. (50) does not constrain the parameter space at the 68.27% CL. We thus include only the theory uncertainties associated to short-distance dynamics in our bounds. For the bo om case we nd the allowed 68.27% CL regions to be: [with short-distance theory error] . (52) For the charm case we nd We see that even if we neglect the present theory uncertainties the constraints from mercury EDM cannot compete with the ones from the neutron EDM. It is instructive to compare the constraints obtained from hadronic EDMs to the constraints from the bound on the electron EDM, obtained recently by the ACME collaboration. e contribution of a modi ed bo om Yukawa to the electron EDM can be easily obtained by the substitutions Q q → Q e and m q → m e in Eq. (4), and similarly for a modi ed charm Yukawa. Using the ACME result, |d e | < 1.1×10 −29 e cm (90% CL) [45], we compute the corresponding allowed intervals. e electron EDM depends solely on the combination κ q sin φ q , which is thus constrained by a one-parameter χ 2 function to be κ b |s b | ≤ 0.32 and κ c |s c | ≤ 0.82 at the 68.26% CL for the bo om-and charm-quark case, respectively. To compare with the neutron and mercury EDM allowed two-parameter 68.26% CL intervals we also list the corresponding ones for the electron EDM: Since for the electron EDM there is no hadronic input, the theory uncertainties originate solely from higher electromagnetic corrections and are small. 2 We see that currently the bound from the electron EDM is stronger than the one from the neutron EDM. However, both experimental progress and the anticipated la ice calculations will strengthen the bounds from neutron and other hadronic EDMs. e combination of leptonic and hadronic EDMs can also be used as a strategy to disentangle e ect of having multiple CP-violating Yukawas.
We illustrate the results of this section for the bo om-and charm-quark cases in Figs. 8 and 9, respectively. We show in colour the allowed 68.26% CL regions of the two-parameter space for di erent EDMs. In the plots on the le we take the two parameters to be κ b(c) sin φ b(c) and κ b(c) cos φ b(c) ; in the plots on the right we choose the parameters to be κ b(c) and φ b(c) . In the upper plots we have included no theory uncertainties; in the lower ones we folded in the theory uncertainties associated to short-distance dynamics.

Conclusions
We presented the complete two-loop QCD anomalous-dimension matrix for the mixing of CP-odd scalar and tensor operators in an EFT valid below the electroweak-symmetry breaking scale. We used the results to perform a next-to-leading-logarithmic RG analysis of the Wilson coe cients from the weak scale to the hadronic scale of 2 GeV, calculating also the requisite nite matching corrections at the heavy-avor thresholds. We applied our calculation to a new-physics scenario where new, avor-conserving, CP-violating phases appear in the Higgs Yukawa couplings to the bo om or charm quark. We calculated the initial conditions at the weak scale up to NLO, and solved the RG equations to compute the induced coe cients of the CP-violating electric and chromoelectric dipole operators and the Weinberg operator in the three-avor EFT at the hadronic scale 2 GeV.
We nd large shi s, as well as a large residual scale and scheme dependence of the dipole Wilson coe cients at next-to-leading-logarithmic order. We interpret this as a hint of a slowly converging perturbation series. e dipole and Weinberg operators contribute to the electric dipole moment of the neutron and mercury. Assuming a Peccei-inn-type solution to the strong CP problem we can then derive constraints on the modi ed Yukawa couplings from the experimental bound on the neutron EDM.
ese constraints are currently not as strong as those derived from measurements of the electron EDM, but will play an important role in future global ts to modi ed Higgs Yukawas. is is true in particular in light of the progress expected in la ice calculations of the hadronic matrix elements, and future improvements regarding the experimental bounds on various hadronic EDMs.

A Unphysical operators
Additionally, the following operators, which are not gauge-invariant, are also required in intermediate steps to determine all counterterms and project the o -shell amplitudes e covariant derivative acting on quarks is de ned as with Q q the quark electrical charge. Accordingly, the gluon eld-strength tensor is e covariant derivative acting on color octets is given by

A.2 Evanescent operators
Next we list the evanescent operators that enter our computation at one and two-loop order. e leading-order anomalous dimension does not depend on their de nition, but the next-to-leading-order one does.
In the q-q-A sector we require the operator Analogously, in the q-q-G sector we require the operator Notice that in Eqs. (62) and (63), as well as in those that follow, we have de ned the γ-algebra structure in terms of anticommutators with γ 5 , i.e., {Γ, iγ 5 }.
is ensures that the operators are self-conjugate not just in d = 4 dimensions but also in d = 4 − 2 dimensions. e evanescent operators required in the q-q sector read e square brackets denote antisymmetrisation normalized as In the q-q sector they read In simplifying the color algebra, we use the following standard relation for the generators of SU(n c ): Consequently, Fierz relations on the Lorentz structures, valid in d = 4, need to be applied on the operators with T a 's to show that they are evanescent.

A.3 Operators related to the infrared rearrangement
e last class of unphysical operators arises because our infrared rearrangement breaks gauge invariance in intermediate steps of the calculation. At the renormalizable level this method generates one gauge-variant operator corresponding to a gluon-mass term, i.e., e "gluon mass", m IRA , is completely arti cial and drops out of all physical results, and Z IRA is an additional renormalisation constant [32]. At the non-renormalizable level the one-loop insertions of the dimension-ve and dimension-six operators can also induce gauge-invariant operators that are relics of the infrared rearrangement. For our calculation, the only relevant ones are the following three

B Renormalisation constants
To obtain the two-loop anomalous dimension of the physical sector we need certain one-loop renormalisation constants involving unphysical operators. We collect them in this appendix. We use the following standard notation for their expansion in α s and In MS, the mixing of evanescent operators into physical also includes nite terms [47], thus in this case the expansion starts with l = 0. e subscripts x and y symbolize sets of Wilson coe cients, for which we use the following notation and standard ordering: e rst necessary input is the mixing of the physical operators into all the evanescent operators that are generated at one-loop. Using the same subscript notation as above, the renormalisation constants (72) e remaining nite mixings of evanescent into physical operators are zero. Furthermore, we need the mixing constants of the physical operators into the operators arising from infrared rearrangement; they are found to be All other mixing constants of physical into the IRA operators are zero. Special care must be taken to obtain the mixings like Z (1,1) qq →P q . Apart from the obvious q ↔ q interchange also the ordering of the operators in the collective blocks is relevant. For example Finally, we need the mixing constants of the physical operators into the e.o.m.-vanishing operators. ey are uniquely xed by the q → q, q → qγ, and q → qg Greens functions and we have veri ed that their values are consistent with the renormalisation of the qg → qγ and qg → qg Greens functions. We nd e two-loop anomalous dimension matrix is given in terms of the one-and two-loop renormalisation constants by e quadratic poles of the two-loop diagrams are xed by the poles of the one-loop diagrams via where β 0 = 11 3 n c − 2 3 N f . As a check of our calculation, we computed these poles directly and veri ed that they satisfy Eq. (76).
In our calculation we needed various eld and mass renormalisation constants up to two-loop level, and we have calculated them explicitly. Writing the expansion with r = q, m, g s , G, u, m IRA denoting the quark eld, quark mass, strong coupling, gluon eld, ghost eld, and arti cial gluon mass renormalisation, respectively, we nd e explicit results read All other physical subblocks are zero at one-loop. Our one-loop results for the physical sector agree with the results in the literature [29,49,50], a er accounting for the di erent normalisation of the operators and the di erent conventions in the covariant derivative.
D Change of renormalisation scheme e (three-loop) mixing among the dipole operators O q 3 and O q 4 has already been calculated in the literature [31] (see Ref. [30] for the earlier two-loop result for the CP-even dipole operators). In this section we show that our results are consistent with these calculations.
ere are three main di erences between our calculation and the earlier ones: (i) the de nition of the physical operators di ers in d space-time dimensions, (ii) the calculation in Ref. [31] has e ectively been performed in the NDR scheme whereas ours in the "Larin" scheme, and (iii) a di erent normalisation of the operators with factors of quark masses and the strong coupling has been chosen. Here we show that there is a unique change of the renormalisation scheme that transforms our result into that of Refs. [30,31].
To address (i), we perform a rede nition of our basis of physical operators as follows: 3 3 e la ice results used in our numerics have been converted to the NDR-MS scheme at the one-loop level [40]. If, in the future, the two-loop conversion of the la ice results becomes available, one should perform the change of scheme discussed in this appendix for our analysis.
All other physical operators are unchanged, i.e. Q i ≡ Q i . e resulting operators correspond, in the NDR scheme, to the operators used in Refs. [30,31], up to a normalisation (see below).
To address (ii), we mimick the results obtained in the NDR scheme by rede ning the evanescent operator E q g to be: with all other operators unchanged. e coe cients in front of the physical operators are uniquely determined by the requirement to reproduce the anomalous dimensions of the dipole operators in the literature. Finally, to address (iii) we take care of the di erent overall normalisation of the operators as follows. If we write the shi ed renormalisation constants as Z = ρZ, where then the shi ed ADMs are given by γ (0) = γ (0) + 2ρ (1,1) , γ (1) = γ (1) + 4ρ (2,1) .
Combining this shi with the change of physical and evanescent operators speci ed above (see, e.g., Ref. [51] for the general expressions), we can reproduce both the one-and two-loop mixing of the dipole operators in Ref. [30,31]. e mixing of four-fermion into dipole operators, on the other hand, cannot be calculated in the NDR scheme. For completeness, we provide the full physical two-loop anomalous dimension matrix in the primed basis:

E Expanding the renormalisation group
To gain a be er understanding of our results, and as an additional check of our calculation, we expand the full solution of the RG equations about the bo om-quark threshold (the procedure for the charm quark is analogous).
Keeping only the leading nonvanishing terms and including the QED contribution, Eq. (4), we nd for the electric dipole where the superscripts "(0)" and "(1)" denote the tree-level and one-loop contributions to the initial conditions of the Wilson coe cients at the weak scale (we omit here the small logarithmic contributions ∝ log(µ 2 ew /M 2 h )). Using α s (m b ) ∼ 0.22, the ratio of QED to LL to NLL is roughly 1 : −9 : −4. We see (as observed already in Ref. [3]) that the contribution of the photonic Barr-Zee diagram is negligible. We also see that the NLL correction is quite large, about half the size of the LL contribution. e leading terms in the expansion of the solution for the RG equations for the chromoelectric dipole, on the other hand, should exactly reproduce the logarithmic parts of the result in Eq. (4) (the constant term is of next-to-next-to-leading logarithmic order and can only be reproduced, within EFT, by a three-loop calculation). Indeed, we nd 1,qq →4,qq γ (0) 4,qq →4,q As expected, the leading contribution to the LL reproduces the quadratic logarithm in Eq. (4), while leading contribution to the NLL result vanishes. is means, in turn, that the NLL corrections tod q start at relative order α s /(4π), with large anomalous-dimension prefactors of (γ  4,qq →4,q )/8 = 56/3, multiplied by the NLO initial condition 3/2. erefore, these constitute sizeable relative corrections, with ratio 1 : −0.6 between the LL and the NLL contributions. Needless to say that in our numerics we use the full solution of the RG equations, where all large logarithms are resummed to leading and next-to-leading order.
For completeness we give the also expansion of the Dicus result [52]: is reproduces the logarithmic term in Eq. (4).