CP-Violation in the Two Higgs Doublet Model: from the LHC to EDMs

We study the prospective sensitivity to CP-violating Two Higgs Doublet Models from the 14 TeV LHC and future electric dipole moment (EDM) experiments. We concentrate on the search for a resonant heavy Higgs that decays to a $Z$ boson and a SM-like Higgs h, leading to the $Z(\ell\ell)h(b\bar{b})$ final state. The prospective LHC reach is analyzed using the Boosted Decision Tree method. We illustrate the complementarity between the LHC and low energy EDM measurements and study the dependence of the physics reach on the degree of deviation from the alignment limit. In all cases, we find that there exists a large part of parameter space that is sensitive to both EDMs and LHC searches.


I. INTRODUCTION
With the discovery of the Higgs-like boson at the LHC [1,2], the remaining particle predicted by the Standard Model (SM) has been found. Up to now, the measured properties of this new resonance show no significant deviation from the SM predictions. Nevertheless, the new boson could reside in a larger structure with an extended scalar sector that incorporates the SM. The possibilities for such extended scalar sectors abound. Among the most widely considered and theoretically well-motivated are Two Higgs Doublet Models (2HDMs). Even with the rather minimal introduction of a second SU(2) L scalar doublet, the possible phenomenological consequences of 2HDMs are rich and diverse. The possibility of new sources of CP-violation is one of the most interesting but, perhaps, less extensively studied.
Explaining the cosmic matter and anti-matter asymmetry requires the existence of additional CP-violation (CPV) beyond that of the SM. Electroweak Baryogenesis (EWBG) is one of the most compelling solutions to this problem [3][4][5]. EWBG fulfills the Sakharov conditions for successful baryogenesis [6] (B violation, out-of-equilibrium dynamics, and both C and CP-violation) through B + L violating sphaleron transitions, a strong first order electroweak phase transition that proceeds through bubble nucleation, and CPV interactions at the bubble wall. While the SM would in principle provide these ingredients, it is known that the CPV effects generated by the Cabibbo-Kobayashi-Maskawa matrix and QCD θ term are too feeble and that the SM-like Higgs scalar is too heavy for a strongly first order electroweak phase transition [7][8][9].
The 2HDMs provide possible solutions to these shortcomings. The viability of a strong first order electroweak phase transition and the favored parameter space of the 2HDMs have been studied in Refs. [10][11][12]. In the CPV sector, the LHC has already excluded the new boson as a pure CP odd scalar at 99.98% CL and 97.8% CL in Ref. [13] and Ref. [14] respectively.
If the boson is a part of the 2HDM, it could nevertheless receive a small CP-odd admixture from CP-violating terms in the scalar potential. This possibility for 2HDM CP-violation is strongly bounded by the non-observation of permanent electric dipole moments (EDMs) of the neutron, electron, and diamagnetic atoms, including mercury and radium [15][16][17][18], as analyzed recently in Refs. [19][20][21][22][23]. The authors of Refs. [24][25][26] also pointed out that LHC searches for additional, heavy scalars can be complementary to EDM searches, especially in regions of 2HDM parameter space where strong cancellations between Barr-Zee EDM diagrams occur. Nonetheless, there exists a window for sufficient CPV to generate the matter-antimatter asymmetry, as shown in Ref. [27].
In what follows, we analyze the prospects for future LHC probes of the CPV 2HDM, building on the previous studies in Ref. [19] and Ref. [24], where EDMs constraints and 8 TeV LHC results in CPV 2HDMs are analyzed in detail. We adopt the framework of CPV 2HDMs with a softly-broken Z 2 symmetry to avoid a problematic tree level flavor changing neutral currents (FCNCs). We consider future LHC searches for a heavy Higgs of mixed CP (denoted h i=2, 3 ) which decays to a Z boson and a SM-like Higgs (h 1 ), and obtain the prospective reach for Run II and the high luminosity phase (HL-LHC). We concentrate on the llbb final state, where the Z boson decays to a pair of leptons (e or µ), and the SM-like Higgs decays to a pair of b quarks, because it is one of the most sensitive channels and because the final state particles allow for a relatively high reconstruction efficiency. We first follow the cut-based analysis procedure described in Ref. [28] to reproduce the ATLAS 8 TeV results and validate our Monte Carlo signal and background generation, then use the Boosted Decision Tree (BDT) [29] method to obtain the 95% CLs exclusion limit for future 14 TeV experiments with integrated luminosities equal to 300 fb −1 , and 3000 fb −1 , respectively. We subsequently translate the prospective exclusion limits into constraints on the parameter space, and find that a large portion of parameter space can be tested with both future LHC and EDMs experiments.
From the global fit of Higgs coupling measurements [31,32], one find that the current data favor the 2HDMs to be close to the alignment limit: β − α = π/2 where α and β are defined in Sec. II B and Sec. II A respectively. Therefore, we summarize our finding in the following two categories with the combined analysis of the future EDMs and LHC exclusion bounds shown in Figs. 6 and 7.
• 2HDMs in alignment limit: With a discovery at the future LHC, the Type-I 2HDM would imply observation of non-zero radium and electron EDMs in the next generation searches, while the Type-II 2HDM would imply non-zero neutron and radium EDMs. A null result at the future LHC will still allow for the CPV 2HDMs if the CPV effect in the model is sufficiently small. Future EDM may still yield non-vanishing results in this case.
• 2HDMs away from the alignment limit: With a discovery at the future LHC, one may or may not expect nonzero EDM results depending on the level of deviation. This is due to the fact that the production of the mostly CP odd Higgs in the model (h 3 defined in Sec. II B) is sensitive to the deviation from the alignment limit which is not suppressed by a small CPV effect. As a result, the discovery at LHC may indicate a relatively large deviation from the alignment limit instead of a large CPV effect. A null result at future LHC may not exclude the CPV 2HDMs if the CPV effect and the deviation from the alignment limit are sufficiently small. For a relatively large deviation from the alignment limit, any non-zero EDM results would disfavor the CPV 2HDMs.
The above conclusions are based on the detailed analysis discussed in Sec. V. We also point out that our analysis will break down in some regions of parameter space that have both small values of tan β (ratio of the vacuum expectation values of the two neutral scalars) and the CPV Higgs mixing angle α b , where the interference between the resonant amplitude ( gg → h 2,3 → Zh 1 ) and non-resonant amplitude (box diagram gg → Zh 1 ) for Zh 1 production may become significant. We do not perform a full analysis of this effect, but rather give a qualitative estimate, as this region does not appear to significantly impact the prospective Run II exclusion reach.
The organization of our paper is the following. In Sec. II, we describe our set-up for CP violation 2HDMs. In Sec. III, we show the analytical formulas used to derive constraints on the parameter space. In Sec. IV we describe details of our simulation and analyses. In Sec. V, we exhibit future LHC constraints and discuss possible issues arising from the interference between the resonant and non-resonant diagrams. Finally, we conclude in Sec. VI. The distributions of kinematic variables used in BDT analysis are listed in Appendix A. The formulas for two-body decay rates of heavy Higgses are given in Appendix B.

