The 2-Higgs-Doublet Model with Soft CP-violation Confronting Electric Dipole Moments and Colliders

We analyze CP-violating effects in both Electric Dipole Moment (EDM) measurements and future analyses at the Large Hadron Collider (LHC) assuming a 2-Higgs-Doublet Model (2HDM) with"soft"CP-violation. Our analysis of EDMs and current LHC constraints shows that, in the case of a 2HDM Type II and Type III, an $\mathcal{O}(0.1)$ CP-violating phase in the Yukawa interaction between $H_1$ (the $125~\textrm{GeV}$ Higgs boson) and fermions is still allowed. For these scenarios, we study CP-violating effects in the neutron EDM and $t\bar{t}H_1$ production at the LHC. Our analysis shows that such an $\mathcal{O}(0.1)$ CP-violating phase can be easily confirmed or excluded by future neutron EDM tests with LHC data providing a complementary cross-check.


Contents
1 Introduction CP-violation was first discovered in 1964 through the K L → ππ rare decay channel [1]. Later, more CP-violation effects were discovered in the K-, B-, and D-meson sectors [2,3] and all the discovered effects are consistent with the explanation given by the Kobayashi-Maskawa (KM) mechanism [4]. However, the KM mechanism itself cannot generate a large enough matter-antimatter asymmetry in the Universe. Therefore, new CP-violation sources beyond the KM mechanism are needed to explain the latter [5][6][7][8].
Experimentally, all the discovered effects of CP-violation till now have appeared in flavor physics measurements, yet they can also be tested through other methods. These can generally be divided into two different categories: (a) indirect tests, which can merely probe the existence of CP-violation but cannot confirm the source(s) behind it; (b) direct tests, which can directly lead us to the actual CP-violation interaction(s).
For indirect tests, there is a typical example that one most often uses, the Electric Dipole Moment (EDM) measurements [9][10][11][12][13]. The reason is that the EDM effective interaction of a fermion is wherein d f is the EDM of such a fermion f , which leads to P-and CP-violation simultaneously [10]. It is a pure quantum effect, i.e., emerging at loop level and, in the Standard Model (SM), the electron and neutron EDMs are predicted to be extremely small [10], |d SM e | ∼ 10 −38 e · cm, |d SM n | ∼ 10 −38 e · cm, (1.2) because they are generated at four-or three-loop level, respectively. Thus, since the SM predictions for these are still far below the recent experimental limits [14][15][16][17] |d e | < 1.1 × 10 −29 e · cm, |d n | < 1.8 × 10 −26 e · cm, (1.3) both given at 90% Confidence Level (C.L.) 1 , these EDMs provide a fertile ground to test the possibility of CP-violation due to new physics. In fact, in some Beyond the SM (BSM) scenarios, the EDMs of the electron and neutron can be generated already at one-or two-loop level, thus these constructs may be already strictly constrained or excluded. In measurements of d e and d n , though, even if we discover that either or both EDMs are far above the SM predictions, we cannot determine the exact interaction which constitutes such a CP-violation. For direct tests, there are several typical channels to test CP-violation at colliders. For instance, measuring the final state distributions from top pair [18][19][20][21][22][23][24][25][26][27] or τ pair [28][29][30][31][32][33][34] production enables one to test CP-violating effects entering the interactions of the fermions with one or more Higgs bosons. The discovery of the 125 GeV Higgs boson [35][36][37] makes such experiments feasible. Indeed, if more (pseudo)scalar (or else new vector) states are discovered, one can also try to measure the couplings amongst (old and new) scalars and vectors themselves to probe CP-violation entirely from the bosonic sector [38][39][40]. At high energy colliders, the discovery of some CP-violation effects can lead us directly to the CPviolating interaction(s), essentially because herein one can produce final states that can be studied at the differential level, thanks to the ability of the detectors to reconstruct their (at times, full) kinematics, which can then be mapped to both cross section and (charge/spin) asymmetry observables.
Theoretically, new CP-violation can appear in many new physics models, for example, those with an extended Higgs sector [41][42][43][44][45][46]. Among these, we choose to deal here with the widely studied 2-Higgs-Doublet Model (2HDM) [45], which we use as a prototypical source of CP-violation entertaining both direct and indirect tests of it. Specifically, the 2HDM with a Z 2 symmetry is used here, in order to avoid large Flavor Changing Neutral Currents (FCNCs) [45], yet such a symmetry must be softly broken if one wants CP-violation to arise in this scenario. We will therefore study the effects of such a CP-violating 2HDM onto (electron and neutron) EDMs as well as processes entering the Large Hadron Collider (LHC), specifically, those involving the production of a top-antitop pair in association with the 125 GeV Higgs boson. This paper is organized as follows. In section 2, we review the construction of the 2HDM with so-called "soft" CP-violation with the four standard types of Yukawa interactions. Then, in section 3, we discuss the current constraints from electron and neutron EDMs, show the reason why we eventually choose to pursue phenomenologically only the 2HDM Type II and Type III for our collider analysis and discuss the importance of future neutron EDM tests. In section 4, we discuss the current constraints from collider experiments on these two realizations of a 2HDM. In section 5, we discuss LHC phenomenology studies on CP-violation effects in the ttH 1 associated production process. Finally, we summarize and conclude in section 6. There are also several appendices which we use to collect technical details.

