The C2HDM revisited

The complex two-Higgs doublet model is one of the simplest ways to extend the scalar sector of the Standard Model to include a new source of CP-violation. The model has been used as a benchmark model to search for CP-violation at the LHC and as a possible explanation for the matter-antimatter asymmetry of the Universe. In this work, we re-analyse in full detail the softly broken $\mathbb{Z}_2$ symmetric complex two-Higgs doublet model (C2HDM). We provide the code C2HDM_HDECAY implementing the C2HDM in the well-known HDECAY program which calculates the decay widths including the state-of-the-art higher order QCD corrections and the relevant off-shell decays. Using C2HDM_HDECAY together with the most relevant theoretical and experimental constraints, including electric dipole moments (EDMs), we review the parameter space of the model and discuss its phenomenology. In particular, we find cases where large CP-odd couplings to fermions are still allowed and provide benchmark points for these scenarios. We examine the prospects of discovering CP-violation at the LHC and show how theoretically motivated measures of CP-violation correlate with observables.


Introduction
After the discovery of the Higgs boson by the ATLAS [1] and CMS [2] collaborations the Large Hadron Collider (LHC) has now entered the Run 2 stage with a 13 TeV centre-ofmass energy. The search for physics beyond the Standard Model (BSM) is now the main goal of the LHC experiments. An important motivation of BSM physics is to provide new sources of CP-violation to fulfil the three Sakharov criteria for baryogenesis [3]. The simplest extension of the scalar sector that can provide a new source of CP-violation, the two-Higgs doublet model (2HDM) [4], can address this issue. Reviews of the 2HDM can be found in [5][6][7]. One of the simplest ways of extending the SM with a CP-violating scalar sector without the addition of new fermions is to add a scalar doublet to the SM content and build a 2HDM softly broken Z 2 symmetric potential with two complex parameters. This model, now known as C2HDM, was first discussed in [8] and has only one independent CP-phase and a simple limit leading to its CP-conserving version. The model has been the subject of many studies [9][10][11][12][13][14][15][16][17][18]. More recently a comparison with a number of other extensions of the SM was performed [19]. The NLO QCD corrections to double Higgs decays were also recently calculated in [20].
One of the primary goals of the LHC Run 2 is to look for new particles but also to probe the CP parity of both the discovered Higgs boson and of any further scalars yet undiscovered. Due to its simplicity, the C2HDM is the ideal benchmark model to test the CP quantum numbers of the scalars at the LHC. The CP-nature of the scalars, in their Yukawa couplings, can be probed directly either in the production or the decays of the Higgs bosons. In this case, it is the relation between the scalar and the pseudoscalar components in the Yukawa couplings that is probed. Some examples of the use of asymmetries to probe the CP-nature of the Higgs boson in the top Yukawa coupling were discussed in [21][22][23] while the decays of the tau leptons were used to probe the tau Yukawa coupling [24][25][26][27][28]. Correlations in the momentum distributions of leptons produced in the decays of the Higgs boson to gauge bosons were used to probe the CP-nature of the Higgs boson couplings to gauge bosons [29][30][31] by the ATLAS [32] and CMS [33] collaborations. The most general CP-violating HV V coupling was used, and limits were set on the anomalous couplings. However, the C2HDM has SM-like couplings to the gauge bosons -it is just the SM coupling multiplied by a number. The anomalous couplings only appear at loop level and are consequently very small. Therefore, in this model, only the Yukawa couplings can lead to direct observations of CP-violation.
There are, however, other ways to probe CP-violation using only inclusive observables. As proposed in [34,35] several combinations of three simultaneously observed Higgs decay modes can constitute an undoubtable sign of CP-violation. In a CP-conserving model, a decay of the type H i → H j Z would imply opposite CP parities for H i and H j . In a renormalisable theory, a Higgs boson decaying to a pair of gauge bosons has to be CPeven 1 .
Hence, the combination of the decays H i → H j Z, H i → ZZ and H j → ZZ is a clear sign of CP-violation. In [35] seven classes of decays were defined, some of which signal CP-violation for any extension of the SM, while others are not possible in a CP-conserving 2HDM but can occur in the C2HDM or even in models with 3 CP-even states. An example of the latter is H i → ZZ (with i = 1, 2, 3) that could signal CP-violation but could also happen in a model with at least 3 CP-even states. In some extensions of the SM, even the ones with just two extra scalars, as is the case of the SM plus a complex singlet, all H i are CP-even, and therefore the decays H i → ZZ all happen at tree-level. From the above classes, we are especially interested in the ones that include the already observed decay of the SM-like Higgs boson, denoted by h 125 , to gauge bosons, h 125 → ZZ. Furthermore, decays of the type H i → ZZ and H i → H j Z were the subject of searches during Run 1 which will proceed during Run 2. Therefore, if a new Higgs boson is observed in the final 1 A CP-odd scalar decays to a pair of gauge bosons at one-loop. It was shown for the CP-conserving 2HDM that the corresponding width is several orders of magnitude smaller than the corresponding tree-level one, H → ZZ [36,37]. states with gauge bosons and also in the final state h 125 Z the scalar sector is immediately established to be CP-violating.
In the C2HDM there is only one independent CP-violating phase. Hence, the only three basis-invariant quantities that signal CP-violation [38,39] are all related and proportional to that phase. In this work, we are interested in finding a relation between the production rates in each of the CP-violating classes and variables that signal CP-violation. Therefore we will test a number of variables in the most interesting classes that include the already observed h 125 → ZZ. Finally, decays of the type H j → H i h 125 can only happen in very specific models and searches for this particular channel are important and should be a priority for Run 2.
The LHC Run 2 will bring increasing precision in the measurement of Higgs production and decay rates. The ATLAS and CMS collaborations have started testing new models such as the singlet extension of the SM and the CP-conserving 2HDM. To probe other extensions of the SM new tools are needed to increase the precision both in the Higgs production rates and also in the Higgs decay widths. Hence, one of the main purposes of this work is the release of the C2HDM_HDECAY code, an implementation of the C2HDM in HDECAY v6.51 [40,41]. The code is entirely self-contained, and the widths include the most important state-of-the-art higher order QCD corrections and the relevant off-shell decays.
The paper is organised as follows. In section 2 we describe the model and in section 3 we discuss the parameter space of the model given the most relevant theoretical and up to date experimental constraints. We identify situations where h 125 could have remarkable properties and give benchmark points for these scenarios. Namely, situations where its couplings to fermions have a large CP-odd component and the possibilities that h 125 is the heaviest of the three scalars in the C2HDM. In section 4, we discuss the production rates of processes that constitute classes of CP-violating decays and relate them with a number of variables that can probe CP-violation. In section 5 we discuss Higgs-to-Higgs decays in the framework of the C2HDM. Our conclusions are presented in section 6. In appendix A we write the Feynman rules for the C2HDM in the unitary gauge. In appendix B we describe the C2HDM_HDECAY code.

