The N2HDM under Theoretical and Experimental Scrutiny

The N2HDM is based on the CP-conserving 2HDM extended by a real scalar singlet field. Its enlarged parameter space and its fewer symmetry conditions as compared to supersymmetric models allow for an interesting phenomenology compatible with current experimental constraints, while adding to the 2HDM sector the possibility of Higgs-to-Higgs decays with three different Higgs bosons. In this paper the N2HDM is subjected to detailed scrutiny. Regarding the theoretical constraints we implement tests of tree-level perturbativity and vacuum stability. Moreover, we present, for the first time, a thorough analysis of the global minimum of the N2HDM. The model and the theoretical constraints have been implemented in ScannerS, and we provide N2HDECAY, a code based on HDECAY, for the computation of the N2HDM branching ratios and total widths including the state-of-the-art higher order QCD corrections and off-shell decays. We then perform an extensive parameter scan in the N2HDM parameter space, with all theoretical and experimental constraints applied, and analyse its allowed regions. We find that large singlet admixtures are still compatible with the Higgs data and investigate which observables will allow to restrict the singlet nature most effectively in the next runs of the LHC. Similarly to the 2HDM, the N2HDM exhibits a wrong-sign parameter regime, which will be constrained by future Higgs precision measurements.


Introduction
The discovery of the Higgs boson by the LHC experiments ATLAS [1] and CMS [2] not only marked a milestone for elementary particle physics but also opened the possiblity to search for new physics (NP) in the Higgs sector itself. Since, so far, a direct discovery of NP in the form of new particles is missing, the Higgs sector plays an increasingly important role. The manifestations of NP in the Higgs sector can be manifold [3]. An immediate direct signal of NP acting in the Higgs sector would be the discovery of additional Higgs bosons, which can be lighter or heavier than the currently known one that has a mass of 125.09 GeV [4]. Indirect signs may appear through modifications in the Higgs couplings to the Standard Model (SM) particles and hence through the observables of the 125 GeV Higgs boson when compared to the SM values. The modifications can be due to strong dynamics behind electroweak symmetry breaking (EWSB) like in composite Higgs models [5][6][7][8][9][10][11][12][13][14][15][16]. In the case of weakly coupled models with extended Higgs sectors, the SM-like Higgs boson mixes with the other Higgs bosons thus changing the couplings to the SM particles. Additionally, new non-SM particles, like e.g. the superpartners in supersymmetric extensions, can contribute to the loop-induced couplings to gluons or photons. Furthermore, the different particle content and the modified couplings induce higher order corrections to the Higgs couplings that can be substantially different from the SM case. Finally, the additional Higgs bosons open the possibility of Higgs-to-Higgs decays. These, and possibly invisible decays due to additional lighter Higgs or other particles that are stable, modify the total width and hence the branching ratios of the SM-like Higgs boson.
With the observed Higgs boson behaving very SM-like [17][18][19][20] it is clear that any extension of the Higgs sector beyond the SM (BSM) has to provide at least one CP-even Higgs boson with a mass of 125 GeV that reproduces the LHC rates. Additional Higgs bosons predicted by the model have to be compatible with the LHC exclusion bounds. A strong constraint on NP models is given by the ρ parameter. This singles out models with singlet or doublet extended Higgs sectors when some simplicity is required. 1 Doublet extended models are particularly interesting due to their relation to supersymmetry. In particular, the 2-Higgs-doublet model (2HDM) [21][22][23] has been extensively studied and considered as a possible benchmark model in experimental analyses. The 2HDM features 5 physical Higgs bosons that, in the CP-conserving version of the model, are given by 2 CP-even, 1 CP-odd and 2 charged Higgs bosons. Upon extending the model by a real scalar singlet field with a Z 2 parity symmetry, there is a symmetric phase containing a viable Dark Matter (DM) candidate. This version of the Next-to-Minimal 2HDM (N2HDM) has been subject to numerous investigations, see e.g. [24][25][26][27][28][29][30][31][32][33][34][35][36][37], while in [38] the phenomenology of the N2HDM with non-vanishing vacuum expectation value (VEV) for the singlet field (Z 2 broken phase) has been discussed. In the Z 2 broken phase, after EWSB the N2HDM Higgs sector consists of 3 neutral CP-even scalars, 1 CP-odd and 2 charged Higgs bosons. The Higgs mass eigenstates, which are now superpositions of the singlet and doublet fields, have an interesting phenomenology that is not only governed by the mixing properties of the doublet fields but also by the amount of singlet admixture to the Higgs mass eigenstates. Thus, the couplings to SM particles can be diluted to such an extent that light Higgs bosons are not excluded by Higgs boson searches at LEP, Tevatron and the LHC in the low-mass range. Such light Higgs bosons then allow for Higgs decays of the heavier Higgs bosons into a pair of light Higgs states. Higgs-to-Higgs decays provide alternative production channels for the heavier Higgs bosons and give access to the trilinear Higgs self-couplings. Their measurement is crucial for our understanding of the Higgs mechanism [39][40][41]. Furthermore, the larger number of parameters, as compared e.g. to the 2HDM, allows for more flexibility in the Higgs sector while being simultaneously in accordance with the experimental and theoretical constraints. This is also the case for the Next-to-Minimal Supersymmetric Model (NMSSM) [42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57] whose Higgs sector is based on two doublets and one complex singlet field. Through the SUSY sector the NMSSM encounters even more parameters. The NMSSM Higgs potential, however, is subject to supersymmetric relations. In particular, the Higgs potential parameters of the two Higgs doublets in the NMSSM are given in terms of the gauge boson couplings, so that neither the NMSSM Higgs masses nor the trilinear Higgs self-couplings can become arbitrarily large. The larger Higgs spectrum of the NMSSM, with an additional pseudoscalar Higgs, and the different Higgs self-couplings induce differences in Higgs-to-Higgs decays as compared to the N2HDM. The supersymmetric relations furthermore lead to constraints in the Higgs boson couplings to the SM particles. Therefore differences in the Higgs rates and also in the coupling patterns, namely in the coupling sum rules, are to be expected. The N2HDM, on the other hand, does not have to respect supersymmetry relations among the masses and couplings. This leads to much more freedom in the choice of parameters of the model and to very different patterns in the couplings of the SM-like Higgs boson. Another class of models that can also provide such extra freedom is given by the scalar singlet framework where one adds hypyercharge-zero singlet scalar fields to the SM, i.e. without introducing any extra doublets [58][59][60][61][62][63][64][65][66][67][68][69][70][71][72][73][74]. Though in some of these models [74] one can still obtain a rich phenomenology of Higgs-to-Higgs decays with several Higgs bosons, the coupling structure of the new Higgs bosons to other SM particles is typically controlled by only a global suppression factor (relative to an SM-like Higgs boson). Thus, this provides still less structure than the N2HDM. Finally, being a model with a 2HDM-like sector, the N2HDM also contains a richer spectrum with a charged and a CP-odd Higgs boson that induce different signatures when compared with scalar singlet models. The N2HDM therefore provides an important benchmark model with a very distinct Higgs boson phenomenology as compared to other commonly studied beyond the SM extensions. The potential of various observables to distinguish between all these models requires a detailed comparison, which is beyond the scope of this paper and deferred to future work. The LHC Higgs data constrain possible deviations induced by NP to be close to the SM case so that only precision measurements allow to reveal BSM signals in the Higgs sector. This calls not only for advanced experimental techniques but also for very precise predictions from the theory side. Thus parameters and observables have to be computed including higher order corrections. Moreover the allowed parameter space of the model has to be evaluated very carefully by checking for consistency with the relevant theoretical and experimental constraints. Only for these parameter regions predictions become meaningful and can be used as guidelines for the experiments. Constraints from the experimental side arise from the Higgs data. The N2HDM has to provide at least one SM-like Higgs boson with a mass of 125 GeV. The additional Higgs bosons must comply with the LHC exclusion limits. Furthermore B-physics and lowenergy physics constraints have to be respected as well as the compatibility with the electroweak precision data. Finally, the symmetric N2HDM, which features a Dark Matter (DM) candidate has to comply with the measured value of the relic density 2 . Theoretical constraints that have to be fulfilled are: that the Higgs potential is bounded from below, that the chosen vacuum is a global minimum and that perturbative unitarity holds. In the N2HDM, the conditions for the first two requirements can be derived from the literature. There exists, however, no analysis so far of all the minima of the N2HDM.
In this work we will determine the allowed parameter space of the N2HDM in the broken phase without applying any approximations on the singlet admixture to the SM-like Higgs boson. Besides taking into account the experimental constraints, we will, in particular, investigate for the first time in great detail the conditions on the N2HDM potential that guarantee tree-level perturbative unitarity, that the vacuum is stable and that the minimum is the global one. We will present the full analysis of the global minimum of the N2HDM potential, which was performed for the first time in [75] where more details can be found. We have implemented the N2HDM in HDECAY [76,77]. This code, called N2HDECAY, computes the N2HDM Higgs decay widths and branching ratios including the state-of-the-art higher order QCD corrections and off-shell decays. Furthermore, the model has been included in ScannerS [70,78] along with the theoretical conditions and the available experimental constraints. Then, this allowed us to perform extensive scans in the parameter space of this model taking into account the experimental and theoretical constraints. We will subsequently investigate the features of the surviving parameter space and the implications for LHC phenomenology.
The outline of the paper is as follows. In section 2 we will introduce the N2HDM together with our notation. Section 3 is dedicated to the description of the theoretical constraints that will be applied here for the first time in full scrutiny without any approximations on the N2HDM Higgs potential. Section 4 describes the parameter scan with the applied constraints. Section 5 is dedicated to the phenomenological analysis. Our conclusions are collected in section 6.

