A Precise Electron EDM Constraint on CP-odd Heavy-Quark Yukawas

CP-odd Higgs couplings to bottom and charm quarks arise in many extensions of the standard model and are of potential interest for electroweak baryogenesis. These couplings induce a contribution to the electron EDM. The experimental limit on the latter then leads to a strong bound on the CP-odd Higgs couplings. We point out that this bound receives large QCD corrections, even though it arises from a leptonic observable. We calculate the contribution of CP-odd Higgs couplings to the bottom and charm quarks in renormalisation-group improved perturbation theory at next-to-leading order in the strong interaction, thereby reducing the uncertainty to a few percent.


Introduction
The precise determination of the Yukawa couplings of the Higgs boson to all fermions has been a focus of particle physics since the discovery of the Higgs boson in 2012 [1,2].In the Standard Model (SM) all Yukawa couplings are aligned with the fermion masses and thus real, but multiple extensions of the SM induce non-trivial phases.This is of particular interest as these phases (mainly in the Yukawa couplings to the third fermion generation) play a leading role in models of electroweak baryogenesis [3][4][5].In this article, we focus on the bottom-and charm-quark Yukawa couplings.
It is well-known that the present bounds on the electric dipole moment (EDM) of the electron place strong constraints on CP-violating phases in the quark Yukawa couplings [6][7][8][9][10][11].However, it is less well-appreciated (and maybe somewhat surprising) that the heavy-quark contributions to the electron EDM receive large QCD corrections, leading to a large implicit uncertainty in the current constraints.In this section, we briefly review the current situation.In the remainder of this article, we calculate the leading logarithmic (LL) and next-to-leading logarithmic (NLL) QCD corrections, thereby reducing the presently O(1) uncertainty to a few percent.
For the purpose of this work, we assume a modification of the SM heavy-quark Yukawa couplings of the form 1   L hq h q h = − y SM q h √ 2 κ q h qh (cos ϕ q h + iγ 5 sin ϕ q h ) q h h .
Here, q h = b, c denotes the bottom-or charm-quark field and h the physical Higgs field.Moreover, y SM q h ≡ m q h e/ √ 2s w M W is the SM Yukawa, with e the positron charge, s w the sine of the weak mixing angle, and m q and M W the heavy-quark and W -boson masses, respectively.The real parameter κ q h ≥ 0 parameterises modifications to the absolute value of the Yukawa coupling, while the phase ϕ q h ∈ [0, 2π) parameterises CP violation.The SM corresponds to the values κ q h = 1 and ϕ q h = 0.
Virtual heavy quarks with CP-odd Higgs couplings induce an electron EDM d e , defined by with σ µν ≡ i 2 [γ µ , γ ν ], via the well-known Barr-Zee diagrams [14].Calculating the two-loop Barr-Zee diagrams (see Figure 1) keeping the non-zero heavy-quark mass and expanding the loop functions for small quark mass and small external momenta leads to [6] d e ≃ −12eQ e Q 2 q h α e (4π) 3  √ 2G F m e κ q h sin ϕ q h x q h log 2 x q h + π 2 3 , up to higher orders in the ratio x q h ≡ m 2 q h /M 2 h (with M h ∼ 125 GeV the Higgs-boson mass).Here, α e is the fine-structure constant, Q q h is the charge of the heavy quark with Q b = −1/3 and Q c = +2/3, and Q e = −1 is the electron charge.Barr-Zee diagrams with an internal Z boson also affect the electron EDM, however, they lead to a contribution that is suppressed with respect to those with an internal photon by the small coupling of the Z boson to electrons [6].
Using the recent bound |d e | < 4.1 × 10 −30 e cm (at 90% CL) [15], Eq. ( 3) implies the bounds 15 and κ c | sin ϕ c | ≤ 0.41 at the 90% CL for the bottom-and charm-quark case, respectively.However, to obtain these bounds, numerical values for the heavy-quark masses must be chosen.It is not clear a priori at which scale the quark masses should be evaluated; µ = M h and µ = m q h would be obvious choices.For the bounds above, we used the values m b (M h ) and m c (M h ); choosing m b (m b ) and m c (m c ) instead leads to the significantly stronger bounds κ b | sin ϕ b | ≤ 0.086 and κ c | sin ϕ c | ≤ 0.13 at the 90% CL.The differences arise from the large QCD running of the quark masses between the two scales µ = M h and µ = m q h .This indicates that the QCD corrections to Eq. ( 3) are large, even though the electron EDM is a leptonic observable.By our explicit calculation, we will show that the predicted value for d e after resolving the ambiguity lies somewhere in between the values obtained by using the two scales M h and m q h .We give here a brief outline of the main ideas of the calculation.First, notice that the result in Eq. ( 3), which was obtained by a two-loop, fixed-order calculation, is numerically dominated by the large quadratic logarithms log 2 x q h .This logarithmic contribution can be reproduced by the one-loop QED renormalisation-group (RG) evolution in an appropriate effective theory (see Section 2), truncated at order α 2 e .The second term in Eq. ( 3), π 2 /3, has no logarithmic dependence on the mass ratio and is formally of next-to-next-to-leading-logarithmic (NNLL) order in RG-improved perturbation theory.Numerically, this term gives a ∼ 6% correction to the logarithmic contribution to d e /e in the bottom case and a ∼ 3.5% correction in the charm case.On the other hand, we have seen above that the choice of different renormalisation scales for the quark masses represents a much larger uncertainty, of the order of 100%.
We, thus, conclude that the QCD corrections to the Barr-Zee diagrams are large, and we expect that these corrections are dominated by leading QCD logarithms, since the product α s log x q h is large (α s denotes the strong coupling constant).These logarithms can be reliably calculated using RG-improved perturbation theory.The result with resummed leading QCD logarithms has the schematic form Here, the first term reflects the quadratic logarithm in Eq. ( 3), while all other terms correspond to the leading logarithms of the diagrams obtained by dressing the Barr-Zee diagrams, Figure 1, with an arbitrary number of gluons.
It is, however, well-known that one can only consistently fix the scale and scheme dependence of the input parameters by going beyond the LL approximation.Hence, we will also perform the NLL calculation, reducing the the uncertainty down to the percent level.This corresponds roughly to the size of the correction that we have estimated above using the π 2 /3 term in the fixed-order result, which can be viewed as part of the NNLL order result in RG-improved perturbation theory.This term is thus not included in our NLL calculation.
Our calculation will show that the QCD perturbation series converges well, as might be expected for a leptonic observable.In addition, we emphasize that this is a complete calculation -no further hadronic input (such as the lattice matrix elements required for hadronic EDMs) is needed.The final result for d e /e lies between the two values obtained from the naive computation, determining the heavy quark masses at the different scales.All these results are illustrated in Figure 4 of Section 4.
This paper is organised as follows: in Section 2, we introduce the effective theories used in the computation.In Section 3, we present the detailed results of our calculation, namely, the initial conditions at the electroweak scale, the calculation of the anomalous dimensions, and the threshold corrections at the heavy-quark thresholds.We also show the final analytical results in a compact form.In Section 4, we present the numerical results of the calculation including updated bounds on the CP-odd heavy-quark Yukawa couplings.We conclude in Section 5. Additional information is presented in two appendices: In App.A, we collect the unphysical operators used in the computation of the anomalous-dimension matrix and in App.B, we show the results for renormalisation constants entering the calculation.