The complex two-Higgs doublet model
The version of the complex two-Higgs doublet model we discuss in this work has an explicitly CP-violating scalar potential, with a softly broken Due to the hermiticity of the Lagrangian, all couplings are real except for m 2 12 and λ 5 . We write each of the doublets Φ i (i = 1, 2) as an expansion around the real vacuum expectation values (VEVs) v 1 and v 2 , in terms of the charged complex fields (φ + i ) and the real neutral fields (ρ i and η i ). The doublets then read The minimum conditions for the potential are where we have introduced We define two CP-violating phases φ(m 2 12 ) and φ(λ 5 ) as Equation (2.5) shows us that these two phases are not independent. We can re-write eq. (2.5) as We choose both vacuum expectation values v 1 and v 2 to be real which together with the condition φ(λ 5 ) = 2 φ(m 2 12 ) [8] ensures that the two phases cannot be removed simultaneously. Otherwise we are in the CP-conserving limit of the model.
The Higgs basis [38,39] {H 1 , H 2 } is defined by the rotation The model has three neutral particles with no definite CP quantum numbers, H 1 , H 2 and H 3 , and two charged scalars H ± . The mass matrix of the neutral scalar states is diagonalised via the orthogonal matrix R [10]. That is, for which we choose the form The Higgs boson masses are ordered such that Note that the mass basis and the Higgs basis are related through    where the matrix T used in ref. [6] for the expression of the oblique radiative corrections is defined as and R H is given by Our choice of the 9 independent parameters of the C2HDM is: and Re(m 2 12 ). With this choice, the mass of the heaviest neutral scalar is a dependent parameter, given by , (2.20) and the parameter space points will have to comply with m H 3 > m H 2 . We will briefly describe here the couplings of the Higgs bosons with the remaining SM fields. A longer list is provided in appendix A, and the full set is contained in the web page [42]. The Higgs couplings to the massive gauge bosons V = W, Z are given by  where g H SM V V denotes the corresponding SM Higgs coupling, given by where g is the SU (2) gauge coupling and θ W is the Weinberg angle. The effective couplings can be written as The Yukawa sector is built by extending the Z 2 symmetry to the fermion fields such that flavour changing neutral currents (FCNC) are absent at tree-level [43,44]. There are four possible Z 2 charge assignments and therefore four different types of 2HDMs described in table 1. The Yukawa Lagrangian has the form These couplings were implemented in the code HDECAY [40,41] which provides all Higgs decay widths and branching ratios including the state-of-the-art higher order QCD corrections and off-shell decays. The description of the new C2HDM_HDECAY code is presented in appendix B.