The N2HDM Higgs Sector
The N2HDM is based on the CP-conserving (or real) 2HDM with a softly broken Z 2 symmetry extended by a real singlet field Φ S . The extension of the 2HDM by real scalar singlet that does not acquire a VEV provides a viable DM candidate [24][25][26][27][28][29][30][31][33][34][35]37]. In [38] the phenomenology of the N2HDM with non-vanishing VEV for the singlet field has been discussed by applying some approximations. In particular, the possibility of a singlet admixture to the 125 GeV Higgs boson has been neglected. In the following no such assumptions will be imposed on the N2HDM potential. In terms of the two SU (2) L Higgs doublets Φ 1 and Φ 2 and the singlet field Φ S , the N2HDM potential is given by The first two lines describe the 2HDM part of the N2HDM while the last line contains the contribution of the singlet field Φ S . This potential is obtained by imposing two Z 2 symmetries on the scalar potential. The first one, called Z 2 , is the trivial generalisation of the usual 2HDM Z 2 symmetry to the N2HDM. It is softly broken by the term involving m 2 12 and can be extended to the Yukawa sector to guarantee the absence of tree-level Flavour Changing Neutral Currents (FCNC). The second Z 2 symmetry, Z 2 , is and is not explicitly broken. If Φ S does not acquire a VEV this second Z 2 symmetry will give rise to a conserved "darkness" quantum number and to the appearance of a dark matter candidate. If Φ S acquires a VEV this quantum number is no longer conserved and there is mixing among all CP-even neutral particles. This same behaviour is still possible if m 2 12 = 0, where both Z 2 symmetries in the potential are exact but need to be spontaneously broken. We will not consider such model further in this study. The two Z 2 quantum numbers assigned to the scalars in the model are shown in table 1. After EWSB the two doublet fields acquire the real VEVs v 1 and v 2 and the singlet field a real VEV v S . They can be parametrised as in terms of the charged complex fields φ + i (i = 1, 2) and the real neutral CP-even and CP-odd fields ρ I (I = 1, 2, S) and η i , respectively. Requiring the potential to be minimized at the VEV leads to three minimum conditions given by with Replacing the doublet and singlet fields in the Higgs potential by the parametrisations (2.4) the mass matrices in the gauge basis are obtained from the second derivatives with respect to the fields in the gauge basis. Due to charge and CP conservation the 7 × 7 mass matrix decomposes into three blocks: the 2 × 2 matrix for the charged Higgs bosons, the 2 × 2 matrix for the CP-odd fields and the 3 × 3 matrix for the CP-even states. Introducing a real singlet field with a VEV, the charged and pseudoscalar sectors of the model remain unchanged with respect to the 2HDM. Consequently, as in the 2HDM, the charged and pseudoscalar mass matrices can be diagonalised by the rotation matrix with t β defined as We have introduced the SM VEV v ≈ 246 GeV and the abbreviations sin x ≡ s x , cos x ≡ c x and tan x ≡ t x . This yields the massless charged and neutral would-be Goldstone bosons G ± and G 0 , the charged Higgs mass eigenstates H ± and the pseudoscalar mass eigenstate A.
Due to the additional real singlet field, the CP-even neutral sector of the N2HDM is changed with respect to the 2HDM. Instead of a 2 × 2 mass matrix we now have a 3 × 3 matrix. In the basis (ρ 1 , ρ 2 , ρ S ) it can be cast into the form for v, t β and v S . The mass matrix can be diagonalised by an orthogonal matrix R which we parametrise as in terms of the mixing angles α 1 to α 3 . Without loss of generality the angles can be chosen in the range (2.14) The matrix R rotates the interaction basis (ρ 1 , ρ 2 , ρ S ) into the physical mass eigenstates H 1 , H 2 and H 3 , and diagonalises the mass matrix M 2 scalar , We use a convention where the mass eigenstates are ordered by ascending mass as In total, the N2HDM is described by 12 independent real parameters. We choose as many parameters with physical meaning as possible. We use the minimisation conditions to trade m 2 11 , m 2 22 and m 2 S for the SM VEV v, t β and v S and replace the quartic couplings by the physical masses and the mixing angles. The soft Z 2 breaking parameter m 2 12 is kept as an independent parameter. Thus, we use the following set of input parameters In appendix A.1 we provide expressions for the quartic couplings in terms of these input parameters.
The singlet field ρ S does not directly couple to the SM particles. Therefore, any change in the tree-level Higgs couplings compared to the 2HDM is due to the mixing of the three neutral fields ρ I . Any coupling not involving the CP-even neutral Higgs bosons remains unchanged compared to the 2HDM and can be found e.g. in [23]. We now provide the couplings of the N2HDM Higgs bosons H i relevant for Higgs decays. We introduce the Feynman rules for the Higgs couplings H i to the massive gauge bosons V ≡ W, Z as where g H SM V V denotes the SM Higgs coupling factor. In terms of the gauge boson masses M W and M Z , the SU (2) L gauge coupling g and the Weinberg angle θ W it is given by We obtain for the effective couplings defined by Eq. (2.19). Replacing the R ij by their parametrisation in terms of the mixing angles yields the effective couplings in table 2.
In order to avoid tree-level FCNCs we extend the Z 2 symmetry (2.2) to the Yukawa sector. This leads to the same four types of doublet couplings to the fermions as in the 2HDM. We show these types in table 3. Consequently, the CP-even H i Yukawa couplings take the same form as the Yukawa couplings of the 2HDM. With the N2HDM Yukawa Lagrangian we obtain the effective coupling factors c(H i f f ) in terms of the mixing matrix elements R ij and the mixing angle β provided in table 4. Replacing the R ij by their parametrisation in terms of the α i this results in the effective coupling expressions given for type I and II in table 5.  The Feynman rule for the H i coupling to the pseudoscalar A and the Z boson is given by where g denotes the U (1) Y gauge coupling and p A and p H i , the four-momenta of the pseudoscalar and the H i , are both taken as incoming. The tilde over the coupling factor indicates that it is not an effective coupling in the sense that it is not normalized to a corresponding SM coupling, since there is no SM counterpart. The corresponding Feynman rule for the H i coupling to the charged pair H ± and W ∓ reads where p H ± denotes the four-momentum of H ± and again all momenta are taken as incoming. The coupling factorsc(H i V H) are provided in table 6.
The trilinear Higgs self-couplings relevant for the Higgs decays into a pair of lighter Higgs bosons are quite lengthy and deferred to appendix A.2. We have used these Feynman rules to implement the N2HDM in HDECAY v6.51 [23]. The resulting code N2HDECAY calculates all N2HDM Higgs boson decay widths and branching ratios including state-of-the-art higher order QCD corrections and off-shell decays. Electroweak corrections, which in contrast to the QCD corrections cannot be taken over from the SM, have been consistently neglected. The program can be downloaded from the url: https://itp.kit.edu/∼maggie/N2HDECAY In appendix B a short description of the program with a sample input and output file is provided.   We finally note that by letting α 2,3 → 0 and α 1 → α + π/2 the N2HDM approaches the limit of a 2HDM with an added decoupled singlet. In the 2DHM, the mixing angle α diagonalises the 2 × 2 mass matrix in the CP-even Higgs sector leading to the two CP-even mass eigenstates h and H, respectively. The shift by π/2 in the limit is necessary to match the usual 2HDM convention. Hence, (2.25)

