Models with Extended Higgs Sectors at Future $e^+ e^-$ Colliders

We discuss the phenomenology of several Beyond the Standard Model (SM) extensions that include extended Higgs sectors. The models discussed are: the SM extended by a complex singlet field (CxSM), the 2-Higgs-Doublet Model with a CP-conserving (2HDM) and a CP-violating (C2HDM) scalar sector, the singlet extension of the 2-Higgs-Doublet Model (N2HDM), and the Next-to-Minimal Supersymmetric SM extension (NMSSM). All the above models have at least three neutral scalars, with one being the 125 GeV Higgs boson. This common feature allows us to compare the production and decay rates of the other two scalars and therefore to compare their behaviour at future electron-positron colliders. Using predictions on the expected precision of the 125 GeV Higgs boson couplings at these colliders we are able to obtain the allowed admixtures of either a singlet or a pseudoscalar to the observed 125 GeV scalar. Therefore, even if no new scalar is found, the expected precision at future electron-positron colliders, such as CLIC, will certainly contribute to a clearer picture of the nature of the discovered Higgs boson.


Introduction
The discovery of the Higgs boson by the LHC experiments ATLAS [1] and CMS [2] has triggered the search for new scalars as predicted by Beyond the Standard Model (BSM) models with extended Higgs sectors. Although no new scalars were found at the LHC up until now, and no solid hints of new physics have been reported by the LHC collaborations, the increasing precision in the measurement of the Higgs couplings to fermions and gauge bosons has dramatically reduced the parameter space of BSM models. Hence, it could be that at the end of the LHC run we will not discover any new particle and will have to rely on future colliders to further search for new physics.
In this work we discuss the phenomenology of several BSM extensions that include extended Higgs sectors at a future electron-positron collider. The models discussed are: the SM extended by a complex singlet field (CxSM), the 2-Higgs-Doublet Model with a CP-conserving (2HDM) and a CP-violating violating (C2HDM) scalar sector, the singlet extension of the 2-Higgs-Doublet Model (N2HDM), and the Next-to-Minimal Supersymmetric SM extension (NMSSM). All the above models have at least three neutral bosons, with one being the 125 GeV Higgs boson. This common feature allows us to compare the production and decay rates of the other two scalars.
The models are investigated by performing parameter scans that take into account the most relevant theoretical and experimental constraints. Our main goal is to answer two questions. The first one is what can an electron-positron collider tell us about the nature of the discovered Higgs boson -is it just part of a doublet, or two doublets; has it a singlet component or a CP-violating one, and if so how large? The second one is, to what extent can a future electronpositron collider distinguish between the different BSM versions if a new Higgs boson is found? Are we able to disentangle the models based on Higgs rate measurements? We hope that we can shed some light on the relevance of a future electron-positron collider for BSM Higgs searches.
The outline of the paper is as follows. In section 2 we briefly introduce the models under study. In section 3 we describe the constraints on the models and how the scans over the parameter space are performed. In section 4 we discuss what we can learn about the nature of the discovered 125 GeV scalar after CLIC. In section 5 the signal rates of the non-SM-like Higgs bosons are compared within the different models. Our conclusions are given in section 6.

Description of the Models
We start with a very brief description of the models analysed in this work and we refer the reader to [3] for a detailed description. Here we will just set our notation and define the free parameters used in each model.

The Complex Singlet Extension of the SM
The first model we discuss is an extension of the SM by a complex scalar field (CxSM) which is defined by a scalar potential with a softly broken global U (1) symmetry given by where S = S + iA is a hypercharge zero scalar field and the soft breaking terms are written in parenthesis. We write the fields as , (2.2) where v ≈ 246 GeV is the vacuum expectation value (VEV) of the h field and v S and v A are the VEVs of the real and imaginary parts of the complex singlet field, respectively. Except for the soft breaking terms, all parameters are real as required by the hermiticity of the potential.
As we further impose invariance under S → S * (or A → −A), a 1 and b 1 are real. We choose to work in the broken phase (all three VEVs are non-zero) because this phase leads to mixing between the three CP-even scalars. Their mass eigenstates are denoted by H i and are obtained from the gauge eigenstates via the rotation matrix R parametrised as where we have defined s i ≡ sin α i and c i ≡ cos α i , and without loss of generality the angles vary in the range  5) and the remaining parameters are determined internally in ScannerS [4,5] fulfilling the minimum conditions of the vacuum. In the broken phase, the couplings of each Higgs boson, H i , to SM particles are rescaled by a common factor R i1 . The expression for all couplings can be found in the appendix B.1 of [6]. All Higgs branching ratios, including the state-of-the art higher order QCD corrections and possible off-shell decays can be obtained from sHDECAY [6] 1 which is an implementation of the CxSM and also the RxSM both in their symmetric and broken phases in HDECAY [7,8]. A detailed description of the program can be found in appendix A of [6].