Experimental and Theoretical Restrictions
The C2HDM was implemented as a model class in ScannerS [45,46]. The most relevant theoretical and experimental bounds are either built in the code or acessible via interfaces with other codes. We have imposed all available constraints on the model and performed a parameter scan. The resulting viable points are the basis for our phenomenological analyses for the LHC Run 2.
The theoretical bounds included in ScannerS are boundness from below and perturbative unitarity [47][48][49]. Contrary to the SM, in the 2HDM coexisting minima can occur at tree-level. Therefore we also force the minimum to be global [50], precluding the possibility of vacuum decay. The points generated comply with electroweak precision measurements, making use of the oblique parameters S, T and U [6]. We ask for a 2σ compatibility of S, T and U with the SM fit presented in [51]. The full correlation among these parameters is taken into account.
The charged sector of the C2HDM has exactly the same couplings as in the 2HDM. 2 Therefore, the exclusion bounds on the m H ± − t β plane can be imported from the 2HDM. The most constraining bounds on this plane come from the measurements of B → X s γ [52][53][54][55][56]. The latest 2σ bounds on this plane were obtained in [56] and force the charged Higgs mass to be m H ± > 580 GeV for models Type II and Flipped, almost independently of tan β. Due to the structure of the charged Higgs couplings to fermions, in models Type I and Lepton-Specific the bound has a strong dependence on tan β. In fact, for tan β ≈ 1 the bound is about 400 GeV while the LEP bound derived from e + e − → H + H − [57] (approximately 100 GeV) is recovered for tan β ≈ 1.8. We further apply the flavour constraints from R b [52,58]. All the constraints are checked as 2σ exclusion bounds on the m H ± − t β plane.
The SM-like Higgs boson is denoted by h 125 and has a mass of m h 125 = 125.09 GeV [59]. We exclude points of the parameter space with the discovered Higgs signal built by two nearly degenerate Higgs boson states by forcing the non-SM scalar masses to be outside the mass window m h 125 ± 5 GeV. Compatibility with the exclusion bounds from Higgs searches is checked with the HiggsBounds code [60][61][62], while the individual signal strengths for the SM-like Higgs boson are forced to be within 2σ of the fits presented in [63]. Branching ratios and decay widths of all Higgs bosons are calculated with the C2HDM_HDECAY code, described in appendix B, which is an implementation of the C2HDM model into HDECAY v6.51. The code has state-of-the-art QCD corrections and off-shell decays, but off-shell decays of one scalar into two are not included. The Higgs boson production cross sections via gluon fusion (ggF ) and b-quark fusion (bbF ) are calculated with SusHiv1.6.0 [64,65], which is interfaced with ScannerS, at next-to-next-to-leading order (NNLO) QCD. As the neutral scalars have no definite CP, we need to combine the CP-odd and the CP-even contributions by summing them incoherently. That is, µ F is given by where in the denominator we neglected the bbF cross section since in the SM it is much smaller than gluon fusion production. The CP-even Higgs production cross sections in association with a vector boson (V H) and in vector boson fusion (V BF ) give rise to the normalised production strengths because QCD corrections cancel at NLO. The effective couplings are defined in eq. (2.21). Both CP-even and CP-odd components contribute to the cross sections for associated production with fermions. As these have different QCD corrections [66], we opted for using leading order production cross sections, and write the strengths as with the coupling coefficients defined in eq. (2.24). We then use HiggsBounds via the ScannerS interface to ensure compatibility at 2σ with all available collider data. Regarding the SM-like Higgs boson, the parameter space is forced to be in agreement with the fit results from [63]. That is, the values given in [63], with µ xx defined as for H i ≡ h 125 , are within 2σ of the fitted experimental values. Here H SM denotes the SM Higgs boson with a mass of 125.09 GeV. As the C2HDM preserves custodial symmetry, We use this method for simplicity. Note that performing a fit to current Higgs data is likely to give a stronger bound than this approach. The C2HDM is a model with explicit CP-violation in the scalar sector. Therefore, there are a number of experiments that allow us to constrain the amount of CP-violation in the model. The most restrictive bounds on the CP-phase [14] (see also [67][68][69][70][71][72]) originate from the ACME [73] results on the ThO molecule electric dipole moment (EDM). All points in parameter space have to conform to the ACME experimental results. As we will discuss in detail later, the bounds can only be evaded either in the CP-conserving limit of the model or in scenarios where cancellations between diagrams with different neutral scalar particles occur [69,70]. The cancellations are related to the orthogonality of the R matrix in the case of almost degenerate scalars [16]. Our results are required to be compatible with the measured EDM values in [73] at 90% C.L. We finalise with a word of caution regarding the electron EDM. The authors of ref. [69] argue that the EDM could be up to one order of magnitude larger than the bound presented by ACME, due to the large uncertainties in its extraction from the experimental data.
In our scan, one of the Higgs bosons H i is identified with h 125 . One of the other neutral Higgs bosons is varied between 30 GeV ≤ m H i < 1 TeV while the third neutral Higgs boson is not an independent parameter and is calculated by ScannerS, but its mass is forced to be in the same interval. In Type II and Flipped, the charged Higgs boson mass is forced to be in the range while in Type I and Lepton-Specific we choose Taking into account all the constraints, in order to optimise the scan, we have chosen the following regions for the remaining input parameters:

