Restricting the parameter space of Type-II Two Higgs Doublet Models with CP violation

We explore the parameter space of the type-II two-Higgs-doublet model with softly broken Z2 symmetry, allowing for CP violation in the scalar potential. Imposing theory-motivated and experimentally-driven constraints, we show that as the CP-violating phases are increased, only small regions of parameter space survive, including regions slightly away from the alignment limit. Electroweak oblique parameters and electric dipole moments emerge as most restrictive constraints. We show that imposing these constraints (as well as theoretical bounds such as perturbativity) the masses of the charged and heavy neutral Higgs, unlike in the CP conserving case, are bound from above and below. In particular, the heavy neutral Higgs masses are almost degenerate. In this parameter space region we highlight the relevant decay signals of the heavy neutral Higgs, involving both Zh and WW=ZZ, indicative of CP violation in the model.


I. INTRODUCTION AND MOTIVATION
Our present day knowledge of particle physics, which seeks to address the fundamental building blocks of matter and how they interact, is encapsulated in the Standard Model (SM). The Standard Model is a relativistic quantum field theory based on SU (3) C ⊗ SU (2) L ⊗ U (1) Y gauge symmetry, where the group SU (2) L ⊗ U (1) Y is responsible for electroweak unification, and SU (3) C provides the description of strong interactions [1][2][3][4][5]. In 2012, the LHC discovered a spin 0 scalar boson with mass about 125 -126 GeV [6,7], and with properties in remarkable agreement with the Higgs boson predicted by the SM, was especially important, as it confirmed experimentally one of the key mechanisms of the SM, responsible for giving masses to the SM particles. The dynamics of this mechanism, known as electroweak symmetry breaking (EWSB), can have a significant effect on the properties of SM particles. Unfortunately, despite its remarkable success, the SM cannot be the ultimate theory. Experimentally, it lacks a dark matter candidate, a mechanism for neutrino masses and mixing (neutrino oscillations) or for sufficient generation of matter-antimatter asymmetry, to name a few. From a theoretical point of view, despite being the most successful model to date, the SM is inherently an incomplete theory [8]. Thus a lot of effort from phenomenologists and theorists has been invested in models beyond the SM (now interpreted as a low-energy effective theory), all of which call for extending the particle content of the SM and sometimes, of the gauge symmetry of the model.
Extension of the particle content includes, in its simplest realization, extending the Higgs boson content. One of the simplest, and most widely studied extension of the SM Higgs sector is the Two Higgs doublet model (2HDM), where a second scalar Higgs doublet is added, with the same quantum numbers as the existing one [9,10]. There are several versions of the model, depending on the choice of Yukawa couplings of the two doublets. One of these 2HDM is adopted by supersymmetry, where up-and down-type quarks cannot acquire masses by coupling to the same Higgs doublet, see for example [11].
With the addition of the new doublet, the model could provide a new source of charge-parity (CP) violation ( [10,[12][13][14]) which is fundamental to explaining matter-antimatter asymmetry of the Universe [15,16].
A significant difference between the CP conserving (CPC) and the CP violating (CPV) 2HDMs is that, in a CPC 2HDM, one of the three scalar physical states, identified as the CP-odd Higgs boson remains decoupled to the other neutral scalars, whereas in the latter all the three neutral Higgs states mix. One of these states is identified with the 125 GeV Higgs boson, and all three scalar states interact with the gauge bosons. Many detailed studies of CPC 2HDMs have been presented, for example in [10,17,18], and other works such as [19][20][21][22] have looked at CPV 2HDMs using CP-odd weak-basis invariants. Also, recent studies of charged Higgs bosons phenomenology in the context of the CPV 2HDMs have been presented in [23][24][25][26][27] and a CPV 2HDM analysis of viable parameter space regions which pass experimental constraints has been given in [28][29][30][31][32].
The 2HDM Yukawa Lagrangian entails additional features with respect to the SM Lagrangian, which are strongly constrained by experimental data. The main difference is that, in the most general 2HDM Yukawa Lagrangian, flavour a Email: mariana.frank@concordia.ca b Email: ericgyabeng2012@gmail.com c Email: mtoharia@dawsoncollege.qc.ca arXiv:2112.14295v2 [hep-ph] 12 Aug 2022 changing neutral currents (FCNC) arise at the tree-level [33,34], and these are strongly constrained by experimental data. One can eliminate these dangerous FCNCs by imposing a Z 2 symmetry on the scalar potential and assigning Z 2 charges to the fermions. Various versions of 2HDMs were studied, which can be classified into four different classes based on which type of fermions couple to a specific doublet. These classes are usually grouped as Type I, Type II, X and Y [10,[35][36][37] 1 .
In this paper we perform a detailed study of the CPV 2HDMs with Type-II Yukawa couplings. In this model, the couplings of the down-type quarks and charged leptons are proportional to tan β, the ratio of the two neutral Higgs bosons vacuum expectation values (VEVs), while the top quark couplings are proportional to cot β, and thus the model is significantly constrained by flavor physics and direct searches of extra Higgs bosons.
Due to its connection with tree-level MSSM, the Type-II 2HDM is a popular extension of the Higgs sector [9,38]. While at tree-level, the MSSM Higgs sector is strictly CP-conserving, it has been shown that under certain circumstances, CP-violating effective Higgs sectors could arise via loop corrections [13,[39][40][41][42][43]. CPV 2HDMs effects in the Higgs sector have been studied extensively in the MSSM limit of the 2HDM in, for example, [44][45][46][47][48].
Here we analyse the consistency of the Type-II CPV 2HDM with Higgs data, electroweak precision and perturbativity and unitarity constraints, as well as CPV in electric dipole moments. Our aim is to investigate in detail the allowed parameter space in the light of the most recent LHC results, in particular the constraints imposed by LHC-Higgs data. In order to study the impact of CP violation on the Higgs bosons and their respective decay channels, we consider the CPV 2HDM potential parameters [32] in the approximate Z 2 case [49]. We study the various real and reduced couplings in a basis containing the physical states with CP violation turned on by fixing the real and imaginary couplings in the scalar potential. The size of these couplings is constrained by perturbativity, unitarity, as well as by electroweak precision parameters, and the imaginary part of these couplings will be further constrained by Electric Dipole Moment (EDM) measurements.
As we increase the amount of CPV by turning on the phases in the Higgs potential, we show that constraints from LHC Higgs data still allow some regions of parameter space slightly away from the alignment limit (the region where the SM Higgs boson decouples from the rest of the scalar states). In those viable regions, the masses of the heavy neutral Higgs bosons are tightly constrained by the theoretical and experimental bounds. These constraints prove to be very stringent also in regards to the rest of parameters such as tan β, or the mixing angle of the Higgs sector α. Production and decay rates of the heavy neutral Higgs bosons, constrained from these considerations, at the LHC are presented and CPV specific signals are highlighted.
Our paper is organized as follows: in Section II we present the Type-II Two Higgs Doublet Model with CP violation, concentrating on the parameters in the scalar potential and their effects on the neutral Higgs states. In Section III we explore the bounds on the parameters in the context of both theoretical and experimental constraints, from perturbativity and unitarity of the potential, the Higgs data, B-decays, as well as the oblique parameters and the EDMs. In Section IV we present the main results, showing the restrictions on the parameter space on the charged and neutral Higgs masses, and in Section V we analyze their production cross sections and plot their decay patterns. We summarize and conclude in Section VI.

II. TWO HIGGS DOUBLET MODEL WITH CP VIOLATION
In this section we review the 2HDM considered in this study. We first introduce the most general Two Higgs doublet scalar potential which breaks SU (2) L × U (1) to U (1) EM , where m 11 , m 22 and λ 1 , λ 2 , λ 3 , λ 4 are real parameters, while m 2 12 , λ 6 and λ 7 are complex, and where the Higgs doublets are parametrized as To avoid tree-level flavor changing neutral currents, one imposes a Z 2 symmetry with Higgs bosons obeying the transformation properties, Eq. (2.3) implies λ 6 = λ 7 = 0 2 , while we allow a non-zero m 12 to break softly the Z 2 symmetry of Eq. (2.1). After electroweak symmetry breaking, the Higgs doublets can be written in terms of the physical states as, GeV and H + is the physical charged Higgs boson with mass M H + .
In the neutral scalar sector, there is mixing of the three states, h 0 1 , h 0 2 and A 0 due to two independent angles and phases from the the imaginary parts of m 12 and λ 5 . The mixing among the three neutral scalars can be parametrized by an orthogonal matrix R [32], with the convention: This form of parametrization is transparent: α mixes the lightest Higgs with the heavier Higgs (as in the CP-conserving case), α b mixes the lightest Higgs with the CP-odd component, while α c mixes the heavier Higgs with the CP-odd component. The physical mass eigenstates are then defined as (h, with h assumed to the the SM-like boson with M h = 125 GeV. This diagonalization procedure yields seven linearly independent equations that link the physical parameters to the parameters in the scalar potential (see for example [30,32]).
for (α b = 0). (2.14) In the CP conserving version of the 2HDM, α b = α c = 0, R is block diagonal, and h and H 2 have no pseudoscalar component.

II.1. CP Violation and Mixing in the Neutral Higgs Sector
CP violation can occur in the 2HDM in the scalar sector, and the effects can be considerable. Setting λ 6 = λ 7 = 0 terms, to avoid hard breaking of Z 2 symmetry, we allow it to be only softly broken by the m 2 12 term, and using the tadpole equations (obtained by minimizing the Higgs potential Eq. 2.1), we are left with 8 independent physical parameters • Three scalar masses, M h , M H2 , and M H ± ; • Three mixing angle α, α b , α c in the neutral scalar sector (one from the CP-conserving sector, two from allowing CP violation); • tan β = v u /v d , the ratio of VEVs; • (m 2 12 ), or ν ≡ (m 12 ) 2 /(v 2 sin 2β). Here M H2,3 are two of the neutral Higgs masses, M H ± is the charged Higgs mass, the neutral-sector input mixing angles are α, α b , α c , and (m 2 12 ) is the real part of the soft Z 2 breaking parameter. We fix M h = 125GeV . As well the third neutral scalar mass M H3 is calculated from the two other ones using matrix elements of R [50] , (2.15) which reduces the number of independent parameters to 7. In this model, H 3 is, in the limit of no CPV, the CP odd Higgs, while in the same limit H 2 is the heavy CP-even Higgs. We require the lightest Higgs properties to be consistent with those of the SM-like boson observed at the LHC, and in particular, set its mass to M h = 125.35 ± 0.15 GeV [51]. We proceed by listing the restrictions imposed on the parameter space of the model. The imaginary part of λ 5 , which is a source of CP violation is given in Eq. 2.12. We concentrate on 2HDMs with a Z 2 symmetry in the Yukawa sector and where Φ d and Φ u give mass to only up or down type fermions, respectively. This corresponds to a CPV version of Type-II 2HDMs, leading to suppressed tree-level flavor changing processes mediated by the neutral Higgs scalars. In general, the imposed Z 2 symmetry in the Yukawa sector is not preserved by renormalization, as terms from the Higgs potential will induce couplings of Φ d , Φ u to both up and down type quarks. However, this does not reintroduce any tree level flavor changing effects because the induced Yukawa matrices are still proportional to the corresponding fermion mass matrices. Neglecting mixing from the CKM matrix, the Yukawa Lagrangian is where Q T L = (u L , d L ), L T L = (ν l , L ), u, d and stand for up and down-type quarks and leptons, respectively. The couplings between neutral Higgs bosons and the fermions and gauge bosons are where H i = h, H 2 , H 3 and c Hif f ,c Hif f = 0 or a Hi ,c Hif f = 0 indicates CP-violation, that is, the mass eigenstate H i couples to both CP even and CP odd operators. The coefficients c Hif f ,c Hif f and a Hi are obtained using the matrix R defined above in terms of α, α b , α c and tan β and are given in Table I. Thus we can express the couplings presented in Table I in terms of the mixing angles in the scalar sector of the CPV 2HDM, α (CPC), α b and α c (CPV). They represent the 12, 13 and 23 angles from the 3 × 3 orthogonal diagonalization matrix R given in Eq. 2.5. For the SM-like Higgs boson, the mixing elements are For example, the normalized coupling of top quarks to the SM-like Higgs has the simple expression [52] (in terms of α, β and α b ) 3 : When the Higgs boson mixes with both the CP-even and CP-odd states, the Higgs coupling to vector bosons V is [53], where as above a h measures the departure from the SM, a h = 1 for a pure CP-even state with SM-like couplings and a Hi = 0 for a pure CP-odd state. The effects of CP mixing will appear also at the loop level, especially in the gg → h and h → γγ rates. At leading order, the Higgs production rates normalized to the SM expectations can be written as [53,54] with the loop form factors given by [17]

III. CONSTRAINTS
The CP-violating 2HDM is constrained from various theoretical observations and from measurements at collider and non-collider experiments. We first discuss constraints on model parameters from theoretical considerations, then those from the Higgs data, from B-physics measurements, from precision data (the S, T, U oblique parameters) and from measurements of EDMs. We restrict ourselves to a brief discussion, and refer to explicit expressions that have appeared before.

III.1. Perturbativity, vacuum stability and unitarity
Constraining the theory to be perturbative imposes constraints on all the couplings in the scalar potential In addition, vacuum stability requires the scalar potential, Eq. 2.1, to be be positive for large values of Φ d , Φ u Both unitarity and perturbativity requirements. lead to the constraints [10, 55]: 3 Note that the Yukawa couplings for the neutral Higgs bosons H i with the up-type quarks in the CPV 2HDM are the first term from the CP-even Higgs coupling, and the second from the CP-odd coupling.
As the Eqs. (3.1) -(3.3) only depend on the absolute value of λ 5 , they do not constrain the CPV phases. These conditions are very important for the CP-violating model, as they provide bounds on the masses of the heavy neutral and charged Higgs. From Eqs. 2.10 -2.11, requiring −8π ≤ λ 4 , λ 5 ≤ 8π, we obtain the from-above and from-below bounds for M 2 where we have introduced the mass Hi . If we now consider Eq. 2.14 and require −8π ≤ λ 5 ≤ 8π, we obtain an upper bound forM 2 And sinceM 2 is present in the inequality involving M 2 H ± , we can finally obtain the upper bound Note that this bound comes from requiring perturbativity on the quartic couplings λ 4 , λ 5 and λ 5 and so is not necessarily the strongest bound. Similar upper bounds can be obtained for both (M 2 H2 + M 2 H3 ) and ν. From Eq. 2.11 we obtain and finally, using the previous bounds onM 2 we obtain

III.2. SM Higgs data
The global data set on the SM Higgs boson includes results from LEP, the Tevatron and the more recent LHC experiments. In order to properly combine the constraints coming from the LHC-Higgs data, we used Lilith [53,56].
Lilith is an open-source library written in Python, which can be used in any Python script as well as in C [57] and C++ [58]/ROOT [59] codes, with a command-line interface. At present, Lilith includes the complete data on Higgs measurements at the Tevatron and LHC Run II at 36 fb −1 , and only the Higgs production and diboson decay at 139 fb −1 . 4 The SM-like Higgs signal rates and masses are compared with the signal rate measurements, with the additional requirement that the lightest Higgs has a same mass and production and decay properties as the observed Higgs peak in the channels with high mass resolutions h → ZZ → 4 and h → γγ. The signal strength for the 2HDM versus the SM is defined as where ω i are the experimental weights of the measurement which includes the experimental efficiency.

III.3. B Physics constraints
The charged Higgs bosons in the 2HDM contribute to the decays of B mesons and thus B-physics data set can be used to constrain their masses and couplings. Note that, since the couplings of the charged Higgs bosons are not sensitive to the parameters in the neutral sector, these constraints are independent of the amount of CP violation in the model. The most stringent constraints on model parameters emerge from the B → X s γ [61] 11) and the most restrictive conditions on charged Higgs masses are for tan β = 1, where M H ± ≥ 580 GeV [62].

III.4. S, T, U parameters
This theoretical constraint comes from the oblique parameters, constrained by the global fit [63][64][65] to be The oblique parameters S, T , and U are observables that combine electroweak precision data to quantify deviation from the SM and thus are used in any BSM to ensure that the model is consistent with the data. An alternative way to evaluate the oblique parameters is to fix U = 0, motivated by the fact that U is suppressed by an additional factor M 2 Z /M 2 H compared S and T , which improves the precision on S and particularly on the T parameter, yielding S = 0.00 ± 0.07, T = 0.05 ± 0.06.
The implication of electroweak precision on the mass spectrum of heavier neutral and on the charged Higgs was thoroughly investigated before [66], in the context of CPV 2HDM. It was shown that, except for the T parameter, S and U fall well within the experimental bounds. Thus working in the limit U =0 is more conservative, as it restricts the T parameter more. Allowing U = 0 indicates a strong correlation between T and U parameters, with values of T > 1 corresponding to values of U > 0.01. In general, the 2HDM can generically produce values for T that exceed the experimental bounds, limiting much of the parameter space.
However, imposing mass splitings [32]: insures that the electroweak precision lie in the allowed ranges.
We shall see that our parameter space surviving all constraints is highly degenerate and thus the electroweak corrections are well within experimental bounds. The S, T and U parameters finite and in principle, observable. Their expressions in are found it, e.g., [67], so we do not repeat them here.

III.5. Electric Dipole Moments
The electron EDM measurement places a very strict constraint on the complex Yukawa couplings in most models. In general, the electric dipole moment of a fermion f corresponds to the imaginary parts of the Wilson coefficient d f of the effective operator Complete expressions for the one-loop evaluation of the EDMs can be found in [30,68,69]. In addition to one-loop contributions, EDMs can originate from two-loop Barr Zee contributions [70], which are proportional to a single power of the electron Yukawa coupling [71,72]. These contributions are more important for Type-X (lepton-specific) 2HDM, and still, even in Type-X, agreement with low-energy data for the muon anomalous magnetic moment (and significant EDM contributions) require very light pseudoscalar masses and very large values of tan β [73][74][75][76][77][78], inconsistent with the parameter space analyzed in this model. It has also been shown that, in unlike in the Type-I 2HDM, the electron and mercury EDMs are not able to probe the parameter space when tan β is close to 1, which will be our case, due to the cancellation in Barr-Zee diagrams [79]. We have not considered these contributions here.
A recent electron EDM measurement was performed using the ThO molecule in [80]. The constrained quantity is |d eff e | < 1.1 × 10 −29 e · cm . (3.14) In 2HDM the CP violating phase strongly affects the magnitude of the EDM through the modified couplings of the Higgs bosons. In Type-II 2HDM the EDM is also expected to be enhanced by tan β [61]. By investigating CP violating effects in extended Higgs sector, the electric dipole moment and its effects on produced particles via protons, photons, or electron and positron collisions have all been studied [81][82][83]. We impose these constraints and explore their effects further in our numerical explorations.

IV. ALLOWED PARAMETER SPACE
In this section, we describe the methodology employed to explore the parameter space. We first investigated the parameter space region by performing detailed calculations of the reduced couplings associated with CPV 2HDM and obtained the different signal strengths associated to the SM Higgs. Using Lilith we explored constraints on the model parameters in the cos(β − α) and tan β plane. Once the SM Higgs bounds were considered in this way, we further employed ScannerS [52], a code that performs parameter scans and checks parameter points in BSM theories with extended scalar sectors. ScannerS incorporates theoretical (perturbative unitarity, vacuum stability, boundedness from below) and experimental constraints (electroweak precision, flavour constraints, Higgs searches and measurements, EDMs) from several different sources in order to determine whether a parameter point is allowed or excluded at approximately 95% CL.
The constraints on the cos(β − α) and tan β plane due to LHC-Higgs data are presented in Fig. 1 which shows the allowed parameter space for different choices of the α b angle, a measurement of the amount of CPV admixture within the SM-like Higgs. Note that the composition of the SM-like Higgs, and thus the allowed parameter region plotted, is independent of α c , which affects only the amount of CPV admixture in the other neutral Higgs bosons. The top left-handed panel shows the parameter restrictions for the CP-conserving case α b = 0, where we see that a significant region of the parameter space survives away from the alignment region (for which cos(β − α) = 0). The interesting particularity of this away-from-alignment region is that the bottom quark Yukawa coupling c hbb is negative while the top Yukawa coupling does not flip sign [84,85]. Also note that this plot is consistent with more recent phenomenological studies of 2HDM such as [86]. The parameter cos(β − α) determines the c hV V and c HV V couplings and is a measure of deviations from alignment. As usual, tan β, the ratio of VEVs, is important in determining the relative strength of the couplings c hbb and c htt (and similarly for c H2,3bb and c H2,3tt ). Increasing α b changes significantly the parameter space, in particular, it reduces the allowed values for tan β, which can be large at α In all graphs the regions of the allowed parameter space shrink considerably when increasing α b (as expected) but noticeably, the region away from alignment survives. This is significant because the away-from-alignment region is a harbinger of non-standard behaviour in the Higgs sector. In the figure, the maroon points correspond to 3σ agreement with the data, green points correspond to 2σ and purple points to 1σ agreement. Increasing α b leads to slowly shrinking of the parameter regions consistent with experiment, and for α b = 0.4 almost the whole parameter space is excluded, with only a few points near the alignment region surviving. Still, it is interesting to note that, for α b = 0.1, parameter points consistent with the Higgs data at 2σ survive, in both the alignment and away-from-alignment regions.
As mentioned before, the Lilith code, which is self-contained, combines different Higgs signal constraints from published data at 36 fb −1 , with the addition of the measurement at 139 fb −1 for the process gg → h → ZZ. The combination of all the most recently published signals has not been fully implemented yet in Lilith. Nevertheless it is possible to show the evolution of each individual bound and gain an insight on the type of pressure that each Higgs signal puts on the parameter space considered here. We find that the allowed parameter space is most sensitive to three main Higgs signals, namely the bb, V V , and γγ channels. Sensitivity to the rest of Higgs decay channels is less important as their bounds on the parameter space under consideration are always less restrictive than the mentioned signals.
To display the effect of reduced experimental error bars thanks to improved statistics and analysis, we will focus on three individual studies from the ATLAS collaboration, comparing their older results at 36 fb −1 with their most recent bounds with 139 fb −1 . In Fig. 2 we chose two examples within the same parameter space as Fig. 1, one for α b = 0.1, the other for α b = 0.3. We show contours for the signal strengths of the Higgs channels gg → h → ZZ (red-thick), V h → (h → bb) (green-light) and gg → h → γγ (black-thin), using the bounds set by ATLAS, with 36 fb −1 (dashed/dot-dashed/dotted) and with 139 fb −1 (solid) of integrated luminosity [87][88][89]. In both cases, the white regions are points in parameter space that lie individually within the bounds from bb, ZZ and γγ, using the ATLAS data at 139 fb −1 of luminosity. Note that these regions do not represent a combination of bounds, but are merely representative of the parameter space that will likely still be allowed once the bounds are properly combined 5 . The  improved statistics observed are clearly apparent in the three signal contours with shrinking of the allowed regions as data is improved from 36 fb −1 to 139 fb −1 . This indicates that while our chosen parameter points are not ruled out, they are under more pressure, and may fit the new data with less C.L.
After delineating the allowed and excluded regions of parameter space using the LHC Higgs data, we proceed to investigate the constraints of the CP-violating phases α b and α c on the masses of the heavier neutral Higgs bosons H 2 and H 3 and on the mass of the charged Higgs boson H ± .
In Fig. 3 we show the surviving parameter space in the M H ± − M H2 plane, where H 2 is the lightest of the two non-SM like neutral Higgs bosons. In these plots α b , cos (β − α) and tan β are fixed such that we are within an allowed parameter space point from Fig. 1 The black dots represent the regions allowed by all constraints in Sec. III, except for EDMs and oblique parameters (we use ScannerS). We have chosen a parameter space region allowed by Higgs data (see Fig. 1) and the main bounds on the black dotted region come from B-physics constraints which are responsible for the requirement M H ± ≥ 650 GeV. In addition, perturbativity and unitarity requirements which restrict the parameter space to within the diagonal border-lines at the left, right and top of each plot. In particular one can see that direct unitarity and perturbativity constraints involve the coupling λ 4 , which depends explicitly on both the charged Higgs mass and the neutral Higgs masses (see Eqs. (2.9) and (2.10)). Within the CP-violating parameters considered here, these unitarity and perturbativity requirements make M H ± bound from above, disallowing it to be much larger than M H2 and M H3 , as shown in the discussion in Sec. III.1.
Inside the black dotted region, we show as the rainbow colored region, the physical domain that passes all tests, including electroweak precision measurements constraints from the oblique parameters (ST U ), as in III.4, and CPV constraints from electric dipole moments, as in III.5. The legend at the right shows the deviation from the prescribed ST U and EDM parameters as measured by χ 2 , restricted to be < 7. We find that, within the allowed parameter region of Fig. 1, points that pass ST U and EDM constraints require an almost fixed value for tan β, consistently around tan β ∼ 0.9 and always smaller than 1 6 . The reason is that a relative mass degeneracy among the heavy scalars is necessary to avoid precision tests bounds, and with non-zero CP violation, the masses of the three neutral scalars are all connected via tan β (see Eq. 2.13). One can check that for small values of cos(β − α) (i.e. for β ∼ α + π/2) a small mass splitting between between M H2 and M H3 can only happen when tan β ∼ 1.
The parameter regions chosen in all four graphs correspond to regions slightly away from alignment in which cos(β − α) > 0. In the case of CP conservation, and in the decoupling limit, the mass M h is independent of the other two neutral boson masses (one CP-even, one CP-odd) which are not restricted. This is not the case for the CPV case, and slightly away from the decoupling limit. Here the parameter α c which determines the CP admixture in the heavier neutral bosons H 2 and H 3 is varied but only values within a limited range lead to allowed points, the exact interval depending on the value α b chosen. Increasing α b enhances the amount of CPV admixture in the SM-like Higgs boson resulting in a decrease in the allowed range by EDMs and ST U constraints in the M H ± − M H2 plane. The allowed regions (in rainbow colors) are separated by straight diagonal lines in M H ± − M H2 , and while the allowed window for M H2 is approximately constant, M H2 ∼ 550 → 700 − 1000 GeV. The range for M H ± , while being run into the TeV region, is restricted by unitarity and perturbativity requirements, as explained above. M H ± shrinks considerably with increasing α b , from 650-950 GeV for α b = 0.1 in the top left hand panel, to 650-685 GeV for α b = 0.4, in the bottom right hand panel. We also note that regions slightly more separated from alignment survive for larger values of α b , such as for α b = 0.4, where a small region of parameter space is allowed for cos(β − α) 0.15. In this region, Γ(H 2,3 → V V ) = 0 and hence, the decay of the heavy neutral Higgs into the electroweak gauge bosons is allowed.
In Fig. 4, we show the allowed parameter space in the M H3 versus M H2 plane from imposing all constraints except the EDMs and the ST U parameters (black points), and with EDMS and oblique parameters (rainbow colored points). Here, the mass of the charged Higgs was varied in the allowed ranges from Fig. 3. The relationship between the masses (squared) of the two neutral heavy scalars is linear, and lies in a restricted range. In all cases we find that these two masses must be relatively degenerate in mass in order to satisfy experimental constraints, in particular restrictions from the ST U parameters. The allowed ranges for both masses shrink with increasing α b , but note also that for larger values of α b the degeneracy requirement between the two masses is slightly relaxed.
In the CPV 2HDM model, a main source of CP violation is encapsulated in λ 5 , given in Eq. 2.12. From [30,[90][91][92][93], the unitarity bound on λ 5 should be λ 5 < 4π. The value depends on the CP violating mixing angles α b , α c and on the masses M H2 and M H3 . We show in Fig. 5 below, λ 5 as a function of M H2 . As before the mass of the SM-like Higgs boson h is kept at 125.09 GeV, and we also take the same input parameters as in Fig. 4 above. We show here only the rainbow-colored points, that is, points that survive all constraints, including EDMs and ST U bounds.

V. PRODUCTION AND DECAY OF HEAVY HIGGS BOSONS AT THE LHC
We now proceed to analyze the implications of the restrictions on the parameter space of M H2 , M H3 , M H ± on the production and decay rates of the heavy gauge bosons at the LHC. We analyse regions of parameter space that pass all restrictions, as shown in Figs. 3 and 4.

V.1. Neutral Higgs Bosons Production
The heavier Higgs bosons are expected to be produced at the LHC dominantly via the gluon fusion process gg → H 2,3 . Additional production processes include vector boson fusion pp → H 2,3 qq , vector boson associated production pp → H 2,3 V , and top associated production gg → H 2,3 tt.
In this work we employ SusHi [94] to calculate cross sections for gluon fusion to next to leading order (NLO), At one loop, for the gluon fusion process, the ratio of the heavy Higgs boson production cross section in a CPV 2HDM to that of a SM-like Higgs is where τ i f = M 2 Hi /(4m 2 f ), τ f = M 2 h /(4m 2 f ) and i = 2, 3, f = t, b. The functions A H 1/2 , A A 1/2 are given by The above expressions are given at LO. In our numerical evaluation, the production cross sections via gluon fusion are calculated at NLO using SusHiv1.6.0 [94] and interfaced with ScannerS [52]. As neutral Higgs bosons in our model have no definite CP assignment, CP-even and CP-odd contributions are summed over [50,93]. In addition, all decay channels including the dominant higher-order effects such as radiative corrections and multi-body channels are included as in HDECAY [95,96].

V.2. Neutral heavy Higgs Boson decay rates
The heavy neutral scalar Higgs bosons decay rates into electroweak gauge bosons are where V = W, Z, δ W = 2, δ Z = 1, and i = 2, 3. We note that in the alignment limit, Γ(H 2,3 → V V ) = 0 when sin α b = sin α c = 0, or away from alignment. These channels open up for non-zero CP violation. The decay rates of the neutral Higgs bosons to SM fermions are 6) where N c = 3 for quarks and 1 for charged leptons.
The heavy scalars can also decay to a pair of gluons via a loop of top or bottom quarks, and the rates are As the decay rate is a CP even quantity, in all the above decay rates, the CP even coefficient c Hif f and the CP odd onec Hif f always contribute incoherently. In addition the heavy neutral scalars, H 2 , H 3 can decay into the Z boson and the SM-like Higgs boson h, To obtain the branching ratios, we include the total width of the heavy Higgs, The branching ratios of the various physical scalars in CPV 2HDM were calculated using C2HDM HDECAY [50] through the HDECAY interface [95,96,99]. Predictions for gluon-fusion Higgs production at hadron colliders for the scalars were obtained using tabulated results from SusHi [94] at NLO. We insert the CPV 2HDM model by turning on the CP violating angles, α b and α c . From our parameter space investigations, we find that the maximum CP violating angle α b allowed by the Higgs signal strength measurements from the LHC is around α b 0.4 as shown on the bottom right panel in Fig. 1. In Fig. 6 and Fig. 7, we present the results for the the production cross sections times the branching ratios of H 2 and H 3 , respectively, for various mixing angles cos(β − α) and CPV angles. For consistency, we maintain the same parameter space analysed in Fig. 3.
In our analysis, the second and third neutral Higgs are close in mass, as a result of the ST U constraints, and because of this, the heavy neutral Higgs bosons decay only into SM particles. The mass of the charged Higgs boson is initially randomly scanned although its value is restricted to lie within the allowed range (note that the third neutral Higgs mass is not an independent parameter, but computed from the mass of the lighter neutral Higgs bosons and the various angles as shown in Eq. 2.15).
We show results for the production and decay of the heavy Higgs slightly away from the alignment limit (in which cos (β − α) = 0). The results shown in Fig. 6 and Fig. 7 are such that cos(β − α) = 0.05 on the top left, cos(β − α) = 0.0175 on the top right, cos(β − α) = 0.1 on the bottom left and cos(β − α) = 0.15 on the bottom right, and are all in a parameter space region allowed by LHC Higgs data from Fig. 1.
Note that the mass regions in each plot are different. This is because we show, in each panel, the M H2 and M H3 regions and α c range surviving the constraints for the chosen values of α b and cos(β − α), consistent with the allowed ones in Fig. 3.
In this regime, the decays into gauge bosons are open, in both the CP-even final states (ZZ and W W ) and in the CP-odd final state (Zh). We plot as solid (dotted) lines, the production times branching ratios for H 2 for the minimum and maximum allowed values for α c from Fig. 3, while keeping the same colour convention for each channel.
In all plots, H 2 → tt is the dominant decay mode, owing to small values of tan β. Production and decay rates into gauge bosons can reach O(1) pb depending on the values of α b and α c ; similar cross-section values can be reached for the CP-odd final state Zh, for different values of the CPV mixing angles. In fact, in all panels we observe how the cross sections between CP-even gauge boson final states and CP-odd final states can flip by changing the value of α c from its minimum value (solid lines) to its maximum value (dotted lines).
In the two top panels, with α b = 0.1, α b = 0.2, when the value of α c is close to 0 (solid lines), the decays into Zh represent the main decay channel involving gauge bosons, with the CP-even decay channels into W W and ZZ being between one to two orders of magnitude smaller. This shows that in this limit of small CPV mixing angles, the state H 2 is "mostly" CP-odd. However, as we increase the value of α c (dotted lines), the decays into CP-even gauge boson channels become larger than into the CP-odd channel Zh, and so the state H 2 can be considered as "mostly" CP-even.
In the lower panels the situation is similar, even though the value of α b is larger, with α b = 0.3, α c = 0.4. The allowed values of α c are constrained to be negative, but again when the absolute value of α c is small (solid lines) the decays in to Zh represent the main decay channel with a gauge boson, so that again, for smaller CPV angles, the state H 2 is "mostly" CP-odd. For a larger absolute value of α c (dotted lines), the CP-even gauge boson decays become larger than the CP-odd channel Zh, and so the state H 2 is now "mostly" CP-even.
In the top panels and the bottom left panel, intermediate values of α c do allow for relatively similar windows for the mass of H 2 . However for the case α b = 0.4 only a limited range of mass values for H 2 satisfies the constraints, and this is distinct for every value of α c . To showcase this effect in the bottom right panel, we also plot the production and decay cross section of H 2 for an intermediate value of α c , shown by dashed lines. The effect of changing α c is entirely non-linear: while the minimum to maximum are slightly shifted from the left (lower masses) to the right (higher masses), the plots for intermediate values of α c do not fall in-between but are shifted further to the right. However both the dashed and solid lines correspond to a "mostly" CP-odd scalar, whereas the "mostly" CP-even scalar corresponds to the dotted line, for the lowest allowed scalar masses.
We expect that plotting all allowed values of α c will result in uninterrupted lines, each segment corresponding to different α c 's. The slight shifting in masses is noticeable in the other panels, but not as pronounced as for α b = 0.4. Fig. 7 analyses the same production and decay cross sections for H 3 , the heaviest neutral boson. We vary the parameters over the same allowed region as in Fig. 6. The plots show the expected behavior of the heaviest neutral scalar relative to the scalar H 2 with respect to its couplings to gauge bosons. For the parameter values where the scalar H 2 was found to be "mostly" CP-odd (small values of α c ), the scalar H 3 is "mostly" CP-even, meaning its decay cross section into W W or ZZ is dominant over the decay into Zh. And conversely when H 2 is "mostly" CP-even, H 3 then behaves as "mostly" CP-odd.
In all cases the dominant decay channels are tt, W W , ZZ and Zh. The approximate degeneracy of the masses M H3 and M H2 has attracted some interest from phenomenologists. At the LHC, experiments looked for the four lepton final states, resulting from decays into ZZ, often referred to as the golden channel when searching for additional heavy scalar resonances. This is motivated by the fact that the corresponding SM backgrounds are small and controllable. Encouragingly, the four-lepton analyses from ATLAS [100] and CMS [101,102] indicate some enhancement in the event rates in final states with high invariant masses. These analyses were interpreted to be consistent with a broad resonance structure around 700 GeV [103], and most recently, also interpreted as a double peak of degenerate states around 680 GeV in the 2HDM [104]. In their case, the degeneracy is imposed to fit the data, while in our analysis, this scenario is completely consistent and imposed by the experimental and theoretical constraints.

VI. SUMMARY AND CONCLUSION
In the present study, we have explored the allowed parameter space of the Type-II Two Higgs Doublet Model with CP violation, and the implications for the production and decays of the additional neutral Higgs bosons at the LHC. After implementing the model with approximate Z 2 symmetry, we set some parameters to zero to forbid flavorchanging neutral currents, but still allowing for complex parameters in the potential. We then proceed by restricting the CP-violating parameter space using both theoretical considerations and experimental data constraints.
First, we showed that after imposing constraints from the LHC-Higgs data some regions of parameter space very close to alignment as well as relatively away from it survive (the alignment region is such that the heavy Higgs fields decouple from the rest of the spectrum leaving only a single SM-like Higgs within the scalar sector).
Within the allowed regions of parameter space, we imposed theoretical constraints involving vacuum stability, unitarity, and perturbativity, electroweak precision constraints, as well as experimental constraints from B-physics and electric dipole moments. We compute the relevant couplings and perform a random scan over two of the 2HDM physical masses, M H2 (the second lightest neutral Higgs) and M H ± , within parameter regions which satisfy the Higgs data and all other constraints.
The main findings presented in this work, representing new contributions to the study of CVP 2HDM can be listed as follows.
• We clearly presented restrictions on the parameter space of the CPV 2HDM coming from LHC Higgs data in the plane tan β -cos(β − α) (the mixing angle in the neutral Higgs boson sector), for fixed values of α b , the amount of CPV admixture in the lightest (SM-like) Higgs. Our plots, which for α b = 0 coincide with available analyses for the CP-conserving 2HDM, also show how the parameter space allowed by Higgs data shrinks with increasing α b .
• Throughout the CPV parameter space explored, tan β ∼ 0.9 emerges as a solid constraint caused by the necessity of having close to degenerate neutral heavy scalar masses in order to pass precision electroweak tests and EDM's.
• We proved analytically (and confirmed numerically) that in the CPV 2HDM the masses of the heavy neutral Higgs, M H2 and M H3 , as well as that of the charged Higgs, M H ± , are bounded from above. We demonstrated that this is a feature specific to CPV only, and is absent in the CPC 2HDMs.
• We also show that, for most of the parameter region allowed, and except for small corners of the parameter space (for the largest α b allowed, and away from alignment), the two heavier neutral Higgs bosons H 2 and H 3 are constrained to be approximately degenerate in mass. This is a condition imposed by S, T parameters. The latter also affect λ 5 , a measure of CP violation in the potential, which is quite large in these scenarios.
• While the masses of the charged and neutral Higgs bosons are varied in a large parameter range, they are constrained to lie in a small mass window. We also separate clearly restrictions coming from oblique S, T parameters and electric dipole moments, from the region which obeys all other experimental and theoretical bounds (such as Higgs data, and B-physics constraints).
• While the exact value of the mass window M H ± − M H2 depends on cos(β − α), α b and α c (the amount of CPV admixture in H 2 and H 3 ), the window is usually only 80-250 GeV wide. There is lower bound on M H ± > 650 GeV as required by B physics. Thus our analysis restricts M H ± , for a given choice of parameters, both from above and below.
• Our plots also show how, increasing the value for α c the two heavier Higgs boson flip from "mostly CP-even" to "mostly CP-odd". And while the dominant decay is into tt, the decays into bosonic states are next largest and can be used to distinguish this scenario.
In conclusion, our work shows that the parameter space of the CPV 2HDM is very constrained, and the neutral Higgs bosons in the CPV Type-II 2HDM show some interesting features, which could distinguish them from Higgs bosons in other models at the colliders.