Direct and indirect constraints on CP-violating Higgs-quark and Higgs-gluon interactions

We investigate direct and indirect constraints on the complete set of anomalous CP-violating Higgs couplings to quarks and gluons originating from dimension-6 operators, by studying their signatures at the LHC and in electric dipole moments (EDMs). We show that existing uncertainties in hadronic and nuclear matrix elements have a significant impact on the interpretation of EDM experiments, and we quantify the improvements needed to fully exploit the power of EDM searches. Currently, the best bounds on the anomalous CP-violating Higgs interactions come from a combination of EDM measurements and the data from LHC Run 1. We argue that Higgs production cross section and branching ratios measurements at the LHC Run 2 will not improve the constraints significantly. On the other hand, the bounds on the couplings scale roughly linearly with EDM limits, so that future theoretical and experimental EDM developments can have a major impact in pinning down interactions of the Higgs.


Introduction
The discovery of a 125 GeV boson at the Large Hadron Collider (LHC) is a breakthrough towards a deeper understanding of the mechanism of electroweak symmetry breaking [1,2]. Current data are consistent with the spin-parity assignment J P = 0 + and indicate that the couplings of this boson to the gauge vector bosons (γ, g, W , Z) and the third family of fermions (t, b, τ ) are consistent with those of the standard model (SM) Higgs boson [3].
The current level of accuracy, however, leaves room for possible deviations from the SM picture. In fact, the Higgs couplings to gauge bosons and t, b quarks are known with an uncertainty of O(20 − 30%) [3], while the couplings of the Higgs to first and second generation fermions are much less constrained [4,5]. Clearly, better knowledge of the Higgs couplings will shed light on the nature of EWSB mechanism and will also have non-trivial implications for other aspects of Higgs phenomenology (such as Higgs portal Dark Matter [6]). Improving the sensitivity and constraints to the Higgs couplings is a major goal of Run 2 at the LHC and is becoming an increasingly important target for low-energy indirect probes.
The analysis of non-standard Higgs couplings can be conveniently performed within an effective field theory (EFT) framework. There are at least two scenarios that can be used to describe current data: (i) Linear realization, in which the observed Higgs forms an electroweak (EW) doublet with the would-be Goldstone modes associated with spontaneous breaking of the EW group (that manifest themselves as longitudinal degrees of freedom of the massive gauge bosons W ± and Z). In this framework the leading dimension-6 operators describing new physics and in particular new Higgs and EW dynamics have been classified in Refs. [7,8]. (ii) The other option is that the boson discovered at the LHC is actually a light composite state associated to new strong dynamics. Explicit models of composite Higgs have been put forward starting with the pioneering work of Refs. [9,10]. This class of models can be best analyzed within the framework of the electroweak chiral Lagrangian with a light singlet Higgs state [11][12][13][14].
In both scenarios (i) and (ii) outlined above, there already exist EFT analyses of non-standard Higgs couplings. Most global analyses (see [3] and references therein) make assumptions about the flavor and CP structure of the Higgs couplings, such as minimal flavor violation [15,16], that reduce the number of operators considered. While systematic studies of flavor-violating Higgs couplings exist in the literature [17,18], analyses of CP-violating (CPV) Higgs couplings have typically focused on subsets of operators [19][20][21][22][23]. Here we wish to initiate a systematic study of the flavor-diagonal CPV couplings of the Higgs, starting with its couplings to quarks and gluons and leaving the discussion of couplings to weak gauge bosons and fermions to future work. Our study is primarily motivated by the need to learn as much as possible in a model-independent way about the recently discovered Higgs, including its CP properties (for recent discussions of CP violation in the Higgs sector in the context of the Two-Higgs Doublet Model see Refs. [24][25][26][27][28]). Moreover, CPV in the Higgs sector might have implications for weak scale baryogenesis in a number of scenarios beyond the SM (BSM). And finally, we expect strong bounds on nonstandard CPV Higgs couplings from permanent electric dipole moments, somewhat in contrast to the CP-conserving couplings, which are harder to constrain.
In this work we focus on the linear EFT realization for the Higgs sector and leave the discussion of strongly interacting light Higgs to a future study. Within this setup, our analysis involves both indirect and direct constraints, along the lines described below: renormalization group evolution from the scale of new physics down to the hadronic scale, including all the relevant SM heavy particle thresholds (Section 2).
• In Section 3 we study in detail the indirect constraints coming from electric dipole moments (EDMs). All bounds are derived assuming that the Peccei-Quinn mechanism [29] is at work. We pay special attention to the role of hadronic and nuclear uncertainties. We present bounds corresponding to current and prospective experimental sensitivities and we assess the impact of improving the theoretical uncertainties on hadronic and nuclear matrix elements.
• In Section 4 we study the direct constraints from LHC Higgs production and decay as well as tt and tth production, presenting bounds from current data and prospective sensitivities at LHC Run 2. We focus here on CP-conserving observables that depend on the square of the CP-violating couplings as these observables currently give the strongest constraints.
• In our analysis we first obtain bounds on the effective couplings by "turning on" one coupling at the time at the high scale. We subsequently study the case in which two operators are switched on simultaneously (Section 5).
• In our concluding discussion (Section 6) we compare the strength of the indirect and direct bounds for the various couplings. We summarize the current status and describe the impact of prospective sensitivities in both planned EDM searches and Run 2 at the LHC.

The set of operators and its renormalization-group evolution
Our analysis assumes the existence of new physics involving heavy degrees of freedom, that modify the low-energy dynamics via a number of SU (3) C × SU (2) W × U (1) Y -invariant local operators of dimension 5 and higher [7,8]. Here we are interested in CPV operators involving the Higgs doublet, quarks, and gluons, so at some scale M / T v (v 246 GeV is the Higgs vacuum expectation value (vev)) we consider the following effective Lagrangian, where L SM denotes the SM Lagrangian. The operators in L 6 are written in terms of the Higgs doublet ϕ, the left-handed quark doublet q L , the right-handed quark singlets u R and d R , and the gluon field strength G a µν . We have introduced the notationφ = iσ 2 ϕ * and σ · G ≡ σ µν G a µν t a , and we are suppressing generation indices. The 3 × 3 matrices Y u,d andΓ u,d induce anomalous Yukawa interactions and quark color dipole moments, respectively. The θ term represents a CPV interaction between the Higgs field and gluons. Note that the couplings Y u,d and θ have massdimension −2, whileΓ u,d has mass dimension −1, due to the explicit factor of 1/v associated withΓ u,d . Finally, ε µναβ denotes the completely antisymmetric tensor with ε 0123 = +1.
In the unitary gauge, we can write the Higgs doublet as ϕ = (0, v + h) T / √ 2. To O(h 0 ) the couplings Y u,d then contribute to the quark mass matrices, while θ produces a shift in the SM QCDθ term. 1 The remaining O(h) terms in Eq. (1) give rise to effects that are not described by the SM, and in particular induce anomalousqqh and CPV Higgs-gluon interactions.
Working in the basis in which the full quark mass matrices (including SM and BSM effect from Y u,d ) are diagonal, the operators of Eq. (1) in combination with the SM Yukawa interactions (L Y ), give the following contributions to L eff , 2 where we use the compact matrix notation Re A ≡ 1/2(A + A † ) and Im A ≡ 1/(2i)(A − A † ). m u,d are the real and diagonal quark mass matrices, while Re Y u,d , Im Y u,d , andΓ u,d are not necessarily diagonal. Studies of the flavor-violating couplings induced by Y u,d have appeared in the literature [17,18], while the CPV third generation couplings have been studied in [19,20]. As for the gluon dipole operators, the EDM constraints on light quark diagonal couplings have been studied in Ref. [21,30], and the top chromo-EDM (CEDM) has been studied in several papers, see for example Refs. [22,23,31]. In this work we focus on the CPV flavor-diagonal couplings arising from Eq. (2) and we ignore the real and the flavor-violating parts of Y u,d . In the dipole operator sector, we focus on the top CEDM as it strongly mixes with Im Y t and θ (in Appendix A we summarize the bounds on the light quark CEDMsd q , q = t). That is, we take at the high scalẽ Γ d = 0 andΓ ij u = δ i3 δ j3Γt , and extend the results of Ref. [22] by taking into account hadronic uncertainties, including additional mixing effects, and considering additional collider constraints from Higgs production.
In summary, our starting point high-scale CPV and flavor-diagonal effective Lagrangian reads: (3) where nowd t ≡ ImΓ t /g s , Im Y q denotes the diagonal entries of the matrices in Eq. (2), and the dots stand for interactions involving two or more Higgs fields.