Constraints on the Parameter Space: c e versus c o
In this subsection, we confront the C2HDM parameter space with all the restrictions presented above. Our aim is to see what is the structure of the remaining parameter space and, in particular, to study the CP-nature of the 125 GeV scalar, encoded in the couplings . These test the CP content of h 125 . We know that h 125 must have some CP-even content because it couples at tree level to ZZ. However, in a theory with CP-violation in the scalar sector (such as the C2HDM), h 125 could have a mixed CP nature. This possibility can be probed in the couplings to fermions in a variety of ways. The simplest case occurs if c e f c o f = 0, meaning that in the coupling to some fermion there are both CP-even and CP-odd components, thus establishing CPviolation. A more interesting case occurs if h 125 (or some other mass eigenstate) has a pure scalar component to a given type of fermion (f ) and a pure pseudoscalar component to a different type of fermion (g). This would render c e f c o We will now show that in the case of the C2HDM this possibility is no longer available for Type I and it is only available in Type II if h 125 = H 2 . In contrast, the possibility still exists in the Lepton-Specific and Flipped models. We will also give several benchmark points to allow a more detailed study of these scenarios.
Applying all the above constraints on the parameter space, we have obtained the points in parameter space that are still allowed. We call this sample 1. We have also performed a scan with the EDM constraints turned off, which we call sample 2. The left plot of figure 1 displays the mixing angles α 1 , which (in the displayed case of H 1 = h 125 ) mixes the CP-even parts of the two Higgs doublets and α 2 , which parametrises the amount of CP-violating pseudoscalar admixture, for Type I. In the right plot we present the CP-odd and the CPeven components of the Yukawa couplings for the Type I model with and without the EDM constraints. Both plots clearly show that the EDM constraints have little effect on the mixing angle |α 2 |, which can go up to 25 • when all constraints are taken into account. The maximum value of this angle can be understood from the bound 0.79 < µ V V < 1.48. In fact, as previously shown in [17] the fact alone that µ V V > 0.79 forces the angle |α 2 | to be below ≈ 27 o . Coming from the bound on µ V V , this constraint will be approximately the same for all types (before imposing EDM constraints), as will become clear in the next plots.
We are also interested in the wrong-sign regime, defined by a relative sign of the Yukawa coupling compared to the Higgs-gauge coupling, realized for c e b < 0. As shown previously in [74,75], the right plot again demonstrates that the wrong-sign regime is in conflict with the Type I constraints because the Yukawa couplings cannot be varied independently.
In figure 2 we present the distributions of the angle α 1 and α 2 for samples 1 and 2 and for a Type II model. The EDM constraints, applied in our sample 1, strongly reduce |α 2 | to small values. Only for scenarios around the maximal doublet mixing case with α 1 ≈ π/4, α 2 can reach values of up to ∼ ±20 • .
The , as discussed in [17]. The inclusion of the EDM constraints, however, clearly rules out this possibility when h 125 = H 1 . Nevertheless, the wrong-sign  regime (c e b < 0) is still possible in the C2HDM for down-type Yukawa couplings. The electron EDM has no discernable effect on the allowed coupling to up-type quarks, as can be read off from the right plot.
The situation changes when we take Type II with h 125 = H 2 , as shown in figure 4. One can still find scenarios where the top coupling is mostly CP-even (c e t 1), while the bottom coupling is mostly CP-odd (c o b 1). It is noteworthy that the electron EDM kills all such points in Type II when h 125 = H 1 , but that they are still allowed in Type II when In table 3 we present three benchmark scenarios in Type II with large CP-violation in the Yukawa sector. The first scenario, BP2m, has maximal c o b with nearly vanishing c e b . Since c e t is always 1 this means that c e t c o b is maximal here. The other two scenarios BP2c and BP2w both have maximal c e b c o b but are in the correct sign and wrong-sign regime, respectively. As discussed above all of the scenarios with large CP-violation in the Yukawa couplings of h 125 require H 2 = h 125 in Type II models.   The situation is even more interesting in the other two Yukawa types. Figure 5 displays the Yukawa couplings for the Lepton-Specific model with and without the EDM constraints.  The down-type quark couplings are tied to the up-type couplings and, thus, heavily constrained to lie close to the SM (fully CP-even) solution. However, figure 5 (left) shows that the charged lepton couplings can still be mostly (or even fully) CP-odd, despite the current EDM constraints. Similarly, in the Flipped model the bottom quark can couple to h 125 in a fully CP-odd fashion, as shown in the left plot of figure 6. In it, we display the Yukawa couplings for the Flipped model with and without the EDM constraints. As will be discussed below, such large CP-odd components are still viable in both Lepton-Specific and Flipped models due to cancellations between the various diagrams entering the EDM calculation. This is also true for Type II when H 2 = h 125 . But it is important to stress that they are not due to large α 2 values, as illustrated in figure 7. The values for α 2 are small, but c o (h 125 bb) grows very fast as α 2 departs from α 2 = 0 for large tan β. It grows roughly as c o (h 125 bb) ∼ s 2 tan β.   Table 4: Benchmark points with large pseudoscalar Yukawa couplings in the Lepton-Specific (LS) and Flipped types. Lines 1-8 contain the input parameters; lines 9-11 the derived third Higgs mass and the relevant Yukawa couplings (multiplied by sgn(c(h 125 V V ))) and the last five lines the signal strengths of h 125 .
In table 4 we present further benchmark scenarios in the Lepton-Specific and Flipped types. The BPLSm scenario has a maximal c o τ coupling with tiny c e τ , thus here h 125 appears CP-even in its couplings to quarks and CP-odd in its couplings to leptons. In BPLSc and BPLSw the product c e τ c o τ is maximal. In the Flipped model we provide three benchmark points. The first one, BPFm, again having maximal c o b while BPFc and BPFw both have large c o b c e b but with opposite signs and are therefore close to the correct and wrong-sign limit, respectively. Note that, in contrast to the Type II, all benchmark points except BPLSw have H 1 = h 125 . We chose BPLSw with H 2 = h 125 since maximising c e τ c o τ close to the wrong-sign limit of the Lepton-Specific model leads to this mass ordering.

