Electric Dipole Moments in Two-Higgs-Doublet Models

Electric dipole moments are extremely sensitive probes for additional sources of CP violation in new physics models. Specifically, they have been argued in the past to exclude new CP-violating phases in two-Higgs-doublet models. Since recently models including such phases have been discussed widely, we revisit the available constraints in the presence of mechanisms which are typically invoked to evade flavour-changing neutral currents. To that aim, we start by assessing the necessary calculations on the hadronic, nuclear and atomic/molecular level, deriving expressions with conservative error estimates. Their phenomenological analysis in the context of two-Higgs-doublet models yields strong constraints, in some cases weakened by a cancellation mechanism among contributions from neutral scalars. While the corresponding parameter combinations do not yet have to be unnaturally small, the constraints are likely to preclude large effects in other CP-violating observables. Nevertheless, the generically expected contributions to electric dipole moments in this class of models lie within the projected sensitivity of the next-generation experiments.


Introduction
Despite the tremendous success of the Standard Model (SM), the arguments for the necessity of an extension are compelling. Specifically, Sakharov's conditions [1] require the presence of additional CP violation with respect to the SM. Assuming CPT invariance, electric dipole moments (EDM) are known to be highly sensitive to new CP-violating phases in new physics (NP) models. The contributions in the Standard Model are extremely tiny (e.g. d SM,CKM n (10 −32 − 10 −31 ) e cm, see e.g. [2,3] and references therein), with one exception: the gluonic operator O GG ∝ µνρσ G µν G ρσ gives in principle a contribution many orders of magnitude above the present experimental limits for e.g. the neutron; this is called the strong CP problem. To explain the absence of this contribution, typically symmetries are invoked, involving additional particles. The most famous example is the Peccei-Quinn mechanism [4], predicting the presence of axions [5,6]. While these have not yet been found in experimental searches, we implicitly assume in this work when discussing hadronic EDMs that the strong CP problem is solved by this or some similar mechanism.
The combination of the resulting tiny SM "background" and very strong experimental upper limits makes EDMs a well suited laboratory to search for NP, complementary to direct searches at e.g. the LHC and Tevatron as well as searches involving flavour-changing processes. The strong suppression in the SM is due to its very specific connection between flavour and CP violation, i.e. the Kobayashi-Maskawa mechanism [7]. When new sources of CP violation are included in NP models, usually large contributions are induced, specifically in models which contain flavour-blind phases. Therefore these models include typically an additional mechanism to keep them at bay. This in turn, as realized first by Weinberg [8], leads in a wide class of models to the situation that the dominant contributions actually stem from twoloop diagrams, when the additional loop allows to avoid strong suppression factors like masses of light quarks or small CKM matrix elements.
An attractive option for NP is provided by Two-Higgs-Doublet models (2HDM), due to their simplicity and their being the low-energy limit of various theories with a viable UV completion. In the most general version of the model, the fermionic couplings of the neutral scalars are non-diagonal in flavour and, therefore, generate unwanted flavour-changing neutral-current (FCNC) phenomena. Different ways to suppress FCNCs have been developed, giving rise to a variety of specific implementations of the 2HDM. In the past, mainly 2HDMs without new sources of CP violation have been considered, especially those with a discrete Z 2 symmetry [9,10]. Recently, however, there has been increased interest in models without this restriction, see e.g. [11][12][13][14][15][16][17][18][19][20][21] and also [22] for a recent review. Potentially huge EDMs used to be the main argument to discard these models. The critical reconsideration of this argument is one of the main motivations for the present work. We show in this article that while the present experimental limits impose strong bounds on the CP-violating parameter combinations, in models with an appropriate flavour structure they have not yet to be unnaturally small. However, large enhancements in other CP-violating observables are very strongly restricted by these bounds. Furthermore, the generic size for EDMs lies well within reach of the next-generation experiments, presently planned and some already in progress. These will therefore provide critical tests for this class of models in the coming years.
The direct observation of the EDM of a charged particle is very difficult, due to the presence of a hugely dominating "monopole" contribution, i.e. its charge. Therefore, the most sensitive measurements, at least so far, stem from neutral systems, especially neutrons and atoms/molecules. Relating them to fundamental parameters involves complex calculations at different scales, often implying large uncertainties. Without their careful estimate no reliable constraints on NP parameters can be obtained. We start therefore in the next section by giving model-independent expressions for these observables in terms of Wilson coefficients of the relevant effective operators, taking recent developments into account and estimating the uncertainties of the QCD, nuclear and atomic calculations in a conservative manner. For a subset of systems, this has been done very recently by one of us in [23]; these results are used when appropriate. This is followed in Sec. 3 by a quick description of the experimental situation, after which we proceed in Sec. 4 to discuss the situation of EDMs in 2HDMs with new sources of CP violation. We start by describing the various sources, pointing out their different importance. To be specific, we then calculate the resulting constraints in the Aligned Two-Higgs-Doublet Model (A2HDM), which has been introduced in [15,17] and whose phenomenology has been further discussed in [24][25][26][27][28]. However, the structure of the model is such that the results hold rather generally. In this context, we point out a general cancellation mechanism for neutral scalar contributions, which questions the way they are commonly treated in the literature. We analyze the phenomenological constraints coming from the presently available experimental bounds in Sec. 5, before giving our conclusions in Sec. 6.