Renormalization-group equations
In order to connect the above operators to measurements taking place at energies below M / T , the renormalization-group equations (RGEs) governing the scale dependence of these operators are 1 We assume in this work that the Peccei-Quinn mechanism [29] is at work, so thatθ (including the shift δθ = (1/2)v 2 θ ) relaxes toθ ind = 0 due to the distortion of the axion potential induced by the higher dimensional operators.
2 Denoting the Standard Model Yukawa couplings by LY = − √ 2qLY d dRϕ − √ 2qLYuuRφ, the quark mass matrices are given by  required. Relating these interactions to EDM experiments necessitates evolving them down to the QCD scale, Λ χ 1 GeV, below which QCD becomes strongly coupled and non-perturbative techniques are required, see Sect. 3.2. As the operators in Eq. (3) do not contribute to EDMs directly, but only through their mixing contributions to other operators, we require an extended basis of operators that includes the light fermion (C)EDMs and the Weinberg operator. Accordingly, we extend the effective Lagrangian as follows: Taking the basis, C q = (d q /eQ q m q ,d q /m q , d W /g s , Im Y q , θ ) T , the one-loop QCD RGEs can be written as [32][33][34][35][36][37][38], where , and N (n f ) is the number of colors (flavors). As we are only interested in the operators of Eq. 3 Note that the electron EDM introduced in Eq. (4), d e , does not appear in the RGEs since it is not affected by one-loop QCD renormalization. However, it is generated by threshold corrections that are discussed below.

Evolution to µ = m t
When considering the contribution of the dimension-6 operators in Eq. (3) to collider observables, it is mainly the mixing among the operators,d t , Y q , and θ , themselves that is relevant. The above RGEs can be used to first run the couplings down to µ = m t , where the top CEDM and Yukawa coupling are integrated out. At this scale, the top Yukawa induces a contribution to θ through a top loop, We present the numerical results of this procedure in Table 1, where we employ the following values [3], All quark masses are given in the MS scheme. A fixed-order perturbation-theory solution of Eq. First the RGEs of Eq. (5) are used to run the operators from µ = M / T to µ = m t , where we integrate out the top quark and the Higgs boson. This implies that the couplings θ , Y t , andd t and their corresponding operators are removed from the EFT below µ = m t . Eliminating these operators gives rise to several threshold corrections to the operators in Eq. (4). The Yukawa interactions contribute to the (C)EDMs [39][40][41][42] and the Weinberg operator [32,43] through Barr-Zee diagrams, shown in Fig. 1. The quark (C)EDMs receive additional contributions from the one-loop diagrams shown in Fig. 2 [44]. The top CEDM gives rise to a one-loop threshold contribution to the Weinberg operator [34,45]. In total we have the following matching conditions, where , and the functions f , g, and h are given by, The first loop terms, contributing to the quark (C)EDMs in Eq. (8), are due to the Barr-Zee diagrams involving quark loops. 5 The second term, contributing to quark EDMs, originates in Barr-Zee diagrams involving an internal W ± loop. The remaining terms are from the one-loop graphs in Fig. 2. The loop terms contributing to d W arise from, respectively, the top qCEDM threshold correction and the fourth Barr-Zee diagram in Fig. 1. The contribution ofd t to the Weinberg operator was also considered in Ref. [22], however, due to its mixing with θ , we obtain somewhat larger contributions ofd t to the operators at 1 GeV (shown in Table 2).
Below µ = m t , our basis consists of the operators explicitly listed in Eq. (4). The RGEs can then be used to run down to µ = m b and subsequently to µ = m c . At these thresholds the After the charm threshold the remaining operators can be evolved to Λ χ using Eq. (5). The numerical result of this analysis is presented in Table 2 for M / T = 1 TeV. A fixed-order perturbation approximation, as the one mentioned in Sec. 2.1.1, is less accurate below m t as α s runs faster in this regime. However, the solution to the RGE below m t is simpler, as only the quark (C)EDMs and the Weinberg operator are involved, and is explicitly given in Ref. [35].

Constraints from electric dipole moments
Strong constraints on the CPV higher-dimensional operators can be derived from experimental upper bounds on EDMs. The strongest constraints arise from measurements on the neutron [46], the 199 Hg atom [47], and the ThO molecule [48]. To interpret the experimental upper bounds, it is necessary to express the observables in terms of the Wilson coefficients of the dimension-6 operators and the corresponding hadronic and nuclear matrix elements. In particular, for the operators involving quarks and gluons this is problematic due to the non-perturbative nature of QCD at low energies. Nevertheless, various techniques have been applied to calculate EDMs directly in terms of CPV quark-gluon operators. Depending on the operator under investigation, the techniques vary in their sophistication and accuracy. Typically in the literature, the uncertainty of the calculations is not taken into account and only the central values of the results are considered. In this work we take into account the theoretical uncertainty and show that in some cases this drastically weakens the constraints on possible CPV from BSM physics.

Experimental status and prospects
We briefly summarize the current experimental status and the outlook on future possibilities. At the moment, the strongest constraints have been set on the EDMs of the neutron [46], d n , the 199 Hg atom [47], d Hg , the 129 Xe atom [49], d Xe , and on the energy shift indicating T violation   in the ThO molecule [48]. As discussed below, for our purposes the latter can be interpreted as a constraint on the EDM of the electron, d e . Recently, a first measurement of the EDM of the atom 225 Ra, d Ra , has been reported [50], but the experiment is not precise enough to impact the constraints discussed below. The outlook of EDM experiments is very positive. Measurements on d n and d e are expected to improve by one to two orders of magnitude, while the limits on d Xe and d Ra will be improved by several orders of magnitude. On the longer time-scale, experiments are being developed to measure the EDMs of light nuclei (proton and deuteron and perhaps helion) in electromagnetic storage rings [51,52]. These experiments have a projected sensitivity of 10 −29 e cm. In Table 3 we summarize the current limits and expected sensitivities for a variety of EDMs. The future sensitivities are meant to be only indicative at this stage (see [53,54] and references therein).

Nucleon EDMs
For the dimension-6 operators under investigation the low-energy operators relevant for hadronic and nuclear EDMs are the light quark (C)EDMs and the Weinberg operator. Theoretically, by far the best understood operators are the quark EDMs whose contributions to the nucleon EDMs have been recently calculated with lattice-QCD techniques [55,56], see Table 4. The uncertainties on the up and down qEDM contributions are 10% − 15%, whereas the strange qEDM contribution is consistent with zero and thus highly uncertain. Although the matrix element is smaller than for the up and down qEDMs, the Wilson coefficients typically scale with the quark mass which means that the largest uncertainty arises from the strange EDM contribution.
Unfortunately no lattice-QCD calculations exist for the qCEDM contributions (see Ref. [38] for preliminary steps towards such a calculation). Instead, the most-used results are obtained with QCD sum rules [57][58][59][60] which are consistent with a chiral perturbation theory (χPT) calculation combined with naive dimensional analysis (NDA) [61]. The matrix elements are shown in Table 4, where it must be stressed that these results apply only if a Peccei-Quinn (PQ) mechanism is invoked to remove the QCDθ term [62]. 6 The uncertainty is estimated to be 6 Note that in the presence of BSM sources of CP violation, the PQ mechanism relaxesθ to a finite valueθ ind significant, O(50%), for the light qCEDM contributions. More problematic is the dependence of the nucleon EDMs on the strange CEDM. Typically, in the PQ scenario, the contribution from the strange CEDM is taken to vanish [57,60]. However, a recent calculation based on SU (3) χPT found a much larger dependence [66]. Here we assume no dependence on the strange qCEDM, but stress that this issue has not been resolved. 7 The least is known about the Weinberg operator. No systematic calculation exists and we must rely on estimates. An estimate based on QCD sum rules [67] gives a somewhat smaller estimate than NDA [32]. Here we take a range, see Table 4, which covers both estimates and also vary the sign of the matrix element. In principle, the matrix elements of d n and d p have an independent sign and magnitude. However, because the coupling to photons goes via the electromagnetic quark current, assuming that the (larger) isovector component dominates we take as the benchmark case that d p has a relative sign with respect to d n , but we vary their magnitude independently. We comment later on the importance of fixing the relative sign.

EDMs of light nuclei
EDM of light nuclei receive two main contributions. The one-body component is determined by the EDM of the constituent nucleons, d n and d p . A second contribution is due to modifications to the nuclear wavefunction induced by the CPV nucleon-nucleon potential. For the operators under consideration, an analysis based on chiral EFT indicates that the CPV potential is dominated by two CPV pion-nucleon interactions 8 [70,71] where N = (p n) T is the nucleon doublet, π the pion triplet, τ the Pauli matrices, andḡ 0,1 two low-energy constants (LECs). Because the quark EDMs contain an explicit photon, their contribution toḡ 0,1 is suppressed by α em /π and can therefore be neglected. The Weinberg operator is chiral invariant and therefore its contribution toḡ 0,1 is suppressed by m 2 π /Λ 2 χ where Λ χ ∼ 1 GeV [58,68]. Nevertheless, power counting indicates that nuclear EDMs can still significantly depend onḡ 0,1 induced by the Weinberg operator [68,72], but explicit calculations show that the largest contributions arise from the constituent nucleon EDMs [73,74]. We will therefore neglectḡ 0,1 from the Weinberg operator. That leaves us with the quark CEDMs, that do induce large values ofḡ 0,1 as indicated by QCD sum rules [75] These values are consistent with an SU (2) χPT analysis [61]. So far, no EDM measurements have been performed on charged particles. As we show in this work, measurements on different systems are crucial to isolate or constrain possible new (induced theta-term), proportional to the coefficient of the new physics operator. Currently the effect ofθ ind is taken into account within the QCD sum rule approach. Progress in lattice-QCD evaluations of dn(θ) (for recent results see [63][64][65]) will also improve the contribution to the nucleon EDM proportional toθ ind . 7 If the PQ mechanism is not invoked, but the strong CP problem is solved via other ways, for example via extreme fine-tuning, then the matrix elements of the up and down CEDMs shift by O(1) factors [57], while the strange CEDM matrix elements are not expected to vanish. In this work we do not pursue this scenario, and assume that the PQ mechanism is at work. 8 For the Weinberg operator, important contributions can arise from CPV nucleon-nucleon interactions, but these vanish for dD [68,69] which is the main focus here.