The case H 3 = h 125
One interesting possibility in the C2HDM is to have the 125 GeV scalar discovered at LHC (h 125 ) coincide with the heaviest Higgs boson (H 3 ). This possibility is excluded for Type II and Flipped, as the B-physics constraints impose that the charged Higgs boson must be quite heavy, m H + > 580 GeV. This poses a problem with the electroweak precision tests, specially with the T parameter. Indeed, the way to accommodate the experimental bounds on T is to have a spectrum that has some degree of degeneracy. Requiring m H + > 580 GeV implies that the other Higgs boson masses cannot all be below 125 GeV. Thus, m H + > 580 GeV is not compatible with m H 2 < 125 GeV.
However, H 3 = h 125 is feasible for Type I and Lepton-Specific, as in these cases Bphysics constraints only impose m H + > 100 GeV (as explained, for low tan β this bound could be slightly higher). The situation here is quite fascinating, because it highlights an interesting complementarity between LHC and the old LEP results. The relevant Feynman rules are:

Measures of CP-violation
Throughout this section we will use the notation H ↓ (H ↑ ) to designate the lightest (heaviest) of the non-h 125 neutral Higgs bosons. Their mass can be below or above 125 GeV. The  manifestation of CP-violation in models with two Higgs doublets can be probed with a number of variables even if there is only one independent CP-violating phase in the scalar sector. The most obvious variable is the phase in one of the complex parameters of the potential, that is, either λ 5 or m 2 12 . As the two phases are not independent, we have opted to use φ(λ 5 ). As discussed in [34,35], there are several combinations of Higgs decays that are a clear signal of CP-violation in any extension of the SM. Others, like the simultaneous observation of the three decays H i → ZZ, i = 1, 2, 3, enable us to distinguish the 2HDM from the C2HDM but are not unequivocal signs of CP-violation in general extensions of the SM. The question we are trying to address now is: are there variables that quantify CPviolation from a theoretical point of view and also correlate with the rates of a combination of decays which would establish CP-violation experimentally?
Starting with φ(λ 5 ), we present in figure 9 three classes of CP-violating processes as a function of the CP-violation phase |φ(λ 5 )| for Type I (left column) and Type II (right column). In the first row we show pp → H ↓ → Zh 125 against pp → H ↓ → ZZ, in the second row we have pp → H ↑ → Zh 125 against pp → H ↑ → ZZ and in the third row we plot pp → H ↓ → ZZ against pp → H ↑ → ZZ. These classes of decays were chosen because together with the already observed process h 125 → ZZ they can be used to identify CP violation, and also because searches for the other processes were performed for Run 1 and will continue during Run 2. There are two striking features in the plots. First, the production rates in Type I are almost one order of magnitude above the ones for Type II. This is because there are constraints that act more strongly on Type II like b → sγ or the EDM constraints as we will see later. Ultimately, it is due to the different structure of the Yukawa couplings: in Type I all Yukawa couplings are equal, making the model harder to constrain. Second, there is no correlation between the magnitude of φ(λ 5 ) and the  production rates, because the points with larger φ(λ 5 ) are almost evenly spread throughout the plot. It is not that we expected that large values of the CP-violating phase would correspond to large production rates for any of the processes in the plot. In fact, maximal CP-violation is attained for specific sets of values of the angles but the production rates are a complicated combination of all parameters of the model. Still some colour structure could have emerged in the plots, but this was not the case.
In view of the negative result for the complex phase, we have looked for other variables [77,78] that could be used as a measure of CPV in the production and decay of the three neutral Higgs bosons [34,35]. The sets of variables proposed in the literature are basically of two types: the ones where the CP-violating variables appear in a sum of squares, and the ones where they appear in a product. The important difference between them is that while the former is zero only when there is no CP-violation in the model, the latter can be zero even if the model is CP-violating. However, if CP is conserved, both variables are zero.
We start by defining a multiplicative variable first proposed in [77] that allows us to distinguish a CP-conserving from a CP-violating 2HDM, which, when the couplings g V V H i are normalised to the SM one, satisfies With the purpose of finding quantities invariant under a basis transformation that change sign under a CP transformation, it was shown in [38] that the simplest CP-odd invariant that can be built from the mass matrix is Furthermore, any other CP-odd invariant built from the mass matrix alone has to be proportional to J 1 [38]. There is a relation between ξ V and J 2 1 that can be written as It should be noted that even if J 1 = 0, CP-violation could occur in the scalar sector (see [38] for details). This measure of CP-violation is not applicable to the fermion-Higgs sector, where the invariants of ref. [39] apply instead. Variables that are a clear signal of CP-violation can be built with the scalar and pseudoscalar components of the Yukawa couplings. In fact, if a model has CP-violating scalars at tree level its Yukawa couplings have the general form c e f + ic o f γ 5 . Thus, as discussed at the beginning of subsection 3.2, variables of the type c e c o clearly signal CP-violation in the model. Therefore, we define the normalised multiplicative variables [78] as measures of CP-violation in the up-and down-quark sectors, respectively. The corre-sponding normalised sum variables [78] are defined as again for the up-and down-quark sectors, respectively. We will use the same combinations of decays as in figure 9 to see if any of the variables proposed can provide a relation between the amount of CP-violation and the production rates in case these processes are observed at the LHC. In figure 10 we present the three classes that signal CP-violation as a function of the CP-violation variables for the Type I C2HDM. Since h 125 → ZZ (not shown) was already observed, any pair of two processes appearing in the plots combined with the latter form a signal of CP-violation. In the first row we show pp → H ↓ → Zh 125 against pp → H ↓ → ZZ, in the second row we have pp → H ↑ → Zh 125 against pp → H ↑ → ZZ and in the third row we plot pp → H ↓ → ZZ against pp → H ↑ → ZZ. In each column, we show the variable that is being probed. For the case of Type I, the top and bottom Yukawa couplings are the same. The general picture is that there is no striking correlation between the large values for the variables (more yellow points) and the large production cross sections for each process. There is a quite even spread of yellow points for the sum variable, ζ f . This is because even if the product of scalar/pseudoscalar components of the SM-like Higgs boson is very constrained, any of the products of the other two Higgs bosons can be very large, yielding a large value for the sum in ζ f . Therefore, variables of this type can always be large, and we will not show them in the remaining plots. Regarding the variables γ f , we can see some structure in the plots as there are cases where more yellow points are clustered closer to the maximum values of the production rates. This trend is clearer in the last row where kinematics play a less important role. In fact, the first row is the most constrained regarding the phase space available while the last row is almost symmetric regarding the reduction of the phase space. The difference between the first and second row regarding kinematics is just that the first row deals with the lightest non-125 Higgs boson while the second row deals with the heaviest one.
In figure 11 we present the same three classes of CP-violating processes as a function of CP-violating variables for the Type II C2HDM. What we see for this Yukawa type is that the distribution of the yellow points is more structured and more clustered in specific regions. In fact, the yellow points tend to cluster more in the regions where the production rates are larger. This behaviour is more striking in the first column, for the variable ξ V , where all yellow points are in the parameter region where the production rates are maximal. For the remaining variables, the distribution of yellow points is again not so structured. Still, all variables are larger when the production rates are both large and are much smaller when both production rates are smaller. However, in all plots, there are always points with small values of the CP-violating variables and large values of the production rates. Hence,  although the variables show us a trend, they are not conclusive as a measure of CP-violation in the scalar sector.
In figure 12 we present the same classes in the Type II C2HDM as in figure 11, but with signal strengths (see eq. (3.4)) within 5% of the SM values. This gives us a hint on what to expect at the end of the LHC Run2, or at the high luminosity LHC. There is a clear effect in reducing the production rates but not in the distribution of yellow points. The main difference is that now no yellow points appear in the first column which means that the points with very large rates were excluded. The distribution of points in the other columns did not change significantly, but the points with the higher rates were also excluded as for the first row. In conclusion, there seems to be an overall reduction in the parameter space of the model leading to smaller production rates. In figure 13 we show the CP-violating parameter ξ V as a function of tan β for Type I (top left), Type II (top right), Lepton-Specific (bottom left) and Flipped (bottom right). The lighter points have passed all the constraints except for the EDM bounds, while the darker points have passed all constraints. In Type I there is not much difference between the two sets of points, and there are no special regions regarding the allowed values of tan β. Also, the maximum value for ξ V is around 0.2 almost independently of tan β. For Type II, the results are much more striking. After imposing the EDM constraints, we end with two almost straight lines (one for tan β ≈ 1 and the other for ξ V ≈ 0), as well as a region around tan β ≈ 3 permitting values of ξ V up to 0.6. This means that tan β can only be large when we approach the CP-conserving limit except for a few points, which lie in the wrong-sign regime. Hence, in a Type II model, points with significant CP-violation can occur for tan β ≈ 1 in the alignment limit or for large tan β for the wrong sign limit. The situation in Flipped is similar to Type II, with a lower maximum value of ξ V ∼ 0.2 after imposing the EDM constraints.
In figure 14 we show the individual contributions to the EDM coming from W -loops, fermion-loops, charged Higgs loops and charged Higgs plus W -loops. For each C2HDM type, we have grouped the contributions to the EDM according to their relative sign. For example, in Type II the contributions of the W -loops (y-axis) and the sum of the contributions of the fermion loops, charged Higgs loops and charged Higgs plus W -loops (x-axis) have opposite sign. The grey shaded region represents the parameter space excluded by the EDM constraints only. The colour code represents the absolute value of the magnitude of φ(λ 5 ). The first important difference between the models is that the maximum value of the individual contributions is around two orders of magnitude smaller in Type I than in Type II. This implies that Type I is less constrained by the EDM bounds. Therefore the remaining constraints play a much more important role in Type I than in Type II. This also leads to the distribution of values for the CP-violating phase in the figure. In Type II large values of φ(λ 5 ) prefer regions where either the EDM contributions are very small, i.e. there are cancellations between loops in the individual contributions, or rely on huge cancellations between different contributions. This is in contrast to Type I, where large values of the CP-violating phase can be found all over the allowed region. The Flipped case behaves roughly like Type I, while the Lepton-Specific case behaves like Type II. This is very different from the usual behaviour of the Yukawa types. The reason is that most observables are dominated by quark effects giving similar behaviour to Type I/Lepton-Specific and Type II/Flipped. Since we are, however, looking at the EDM of the electron the lepton Yukawa couplings are the most important ones and those are equal in Type I/Flipped and Type II/Lepton-Specific.
Finally, it is important to comment on the different impact of the EDM constraints on the different models, regarding the possibility of having large pseudoscalar Yukawa couplings. For simplicity, we focus our discussion on the H 1 = h 125 case. In Type I the pseudoscalar Yukawa couplings c o f are proportional to sin(α 2 )/ tan β for any fermion type f . Since α 2 is at most of the order 20 • and tan β is constrained to be above 1, c o f will be below about 0.4. For the other Yukawa types either c o b or c o τ are proportional to sin(α 2 ) tan β. In this case, even if α 2 is small the pseudoscalar coupling can be substantial for large tan β. We have observed that in Type II the values of α 2 can go up to 20 • while in the Flipped model they barely reach 5 • . In both cases, c o b could still be large because large values of tan β are allowed in both models. So what is the reason for having large c o b for the Flipped model but not for Type II when H 1 = h 125 ? The main reason is that the EDM constraints are less stringent in the Flipped model than in Type II due to relative signs between the CP-odd lepton Yukawa couplings and the remaining CP-odd Yukawa couplings. This ends up flipping several signs in one model relative to the other leading to cancellations between loops and much smaller individual EDM contributions. In Type II we found that this leads to the result that we can have large tan β only when α 2 is very close to zero which is not the case for the Flipped model.