Model Set-up
In this section, we briefly review the 2HDM with a softly broken Z 2 symmetry and how CP-violation arises in such a model. We mainly follow the conventions in [47][48][49]. The Lagrangian of the scalar sector can be written as Under a Z 2 transformation, we can have φ 1 → φ 1 , φ 2 → −φ 2 , thus, in the scalar potential, all terms must contain even numbers of φ i . However, if the Z 2 symmetry is softly broken, a term ∝ φ † 1 φ 2 is allowed, thus the scalar potential becomes Here φ 1,2 are SU(2) scalar doublets, which are defined as (2. 3) The parameters m 2 1,2 and λ 1,2,3,4 must be real, while m 2 12 and λ 5 can be complex. Further, v 1,2 are the Vacuum Expectation Values (VEVs) of the scalar doublets with the relation |v 1 | 2 + |v 2 | 2 = 246 GeV. The ratio v 2 /v 1 may also be complex 2 , and we define t β ≡ |v 2 /v 1 | as usual 3 .
As was shown in [45], CP-violation in the scalar sector requires a nonzero m 2 12 . For the three possible complex parameters m 2 12 , λ 5 , and v 2 /v 1 , we can always perform a field rotation to keep at least one of them real. In this paper, we choose v 2 /v 1 to be real (thus both v 1,2 are real) like in [47][48][49], and have the relation following the minimization conditions for the scalar potential. If Im(m 2 12 ) and Im(λ 5 ) are non-zero, CP-violation occurs in the scalar sector.
We diagonalize the charged components as where H + is the charged Higgs boson and G + is the charged Goldstone. Similarly, for the CP-odd neutral components, where A is the physical CP-odd degree of freedom and G 0 is the neutral Goldstone. In the CP-conserved case, A is a pseudoscalar boson while, in the CP-violating case, A has further mixing with the CP-even degrees of freedom as    (2.7) 2 We can always fix v1 real through gauge transformation and v2 may be complex at the same time. 3 In this paper, we denote sα ≡ sin α, cα ≡ cos α, and tα ≡ tan α.
Here H 1,2,3 are mass eigenstates and we choose H 1 as the lightest one with mass m 1 = 125 GeV, so that it is the discovered SM-like Higgs boson. The rotation matrix R can be parameterized as When α 1,2 → 0, H 1 becomes the SM Higgs boson. If m 1,2 , α 1,2,3 and β are known, m 3 can be expressed as In the mass eigenstates, the couplings between neutral scalars and gauge bosons can be parameterized via The coefficients are then Next we turn to the Yukawa sector. Due to the Z 2 symmetry, a fermion bilinear can couple to only one scalar doublet, with the formQ L φ i D R ,Q Lφi U R , orL L φ i R , thus it is helpful to avoid the FCNC problem [45]. Hereφ i ≡ iσ 2 φ * i and left-handed fermion doublets are defined as Q i,L ≡ (U i , D i ) T L and L L ≡ (ν i , i ) T L , for the i-th generation. Since the scalar potential contains a φ 1 ↔ φ 2 exchange symmetry, we can set the convention in whichQ L U R always couple to φ 2 so that there are four standard types of Yukawa couplings [45,49]: (2.14) The fermion mass matrix is We parameterize the Yukawa couplings of mass eigenstates as For CP-violating models, c f,i are complex numbers and we list them in Appendix A for all four types of Yukawas. In all these models, Im(c f,1 ) ∝ s α 2 , thus α 2 is an important mixing angle which measures the CP-violating phase in the Yukawa couplings of H 1 .

Current EDM Constraints and Future Tests
In this section, we analyze the EDM constraints of the electron and neutron for the four types of 2HDM in some detail. The b → sγ decay requires the charged Higgs mass to be m H ± 600 GeV for all the four types of Yukawa couplings when t β ∼ 1 [50,51]. If t β gets larger, the constraints will become weaker for Type I and III Yukawa couplings. The oblique parameters [52,53] will then favor the case m H 2,3 500 GeV [54][55][56][57] 4 . With such choices for the scalar masses, the vacuum stability condition favors µ 2 ≡ Re(m 2 12 )/s 2β (450 GeV) 2 [48]. Notice that µ 2 will modify the charged Higgs couplings a little, but it is not numerically important to the EDM calculation, so we fix it at µ 2 = (450 GeV) 2 in the rest of this work. More discussions about the scalar couplings will appear in Appendix B.
An electron EDM measurement places a very strict constraint on the complex Yukawa couplings in most models. As a rough estimation, if we consider CP-violation only in the 125 GeV Higgs interaction with the top quark, the typical constraint is arg(c t,1 ) 10 −3 [58]. However, some models (including the 2HDM) allow for the accidental cancellation among various contributions, so that larger arg(c t,1 ) may still be allowed [59][60][61][62][63][64]. In such cases, neutron EDM constraints will also become important, as shown in the analysis later in this section.