II. CPV 2HDM MODEL DESCRIPTION
In this section, we describe details of the CPV 2HDM framework that will be used in the following discussions.

A. General 2HDM Scalar Potential
The most general 2HDM scalar potential containing two Higgs doublets φ 1 and φ 2 can be expressed in the following form: Two fields φ 1 and φ 2 can be expressed as with in general v 1 and v 2 complex and v = |v 1 | 2 + |v 2 | 2 = 246 GeV. We also denote that tan β = |v 2 |/|v 1 |. One can always perform a SU (2) L × U (1) Y gauge transformation to go into a basis where v 1 is real while v 2 = |v 2 |e iξ is still complex.
To guarantee that there are no FCNCs at tree level, one can assign Z 2 charges to the two Higgs doublets as well as the fermion fields such that each fermion can only couple to one of the Higgs doublets. Depending on the transformation of the fermion fields under the Z 2 symmetry, there can be various types of 2HDMs that we will introduce in Sec. II C. The Z 2 symmetry implies the potential parameters m 2 12 and λ 6,7 vanish, which in turn forbids the presence of CP phases in the potential. Therefore, we retain the m 2 12 term which only softly breaks the Z 2 symmetry. In general, this soft Z 2 symmetry breaking term together with quartic Z 2 conserving term would induce new quartic Z 2 breaking terms by renormalization, but they are at one-loop level and thus do not induce new FCNC at tree level.
Hermicity implies that there are only two complex parameters, m 2 12 and λ 5 , in the potential. With the global phase redefinition of the fields φ j → e iθj φ j , one may define two rephasing invariant phases as in Ref [19], The minimization of the potential yields that: Im(m 2 12 e iξ ) = v 2 sin β cos βIm(λ 5 e 2iξ ) .
Eq. 6 above indicates that the value of ξ is determined by given m 2 12 and λ 5 . Expressing this equation with rephasing invariant phases implies: In short, there is only one CP independent phase in the potential after electroweak symmetry breaking(EWSB). Using this rephasing freedom of the fields, we will work in a basis where ξ = 0 and encode this invariant CPV phase into a CPV angle in the diagonalization matrix for the neutral Higgs sector.