Higgs-to-Higgs decays
In the last section we have discussed the relation between some of the classes of decays that probe CP-violation with a number of variables proposed in the literature to measure CP-violation in the scalar sector. These particular classes were chosen not only because they can yield large production rates but also because there are currently searches being performed for these channels at the LHC that have already started during Run 1. There are, however, other classes of decays that can probe CP-violation and were the subject of a study that led to the production of benchmarks for Run 2 [35,79] Lepton-Specific Flipped  σ(pp → H ↑ → ZZ) < 1 fb for the lower plots. It is clear from the plots that, with the extra restriction on the ZZ final state, the cross sections now barely reach 10 fb for the two decay scenarios and for all types. Hence, although possible, it will be very hard to detect the new scalars in the h 125 h 125 final state if they are not detected in the ZZ final state. One should note that the cross section for di-Higgs production in the SM is about 33 fb. Consequently, a resonant di-Higgs final state such as the one presented in figure 15 would easily be detected because the cross sections can reach the pb level. However, it is also clear that once we force σ(pp → H ↑ → ZZ) < 1 fb it is no longer possible to detect these di-Higgs states even at the High Luminosity LHC.
In figure 17 we show the production rates for the process pp → H ↑ → H ↓ h 125 as a function of the heavier Higgs mass, for all C2HDM types. For this channel the rates can reach at most about 100 fb, and only for Type I and Flipped. In Type II the rates are at most at the fb level. The rates for the H ↓ h 125 final state with the extra condition σ(pp → H ↑ → ZZ) < 1 fb are shown in figure 18. The maximum rates (for low masses) are now reduced by about a factor of 5 for Type I. However, the rates do not decrease much for the Flipped C2HDM, and some signal at LHC Run 2 could point to this C2HDM type. Finally, although H ↑ → H ↓ h 125 appears hard to detect in these models it is nevertheless a clear signal of non-minimal models and should therefore be a priority for the LHC Run 2.
We end this section with the production rate for the process pp → h 125 → H ↓ H ↓ as a function of the lighter Higgs mass for the various C2HDM types, which are shown in figure 19. Most points correspond to a mass of the heavier state above 125 GeV. But, as shown in figure 8 Lepton-Specific Flipped  (for Type II and Flipped). In order to understand what would be the behaviour when choosing a definite final state, we have checked that the rates pp → h 125 → H ↓ H ↓ → bbτ + τ − are still above the pb level for all model types.