Electron EDM
A recent electron EDM measurement was performed using the ThO molecule [14]. The exact constrained quantity is |d eff e | ≡ |d e + kC| < 1.1 × 10 −29 e · cm. (3.1) The second term measures the contribution from CP-violating electron-nucleon interactions via L ⊃ C N N ēiγ 5 e , (3.2) where the coefficient C is almost the same for proton and neutron. Here, k ≈ 1.6 × 10 −15 GeV 2 e · cm, which was obtained for ThO [65,66], however, for most other materials with heavy atoms, this quantity appears to be of the same order [10,67]. The contribution from electron-nucleon interactions is usually sub-leading, though it can also become important.
The typical Feynman diagrams contributing to the electron EDM in the 2HDM are listed in Figure 1. Diagrams (a)-(e) are Barr-Zee type diagrams [68] with the top quark t, W ± -boson, or charged Higgs H ± in the upper loop, while diagrams (f) and (g) are non-Barr-Zee type. Such seven diagrams contribute directly to d e . Diagram (h) shows the Figure 1. Typical Feynman diagrams contributing to the electron EDM in the 2HDM. The blue lines can be γ or Z while red lines are neutral Higgses H 1,2,3 . Diagrams (a)-(g) will contribute to d e directly while diagrams (h)-(i) will contribute to the electron-nucleon interaction term.
from now on. For the fermion-loop contribution in which the top quark is dominant, we The analytical expressions are too lengthy to present them here so that we list all of them in Appendix C. The interaction induced effective EDM terms are [66,[75][76][77] d int e,q,i = (3.14) The nucleon matrix elements O ≡ N |O|N and their values are similar for proton and neutron. Thus we choose the average values of proton and neutron considering three active quarks (u, d, s) at the hadron scale ∼ 1 GeV [77][78][79][80][81], as listed in Table 2. Summing all For each part above, d j e ∝ m e thus it is suppressed by the small electron mass. We can extract C j e ≡ d j e /(−m e ), which is independent of the fermion mass. This coefficient is not useful in the electron EDM calculation, but it will be helpful in order to map the corresponding part into the quark EDM, which is important in the neutron EDM calculation below.

Neutron EDM
The neutron EDM calculation is more complex as it involves more contributions and QCD effects. As shown in Figure 2, there are three types of operators contributing to the neutron where g s is the QCD coupling constant, t a is a generator of the QCD group and f abc denotes a QCD structure constant. At a scale µ, and For convenience we also redefine w(µ) ≡ g s (µ)C g (µ). Notice that these EDMs should be first calculated at the weak scale µ W ∼ m t .
The calculation methods of C u and C d are the same as those for d e through diagrams (a)-(g) in Figure 1. For the quark EDM, we perform the calculation at the weak scale µ W ≈ m t and list the results of the C j q evaluation [71] as follows: Here, eachC j e means C j e with a replacement c e,i → c q,i in the Yukawa couplings. The contributions including the Z boson in the Bar-Zee diagram become important in the quark EDM calculation, because there is no accidental suppression like that in the electron EDM calculation. For the CEDM terms, only Barr-Zee diagrams with a top loop contribute. The result at the weak scale µ W ∼ m t is then [71,72] (3.26) The coefficient of the Weinberg operator at weak scale is [10,72] 27) and the loop integration W (z) is listed in Appendix C.
To calculate the EDM of the neutron, we must consider the QCD running effects to evolve these to the hadron scale µ H ∼ 1 GeV. The one-loop running gives [72,[82][83][84][85]  There is no quark mass dependence in C q orC q and the evolution of C g is equivalent to w(µ H ) = 0.41w(µ W ). According to Equation 3.20, we only need the quark mass parameters at µ H ∼ 1 GeV in the final calculation. The one-loop running mass effect is [2] m q (1 GeV)/m q (2 GeV) = 1.38 (3.29) and, with the lattice results at 2 GeV [2,86,87], we have The hadron scale estimation was performed based on QCD sum rules [10,72,88,89] with an uncertainty of about 50%. Combining all these results above, we have

Numerical Analysis for the 2HDM
In this section we analyze the 2HDM with soft CP-violation, including all the four types of Yukawa interactions. For the electron EDM, the Type I and IV models give the same results, while the Type II and III models give the same results 7 . For Type I and IV models, numerical results show that there is no cancellation among various contributions to the electron EDM, thus the CP-violating phase is strictly constrained. This behavior is consistent with the results in which only the contribution from H 1 is considered [58]. We take m 2,3 ∼ 500 GeV and m ± ∼ 600 GeV as a benchmark point and find d I,IV e −1.3 × 10 −26 s α 2 /t β (3.33) in the region t β ∼ (1−10). This result is not sensitive to α 1,3 and gives |s α 2 /t β | 8.5×10 −4 , which means the CP-phase | arg(c f,1 )| 8.5 × 10 −4 for f = i , U i . This is extremely small and would not be able to produce interesting CP-violating effects, so in the rest of this work, we do not discuss further these two 2HDM realizations. For Type II and III models, in contrast, numerical results show significant cancellation behavior for some parameter regions in the electron EDM calculation and thus α 2 is allowed to reach O(0.1). For these two models, we can discuss two different scenarios: (a) the heavy neutral scalars H 2,3 are close in mass and α 3 can be changed in a wide range; (b) H 2 and H 3 have large mass splitting, and thus α 3 must be close to 0 or π/2. 6 The light quark condenzation is chosen as qq (1 GeV) = −(254 MeV) 3 [90], which is a bit larger than that from [10,88]. Ref. [90] presents the lattice result qq (2 GeV) = −(283 MeV) 3 and also shows the Renormalization Group Equation (RGE) running effect as d mqqq (µ)/d ln µ ∝ m 4 q , which is negiglible for u and d quarks. Thus we have qq (1 GeV)/ qq (2 GeV) = mq(2 GeV)/mq(1 GeV) = 0.73. 7 During the calculation of diagram (a) in Figure 1, we consider only top quark in the upper loop and ignore the small contributions from other fermions. Such approximation is good enough when t β is not too large, for example, 10. . Cancellation behavior between β and α 1 in scenario (a) of a Type II and III 2HDM. We choose m 2 = 500 GeV, m ± = 600 GeV and α 3 = 0.8 as an example. In both plots, the solid lines are the boundaries with |d e | = 1.1 × 10 −29 e · cm and the regions between solid lines are allowed by the ACME experiment while the dashed lines mean d e = 0. In the left plot, we choose a Type II model. The blue, orange, and red lines are shown for α 2 = 0.05, 0.1, 0.15, respectively. In the right plot, we fix α 2 = 0.1 and show the comparison between the Type II and Type III models. The orange lines are for the Type II model while the cyan lines are for the Type III model. We first consider Scenario (a). For Type II and III, we show the cancellation behavior of the electron EDM in Figure 3. In this scenario, the cancellation behavior is not sensitive to α 3 in a wide region (for example, 0.2 α 3 1.4). Thus, we choose α 3 = 0.8 as an example and show the cancellation behavior in the β-α 1 plane. The electron EDM sets a strict constraint which behaves as a strong correlation between β and α 1 . Numerical analysis shows that, with fixed heavy scalar masses, the location where the cancellation happens is not sensitive to α 2 , but the width of the allowed region is almost proportional to 1/s 2α 2 . We show this behavior for the Type II model in the left plot of Figure 3, for α 2 = 0.05, 0.1, 0.15, using blue, orange and red lines, respectively. The cancellation behavior in the Type III model is similar to that in the Type II model, because the Barr-Zee diagram with a bottom quark loop is negligible and thus the only difference comes from the electron-nucleon interaction part. In the right plot of Figure 3, with fixed α 2 = 0.1, we show the comparison results between the Type II model (orange lines) and Type III model (cyan lines), finding that they are almost the same. When m 2,3 increases, the location where the cancellation happens will also change slowly and we show the corresponding results in Figure 4. When m 2 increases from 500 GeV to 900 GeV, the cancellation location also moves slowly from about β 0.76 to β 0.84. The width of the allowed region is almost independent of the heavy scalar masses, as it is sensitive only to α 2 . The cancellation behavior leads to the conclusion that there is always a narrow region which is allowed by the electron EDM measurement, thus we cannot set a definite constraint on the CP-violation mixing angle α 2 only through the electron EDM, such as in the ACME experiment.
In contrast, the neutron EDM calculation does not involve such a cancellation behavior in the same region as the electron one, thus it can be used to set direct constraints on the CP- . Mass dependence in the cancellation region in the Type II model. Choosing m ± − m 2 = 100 GeV, α 3 = 0.8, α 2 = 0.1 and α 1 = 0 as an example, the black line shows the value of β satisfying d e = 0 while the dark blue region satisfies |d e | < 1.1 × 10 −29 e · cm, which is allowed by the ACME experiment at 90% C.L. If we set |α 1 | < 0.1, the light blue region is allowed. Results in the Type III model are almost the same and thus we do not show these. In the parameter region allowed by the electron EDM constraints, the CEDM of the d quark contributes dominantly to the neutron EDM. Numerical analysis shows that the neutron EDM d n ∝ s 2α 2 and it is not sensitive to α 1,3 . We calculate its .32 in the parameter region allowed by ACME experiment. We choose α 1 = 0 and α 3 = 0.8 as an example, but the modification due to these two angles is less than percent level, which is far smaller than the uncertainty in the theoretical estimation (about 50% level). In the right plot, we show the limit on α 2 in the Type II (blue) and III (orange) models. The solid lines are obtained through the estimation of central value and the dashed lines are the boundaries considering the theoretical uncertainty. If theoretical uncertainties are taken into account, we cannot set any limit on α 2 in the Type III model through neutron EDM measurements.  Type II and III models using the central value estimated in  Equation 3.32 and show the results in the left plot of Figure 5. In the Type II model, α 2 is constrained by the neutron EDM (the latest result is |d n | < 2.2 × 10 −26 e · cm at 95% C.L. [17]). Using the central value estimation in Equation 3.32, |α 2 | (0.073 − 0.088) if m 2 changes in the range (500−900) GeV, as shown in the right plot of Figure 5. Considering the uncertainty in the neutron EDM estimation [88], a larger |α 2 | ∼ 0.15 can also be allowed.
Here, we do not consider the region α 2 close to π/2 since it corresponds to the case in which H 1 is pseudoscalar component dominated, which can be excluded by other experiments, see the next section. In the Type III model, there is almost no constraint on α 2 from the neutron EDM 8 , because there is an accidental (partial) cancellation between the two terms (see Equation 3.26) in the d quark CEDM contribution, which dominates the neutron EDM calculation. Next, we discuss Scenario (b), in which a large mass splitting exists in m 2,3 , corresponding to the cases in which α 3 is close to either π/2 or 0. From Equation 2.9, we can find two solutions for t α 3 as In the large mass splitting scenario, α + 3 is close to π/2, and α − 3 is close to 0. We choose as an example m 2 = 500 GeV, m 3 = 650 GeV and m ± = 700 GeV. In the α + 3 case, H 2 is a CP-mixed state in which the pseudoscalar component is dominant, while H 3 is almost a pure scalar. Conversely, in the α − 3 case, H 3 is a CP-mixed state while H 2 is almost a pure scalar. In this scenario, the large mass splitting between H 2,3 leads to a significant H 3 → H 2 Z decay, because the coupling is just c V,1 , which is not suppressed by mixing angles. Numerical analysis shows a similar cancellation behavior as Scenario (a) in both α ± 3 cases. We show the results of the Type II model in the upper two plots in Figure 6. Similar to Scenario (a), the cancellation behavior in the Type III model is almost the same as that in the Type II model and we show the comparison in the lower two plots in Figure 6.
The behavior of the neutron EDM is also similar to that of Scenario (a). In the regions allowed by electron EDM constraint, d n is only sensitive to α 2 and is almost independent of α 1 . With the mass parameters chosen above, and using the indices II/III and +/− to denote Type II/III models and α +/− cases, we have  Figure 6. Similar to Scenario (a), the electron EDM sets a strict constraint which behaves as a strong correlation between β and α 1 . We show the cancellation behavior of the Type II model in the upper two plots and present the comparison between the Type II and III models in the lower two plots. The notation is the same as in Figure 3. The left two plots correspond to the case α + 3 in the allowed region α + 3 π/2 − 1.5 × 10 −2 α 2 while the right two plots correspond to the case α − .32. Thus, we can obtain the upper limit on α 2 in the Type II model as There is no constraint on α 2 fromthe neutron EDM in the Type III model, due to the same reason as discussed above for Scenario (a).

Future Neutron EDM Tests
Several groups are currently planning new measurements on neutron EDM, to the accuracy of O(10 −27 e · cm) or even better [11,13,[91][92][93][94]. Such an order of magnitude improvement in accuracy would be very helpful to perform further tests on the 2HDM Type II and III scenarios considered here. Figure 7. Upper limit on α 2 in the Type II and III models when the future limit decreases to |d n | < 10 −27 e · cm. The color scheme is the same as above: blue for the Type II model and orange for the Type III model. The solid lines are obtained using the central value estimation and, if we consider the theoretical uncertainty estimation of [88], the boundaries of the limits on α 2 are the dashed lines.
If no anomaly is discovered in future neutron EDM measurements, the upper limit on d n would improve to about 10 −27 e · cm, and there would be more stringent limits on α 2 in both Type II and III models, as shown in Figure 7 for Scenario (a). Further, α 2 can be constrained to O(10 −2 ) in the Type III model and to O(10 −3 ) in Type II model. Similar constraints can be placed in Scenario (b). In contrast, if α 2 ∼ O(0.1), there will be significant BSM evidence in future neutron EDM measurements. In the models which contain a similar cancellation mechanism in electron EDM, the neutron EDM experiments may be used to find the first evidence of CP-violation or set the strictest limit directly on the CP-violating phase α 2 .

Summary on EDM Tests
In the previous subsections, we have discussed the electron and neutron EDM tests in the 2HDM with soft CP-violation. There is no cancellation mechanism in the Type I and IV models and thus the electron EDM can set strict constraints on the CP-violation angle as arg(c f,1 ) s α 2 /t β 8.5 × 10 −4 . However, this value is too small to give any observable CP effects in other experiments, thus we decided not to have further discussions on these two 2HDM realizations. In contrast, cancellations among various contributions to the electron EDM can occur in the Type II and III models. Here, we still face stringent constraints but these will induce a strong correlation between β and α 1 . We cannot set constraints directly on the CP-violation mixing angle α 2 though. The behavior is the same in the Type II and III models. In fact, it is also the same in both Scenario (a), in which m 2,3 are close to each other, and in Scenario (b), in which m 2,3 have large splitting. A cancellation generally happens around t β ∼ 1 with the exact location depending weakly on the masses of the heavy (pseudo)scalars.
Current measurements of the neutron EDM can set an upper limit on |α 2 | (0.073 − 0.088) in the Type II model, depending on different scenarios and masses, if we take the central value of the neutron EDM estimation. Such limits can be weakened to about 0.15 if we consider the theoretical uncertainty. But one cannot set limits on α 2 in the Type III model, because the CEDM of the d quark in this model is suppressed by a partial cancellation. However, α 2 in the Type III model is constrained by collider tests, which will be discussed in the next section.
Finally, we showed the importance of future neutron EDM measurements in our models relying on the cancellation mechanism in the electron EDM. For α 2 ∼ O(0.1), there would be significant evidence in future neutron EDM experiments, which will be more sensitive than any other experiments. And if there is no evidence of non-zero neutron EDM, the improved limit on the neutron EDM will set strict constraints on the CP-violation mixing angle: the upper limit of |α 2 | will reach O(10 −2 ) in the Type III model and O(10 −3 ) in the Type II model.

Current Collider Constraints
Any BSM model must face LHC tests. In our 2HDM with soft CP-violation, as mentioned, we treat H 1 as the 125 GeV Higgs boson. In this scenario then, the latter mixes with the other (pseudo)scalar states and its couplings will be modified from the corresponding SM values. However, these modified couplings are constrained by global fits on the so-called Higgs signal-strength measurements. In addition, the scalar sector is extended in a 2HDM, so that direct searches for these new particles at the LHC will also set further constraints on this BSM scenario. In this respect, we discuss only the 2HDM Type II and III, in which the cancellation behavior in the electron EDM requires t β close to 1.

Global Fit on Higgs Signal Strengths
The Higgs boson H 1 can be mainly produced at the LHC through four channels: gluon fusion (ggF), vector boson fusion (VBF), associated production with vector boson (V + H 1 , here V = W, Z) or a top quark pair (tt + H 1 ) [95][96][97][98]. The decay channels H → bb, τ + τ − , γγ, W W * and ZZ * have already been discovered. Define the signal strength µ i,f corresponding to production channel i and decay channel f as follows: where σ i denotes the production cross section of the production channel i amongst those listed above, Γ f denotes the decay width of channel f and Γ tot denotes the total decay width of H 1 . A quantity with index "SM" denotes the value predicted by the SM. Such signal streengths for different channels have been measured by the ATLAS [99][100][101][102] and CMS [103][104][105] collaborations: we list them in Table 3. As intimated, in the 2HDM, H 1 couplings to SM particles are modified due to the mixing with other (pseudo)scalars and thus the aforementioned signal strengths are modified. The production cross sections satisfy [106][107][108]     We perform χ 2 -fits where where µ th i,f is the theoretically predicted signal strength, µ exp i,f is the experimentally measured one and δµ i,f is the associated uncertainty. The possible small correlations across production and decay channels are ignored. For a 2HDM, χ 2 depends only on β, α 1,2 . We perform global fits for the Type II and III models, in which α 2 ∼ O(0.1) is still allowed. The minimal χ 2 (denoted by χ 2 min ) obtained from ATLAS and CMS data as well as the combined one are listed in Table 5. The fitting, normalized to the degrees of freedom , is good enough because the models approach the SM limit when α 1,2 → 0. If one then defines δχ 2 ≡ χ 2 − χ 2 min , this is useful to find the allowed parameter regions of the two 2HDM realizations considered. Our numerical study shows that the results depend weakly on β. We choose β = 0.76 (corresponding to m 2,3 ∼ 500 GeV in Scenario (a)) as an example and show the allowed region from combined ATLAS and CMS results in the α 2 − α 1 plane in Figure 8. For both Type II and III, the global fit requires |α 2 | 0.33 in the region β ∼ (0.7 − 1). Forthe Type II model, this constraint is weaker when compared with that from the neutron EDM. Nevertheless, it can set a new constraint on |α 2 | for the Type III model. The allowed range for |α 1 | in the latter is wider than the one in the Type II model, in fact. In both models, α 1 is favored when close to 0, thus, in the following discussion, we usually fix α 1 = 0.02, a value which is not far from the best fit points in most cases. In Figure 9, we show instead the allowed regions in the α 1 − β plane for fixed α 2 = 0.1, 0.2 in the Type III model. The dependence on β is indeed weak, but it increases somewhat when α 2 gets larger, as shown in the figure.

LHC Direct Searches for Heavy Scalars
In the 2HDM, there are four additional scalars, H 2,3 and H ± , beyond the SM-like one H 1 . Thus, we must also check the direct searches for these (pseudo)scalars at the LHC.
Notice that H 2,3 decay to tt dominantly and we show their decay widths and BRs in Appendix E. The H 2,3 → 2H 1 decays are ignored because such channels are suppressed in the allowed parameter region isolated so far. In Scenario (b), H 3 → ZH 2 decay is also open if m 3 − m 2 > m Z . In addition, H − decays totb dominantly.
We first consider the process gg → H 2,3 → ZZ. Theoretically, this process is sensitive to the couplings between H 2,3 and the gauge vector bosons, hence sensitive to α 2 . Experimentally, this process is the most sensitive channel in searching for heavy neutral scalars. The current LHC limit for m 2 = 500 GeV is σ gg→H 2,3 →ZZ 0.1 pb at 95% C.L. [109] at √ s = 13 TeV with about 40 fb −1 of luminosity. In such processes, our numerical study shows that the interference between H 2,3 production and the SM background is small, thus we can safely consider only the resonant production. However, for Scenario (a), in which H 2,3 are close in mass such that |m 2 − m 3 | O(GeV) Γ 2,3 20 GeV for m 2 500 GeV (where we have denoted by Γ 2,3 the widths of the two heavy Higgs states), we must consider the interference between the H 2 and H 3 production processes. To the one-loop order, we have σ gg→H 2,3 →ZZ = σ S + σ P , (4.11) where σ S is the contribution from Re(c t,2,3 ) and σ P is the contribution from Im(c t,2,3 ).
Their ZZ invariant mass distributions are then separately given by (4.13) In the equations above, f g (x) denotes the gluon Parton Distribution Function (PDF), which, in our numerical study, is chosen to be the MSTW2008 set [110]. The function [106,107] Γ 0 (q) = q 3 α em is the decay width to the ZZ final state of a would-be SM Higgs boson with mass q. The functions [106,107]σ are the parton-level cross sections of a pure scalar(pseudoscalar) state with couplings c t = 1(i). The loop functions A 1 and B 1 are listed in Appendix D. Thus, the total cross section is For m 2 500 GeV, we choose ∆q = 50 GeV as the mass window where interference is accounted for.
Numerically, we show the cross sections depending on the mixing angles in Figure 10 by fixing m 2 = 500 GeV in the Type III model. The left plot is for Scenario (a) and the right plot is for Scenario (b) for the α + 3 case. In both scenarios, we can see that α 2 0.27 is favored when m 2 = 500 GeV. Thus, in the following analysis, we generally choose α 2 = 0.27 (unless stated otherwise) as a benchmark point, corresponding to the largest allowed CPviolation effects. For Scenario (a), when we choose α 2 = 0.27, α 3 0.4 or 1.2 is favored, which still keeps H 2,3 nearly degenerate in mass. For Scenario (b) and the α − 3 case, c 2,V is suppressed (close to α 1 ), thus it faces no further constraints here. In the Type II model, we can obtain the same cross section as that in the Type III model with the same parameters.
In the Type II model, due to the stricter neutron EDM constraint, the considered parameter space is always allowed.
As mentioned, H 2,3 decay dominantly to a tt final state and the current LHC limit for m 2 = 500 GeV is about σ pp→H 2,3 →tt 7 pb at 95% C.L. [111] at √ s = 13 TeV and 36 fb −1 of luminosity. In contrast to the ZZ channel, the interference with SM background is very important in the tt channel [112,113], which strongly decreases the signal cross section compared with the pure resonance production cross section, so long that non-resonant Higgs diagrams can be subtracted [? ]. The total cross section can be divided into Here, σ SM denotes the SM cross section while σ res and σ int denote the resonant and interference cross section, separately. Furthermore, δσ tt is the cross section difference between 2HDM and SM, i.e., whereσ denotes the parton-level cross section as a function of the tt invariant mass q.
Following the results in [112,113], we havê σ res =σ res,S +σ res,P Here, q 2 = x 1 x 2 s, β t = 1 − 4m 2 t /q 2 is the velocity of the top quark in the tt center-ofmass frame. In our numerical study, we set q in the range m 2 − ∆ q/2 < q < m 2 + ∆ q/2, where we choose the mass window ∆ q = 100 GeV for m 2 = 500 GeV. We choose the MSTW2008 PDF [110] as above. We show the cross sections for some benchmark points in both Scenario (a) and Scenario (b) in Table 6. The numerical results show that, for all benchmark points we consider, the interference with the SM background significantly breaks the resonance structure of H 2,3 and decreases the cross sections to around (even below) 0, which means the tt resonant search at the LHC cannot set limits on this model 9 .
The H 2,3 states can also be produced in association with a tt pair at the LHC, thus we should also check this constraint for our favored benchmark points. Since H 2,3 mainly decay into a tt pair, the whole production and decay process will modify the cross section of the pp → tttt process (which we denote by σ 4t ), which current LHC limit is about 22.5 fb at 95% C.L. [116] at √ s = 13 TeV with 137 fb −1 of integrated luminosity. The interference effects between SM and BSM contributions are expected to be significant [117]. We estimate this cross section in the 2HDM considering all interference effects by using Madgraph5_aMC@NLO [118,119]. We then show the numerical results in Table 7 for some benchmark points, all allowed by current LHC limits. Finally, we should also check the direct LHC limits on the charged Higgs boson H ± . As mentioned above, b → sγ decay favors a heavy H ± state with mass m ± 600 GeV [50,51]. For m ± = 600 GeV, the current LHC limit is about 0.1 pb at 95% C.L. [120,121] at √ s = 13 with some 36 fb −1 of luminosity TeV. For large t β , the interference effect is negligible [122]. However, in the Type II and Type III models with CP-violation as considered above, t β ∼ 1 is favored. For m ± 600 GeV, its width Γ ± 30 GeV, which leads to significant interference effects. Again, we estimate the cross section considering all interference effects using Madgraph5_aMC@NLO [118,119]. If we denote by δσ ± the cross section modification (including both the resonant and interference effects) to SM ttbb process, our numerical estimation show that for m ± = 600 GeV and β = 0.76. That means that the interference effect significantly decreases the H ± production cross section in this parameter region, thus the latter is not constrained by current LHC experiments.

Summary on Collider Constraints
The 125 GeV Higgs (H 1 ) signal strength measurements lead to a constraint |α 2 | 0.33, which depends weakly on β. The LHC direct searches for heavy neutral scalars decaying to the ZZ final state set a stricter constraint |α 2 | 0.27 for m 2 = 500 GeV in both Scenario (a) and (b). When m 2 (550 − 600) GeV, the constraint from direct searches becomes weaker than that from the global fit ton the H 1 signal strengths. In further analysis, we prefer to choose α 2 = 0.27, which is the largest allowed value for m 2 = 500 GeV. We have also checked the constraints from tt, tttt and charged Higgs boson searches, in which the interference effects are very important. All benchmark points that we have considered are allowed by current LHC measurements. In the remainder of this work, we focus on the phenomenology of CP-violation in ttH 1 associate production. We will instead consider the production and decay phenomenology of the heavy (pseudo)scalars H 2,3 in a forthcoming paper.

LHC Phenomenology of CP-violation in ttH 1 Production
In this section, we study the production of the neutral Higgs bosons H 1 in association with a tt pair at the LHC. We start by discussing the setup used in the calculation and finish by highlighting the results for the ttH 1 final state.

Phenomenological Setup
Events are generated at Leading Order (LO) using Madgraph5_aMC@NLO [118,119]. Cross sections of signal processes are calculated using a UFO model file [123] corresponding to the general 2HDM [124] slightly modified to account for CP-violation effects in vertices involving both the neutral (H i , with i = 1, 2, 3) and charged (H ± ) Higgs boson states. Here, we employ the LO version of the Mmhtlo68cl PDF sets [125]. For both the signal and background processes, we have used the nominal value for the (identical) renormalization and factorization scales to be equal to half the scalar sum of the transverse mass of all final state particles on an event-by-event basis, i.e.: In the computation of the parton level cross sections, we have employed the G µ scheme, where the input parameters are G F , α em and m Z , the numerical values of which are given by Uncertainties due to the scale and PDF variations are computed using SysCalc [126]. In order to keep full spin correlations at both the production and decay stages of the top quarks, we have employed MadSpin [127]. Pythia8 [128] is used to perform parton showering and hadronization -albeit without including Multiple Parton Interactions (MPIs) -to the events, eventually producing a set of event files in HepMC format [129]. The HepMC files are passed to Rivet (version 2.7.1) [130] for a particle level analysis. In the latter, jets are clustered using the anti-k T algorithm using FastJets [131,132] 10 .
The particle level events are selected if they contain charged leptons, high jet multiplicity of 4-6 jets where some of these are b-tagged (see below) and missing transverse energy which corresponds to the SM neutrino from W boson decays. Only prompt electrons and muons directly connected to the W boson are accepted, i.e., we do not select those coming from τ decays. Electrons are selected if they pass the basic selection requirement of p e T > 30 GeV and |η e | < 2.5 while muons are selected if they satisfy the conditions p µ T > 27 GeV and |η µ | < 2.4. Jets are clustered with jet radius ∆R = 0.4 and selected if they satisfy p j T > 30 GeV and |η j | < 2.4. For b-tagging, we use the so-called ghost-association technique [139,140]. In this method, a jet is b-tagged if all the jet particles i within ∆R(jet, i) < 0.3 of a given anti-k T jet satisfy p i T > 5 GeV. Furthermore, in our analysis, we select events if they contain two charged leptons (with opposite electric charge), at least four jets (where at least two of them are b-tagged) and missing energy (coming from neutrinos in the leptonic decays of both top (anti)quarks).

Inclusive ttH 1 Cross Section
The LO parton level Feynman diagrams are depicted in Figure 11. The cross section has two contributions: one from qq annihilation (diagram (a) in Figure 11), which is expected to dominate in the region of medium and large x =p i /P (with P being the longitudinal momentum of the incoming proton) and one from gg fusion (diagrams (b) and (c) in Figure 11) dominating at low x. For the calculation of the cross section, we employ Mad-graph5_amc@nlo [118,119] with the Mmhtlo68cl and Mmhtnlo68cl PDF sets [125] in the 4-flavor scheme. Systematic uncertainties are divided into two categories: scale and Figure 11. Representative Feynman diagrams corresponding to ttH 1 production at LO. They consist of production through qq annihilation (diagram (a)) and through gg fusion (diagrams (b)-(c)).
Furthermore, PDF uncertainties are estimated using the Hessian method [141].
In Table 8, we show the results of the cross section both at LO and the Next-to-LO (NLO) in the SM. We can see that the NLO corrections imply a K-factor of about 1.17 in the case when no cuts are applied on the Higgs boson transverse momentum and for the case where p H T > 50 GeV. The K-factor slightly increase to 1.25 when a more stringent cut (p H T > 200 GeV) is applied. Furthermore, the theoretical uncertainties are dominated by those associated to scale variations which significantly decrease when we go from LO to NLO. PDF uncertainties are subleading and mildly dependent on the Higgs p T cut. Finally, we notice that the ggF contribution is dominant accounting for 68( 71.5%), at LO(NLO), of the total cross section in the case of p H T > 50 GeV and slightly decreasing to 59%( 67%) for the p H T > 200 GeV case. In the 2HDM, the ttH 1 coupling is given by with c t,1 = c α 2 s β+α 1 /s β − is α 2 /t β (see Equation A.1) for all four types of models. The ttH 1 production cross section behaves as shown in Equation 4.4. The presence of the pseudoscalar part in the ttH 1 coupling can drastically changes the value of the cross section as can be seen in Figure 12.

Results
In this subsection, we show the results of the sensitivity of certain spin observables to the nature of the ttH 1 vertex. Description of the observables is given in Appendix F. In Figure 13, we show the results of the normalized distributions corresponding to spin observables in the SM (red lines) and the 2HDM with α 2 = 0.27 (blue lines). The pure pseudoscalar case α 2 = π/2 is also shown in green lines as a comparison in the figures.
To avoid clutter, we show only the most sensitive observables, i.e., the difference in the azimuthal angles of the charged leptons and the cos θ n + cos θ n − spectrum in the transverse basis.
In order to quantify the sensitivity of the various spin observables to the benchmark points, we compute forward-backward asymmetries. An asymmetry A O on the observable O is defined by 7) Figure 13. Left: Normalized differential cross section 1/σ dσ/d|∆φ + − | versus |∆φ + − | in the laboratory frame. Right: Normalized differential distribution (1/σ)dσ/d(cos θ n a cos θ n b ) versus cos θ a cos θ b in the transverse basis. The distributions are shown for the SM (red) and for the selected benchmark point, α 2 = 0.27, of the Type III 2HDM (blue). We also show the distribution for the pure pseudoscalar case α 2 = π/2 in green lines as a comparison, which is of course excluded. To quantify deviations from the SM expectations, we compute the χ 2 as with σ O the uncertainty on the measurement of the asymmetry in the SM. We assume that the asymmetry follows a Gaussian distribution, in which case the uncertainty is given by [142] 9) where N = A × σ × L. Here, A × is the acceptance times the efficiency of the signal process after full selection. In our case, we find that the efficiency is about 14% for both the SM and the 2HDM, at L = 3000 fb −1 , and σ is the cross section times the BRs, i.e., In this calculation, we assume a b-tagging efficiency of about 80%. In Table 9, we show the expected deviations from the SM expectation at L = 3000 fb −1 . We can see that, unfortunately, for α 2 = 0.27, the χ 2 cannot be larger than 0.93 considering only a single observable. After combining all the observables in Table 9, the χ 2 can reach about 4.4.
The results depends weakly on β and α 1 in our favored region (t β ∼ 1 and α 1 ∼ 0), because the observables are sensitive only to the ttH 1 CP-violating phase s α 2 /t β in this region. Possible improvements may be made by using, e.g., the cos θ in the r-basis (see Appendix F) in the "single lepton plus jets" final state.