Effective Theories Below the Weak Scale
A precise determination of the electron EDM in the presence of CP-odd Higgs couplings to the bottom and charm quarks requires an effective field theory that allows to sum large QCD logarithms to all orders in the strong coupling constant and includes the effect of the combined α s and α e RG evolution.Based on the "full" Lagrangian in Eq. ( 1), we construct the effective Lagrangian below the electroweak scale, µ ew , by integrating out the heavy degrees of freedom of the SM (the top quark, the weak gauge bosons, and the Higgs).EDMs are then induced by non-renormalisable operators that are CP odd.In the current work, we focus on the part of the effective Lagrangian that is relevant for predicting the electron EDM in the presence of CP-odd Yukawas couplings to the bottom and charm quarks.In this case, the relevant effective Lagrangian reads where the four linearly independent operators are and C eq h 1 , C q h e 1 , C eq h 2 , and C e 3 are the corresponding Wilson coefficients.Additional CP-odd operators that are suppressed by additional factors of m light /m q h (where m light corresponds to a light quark mass) are denoted by the ellipsis.We defined the electron dipole operator with a factor of the running quark mass m q h ≡ m q h (µ), to avoid awkward ratios of quark and lepton masses in the anomalous dimensions.Throughout this work the γ 5 matrix is defined by where ϵ µνρσ is the totally antisymmetric Levi-Civita tensor in four space-time dimensions with ϵ 0123 = −ϵ 0123 = 1, and we use the notation F µν = 1 2 ϵ µνρσ F ρσ .We treat γ 5 within the "Larin" scheme; for the details and subtleties we refer to Ref. [9].The non-standard sign convention for O e 3 is related to our definition of the covariant derivative acting on fermion fields f , The basis of all flavour-diagonal, CP-odd operators is closed under the QED and QCD RG flow, as both interactions conserve CP and flavour.Below the electroweak scale, there is a tower of effective theories relevant for predicting the electron EDM.For the bottom case, q h = b, we employ the effective Lagrangian in Eq. ( 4) for the five-flavour theory, while for the charm case, q h = c, we employ Eq. ( 4) for both the five-flavour and the four-flavour theory (see discussion in Section 3.3 for the threshold corrections on couplings and Wilson coefficients at the bottom-quark scale).
The effective theory below the heavy-quark scale, µ q h , does not contain four-fermion operators with heavy quarks, and we use the modified effective Lagrangian in which -as opposed to Eq. (4) -the dipole operator is defined with the conventional factor m e : m e e (ēσ µν e) Fµν .
Note that for q h = b the Lagrangian in Eq. ( 8) refers to the four-flavour Lagrangian, while for q h = c to the three-flavour one.This definition of L eff,q h implies that (cf.Eq. ( 2)) In the next section, we describe how the RG evolution within the effective field theories relates d e to the parameters κ q h and ϕ q h of the "full" theory in Eq. (1).