B. Higgs Mass Eigenstates
After EWSB, we can use the following relations to diagonalize the mass matrix for the charged Higgs sector, which separates the physical charged Higgs and would-be Goldstone bosons: This leads to a relationship between the mass of the charged Higgs and parameters in the scalar potential: where the parameter ν sets the hierarchy between the SM-like Higgs and charged Higgs. The mass term in the Lagrangian is given by L mass where λ 345 represents λ 3 + λ 4 + Re(λ 5 ). A rotation matrix R defined below can be used to diagonalize the mass matrix: where s α and c α are short hands for sin α and cos α. Under this rotation matrix, we have We demand that three rotation angles are in the following range: With this diagonalization procedure, one can obtain six linearly independent equations which can be solved for the parameters in the scalar potential in terms of the physical parameters, as shown below [19], The last equation relates the two CPV angles, α c and α b , and indicates that there exists only one independent CPV phase in our model. Using Eq. (9) and the minimization condition Eq. (4) we obtain the full relationships between model parameters (λ 1 , λ 2 , λ 3 , λ 4 , Reλ 5 , Imλ 5 , m 2 11 , m 2 22 , Rem 2 12 , Imm 2 12 ) and phenomenological parameters (v, tan β, ν, α, α b , α c , m h1 , m h2 , m h3 , m H + ). Through Eq. (18), one can solve for the angle α b in terms of α c , Conversely, one could obtain the formula for α c in terms of α b . However, two solutions will be generated when solving the second order equation for tan α c . Here we adopt the convention in Ref. [24], where sin α max b sets a theoretical bound on the CPV angle α b which comes from the requirement of the existence of a real solution for tan α c : To eliminate the tree level FCNCs, one can assign Z 2 charges to different fermion fields. In general, this would lead to four possible arrangements in the Yukawa sector, which are often dubbed Type-I, Type-II, Lepton-specific and Flipped 2HDMs [30,47,48]. In this work, we only concentrate on the first two, since Type-I (Type-II) differs from Lepton-specific (Flipped) only in the lepton sector and they should behave similarly to the first two in our collider and EDMs experiments. Under the Z 2 symmetry fermion fields transform as The corresponding Yukawa interactions invariant under the Z 2 symmetry are: The interaction of the physical Higgs with fermions and with vector bosons can be parametrized as where c f,i (c f,i ) represents the scalar (pseudo-scalar) component of the physical Higgs h i coupling to fermions while a i stands for the coefficient of h i coupling to the vector bosons. Analytic expressions for these coefficients are given in terms of the phenomenological parameters in Table I. Higgs global fits to the CP conserving 2HDM from current LHC measurements indicate that the couplings are close to the alignment limit: β − α = π/2 [31,32].
Hence, we concentrate on the region having only small deviations from this limit in our study. The interaction between the heavy Higgses, SM Higgs and Z bosons can be parametrized in the following form: with the coefficient g iz1 expressed as: We parametrize the deviation from the alignment limit by a small variable θ where β − α = π/2 + θ. Then we expand coupling g iz1 in the limits of small α b (CPV angle) and θ , which gives, Thus, near the alignment limit, the decay h 2 → Zh 1 could occur only if α b = 0, assuming it is kinematically allowed. In contrast, the decay h 3 → Zh 1 could arise even in the α b = 0 limit so long as there exists a departure from exact alignment. Consequently, one may interpret null results of any search for a heavy scalar decaying to a Z-boson and a SM-like Higgs boson in terms of constraints on either α b or θ. In what follows we will, thus, consider the present and prospective constraints on α b in two cases: θ = 0 and θ = 0.

A. Production of Heavy Higgs
At the LHC, the dominant production mode for a heavy Higgs is via gluon fusion. Therefore, we restrict our study on this specific production mode. The one loop gluon fusion production cross-section of a heavy Higgs is obtained by rescaling the value of the production cross-section for the SM-like Higgs: with τ i f = m 2 hi /(4m 2 f ), the ratio of the mass squared of the heavy Higgs (h i ) to 4 times the mass squared of the fermion running in the loop. Here, σ(gg → H SM ) represents the gluon fusion production cross-section of a heavy Higgs with SM couplings. The functions F H 1/2 and F A 1/2 are defined in the following: As one can see from Eq. (31), the numerator involves the sum of two contributions arising from the CP-odd and CP-even components of the physical Higgs boson, respectively. Denoting M gg→hi CP−odd and M * gg→hi CP−even as the CP-odd and CP-even parts of the gluon fusion matrix elements, we see that the interference term M gg→hi CP−odd M * gg→hi CP−even vanishes after integrating over final state phase space due to parity. The heavy Higgs production cross-section in this form automatically takes into account the K-factor, if one uses the production σ(gg → H SM ) with higher order corrections.
Here we obtain the values of σ(gg → H SM ) from the website [33].