Conclusions
In this work we have analysed in detail a minimal complex version of the 2HDM, known as the C2HDM. The inclusion of all available experimental and theoretical constraints allowed us to present an up to date status of the model. We have shown that large CP-odd Yukawa couplings of h 125 are still possible in all Yukawa types except Type I. However, in Type II this is only a possibility if H 2 = h 125 . We provided two different interesting kinds of benchmark points which are in agreement with all current observations: for h 125 coupling like a scalar to some fermions and like a pseudoscalar to others, and for scenarios where the scalar and pseudoscalar Yukawa couplings of h 125 to some fermion are of similar size.
The model has only one CP-violating parameter. In previous works we have presented classes of three decays that are a clear sign of CP-violation, when all three are observed. In this work we looked for correlations between the production rates of the processes in each class and the phase φ(λ 5 ) that measures CP-violation. We have shown that there is no correlation for the most relevant classes. We have then tested the correlation with other CPviolating variables proposed in the literature. The conclusion is that some correlation can be seen between the production rates and the variables. This is particularly true for the variable ξ V and even more for Type II. However, in most cases there is almost no correlation between the high rates of the CP-violating processes and the proposed CP-violating variables. The results also tell us that measuring small rates should not be interpreted as a sign of a small CP-violating angle.
Finally, we presented in the last section the production rates for the scenarios where the Higgs bosons decay to two other scalars. The search for a scalar decaying into two h 125 Higgs bosons was already performed during Run 1. However, the search for H i → H j h 125 has not started yet. It is clear from the results that it will be much harder to probe CP-violation with classes of decays that involve scalar to scalar decays. However, on one hand it could be that all other decays would be very constrained. On the other hand, the observation of scalar to scalar decays would be a first step to reconstruct the Higgs potential. Hence, these searches should be a priority for Run 2 and especially for the high luminosity LHC.
As part of a common effort to a proper interpretation of the LHC results we release the code C2HDM_HDECAY that calculates the decay widths and branching ratios of all the C2HDM scalars including state-of-the-art QCD corrections.