Renormalisation Group Evolution
Our goal is the summation of all leading and next-to-leading logarithms via RG-improved perturbation theory.The calculation proceeds in the following steps.First, we integrate out the Higgs and weak gauge bosons together with the top quark at the electroweak scale, µ ew ∼ M h .This matching calculation at µ ew induces the initial conditions for the five-flavour Wilson coefficients appearing in Eq. ( 4).We collect them in Section 3.1.Subsequently, we perform the RG evolution from µ ew down to the bottom-quark threshold, µ b ∼ m b (m b ).The anomalous dimensions relevant for the mixed QCD-QED RG evolution at NLL accuracy are computed here for the first time, see Section 3.2.The next step depends on whether q h = b or q h = c.For the bottom case, q h = b, we match directly to the four-flavour version of the Lagrangian in Eq. ( 8) to obtain the prediction of the electron EDM.The relevant threshold corrections at µ b are discussed in Section 3.3.For the charm case, q h = c, we must instead match at µ b to the four-flavour version of the Lagrangian with four-fermion operators in Eq. ( 4) and additionally perform the RG evolution in the four-flavour theory from µ b down to µ c ∼ m c (m c ).The corresponding anomalous dimensions are also given in section 3.2.Finally, we match to the three-flavour version of the Lagrangian in Eq. ( 8) to obtain the prediction of the electron EDM.
The calculations of the amplitudes relevant for computing the initial conditions of Wilson coefficients and the hitherto unknown anomalous dimensions have been performed with the package MaRTIn [16] which is based on FORM [17] and implements the two-loop recursion algorithms presented in Refs.[18,19].The amplitudes were generated using QGRAF [20].
The RG evolution between the different matching scales is significantly more involved than in applications in which the electromagnetic coupling, α e , can be neglected.The reason is that the leading contribution to the electron EDM contains a term with two powers of the large logarithm, i.e, log 2 x q h , see Eq. ( 3).Within the effective theory, this term is obtained by a LL QED calculation, truncated at order α 2 e .However, the numerically relevant corrections to this result are not further electromagnetic α n e log n corrections, but large, logarithmically enhanced QCD corrections which must be summed to all orders to obtain an accurate prediction.Therefore, to properly account for the numerically relevant corrections we must consistently solve the mixed QCD-QED RG equations.
This can be achieved consistently using the general formalism developed in Ref. [21].The main idea is that, since α s log x q h is O(1), such products of QCD coupling times large logarithms must be resummed for an accurate prediction.By contrast, the product α e log x q h is small and does not require resummation.The large logarithms appearing in α e log x q h are expressed in terms of the resummed α s log x q h , i.e. α e log x q h = (α e /α s ) × α s log x q h .As a consequence, the conventional expansion in terms of α s and α e is replaced by an expansion in α s and κ ≡ α e /α s .Therefore, in this setup the Wilson coefficients at some scale µ have the expansion with αs ≡ α s /4π.By implementing the mixed RG evolution we will compute the C (nm) i (µ) coefficients relevant for the electron EDM; for details we refer to Ref. [21].The analytical results are presented in Section 3.4.