B. Decay of Heavy Higgs
The dominant two body decay modes of the heavy Higgses are taken into account with Γ tot expressed in the following form: where the "+ · · · denote the tiny decay rates to a pair of light fermions and photons, and Z boson and photon, which we have neglected. In addition, we ignore decay rate of a heavy Higgs to one SM-like Higgs and another heavy Higgs, as well as a pair of heavy Higgses because they are forbidden by kinematics due to the mass hierarchy we choose in our benchmark model. The analytical expression for each two-body decay rate can be found in the Appendix B.

IV. SIMULATION DETAIL
In this section, we will discuss details of our collider simulation. We first reproduce the result of 8 TeV ATLAS exclusion limit on σ(gg → h i ) × Br(h i → Zh 1 ) × Br(h 1 → bb) obtained by searching for a heavy Higgs h i=2,3 decaying to Z( + − )h 1 (bb) [28] (As in Ref. [28] we do not include a Br(Z → + − ) factor because it is assumed to have the SM value) . We then use a BDT method to perform events classification and derive the projected exclusion limit for a future 14 TeV search. Events are generated by MadGraph 5 aMC@NLO [34] and then passed through PYTHIA6 [35] for parton showering. Finally Delphes3 [36] is used for fast detector simulation.

A. 8 TeV Result Reproduction
We use the cuts described in Ref. [28] as follow: • The events must have 2 electrons or 2 opposite charged muons with p e,µ T > 7 GeV and |η e |(|η µ |) < 2.5(2.7) • The leptons must have p e,µ T,lead > 25 GeV, and if the leptons are µ + µ − pairs, then one of the µ must satisfy |η µ | < 2.5 • The events must have exactly 2 tagged b-jets with p lead b,T > 45 GeV and p sub b,T > 20 GeV • The reconstructed invariant mass for dilepton and dijet systems should satisfy: 83 < m < 99 GeV and 95 < m bb < 135 GeV where H T is defined as the scalar sum of all jets and leptons in the events 3 is the reconstructed mass of heavy Higgs. For the detector simulation, we use the default Delphes ATLAS cards with following modifications. The values are modified to be consistent with those used in the ATLAS analysis [28]: • The isolation conditions for leptons: Change DeltaRMax from 0.5(default) to 0.2; Change PTMin from 0.5(default) to 0.4(1) for electron(muon); Change PTRatio from 0.1(default) to 0.15. These changes will increase the lepton identification in the boosted regime.
For the signal process, we only take into account the gluon fusion production mode of the heavy Higgs. As for the background processes, we consider the two major backgrounds Zbb and tt as well as to sub-leading backgrounds SM Zh and diboson ZZ backgrounds. For all the backgrounds, we generate events with one additional jet with jet matching. The numbers of events generated and the corresponding acceptance times efficiency are given in Table II. The cross-sections are normalized to the values with higher order corrections. The K-factors for Zbb , tt, Zh, ZZ are calculated based on the result in Ref. [37][38][39][40]. One can observe that the Z( )bb background is a bit larger than the ATLAS result in Table II. This maybe due to the fact that ATLAS used a data-driven method to estimate the number of Z( )bb background events, which may include some effects that our fast detector simulation cannot fully replicate. However, one can also see that these kinds of effects are at a controllable level; our simulation result agrees with ATLAS results within at most 20% uncertainty. Since we may also expect the same kind of effect in 14 TeV simulations, our projected exclusion limit result will be conservative.

Backgrounds/
Signal σ(pb) σ × L simulated # of events after cuts # of expected event in Ref. [28] A ×  [28]. The third column is the total number of events produced at 8TeV LHC with the integrated luminosity equal to 20.3 fb −1 . The fourth column is the number of events left for each background after all the cuts with the integrated luminosity equal to 20.3 fb −1 . The fifth column gives the number of events left with the same cuts estimated by ALTAS in Ref. [28]. The last column gives the acceptance times the efficiency after all the cuts obtained by our simulation.
We present the reconstructed invariant mass of the heavy Higgs in Fig. 1(a) which can be compared with the ATLAS result in Fig.3(b) in Ref. [28]. With this binned distribution we use a profile likelihood method as used in the ATLAS paper to reproduce the 95% CLs exclusion limit. A comparison with ATLAS result is given in Fig. 1(b), the red curve is our reproduced exclusion limit, and the blue curve is the ratio of the ATLAS results to our reproduced values. One can see that, the ratio is generally less than one which corresponds to the excess of Z( )bb background in our simulation. The peak at 800 GeV is due to the lack of background statistics and downward fluctuation near m hi = 800 GeV, but for our benchmark model where m h1 = 550, m h2 = 600 GeV, the ratio seems reasonably close to one.  1: Fig.(a) shows the reconstructed invariant mass distributions with + − bb final state. The signal comes from a heavy Higgs of mass 500 GeV and production cross-section 0.03 pb with an integrated luminosity 20.3 fb −1 . Fig.(b) demonstrates the 95% exclusion limit on the signal σ(gg → A)Br(A → Zh)Br(h → bb). The red curve is our result using the distribution in Fig.(a) with profile likelihood method while the blue curve is the ratio of the ATLAS result (in Fig.3(b) of Ref. [28]) to our reproduced expected exclusion limit.