Theoretical Constraints
In this section we investigate the conditions on the N2HDM imposed by theoretical considerations. These are the requirements for tree-level perturbative unitarity, that the vacuum is stable and that it is the global minimum of the scalar potential. For the first two requirements the corresponding conditions in the N2HDM can be derived from the literature. In the following we will summarize them before presenting the analysis of the stationary points of the N2HDM (cf. also [75]). The model, along with the theoretical conditions, has been implemented in ScannerS. This allows us to perform extensive scans in the parameter space of the N2HDM taking into account both the experimental and the theoretical constraints, as will be described in detail in section 4.

Tree-level Perturbative Unitarity
Tree-level perturbative unitarity is ensured by requiring that the eigenvalues of the 2 → 2 scalar scattering matrix are below an absolute upper value given by 8π [79]. It can be useful to impose a limit smaller than 8π at tree level to safeguard for possible large enhancements of the scalar couplings through higher order corrections. Following the procedure and notation of [79], we have calculated the full 2 → 2 scattering matrix of the fields in the gauge basis, The resulting matrix is block diagonal. The block matrices not containing ρ S have the eigenvalues (3.28) These are the same eigenvalues as found in the 2HDM. The new contributions due to the singlet field yield the eigenvalues and the eigenvalues a 1,2,3 , which are the real roots of the cubic polynomial Not all of the above eigenvalues are independent. As noted in [79], we have The conditions on f 1 , e 2 and f − can therefore be dropped, as they follow from the conditions on p 1 , e 1 and f + . Since λ 1 , λ 2 > 0 is necessary for the potential to be bounded from below (see Eqs. (3.50) and (3.51)) we obtain The resulting conditions for tree-level perturbative unitarity are thus given by where a 1,2,3 are the real roots of Eq. (3.37).