Conclusions
In this work, we have analyzed soft CP-violating effects in both EDMs and LHC phenomenology in a 2HDM with soft CP-violation. In this scenario, the mixing angle α 2 is the key parameter measuring the size of CP-violation since the CP-violating phases in H 1 ff Yukawa vertices are proportional to s α 2 .
We have considered all four standard types of Yukawa couplings, named Type I-IV models, in our analysis. In Type I and IV models, there is no cancellation mechanism in electron EDM calculations, leading to a very strict constraint on the CP-violating phase | arg c t/τ,1 | 8.5 × 10 −4 , which renders all CP-violating effects unobservable in further collider studies for these two models.
In Type II and III models, we have discussed two scenarios: (a) H 2,3 are closed in mass while α 3 is away from 0 or π/2; and (b) H 2,3 have a large mass splitting while α 3 must appear close to 0 or π/2. The cancellation behavior in the electron EDM leads to a larger allowed region for α 2 in both scenarios. In such two models, t β is favored to be close to 1, whose location depends on the masses of the heavy (pseudo)scalars, with a strong correlation with α 1 . The electron EDM alone cannot set constraints on α 2 directly. In the Type II model, |α 2 | 0.09 is estimated from the neutron EDM constraint if we consider only the central value estimation and this constraint can be as weak as 0.15 if theoretical uncertainty in neutron EDM estimation is also considered. In the Type III model, no constraint can be drawn from the neutron EDM and |α 2 | 0.27 is estimated from LHC constraints if m 2 500 GeV. Other LHC direct searches do not set further limits for the 2HDM.
Our analysis shows the importance of further neutron EDM measurements to an accuracy of O(10 −27 e · cm). An α 2 of the size ∼ O(0.1) will lead to significantly non-zero results in such experiments. If CP-violation in the Higgs sector exists, as we have discussed, first evidence of it is expected to appear in the neutron EDM measurements. Conversely, if there is still a null result for the neutron EDM, constraints on |α 2 | can be pushed to about 4 × 10 −3 in the Type II model and 2 × 10 −2 in the Type III model. Thus, we conclude that, for models in which a cancellation mechanism can appear in the electron EDM, the neutron EDM measurements are good supplements to find evidence of CP-violation or set constraints on the CP-violating angle directly.
We have also performed a phenomenological study of soft CP-violation in the 2HDM for the case of ttH 1 associate production at the LHC with a luminosity of 3000 fb −1 . With fixed β and α 1,2 , its properties are independent of the mixing angle α 3 and the masses of the heavy (pseudo)scalars H 2,3 and H ± . Upon choosing the bencdhmark point β = 0.76, α 1 = 0.02 and α 2 = 0.27, we constructed top (anti)quark spin dependent observables and tested their deviations from the SM. Amongst these, a single observable, the azimuthal angle between the two leptons from fully leptonic tt decays, ∆φ + − , is the most sensitive one, with χ 2 = 0.93, meaning that we can hardly achieve any higher signal significance using any other single observable. After combining all the observables, χ 2 can reach 4.4.
In the light of this, then, we conclude that future neutron EDM experiments could provide more useful tests of soft CP-violation in the 2HDM than the LHC experiments.
In this paper, we did not perform phenomenological studies of the heavy (pseudo)scalars (H 2,3 or H ± ), for which interference effects with the SM background are very important and thus need a dedicated treatment. We will turn to them in a forthcoming paper. c X