The 2HDM and the C2HDM
In this section we introduce the real (2HDM) and complex (C2HDM) versions of a particular 2-Higgs-Doublet model, where we add a second doublet to the SM scalar sector. The Higgs potential is invariant under the Z 2 transformations Φ 1 → Φ 1 and Φ 2 → −Φ 2 and is written as By extending the Z 2 symmetry to the fermions we guarantee the absence of tree-level Flavour Changing Neutral Currents (FCNC). If all parameters of the potential are real and the VEVs in each doublet are also real the potential is CP-conserving and we call the model 2HDM; if the VEVs are real but m 2 12 and λ 5 are complex, with different phases, the model is CP-violating and we call it C2HDM [9]. Both the 2HDM and the C2HDM have two charged Higgs bosons and three neutral scalars. In the 2HDM the neutral scalars are h and H, the lighter and the heavier CP-even states, while A is the CP-odd state. In the C2HDM we have three Higgs mass eigenstates H i (i = 1, 2, 3) with no definite CP and that are ordered by ascending mass according to m H 1 ≤ m H 2 ≤ m H 3 . The rotation matrix, R, that diagonalises the mass matrix is parametrised as defined for the complex singlet extension case in Eq. (2.3) and with the same range as in Eq. (2.4) for the mixing angles. The CP-conserving 2HDM is obtained from the C2HDM by setting α 2 = α 3 = 0 and α 1 = α + π/2 [10]. In this case the CP-even mass eigenstates h and H are obtained from the gauge eigenstates through the rotation parametrised in terms of the angle α. The 2HDM has eight independent parameters while the C2HDM has nine independent parameters. We define for both versions of the model v = v 2 1 + v 2 2 ≈ 246 GeV and tan β = v 2 /v 1 . For the 2HDM we choose as independent parameters v , tan β , α , m h , m H , m A , m H ± and m 2 12 , (2.7) while for the C2HDM we choose [11] v , tan β , where m H i and m H j denote any two of the three neutral Higgs bosons but where one of them is the 125 GeV scalar. The remaining mass is obtained from the other parameters [11]. We write the couplings to massive gauge bosons (V = W, Z) of the Higgs boson H i in the C2HDM as where [12] c and g H SM V V denotes the SM Higgs coupling factors. In terms of the gauge boson masses M W and M Z , the SU (2) L gauge coupling g and the Weinberg angle θ W they are given by Both the 2HDM and C2HDM are free from tree-level FCNCs by extending the global Z 2 symmetry to the Yukawa sector. The four independent Z 2 charge assignments of the fermion fields determine the four types of 2HDMs depicted in Table 1. The Yukawa Lagrangian is defined by where ψ f is the fermion field with mass m f . In Table 2 we present the CP-even and the CPodd components of the Yukawa couplings, c e (H i f f ) and c o (H i f f ), respectively [12]. All Higgs branching ratios can be obtained from C2HDM HDECAY [13] 2 which implements the C2HDM in u-type d-type leptons  HDECAY [7,8]. These include state-of-the art higher order QCD corrections and possible off-shell decays. The complete set of Feynman rules for the C2HDM is available at: http://porthos.tecnico.ulisboa.pt/arXiv/C2HDM/ where for the SM subset the notation for the covariant derivatives is the one in [14] with all η's positive, where the η's define the sign of the covariant derivative (see [14]). Note that the 2HDM branching ratios are part of the HDECAY release (see [7,8,15] for details).

The N2HDM
The version of the N2HDM used in this work was discussed in great detail in [16]. This extension consists of the addition of an extra doublet and an extra real singlet to the SM field content. The potential is invariant under two discrete Z 2 symmetries. The first Z 2 symmetry is just a generalisation of the one used for the 2HDM in order to avoid tree-level FCNCs, and that is softly broken by m 2 12 ; the second one is defined as and it is not explicitly broken. Φ 1 and Φ 2 are doublet fields and Φ S is a singlet field. The most general form of this scalar potential invariant under the above transformations is 3 (2.14) The doublet and singlet fields after electroweak symmetry breaking can be parametrised as where v 1,2 are the VEVs of the doublets Φ 1 and Φ 2 , respectively, and v S is the singlet VEV. The singlet VEV breaks the second Z 2 symmetry, precluding the existence of a dark matter candidate. As this is a CP-conserving model, with no dark matter candidate, we end up with three CP-even scalars, one of which plays the role of the 125 GeV Higgs boson, a CP-odd scalar and two charged scalars. The orthogonal matrix R that diagonalises the mass matrix is again parametrised as in Eq. (2.3) in terms of the mixing angles α i with the same ranges as before, see Eq. (2.4). The physical CP-even eigenstates, denoted by H 1 , H 2 and H 3 , are ordered by ascending mass as (2.16) We choose as the 12 independent parameters the set The expressions of the quartic couplings in terms of the physical parameter set can be found in appendix A.1 of [16]. All Higgs branching ratios, including the state-of-the art higher order QCD corrections and possible off-shell decays can be obtained from N2HDECAY 4 [16,18] which implements the N2HDM in HDECAY [7,8].