Atomic screening
Best values of a 0,1 Estimated ranges of a 0,1  physics, and we therefore investigate the potential impact of these measurements. The EDM of the deuteron is given by [69,73,76] where the small uncertainties are taken from Ref. [72]. The 3 He EDM has been analyzed within the same framework [72,73] and depends onḡ 0 as well. In addition, it depends on CPV nucleon-nucleon interactions induced by the Weinberg operator. We do not consider a 3 He EDM measurement in what follows.

Atomic EDMs
Next we focus on diamagnetic atoms. Schiff's theorem [77] tells us that the EDM of a pointlike nucleus is screened by the electron-cloud, ensuring that the total atomic EDM vanishes. However, in heavy diamagnetic atoms the conditions for Schiff's theorem are violated by the finite size of the nucleus. For the operators we consider, the dominant contributions 9 to diamagnetic atomic EDMs then arises from the nuclear Schiff moment S A [78]. The atomic EDM can then be written in term of an atomic screening factor A A times S A . The latter can be expressed as a function ofḡ 0,1 and the nucleon EDMs d n,p : Whereas the atomic uncertainties are rather minor (O(20%), see Refs. [79][80][81]), the dependence of S A onḡ 0,1 is far more uncertain due to the complicated nuclear many-body problem [82,83] (for a detailed discussion, see Ref. [84]). We give the best value and range for A and a 0,1 in Table 5. In addition, the nuclear Schiff moments depend on the constituent nucleon EDMs.
As far as we are aware, this has only been calculated for d Hg [74], with the result α n = 1.9(1) and α p = 0.20 (6). We neglect possible contributions from CPV short-range nucleon-nucleon interactions but stress that this assumption is untested. Finally, we discuss the constraint on the electron EDM. For the operator set discussed in this paper, there appear no significant contributions to CPV electron-quark interactions 10 , such that paramagnetic EDMs are dominated by d e . The strongest constraint then arises from the ThO measurement which gives [48] d e ≤ 8.7 · 10 −29 e cm , at 90% confidence level (c.l.). The conversion of the ThO measurement into a bound on d e entails a theoretical uncertainty from atomic and molecular dynamics, estimated at the 15% level [85,86]. Since this is substantially below the hadronic and nuclear uncertainties, we neglect it in our analysis.

Analysis strategy: central, conservative, and minimized bounds
In most of the existing literature, when discussing EDM constraints on BSM physics, the theoretical uncertainty of the hadronic and nuclear matrix elements is not taken into account. Bounds are obtained by considering the central values given in the previous section, leading to strong constraints on many BSM models. In this work, we investigate how the constraints are softened if we do consider the range of the matrix elements. To do so, we present bounds obtained by three different choices of matrix elements: 1. Central: Here we take the central value of the hadronic and nuclear matrix elements. This is the usual method of deriving EDM constraints on BSM physics.

2.
Conservative: In this case we minimize the absolute value of each hadronic and nuclear matrix elements within their given range. For example, in case of the qCEDMs we take For ranges which include zero, such as the dependence of d n on the strange qEDM or the dependence of d Hg onḡ 1 , we set, somewhat arbitrarily, the matrix elements to one tenth of the central value. For example, d n = 0.0008 d s .

Minimized:
Here we vary the matrix elements within their allowed range assuming a flat distribution, and minimize the total χ 2 of the set of EDM experiments. This method corresponds to the Range-fit (Rfit) procedure defined in Ref. [87]. It always gives the weakest constraint of the three methods discussed in this work as it allows for cancellations between different contributions. This approach gives the most conservative (perhaps overconservative, but realistic) constraints.
For matrix elements with an uncertain sign, such as the dependence of d n and d p on d W , we calculate the bounds for all permutations of the signs and present the most conservative one.

Single coupling analysis
Following the above strategies, we present the bounds on the CPV operators in Table 6. We assume here that only a single CPV coupling is turned on at the scale M / T . From the first line of the table, it is clear that d e is mainly sensitive to the Yukawa couplings of the heavy quarks, while it does not constrain the up-and down quark-Yukawa couplings and θ at a significant level. Considering the excellent theoretical accuracy in case of d e , we always take the central value of the matrix elements and do not consider the conservative or minimized case.
In contrast, d n and d Hg obtain large contributions from θ , the light-quark Yukawa couplings, and the top CEDM. Compared to d e , these EDMs are less sensitive to Im Y t , even when using central matrix elements. With central values, d Hg gives the strongest constraints on Y u,d , but this does not take into account the considerable theoretical uncertainties. Once these are taken into account, several bounds are changed dramatically. Moving from the central to the conservative strategy, the matrix elements for the Weinberg operator decrease by a factor five, which is reflected in the constraints ond t and Y c,b which mainly induce EDMs via d W . Similarly, the constraint on Y s , which mainly induces d s , is softened by a factor ten due to the uncertainty in the nucleon matrix element. Finally, the d Hg constraints on Y u,d are severely weakened due to the uncertain status of the nuclear matrix elements connecting d Hg andḡ 1 .
Moving to the minimized case, we see that the bounds become softer. In most cases, the bounds from d n are only mildly effected. The main exceptions are the bounds on Y b and, to lesser extent,d t for which several contributions of similar size contribute to d n . These contributions can mutually cancel within the minimization strategy, such that no significant constraint on Y b remains. The dramatic change in the constraint on Y s arises because the allowed range of the matrix element connecting d s to d n includes zero. The minimizing strategy has more severe consequences for the bounds from d Hg . For all operators, the uncertainties in the matrix elements are large enough to kill the constraints. This clearly reflects the additional uncertainty due to the nuclear many-body problem. Although this might sound as an extremely conservative conclusion, we show in the next section that modest theory improvements could drastically change the impact of diamagnetic measurements.
When combining the constraints from d n , d Hg , and d e , we obtain a significant constraint in all cases, even when using the minimizing strategy. Within the context of constraining non-standard Higgs couplings, this shows once more the importance of complementary EDM probes.
Finally, we briefly discuss the dependence of the constraints on the scale of new physics, M / T . This dependence enters in two ways. First of all, the couplings scale as M −2 / T , such that the constraints on the dimensionless couplings are less stringent for higher values of M / T . Second, the value of M / T affects logarithmically the evolution to lower energies. To illustrate this effect, we show the resulting constraints assuming three values for the scale of new physics M / T = 1, 10, 100 TeV in Table 7. As might be expected, these constraints differ by

Impact of more accurate hadronic and nuclear matrix elements
It is extremely instructive to study the impact of better theoretical control on the hadronic and nuclear matrix elements appearing in the EDM expressions. There are a number of matrix elements which have the largest uncertainty: • The dependence of d n,p on d s and d W . We investigate what happens if these matrix elements were known with 50% accuracy. That is, we take d n = (0.008 ± 0.004)d s and d n = (50 ± 25) MeV d W . Similarly for d p , but with a relative sign on the d W element.
• The dependence of d n,p ond u,d has an uncertainty of 50%. We reduce this to 25%.
That is a 0 = 0.13 ± 0.065 and a 1 = 0.25 ± 0.125. In the next section, we do the same for S Xe and S Ra .
In the bottom row of Table 6, we present the bounds on the CPV operators assuming these improved matrix elements. We see that the bounds on Y u and Y d , are only slightly improved, while Y t is unaffected. The consequences for the limits on the other couplings are larger, with improvements of a factor 3 to 25 depending on the coupling. The bound on Y s would be improved by three orders of magnitude. An important observation is that once we include the improved matrix elements, the minimized constraints come close to the central values constraints. That is, a comparison of the dn(10 −28 ) 5.4 · 10 −9 2.8 · 10 −9 4.1 · 10 −6 1.8 · 10 −6 1.4 · 10 −4 3.2 · 10 −5 1.9 · 10 −4 3.7 · 10 −6 dp(10 −29 ) 2.2 · 10 −10 5.5 · 10 −10 4.1 · 10 −7 rows "Comb. Cen." and "Future Min." tells us that almost all constraints only differ by a factor 2. The exceptions are the bounds on Y u,d which differ by a factor of 5 and 3. This indicates that once the hadronic/nuclear theory is at this level of precision, there is very little room for mutual cancellations between contributions. At this point, we can exploit the full power of the impressive experimental constraints on EDMs.