B. 14 TeV Prediction
We use the same Delphes card when generating events for the 14 TeV case. The preselection cuts we use are almost the same as those for the 8 TeV case. In order to get a sufficiently large sample for BDT analysis, we expand the mass window for m bb from 95 ∼ 135 GeV to 60 ∼ 140 GeV. Also, rather than implementing the E miss T / √ H T and p Z T cuts, we allow the BDT to optimize them. The numbers of events generated and the acceptance times efficiency after preselection for signal and backgrounds are given in Table III   ROOT, Toolkit for Multivariate Data Analysis (TMVA) [41] and the BDT method for the classification of signal and background events. The variables used for the classification are listed below: where p lead,sub T,(j, ) represent the leading and subleading p T of leptons and jets; m and m bb are the invariant masses of dijet and dilepton systems, respectively; p h,Z T stands for the reconstructed p T of the Z boson and the SM Higgs; E miss T / √ H T is the ratio of the missing transverse energy to √ H T defined in the previous subsection; ∆R ,bb,Zh are the angular separations of two leptons, two bjets and reconstructed Zh, respectively, with ∆R ab = (η a − η b ) 2 + (φ a − φ b ) 2 . ∆φ Zh is the separation of the azimuthal angles between Z and h. The distributions of these variables are shown in Figures 9 to 12 in Appendix A.
We select a representative point with M h2 = 550 and M h3 = 600 GeV as the signal to train the BDT. The BDT The distributions of the BDT ouput for two heavy Higgs masses are shown in figs. 2(a) and 2(b). One could find that the discriminating power is a bit better for the heavier Higgs as we expected. The next step is to select a cut on the BDT output to obtain the most stringent 95% exclusion limits. After applying the BDT cuts shown in Table. IV and V, we use the reconstructed heavy Higgs mass distribution of the remaining events to derive the 95% exclusion limit on σ(gg → h 2,3 )Br(h 2,3 → Zh 1 )Br(h 1 → bb). We show the resulting prospective exclusion limits in Tables IV and V. We also perform a cut-based analysis with the same ATLAS cuts described in IV A, and the results are shown in the "cut-based result" column in Tables IV and V. One can see that the exclusion limits of our BDT analysis are 30% to 50% better (lower) than the cut-based analysis results.