The NMSSM
Supersymmetric models require the introduction of at least two Higgs doublets. The NMSSM extends the two Higgs doublet superfieldsĤ u andĤ d of the Minimal Supersymmetric extension (MSSM) by a complex superfieldŜ. The µ problem of the MSSM is thus solved dynamically when the singlet field acquires a non-vanishing VEV. The NMSSM Higgs sector consists of seven physical Higgs states after EWSB. These are, in the CP-conserving case investigated in this work, three neutral CP-even, two neutral CP-odd ones and a pair of charged Higgs bosons. The NMSSM Higgs potential is derived from the superpotential, the soft SUSY breaking Lagrangian and the D-term contributions. The scale-invariant NMSSM superpotential reads in terms of the hatted superfields For simplicity, we have only included the third generation fermion superfields here. They are given by the left-handed doublet quark ( Q 3 ) and lepton ( L 3 ) superfields and the right-handed singlet quark ( t c R , b c R ) and lepton ( τ c R ) superfields. The first term in Eq. (2.18) takes the role of the µ-term µĤ dĤu of the MSSM superpotential, the term cubic in the singlet superfield breaks the Peccei-Quinn symmetry thus avoiding a massless axion and the last three terms represent the Yukawa interactions. The soft SUSY breaking Lagrangian consists of the mass terms for the Higgs and the sfermion fields, that are built from the complex scalar components of the superfields, The contribution to the soft SUSY breaking part from the trilinear soft SUSY breaking interactions between the sfermions and the Higgs fields reads where the A's denote the soft SUSY breaking trilinear couplings. The gaugino mass parameters M 1,2,3 of the bino (B), winos (W ) and gluinos (G), respectively, that contribute to the soft SUSY breaking are summarised in We will allow for non-universal soft terms at the GUT scale. The expansion of the tree-level scalar potential around the non-vanishing VEVs of the Higgs doublet and singlet fields, leads to the Higgs mass matrices for the three scalars (h d , h u , h s ), the three pseudoscalars (a d , a u , a s ) and the charged Higgs states (h ± u , h ∓ d ) that are obtained from the second derivative of the scalar potential. The VEVs v u , v d and v s are chosen to be real and positive. Rotation with the orthogonal matrix R S that diagonalises the 3 × 3 mass matrix squared, M 2 S , of the CP-even fields, yields the CP-even mass eigenstates H i (i = 1, 2, 3), (2.23) They are ordered by ascending mass, The CP-odd mass eigenstates A 1 and A 2 result from a rotation R G separating the massless Goldstone boson followed by a rotation R P into the mass eigenstates, which are ordered by ascending mass, The three minimisation conditions of the scalar potential are used to replace the soft SUSY breaking masses squared for H u , H d and S in L mass by the remaining parameters of the tree-level scalar potential. This leads to the following six parameters parametrising the tree-level NMSSM Higgs sector, We have chosen the sign conventions such that λ and tan β are positive, whereas κ, A λ , A κ and µ eff are allowed to have both signs. Contrary to the non-SUSY Higgs sector extensions introduced in the previous sections, the Higgs boson masses are not input parameters. They are instead calculated from these, including higher order corrections. These are crucial to shift the mass of the SM-like Higgs boson to the observed value of 125 GeV. Due to these corrections also the soft SUSY breaking mass terms for the scalars and the gauginos as well as the trilinear soft SUSY breaking couplings contribute to the Higgs sector.