Impact of improved experimental bounds and additional probes
We now study how the constraints would change with improved measurements of d e , d n , and d Hg or with measurements of d Xe and d Ra that are currently not competitive. The impact of improved limits on d n , d Hg , or d e , can be simply obtained by rescaling the bounds in Table 6. In Table 8 we show constraints using expected future experimental sensitivities and central values of the matrix elements. For all couplings a measurement of d n at 10 −28 e cm would be more constraining than d e at 5 · 10 −30 e cm. This observation, however, does not take into account the hadronic uncertainties which we discuss below. Future experiments on light-nuclear EDMs such as d p and d D can have a large impact as well. The projected d p measurement at 10 −29 e cm, would be roughly 10 times better than the d n measurement. This factor is not surprising considering that the matrix elements for d n and d p are very similar. For several couplings, a measurement of d D at 10 −29 e cm would even be more constraining than a d p measurement with the same accuracy. In particular, couplings such as Y u,d that induce light-quark CEDMs give relatively large contributions to d D . This behavior illustrates the complementarity of a d D measurement [58,68]. On the other hand, d D is less sensitive than d p to couplings such as Y c,b andd t which induce relatively large contributions to d W . The central values of the matrix elements linking d W to d n and d p have opposite signs, and therefore the sum of nucleon EDMs, d n + d p , that enters in d D , see Eq. (13), vanishes. This conclusion strongly depends on the relative sign of the nucleon matrix elements which is highly uncertain.
In Table 9 we perform the same analysis using the minimization strategy. In the first two rows, we repeat the combined constraints with current and future matrix elements. In the next two rows, we include the expected increase in sensitivity of future d n and d e experiments. In rows 5 and 6 we add the prospected d Xe and d Ra measurements with and without improved matrix elements. From row 5 we see that d Xe and d Ra mainly improve the constraints on Im Y u,d , but at most a factor 4. From the comparison between rows 2 and 5, we conclude that theory improvements can have as much impact as additional experimental probes. Once improved matrix elements are added, measurements of d Xe and d Ra improve the constraints on Im Y u,d by roughly an order of magnitude and on θ by a factor of 5 over current constraints with improved matrix elements (i.e. row 2 of the table). This observation reflects that improvements in experiments and theory must go hand in hand.
In the last two rows we study the impact of d p and d D measurements. Due to the very high accuracy (10 −29 e cm) the bounds are strongly improved over the current constraints, but it must be said that these experiments have a longer time-scale. Interestingly, the bounds on Im Y s,b are dramatically improved if more accurate matrix elements are used, once more underscoring the strong impact of hadronic and nuclear theory.

Constraints from colliders
In this section we discuss the constraints that collider observables impose on the couplings θ , the pseudoscalar Yukawa couplings of the light quarks, Im Y u , . . ., Im Y b , and the top pseudoscalar Yukawa and CEDM, Im Y t andd t . We focus on total production cross sections and on branching ratios, that are sensitive to the square of the coefficients of the CPV operators. Additional information could be obtained by a study of observables that depend linearly on the CPV coefficients.
The operators we study are all constrained by the Higgs signal strengths, which are observed to be compatible with the SM [88,89]. For a given Higgs production mechanism, i → h, followed by the decay of the Higgs to the final state f , the signal strength in the presence of the dimension-6 operator O is defined as where σ SM and σ O are, respectively, the production cross sections in the SM and the correction induced by O. Γ SM,O h→f are the decay widths in the channel f and Γ SM,O tot the Higgs total width. In Sections 4.1 and 4.2 we discuss how the operators we consider affect the production and decay signal strength, µ i and µ f , and extract bounds from the LHC Run 1 [88,89].
In addition, we discuss the bounds on the top CEDM from the tt cross section, and bounds ond t and Im Y t from the tth cross section.