B Scalar Couplings
The scalar couplings in the potential can be expressed using the physical parameters as [47][48][49] Consider the-bounded-from-below conditions as [45] then µ 2 (450 GeV) 2 is favored and thus we choose µ 2 = (450 GeV) 2 in the analysis.
The couplings between neutron and charged scalars c i,± are [49] c where R is the matrix in Equation 2.8. These couplings are useful in the calculations of fermionic EDMs seeing the contribution of a charged Higgs boson.

C Loop Integrations for EDM
The loop functions in the calculation of the Barr-Zee diagrams are [68][69][70][71][72][73][74]: Denoting the loop functions in the non-Barr-Zee type diagrams with a W boson are [70] (D a W ) i = − The loop functions in the non-Barr-Zee type diagrams with a Z boson are instead [70] (D a In the functions (D p W ) i , we have z a ≡ z W H i while, in the functions (D p Z ) i , we have z a ≡ z ZH i . Last, the loop function for the Weinberg operator is [72] (C.16)

D Loop Integrations for Higgs Production and Decay
The loop functions for Higgs production and decay are [106,107] A 0 (x) = x − I(x) x 2 , (D.1)

E Decay of Heavy (Pseudo)scalars
For heavy neutral (pseudo)scalars, we consider the decay channels H 2,3 → tt, W W, ZZ and ZH 1 . The partial decay widths are given by Here k = i or 1, and the functions In Scenario (b), since H 2,3 have large mass splitting, we should also consider the H 3 → ZH 2 decay. Its partial width is The charged Higgs boson H + decays mainly to tb in the small t β region. Ignoring the coupling term proportional to m b , we have (E.7) Besides this, H + also have subdominant decay channels, like W + H i [49], yielding For β = 0.76 and m ± = 600 GeV, Γ H + →tb = 33 GeV while the sum for all three neutral scalars i Γ H + →W + H i 5 GeV for |α 2 | 0.27.

F CP-violation Observables in the ttH 1 Channel
In this section, we give an overview of the different observables that we have used in this study to pin-down the spin and CP properties of the SM-like Higgs boson produced in association with a tt pair. First, one can study directly the spin-spin correlations of the tt pair by measuring the differential distribution in cos θ a cos θ b of the emerging leptons, where α is the spin analyzing power of the charged lepton and θ a,b = (ˆ a,b ,Ŝ a,b ), withˆ a,b being the direction of flight of the charged lepton in the top quark rest frame andŜ a,b the spin quantization axis in the basis a. Furthermore, C ab is the correlation coefficient which is related to the expectation value of cos θ a cos θ b using Equation F.1. In the following, we consider three different bases: the helicity basis (a = k), the transverse basis (a = n) and the r-basis, see, e.g., [20,143] for more details about the definitions of the spin bases and [144,145] for reported measurements of these observables in tt production. It was found that the tt spin-spin correlations in the transverse and r-bases are good probes of CP-violation, e.g., through the anomalous chromomagnetic and chromoelectric top quark couplings [143] 11 . Furthermore, we consider the opening angle between the two oppositely charged leptons produced in the decays of the top (anti)quarks which is defined by wherep + (p − ) is the direction of the flight of the charged lepton + ( − ) in the parent top (anti)quark rest frame. The azimuthal angle ∆φ + − = |φ + − φ − | is a clean observable to measure the spinspin correlations between the top and the antitop quarks. The momenta of the charged leptons are usually measured in the laboratory frame [146,147]. This observable shows a high sensitivity to the degree of correlations between the top (anti)quarks in tt production. However, since we are considering the ttH 1 production mode, the presence of the Higgs boson may wash out the sensitivity of ∆φ to the correlations, though we have found this not to be the case.
In addition to the aforementioned observables, we also study the sensitivity of the following angle [21] cos θ H 1 = (p + ×p H 1 ) · (p − ×p H 1 ) |(p + ×p H 1 )||(p − ×p H 1 )| , (F.3) 11 In ttH1 production, the contribution of ggF is about 70% of the total cross section. Hence, the initial state is mostly Bose-symmetric. Following the recommendations of [143], the value of cos θ is multiplied by the sign of the scattering angle ϑ =p ·p t withp t = p t /|p t | the top quark direction of flight in the tt rest frame andp = (0, 0, 1).
wherep + ,p − andp H 1 are the directions of flight of the postively-, negatively charged lepton and of the reconstructed Higgs boson in the laboratory frame. The θ H 1 angle defines the angle spanned by the charged lepton momenta projected onto the plane perpendicular to the Higgs boson direction of flight.
One can obtain the polarization of the (anti)top quark by integrating Equation F.1 over the angle θ a (or θ b ): 1 σ dσ d cos θ a ± = 1 2 1 + α ± P a t,t cos θ a ± , (F. 4) which applies to all the spin quantization axses used here. It was also found that the energy distributions of the top quark decay products carry some information on the polarization state of the top (anti)quark [148][149][150][151][152][153][154][155][156]. We follow the same definitions used by [150,152] and study the ratios of the different energies. We give the first two observables as follows where E , E b and E t are the energies of the charged lepton, b-jet and top quark in the laboratory frame. Finally, we consider the energy of the charged leptons in the laboratory frame where m t = 172.5 GeV is the pole mass of the top quark. Figure 14. The double differential cross section in cos θ n cos θ n in the transverse basis (left) and in cos θ k cos θ r − cos θ r cos θ k in both the helicity and r-basis (right). Red lines are for the SM, blue lines are for the signal benchmark point with α 2 = 0.27 while green lines are for the pure pseudoscalar case α 2 = π/2 as a comparison, which is of course excluded. In Figure 14, we show the cos θ k + in the helicity basis (left panel) and an asymmetric cos θ ℓH Ratio to SM combination of the same double-angle distribution in both the helicity and r-basis (right panel). The two figures show clearly significant sensitivity for α 2 = π/2, while for α 2 = 0.27 the sensitivity is rather mild. In Figure 15, we display the spectrum in the opening angle between the leptons from the fully leptonic decay of the tt system (left) and the cos θ H 1 distribution (right). The same conclusions as in the previous case apply for the case of α 2 = 0.27 while differences of about 10% with respect to the SM case can be reached for α 2 = π/2.

G Top Quark Reconstruction
For tt spin-spin correlation and polarization observables in the top quark rest frame, it is mandatory to fully reconstruct the top (anti)quark four-momentum. In this regard, we employ the PseudoTop definition [157] widely used by the ATLAS and the CMS collaborations for, e.g., validation of MC event generators. We slightly modify the Rivet implementation of the CMS measurement of the tt differential cross section at √ s = 8 TeV [158]. We minimize the following quantity to select the hadronic, leptonic (anti)top quarks and SM-like Higgs boson decaying into bb. In Equation G.1, m t , m W and m H are the masses of the top quark, W boson and the Higgs boson, respectively, whilet (t h ) is the momentum of the (anti)top constructed in the leptonic(hadronic) decays of the W boson, withp H 1 the four-momentum of the Higgs boson candidate. In the reconstruction procedure, all jets and leptons in the event are considered provided they satisfy the selection criteria which was highlighted in subsection 5.1. Validation plots for the PseudoTop reconstruction method in ttH 1 (→ bb) (green) and the QCD-mediated ttbb (red) are shown in Figure 16.