Parameter Scans
The analyses are performed with points, each corresponding to a set of the parameters chosen for a given model, that are in agreement with the theoretical and experimental constraints. The discovered SM-like Higgs boson mass is taken to be [19] m h 125 = 125.09 GeV , (3.26) and we suppress interfering Higgs signals by forcing any other neutral scalar to be outside the m h 125 ± 5 GeV mass window. Any of the Higgs bosons is allowed to be the discovered one except for charged and pure pseudoscalar particles. The vacuum expectation value v is fixed by the W boson mass and all calculations of cross sections and branching ratios do not include electroweak corrections as they are not fully available for all models. All models except for the NMSSM, the scan of which will be described below, have been implemented as ScannerS model classes. This allowed us to perform a full parameter scan that simultaneously applies the constraints we will now briefly describe. The theoretical bounds are common to all models although with different expressions. We force all potentials to be bounded from below, we require that perturbative unitarity holds and that the electroweak vacuum is the global minimum (using the discriminant from [20] for the C2HDM). Compatibility with electroweak precision data for the CxSM was imposed by a 95% C.L. exclusion limit from the electroweak precision observables S, T and U [21,22] -see [23] for more details. The same constraints for the C2HDM use the expressions in [24] while for the N2HDM we use the formulae in [25,26]. For the computed values of S, T and U we ask for a 2σ compatibility with the SM fit [27] taking into account the full correlation among the three parameters.
95% C.L. exclusion limits on non-observed scalars have been applied by using HiggsBounds [28] which include LEP, Tevatron and up-to-date LHC experimental data. Compatibility with the Higgs data is enforced using the individual signal strengths fit [29] for the h 125 . The branching ratios for the different models were calculated using the modified versions of HDECAY as described in the previous sections. All scalar production cross sections can be easily obtained from the corresponding SM one except for the gluon fusion (ggF ) and b-quark fusion (bbF ) which were determined using SusHi v1.6.0 [30,31]. For the C2HDM, the CP-even and the CP-odd Yukawa coupling contributions are calculated separately and then added incoherently, giving where we neglected the bbF contribution for the SM in the denominator. Analogous expressions were used for the other models which do not have a CP-odd component. Models with two doublets with or without extra neutral singlets always have a pair of charged Higgs bosons. In this study the charged Higgs Yukawa couplings are always proportional to two parameters only: the charged Higgs mass and tan β. These couplings are constrained by the measurements of R b [32,33] and B → X s γ [33][34][35][36][37], which yields 2σ exclusion bounds on the m H ± − t β plane. The latest calculation of [37] enforces, almost independently of the value of tan β, m H ± > 580 GeV (3.28) in the Type II and Flipped models while in Type I and Lepton Specific models this bound is not only much weaker but it has a much stronger dependence on tan β.
Finally there are bounds that apply only to the C2HDM because constraints on CP violation in the Higgs sector arise from electric dipole moment (EDM) measurements. Among these the EDM of the electron imposes the strongest constraints [38], with the experimental limit given by the ACME collaboration [39]. We require our results to be compatible with the values given in [39] at 90% C.L. A detailed discussion of the constraints specific to the C2HDM can be found in [13]. With all the above constraints taken into account, the initial range of parameters chosen for each model is as follows:

• The CxSM Parameter Range Scan
The non-125 GeV Higgs bosons are chosen to be in the range (3.29) The VEVs v A and v S are varied in the range The mixing angles α 1,2,3 vary within the limits (3.31) • The (C)2HDM Parameter Range Scan The angles vary in the range and − π 2 ≤ α 1,2,3 < π 2 . (3.33) The value of Re(m 2 12 ) is in the range 0 GeV 2 ≤ Re(m 2 12 ) < 500000 GeV 2 . (3.34) In type II, the charged Higgs mass is chosen in the range • The N2HDM Parameter Range Scan In view of what was discussed for the previous models, the ranges for the parameters of the N2HDM are 580 GeV ≤ m H ± < 1 TeV (type II) .

(3.39)
Note that the 125 GeV Higgs boson can be the lighter as well as the heavier scalar. This possibility is not excluded in any of the models.