Limits on θ and Im Y q from Higgs production and decay
The most important manifestations of the operator θ and of the non-standard Yukawa couplings of the light quarks at colliders are the modification of the Higgs production cross section and decay width. In the SM, the dominant mechanism of Higgs production is gluon fusion through a top loop. The gluon fusion cross section has been computed at N 2 LO in α s [120][121][122], with the inclusion of top quark mass effects and of electroweak corrections. Recently, the inclusion of N 3 LO corrections has been completed [123]. Here, for consistency with the calculation of the production cross section induced by θ and Im Y q , we use the N 2 LO expression, and, at this order, for a Higgs mass of m h = 125 GeV, the SM gluon fusion cross section is [124] σ SM ggF = 19.2 +1.4 The first uncertainty takes into account the effects of missing terms in the perturbative expansion of the cross section, and it is obtained by varying the renormalization and factorization scales. The second uncertainty is the combination of the uncertainties due to the parton distribution functions (PDFs) and α s . Here and below, we always quote PDF and α s uncertainties at the 90% c.l. Higgs production through the SM light-quark Yukawa couplings is negligible, with the exception, to some extent, of the b quark [125]. Since the b quark contribution to the Higgs production cross section is only a few percent of Eq. (18) [125,126], we neglect it in our discussion.
The tree-level contributions of θ and the pseudoscalar Yukawa couplings to the Higgs production cross section are shown in Fig. 3(a,b). The cross section induced by θ and the pseudoscalar Yukawa couplings has been computed at N 2 LO in Refs. [127,128] and [125], respectively. The calculation of Higgs production through the coupling to GG was performed in the framework of supersymmetric models, where a neutral pseudoscalar Higgs boson A couples to the top quark with a pseudoscalar Yukawa-type coupling. The coupling AGG is then induced by integrating out the top quark, and, in the MSSM, its coefficient is v 2 θ = cot β. The calculation of (a) Refs. [127,128] can be simply adapted to the CPV coupling of the scalar Higgs boson to GG, by not fixing the coupling to cot β, and by neglecting higher-order corrections to the matching coefficients of the effective operators that are specific to the model of Refs. [127,128]. The Higgs production cross section induced by scalar and pseudoscalar Yukawa couplings was computed at N 2 LO in Ref. [125]. While Ref. [125] focuses on the b quark Yukawa coupling, the results can be used for any massless quark. At 8 TeV, the gluon fusion cross section induced by θ is where we neglected electroweak corrections. The first error in Eq. (19) is the scale uncertainty, obtained by varying the factorization scale µ between m h /2 and 2m h . The second error is the combination of the PDF and α s uncertainties. To estimate the central value and the PDF error, we followed the recipe of the PDF4LHC working group [129]. We evaluated the cross section using three N 2 LO PDF sets, CT10 [130], MSTW08 [131] and NNPDF2.3 [132], following the prescriptions of each collaboration to extract the 90% c.l. PDF and α s uncertainties. As our central value and PDF error we quote the midpoint and the width of the envelope provided by the central values and PDF and α s errors obtained with these three different PDF sets. Eq. (19) shows that scale variations and PDF errors have approximately the same importance. The coupling θ in Eq. (19) is evaluated at the renormalization scale µ, which we set to µ = 125 GeV. The RGE of θ is discussed in Eq. (5). If only θ is turned on at the scale M / T , to a very good approximation we can neglect the running of θ and interpret the coupling in Eq. (19) as the coupling at the scale M / T . We summarize the Higgs production cross section induced by Im Y q in Table 10. The couplings Im Y q are scale dependent. Since the calculation of the cross section neglects all mass effects, we also neglected the mixing of the light quark Yukawas to θ and to the light quark CEDM. In this approximation, the RGE of Im Y q is diagonal, and, by an appropriate choice of scheme, can be made identical to that of the quark masses [133]. For consistency with the calculation of Ref. [125], we used the three-loop anomalous dimension to run the couplings from the reference scale µ 0 = 1 TeV to the renormalization scale µ = m h . For the u, d, and s quark, once N 2 LO corrections are included, the uncertainty is dominated by PDF errors. For the c and b quarks, PDF errors and scale variations are comparable. Notice that the PDF error is particularly large for Im Y s , reflecting some issues in the determination of the strange quark PDF.
We can use the cross sections in Eq. (19) and Table 10 to construct the production signal strength µ ggF appropriate for our scenario, given by the ratio of the single Higgs production cross section in the presence of the CPV operators in Eq. (3) and the SM cross section. For θ , the signal strength is where we neglected electroweak corrections to the SM and θ cross sections, and worked in the m t → ∞ limit. The Higgs production cross section induced by the hGG operator is very similar to the SM cross section, that proceeds via the hGG operator [127,128]. As a consequence, the signal strength is very close to the tree level value of 9/4, and the PDF, α s and scale errors cancel almost completely in the ratio, leaving a negligible error on µ ggF . The production signal strengths for the pseudoscalar Yukawa couplings are summarized in Table 11. In this case, the scale variations are almost completely determined by the error on the SM cross section, while PDF errors of the same size appear in both SM and ImY q cross sections.   TeV, with theoretical uncertainties.   The operator θ does not significantly affect the decay channels that are relevant at the LHC, γγ, ZZ * , W W * and bb. It contributes, however, to the Higgs decay into gluons, thus affecting the total width. The contribution of θ to the width can be extracted from the review [136], to which we refer for references to the original calculations. In Ref. [136] the decay width of a pseudoscalar Higgs boson A into gluons via the AGG operator was considered. As for the production cross section, the contribution of the CPV coupling hGG to the width can be obtained by replacing cot β with v 2 θ , and sending the top mass to infinity. We find The contribution of θ to the width is thus less than 10% of the contribution to the production cross section. We can approximate µ θ with µ θ ggF and use the production signal strengths extracted by the ATLAS and CMS collaborations [88,89], to derive bounds. Here, and in the rest of the section, when citing experimental results we quote the uncertainty reported in the original publication, usually at the 68% c.l., and we rescale it to ATLAS 1.2 · 10 −2 1.0 · 10 −2 0.6 · 10 −2 0.7 · 10 −2 0.5 · 10 −2 CMS 1.4 · 10 −2 1.3 · 10 −2 1.0 · 10 −2 1.1 · 10 −2 1.0 · 10 −2 future 0.6 · 10 −2 0.5 · 10 −2 0.4 · 10 −2 0.4 · 10 −2 0.4 · 10 −2 The corrections to the width are more important in the case of the pseudoscalar Yukawa couplings. For the u, d, s and c quark, the pseudoscalar Yukawa contribute mainly to the total width, while Im Y b also affects the bb decay channel. Im Y b and Im Y c contribute to the γγ width, but the effect is negligible with respect to the correction to the total width. The decay width is related to the decay width of a pseudoscalar A boson into quarks, photons, and gluons, with appropriate replacement of the couplings. Using the expressions of Ref. [136], we obtain We performed a fit to the signal strengths in the various production and decay channels observed by the ATLAS and CMS collaboration [88,89,[137][138][139][140][141][142][143]. In Table 12 we show the 90% bounds we obtain by turning on one CPV coupling at a time. We find that current LHC data exclude pseudoscalar Yukawa couplings greater than about 1%. For the lightest quarks, u and d, the correction to the production cross section is compensated by the dilution of the γγ, W W * , and ZZ * decay channels, resulting in weaker bounds. For heavier quarks, s, c, and b, the smaller PDFs suppress the production cross section, and the most important effect is the correction to the total width.
In Table 13 we show the dependence of the bound on the scale M / T . In the case of the pseudoscalar Yukawa couplings, the running of Im Y q is such that the bounds get slightly stronger as the new physics scale is increased. The effect is mild, about 20%-25% in going from M / T = 1 TeV to M / T = 100 TeV. For θ , the effect of the running is negligible. The pseudoscalar Yukawa couplings modify another important Higgs production mechanism namely that of associated production with a W or Z boson. However, we find that the associated production cross section is a factor of 10 3 smaller than the single Higgs production cross section, yielding significantly weaker bounds on Im Y q .
At 14 TeV, the ratio σ θ ggF /σ SM ggF remains substantially identical to Eq. (20), due to the fact that the SM cross section, induced by the hGG effective operator, and the θ cross section, induced by hGG, have the same scaling with the center-of-mass energy. In the case of the t/mt M / T = 1 TeV 1.2 · 10 −2 1.0 · 10 −2 0.6 · 10 −2 0.7 · 10 −2 15 · 10 −2 0.5 · 10 −2 0.27 5.2 · 10 −2 M / T = 10 TeV 1.1 · 10 −2 0.8 · 10 −2 0.5 · 10 −2 0.6 · 10 −2 14 · 10 −2 0.4 · 10 −2 0.27 2.8 · 10 −2 M / T = 100 TeV 1.0 · 10 −2 0.8 · 10 −2 0.5 · 10 −2 0.5 · 10 −2 13 · 10 −2 0.4 · 10 −2 0.27 2.1 · 10 −2 pseudoscalar couplings to the u, d and s quarks, the cross section grows more slowly than the gluon fusion cross section, leading to smaller corrections to the production signal strength at 14 TeV. For c and b quarks, the PDFs are obtained perturbatively, from the splitting of gluons in heavy quark pairs. Thus, the ratio of the Yukawa and gluon fusion cross section remains approximately constant at higher center-of-mass energy. The ATLAS and CMS collaborations have released projections on the fractional error on the signal strength [134,135]. With the integrated luminosity of 300 fb −1 that will be reached at the LHC Run 2, the error on µ ggF will be dominated by the theoretical uncertainties on the gluon fusion cross section. After the inclusion of the recently completed N 3 LO corrections [123], these are in turn dominated by PDF and α s uncertainties. As far as the experimental errors are concerned, ATLAS projects a reduction of the error to 6% at 300 fb −1 , and 4% at 3000 fb −1 . Assuming a central value µ ggF = 1 and a combined theoretical and experimental error of 10%, the bound on θ improves only slightly, v 2 θ < 0.21.
We obtain the projected limits on the Yukawa couplings by assuming a 10% error on the gluon fusion signal strength with Higgs decaying in γγ, ZZ * and W W * , a 20% error on the vector boson fusion signal strength with Higgs decaying in W W * , and a 30% error on the Higgs to bb signal strength [134,135]. The projected limits on the Yukawa are given in Table 12, which shows a possible improvement to the 0.5% level for all the pseudoscalar Yukawa couplings.

Limits on top CPV couplings
The top pseudoscalar Yukawa coupling and CEDM can be probed in two ways. First of all, these couplings contribute to the Higgs gluon fusion production cross section at one loop. Because the dominant Higgs production mechanism in the SM also proceeds via top loops, we can expect it to be extremely sensitive to anomalous top-Higgs and top-gluon couplings. Secondly, these couplings affect processes with top quark pairs in the final state. In particular, we focus on the tt total cross section, for which the SM prediction is known very precisely, at the N 2 LO accuracy, and on the associated production of the Higgs and a tt pair (tth production). As discussed in Sec. 2, θ receives a threshold correction from Im Y t , and mixes with the top CEDM. The relevant diagrams are shown in Fig. 3 (c)-(e). The bound on θ can therefore be used to constrain Im Y t andd t . The gluon fusion cross section induced by Im Y t is very similar to Eq. (20), with the only difference that O(α 2 s ) corrections to the matching coefficient of Im Y t to θ need to be considered, as done in Refs. [127,128]. Their effect is to shift the signal strength in Eq. (20) by 0.005, which has no consequence on the constraints. Im Y t modifies the γγ and gg decay widths. In the case of the γγ branching ratio, we find that the corrections to γγ and to the total width are very similar, and µ γγ is, accidentally, very close to one The W W * and ZZ * branching ratios are affected by the contribution to the total width. As for θ , the corrections to the decay signal strength are about 10% of the correction to production, and, with the current experimental accuracy, can be neglected. Thus we can extract a bound on Im Y t from the total gluon fusion signal strength, Eq. (22), obtaining where Im Y t and the top mass are evaluated at µ = 1 TeV. Notice that the ratio Im Y t /m t is RG invariant, since scalar and pseudoscalar currents have the same anomalous dimension up to three loops [133]. The top CEDM mixes with θ , with anomalous dimension given in Eq. (5). We solve the RGE and run θ andd t to the top threshold, where we integrate out the top quark and stop the running. In this case, the RGE is only known at LO, and O(α s ) and O(α 2 s ) corrections to the evolution and to the matching coefficient ofd t onto θ are not known. Therefore, to put bounds ond t we use the tree-level value of the signal strength. We obtain The strong limit is a consequence of the large mixing of the top CEDM and θ . In similar fashion, the gluon fusion cross section can be used to constrain the top chromo-magnetic dipole moment through its mixing onto hGG. We discuss this in more detail in the next section. Looking to the future, a 10% accuracy in the measurement of µ ggF would allow to slightly improve the constraint to 4%. In light of the strength of the constraints, it will be interesting to include higher-order corrections to the mixing ofd t and θ . In addition,d t and Im Y t contribute to processes involving the production of top quarks. The top CEDM affects the tt and tth cross sections, while Im Y t only contributes to tth.
The top CEDM contributions to the tt cross section are shown in Fig. 3 (f)-(i). The cross section induced by the top chromomagnetic and chromoelectric dipole moments was computed in Refs. [144,145], and we consider here terms that are at most quadratic ind t . The SM tt cross section is known at N 2 LO accuracy [146]. Combining these results, the tt cross section in the presence of a top CEDM at √ S = 8 TeV is where the SM cross section was computed using the program TOP++ [147], and includes N 2 LO corrections and soft gluon resummation. The contribution ofd t was computed at LO, and we included only PDF errors. The cross section, andd t , are evaluated at the renormalization scale µ = m t . In the SM, NLO and N 2 LO corrections to the tt cross section are large [146], suggesting the need to include NLO corrections for the dipole operators as well [148].
Some of the Feynman diagrams showing the contribution of the top CEDM, pseudoscalar Yukawa and their interference to the associated production of a Higgs boson and a tt pair are shown in Fig. 3 (j)-(r). We computed the cross section at LO, retaining terms at most quadratic in the coefficients of the CPV operators. Taking the ratio with the SM cross section, we obtain the production signal strength where we evaluated the SM and the cross section induced byd t and Im Y t at LO in α s . We chose m t as factorization scale and Eq. (33) is expressed in terms of couplings at µ = m t . We performed the calculation with the CT10, NNPDF2.3 and MSTW08 NLO PDF sets, and included only PDF and α s error. Corrections to the Higgs branching ratios are small and can be neglected. The tt production cross section has been measured both at CMS and ATLAS. At 8 TeV [149,150] σ AT LAS where the first uncertainty is due to statistics, the second to systematics, and the third to the limited knowledge of the integrated luminosity. Current Higgs measurements can be used to infer the tth signal strength, although with large uncertainties. ATLAS and CMS reported We can get a bound ond t by demanding that the BSM cross section is less than the difference between the observations and SM predictions. At the 90% c.l., we find v 2 m t |d t | < 0.23, both from ATLAS and CMS. This bound is in agreement with the analysis of Ref. [23]. Similar bounds can be extracted from the first 13 TeV data. From the tth signal strength, we obtain Interestingly, already with current data, and notwithstanding the large uncertainties on the tth cross section, the tt and tth processes show comparable sensitivities to a top CEDM. The limit on the Yukawa is a factor of 10 weaker, v 2 Im Y t AT LAS < 1.6, 0.6 < v 2 Im Y t CM S < 2.1.  In Table 13 we show the dependence of the bounds on Im Y t andd t on M / T . While the bound on Im Y t depends mildly on the new physics scale, the strong mixing ofd t and θ causes the bound ond t to get stronger by a factor of 2.5 as M / T is increased to 100 TeV. The bounds will be improved by the LHC Run 2. Ref. [23] discusses how the bounds on the top CEDM from the tt cross section could be reduced to about 8%, in particular by studying differential distributions, and focusing on events with large tt invariant mass.
ATLAS and CMS project to reach a 30% (15%) uncertainty on the tth signal strength with 300 fb −1 ( 3000 fb −1 ) data. Assuming an observed signal strength compatible with the SM, a 30% accuracy allows to improve the bounds to better than 10% level v 2 m t |d t | < 0.06.
At 14 TeV, more information could be gained by looking at differential distributions. In Fig.  4 we show the differential cross section with respect to the invariant mass of the tt pair, m tt , induced by the SM and the top CEDM. We work at LO, and setd t to the projected maximum value allowed by the gluon fusion cross section at 14 TeV, |v 2d t /m t | < 0.04. While the total tth cross section will not be able to improve this limit, at large invariant mass the contribution of the top CEDM increases, being 60% of the SM at 1.5 TeV, and overtaking the SM for m tt > 2.5 TeV. Thus, the study of events at large m tt could provide a route to further improve the bound ond t . These considerations, of course, are valid only under the assumption that new degrees of freedom generating the non-standard chromo-electric top couplings are sufficiently heavy so that a local operator analysis provides a good description of the process pp → tth. Assuming that the top CEDM operator is generated at loop level by new particles with a common mass m * this criterion roughly speaking implies that m tt + m h m * . This should be kept in mind when analyzing the range of applicability of Fig. 4. Since current and prospective EDM bounds on v 2d t /m t are consistent with the mass scale m * being in the multi-TeV range (depending on coupling strengths), there are classes of models in which Fig. 4 remains valid all the way to Figure 5: Bounds ond t and Im Y t from gluon fusion, tt production, and tth production. Solid lines denote current bounds. For each observables we show the most stringent limit, namely the CMS limit for gluon fusion, and the ATLAS limit for tth. The tt bound is approximately the same for the two experiments. Dashed lines denote the projected bounds from the LHC Run 2, assuming that µ ggF = µ tth = 1 with 10% and 30% uncertainty, respectively. The projected tt bound relies on the analysis of Ref. [23] . m tt = 3 TeV. Fig. 5 summarizes the limits on Im Y t andd t set by gluon fusion, tt and tth. Solid lines denote current limits, and we show the most stringent bounds, that is the CMS bound for gluon fusion, and the ATLAS bound for tth. Dashed lines denote the projected bounds from LHC Run 2. An interesting feature of Fig. 5 is that the three observables we considered show comparable sensitivities tod t and, to a lesser extent, Im Y t . Thus, in the presence of a significant deviation from the SM in any of these three observables, it will be possible to look for signals in the remaining two, and gain more insight on the origin of the deviation. EDM constraints on this combination of couplings are discussed in Sect. 5.4.

The top chromo-magnetic dipole moment
In Sect. 4.2 we found that the strongest collider limits ond t arise from the Higgs gluon fusion production cross section. Here we briefly step aside from the main focus of the paper on CPV operators, to remark that a similar observation applies to the CP-conserving top chromo-magnetic dipole moment (CMDM) [91].
We consider the effective Lagrangian involving the top quark and the Higgs.
Eq. (40) contains the SM Yukawa coupling, and three dimension-6 operators, closely related to those discussed in Eq. (3). The first is a correction to the top Yukawa, the second is a coupling of the Higgs to the gluon field strength, andc t is the top CMDM. In the SM, Re Y t ,c t and c H vanish, and the hGG operator is generated below the top threshold, with coefficient 1/v 2 in the normalization of Eq. (40). Thus, the effective Lagrangian (40) leads to a signal strength where this expression is valid at LO in α s . The dimension-6 operators in Eq. (40) mix, in exactly the same way as their CPV analogs.
where the differences with Eq.
At LO, the signal strength induced by the top CMDM is As was the case for the top CEDM, the corrections to the Higgs width induced by the top CMDM are less important, and can be neglected. The requirement that the signal strength is in agreement with the observed µ ggF can be satisfied in two ways. The CMDM can be negative, and large enough to cancel the SM contribution to gluon fusion. This is achieved for a small interval around v 2c t /m t = −0.25, a value already excluded by the LHC and Tevatron measurements of the tt inclusive cross section [23]. The other solution is where the limit is on the coefficient at the scale M / T = 1 TeV. The limits are stronger than those ond t , because µ ggF has a linear dependence onc t , whiled t only contributes quadratically.
The bounds in Eq. (45) cut significantly into the region allowed by the tt cross section, v 2c t /m t ∈ (−0.10, 0.05), and are already competitive with the projected bounds from the LHC Run 2, studied in Ref. [23]. Furthermore, the reduction of the uncertainty on µ ggF to the 10% level at Run 2 would improve the limits in Eq. (45) to |v 2c t /m t | < 0.006. The Higgs gluon fusion production cross section thus appears to be the ideal place to look for anomalous top couplings.
This observation deserves two specifications. First of all, the contribution ofc t to c H is generated by running between the scale M / T , where we assume any BSM degrees of freedom that t/mt Current 1.2 · 10 −2 1.0 · 10 −2 0.6 · 10 −2 0.7 · 10 −2 15 · 10 −2 0.5 · 10 −2 0.27 5.2 · 10 −2 LHC Run 2 0.6 · 10 −2 0.5 · 10 −2 0.4 · 10 −2 0.4 · 10 −2 12 · 10 −2 0.4 · 10 −2 0.21 4.0 · 10 −2 generate the top CMDM to be integrated out, and m t . If the scale separation between M / T and m t is small, we cannot rely only on logarithmically enhanced terms, and need to consider also the finite contributions ofc t to c H with β 2 = 1 − Secondly, the strong bounds in Eq. (45) assume that only one coupling,c t , is turned on at M / T , so that the contribution ofc t to µ ggF is not influenced by other terms. It therefore remains important to look for direct effects ofc t . In particular, since the sensitivity of the tt and tth cross sections onc t is only moderately weaker than µ ggF , these three observables constitute ideal orthogonal probes to pin down a top CMDM.

Summary of collider bounds
In Table 14, we summarize current bounds on θ , Im Y q , Im Y t andd t that can be extracted from measurements of the Higgs production and decay processes, of the tt cross section, and of the tth signal strength performed by the ATLAS and CMS collaborations during the LHC Run 1. For each coupling, we listed the strongest bound. In the second row, we summarize the projected bound from the LHC Run 2, at 14 TeV center-of-mass energy and with integrated luminosity of 300 fb −1 . The projected bounds are obtained by assuming µ ggF and µ tth to be in agreement with the SM, with uncertainties of 10% and 30% respectively [134,135].
The current measurements of the Higgs production cross section and branching ratios allow one to bound the pseudoscalar Yukawa couplings of the light quarks at the level of 0.5 to 1%, that is, they exclude pseudoscalar couplings much bigger than the SM bottom Yukawa. The higher luminosity of the LHC Run 2, and the consequent reduction of the uncertainties on the signal strength to the 10%-20% level, will allow to improve these bounds, especially for lighter quarks.
A comparison with the EDM constraints in Table 6 shows that collider cannot compete with EDM constraints on the pseudoscalar Yukawa of the first generation. Indeed, EDM bounds already forbid Im Y u,d larger than the SM u and d Yukawas. The EDM bound on Im Y c is very close to the collider bound. While the LHC Run 2 will improve the bound on Im Y c by at most a factor of two, the next generation of EDM experiments will probe this coupling at the 10 −5 level, out of the reach of collider experiments. It is nonetheless important to pursue direct probes of this coupling, for example by studying decays of the Higgs to cc [151]. Current collider and EDM bounds on Im Y s and Im Y b are comparable, with the LHC having a slight edge. In the case of Im Y s , EDMs are not very constraining because of the poor knowledge of the nucleon matrix element of the strange EDM (and CEDM). A modest improvement on the theory, coupled with the next generation of EDM experiments, will put Im Y s out of the reach of collider experiments. Im Y b is more interesting because even with improved theory, the collider and EDM bounds are of approximately the same size. It is therefore important to get as many handles on Im Y b as possible, by studying inclusive Higgs production, the decay h → bb [20], and the associated production of h and a bb pair, with tagged b jets [126,152,153].
The coupling of the Higgs to GG is constrained by µ ggF at the 30% level, with possibility to improve to 20% in Run 2. Also in this case, collider experiments are competitive with EDM bounds.
Finally we discuss the CPV couplings of the top quark, Im Y t andd t . It is interesting that the current collider bounds on Im Y t andd t are dominated by the contribution of top loops to hGG. In the case of Im Y t , the bound from gluon fusion is a factor of ten stronger than the direct bound via µ tth . In the case of the top CEDM tt and tth probed t at the same level, a factor of 3 weaker than the bound from µ ggF . However, especially in the presence of a CEDM, tt and tth have a greater chance of improvement at Run 2, getting much closer to the gluon fusion bound. Furthermore, were significant deviations from the SM to be observed in µ ggF , µ tth or in the tt cross section, the fact that these three observables have roughly the same sensitivity tõ d t would offer the exciting possibility to prove or exclude that the origin of the signal is a top CEDM. Comparing to EDM bounds, Table 6, we see that Im Y t is strongly constrained by the electron EDM (although this constraint strongly depends on the SM prediction of the electron Yukawa coupling, see the discussion in Ref. [20]). On the other hand, the bounds ond t are very close to the LHC bounds, which makes the study of the top CEDM (and CMDM) even more interesting.
We conclude by noting that our analysis of collider observables has focused on CP-even observables that are sensitive to the square of CPV couplings. More information could be gained by studying differential observables, see Fig. 4, or observables such as spin correlations that are linear in the CPV couplings (see for example [94,119,154]).

Direct vs indirect constraints: interplay of couplings
Although it is interesting to study constraints on the individual CPV dimension-6 operators, in most BSM realizations several will be generated at the same time. Clearly a single EDM experiment can only constrain or identify a single combination of operators and several measurements are needed to isolate the individual couplings. In this section we study how EDM and collider experiments can constrain or identify combinations of CPV couplings and to what extent various experiments are complementary. We also focus on the role of the uncertainties in matrix elements which have a strong impact on the EDM analysis. As there are many combinations of CPV operators that can be studied, we focus here on a subset of cases which, in our opinion, are most interesting. Other combinations or specific BSM scenarios can be studied in similar fashion.

Im
We begin the analysis by studying the case where CPV predominantly occurs in the interactions between the lightest two quarks and the Higgs field. As shown in the previous sections, there are no significant constraint from d e or collider experiment and we therefore focus on hadronic and nuclear EDMs. In the left panel of Fig. 6, we show 90% c.l. contours in the Im Y u -Im Y d plane, arising from the current d n and d Hg bounds. The solid lines correspond to the central values of the matrix elements. In this case, the two EDM experiments are very complementary because d Hg is dominated by CPV pion exchange proportional toḡ 1 . In the conservative case (dashed lines) the complementarity is reduced leading to a significantly larger contour. The loss of complementarity is amplified once we apply the minimization strategy as can be seen in the right panel of Fig. 6. A free direction emerges indicating that large values of Im Y u and Im Y d cannot be excluded. Note that, as expected, the central and conservative contours always lie inside the minimized contour.
Additional information is needed to eliminate the free direction in the Im Y u -Im Y d plane. On the theoretical side, improved matrix elements definitely help. Using the benchmark matrix elements given in Sec. 3.4.1, we obtain the thick dashed contour in the right panel. We see that a modest improvement on the matrix elements (in most cases 50% uncertainty is sufficient) would already greatly improve the bounds and remove the free direction.
Alternatively, we can study additional EDM measurements. This could be achieved by improving d Hg by orders of magnitude, but this is unlikely to happen. Improving d n alone would still allow for a free direction. Instead we study the impact of EDM measurements on different systems. In the left panel of Fig. 7, we show constraints from d p , d D , and d Ra using central

Im
We now focus on the CPV Yukawa couplings of down-type s and b quarks. We show the central constraints in the left panel of Fig. 8, for the case with a positive Weinberg matrix element. It is clear that the electron EDM mainly constrains one of the couplings, Im Y b , while the neutron and mercury EDMs probe nearly the same combination of couplings. In case of a negative d W matrix element, the neutron and mercury constraints are somewhat more complementary, leading to slightly stronger constraints. We show constraints resulting from the minimized strategy in the right panel of Fig. 8. In this case the Im Y s direction becomes unconstrained due to the large uncertainties related to the matrix element of the strange EDM, while the constraint in the Im Y b direction is hardly affected. We find that the free direction would not be eliminated by a measurement of d p (or d D , d Ra , d Xe ) at the current d n sensitivity, although this would improve the constraint on Im Y b . In contrast, better knowledge of the strange EDM matrix element does eliminate the free direction.  Fig. 6). The right panel shows in black the current combined minimized EDM constraints. The thick-dashed contour shows the combined limit with improved matrix elements, while a future d p measurement is shown in green. The red ellipse represents the collider constraints.
In this case, current collider bounds, denoted by the red ellipse on the right panel of Fig.  8, are stronger than EDM bounds in the minimized strategy. In particular, they constrain the free Im Y s direction. Even in the presence of better theoretical handling of the strange matrix elements, the bounds from collider still play an important, complementary role. The study of the Higgs signal strengths at the LHC Run 2 will improve the bounds on Im Y b and Im Y s by 20% and 40%, respectively. Further improvement could come from the study of exclusive decays of the Higgs into bb or ss mesons [4,155,156].