V. RESULTS AND DISCUSSION
We now translate our simulated exclusion limit into constraints on the parameter space of CPV 2HDMs. We use a benchmark point below which is consistent with the electroweak precision measurements and muon g − 2 data as discussed in Ref. [24]: m h2 = 550 GeV, m h3 = 600 GeV, m H ± = 620 GeV, ν = 1 (37) and show the constraints on sin α b vs tan β plane. The 95% CLs exclusion limit is given by: where σ L is the exclusion limit listed in Tables IV and V. We assume that the resonance process is dominant when the invariant mass of two gluons is approaching the mass of the heavy Higgs. This is not always true in the parameter space we consider, especially in the limit of small θ. The gluon fusion to Zh 1 box diagram may become important and interfere with the resonant triangle diagram. This may change the distribution of the invariant mass of Zh 1 . Here, we will simply identify the region of parameter space that may suffer from this effect, leaving a detailed analysis for future study. To proceed, we will compare the relative scale of the amplitude squared of the resonant and non-resonant gg → Zh 1 processes at the center of mass energy √ s = m h2 , and m h3 . For the resonance contribution we use the following formula where G F and α s are the Fermi constant and strong coupling constant, respectively, while g iz1 , F We now consider the prospective future reach at the LHC. In the alignment limit, the sensitivity comes primarily from the resonant production of h 2 , as expected from Eqs. 29 and 30. We show the prospective exclusion regions associated with h 2,3 production separately in Fig. 3. The pink region is forbidden by the requirement of the electroweak symmetry breaking. The green, blue, and magenta regions represent the prospective exclusion limits for the LHC integrated luminosities equal to 100 fb −1 , 300 fb −1 and 3000 fb −1 , respectively. The black contours correspond to log 10 |M | 2 i in Eq. (39) with s = m h2,3 for i = 2, 3. If we require |M | 2 i > 10 −4 to guarantee the dominance of the resonant production, then there will be some parts of the prospective exclusion region for 3000 fb −1 at low tan β that may not be valid for our analysis. One could observe that from Fig. 3(b) and Fig. 3(d) there is a loss of sensitivity for h 3 near tan β ≈ 1 in the alignment limit for both Type-I and Type-II models. This is due to the cancellation effect in the coupling g 3z1 when −α = β = π/4. The situation is similar for non-vanishing but small θ. An illustration for cos(β − α) = 0.02 is shown in Fig. 4. However, for a large deviation, for example cos(β − α) = 0.05 in the Type-II model and cos(β − α) = 0.1 in the Type-I model, the constraints from the resonance production of h 3 become important, and can cover large part of the parameter space. In this situation, the effect of the non-resonant production is negligible due to a relatively large |M | 2 3 . This can be seen in Fig. 5. In addition, we take into account the constraint from the Higgs coupling measurements at 7 and 8 TeV LHC as what was done in Ref. [24]. The channels included in the χ 2 analysis are: h 1 → W W , ZZ, γγ, bb, τ τ . We also include the present and prospective constraints from EDM searches given in Ref. [19], which are summarized in Table VI. We find that the constraints from LHC and low energy experiments are complementary. Figures 6 and 7 demonstrate the exclusion limits from both LHC and EDM searches. In each figure, the orange region gives the current LHC exclusion limit. The blue, and magenta regions represent prospective future LHC limits for integrated luminosities equal to 300 fb −1 and 3000 fb −1 , respectively. The light green and light blue regions are excluded by the neutron EDM and electron EDM searches, respectively. The light red and light yellow represents current constraints from the mercury and prospective radium atomic EDM searches. The gray regions are excluded by the Higgs coupling measurements. The pink region is again as described above theory-inaccessible. There are also constraints that we do not show in Fig. 6 and Fig. 7 from heavy flavor physics [43], which exclude the regions of parameter space with tan β less than 0.9 in both Type-I and Type-II models for the benchmark point we choose. These constraints can be relaxed if other new particles are introdueced in addition to the 2HDMs or some non-trivial flavor structure [44,45]. Figure 3 of Ref. [43] demonstrates the constraints on the tan β vs m H + plane for the Type-I and Type-II 2HDMs. The most stringent bounds on tan β come from B s −B s mixing and B 0 s → µ + µ − for the Type-I model, and from B s −B s and B d −B d mixing for the Type-II model.
In Figs. 6(a) and 6(b), we show the current and prospective exclusion regions for the Type-I model in the alignment limit. One can see that the reach of the collider search is not competitive with that of the electron EDM search even at the end of the HL-LHC phase, especially in the low tan β region. This is due to the fact that the collider search is sensitive to Br(h 2,3 → Zh 1 ) in addition to the h 2,3 production cross sections. In the alignment limit, the Zh 1 channel is fed mainly from the decay of h 2 , and the coupling g 2z1 is suppressed by the CPV angle α b as shown in Eq. (29). Moreover, for low tan β, the couplings of h 2 to quarks are enhanced, which leads to a suppression on Br(h 2 → Zh 1 ) and an increasing gluon fusion h 2 production cross section. However, the increase of the latter cannot compensate for the decreasing Br(h 2 → Zh 1 ). The net effect is a reduced Zh 1 signal strength. In contrast, the electron EDM is sensitive to the pseudo-scalar couplings that are enhanced at low tan β in the Type-I model, so electron EDM searches exclude a large part of parameter space in the low tan β region in Fig. 6(b). In Figs. 6(c) to 6(f), we present the current and prospective exclusion regions away from the alignment limit. We can observe that the current collider constraints are not as strong as those from EDMs . However, the future LHC reach can be comparable to that of the EDMs searches, and even better at moderate tan β. One can observe this feature from Eq. (30), where g 3z1 is proportional to θ which describes the level of deviation from the alignment limit, in contrast to Eq. (29) where g 2z1 is suppressed by the small CPV angle α b . One can also observe that in the large tan β region the LHC loses sensitivity. The reason is that at large tan β, pseudo-scalar couplings of both t and b quarks to h 3 are suppressed, so the total partonic production cross-sectionσ(gg → h 3 ) decreases as tan β increases. Despite the possible increase in Br(h 3 → Zh 1 ), the over all effect is a decreasing trend of signal rate towards large tan β resulting in an untestable region for the LHC search.
In the Type-II model, the results are shown in Fig. 7. In contrast to the Type-I model, the electron and mercury EDMs are not able to probe the parameter space when tan β is close to one due to the cancellation in Barr-Zee diagrams indicated in [19] and [23], whereas the neutron and radium EDMs retain sensitivity in this region. In the situation that is close to the alignment limit, as one can observe from Fig. 7(a) to 7(d), the future LHC reach can help to test the region where tan β is close to one. However, the reach of future neutron and radium EDM constraints still exceeds that of the LHC. When the deviation from the alignment limit is as large as cos(β − α) = 0.05, the future LHC may probe a large portion of the parameter space for reasons similar to those for the Type-I model: g 3z1 is sensitive to this deviation and is not suppressed by the CPV angle. Moreover, some portions of the large tan β region cannot be accessed, but for reasons different from the Type-I case. In the Type-II model, the pseudo-scalar coupling of the t-quark to h 3 is suppressed at large tan β, while the pseudo-scalar coupling of the b-quark to h 3 is enhanced. However, for the range of tan β we are interested in, the enhancement of the b-quark loop contribution tô σ(gg → h 3 ) cannot compensate for the suppression of the t-quark loop effect as one can see from Fig. 8(b). Thus, σ(gg → h 3 ) decreases with increasing tan β. As for Br(h 3 → Zh 1 ), due to the increasing Br(h 3 →bb) and decreasing Br(h 3 →tt), the overall effect leads to a decreasing Br(h 3 → Zh 1 ). A decreasing production cross-section combined with a decreasing decay branching ratio makes the large tan β region relatively inaccessible for the LHC in the Type-II model. Now we argue that the future result of EDM and LHC experiments are expected to be complementary to each other and combining information from two kinds of experiments would help us better determine if the 2HDMs are realized in the nature. Since the global fit of the Higgs coupling measurements constrains 2HDMs in the parameter space that are close to the alignment limit, we summarize our results in two categories: 2HDMs are in the exact alignment limit (cos(β − α) = 0), 2HDMs deviate from the alignment limit (i.e. cos(β − α) = 0).
• 2HDMs in the alignment limit: -Future LHC makes a discovery As discussed above, in the alignment limit, the productions of the heavy Higgses h 2 and h 3 are purely determined by the size of CPV angle α b , so the reach of future LHC is merely sensitive to the CPV effect in the model. From Fig. 6(b), one can observe that, in the Type-I model, the reach of future LHC is entirely inside the reach of future radium and electron EDM experiments. Thus, one can conclude that if Type-I model is true, with a discovery at the future LHC one should also observe non-zero radium and electron EDMs, otherwise the null results of radium and electron EDMs will veto the Type-I CPV 2HDM. A similar conclusion can be drawn for the Type-II model by observing Fig. 7(b), where one can find that the LHC sensitive region is well within the reaches of radium and neutron EDMs. Hence, if the Type-II model is true, the discovery of the future LHC should lead to the observations of non-zero radium and neutron EDMs. The left (right) column is for the current (future) exclusion limit. The orange region is excluded by the current LHC data. The blue and magenta regions represent the future LHC limit with integrated luminosities equal to 300 fb −1 and 3000 fb −1 , respectively. Light transparent red represents the constraint from mercury EDM, light blue denotes electron EDM, light transparent green stands for neutron EDM, and light yellow signifies future prospective radium EDM. The gray region is excluded by the coupling measurement of the SM-like Higgs and the pink region is theoretically inaccessible due to the absence of a real solution for αc. The benchmark point used here is m h 2 = 550 GeV, m h 3 = 600 GeV, m H ± = 620 GeV, ν = 1.
-Future LHC gives a null result For both Type-I and Type-II models, a null result from future LHC does not exclude the possibility of CPV in 2HDM as long as CPV angle α b is sufficiently small. Meanwhile, any non-zero EDMs that correspond to the regions of parameter space within the reach of LHC would disfavor the CPV 2HDMs.
• 2HDMs away from the alignment limit: The left (right) column is for the current (future) exclusion limit. The orange region is excluded by the current LHC data. The blue and magenta regions represent the future LHC limit with integrated luminosities equal to 300 fb −1 and 3000 fb −1 , respectively. Light transparent red represents the constraint from mercury EDM, light blue denotes electron EDM, light transparent green stands for neutron EDM, and light yellow signifies future prospective radium EDM. The gray region is excluded by the coupling measurement of the SM-like Higgs and the pink region is theoretically inaccessible due to the absence of a real solution for αc. The benchmark point used here is m h 2 = 550 GeV, m h 3 = 600 GeV, m H ± = 620 GeV, ν = 1.
-Future LHC makes a discovery In this case, the result from future LHC is not purely sensitive to the CPV effect since the coupling of h 3 to Zh 1 is proportional to the level of deviation from the alignment limit and is not suppressed by the CPV angle α b . A discovery at the future LHC may or may not imply a non-zero EDM result, a situation that depends largely on the magnitude of deviation from the alignment limit. Moreover, if the deviation is relatively small, such as cos(β − α) = 0.02, the exclusion limit would mainly come from h 2 , for which production is purely sensitive to the CPV angle α b , as shown in Fig. 6(d) and Fig. 7(d) for the Type-I and Type-II models, respectively. In addition, one would expect the non-zero radium and electron EDMs for the Type-I model and non-zero radium and neutron EDMs for the Type-II model. Thus, a null EDM result would disfavor both Type-I and Type-II CPV 2HDMs. On the other hand, if the deviation is relatively large, such as cos(β − α) = 0.1 in the Type-I model ( Fig. 6(f)) and cos(β − α) = 0.05 in the Type-II model ( Fig. 7(f)), the exclusion power would be dominated by the h 3 decay. As mentioned above, the discovery of h 3 does not necessarily lead to sizable EDMs, and therefore a more detailed study of the CP properties of the newly discovered particle would be needed.
-Future LHC gives a null result In this case, for sufficiently large deviations as shown in Fig. 6(f) and Fig. 7(f), the discovery of any EDM results would indicate that the CPV source is not consistent with the CPV 2HDMs. On the other hand, for relatively a small deviation as shown in Fig. 6(d) and Fig. 7(d), the CPV 2HDMs is still available if CPV angle α b is sufficiently small. Moreover, any non-zero EDMs that correspond to the regions of parameter space within the reach of LHC would disfavor the CPV 2HDMs.
Finally, we comment on the potential constraints from the viability of successful EWBG in 2HDMs. One can potentially include the allowed regions where EWBG is viable in Figs. 6 and 7. For example, the authors of Ref. [27] studied the CPV for EWBG and identified some regions of parameter space that seem favorable. They pointed out that the CP violating phase necessary for successful baryogenesis is sensitive to tan β and the masses of the heavy Higgses. That work also concentrated on parameter space region where the dominant decay channel of h 3 is to the Zh 2 final state (A → ZH in the CP-conserving limit), based on earlier studies of the electroweak phase transition [10,49]. A strong first order electroweak phase transition favors -but does not absolutely require -regions of parameter space leading to dominance of this decay mode. In this paper our main focus is on constraints of the CP violating phases from the LHC and EDMs in 2HDMs. It is possible that for the spectra considered here, the CPV 2HDMs can accommodate a strong first order electroweak phase transition and give rise to the CPV asymmetries needed for successful baryogenesis. A detailed and general analysis of this possibility is beyond the scope of the present paper and will be left for future study.