The NMSSM Parameter Scan
For the NMSSM parameter scan we proceed as described in [6,40] and shortly summarise the main features. We use the NMSSMTools package [41][42][43][44][45][46] to calculate the spectrum of the Higgs and SUSY particles with higher order corrections included. The package also checks for the constraints from low-energy observables. It provides the input required by HiggsBounds which verifies compatibility with the exclusion bounds from Higgs searches. The relic density obtained through an interface with micrOMEGAS [46] is required not to exceed the value measured by the PLANCK collaboration [47]. The spin-independent nucleon-dark matter direct detection cross section, that is also obtained from micrOMEGAS, is required not to violate the upper bound from the LUX experiment [48]. We furthermore test for compatibility with the direct detection limits from XENON1T [49]. The mass of one of the neutral CP-even Higgs bosons has to lie between 124 and 126 GeV. The signal strengths of this Higgs boson have to be in agreement with the signal strength fit of [29] at the 2 × 1σ level. For the production cross sections, gluon fusion and bb annihilation, we take the SM cross sections and multiply them with the effective couplings obtained from NMSSMTools. The SM cross section values are obtained from SusHi [30,31]. In gluon fusion the next-to-leading order (NLO) corrections are included with the full top quark mass dependence [50] and the next-to-next-to-leading order (NNLO) corrections in the heavy  quark effective theory [51][52][53][54][55]. For Higgs masses below 300 GeV the next-to-next-to-next-toleading order (N 3 LO) corrections are taken into account in a threshold expansion [56][57][58][59]. For masses above 50 GeV bb annihilation cross sections that match between the five-and four-flavor scheme are used obtained in the soft-collinear effective theory [60,61]. They equal the results from [62,63]. For masses below 50 GeV, cross sections obtained in the Santander matching [64] are used, with the five-flavor scheme cross sections from [65] and the four-flavor scheme ones from [66][67][68]. The branching ratios are obtained from NMSSMTools. We cross-checked the Higgs branching ratios of NMSSMTools against NMSSMCALC [69]. We demand the masses of all Higgs bosons to be separated by at least 1 GeV in order to avoid overlapping signals. The obtained parameter points are also checked for compatibility with the SUSY searches at LHC. We require the gluino mass and the lightest squark mass of the second generation to be above 1.85 TeV, respectively, [70]. The stops have to be heavier than 800 GeV [71] and the slepton masses heavier than 400 GeV [72]. The absolute value of the chargino mass must not be lighter than 300 GeV [73]. The scan ranges applied for the various parameters are summarised in Table 3. Perturbativity is ensured by applying the rough constraint λ 2 + κ 2 < 0.7 2 . Over the last years, predictions for the measurement of the Higgs couplings to fermions and gauge bosons were performed for CLIC for some benchmark energies and luminosities. Table 4 shows the expected precision in the measurement of the Higgs couplings and was taken from [74] (see [74,75] Table 4 have a central value of 1. Small deviations from the central value will not have a significant effect on our results because the errors are very small. If significant deviations from the SM predicted values are found the data has to be reinterpreted for each model. Starting with the simplest extension, the CxSM, there are either one or two singlet components that mix with the real neutral part of the Higgs doublet. In the broken phase, where there are no dark matter candidates, the admixture is given by the sum of the squared mixing matrix elements corresponding to the real and complex singlet parts, i.e. with the matrix R defined in Eq. (2.3). If a dark matter candidate is present one of the R ij , j = 2, 3, is zero. In any case the Higgs couplings to SM particles are all rescaled by a common factor. Therefore, we just need to consider the most accurate Higgs coupling measurement to get the best constraints on the Higgs admixture. The maximum allowed singlet admixture is given by the lower bound on the global signal strength µ which at present is Therefore we expect to constrain the admixture of a singlet to close to 0.1% at the end of the CLIC operation. This implies, for this particular kind of extensions, that the chances of finding a new scalar are reduced due to the orthogonality of the R matrix. Note that in the limit of exact zero singlet component the singlet fields do not interact with the SM particles. The results for a real singlet are similar, with the bound being exactly the same but with a two by two orthogonal matrix replacing R. In this case it is exactly the value 0.11% that multiplies all production cross sections of the non-SM Higgs boson, after CLIC@3TeV. We now discuss the C2HDM as this is the model with a CP-violating scalar and one that shows a quite different behaviour in the four independent Yukawa versions of the model. In fact, the constraints act very differently in the four Yukawa versions of the model as shown in [13]. This is particularly so for the EDMs [13] -while for Type II the electron EDM constraint almost kills the pseudoscalar component of the the bbH coupling, the same is not true for the Flipped model and for the pseudoscalar component of the Higgs couplings to leptons in the Lepton Specific model. Since different Yukawa couplings enter the two-loop Barr-Zee diagrams, a small EDM can either be the result of small CP-violating Yukawa couplings or come from cancellations between diagrams. This can even allow for maximally CP-violating Yukawa couplings of the h 125 in some cases [13]. So now the question is: in the long run, can CLIC give us relevant information that complements the one from EDMs? How far can one expect to go in the knowledge of the Higgs nature by putting together CLIC and EDM results, how well can one constrain the CP-violating component of the 125 GeV Higgs boson?
In Fig. 1 (left) we present the mixing angles α 2 versus α 1 for the C2HDM Type I. The blue points are for Sc1 but without the constraints from κ Hgg and κ Hγγ ; the green points are for Sc1 including κ Hgg (the measurement of κ Hγγ was not included because it is not available) and the red points are for Sc3 including κ Hgg and κ Hγγ . Note that the κ Hgg and κ Hγγ are the only measurements of couplings that can probe the interference between Yukawa couplings (in the case of κ Hgg ) and between Yukawa and Higgs gauge couplings (in the case of κ Hγγ ). In the right panel of Fig. 1 we show the pseudoscalar component of the b-quark Yukawa coupling c o b versus its scalar component c e b . Because in Type I all Yukawa couplings are equal, this plot is valid for all Type I Yukawa couplings. One can then expect, by the end of the CLIC operation, all pseudoscalar (scalar) Type I Yukawa couplings to be less than roughly 5% (0.5 %) away from the SM expectation. We again stress that this result assumes that experiments will not see deviations from the SM.
Recently, in [76] a study was performed for a 250 GeV electron-positron collider for Higgsstrahlung events in which the Z boson decays into electrons, muons, or hadrons, and the Higgs boson decays into τ leptons, which subsequently decay into pions. The authors found that for an integrated luminosity of 2 ab −1 , the mixing angle between the CP-odd and CP-even components, defined as could be measured to a precision of 4.3 o which means that this is the best bound if the central measured value of the angle is zero. Their result is translated into our notation via Taking into account the values in Fig. 1 (right)   In Fig. 2 (left) we present the mixing angles α 2 vs. α 1 for the C2HDM Type II. In the right panel we again show the pseudoscalar component of the b-quark Yukawa coupling c o b vs. its scalar component c e b . The blue points are for Sc1 without the constraints from κ Hgg and κ Hγγ . These loop induced couplings are the only ones where interference between Yukawa couplings and Higgs gauge couplings occur. Therefore, whatever the precision on the measurement of treelevel couplings is, the result will always be a ring in that plane, that will become increasingly thiner with growing precision. However, even for CLIC@350GeV, if the constraint for κ Hgg is included, the ring is reduced to the green arch shown in the figure. By the end of the CLIC operation the arch will be further reduced to the red one. As discussed in previous works, a very precise measurement of κ Hgg or κ Hγγ will kill the wrong-sign limit 5 , which corresponds in the figure to c e b = −1. Now, how do these bounds compare do the direct ones from h 125 → τ + τ − ? In Type I the same bounds apply to all ψ CP . At the same time the bound on ψ top CP is the same in all models and it was already discussed for Type I. In Type II ψ bottom CP = ψ τ CP and from Fig. 2  (right) we obtain bounds on ψ bottom CP that are of the order of 30 o for CLIC@350GeV and 15 o for CLIC@3TeV. Therefore, we conclude that for Type II the indirect bounds cannot compete with the direct ones. The EDM constraints also play a very important role in probing the CP-odd components of the couplings. In fact, in the particular scenario of the Type II C2HDM in which the lightest Higgs boson is the 125 GeV scalar, the bound is already constraining ψ bottom CP to be below 20 o [13] clearly competing with the expectations for CLIC.
The present best measurement for the electron EDM was obtained by the ACME collaboration, with an upper bound of |d e | < 9.3 × 10 −29 e cm (90% confidence) [39] and by the JILA collaboration with an upper bound of |d e | < 1.3 × 10 −28 e cm (90% confidence) [79]. ACME II is expected to increase the statistical sensitivity by an order of magnitude [80] relative to the ACME I result. There are several other planned experiments that could result in an increase in sensitivity by two to three orders of magnitude [81,82]. These experiments together with the input from CLIC would certainly improve our knowledge on the nature of the Higgs boson. Figure 3: tan β as a function of sin(α1 − π 2 ) for Type I in Sc1 (left) and Sc3 (right). The factor − π 2 is due to a different definition of the rotation angles relative to the 2HDM. Also shown in the colour code is the amount of singlet admixture present in h125.
The predictions for the N2HDM are very similar to the ones for the 2HDM and we will discuss them together. Although the N2HDM has an extra singlet field relative to the 2HDM, the couplings to gauge bosons and fermions are very similar. For instance, for the lightest Higgs boson the couplings to massive gauge bosons are related via g N 2HDM which results in some extra freedom for the N2HDM parameter space. In Fig. 3 we show tan β as a function of sin(α 1 − π 2 ) for Type I in Sc1 (left) and Sc3 (right). These are typical plots not only for a Type I N2HDM but also for a Type I 2HDM (and very similar plots are obtained for the Lepton Specific versions of both models). The only really notable difference between the N2HDM and the 2HDM is the colour bar where we show the percentage of the singlet component in the 125 GeV Higgs boson, Σ 125 = (R i3 ) 2 . In a previous work [3] we have shown that before the LHC run 2 the percentage of the singlet admixture was below 25% for Type I and the predictions for CLIC@350GeV and CLIC@3TeV are below about 0.8% and 0.2%, respectively. As expected, the allowed parameter space gets closer and closer to the SM line, that is the line sin(β − α) = 1 (alignment limit). Note that unless one detects a new particle there is no way to find the value of tan β if the models are in the alignment limit. In fact, considering that the lightest Higgs boson is the 125 GeV one, if we are in the alignment limit, sin(β − α) = 1 in the 2HDM, 6 all couplings of the 125 GeV Higgs boson to the other SM particles are independent of the value of tan β (including the triple Higgs coupling). If the 125 GeV Higgs boson is not the lightest scalar in the model, the limits change but the physics is the same. 2 ) for Type II in Sc1 (left) and Sc3 (right). The factor − π 2 is due to a different definition of the rotation angles relative to the 2HDM. Also shown in the colour code is the amount of singlet present in h125.
In Fig. 4 we show tan β as a function of sin(α 1 − π 2 ) for Type II in Sc1 (left) and Sc3 (right). These are typical plots not only for a Type II N2HDM but also for a Type II 2HDM (and very similar plots are obtained for the Flipped versions of both models). As previously discussed we see that the right leg, corresponding to the wrong-sign limit, is very dim in the left plot and vanishes in the right plot. Again, this is true for both the 2HDM and the N2HDM. As for the percentage of the singlet component, it was constrained to 55% for Type II N2HDM at the end of run 1 [3] and the predictions for CLIC@350GeV and CLIC@3TeV are below about 0.8% and 0.2%, respectively.
We end this section with a discussion on the correlations between different cross section measurements for the different models. In Fig. 5 we present µ t = σ BSM tth /σ SM tth as a function of   could hint on a specific model. Take for instance the plot on the right and let us assume that the µ's could be measured with 5% precision. In this case a measurement (µ t , µ V ) = (1, 0.85) indicates that the model cannot be the C2HDM Type II nor the NMSSM. A measurement (µ t , µ V ) = (1.2, 1.0) excludes the NMSSM but not the remaining two models, in their Type II versions. Finally, we present in Fig. 6 µ t = σ BSM tth /σ SM tth as a function of g BSM tth /g SM tth for the 2HDM and N2HDM Type I and the CxSM (left) and for the 2HDM and N2HDM Type II and the NMSSM (right) for 1.4 TeV. The behaviour is very similar for all models, the cross sections grow linearly with the top Yukawa couplings. So in this case a deviation from the SM expectation could exclude some models and at the same time make a prediction for the change in the coupling relative to the SM for a given model. Note that because e + e − →tth (for which both Yukawa couplings and Higgs gauge couplings contribute) is not kinematically allowed for 350 GeV, the study of the correlations between this process and associated or W -fusion cross sections (for which only Higgs gauge couplings contribute) can only be performed for 1.4 TeV.