Im
We now turn to the anomalous Yukawa couplings of the third generation. In this case, the constraints depend strongly on the sign of the nucleon matrix element for the Weinberg operator d W . In Fig. 9 we show the results for the least constrained case, the case of a negative neutron matrix element. We present constraints from d e , d n , and d Hg using central matrix elements in the left panel of the figure. The constraints originate mainly from the interplay between d n and d e : the anomalous top-Yukawa couplings is strongly constrained by d e , while Im Y b is mainly constrained by d n .
In the right panel of Fig. 9 we show constraints using the minimization strategy. In this case the constraints are significantly weaker leading to an almost free direction. More precise matrix elements could significantly improve the constraints as can be seen from the black dashed contour. A d D measurement would give a similar, but weaker, constraint (orange dashed ellipse).
Also in this case, the collider bound on Im Y b plays an important role, by eliminating the free direction in the Im Y b -Im Y t plane. The combined current LHC-EDM bound is indeed better than the projected EDM bound with improved matrix elements. The LHC Run 2 is likely to probe the pseudoscalar top Yukawa at the 10% level, still too far from EDM bounds to be |d D | ≤ 2.9·10 -26 e cm ATLAS Figure 9: The left panel shows the 90% c.l. contours using central matrix elements and current EDM measurements (notation as in Fig. 6). The right panel shows in black the combined minimized EDM constraints with current matrix elements (solid) and improved matrix elements (dashed). The impact of a future d D measurement is shown in orange. The horizontal band denotes the collider constraints.
relevant. As discussed in the Im Y b -Im Y s case, it will be important to get as many handles as possible on the bottom quark Yukawa.