Comparison to existing work
There is a huge amount of literature on EDMs, and there is no hope of reviewing it here; instead, we refer the reader to [2,[29][30][31][32] for recent reviews. Generally, most of the analyses in the literature are performed within the framework of supersymmetric models (SUSY) (for recent examples, see [33][34][35][36][37][38][39][40][41] and also the phenomenological analysis in [30]). While in principle the 2HDM contributions are present in these models as well, they are usually subdominant, which is why they do not receive much attention. Especially the charged Higgs exchange is usually negligible in these models, as it does not exhibit the strong tan β-enhancement of other terms, which is why some of the corresponding contributions discussed below are not incorporated at all in these analyses.
Recent studies more closely related to our work include [42][43][44]. In the first of these, the authors discuss one contribution discussed below, namely the charged Higgs contribution to the neutron EDM. The results are similar to ours 1 , apart from a different treatment of the hadronic matrix element, which yields weaker constraints in our case. The second article discusses EDM contributions in the context of Minimal Flavour Violation (MFV), including complex phases in that framework. The authors perform the analysis in the decoupling limit and assume a small breaking of the Z 2 symmetry, as was assumed already for the 2HDM analysis in [45]. Their results are therefore relevant for a subset of our parameter space. They conclude, as we will below, that one-loop contributions are generally not exceeding the experimental limits. In addition, they consider a subset of the two-loop contributions we discuss below, corresponding to the more restrictive assumptions they make. Finally, in [44] the authors discuss a subset of MFV operators which might generate a new phase in B s mixing; the corresponding operators are not relevant in our context.

Model-independent expressions for EDMs
From the point of view of particle physics, the proper starting point for a model-independent analysis is the following effective Lagrangian at the hadronic scale (here up to dimension six, see e.g. [2]): with the operator basis The operators in Eqs. (2) are the (colour-)EDM operators O γ,C f for light fermions (f = e, d, u, s), the Weinberg operator O W and T-and P-violating four-fermion operators O f f without derivatives (see, e.g., [46]). The factors of 1/2 for the (C)EDM operators are included to identify the coefficients d γ,C f with the classical electric/gluonic dipole moment in the corresponding limit. The analysis of their influence on experimental observables is divided into two steps: first, the observables have to be expressed in terms of the coefficients of this effective lagrangian. This step can be done independently from the NP model considered and is performed in this section. The necessary calculations are on the QCD, nuclear, and/or atomic/molecular level. They typically involve relatively large uncertainties. Their careful assessment is essential to obtain reliable bounds, which is why we will pay close attention to this. In the second step, performed exemplarily for the A2HDM later in this paper, the coefficients have to be calculated in terms of parameters of the NP model considered, allowing to obtain the constraints on the latter. In calculations on the QCD level, the corresponding matrix elements are often known only up to a factor of a few, sometimes without a definite sign. There are different methods to calculate/estimate them; while Naive dimensional analysis (NDA) [47] is still used occasionally, mostly due to its simplicity, its estimates are known to be uncertain e.g. by arbitrary powers of 4π (see e.g. [48]), which is why we do not consider these estimates here. Instead we are going to use QCD sum rule estimates, where such factors are absent and which are supposed to be uncertain "only" by the aforementioned factor of a few (depending on the operator). The main reason for this limited precision is that for sum rules with baryons the suppression of excited states does not work as well as for mesons. For a review on these issues, see [2]. While ultimately progress may come from Lattice QCD, there are severe difficulties obtaining reliable results at the moment, such that we are not aware of available results competitive to the ones used here. Note also available calculations in the framework of Baryon Chiral Perturbation Theory [49] (see e.g. [32,50] for recent analyses and references therein). There, the scaling of the various matrix elements can be analyzed from the chiral properties of the operators, leading to a systematic classification scheme. However, it is typically accompanied by NDA estimates, as unknown low-energy constants prevent quantitative estimates. We therefore do not use their results here quantitatively.
The involved calculations on the nuclear and atomic/molecular level are in very different shape, uncertainties ranging from a few to several hundred percent; these are commented upon in the appropriate subsections. These relatively large uncertainties are in a sense of minor importance for the experimental searches, since NP contributions can easily be larger than the SM ones by several orders of magnitude. However, they are essential in obtaining bounds for NP parameters from the experimental limits. Readers not interested in their detailed discussion find the final expressions for the corresponding EDMs in Eqs. (3),(5) and (31). For paramagnetic systems, we use the results from [23], but perform an update to include the very recent measurement with thorium monoxide (ThO) [51]; its results are summarized in Table 1.

The neutron EDM
The neutron EDM can be related to the coefficients in Eq. (1) by QCD calculations alone. Here we collect the necessary formulae, for details see again e.g. the review [2]. This EDM is dominated by contributions from the (C)EDMs of its constituents and the Weinberg operator, while four-quark operators play a minor role.
The QCD sum rule calculation for the contribution from the quark (C)EDMs yields [52,53] where µ h ∼ 1 GeV denotes a hadronic scale. In the following we suppress the scale dependence in the notation for brevity and evaluate at µ h = 1 GeV unless stated explicitly. The uncertainty given here for these matrix elements is similar to the estimate given in [52]. However, given the results in [53], we extended the range to include lower values. 2 This incorporates larger values for the normalization factor λ n , determined by the matrix element of the nucleon and its interpolating current, see [54][55][56]. 3 We note that alternative treatments are compatible within the estimated level of precision, however indicating in some cases higher sensitivity, see e.g. [59]. Note furthermore that the quark condensate qq in this formula combines with the light quark masses in the Wilson coefficients as [60] ( with the sign left undetermined. This expression is based on several estimates that all lead to similar results, but is not a direct calculation. Finally, for the sake of completeness, for an exemplary four-quark contribution to the neutron EDM the sum rule estimate results in [62] |d n (C bd )/e| = 2. 6 again with an unspecified sign. Note that here four-fermion operators involving the beauty quark, defined analogously to their equivalents with light fermions, contribute below the b-quark mass scale effectively via an effective two-gluon coupling of the down quark, which is also why the coupling is to be evaluated at µ b ∼ m b . The contribution from up-type quarks is ignored, as enhanced couplings in that sector (corresponding e.g. to tan β 1 in a Type II model or to |ς u | 1 in the A2HDM) are usually excluded.
2 Note that the analytical differences have minor numerical impact. 3 Note, however, that the rather large central value in [54] does lead to a too small value for the nucleon sigma term σ πN when using the sum rule in [57] [58]. 4 Here and in Eq. (6) the authors state a 100% uncertainty for the result, which we incorporate as allowing for twice and half the computed value.

EDMs of atoms
For atoms, Schiff's theorem [63] implies a vanishing EDM in the non-relativistic limit for systems of particles whose charge distribution is identical to their EDM distribution. The limits from the nonobservation of these EDMs are then related to violations of the conditions for this theorem, and separated into two classes, depending on which of the approximations is more strongly violated. For reviews on atomic calculations, see e.g. Refs. [29,46].
In paramagnetic atoms, i.e. atoms with non-vanishing total angular momentum, relativistic effects are important, which are largely enhanced for atoms with a large proton number [64][65][66], scaling at least like d ∼ Z 3 . This implies a sensitivity mainly to the electron EDM, but also electron-nucleon interactions are enhanced, described by The coefficients of both classes of contributions are estimated in atomic multi-body calculations. In some publications, these operators are classified instead according to their isospin, where X = S, P, T and the Dirac structures Γ 1,2 X can be read off from Eq. (7). In diamagnetic atoms, where the total angular momentum vanishes, the finite size of the nucleus is the main source for the violation of Schiff's theorem. The dominant contribution to the corresponding EDM stems from its nuclear Schiff moment, which can be expressed in terms of the nucleon EDMs and pion-nucleon couplings, which are in turn related to the basic terms in Eq. (1). Although the quark CEDMs typically give the dominant contribution, the above electron-nucleon interaction is relevant as well.
For an atom with proton number Z p , neutron number Z n and consequently nucleon number A = Z n + Z p , the parameter combinations effectively contributing to the EDMs read (see again e.g. [29,46]) where σ N at denotes the sum over the spin of the indicated nucleon species in the corresponding nuclear state, and we used σ at = σ n at + σ p at and σ i at = σ i at I/I, where I denotes the total nuclear spin. The spin sums stem from the quantum-mechanical expressions derived from the pseudoscalar operator N γ 5 N . In this equation, the contribution from the first term in Eq. (7) is seen to be additionally enhanced, because the contributions from neutrons and protons enter spin-independently. This renders this term dominant for paramagnetic systems, as for the other two coefficients closed shells in the nucleus barely contribute. In diamagnetic atoms, it does not contribute at leading order, however, which is why the relative influence of the other two terms is relatively enhanced. In fact, if present, among the electron-nucleon interactions the third term is typically dominant in this case. In general, the definitions forC (at) X imply a dependence of these coefficients on the system considered. However, because of (Z n + Z p )/A = 1 andC n S ≈C p S , this is usually neglected in the case ofC S . More importantly, the ratios Z N /A are approximately universal for the systems considered here, leading to a universalC S even forC n S =C p S [23]. However, for the spin-dependent terms the relative weights are not atom-independent, such thatC at P,T depend on the atom ifC n P,T =C p P,T . To remind the reader of that fact, we added the label 'at' on the corresponding quantities.
Expressed in terms of the isospin coefficients, the effective contributions correspond to Note again that the coefficient of the triplet contribution is neither atom-independent nor generally small in the latter case; for example, σ p Xe ≈ σ n Xe /3 and σ p Hg ≈ σ n Hg /10 [67], implying ( σ n at − σ p at )/ σ at ∼ 1 for the latter. Note furthermore that the coefficient for C P is sometimes mistakenly given as (Z n − Z p )/A. ThO only,C S = 0, 90% CL [51] 0.089 × 10 −27 e cm †, ‡ 0.6 × 10 −8, ‡ Table 1: New limits on the electron EDM andC S , including the measurement in the ThO system [51], see text. † : Using W d from [69]. ‡ : Theory errors neglected.

The EDM of paramagnetic systems
The EDM of paramagnetic systems is dominated to very good approximation by the contributions from d e andC S , as explained above. The presently most constraining measurement from this class is performed with ThO [51]. Their result is given in terms of an angular frequency, corresponding to an energy shift, which can be parametrized as using the conventions from [23] for the atomic constants W d,c . For these, we obtain W d = −(3.67±0.18)× 10 25 Hz/e cm from [68][69][70] and W c = −(598 ± 90) kHz [69]. Note that the calculations for the former are consistent; we use the value given in [70] (corresponding to E eff = 75.6 GV/cm), and enlarge only slightly the uncertainty to 5% due to the rather large difference to the central value in [69]. Note furthermore the consistency within uncertainties between the explicit calculations and the analytical estimate in [71].
Since each measurement only constrains a combination of the two contributions, conservatively no constraint on the electron EDM alone can be obtained from any single measurement. The combination with previous measurements, performed with thallium (Tl) atoms and ytterbium fluoride (YbF) molecules [72,73] allows for a model-independent determination of the electron EDM, which improves significantly using in addition information from the mercury (Hg) system [23], see Fig. 1. However, in contrast to the situation in [23], the Hg measurement does not provide a competitive bound onC S compared to the ThO one when setting d e to zero. Therefore, this procedure results in an extremely conservative bound, which is allowing for arbitrarily large cancellations between the different contributions and includes conservative estimates for the uncertainties of the various coefficients. In addition to this value, we obtain additional ones using assumptions on the maximal amount of fine-tuning: we restrict the contribution fromC S alone not to exceed n = 1, 2, 3 times the measured limit for ThO and use this as an additional constraint, thereby using effectively the ThO result twice. While this is clearly not as rigorous as the above limit, it is still more conservative than the common procedure to set the contribution from the electron-nucleon interaction simply to zero. This yields the inner solutions in Fig. 1; the corresponding upper limits for d e andC S are given in Table 1, together with the values quoted in [51], which are obtained by setting the other contribution to zero and neglecting theory uncertainties. Note that, with a second competitive measurement, d e andC S can be extracted again without additional assumptions, see again [23]. In the phenomenological section below, we use all values presented in Table 1, in order to demonstrate the progress due to the new measurement and to compare the various upper limits. We consider n = 2 already a conservative choice, since there is no dynamical relation between the two contributions, rendering large cancellations unlikely. Nevertheless, the necessity to introduce this kind of assumption demonstrates the importance of independent measurements in other systems, ideally with strongly differing values for the ratio W d /W c like, e.g., rubidium.

The EDM of diamagnetic atoms
For diamagnetic atoms mainly finite-size effects of the nucleus determine the EDM. More specifically, its main source is the CP-odd nuclear Schiff moment 5 [63]. Although contributions from the nucleon EDMs are present as well, it is dominated by T, P -odd nuclear forces. These are represented by the interference of CP-even and -odd pion-nucleon interactions, the latter of which depend on the CEDMs of the up and down quark and four-quark operators (see again e.g. [29]). All of the necessary calculations are very involved, and the wide range of results indicates that the related theoretical uncertainties are large; for recent discussions see,e.g., [32,34]. The first step, namely relating the atom EDM to the Schiff moment, is parametrized as with the constant C at Schiff being the result of multi-particle computations, modeling the electron-nucleon configurations in the corresponding atom. Due to the recent measurement in [78], the interest in these calculations has increased especially for Hg, leading to two recent results [67,79], 6 from which we infer which is now in agreement with the updated value of [79] (the preliminary result reads C Hg Schiff = −2.46 [80]), strengthening the confidence in these calculations. The value also agrees with the earlier calculation [81] and is reasonably close to an old estimate [82].
In the next step, the Schiff moment is related to the CP-odd and -even πN N coupling constants [83], parametrized as [84] (note the different sign conventions for these constants used in the literature) The isotensor coefficient is set to zero in the following, as its effect is suppressed by an additional factor of the mass difference of light quarks [85]. The CP-even coefficient is given by g πN N = 13.17 ± 0.06 [86], the uncertainty of which is negligible in this context. The corresponding nuclear calculations for mercury span a wide range and have in the case of a 1 also different signs in some of the calculations, see Table 2. While in principle the calculations in [84] are more advanced than the previous ones, for mercury at some stage all the interactions used show problems, and the differences between the calculations are not well understood [84]; in absence of errors in one or several of the calculations, the problem might stem from the fact that mercury is a soft nucleus [84]. We therefore estimate conservatively the following ranges: covering the full range of results shown in Table 2. We note that the possibility of vanishing a 1 implies that no constraint can be obtained conservatively on the corresponding isovector combination of CEDMs. Ref.
a 0 e fm 3 a 1 e fm 3 b e fm 3 [87] 0.00004 0.055 - [88] 0.010 (0.002-0.010) 0.074 (0.057 -0.090) -  [34] in that the SIII Skyrme interaction results in [84] were the most trustworthy of their calculations. In fact, it is shown in Ref. [88] that this interaction yields the worst results in reproducing the observables which can be used as experimental crosschecks, which is why the authors of [84] themselves regard it as critical.
Below we will show results for a representative value, in order to illustrate the potential of this observable, given a more reliable theoretical situation. Regarding the coefficient a 0 , we point out that the tiny value obtained in [87] might be the result of accidental cancellations, see the discussion in [88]. Finally, the parameter b has so far only been calculated by one group; given the unclear situation, an additional independent calculation would be worthwhile. In the last step, the CP-odd πN N -couplings have to be related to the (C)EDMs of quarks. For this, typically a relation from the partially conserved axial current is used for the pion and QCD sum rules for the remaining nucleon matrix elements of quark currents. The main difficulty in this case is that for baryon sum rules in external fields the Borel transform does not exponentially suppress the contributions from excited states, leading to a large uncertainty. For the isovector coupling, this can be improved by tuning an unphysical parameter to suppress these higher order terms, leading to [85] g (1) πN In the isoscalar sector a similar tuning is not possible, allowing for [85] g (0) thereby also questioning the sensitivity to this combination of CEDMs. An additional contribution to g (1) πN N comes from four-quark operators, reading [33,62,89,90] where naive factorization has been used for the four-quark matrix elements. Here κ parametrizes the uncertainty in the strange quark content of the neutron, we use σ πN = N |m uū u + m dd d|N = (59 ± 7) MeV [92], andβ = 11 − 2n l /3 for n l light quarks [90]. Recent lattice studies [93][94][95][96] (see also [97][98][99][100][101]) indicate a smaller value for κ than assumed previously (see e.g. [102] and references therein), thereby reducing the influence of the strange quark on EDMs. However, while agreeing on a smaller order of magnitude, the range implied by these calculations is still relatively large. We combine them to arrive at where we again chose a conservative range for the central value, reflected by the second uncertainty, while the first one is of statistical origin. However, as for the neutron, the four-quark contributions are subleading in 2HDMs, see the discussion in Sec. 4.
The Schiff moment receives contributions from the nucleon EDMs as well. While this contribution is not expected to be dominant, the resulting constraint for the neutron EDM is actually of the same order like the one from the dedicated experiments; using e.g. Eq. (14), the range for C Hg Schiff given in Eq. (15) and the expression S(d n ) = 1.9 fm 2 d n [103] (for simplicity with its central value), we obtain |d n | ≤ 7.8 × 10 −26 e cm, which is only about a factor of two weaker than the present direct limit [104]. However, there is no way to combine these limits, therefore we just consider the latter in the following.
Additional sources from electron-nucleon interactions and the electron EDM are present. Regarding the latter, the value usually used in the literature for mercury reads d Hg (d e ) = 1.16 × 10 −2 d e [105]. The corresponding calculation, however, shows a high sensitivity to higher order effects; the "corrections" to a previous estimate [106] amount to ∼ 200% and change the sign. The authors point out the sensitivity to correlation effects (which have been found to be large for mercury for the coefficients discussed above), making a new calculation mandatory. In light of this situation we do not see a way to extract a meaningful upper limit on the electron EDM from mercury until the theoretical situation improves. However, even taking the central value quoted above, the bound would be weaker than the one from paramagnetic systems.
The electron-nucleon interactions are induced via the three operators in Eq. (7). In this case, thẽ C S contribution is suppressed, as to leading order its contribution from closed electron shells vanishes; generically this leads to a dominance of the term proportional toC T , if it is present. However, for the 2HDMs considered here, only the scalar-pseudoscalar operators are present. The contributions proportional toC P are often neglected, as its coefficient is one order of magnitude smaller than the one ofC S , even in this case, due to the suppression by the nucleon mass. However, expressing the corresponding matrix elements in terms of coefficients of the four-fermion ones shows basically the opposite behaviour, rendering the sensitivity to fundamental parameters similar. All types of contributions are relevant in some part of parameter space [33].
Given the large theoretical uncertainties in the contributions to the Schiff moment, the constraints on the electron-nucleon interaction might be the most important one at the moment. The coefficients in the relation d Hg (C S,P ) are obtained again in atomic calculations; usually only the coefficient of the tensor operator is calculated and approximate analytic relations are used to obtain the others 8 [29,67,106,107]: where µ denotes the magnetic moment of the nucleus in terms of nuclear magnetons µ N . We expect the uncertainty for these relations to be relatively small, as also indicated by a recent explicit calculation for a variety of atoms [67], which is why we neglect it in the following. For the tensor coefficient, parametrized by recent results read C Hg C T = −5.1 [67] and C Hg C T = −4.3 [79]. Thus, using Eqs. (24) and (25), we obtain where we used µ Hg = 0.506 µ N 9 and σ = −1/3 I/I, the estimate from a simple shell model for the nucleus, and the common convention d = dI/I. The next step is to relate the coefficientsC S,P to the effective operators discussed above. The contributing operators are four-fermion operators with electrons and light quarks, and an effective electron-2-gluon vertex from integrating out the heavy quarks. Again neglecting the up-type quark contributions, they can be parametrized as follows [33,62,89,90]: where the same matrix element appears as for C qq , see Eq. (20). The missing ingredients are the expressions forC P in terms of the coefficients of four-fermion operators. We use the estimates for the isospin coefficients (cf. Eqs. (8) and (11) again neglecting up-type quark contributions. Finally, from these considerations we obtain the following result for mercury: with the expressions forḡ (1,0) πN N given in Eqs. (18), (19) and (22). The possible vanishing of the coefficients of the isoscalar and -vector CEDM contributions implies that conservatively no bound can be obtained for them. Usually these contributions are assumed to be the dominant ones in this system, underlining the importance of theoretical developments to clarify the situation. Below, we will show the limits that would result for the central values in Eq. (31), however only for illustration purposes.

Renormalization of the effective operators
To connect the relevant Wilson coefficients at the hadronic scale with the short-distance calculation at the electroweak one, the renormalization group running has to be taken into account. In general, QCD effects tend to reduce the value of the different coefficients (see e.g. [108]), 10 apart from the four-quark one [109,110]. We neglect its mixing into the CEDM because of its smallness; however, we take its enhancement into account in the estimate below. As pointed out in [108], the mixing of the CEDM-into the EDM operators constitutes a large effect. On the other hand, we consider the NLO running of minor importance at present, given the large uncertainties in the hadronic matrix elements. Furthermore, the mixing of the Weinberg operator into the CEDM ones is of higher order in α s and therefore neglected in the following. For models, in which this is not the case, the operator mixing has recently been discussed in [110]. Denoting C = (d γ f /2, d C f /2, C W ) (see Eq. (1)), this results in the following leading order expressions [108,[111][112][113]: where β 0 = (11N C − 2n f )/3, N C = 3, C F = 4/3, n f denotes the number of active flavours and q q = 2/3, −1/3 the charge for up-and down-quarks, respectively. As we expect the Higgs masses to be of the order of m t (as is the mass of the already observed scalar), we choose µ tH ∼ m t as the common matching scale where top quark and scalars are integrated out. We use the solution to Eq. (32) to scale down to µ ∼ m b , where in addition the beauty quark is integrated out, thereby matching O C b onto O W . The matching condition reads [112,114]   where µ + b (µ − b ) refers to the same scale µ b , but in the n f = 5(4) theory, respectively. We emphasize that this matching, together with the larger anomalous dimension of the Weinberg operator, implies a relative enhancement of the contribution involving charged-scalar exchange compared to the one involving neutral scalars, as the suppression from the running is weaker for the CEDM contribution. When going to the n f = 3 theory, the charm contribution to C W becomes local, which is however severely suppressed because of m c m t , and therefore neglected. The solution of Eq. (32) reads where we introduced η i−j = α s (µ i )/α s (µ j ), η = η t−h , and κ i = γ i /(2β 0 ). For the sake of simplicity, Eq. (35) is displayed for constant n f throughout the integration, but its change is taken into account in the numerical analysis. Finally, regarding the Wilson coefficients of the semileptonic four-fermion operators, we note that they scale like the quark masses, therefore the combinations C qe /m q and C eq /m q are scale-independent.

Experimental status
At present, the limits most sensitive on the various sources discussed above stem from searches for EDMs of Tl [72], YbF [73], ThO [51], Hg [78] (see also [115] for a more detailed discussion) and the neutron [104], see Table 3. The physical origin of their EDMs is quite different, making them complementary sources of information. Although these limits have different orders of magnitude, their different dependences on the fundamental parameters of the theory actually lead to similar sensitivities.
Several developments allow to expect significantly improved bounds or a non-zero measurement in the near future, see also e.g. [2,[29][30][31][32]. The first option is to reduce the uncertainties within the established methods, but in the longer term techniques exploiting different features like octupole deformation hold the promise of qualitatively improving the sensitivities further. Regarding octupole deformation, important experimental progress has been reported recently in [133].
There are several experiments for the neutron EDM planned and running or under construction (see [134] for a recent result and again Table 3) using different techniques to obtain higher neutron densities to achieve an up to two orders of magnitude improved bound.
Regarding the electron EDM and the electron-nucleon coefficientC S , with the experiments for thallium completely dominated by their systematic uncertainties, significant advancement seems difficult within this system. An improvement, again up to two orders of magnitude, might come instead from the cesium, rubidium and francium systems [122][123][124][125][126], which can be stored to obtain longer oscillation times. The expected limits correspond to probing the electron EDM down to 10 −29 e cm in the midterm future (2-3 years), and even 10 −31 e cm has been envisaged for the farther future in [125]. Further measurements with YbF are expected to strengthen the present limit in the short term [73] and various experiments are underway to gain sensitivity down to ∼ 10 −30 e cm or further [127,128] (see e.g. [31,32] for a more complete list). A key technique is the rejection of systematic errors by using the so-called Ω-doublet structure of a subset of paramagnetic molecules (characterized by two very closely lying states of opposite parity, leading to an extremely high polarizability), as demonstrated in the recent experiments [135] -so far obtaining a less stringent limit than the one from YbF -and [51].
In the context of the analysis in [23], the expected presence of several measurements with similar sensitivities will allow to improve model-independently the limits on the electron EDM as well as the constantC S , taking into account possible cancellations and at the same time removing the input from the Hg system and assumptions on fine-tuning.
In the future, trapped molecular ions might also be used as sensitive probes for EDMs, however, at the moment there are still severe experimental and theoretical challenges to overcome. Finally, also solid state systems are being explored as sensitive probes for the electron EDM [136,137]. While again some experimental as well as theoretical progress is necessary before competitive results can be achieved, recent results show the progress in this field [138].
For diamagnetic systems, apart from some improvement from the Hg system itself [78,129], significant improvement is aimed at using xenon ( 129 Xe), for which the theoretical treatment is similar to the one described above. Further progress is expected with different enhancement mechanisms like intrinsic octupole moments and, related to that, closely neighboured parity doublets, which allow for large enhancement factors for the corresponding Schiff moments. Prominent examples are radium and radon; however, the calculation of the corresponding matrix elements is more complicated, making again theoretical uncertainties a critical issue. Furthermore, also diamagnetic molecules are under investigation. A first measurement exists in the TlF system [139], but the planned experiments are expected to improve greatly on the present limits, see again e.g. [31] for a recent list of experiments. Generally, due to the various possible contributions to the EDMs, measurements in different diamagnetic systems are even more necessary to disentangle the sources and potentially differentiate between NP models. Ultimately, this could be done in an analysis similar to the one in paramagnetic systems [23], but for that a lot more information than presently available is necessary.
Finally, new techniques are being used for measuring the EDMs of charged particles directly by using a storage ring [140][141][142][143], e.g. for muons, where the present limit stems from a storage ring experiment already [144], the proton [145], which is supposed to be tested down to 10 −29 e cm, or the Deuterium nucleus, which has the advantage of being lightly bound and allowing thereby to circumvent the large uncertainties present e.g. in the nuclear calculations for mercury. There are also proposals to use the technique for molecular ions, see e.g. [146].

EDMs in 2HDMs
We now address the model-dependent second step in relating EDMs to model parameters, i.e. calculating the relevant effective coefficients in specific models. The model dependence is in some sense more severe in EDMs than for other observables, for the following reason: as generic one-loop contributions are excluded already, an additional mechanism is necessary to render them small. This implies that the usual power counting is not sufficient, but that this suppression mechanism has to be incorporated. As a result, even if a NP model has a 2HDM as an intermediate effective theory, this does not necessarily imply that limits calculated at that level hold for the full theory, as can be seen e.g. in many SUSY models.
In this section we limit the discussion to 2HDMs. While generally even for these, limits cannot be given model-independently, we hold the discussion as general as possible, and the results, while given in the parametrization of the A2HDM, can be easily transferred to other frameworks.

Contributions to EDMs in 2HDMs
We start by listing the contributions to the different effective operators in Eq. (1) within a 2HDM. As most recent analyses have been done within a SUSY framework, we will comment on the differences to the situation there when appropriate.
• Four-fermion operators: They induce the leading contributions in the SM [2], but there their effects remain extremely small. In 2HDMs, they are induced by CP-violating Higgs exchanges. While they can have contributions at tree level, in that case a further suppression by two light-fermion Yukawa couplings applies. If these are proportional to (or of the order of) the corresponding masses (as e.g. in models with a Z 2 symmetry, the A2HDM, MFV, Type III, . . .), the ones with light fermions are suppressed to an acceptable level. The proportionality implies also that the induced coefficients divided by the corresponding masses are family-universal. Those involving heavy fermions do not contribute directly, but induce higher-dimensional operators like (f f )GG, again on an acceptable level, cf. Eqs. (6) and (22).
There are two categories: CP-violating four-quark operators contribute to the nucleon EDM directly, or to the Schiff moments of nuclei by inducing CP-violating meson-nucleon interactions. As we will show below, for Higgs couplings of the order of the fermion masses, in the 2HDM both contributions are subleading and can be neglected. In SUSY, they can receive contributions proportional to tan 3 β from threshold corrections, rendering them more important there and even dominant for very large values of tan β [62,147].
The second category consists of semileptonic operators. These induce CP-violating electron-nucleon couplings in atoms, as discussed in Sec. 2.2. While in principle they are as suppressed as their fourquark equivalents, they receive very strong enhancement in heavy atoms due to the large number of nucleons and electrons, making their inclusion mandatory.
• Weinberg operator: The contribution to this operator starts at the two-loop-level in 2HDMs, schematically shown in Fig. 2(a). It is neither suppressed by small quark masses nor by small CKM elements, and therefore is expected to be large. However, for two reasons it is not completely dominating: first, the matrix element given in Eq. (5) is of the order of a light quark mass instead of a typical hadronic scale, and second the RGE running yields a strong suppression, see Eq. (32). As mentioned before, the second point is also the reason why, contrary to naive expectations, the neutral Higgs contribution is generally suppressed compared to the charged one, cf. Sec. 2.3. In SUSY, the graphs discussed here are typically subleading, which is why they are often ignored.
• (C)EDMs of light quarks: In the SM they vanish at the one-and even the two-loop level, leading to a tiny result [148]. In a general 2HDM, however, they can be generated at the one-loop level and are by far the leading contributions, which is why an additional mechanism for their suppression is necessary. An example are models with a Z 2 symmetry, where these loops are CP-conserving like in the SM. In the A2HDM, but also more generally in models where the Higgs couplings are related to CKM-matrix elements and quark masses, the one-loop contribution for the light fermions is suppressed by at least one corresponding mass factor, together with factors like m 2 , rendering them one to two orders of magnitude smaller than the contributions discussed in the following. The reason is that the latter factors are circumvented in Barr-Zee(-type) diagrams [149][150][151], see Fig. 2(b), which is why these two-loop contributions dominate in this class of models. The corresponding contributions are given later in this section. They receive contributions from neutral scalars, discussed in the above papers, but also from charged ones. These contributions have been discussed partly e.g. in [152][153][154], but to the best of our knowledge e.g. the ones for the down-quark EDM with a top and beauty quark in the loop are still missing. To construct a Barr-Zee diagram with a charged Higgs, a second charged current is necessary; therefore there are no contributions to the CEDMs from these graphs.
In SUSY, there are more one-loop contributions from loops with gauginos and sfermions, generally leading to strong bounds on the imaginary parts of the corresponding couplings. From the Higgs sector, the two-loop contributions again dominate, due to the arguments given above.
• Electron EDM: The SM contribution to this is tiny, as for m ν → 0 it vanishes even on the three-loop level [155]. In 2HDMs, the one-loop contributions are real unless a neutrino coupling is involved, which is why the dominant contributions are again on the two-loop level, from Barr-Zee diagrams. In SUSY, already on the one-loop level sizable contributions appear from gaugino-slepton loops, therefore again the Higgs contributions do not receive that much attention.
Because of the arguments given above, we will explicitly consider only the contributions stemming from the two-loop diagrams and the semileptonic four-fermion operators important for atoms and molecules. It should be emphasized again that the limits obtained within 2HDMs are sensitive to the UV completion of the model, as their sensitivity to two-loop contributions already shows. Especially in SUSY there are usually large one-loop contributions dominating, which are not included here.
The contributions listed above are related to different sources of CP violation in 2HDMs: while e.g. the charged Higgs contribution to the Weinberg operator stems only from CP violation in the Yukawa couplings of the model, diagrams involving neutral scalars in general receive contributions from the Higgs potential as well. Before providing results for specific diagrams, we discuss the different classes of contributions, pointing out their general features.

Charged Higgs contributions in 2HDMs
The Lagrangian for charged Higgs exchange can for vanishing neutrino masses be parametrized as where V is the CKM mixing matrix and the form reflects the fact that we will be mostly concerned with the A2HDM, where ς u,d,l are complex numbers of O(1); for a general 2HDM, they are arbitrary matrices and the dependences on the quark mass matrices are artificial, i.e. just a possible normalization. Note that we consider the latter form simply as a phenomenological parametrization. For the contributions calculated explicitly below, this implies only the generalization of the factors ς u,d,l . However, when these elements are indeed arbitrary, other contributions are possibly dominant, since especially the suppression for the one-loop contributions to (C)EDMs explained above can be spoiled. If the scaling does approximately hold, the bounds obtained below are valid for the corresponding generalized couplings.

Neutral Higgs contributions in 2HDMs
The flavour-diagonal Higgs couplings are parametrized analogously as with the fields ϕ 0 i = {h, H, A} denoting the neutral scalar mass eigenstates. Introducing the notation F (f ) for the species of a fermion, e.g. F (u) = F (c) = F (t) = u, we write the fermion couplings to neutral scalars as to allow for the general form of ς u,d,l . Here, R is an orthogonal matrix defined by M 2 diag = RM 2 R T , relating the mass eigenstates to the neutral scalar fields S i in the Higgs basis, where Φ T , and the fields G i are the Goldstone bosons absorbed by the gauge bosons: ϕ 0 i = R ij S j . In a general 2HDM, the ς u,d,l are the matrices introduced in Eq. (38), only the diagonal elements of which are relevant here.
The fact that these interactions involve three neutral bosons, two of which have unknown masses, and that the matrix R depends on the scalar potential, which is largely unknown so far as well, renders these contributions very hard to deal with phenomenologically, even when a specific model like the A2HDM is assumed. To avoid these difficulties, in the literature typically the dominance of the contribution from the lightest scalar is assumed; this is however problematic, as we will discuss below. While we will still apply this assumption occasionally to obtain indicative numbers for the neutral couplings, we can also include new information compared to earlier analyses: thanks to the huge amount of data collected recently by the LHC experiments and to lesser extend the Tevatron ones, we have some information already on the matrix R. The collider data shows that the Higgs-like state discovered at the LHC couples to W + W − and ZZ with a strength close to the SM one; assuming that it corresponds to the lightest neutral scalar h, one gets |R 11 | > 0.80 at 90% CL [27,156]. The orthogonality of R implies then |R 21 | 2 + |R 31 | 2 < 0.60 at 90% CL.
We now return to the assumption of dominating contributions from the lightest Higgs. Since EDMs are T and therefore CP violating, the contributions from neutral scalars typically involve the combinations Re y While the rotation matrix R is unknown (to the extend discussed above), we do know that it is orthogonal. This property yields one central relation for these couplings (ξ d,l = 1, ξ u = −1): This sum is therefore independent of the scalar potential and obviously vanishes for f = f , couplings ς f,f with identical phases (e.g. real couplings, as for example present in Z 2 models or MFV as defined in [45]), and also for F (f ) = F (f ) when the ς f,f are family-universal (as in the A2HDM). In the expressions below, the terms are weighted typically by some function of the neutral Higgs masses, making it most relevant for degenerate Higgs masses. However, the expression implies that all contributions stemming from CP violation in the potential involve mass differences of the neutral scalars, and that generally large cancellations can be expected in the neutral sector. The precise form of the matrix R depends on the potential for the scalar fields; note that in general the mass eigenstates do not correspond to CP eigenstates. For a CP-invariant potential, specifically, the rotation takes the simple form whereα is often denoted α − β in Z 2 models, leading to with similar expressions for the remaining combinations. Note that in this case all contributions vanish for real ς u,d , while in general mixing between the CP-odd and -even components can induce CP violation as well. The general argument above is strengthened by a second important observation, namely that the scalar mixing angles are not independent of their masses. To be specific, let us consider the limit where the second scalar doublet Φ 2 receives a very large mass and decouples from the low-energy effective theory. In this limit, the (SM-like) light Higgs has a mass M 2 h ∼ O(v 2 ), while all the other scalars become heavy and degenerate, i.e., M 2 . If the potential is CP symmetric, the mixing angle in Eq. (43) vanishes in the decoupling limit: tanα ∼ O(v 2 /M 2 H ± ). More generally, allowing for CP violation in the scalar potential, this limit yields the following form for the scalar mixing matrix: with some potential-dependent angle θ CP which vanishes if CP is conserved. 11 This implies Imy ϕ 0 1 f = 0 and, therefore, the cancellation of Eq. (42) takes place only between ϕ 0 2 and ϕ 0 3 , which in addition have equal masses in this limit. Thus, in the absence of complex Yukawa couplings, the sum of scalar contributions would vanish even with mass-dependent weight factors. This fact is sometimes overlooked in the literature, leading to claims of non-vanishing contributions in the decoupling limit which are not correct in this context.
Together, these observations imply two strong statements: 2. Even with the right-hand side present, the only contribution not suffering this suppression stems from the factors ς u,d,l which determine also the charged Higgs interactions.
In both cases, the approximation of simply taking the contribution from the lightest Higgs is not a good one; specifically, it is not conservative.

The Aligned Two-Higgs-Doublet Model
We are now prepared to proceed and give the expressions for the relevant coefficients within specific models. We do this exemplarily for the A2HDM [15,17]. We discuss here only the constraints from EDMs; for other phenomenological constraints, see [24][25][26][27].
In the A2HDM, the problem of FCNCs is circumvented by assuming at some scale Λ A alignment of the two Yukawa matrices present for each fermion species. The flavour-changing Higgs couplings are then determined by the CKM matrix and the three complex parameters ς u,d,l mentioned above, constituting new sources for CP violation. The various models with Z 2 symmetry appear as limiting (CP-conserving) cases of these parameters, see [15]. While renormalization induces some misalignment, the structure of the model prevents these effects from becoming sizable [15,17,157,158]. Note, however, that the operators additionally generated by the misalignment are not relevant in this context, since they are not flavour-diagonal.
The resulting Yukawa couplings have the form given in Eqs. (38) and (39), where now ς u,d,l are complex numbers of O(1) instead of matrices. Specifically, as mentioned above, the right-hand side of Eq. (42) vanishes in this case for two fermions of the same electric charge.
We now turn to calculating the expressions for the different classes of diagrams in the A2HDM, contributing to the effective coefficients in Eq. (1). The phenomenological analysis of these expressions is postponed to the next section.

Four-fermion operators
For the four-fermion interactions, cf. Fig. 2(c), we obtain for the A2HDM Because the neutral Higgs coupling y ϕ 0 i f is identical for fermions of the same charge, the ratio C f f /(m f m f ) is rendered family-independent. As noted above, the electron-quark couplings are important for the EDMs of atoms and molecules. An estimate of the contributions for the neutron from four-quark operators on the other hand reads where we included the RGE enhancement by an approximate factor of five, cf. [110]. Note that additionally the cancellation discussed above has to be considered. This implies at most a moderate contribution, which is well below the two-loop contributions discussed later. Therefore, we neglect it in the following. An analogous statement holds for mercury.

The Weinberg operator
As mentioned before, the Weinberg operator is of special importance, as its contribution is neither suppressed by light quark masses nor by small CKM elements. Here we have calculated the different contributions in the A2HDM; our results agree with the results obtained in [8,159] when translating them into the language of complex propagators used there.

Charged Higgs contribution
As described in Sec. 2.3, we perform the analysis of the charged Higgs contribution in an effective field theory framework [112,114], which simplifies the problem to the calculation of one-loop diagrams. The corresponding amputated diagram contains contributions from two operators; the correct coefficient can be read off from the Dirac structure γ µ / qγ 5 , for which the additional contribution is absent [160]. Our result reads where x tH = m 2 t /M 2 H ± , which is to be used in Eq. (37) to obtain its contribution to C W . We have checked that this result agrees with [114], noting that their Lagrangian for charged Higgs exchange corresponds to ours for n = 2 scalar doublets, −Y 12 /Y 11 = ς u and Y 22 /Y 21 = −ς d . Our common result in turn corresponds at the matching scale to the one obtained in [159].

Neutral Higgs contribution
For the neutral Higgs contribution the full two-loop diagram has to be calculated, as for a top quark in the loop internally only heavy degrees of freedom appear. The calculation proceeds via the following steps: the three-gluon matrix element is obtained by using from every field-strength tensor only the part containing derivatives, and summing over all possible permutations, leading to ( Here we ignored higher orders in p 2 i /M 2 (M ∈ {M H , m t }) and used µ a (p 1 ) b µ (p 2 ) = µ b (p 1 ) a µ (p 2 ) as well as p µ µ a (p) = 0. The other side of the matching condition is calculated by again summing over the different momentum configurations for the two-loop diagram, identifying the part proportional to the same Dirac structure in the corresponding expression, expanding carefully in the external momenta, and using the Feynman parametrization for the remaining integrals. The different integrals combine to give the result 12 which is again in agreement with [8] for the top-quark contribution, where however the b-quark one was neglected. Here, h(m, M ) is defined by 13 implying a much weaker suppression of the corresponding contribution, which might be compensated if |ς d | |ς u |.

Barr-Zee diagrams
The diagrams for EDMs introduced by Barr and Zee [149] (and later generalized for the gluonic dipole moment [150,151], see furthermore [161][162][163][164][165]) are proportional to the light quark mass, which at first sight leads to the conclusion that they should be tiny compared to the contribution from the Weinberg operator, which does not suffer this suppression. However, the following arguments show that their contributions are in fact comparable (cf. also [150]): • e 3 and e g 3 s are of similar size at µ tH .
• The anomalous dimension for the Weinberg operator is larger, implying a stronger suppression from the running.
• The parametric integral of the Weinberg operator is smaller. 12 Note that in principle the correct procedure for the b-quark contribution would be analogous as for the charged Higgs contribution, i.e. integrating out the Higgs, running the resulting 4-quark operator down to µ ∼ m b and matching it on the Weinberg operator. This produces a potential enhancement from a smaller anomalous dimension. However, considering the enhancement for the charged Higgs, the resulting contribution would be at most on the level of the one from Barr-Zee diagrams discussed below. As their relative sign is unknown, it would therefore not improve the limit on Im(y 2 d ) given later, which is why we use this simplified treatment. 13 The inner integral can be done analytically, simplifying the numerical analysis. Note the factor of 2 between the definition of h(m, M ) in [8] and [159], the latter of which we are using here. discussed here. As in [149] we assume that they do not exhibit strong cancellations with the ones included in our calculation. 15 Furthermore we observe that in general the CEDM contribution dominates over the EDM one. Therefore our calculation is still expected to give reasonable upper limits on the CP-violating parameters in the Yukawa sector, to the extent discussed further in Secs. 4.3 and 5.2.

Charged Higgs contribution
As mentioned before, for CP-violating charged Higgs couplings there exist a number of corresponding diagrams not calculated in [149]. For the electron EDM, the contribution with CP violation stemming from the Yukawa couplings has been calculated in [152]. In general, the translation to quarks is non-trivial, as the authors give the result only for m ν (= m b ) = 0, while for the quark EDM e.g. contributions with relative weight m b ς d /(m t ς u ) could exist. However, our analysis yields that all additional contributions are either of the order , or CP-conserving, implying that the translation is in this case possible without evaluating new diagrams. We therefore start from the result from [152], 16 identifying c t = ς u and c e = −ς l , and use 17

Charged Higgs contributions
The contributions from charged-Higgs exchange are relevant for the neutron and electron EDMs, only; in diamagnetic systems, they are usually negligible, since they contribute neither to CEDMs nor electronnucleon interactions. They vanish whenever the relevant factors ς u,d,l lack a phase difference, similarly to the SM contributions. We start by analyzing the constraint from the electron EDM as obtained in Sec. 2.2. The charged Higgs contributes via Barr-Zee diagrams, cf. Eq. (59); the resulting constraint is shown in Fig. 3 on the left, and implies |Im(ς u ς * l )| 0.02−0.34 (ς u,33 ς * l,11 ), depending on the charged scalar mass and the choice for |d e | in Table 1, together with |Im(ς u ς * l )|/M 2 H ± ≤ 10 −5 GeV −2 , to be compared with |ς u ς * l |/M 2 H ± ≤ 10 −2 GeV −2 obtained in [17]. This demonstrates already the strength of EDMs in constraining CP-violating parameter combinations. The main contribution to the neutron EDM stems from the Weinberg operator, especially since there are no sizable contributions to the CEDM involving the charged Higgs. For the considered range of charged Higgs masses, the relative contribution from the corresponding Barr-Zee diagrams, cf. Eq. (61), is about 15% of the one from the Weinberg operator.
Using Eqs. (37) and (48), we plot the resulting constraint in the Im(ς * u ς d ) − M H ± plane (ς * u,33 ς d,33 ) in Fig. 4. For a charged-Higgs mass of ∼ 500 GeV, Im(ς * u ς d ) 1 remains allowed, which is strengthened to ∼ 0.3 for light masses. We emphasize that therefore no fine-tuning is necessary to avoid this bound; however, the next-generation experiments will put this scenario to a non-trivial test, i.e. we would generally expect contributions within the projected sensitivity.
To illustrate the impact of this bound in the A2HDM, we show on the right-hand side the comparison to the one arising from the branching ratio for b → sγ [24] in the complex ς u ς * d plane, an observable known for its high sensitivity to a second Higgs doublet. While an imaginary part of O(1) is still possible, it follows from the discussion in [24] that large effects in other observables like A CP (b → sγ) are excluded by this constraint.

Neutral Higgs contributions
As discussed in Sec. 4.3, the contributions from neutral scalars are more involved phenomenologically. Specifically, cancellations are likely to play an important role, cf. Eq. (42). Since these cancellations take place for both limiting cases, universal Higgs masses and decoupling, and furthermore the mixing into the lightest mass eigenstate is rather small, see the discussion in Sec. 4.3 and in [27,156], we use the right-hand side of Eq. (42) as an approximation of the appearing sums. However, since the two limits imply different patterns for the single contributions, we evaluate the mass-dependent functions for a varying effective mass M ϕ , allowing the corresponding coefficient to take any value between the two The constraints shown can be translated back into the corresponding parameter combinations whenever a specific model with known Higgs masses is discussed. For neutral Higgs exchanges between different families, the resulting constraints allow for a comparison with the charged Higgs contributions, albeit with some caution. Note that the two contributing terms, Re(y , but for f, f = d, l with opposite signs, implying further cancellations, since the coefficient functions given in the previous section have identical signs. If one of the involved fermions is an up-type quark, the two contributions instead strengthen the bound.
For the cases in which the right-hand side vanishes, we provide the value of the contribution from the lightest scalar as a reference, which is not to be understood as a conservative limit of any kind, but as the strongest obtainable limit for the corresponding couplings in a specific model.
For neutral scalars, we do not include the contributions from the Weinberg operator, for two reasons: first they are subject to the cancellations discussed above (even in the most general case), second they are slightly smaller than the contributions from Barr-Zee diagrams with the same coefficients (for universal ς u,d ), which enter now via the chromomagnetic moments.
We start again with the constraints from the electron EDM. In Fig. 3 on the right, the constraint for |Im(ς l ς * u )| (ς l,11 ς * u,33 ) is displayed, plotted against M ϕ . This parameter combination is now bound to be 0.01 − 0.2, depending on M H ± and the choice for |d e |; this is about a factor 2 stronger than the charged-Higgs constraint for M H ± ∼ M ϕ , as can also be deduced directly from Eqs. (58) and (59). While this again does not yet call for severe fine-tuning of the parameters at the moment, the bounds are strong already, especially when accepting the bounds from ThO with restricted fine-tuning. Clearly, the coming experiments, see Table 3, will explore a region of parameter space in which we would generally expect a signal. The contribution with a tau lepton in the loop, proportional to Im(y 2 l ) (Re(y l,11 )Im(y l,33 ) and Re(y l,33 )Im(y l,11 )), is subject to strong cancellations in the A2HDM; there is therefore no conservative limit. The contribution at the lightest Higgs mass yields |Im(y 2 l )|/2 ≤ 2 − 15, depending on |d e |. For a beauty quark in the loop, the constraint is weaker than the one obtained from the bound on the electronnucleon couplingC S ; it is therefore omitted. It is noteworthy, however, that the single contributions are smaller than the ones with the tau in the loop, despite the larger mass of the beauty quark, due to the smaller charge and the occurring cancellation. Finally, the gauge boson loops give potentially large contributions, however again subject to strong cancellations. Furthermore, since the admixture of the lightest mass eigenstate with the second doublet is small, this contribution gets further suppressed. Having this in mind, however, the contribution from the lightest neutral scalar yields |R 11 Im(y ϕ 0 1 l )| ≤ 0.01 − 0.07, again depending on the value for |d e |.
The main constraint from the neutron EDM is again for the combination ς * u ς d (ς u,33 ς d, 11 ), since it is enhanced by the top mass and involves different families. The resulting constraint, shown in Fig. 5, is similar to the one from the charged Higgs exchange via the Weinberg operator; given this situation, the treatment for the hadronic matrix element is decisive for their relative strength and possible cancellations. The other constraints are either again subject to strong cancellations (|Im(y 2 u )|/2 ≤ 1.4, |Im(y 2 d )|/2 ≤ 26 and |R 11 Im(y φ 0 1 d )| ≤ 3.6 for the lightest scalar) or not constraining due to the small masses involved. The final constraints we consider stem from the mercury EDM. As discussed above, relating this observable to fundamental parameters is complicated by large theory uncertainties. However, e.g. the electron-nucleon couplings are not as strongly affected by these uncertainties, providing a more reliable bound. Furthermore, it is a conservative one: this contribution is not expected to be dominating this observable; for that reason, assuming this contribution to saturate the experimental limit is conservative. This fact has been used in [23] to obtain the limitC S ≤ 7×10 −8 (mainly) from the mercury measurement, thereby allowing for a model-independent limit on the electron EDM. Here, as we are expressing both C S andC P by coefficients of four-fermion operators, we make this assumption for their combination appearing in Eq. (31); additionally, we show the bounds from ThO with the fine-tuning assumption. The resulting constraints are shown in Fig. 6 on the left. They do not appear very strong numerically, but constrain a parameter combination which was allowed to be very large before and are therefore relevant. Note that the contribution fromC P weakens slightly the constraint compared to usingC S ≤ 7 × 10 −8 , but not severely.  Figure 6: Actual constraint from the mercury EDM (95% CL) in the |Im(ς d ς * l )| − M ϕ plane (left) for the old bound on |C S | as well as the ones from ThO, and the potential constraint from the same system in the |Im(ς d ς * u )| − M ϕ plane (right), see text.
Further contributions enter via the (Barr-Zee-)CEDM contributions to the Schiff moment, yielding potentially strong bounds; we illustrate their potential impact by using simply central values for the hadronic parameters in the equations above to obtain a constraint. An exemplary result, corresponding to a more reliable theoretical situation, is shown in Fig. 6 on the right. We note that for the assumed situation, it would be the strongest limit on |Im(ς d ς * u )| available. Furthermore, the physical mechanisms are different for the various systems. Specifically, for mercury the charged Higgs plays a minor role, so a possible cancellation for the neutron between these two contributions cannot take place here. Theoretical progress for this observable would therefore be very valuable.

Conclusions
EDMs are very sensitive probes of NP models incorporating additional sources of CP violation. In particular, they strongly constrain possible new flavour-blind phases, as those present in generic 2HDMs without tree-level FCNCs. We have critically analyzed the present experimental limits on EDMs of elementary particles and composite systems (nucleons, nuclei, atoms and molecules), and have derived the resulting phenomenological constraints on the 2HDM parameters.
To be specific, our final results are written in the context of the A2HDM, where the alignment in flavour space of the two Yukawa matrices coupling to a given right-handed fermion guarantees the absence of tree-level FCNCs, while allowing for flavour-blind Yukawa phases. This theoretical framework includes (and generalizes) all particular (CP-conserving) types of 2HDMs based on discrete Z 2 symmetries, usually adopted in the literature. Nevertheless, our findings can be directly applied to even more general Yukawa structures with simple notational changes.
The symmetries of the A2HDM protect in a very efficient way the flavour-blind phases from undesirable phenomenological consequences. Although the present experimental limits impose indeed strong bounds on the CP-violating parameter combinations, O(1) contributions remain allowed. However, large enhancements in other CP-violating observables are already strongly restricted by the present EDM bounds.
A strong caveat to keep in mind is the strong sensitivity of the EDM predictions to the UV completion of the low-energy 2HDM. Since the A2HDM flavour symmetries strongly suppress any possible treelevel or one-loop contribution, the predicted EDMs originate from two-loop diagrams. Therefore, these theoretical results could easily be changed by NP contributions beyond the 2HDM, as happens for instance in supersymmetry, and unexpected cancellations could also take place. The EDM constraints should then be interpreted with a lot of care.
Within the A2HDM, the dominant mechanisms generating non-zero EDMs are charged and neutral scalar exchanges through two-loop diagrams of the Weinberg and Barr-Zee type. While the charged Higgs contributions can be determined unambiguously, the mixing among the three neutral scalars makes their effect much more subtle. We have shown that the neutral scalar contribution for a given fermion species vanishes exactly in two opposite limits: universal Higgs masses and decoupling. The null result is due to the orthogonality of the scalar mixing matrix, which generates exact cancellations among the contributions from the three neutral scalars. This fact has been sometimes overlooked in the literature, leading to claims of non-vanishing contributions in the decoupling limit which are not correct in this context. In particular, simply taking the contribution from the lightest Higgs is not necessarily a good approximation. In order to obtain a phenomenological estimate of the neutral scalar effect we have taken an average scalar mass to evaluate any mass-dependent function and followed the prescription indicated in Eq. (62); we have only provided as a reference the value of the contribution from the lightest scalar in those cases where the right-hand side of Eq. (62) vanishes.
Our final phenomenological results are shown in Figs. 2 to 5. In spite of all previous comments of caution, these plots indicate that interesting signals could be expected within the projected sensitivity of the next-generation of EDM experiments. Experimental progress in this field could then bring a break-through in the search for NP phenomena.