Signal Rates of the non-SM-like Higgs Bosons
In this section we present and compare the rates of the neutral non-SM-like Higgs bosons in the most relevant channels at a linear collider. We denote by H ↓ the lighter and by H ↑ the heavier of the two neutral non-h 125 Higgs bosons. All signal rates are obtained by multiplying the production cross section with the corresponding branching ratio obtained from sHDECAY, C2HDM HDECAY, N2HDECAY and NMSSMCALC. For the particular processes presented in this section, there is no distinction between particles with definite CP-numbers and CP-violating ones and they are therefore treated on equal footing. The main production processes for a Higgs boson at CLIC are associated production with a Z boson, e + e − → ZH i , and W -boson fusion e + e − → ννH i . We will be presenting results for two centre-of-mass energies, √ s = 350 GeV and √ s = 1.4 TeV. In the case of the former the cross sections are comparable in the mass range presented while for the latter the W -boson fusion cross section dominates in the entire Higgs boson mass range. In order to give some meaning to the event rates presented in this section, we will use as a rough reference that at CLIC 10 −1 fb for Sc1 correspond to 50 signal events and 10 −2 fb for Sc2 correspond to 150 signal events. In Fig. 7 we present the total rate for e + e − → ννH i → ννγγ as a function of the Higgs boson mass for the CxSM and for the Type I versions of the N2HDM and C2HDM. Also shown is the line for a SM-like Higgs boson. On the left panel we present the results for the lighter Higgs boson, H ↓ , and on the right we show the results for the heavier Higgs boson, H ↑ . The trend shown in the two plots is the same for all other final states. There is a hierarchy with the points of the N2HDM reaching the largest cross sections followed closely by the C2HDM and finally by the CxSM. This is easy to understand since the CxSM is the model with the least freedom -all couplings of the Higgs boson to SM particles are modified by the same factorwhile the N2HDM is the least constrained model. This means that it is possible to distinguish between the singlet and the Type I doublet versions if a new scalar is found with a large enough rate. The γγ final state is one where the branching ratio decreases very fast with the mass.