Initial Conditions at the Weak Scale
We augment the SM by flavour-conserving, CP-odd Higgs Yukawa couplings to the heavy quarks q h = b, c, as parameterised in Eq. ( 1).At a scale µ ew ≈ M h we integrate out the heavy degrees of freedom of the SM and match to the effective five-flavour theory relevant for the electron EDM in Eq. ( 4).The initial conditions for the Wilson coefficients relevant for our calculation are obtained by evaluating the tree-level and one-loop Feynman diagrams such as those shown in Figure 2; we find All parameters appearing above correspond to parameters in the five-flavour theory evaluated at the scale µ ew , i.e., m q h = m q h ,5fl (µ ew ) and α e = α e,5fl (µ ew ).We treat the heavy quarks and the electron as massless at the electroweak matching scale; therefore, no powers of m e /M h or m q h /M h appear in the results.The explicit factors of the electron and heavy-quark masses arise by expressing the Yukawa couplings in Eq. ( 1) in terms of m e and m q h .We have included all terms up to corrections of quadratic order in the strong and electromagnetic coupling constants, as only these are required for our analysis.
The O(α e ) coefficient, Eq. ( 14), is well-defined only after specifying the basis of evanescent operators which can be found in App. A. In Section 3.4, we will use the framework of Ref. [21] to solve the mixed QCD-QED RG as an expansion in αs and κ.Based on Eqs. ( 12)-( 14) we find the contributing expansion coefficients in Eq. ( 11)