Boundedness from below
We consider the potential to be bounded from below in the strong sense, which means that the potential is required to be strictly positive as the fields approach infinitiy. The corresponding necessary and sufficient conditions have been given in [80] and translated to our notation in [75]. They depend on the discriminant The allowed region is given by with and

Global Minimum Conditions for the N2HDM Potential
In general, the minimum needs not be the global minimum, if the tunnelling time for the vacuum to tunnel into the global minimum [81,82] is larger than the age of universe. As the calculation of the tunnelling time is beyond the scope of this work, we do not discuss such metastable vacua and restrict to the stronger requirement that the vacuum is at the global minimum. While for the 2HDM it has been proven [83] that the existence of a normal minimum of the form Eq. (2.4) precludes the existence of a deeper charge-or CP-breaking minimum, this does not generalise to the N2HDM. Counter-examples that underline this statement can be found in the appendix of [75]. For the analysis of the global N2HDM minimum we therefore have to include the possibility of CP-and charge-breaking minima. We consider the most general constant field configuration, where all fields are real, Here we have already exploited the SU (2) L × U (1) Y gauge symmetry to eliminate four degrees of freedom. Any other possible constant field configuration of the N2HDM can be projected onto this one through a gauge transformation. By v cb we denote the charge-breaking and by v cp the CP-breaking constant fields. In the following we will refer to the constant fields as VEVs although this is technically only correct if the configuration describes a minimum of the scalar potential. By expanding the field configuration Eq. (3.52) in the potential one observes a set of Z 2 symmetries for the real fields we are using, and we can thus choose, without loss of generality, all VEVs except for v 2 to be positive. We then proceed to find all possible stationary points of the N2HDM. This detailed analysis is presented in appendix C.
We want the minimum to conserve both electric charge and CP and to give rise to three CP-even massive scalars. From the vacuum structure point of view this means that v 1 = 0, v 2 = 0 and v S = 0 while v cb = 0 and v cp = 0 at the chosen minimum. In order to ensure that this is the global minimum we proceed as follows: • We choose the model parameters such that there is a minimum with • Using the stationarity conditions presented in appendix C we look for all other possible stationary points of the potential with this set of parameters.
• We discard all sets of parameters for which we find a stationary point below the minimum.
This procedure leads to a global minimum that is CP-conserving, preserves electric charge and allows for the singlet to mix with the CP-even scalars from the doublets. We note here that, in the N2HDM, the classification of the possible types of minima, in terms of preservation or breaking of electric charge or CP, follows that of the 2HDM regardless of the value of v S -see also appendix C.