The 350 GeV CLIC
Still it is clear that there are regions of the parameter space that have large enough production rates to be detected at the 350 GeV CLIC. We would like to stress that the behaviour seen in the plots regarding the event rates for the lighter (left) and for the heavier (right) scalar is the same for the remaining final states and we will only show plots for the lighter Higgs boson in the remainder of this section.
In Fig. 8 we present the total rate for e + e − → ZH ↓ → Zbb (left) and for e + e − → ννH ↓ → ννbb (right) as a function of m H ↓ for √ s = 350 GeV, for the NMSSM and for the Type II versions of the N2HDM and C2HDM. Clearly there is plenty of parameter space to be explored in the NMSSM and even more in the Type II N2HDM. For the Type II C2HDM, as discussed in a previous work [13], the constraints are such that points with masses below about 500 GeV are excluded. Again there are regions where the models can be distinguished but not if the cross sections are too small. As expected, for this centre-of-mass energy there is not much difference between the two production processes (for instance for a 125 GeV scalar σ(e + e − → ZH i ) = σ(e + e − → ννH i ) for √ s ≈ 400 GeV; as the scalar mass grows so does the energy for which the values of the cross sections cross). We have also checked that the behaviour of the total rates does not change significantly when the Higgs boson decays to other SM particles. That is, although the rates are much higher in H i → bb than in H i → γγ, the overall behaviour is the same. The highest rates are obtained in all models for the final states bb, W + W − , ZZ and τ + τ − .