Im Y t -d t
Finally, we consider the case that BSM physics contributes mainly to the top CEDM,d t , and its anomalous Yukawa coupling, Im Y t . In this case, the sign of the d W matrix element does not affect the constraints too much. We use a negative value which gives the weakest constraints. We show EDM constraints in the left panel of Fig. 10 using central matrix elements. The plot looks similar to the Im Y t -Im Y b plot, with the strongest constraint ond t (Im Y t ) arising from d n (d e ). After minimizing over the matrix elements, see the right panel in Fig. 10, the constraints become weaker by roughly an order of magnitude in thed t /m t direction. The constraint in the Im Y t direction is much less affected because it mostly arises from d e where the uncertainties are smaller. The reduced sensitivity could be almost completely overcome with improved matrix elements, as can seen by the dotted contour. Finally, we show the impact of a d p measurement at the current d n level. A d D measurement would be complementary as well, but is sensitive to cancellations in the sum of nucleon EDMs, d n + d p , in case of the Weinberg operator that is induced byd t (see also the discussion in Sect. 3.4.2). Collider constraints from the gluon fusion process, depicted in red, are very close to the minimized contour, but slightly too weak to have an impact. Constraints from tt and tth are at the moment not competitive with gluon fusion or with EDM constraints. However, they have the largest margin of improvement at the LHC Run 2, and are likely to become relevant, especially in the absence of theoretical improvement on the hadronic matrix elements.
|d p | ≤ 2.9·10 -26 e cm gg→h (CMS) Figure 10: The left panel shows the 90% c.l. contours using central matrix elements and current EDM measurements (notation as in Fig. 6). The right panel shows in black the combined minimized EDM constraints for current (solid) and improved (dashed) matrix elements. A future d p measurement is shown in green.