VI. CONCLUSIONS
The CP properties of the SM-like Higgs boson provide a portal to investigate the possibility of CP violation beyond the Standard Model, an important ingredient for successful baryogenesis. A CP violating 2HDM with an approximate Z 2 symmetry may provide this source of CP-violation. Future searches for new scalars arising in the 2HDM, together with next generation EDM searches, may both discover the 2HDM and determine whether or not it contains the CPV interactions necessary for electroweak baryogenesis. In this paper, we have studied how well a future LHC search for a heavy Higgs with Zh 1 decay mode and Z( )h 1 (bb) final state can probe the parameter space of CP violating 2HDMs. We used the BDT method to estimate the future exclusion limits on the observable σ(gg → h 2,3 ) × Br(h 2,3 → Zh 1 ) × Br(h 1 → bb). We then compared this reach with that of next generation searches for the EDMs of the electron, neutron and neutral atoms.
Our results in Figs. 6 and 7, lead to the following conclusions: (1) In the exact alignment limit, a discovery at the LHC would imply observable radium and electron EDMs for the Type-I model and observable radium and neutron EDMs for the Type-II model. Null LHC results would still be consistent with CPV 2HDMs if the CPV angle is sufficiently small, but would not preclude observable EDMs. (2) Away from the alignment limit, a discovery at the LHC may or may not imply non-zero EDM results, depending on the level of deviation from the alignment limit. If the deviation is small, the LHC reach is mainly dominated by the production of the mostly CP-even Higgs (h 2 ). In this case, one can reach the same conclusion as in the alignment limit for both Type-I and Type-II models. However, with a relatively larger deviation, an LHC discovery may not imply non-zero EDM results because the production of the mostly CP odd Higgs (h 3 ) is not purely sensitive to the variation of the CPV phases. In addition, the LHC reach would cover most of the parameter space that can also be probed by the EDM searches. As a consequence, a null LHC result, together with observation of the EDMs, would disfavor the CPV 2HDMs. Finally, we also point out that our analysis may break down in the small tan β and α b regions, due to the non-negligible interference effect between the resonant and non-resonant productions of Zh 1 through the gluon fusion channel. We leave the detailed study of this effect for future work.