Anomalous Dimensions
To solve the RG in the effective theories below the electroweak scale we need to include the running of coupling constants, mass anomalous dimensions, and the mixing of operators.In this section, we collect the results that enter the analysis at NLL accuracy.The running of the Wilson coefficients from the electroweak matching scale down to the relevant quark scale is governed by the RG equation where γ ji are the components of the anomalous dimension matrix (ADM).We choose the ordering of Wilson coefficients as The ADM admits an expansion in powers of2 αs ≡ α s /4π and αe ≡ α e /4π, with αs , αe , and γ (nm) depending on the number of active fermion flavours in the effective theory.Using this expansion, the ADM can be organised by loop order.In general, it depends on the number of active fermion flavours and on the flavour of the heavy quark q h .Below we present our results entering the five-and four-flavour RG evolution required for the case q h = b and q h = c.By explicit calculation, we find at one loop with is the number of quark colors, n u is the number of active up-type quarks, n d is the number of active down-type quarks, n ℓ is the number of active charged leptons and In both the bottom and charm cases, we have n u = 2 and n ℓ = 3.
At two-loop, the pure QCD ADM is the same for the bottom-and charm-quark case and depends on e q h e q h e q h e q h e e γ q h the number of active quark flavours, For the bottom case, we only need the five-flavour theory, i.e., we always fix n d = 3.For the charm case, we must solve the RG in both five-and four-flavour theory in which case n d = 3 and n d = 2, respectively.
For the two-loop mixed QCD-QED and the pure two-loop QED ADMs, explicit factors of the quark charges result in different results for the bottom and charm cases.For notational simplicity, we quote the ADMs for the two cases separately after having substituted the electric charges.For the bottom, we find, after fixing [ For the charm, we find, keeping the n d -dependence explicit γ [ Some representative Feynman diagrams used to compute the two-loop ADMs are shown in Figure 3.The two-loop ADM are well-defined only after specifying the basis of evanescent operators; this is done in App. A. All two-loop results in this section are new, to the best of our knowledge.The one-loop RG evolution has also been calculated in Ref. [22].

Threshold Corrections
When performing the running of the Wilson coefficients, we must also integrate out the heavy quarks at the respective scales.This leads to threshold corrections to the gauge couplings and quark masses, as well as to the Wilson coefficients.Below we describe the origin of the threshold corrections for each case, and collect the corresponding results.Due to the mixed RG evolution, the matching of the effective theories must also be performed as an expansion in αs and κ, cf.Eq. (11).We stress that the results presented below are only applicable for our specific case in which the only non-vanishing initial conditions are the ones in section 3.1.Other UV completions, with more contributions to the initial conditions, can receive additional threshold corrections (for instance, if C e 3 (µ ew ) ̸ = 0 at one-loop).
The Bottom-Quark Case In the case of an anomalous bottom-quark Yukawa coupling, q h = b, the only relevant threshold below the weak scale is at µ b in which we match the five-flavour version of the Lagrangian in Eq. ( 4) to the four-flavour version of the Lagrangian in Eq. ( 8).There are three effects that induce a non-trivial correction to the matching onto C e 3 : the threshold correction for α s when matching from the fiveflavour onto the four-flavour theories (the corresponding one for α e does not contribute because in our case C with γ 0 b,s = 6C F the leading-order QCD mass anomalous dimension, the α s threshold correction δα s = 2/3, and m b = m b (m b ) in the five-flavour theory in the MS scheme.We indicate by explicit subscripts "5fl" and "4fl" in which effective theory the various quantities are defined.The couplings in Eq. ( 27) are evaluated at µ b , i.e., αs,4fl (µ b ) and κ 4fl (µ b ).
The Charm-Quark Case In the case of an anomalous charm-quark Yukawa coupling, q h = b, there are threshold corrections both at the bottom-and the charm-quark thresholds, µ b and µ c , respectively.
At µ b , the effective Lagrangian in the five-and four-flavour theories are the same (see Eq. ( 4)), and at NLL accuracy the only relevant effect is the decoupling of α s .The only non-trivial matching condition at the bottom-quark scale then reads for the charm-quark case, q h = c, Other operators beyond O e 3 also receive threshold corrections.However, none of these terms enter the final three-flavour value of C e 3 (µ c ) at NLL order, so we do not list them here.At the charm-quark threshold, µ c , the threshold correction from matching the four-flavour onto the three-flavour theory with the single operator in Eq. ( 8) is analogous to the bottom case.Accordingly, we find for e,(02) 3,4fl (µ c ) + 24 3,4fl (µ c ) − 16 log where γ 0 c,s = 6C F , δα s = 2/3, and m c = m c (m c ) in the four-flavour theory in the MS scheme.

Analytic Solution of the RG
In this section, we show the final result for the electron EDM Wilson coefficients after implementing the mixed RG evolution (see discussion above Eq.( 11)) using the ADMs from section 3.2, and including the threshold corrections from section 3.3.We obtain obtain the exact analytical result for the contribution to the electron dipole moment up to O(α s κ 2 ) including the resummation of QCD logarithms: e,(02) 3 with α e , αs , and κ ≡ α e /α s evaluated at the scale µ q h in the four-and three-flavour theory for the bottom-and charm-quark case, respectively.The coefficients C e,(mn) 3 (µ q h ) are functions of the initial conditions at µ ew and of ratios of the values of α s at different scales.For brevity we introduce corresponding compact notations In both the bottom and charm-quark cases, the operator O 3 e in Eq. ( 8) below the respective heavy-quark scale is modified only by QED running effects which are negligibly small.

Bottom-Quark Result
In the bottom-quark case, q h = b, we find for the coefficients in Eq. ( 31 where m b = m b (m b ) in the five-flavour theory and in the MS scheme.

Charm-Quark Result
In the charm-quark case, q h = c, we find for the coefficients in Eq. ( 31)

C
e,(02) 3 where m c = m c (m c ) in the four-flavour theory and in the MS scheme.

Numerical Results
In this section, we present the numerical results based on the analytic expressions of the previous section and obtain constraints on CP-odd Higgs Yukawas to the bottom and charm quarks from the electron EDM.Combining the initial conditions in section 3.1 with the RG solution in section 3.4 leads to the electron EDM prediction, cf., Eq. (10), This result demonstrates explicitly how our RG-improved calculation removes the ambiguity of the fixed-order computation by resumming the large α s logarithms and thus defining the scale at which the masses are evaluated: in contrast to the fixed-order result in Eq. ( 3), d e is here proportional to m q h (m q h )m q h (µ ew )/M 2 h and the functions F (N)LL q h , which contain the α s resummation and do not depend on the large logarithms log m q h /M h .Expanding the LL result of Eq. (31) in α s (µ ew ) we recover exactly the log 2 x q term in Eq. ( 3).On the other hand, we cannot reproduce the π 2 /3 term in Eq. ( 3) since it must arise from the higher-order terms in Eq. ( 31) of order κ 2 α2 s , which are counted as NNLL in the QCD resummation.
The full expressions for F (N)LL q h are readily extracted from the results in section 3.For convenience, we give their numerical values for the special case of and µ c = m c (m c ).In this case, the logarithms in F NLL q h vanish and the only dependence is on η 5/4 .We find where we used the values η 5 = 0.506 and η 4 = 0.605 obtained by solving the mixed QCD-QED RG equations at two-loop accuracy using the numerical input in Table 1.
Table 1: Input parameters used in evaluating the low-scale Wilson coefficient C e 3 (µ q h ) and equivalently d e /e.All values are taken from Ref. [23]; running parameters are evaluated in the MS scheme.

Parameter Value
Parameter Value Next, we use the full analytic expressions to estimate the uncertainty in predicting d e , and provide the corresponding bounds on the anomalous CP-odd q h -Yukawa couplings.We estimate the uncertainty due to the truncation of the perturbation series in two ways: • The dependence on the matching scales cancels in our result to the order we calculated (α s κ 2 ).
The residual scale dependence is sensitive to higher-order terms in the perturbation series.Therefore, we evaluate the Wilson coefficient C e 3 at the fixed low scale m q h (m q h ), and separately vary all matching scales (µ ew and µ b for the bottom-quark case; µ ew , µ b , and µ c for the charmquark case).We fix all scales that are not varied to their "central" values (µ ew = M h and • There is a further ambiguity in our result that would only be resolved by a NNLL calculation: we can evaluate the NLL correction to d e , i.e., the whole term proportional to F NLL q h in Eq. ( 31), using either two-or one-loop values for all masses and couplings.We use this numerical difference as a further way to estimate the uncertainty; this difference effectively smears the lines obtained for the NLL scale variation, as described above, into the red bands shown in Figure 4.
As central values for C e 3 or equivalently d e at NLL accuracy we take the average of the maximal and minimal values obtained by all scale variations and differences as described above.Half of that difference is then assigned as the theoretical uncertainty associated with missing higher-order terms.This leads to The corresponding results are further illustrated in the two upper and two lower plots of Figure 4 for the bottom-and charm-quark cases, respectively.The plots on the left show the scale dependence of d e /e on µ q h at LL (dashed, black line) and NLL (red band) accuracy.The plots on the right show the corresponding dependence on µ ew .For the LL result we use one-loop values for all masses and couplings, and both one-and two-loop values for the NLL result, as discussed above.The comparison of LL and NLL transparently shows how the NLL computation presented here drastically reduces the large uncertainties associated to QCD corrections.In the left two panels, the variation of the matching scale at the bottom and charm thresholds is shown, respectively, while the right two panels show the dependence on the electroweak matching scale.In all plots, the dashed lines show the scale variation of the LL result.The scale variation of the NLL results is indicated by the red bands.The boundaries of the red bands are obtained by the evaluating the NLL scale variation in two ways, as discussed in the main text.Finally, the two gray dotted lines show the fixed-order results that have been used in literature so far; they correspond to evaluating the quark masses in Eq. (3) at the electroweak scale (the line marked "high") and at the heavy-quark threshold (the line marked "low"), and neglecting all other QCD corrections.
Furthermore, the gray dotted lines in Figure 4 show the naive result of the fixed-order calculation, Eq. ( 3), that has been used so far in the literature.The two lines correspond to using different values for the heavy-quark mass m q h : the line marked "high" corresponds to using the value evaluated at the electroweak scale (m q h (M h )), while line marked "low" corresponds to using the mass evaluated at the low scale (m q h (m q h )).The spread of these three lines illustrates the level of ambiguity in the fixed-order result.The NLL computation of the current work removes this ambiguity almost entirely.
The electron EDM has been indirectly constrained by measuring the EDM of polar molecules; the currently strongest bound is |d e | < 4.1 × 10 −30 at 90% confidence level [15].To perform a rough combination of this experimental constraint with our derived theory uncertainties we interpret the measurement as a Gaussian centered at zero and the above bound as the corresponding 90% confidence level (CL) interval, i.e., including negative values for d e .Adding the corresponding experimental "1σ" uncertainty in quadrature with the theory uncertainty we find based on an one-parameter χ 2 function in terms of κ q h sin ϕ The experimental bound on the electron EDM [15] (see also Ref. [24]) translates into strong constraints on new CP-violating phases in various extensions of the SM.This is significant because the SM contribution to the electron EDM is estimated to lie nine orders of magnitude below the current experimental sensitivity [25,26].In this context, it is important to point out that the current strongest bounds on the electron EDM are obtained indirectly via the measurement of the EDM of polar molecules; it has recently been shown that the SM contribution to those EDMs lies only about five orders of magnitudes below the current experimental limit [27].Moreover, in addition to the electron EDM operator, CP-odd semileptonic four-fermion operators contribute to the EDM of polar molecules.While such operators are generated in our model, their contributions are suppressed by roughly five orders of magnitude compared to the electron EDM, and can be safely neglected.CP-violating phases such as those in Eq. ( 1) appear in several well-motivated beyond standard model theories (see e.g.Refs.[12,13]), and the electron EDM is capable of placing stringent bounds on these phases; these have mainly been studied in an EFT approach [6,7,[9][10][11][28][29][30].As shown in Ref. [9], analoguous, albeit weaker bounds on CP-odd bottom and charm Yukawas can aso be obtained from hadronic EDMs.
Implicit in the electron EDM bounds is a large O(1) QCD uncertainty that has so far been neglected, even though it leads to sizeable ambiguities in the resulting constraints.Fortunately, since the electron EDM is a leptonic observable, this ambiguity can be removed systematically by a perturbative calculation, without the need for additional non-perturbative information.In this work, we have calculated the contribution to the electron EDM of CP-odd Higgs couplings to the bottom and charm quarks in RG-improved perturbation theory, summing the leading and next-to-leading large logarithms proportional to the strong coupling constant.This calculation has reduced the residual ambiguity in the bound to the level of a few percent, as discussed in detail in Section 4. The perturbation series shows good convergence, as expected for a leptonic observable.If, in the future, a non-zero electron EDM were observed, the error could be further reduced by summing the next-to-next-to-leading QCD logarithms, as well as taking into account the QED evolution of the electric dipole moment below the heavy-quark thresholds.

A Unphysical Operators
For completeness, we list here the "unphysical" operators entering our calculation in intermediate steps.They are called "unphysical" because they vanish either via the equations of motion (e.o.m.) of the quark fields for onshell external states, or via algebraic relations that are valid in d = 4, but not in d ̸ = 4.
A.1 E.o.m.-vanishing Operators These operators have matrix-elements that vanish via the e.o.m. of the quark field.The following two gauge-invariant operators enter our computation at the two-loop level: The covariant derivative acting on electron fields is defined as with Q e = −1 the electron electrical charge.

A.2 Evanescent Operators
Next we list the evanescent operators that enter our computation at one-and two-loop order.The leading-order ADM does not depend on their definition, but the next-to-leading-order ADM does.In the q-q-A sector we need the operator The evanescent operators required in the e-q h sector read ( The square brackets denote antisymmetrisation normalised as

A.3 Operators Related to the Infrared Rearrangement
The last class of unphysical operators arises because our use of infrared rearrangement breaks gauge invariance in intermediate steps of the calculation.At the renormalisable level this method generates two gauge non-invariant operators corresponding to a gluon-mass and a photon-mass term, i.e.
The mass, m IRA , is completely artificial and drops out of all physical results, and Z IRA,g , Z IRA,γ are two additional renormalisation constants [31].One-loop insertions of the dimension-five and dimension-six operators can induce further gauge-invariant, higher dimension operators that are relics of the infrared rearrangement.For our calculation, the only relevant one is

B Renormalization Constants
The following are the SM counterterms necessary for the calculation: Z (1,0) Z (1,1) Z (1,0) Z (1,1) Z (2,0) where we have organised the contributions according to and f denotes either a charged-lepton or a quark field.
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 write Z x→y where the subscripts x and y symbolize sets of Wilson coefficients, for which we use the following notation and standard ordering: The first 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 read Z (0,1) (1,0) P →E = 0.The remaining mixing of physical operators into evanescent operators is zero at one-loop.Furthermore, the finite part of the mixing of evanescent into physical operators is subtracted by finite counterterms [32].They read Z (0,1) The remaining finite mixing of evanescent into physical operators is zero at one-loop.Furthermore, we need the mixing constants of the physical operators into the operators arising from infrared rearrangement; they are found to be Z (0,1) All other mixing constants of physical into the IRA operators are zero.

Figure 1 :
Figure 1: Barr-Zee diagram that contributes to the electron EDM from a CP-odd heavy-quark Yukawa coupling with the Higgs (red square).

Figure 2 :
Figure 2: Examples of leading order (left) and next-to-leading order (right) Feynman diagrams that contribute to the initial conditions of the four-fermion sector in the EFT.CP-odd heavy-quark Yukawa couplings are indicated by red squares.

Figure 3 :
Figure 3: Examples of two-loop QCD (left), mixed QED-QCD (middle), and two-loop QED (right) Feynman diagrams entering the computation of the anomalous dimension matrix.Effective operator insertions are represented by the unshaded boxes.
(µ b ) = 0); the different normalization of the dipole operators O e 3 and O e 3 in the two theories, which leads to a factor of m b /m e ; and a threshold correction from one-loop insertions of O be2 in the five-flavour theory.Details can be found in the analogous discussion in Ref.[9].Taking all of these effects into account, we find for the electron dipole operator in the four-flavour theory at µ b (µ b ) , (µ b ) .
) where the heavy quark masses at the electroweak scale are given by their NLL running relationsm b (µ ew ) m b (m b ) The maximal residual scale dependence then provides our first uncertainty estimate on C e 3 or, equivalently, on d e .The ranges for the scale variations are chosen as µ ew ∈ [60 GeV, 250 GeV], µ b = [2 GeV, 8 GeV], and µ c = [1 GeV, 2 GeV].The scale variations are shown in Figure 4 and further discussed below.(The µ b variation is not explicitly shown for the charm-quark case, as it looks very similar to the µ c variation.)

Figure 4 :
Figure4: Residual scale dependence of the electric dipole moment induced by CP-odd bottom-quark couplings (upper two panels) and CP-odd charm-quark couplings (lower two panels).In the left two panels, the variation of the matching scale at the bottom and charm thresholds is shown, respectively, while the right two panels show the dependence on the electroweak matching scale.In all plots, the dashed lines show the scale variation of the LL result.The scale variation of the NLL results is indicated by the red bands.The boundaries of the red bands are obtained by the evaluating the NLL scale variation in two ways, as discussed in the main text.Finally, the two gray dotted lines show the fixed-order results that have been used in literature so far; they correspond to evaluating the quark masses in Eq. (3) at the electroweak scale (the line marked "high") and at the heavy-quark threshold (the line marked "low"), and neglecting all other QCD corrections.