The Parameter Scan
In order to perform phenomenological analyses we need viable parameter points, i.e. points in agreement with theoretical and experimental constraints. To obtain these points we use the program ScannerS to perform extensive scans in the N2HDM parameter space and check for compatibility with the constraints. We denote the discovered SM-like Higgs boson with a mass of [4] m h 125 = 125.09 GeV (4.53) by h 125 . In the following, we exclude all parameter configurations where this Higgs signal is built up by multiple resonances by demanding the mass window m h 125 ± 5 GeV to be free of any Higgs bosons except for h 125 . Furthermore, we do not include electroweak corrections in the parameter scans nor in the analysis, as they are not (entirely) available for all observables and cannot be taken over from the SM.
We check for the theoretical constraints on the N2HDM at tree level as described in section 3. Tree-level perturbative unitarity is verified by using Eqs. (3.43)- (3.47). This method yields a shorter run-time than the model-independent numeric check implemented in ScannerS. We checked that both methods lead to the same results. Equations (3.50) and (3.51) are used to guarantee that the potential is bounded from below. The vacuum state found by ScannerS is required to be the global minimum, otherwise it is rejected. As described in 3.3, the check is performed by comparing the value of the scalar potential at the ScannerS vacuum with the values at all of the stationary points.
Many of the experimental constraints applied on the 2HDM also hold for the N2HDM. The constraints on R b [84,85] and B → X s γ [85][86][87][88] are only sensitive to the charged Higgs boson so that the 2HDM calculation and the resulting 2σ exclusion bounds in the m H ± − t β plane can be taken over to the N2HDM. Note that the latest calculation [88] enforces m H ± > 480 GeV (4.54) in the type II and lepton specific 2HDM. In the type I model on the other hand the bound is much weaker and more strongly dependent on tan β. The oblique parameters S, T and U are calculated with the general formulae in [89,90], and 2σ compatibility with the SM fit [91] including the full correlations is demanded.
The N2HDM must comply with the LHC Higgs data. This requires one scalar state to match the observed signal rates for a Higgs boson of about 125 GeV. Furthermore, the remaining Higgs bosons must be consistent with the exclusion bounds from the collider searches at Tevatron, LEP and LHC, where the strongest constraints arise from the LHC Run 1 data. ScannerS provides an interface with HiggsBounds v4.3.1 [92][93][94] which we use to check for agreement with all 2σ exclusion limits from LEP, Tevatron and LHC Higgs searches. The required input for HiggsBounds are the cross section ratios (relative to a SM Higgs boson of the same mass) of the different production modes, the branching ratios and the total widths for all scalars. We compute the latter two with the program N2HDECAY. The production cross sections through gluon fusion (ggF) and b-quark fusion (bbF) are obtained at next-to-next-to-leading order (NNLO) QCD from SusHi v1.6.0 [95,96] which is interfaced with ScannerS. 3 The remaining cross section ratios at leading order (and also at higher order in QCD) are given by the effective couplings squared. For example for the production in association with a fermion pair it is c(H i ff ) 2 (Table 4), and for the gauge boson mediated cross sections (vector boson fusion and associated production with Compatibility of the discovered Higgs signal with h 125 is checked by using the individual signal strengths fit of Ref. [97]. The needed decay widths and branching ratios are taken from N2HDECAY. The fermion initiated cross section normalized to the SM, is obtained with the NNLO QCD cross sections taken from SusHi. In the normalization we neglect the bbF cross section, which in the SM is very small compared to gluon fusion. The production through vector boson fusion (VBF) or through associated production with a vector boson (VH) normalized to the SM, µ V , is given by where H i is identified with h 125 . The QCD corrections to massive gauge boson-mediated production cross sections cancel upon normalization to the SM. The properties of the h 125 are checked against the six fit values of given in [97], with µ xx F defined as For H i ≡ h 125 we require agreement with the fit results of [97] at the 2 × 1σ level.
In the numerical analysis we will show results for type I and type II N2HDM models. For the scan with the input parameters from Eq. (2.18) we fix v to the SM value and choose t β in the range 0.25 ≤ t β ≤ 35 . (4.59) As the lower bound on t β from the R b measurement is stronger than the lower bound in Eq. (4.59), the latter has no influence on the physical parameter points. We transform the mixing matrix generated by ScannerS to the parametrisation of Eq. (2.13) such that the mixing angles are allowed to vary in the ranges (4.60) We identify one of the neutral Higgs bosons H i with h 125 and allow the remaining neutral Higgs bosons to have masses within In the type II model, the charged Higgs mass is chosen in the range  The condition m 2 12 > 0 is found to be necessary for the minimum to be the global minimum of the scalar potential.

Phenomenological Analysis
We start with the investigation of the parameter distributions and from now on denote the lighter of the two non-h 125 CP-even Higgs bosons by H ↓ and the heavier one by H ↑ .
The inspection of the mass distributions resulting from our scan for type II shows that the masses can take all values between 30 GeV and 1 TeV. Furthermore, we note that it is possible that both H ↓ and H ↑ are lighter than h 125 and also that H ↓ and A have masses below 125 GeV. Due to the lower bound m H ± ≥ 480 GeV and the constraints from the EWPT, which force at least one of the non-SM-like neutral Higgs bosons to have mass close to the charged Higgs mass, there is no scenario with all non-SM-like neutral Higgs masses below 125 GeV.
In type I, overall we have lighter Higgs spectra because of the much weaker lower bound on the charged Higgs mass of m H ± > ∼ 80 GeV. Consequently, here we can also have situations where h 125 is the heaviest Higgs boson in the spectrum.

The wrong-sign Yukawa Coupling Regime
The wrong-sign Yukawa couplings regime, which was discussed in [98][99][100] for the CP-conserving 2HDM, is the parameter region where the coupling of the h 125 to the massive gauge bosons is of opposite sign with respect to the coupling to fermions, c(h 125 ff ). This region is in contrast to the SM case where both couplings have the same sign, which can have interesting phenomenological consequences such as the non-decoupling of heavy particles [98,101]. This region is not excluded by the experimental data for a type II model with an opposite sign Yukawa coupling to downtype fermions. For the coupling to up-type fermions this wrong-sign scenario is only realized for tan β < 1 which is excluded. In the type I N2HDM the Yukawa couplings to up-and down-type fermions are the same so that the wrong-sign coupling regime cannot be realized for any of the quark types, as it requires tan β < 1. We investigate the extent to which this scenario can be realized in the N2HDM type II. In the 2HDM the wrong-sign regime is obtained for parameter values where sin α > 0, while the correct-sign regime is realized in the opposite case. In order to match the 2HDM description, cf. Eq. (2.25), we use the condition sgn(c(h 125 V V )) · sin(α 1 − π 2 ) > 0 (5.66) for the wrong-sign limit. Figure 1 displays tan β versus sgn(c(h 125 V V )) · sin(α 1 − π 2 ) for all parameter points from our N2HDM type II scan that survive the imposed constraints and where the SM-like Higgs is given by H 1 = h 125 . 4 The colour code quantifies the singlet admixture Σ h 125 of the SM-like Higgs boson h 125 . We define the singlet admixture Σ h 125 through i.e. the absolute value squared of the mixing matrix element describing the mixing of the singlet field with the SM-like Higgs state. In the right plot we have inverted the colour ordering. In the left panel of Fig. 1 we observe a large number of points in the 2HDM limit, i.e. with small singlet admixture. Such points are distributed in two branches with a shape that agrees with the 2HDM. This allows to verify the identification of the left branch with the correct-sign regime and of the right branch with the wrong-sign regime. The inverted colour ordering in Fig. 1 (right) allows us to investigate the repartition of the singlet admixture over the two limiting  cases. Overall we see that the singlet admixture can be considerable. In the wrong-sign regime it reaches up to about 30% while in the correct-sign regime it can even be as large as 55%. The points with the largest singlet admixture can be found for small values of tan β. In the following analysis we will comment further on the interesting phenomenology of the wrong-sign regime.
Constraining the wrong-sign regime An important question to ask is to which extent will the collection of more precise data, obtained at the LHC Run II and in the high-luminosity phase, be able to constrain the N2HDM parameter space and in particular the wrong-sign regime. In Fig. 2 we show again the allowed region in the tan β versus sgn(c(h 125 V V )) · sin(α 1 − π/2) projection of the parameter space, where the grey points respect all theoretical and experimental constraints, in particular the reduced signal strengths from the six-parameter fit of Ref. [97] for the SM-like Higgs (here H 1 ≡ h 125 ). We then successively constrain the µ-values further, by assuming that future more precise measurements are able to achieve a precision of 5%, with a central value of 1. Hence The plot shows that the correct-sign regime given by the left branch is most strongly constrained by the τ τ final state. While the wrong-sign limit is also very sensitive to this observable its compatibility with the data is fundamentally different as there are no black nor blue points in the right branch. This behaviour is represented in different form in Fig. 3, in which we have plotted µ V /µ F versus µ γγ . The depicted yellow areas show the points from the parameter scan compatible with all constraints in the correct-sign regime, for the N2HDM (left plot) and, for comparison, also for the 2HDM (right). 5 The pink points represent the wrong-sign regime. Here and from now on h 125 can be any of the CP-even, neutral Higgs bosons (i.e. any H i in the N2HDM and h or H in the 2HDM). The white triangle denotes the SM result. The yellow and pink regions completely overlap in the N2HDM in contrast to the 2HDM. Here we have less parameters and hence less freedom to reach compatibility with the constraints than in the N2HDM so that the yellow area is more restricted. 6 The allowed area of the wrong-sign regime, on the other hand, is the same in both models. The reason is that in the N2HDM the singlet admixture can at most reduce the Higgs couplings to SM particles and hence the µ-values. As the wrong-sign regime in the 2HDM is already touching the lower bounds in the presented µvalues the N2HDM cannot add anything new to this region. From these figures we immediately infer that in the wrong-sign regime the ratio µ V /µ F cannot reach 1, which explains the missing blue and black points in Fig. 2. The measurement of µ V /µ F is hence a powerful observable to constrain the wrong-sign regime with values of µ V /µ F > ∼ 0.9, excluding this scenario. We note that the pink region in Fig. 3 is rather insensitive to an increase in the precision of µ V V to 5% around 1.
In Fig. 4 the pink points show the wrong-sign regime in the µ γγ versus µ τ τ plane. The left panel is for the N2HDM, with h 125 being any of the H i , whereas the right panel is for the 2HDM, with 7 h 125 ≡ h. We show on top, in green, a further restriction of the sample assuming that future measurements can constrain the µ V V , V = Z, value to 5% around the SM value. Again we observe that the 2HDM area, with a smaller number of parameters, is more constrained than the N2HDM, with the upper bounds in both models being about the same. More importantly, we note that the increase in the precision of µ V V forces the reduced photonic rate to be below about 0.9. 8 Thus, the the wrong-sign regime can be excluded by increasing the precision in the µ V V measurement and observing µ γγ > ∼ 0.9. The outliers in Fig. 4 are points where h 125 has a substantial decay width into an (off-shell) pair of light Higgs bosons. The resulting increase of the total width reduces the branching ratio into τ τ and thus µ τ τ . In the majority of the scenarios the light Higgs boson is the pseudoscalar A. If, however, h 125 is not the lightest CP-even neutral Higgs boson, decay widths of h 125 → H 1 H 1 (h 125 → hh in the 2HDM) can also be substantial.

Phenomenology of the h 125 Singlet Admixture
Type II N2HDM: The large number of parameters in the N2HDM allows for considerably non-standard properties in the phenomenology of the SM-like Higgs boson. In particular, in type II, which we discuss first, significant singlet admixtures of up to 55% are still compatible with the LHC Higgs data. This can be inferred from Figs. 5 and 6. Figure 5 displays the correlation between pairs of the effective couplings squared of h 125 to the SM particles. The left plot shows the coupling to top quarks versus the coupling to bottom quarks, and the right plot shows the coupling to the massive gauge bosons versus the coupling to bottom quarks. We remind that h 125 can be any of the neutral CP-even Higgs bosons. The influence of the singlet admixture is quantified by the colour code. In Fig. 6 (left) we show the reduced signal rate in τ final states, µ τ τ , versus the one into massive vector bosons, µ V V . In Fig. 6 (right) we plot the ratio µ V /µ F of the vector boson induced production over the fermionic production, each normalized to the SM, versus the photon final state signal strength µ γγ . The white triangle indicates the SM values of the signal strengths. The dashed lines are the experimental limits on the respective signal strengths. They show that the N2HDM parameter space is constrained by the upper and lower limits on µ V V and µ γγ and the lower limits on µ τ τ and µ V /µ F , respectively. As can be inferred from Fig. 6 (left) enhanced rates in the τ final state are still allowed by the experimental data where the largest enhancement is reached for small admixture, i.e. in the 2HDM-like regions. The area of the enhanced µ τ τ can be divided in three regions that shall be explained separately. The largest enhancement of up to 40% is obtained for simultaneously enhanced µ V V . The enhancement is due to the production mechanism and corresponds to the enhanced couplings to top quarks in Fig. 5 (left) while the involved decays remain SM-like. This is confirmed by Fig. 7 (left), which shows the value of µ τ τ in the plane of the effective couplings squared, c 2 (h 125 tt) and c 2 (h h125 bb), for parameter points in the correct-sign regime. The largest µ τ τ , given by the yellow points, are found for large effective couplings to top quarks. The 2HDM-like region in Fig. 6, where we have enhanced µ τ τ values but reduced µ V V , is due to enhanced effective couplings to τ leptons and b quarks, i.e. the spikes in Fig. 5 and Fig. 7 (left), respectively, where µ τ τ reaches values of up to about 1.4. On the other hand, the effective coupling to gauge bosons cannot exceed one, so that overall the branching ratio into gauge boson is reduced in favor of BR(h 125 → τ τ ). Finally, the points with µ τ τ > 1 and µ V V ≈ 1 are located in the wrong-sign regime and have reduced couplings to the gauge bosons. In Fig. 5 (right) these are the points below the dashed line and isolated from the bulk of the points. Figure 7 (right) shows the same coupling values squared, but now only for points in the wrong-sign limit, where the area with enhanced µ τ τ can easily be identified. It is the resulting reduced decay width into V V that increases the branching ratio into τ τ and thus µ τ τ . Note that in Fig. 5 (left) the sharp lines departing from the SM point to smaller and larger values of c 2 (h 125 bb) are again due to the unitarity of the mixing matrix.
While the enhanced µ τ τ is a feature of the 2HDM-like regions, singlet admixtures as large as 55% can be compatible with the current data. Interestingly the best measured quantities µ V V and µ γγ are not the ones with the highest constraining power on the singlet admixtures. A value of µ V V = 1 still allows for Σ h 125 of up to 50%, and a value of µ γγ = 1 can be compatible with Σ h 125 values of up to about 40%. A measurement of µ τ τ ≈ 1 on the other hand constrains Σ h 125 to be below about 25%. And a measurement of µ V /µ F ≥ 1 enforces Σ h 125 < ∼ 20%. This behaviour can be understood by inspecting the involved couplings individually. Overall, high singlet admixtures induce reduced couplings. However, Fig. 5 shows that the effective coupling to b-quarks is reduced more strongly than that to top-quarks, where the latter dominates gluon fusion production. Furthermore, among the physical points, a large singlet admixture leads to  the effective gauge coupling being larger than the effective bottom coupling. This means that the total width is reduced due to the strongly reduced Γ(h 125 → bb), which dominates the total width in the SM case. In turn, this increases the BR(h 125 → V V ) and BR(h 125 → γγ) strongly enough to balance the reduced gluon fusion production cross section and the reduced partial widths in these channels, so that overall SM-values of µ V V and µ γγ can be compatible with large singlet admixtures. The decay width into τ leptons, in contrast, gets rescaled by the same coupling as the bb channel in the type II model, so that the µ τ τ is reduced for large singlet admixtures. Since with increasing singlet admixture c(h 125 V V ) is more strongly reduced than c(h 125 tt) the value µ V /µ F is reduced for non-zero singlet admixtures. As the coupling to gauge bosons reaches at most 1, an enhanced µ V /µ F can only be due to a smaller coupling of h 125 to tt. A value of c(h 125 V V ) close to 1 with a simultaneously reduced c(h 125 tt) is only possible for small singlet admixtures. The corresponding bulge in Fig. 6 (right) contains points for which all other µ-values are close to their lower experimental boundaries. They are characterised by large couplings of the SM-like Higgs to a charged Higgs pair, which enters the loop-induced Higgs decay into photons and keeps µ γγ above its lower experimental boundary.
We conclude our discussion of the type II N2HDM by displaying in Fig. 8 the allowed N2HDM parameter region in the tan β versus Σ h 125 plane (grey points) and its subsequent restriction by more precise measurements, with the same colour code as defined in Eq. (5.68). It reflects our findings that the singlet admixture is most powerfully constrained by a precise measurement of µ τ τ while being much less sensitive to the remaining µ-values. The restriction power depends on the value of tan β. For medium values of tan β a singlet admixture of up to 20% is still compatible with a 5% measurement of µ τ τ whereas for small values of tan β even up to 37% is allowed. Only the simultaneous measurement of all µ-values within 5% around the SM value Type I N2HDM: We now turn to the discussion of the N2HDM type I. In this case the doublet Φ 2 couples to both up-and down-type quarks. Consequently, any change in the coupling to the top quark is induced also in the bottom quark coupling and vice versa. This restricts the available freedom in the choice of the parameters of the N2HDM so that, overall, less pronounced deviations from the SM or from the pure 2HDM are expected. This can be inferred from Figs. 9 and 10, which show, respectively, the singlet admixture of h 125 as a function of the reduced signal strengths and the effective couplings. In constrast to the type II N2HDM the singlet admixture can reach at most 25%. Figure 9 shows the most constraining signal strenghts, i.e. µ τ τ versus µ V V in the left plot and µ V /µ F versus µ γγ in the right plot. From this figure it can be inferred that the allowed areas in these two planes are more reduced than in type II. While µ γγ is delimited by the present LHC data, only the lower bound of µ V V is restricted due to the LHC fit value, while the upper µ V V bound and the boundaries of µ τ τ and µ V /µ F of the allowed areas are already well below the restrictions set by the LHC data. The highest singlet admixtures come with reduced signal strengths while the ratio µ V /µ F ≈ 1. This is in accordance with Fig. 10, which shows the effective couplings squared, c(h 125 V V ) 2 versus c 2 (h 125 bb) (= c 2 (h 125 tt)) together with the singlet admixture. Both effective couplings are reduced almost in parallel with rising Σ h 125 so that µ V /µ F ∼ c(h 125 V V ) 2 /c(h 125 f f ) 2 is close to one for large singlet admixtures.
The SM-like dark-blue boundary in Fig. 9 (left) with enhanced µ τ τ corresponds to the right boundary of the area in Fig. 10 where c 2 (h 125 bb) is enhanced and c 2 (h 125 V V ) is reduced. The simultaneous enhancement of the Higgs coupling to the bottom quarks and τ leptons leaves the branching ratio into τ 's unchanged, but the enhanced Higgs couplings to fermions enhance the main production mechanism so that overall µ τ τ is enhanced. Applying the same reasoning it is clear that the dark-blue boundary with reduced µ τ τ values corresponds to the left dark-blue boundary in Fig. 10. This is also the area in Fig. 9 (right) corresponding to enhanced µ V /µ F values for µ γγ < ∼ 1.3 in the 2HDM-like region (blue area). The reduced coupling to the bottom quarks reduces the dominating decay into bb. The stronger reduction of c(h 125 bb) as compared to c(h 125 V V ) overcompensates the reduced Γ(h 125 → V V ) so that overall BR(h 125 → V V ) is enhanced and makes up for the reduction of the production cross section due to the smaller couplings to fermions. Simultaneously, this also ensures that µ V /µ F is enhanced. An analogous reasoning allows to identify the dark blue lower µ V /µ F values for µ γγ < ∼ 1.3 with the right  boundary in Fig. 10. Figure 11 allows to analyse by which measurement the singlet admixture can be most effectively constrained. The colour code has been defined in Eq. (5.68). As can already be inferred from Fig. 9 the singlet admixture is mostly restricted by the precise measurement of µ ZZ , down to about 7.5% for a 5% accuracy in µ ZZ . The simultaneous measurement of all µ-values with a precision of 5% hardly improves on this constraint. When comparing with the type II case, we can conclude that the structure of the Yukawa couplings decides which final state in the Higgs data is the most effective one in constraining the N2HDM, i.e. its singlet admixture.

Conclusions
In this paper we have investigated the N2HDM, which is based on the CP-conserving 2HDM extended by a real scalar singlet field. It combines a parameter space which is larger than the 2HDM with a greater freedom in the choice of the parameters (compared to singlet extended supersymmetric models as e.g. the next-to-minimal supersymmetric model). This allows for an interesting phenomenology that is still compatible with the experimental data. Thus the Higgs couplings can carry a substantial singlet admixture. However, in order to be able to determine the allowed parameter space and thereby perform meaningful phenomenological analyses, the investigation of the constraints on the model had to be put on solid ground. In this paper we have performed a thorough analysis of the theoretical constraints on the N2HDM Higgs potential. First, we have collected from the literature the formulae for the N2HDM that test tree-level perturbative unitarity and stability of the vacuum. Moreover, for the first time, we have presented a detailed analysis of all the stationary points of the potential to obtain its global minimum. The model, together with the theoretical constraints, has been implemented in ScannerS. For the test of the experimental constraints the necessary branching ratios were obtained with the new program N2HDECAY. We have written this code based on HDECAY to provide the N2HDM branching ratios and total widths including the state-of-the-art higher order QCD corrections and off-shell decays. With these preparatory works completed we were then in the position to subject the N2HDM to critical theoretical and experimental scrutiny in the second part of the article.
Taking into account the theoretical and experimental constraints in the N2HDM, we performed a parameter space scan and analysed the properties of the allowed regions. In the type II N2HDM substantial singlet admixtures of up to 55% were found to be compatible with the data. It turned out that the most precisely measured quantities, µ V V and µ γγ , are not the most effective ones in constraining the singlet admixture, but instead it is the µ τ τ observable. In the type I N2HDM, on the other hand, the parameter space is more constrained due to the universal fermion coupling to up-and down-type quarks and hence the allowed singlet admixture remains below about 25%. In the future, in this case, the singlet admixture will be most strongly constrained by the precise measurement of µ ZZ .
Like in the 2HDM we find a wrong-sign regime in the N2HDM. While overall the allowed parameter space is larger in the N2HDM compared to the 2HDM, the wrong-sign regimes in both models are comparable. In this regime µ V /µ F is found to be well below 1, so that a measurement of a value near 1 excludes this scenario. Moreover, a future measurement of µ V V with a precision of 5% around the SM value and the observation of µ γγ > ∼ 0.9 eliminates the wrong-sign regime. These findings are consistent with the observations in the 2HDM.
With the analysis tools provided in this study, the next natural step would be to compare the N2HDM with other extended Higgs sectors with a similar theoretical structure and Higgs spectrum, to investigate which observables allow to distinguish (or to exclude) the models.
factor i can be cast into the form where ijk denotes the totally antisymmetric tensor with 123 = 1. Note, that in Eqs. (A.10) to (A.14) there is no summation over repeated indices. The sums of different powers of R ij arise from simplifications that exploit the orthogonality of the mixing matrix. The employed formula reads The matrix / R ij is the submatrix formed by deleting the i-th row and the j-th column from R. The indices i and j take any values in {1, 2, 3}.

B The Fortran Code N2HDECAY
The code N2HDECAY is the N2HDM implementation in the program HDECAY, written in Fortran77. It is based on HDECAY v6.51. The code is completely self-contained. All changes related to the N2HDM have been implemented in the main file n2hdecay.f. Further linked routines have been taken over from the original HDECAY code. The implemented decay widths include the most important state-of-the-art higher order QCD corrections and the important off-shell decays. They can be taken over from the SM and the minimal supersymmetric extension (MSSM), respectively, for which HDECAY was originally designed. The electroweak corrections have been consistently turned off as they cannot be taken over from the available corrections in the SM and/or MSSM.
The code is compiled with the makefile by typing make. This produces an executable file called run. Typing run executes the program, which calculates the branching ratios and total widths that are written out together with the mass of the decaying Higgs boson. The output files are called br.X N2HDM y. Here X=H1, H2, H3, A, H+ denotes the decaying Higgs particle. 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. In the following we present the example of an output file as obtained from the above input file. The produced output in the four output files br.H3 N2HDM y for the heaviest neutral Higgs boson is given by

MH3
BB TAU TAU  MU MU  SS  CC  TT  ---- All files necessary for the program can be downloaded at the url: http://www.itp.kit.edu/∼maggie/N2HDECAY The webpage contains a short explanation of the program and information on updates and modifications of the program. Furthermore, sample output files can be found for a given input.

C Global Minimum Conditions
We start by recalling that, up to gauge symmetries, the most general constant field configuration, where all fields are real, is The subscripts cp and cb refer to the case where CP or charge are spontaneously broken in scenarios where also both v 1 and v 2 are non-zero. In the case where only v cp and/or v cb are non-zero there is in fact no CP or charge breaking. One way to see this is by noting that such a configuration is continuously connected by a gauge transformation to the case where only v 2 = 0 (where it is more clear that there is no CP or charge breaking). Such a gauge transformation amounts to a redefinition of the charge operator by a rotation.
In order to find all possible minima we consider the stationarity conditions for the VEVs, where λ 345 has been defined in Eq. (2.8) and The derivatives with respect to the degrees of freedom, which have been removed through a gauge transformation in Eq. (C.16), contribute with three further conditions From Eqs. (C. 19) and (C.20) we infer that except for the special case the VEVs v cb and v cp cannot be simultaneously non-zero. Eq. (C.24) is therefore always trivially satisfied. From Eqs. (C.17) and (C.18) we conclude that We further observe from Eqs. (C.23) and (C.25) that The configurations forbidden by non-vanishing m 2 12 are inert stationary points where the whole doublet VEV can be brought into one Higgs doublet through a basis transformation. Models with inert Higgs doublets show a very different phenomenology and we only consider points with

C.1 2HDM-like stationary points
The three cases of the 2HDM-like stationary points are obtained setting v S = 0. These are the stationary points of a 2HDM potential with the same parameters m 2 11 , m 2 22 , m 2 12 and λ 1−5 .
Case I: This case with a CP and charge conserving minimum is the most complicated one. We start by rewriting the minimum conditions in terms of Case I IIa IIb sI sIIa sIIb s v 1 1 1 Without loss of generality all VEVs except v 2 can be chosen positive due to the Z 2 symmetry. (C.31) The resulting system of equations is used to eliminate v leading to a single quartic equation for sin 2 δ, The ∓ signs correspond to the two possible signs of sin δ in the interval Eq. (C.31). The values for v 1 and v 2 are then obtained from Eq. (C.30) where both possible signs of v 2 need to be considered. Altogether this yields up to 16 possible solutions for v 1 and v 2 given by the four solutions for sin 2 δ times two for the sign of sin δ times two for the sign of v 2 . It should be noted only two of these solutions are independent as shown in [102,103].
Case II: The system of minimum conditions obtained here can be solved analytically both for case IIa and IIb. The formula for the value of the potential at these points can be cast into the form (C.41)

C.2 Stationary points with a singlet VEV
The stationary points for a non-vanishing singlet VEV v S = 0, as possible in the N2HDM, are obtained analogously to the previous cases.

Case sI:
Again the case with non-vanishing v 1 and v 2 but zero v cp and v cb is the most complicated one and leads to a quartic equation for sin 2 δ, with δ defined in Eq.   If λ 4 = λ 5 all VEVs can be simultaneously non-zero. In this case the scalar potential takes the stationary value .
(C. 55) In order to check if our minimum is the global one we compare the value of the scalar potential at our minimum with the values of the potential at all the stationary points listed above. In practice this means that we need to compare with the up to five different analytically known values of cases (s)II and case s and with the numerical solutions of case (s)I.