The 1.4 TeV CLIC
As the centre-of-mass energy rises the W -fusion process becomes the dominant one. In Fig. 9 we present the total rate for e + e − → ννH ↓ → ννZZ as a function of the lighter Higgs mass for √ s = 1.4 TeV. In the left panel we show the rates for the CxSM and for the Type I N2HDM and C2HDM while in the right panel plots for the NMSSM and the Type II N2HDM and C2HDM are shown. We can expect that total rates above roughly 10 −2 fb can definitely be explored at CLIC@1.4TeV. Hence, all models can be explored in a very large portion of the parameter space and again there are regions where the models are clearly distinguishable. The plots do not present any major differences when we change the final states as previously discussed.  However, once the 350 GeV run is complete, even if no new scalar is found, the measurement of the 125 GeV Higgs couplings will be increasingly precise which in turn reduces the parameter space of the model. In Fig. 10 we present the total rate for e + e − → ννH ↓ → ννZZ as a function of the lighter Higgs boson mass for √ s = 1.4 TeV (same as Fig. 9) but where we have included the predictions on the Higgs coupling measurements after the end of the 350 GeV run. We see that after imposing the constraints on the Higgs couplings the cross sections decrease by more than one order of magnitude. We find that the models can all be probed but are no longer distinguishable just by looking at the total rates to SM particles. Interestingly, all points from the NMMSM disappear when we impose the constraints from the 350 GeV run. This is of course related to the fact that we have used the SM central values for all predictions but it could very well be that at the end of this run we could be celebrating the discovery of a new NMSSM particle -or from any another model!

Conclusions
In this work we investigated several extensions of the SM: the CxSM, the 2HDM, C2HDM and N2HDM in the Type I and Type II versions and the NMSSM. This analysis was performed for three CLIC benchmarks with centre-of-mass energies of 350 GeV, 1400 GeV and 3000 GeV. Starting by analysing the effect of the precision of Higgs coupling measurements in the different models we have probed the nature of the discovered Higgs boson by looking at its admixture with both singlet and pseudoscalar components. We concluded that a substantial improvement in constraining the admixtures of both the singlet and the pseudoscalar one from tens of percent to values below 1% is attained when going from the LHC to the last stage of CLIC. Pseudoscalar components of CP-violating scalars will also be very constrained by the precise measurements of the loop induced processes h 125 → gg and h 125 → γγ. These components can also be measured directly in physical production and decay processes such as e + e − →tth 125 or h 125 → τ + τ − . The interplay with future measurements of EDMs should also be taken into account. If no new physics is found at CLIC and the measured values for the Higgs couplings are in agreement with the SM prediction, the Higgs boson will probably prove to have origin in a doublet with any admixture below the 1% level.
We have then compared the rates in the different models in order to understand if they can be probed at the first stage of CLIC and then to understand if knowing the results of the first stage if there is still parameter space for the second stage. We also wanted to understand if there are regions of the parameter space where only one or some models are allowed. Indeed this is the case for the models that have more freedom, such as the N2HDM, having the largest rates. We have also found that there are no significant differences between the behaviour of heavy and light scalars, and while different final states for the Higgs bosons result in different event rates they do not change the rates reachable in the models relative to another.