A.1 Couplings of neutral Higgs bosons to fermions
The couplings of neutral Higgs bosons to fermions can be written in general for all the neutral Higgs bosons of the model H i in a compact form with the coefficients presented in table 2.

A.2 Couplings of charged Higgs bosons to fermions
The couplings of the charged Higgs bosons to fermions can be expressed in the following Lagrangian , for quarks and leptons in an obvious notation and P L = (1 − γ 5 )/2 and P R = (1 + γ 5 )/2. The couplings η L,R are given in table 5. In these expressions, we neglect the masses of the neutrinos, so in the last line in table 5 the zeros mean that the corresponding mass in eq. (A.2) is zero.

A.3 Cubic interactions of neutral Higgs bosons
The cubic interactions of neutral Higgs bosons, [H i , H j , H k ], have long expressions. We do not write them here but we collect them in one web page [42]. The expressions there are the Feynman rules without the i.

A.4 Cubic interactions of neutral and charged Higgs bosons
These are (on the right-hand side of these expressions we write the Feynman rule including the i), +R i2 (−(Re(λ 5 ) − λ 2 + λ 4 ) cos(β) 2 sin(β) + λ 3 sin(β) 3 ) where the λ i or g H i H + H − can be read from eq. (A.3). The λ i are in the notation used in ref. [16] and should not be confused with the parameters in the potential.

A.5 Cubic interactions with gauge bosons
One gauge boson Two gauge bosons where in eq. (A.8) and eq. (A.9) we used a notation to make contact with our previous conventions [16].

A.6 Quartic interactions with Higgs bosons
Again these interactions involve quite long expressions and they are given in our web page [42].

B The Fortran Code C2HDM_HDECAY
The code C2HDM_HDECAY is the implementation of the CP-violating 2HDM in the program HDECAY v6.51 [40,41], which is written in Fortran77. All changes with respect to the C2HDM have been included in the main file hdecay.f, which is now called chdecay.f. Further linked routines have been taken over from the original HDECAY program, so that the code is completely self-contained. The decay widths are computed including the most important state-of-the-art higher order QCD corrections and the relevant off-shell decays. Note, that it does not include off-shell Higgs-to-Higgs decays, but only on-shell decays into a lighter Higgs pair. The QCD corrections can be taken over from the SM and the minimal supersymmetric extension (MSSM), respectively, for which HDECAY was originally designed. The electroweak corrections on the other hand cannot be adapted from the available corrections in the SM and/or MSSM so that they have been consistently turned off.
The C2HDM input parameters are specified in the input file c2hdecay.in which has been obtained by extending the original HDECAY input file hdecay.in. The C2HDM branching ratios and total widths are calculated after setting the input value C2HDM= 1 in c2hdecay.in. The required input parameters are set in the block 'complex 2 Higgs Doublet Model'. Here the user specifies the values of two of the neutral Higgs boson masses, the third one is computed from the input values, and the charged Higgs boson mass. Furthermore, the mixing angles α 1,2,3 and tan β have to be set as well as the real part of the mass parameter M 2 12 . The type of the fermion sector is chosen via the input variable TYPE_cp. The values 1, 2, 3 and 4 correspond to the I, II, Lepton-Specific and Flipped types, respectively. For illustration, we display part of an example input file for the C2HDM case.
The code is compiled with the file makefile. By typing make an executable file called run is produced. The program is executed with the command run. It calculates the branching ratios and total widths which are written out together with the mass of the decaying Higgs boson. The names of the output files are br.Xy_C2HDM. Here X=H1, H2, H3, c denotes the decaying Higgs particle, where 'c' refers to the charged Higgs boson. Files with the suffix y=a contain the branching ratios into fermions, with y=b the ones into gauge bosons and the ones with y=c, d the branching ratios into lighter Higgs pairs or a Higgs-gauge boson final state. For illustration, we present the example of an output that has been obtained from the above input file. The produced output in the four output files br.H3y_C2HDM for the heaviest neutral Higgs boson is given by The program files can be downloaded at the url: http://www.itp.kit.edu/∼maggie/C2HDM There one can find a short explanation of the program and information on updates and modifications of the program. Furthermore, sample output files are given for a sample input file.