Discussion
In this paper we have presented a detailed study of both direct and indirect constraints on non-standard CP-violating Higgs couplings to quarks and gluons. Working within a linear EFT framework, we have focused on the leading flavor-conserving dimension-6 operators coupling the Higgs doublet to quarks and gluons, namely the CPV Yukawa couplings, the chromo-electric dipole operators and the ϕ † ϕGG operator (see Section 2). We have first obtained bounds on the effective couplings by assuming that at the high-scale (where we match the EFT to the underlying new physics model) only one coupling at the time dominates (Sections 3 and 4). A summary of current and prospective bounds from EDMs is provided in Table 9. Similarly, in Table 14 we present a summary of current and prospective bounds from single Higgs production and decay, and from tt and tth production at the LHC. In Section 5 we have then studied a few selected cases in which two couplings dominate at the matching scale, summarized in Figs. 6-10. Throughout our analysis, we have payed special attention to the theoretical uncertainties. For the extraction of direct bounds from LHC production cross sections, uncertainties arise from α s , parton distribution functions, and the residual dependence on the renormalization and factorization scales. The uncertainty estimate is straightforward, and typically leads to effects of O(10%) (see Section 4). 11 On the other hand, the non-perturbative matrix elements at the hadronic and nuclear level pose a greater theoretical challenge. The uncertainties corresponding to model calculations are quite large, and in some cases not even the sign of a matrix element is determined reliably. We have obtained bounds with different treatments of theoretical input. These treatments include two extreme cases: (a) Taking the central value of the matrix elements, ignoring the uncertainty. The resulting bounds reflect the maximal physics reach of EDM experiments. (b) Assuming that the matrix elements have a flat distribution in a certain range corresponding to existing calculations: this is the range-fit method used in Ref. [87] in the context of fits to the elements of the Cabibbo-Kobayashi-Maskawa matrix. The resulting bounds account for the theoretical uncertainties in the safest possible way (perhaps over-conservative, but realistic). We have also explored a third option (c): using the range-fit method with reduced theoretical errors, at the 25-50% level (see discussion in Section 3.4.1), anticipating progress in the next few years from both lattice QCD and nuclear many-body theory.
We would like to highlight the following points of our analysis: • Concerning the EDMs, a key result of our study is that hadronic and nuclear uncertainties greatly dilute the nominal constraining power (i.e. the one obtained by using central values for all matrix elements, ignoring their uncertainty). The dilution effect comes about because a given high-energy coupling generates via RGE and threshold corrections a number of operators at low-energy, whose contribution can cancel each other due to the poorly known matrix elements. From Table 6 one sees that when going from central values of the matrix elements to range-fit method, the 199 Hg bounds essentially disappear, while the neutron bounds are weakened by up to an order of magnitude, depending on the coupling, and they are eliminated in the case of Im Y s,b . Nonetheless, when considering all existing EDM constraints (d e , d n , d Hg ), it is still possible to obtain bounds on non-standard Higgs couplings. Using current theoretical uncertainties these bounds are summarized in the next-to-last row of Table 6. The bounds on Im Y s,b,t are currently determined by the ThO EDM limit, while the bounds on Im Y u,d,c , θ , andd t are set by the neutron EDM limit. This complementarity is also quite evident in the two-couplings analysis of Section 5.
• Another noteworthy result is that with the improved matrix element precision advocated in Section 3.4.1 the bounds obtained with the range-fit method come very close to the ones obtained with central value matrix elements. That is, comparing the row "Comb. Cen." and "Future Min." in Table 6, most numbers only differ by a factor 2. The exceptions are the bounds on Im Y u and Im Y d which are different by a factor 3 and 5, respectively. This follows from the fact that once matrix elements are known at the 25-50% level, there is very little room for cancellations and one essentially exploits the full power of experimental constraints. We reiterate here that the desirable target uncertainties for hadronic and nuclear matrix elements are: (i)  12 These targets do not seem unrealistic and further motivate systematic studies of hadronic and nuclear matrix elements with lattice QCD and modern nuclear many-body methods.
• In Table 9 we summarize possible future scenarios, by taking into account (i) improved matrix elements according to the benchmarks of Section 3.4.1; (ii) improved sensitivities on existing systems (d e , d n ) and EDM measurements in additional systems (d Xe , d Ra , d p , d D ); and combinations of both (i) and (ii). This exercise shows that improving the theory can have as much, or even more, impact as additional measurements. Note that this point is also evident from the plots in Section 5. For the couplings under consideration here, the anticipated improvements in the neutron and ThO EDMs (as well as the addition of proton and deuteron EDMs at the level of 10 −29 e cm) will have the largest impact. In any case, regardless of which new experimental probe becomes available, the constraints on couplings dramatically improve by using more accurate matrix elements as discussed above.
• Our analysis of collider observables has focused on (total) production and decay processes, that are sensitive to the square of CP-odd couplings. Additional information could be gained by studying more differential observables, as briefly showed in Fig. 4, or observables, such as spin correlations, that depend linearly on the new physics couplings.
• We noticed that the Higgs gluon fusion production cross section provides a very strong bound on the top CMDM, better than the direct constraint from the tt total cross section.
With the integrated luminosity of the LHC Run 2, in the absence of deviations from the SM, the limit from gluon fusion will significantly improve, to better than 1%. If multiple couplings are generated at the high scale M / T , the gluon fusion, tt and tth cross section provide complementary observables, ideal to pin down a top CMDM (and CEDM).
• Complementarity of EDMs and LHC constraints: Currently, our best knowledge of the non-standard CPV Higgs couplings comes from a combination of EDMs and LHC constraints, summarized in Table 15. 13 The strongest constraints ond q =t and Im Y u,d,t arise (by far for the light flavors) from EDMs, while for Im Y s,c,b ,d t and θ the current bounds from EDM and LHC are comparable, once we take into account the uncertainties in hadronic and nuclear matrix elements. In all cases, except for Im Y b , improved matrix elements would strengthen the current EDM constraints and put these couplings out of the reach of LHC, at least with the observables we considered. Because EDM and LHC experiments probe different combinations of couplings, they complement each other in cases where more than one coupling is simultaneously generated at the high scale. In Sect. 5, we have studied several cases where, with the current status of the hadronic and nuclear theory, only by combining LHC and EDM constraints significant constraints are obtained.
Looking to the future, the prospects for improving bounds on the non-standard couplings from EDMs are excellent, especially if experimental progress will be accompanied by improved matrix elements. On the other hand, the bounds obtained from the LHC will improve little with increased center-of-mass energy and luminosity. Although the constraints on Y u,d,s are expected to become more stringent by up to a factor of two, the expected improvements for Y c,b,t , θ andd t are more modest, see Table 14. The reason for this is that the additional non-standard contributions to Higgs production induced by Im Y q and θ grow with the center-of-mass energy more slowly or at the same rate as the SM gluon fusion cross section. So, significant improvements will be possible only with a substantial reduction of the uncertainties on the SM gluon fusion cross section. Better prospects exists for the top CEDMd t , in which case the tt and tth cross sections grow faster than the SM, and additional information can be extracted from the shape of differential distributions. As a result, as can be seen from Tables 9 and 14, anticipated 13 For ease of comparison with the existing literature, we also quote here the bounds on the non-standard Yukawa couplings in terms of the parametersκq, defined by L = (mq/v)κqqiγ5q h. Multiplying the entries of Table 15 by v/mq we obtain:κu < 0.45,κ d < 0.11,κs < 37,κc < 2.7,κ b < 0.5,κt < 0.01.
In this work we have focused on new CPV couplings of the Higgs to quark and gluons. In light of the upcoming Run 2 at the LHC and EDM searches with improved sensitivities, we think it will be timely to systematically analyze all possible CPV Higgs couplings. In this context, several new directions are worth exploring. First, as evident from our discussion, it would be interesting to study observables involving the Higgs at the LHC, that are linearly sensitive to the non-standard couplings. Second, in a framework in which the observed Higgs is part of an EW doublet, additional CPV operators appear at dimension-6, generating CPV Higgs couplings involving electroweak bosons and fermions [7,8]. We plan to study these in a subsequent work, focusing again on the best information that can be obtained from both direct and indirect probes. Finally, it would be interesting to perform a comparative analysis of the linear EFT versus the more general EFT based on the EW chiral Lagrangian with a light Higgs [13,14]. In this framework, the non-standard (possibly CPV) Yukawa couplings rise to the level of leading order couplings, and some symmetry relations are lost (e.g. in the dipole operators the coefficients of O(h 0 ) and O(h) are independent). In this context it would be very valuable to identify experimental tests involving a combination of EDMs and LHC observables that would discriminate between the two scenarios, and thus shed light on the nature of the Higgs and electroweak symmetry breaking.   (4)) at low energies, Λ χ 1 GeV. Here we assumed the scale of new physics to be M / T = 1 TeV. A dash, " − ", indicates no, or a negligible, contribution.
At the LHC thed q =t operators contribute to single-Higgs production and can be bound at the level of vd q ∼ 4-20%. Much stronger constraints arise from EDMs (four to seven orders of magnitude stronger), and their analysis can be performed in a similar way as for the top CEDM. The evolution ofd q =t to low energies is well described by the RGEs of (d q ,d q , d W ) [35]. This means the up, down, and strange CEDMs only give rise to the up, down, and strange EDMs and CEDMs, respectively, at low energies. Instead, the charm (bottom) CEDM induces a threshold correction to the Weinberg operator at m c (m b ), see Eq. (10). In turn, the induced Weinberg operator generates all the light quark (C)EDMs, d u,d,s andd u,d,s , when evolved to Λ χ [30]. This gives rise to the contributions to the operators at Λ χ shown in Table 16.
The operators at Λ χ can again be related to EDMs as discussed in Sec. 3.2. The resulting EDM constraints are presented in Table 17. A clear difference with the bounds on the Yukawa couplings andd t is that the electron EDM does not constrain any of the operators considered in this appendix. In addition, the constraints from d Hg vanish when applying the minimization procedure, as was the case in Table 6. Nonetheless, there are significant constraints on most quark CEDMs from the neutron EDM. The exception isd s , which remains unconstrained in the minimized case due to the uncertain d s matrix element. It is important to note that the constraints in the "Future Min." row differ from those in the "Comb. Cen." row only by a factor of two for most couplings, and no more than a factor of 5. Thus, an improvement of the matrix elements, as described in Sec. 3.4.1, would again allow one to exploit the full potential of the experimental limits.