Vacuum Instabilities in the N2HDM

The Higgs sector of the Next-to-Minimal Two-Higgs-Doublet Model (N2HDM) is obtained from the Two-Higgs-Doublet Model (2HDM) containing two complex Higgs doublets, by adding a real singlet field. In this paper, we analyse the vacuum structure of the N2HDM with respect to the possibility of vacuum instabilities. We show that while one type of charge- and CP-preserving vacuum cannot coexist with deeper charge- or CP-breaking minima, there is another type of vacuum whose stability is endangered by the possible occurrence of deeper charge- and CP-breaking minima. Analytical expressions relating the depth of different vacua are deduced. Parameter scans of the model are carried out that illustrate the regions of parameter space here the vacuum is either stable or metastable as well as the regions where tunnelling to deeper vacua gives rise to a too short lifetime of the vacuum. Taking other experimental and theoretical constraints into account, we find that the vacuum stability constraints have an important impact on the phenomenology of the N2HDM.


INTRODUCTION
The Higgs mechanism [1][2][3][4][5] has been introduced to generate particle masses without violating gauge symmetries. It is based on a sufficiently stable vacuum with non-zero vacuum expectation value v ≈ 246 GeV. Within the Standard Model (SM), the stability of the electroweak (EW) vacuum is guaranteed at lowest order as a consequence of the postulated form of the Higgs potential. Through higher-order corrections, the stability of the EW vacuum becomes intimately related also to the other Standard Model (SM) parameters [6,7], in particular the top quark mass. When extrapolating the SM to high energy scales it turns out that the EW vacuum is metastable for scales larger than about 10 10 GeV, which means that for these scales the vacuum is no longer absolutely stable but has a lifetime that is longer than the age of the universe.
Extensions beyond the SM (BSM) typically introduce new additional scalar degrees of freedom. While the loop contributions of these scalar particles may counteract the impact of the top quark loop, the presence of the additional scalars modifies the structure of the Higgs potential such that additional vacua can occur that are different from the one related to the correct EW symmetry breaking (EWSB). There can be vacua that break the CP symmetry (CP breaking) or the conservation of electric charge (charge breaking), in supersymmetric models even color breaking minima can occur. Moreover, there is the possibility of a second EW minimum but with a wrong vacuum expectation value (VEV) v = 246 GeV, as for example in the 2-Higgs Doublet Model (2HDM) where this situation was named "panic vacuum" [8][9][10][11][12]. In case an additional vacuum is deeper than the EW one then tunnelling can occur into dangerous non-physical vacuum configurations [13]. If this happens at time scales beyond the lifetime of the universe the EW vacuum is considered as metastable. However, parameter regions giving rise to faster vacuum decays are regarded as unphysical and should be excluded. Hence the requirement of a stable EW vacuum at cosmological time scales has immediate consequences for the allowed parameter space of the models. A thorough analysis of their vacuum structure is therefore crucial to correctly identify the allowed parameter space and consequently make appropriate predictions for observables and signatures for the experimental studies.
The analysis of the 2HDM [8,9,14,15] has shown that if a "normal" vacuum exists, i.e. a vacuum that is EW breaking but charge and CP conserving, any stationary point that is charge or CP breaking is necessarily a saddle point that lies above the normal minimum. There is also the possibility to have a second normal minimum but with the wrong VEV, i.e. a panic vacuum state. The Inert 2HDM, a 2HDM with an exact Z 2 symmetry, can have two types -Inert and Inert-like -of minima which can coexist with one another for certain parameter relations. The one-loop study, however, has shown [16] that the parameter regions where this is the case can change at loop level. The analysis of the possibility of a strong first order phase transition in the context of a CP-conserving and CP-violating 2HDM conducted in [17,18] revealed as a side product that the allowed minima at leading order do not necessarily lead to stable physical configurations at next-to-leading order (NLO) and vice versa. The developed code BSMPT [19] allows for studies of the vacuum structure at NLO (at zero and at finite temperature) of arbitrary user-defined BSM extensions. This is also the case for Vevacious [20,21], designed for general BSM models, including one-loop and temperature effects. Recently, members of this collaboration have presented an approach at leading order for an efficient and reliable evaluation of the constraints from vacuum stability and applied it to the minimal supersymmetric extension of the SM (MSSM) [22]. As shown in [23,24] and also discussed in [22], for calculations of the vacuum decay lifetime the loop-corrected effective potential in general does not correspond to a consistent perturbative expansion. A first analysis of the vacuum structure of the N2HDM has been carried out by some of the present authors in [25]. The N2HDM, which is obtained upon extension of the 2HDM with a real singlet field (which may acquire a VEV), was shown to exhibit a different vacuum structure than the 2HDM. Thus e.g. charge or CP breaking minima deeper than the normal minimum can exist.
In this paper we perform a detailed analysis of the vacuum structure of the N2HDM. We classify the different possible types of vacua and derive analytical expressions for the comparison of the values of the potential at minima of different nature. In contrast to [25], where a general phenomenological analysis of the N2HDM was performed (and where parts of our results have been presented in a numerical approach), we concentrate here on the vacuum structure itself and its implications on the model. By applying the method of [22], we investigate the requirement of a sufficiently stable physical minimum on the allowed parameter range. In particular, we investigate here for the first time the impact of the N2HDM vacuum structure on the phenomenology of the model. Moreover, we discuss the importance of including parameter regions with a metastable vacuum in phenomenological analyses in order to avoid incorrect conclusions on the viability of parameter space regions. Our study thus makes important new contributions to properly constraining the viable parameter space taking into account the theoretical constraints from the requirement of a stable vacuum.
The paper is organised as follows. In Sec. 2 we introduce the model and the different types of possible minima. It contains the detailed analytical analysis of the stability of the different minima. Section 3 is dedicated to the numerical analysis of the vacuum structure of the N2HDM. We describe our parameter scan and the method that we apply in order to identify the regions where the vacuum is stable or metastable. Subsequently we present and discuss our numerical results. We conclude in Section 4. The appendix contains a derivation that is used in our determination of the nature of the stationary points.

THE MODEL AND POSSIBLE MINIMA
The particle content of the N2HDM is identical to the one of the 2HDM in the fermionic and gauge sectors, but includes an extra real singlet scalar field, Φ S . To reduce the large number of parameters of the scalar potential, and to allow for the possibility of interesting phenomenology, such as dark matter, three discrete symmetries are imposed: (a) a Z 2 symmetry in which one of the doublets is affected by a sign change, Φ 1 → Φ 1 , Φ 2 → −Φ 2 and Φ S → Φ S ; (b) another Z 2 symmetry which leaves the doublets unchanged but changes the sign of the singlet, Φ 1 → Φ 1 , Φ 2 → Φ 2 and Φ S → −Φ S ; (c) the standard CP symmetry, Φ 1 → Φ * 1 and Φ 2 → Φ * 2 -since the singlet is real, the CP transformation does not affect it. After imposing these symmetries only terms quadratic and quartic in the fields are allowed and the most general scalar potential is given by where all parameters in the potential are real. We allow for the Z 2 symmetry (a) to be softly broken by the m 2 12 term. The theory obviously also includes fermions, and the Yukawa Lagrangian will depend on the choices made to extend the discrete symmetries imposed upon the scalar sector to the fermion one. Due to gauge invariance the singlet field Φ S couples to neither fermions nor gauge bosons. Therefore, the Yukawa Lagrangian will have four different possible forms, identical to the different types of 2HDM Yukawa Lagrangians. All of the four different possibilities lead to flavour conservation in scalar interactions. One of the possibilities (achieved if all right fermion fields change sign under the first Z 2 symmetry defined above) is a Type-I model, in which all fermions only couple to the doublet Φ 2 , and the Yukawa Lagrangian for the third generation is given by with Q L and L L denoting the left-handed quark and lepton doublets, and t R , b R and τ R the right-handed top, bottom and tau singlets. The remaining three Yukawa types can be defined analogously [25]. The N2HDM contains different phases, depending on the type of symmetry breaking that occurs. Vacuum expectation values for the scalar fields will lead to vacua which may, or may not, preserve the symmetries imposed. Let us now review the different types of vacua possible in the N2HDM. For the purpose of studying the interplay between different possible vacua, it is convenient to introduce a bilinear formalism, similar to that which has been developed for the 2HDM [8-10, 14, 15, 26-34]. This formalism has been applied to models with different scalar content, for instance the 3HDM [35,36] or the complex singlet-doublet model [37]. For the N2HDM let us define five real quantities, Further, we define the vectors X, A and the symmetric matrix B as in terms of which the potential of eq. (2.1) can be rewritten as In what follows we shall also make extensive use of the vector It can easily be shown that, at a given stationary point in which the fields acquire vacuum expectation values such that T , the value of the potential at that stationary point, V SP , is given by As explained in [25], by using the gauge freedom of the model, it is always possible to bring the most generic possible vacuum (in which, in principle, one would have nine different VEVs to consider, since the scalar doublets and singlet have a total of nine real component fields) to a simple form, to wit where all VEVs v X are, without loss of generality, real. The charge breaking VEV v cb breaks electromagnetic symmetry (giving the photon a mass) and the VEV v cp breaks CP conservation. It is easy to verify that these VEVs cannot coexist simultaneously. In other words, the minimisation of the potential implies that, if v cb = 0 then v cp = 0, and vice-versa. Different non-zero VEVs lead to different types of symmetry breaking, originating from minima which preserve, or not, distinct symmetries. The classification of all possible vacua was first made in [25], but here we adopt a different notation better suited for our analysis. There are two possible charge charge breaking vacua; two CP breaking vacua; two normal (electroweak breaking, but charge and CP conserving) vacua; and a single vacuum for which electroweak symmetry is unbroken. 1 Thus a total of seven possible types of vacua, or phases, exists in the model. The two electroweak breaking but charge and CP conserving vacua of the N2HDM most closely resemble a SM-like vacuum, in that they have a CP-even scalar field which can mimic the SM Higgs boson. However, the N2HDM involves extra scalars, including a charged one and several neutral ones with definite CP quantum numbers, and possibly a dark matter candidate.
The first normal stationary point N (denoted I in [25]) occurs when the parameters of the potential are such that the minimisation conditions of the potential allow a solution for which both doublets have neutral, real VEVs and the singlet has none. This vacuum therefore preserves the Z 2 symmetry of the singlet -the singlet has no VEV and does not mix with the remaining neutral scalars. Hence this corresponds to the dark matter phase of the model, with the VEVs This results in the following values for the X and V vectors (defined in eqs. (2.4) and (2.6)): with v 2 = v 2 1 + v 2 2 . The entries of V N are dictated by the N minimisation conditions and by the respective eigenvalues of the scalar mass matrices, where m 2 H ± is the squared charged scalar mass at this stationary point and m 2 D the squared mass of the singlet field. These are given by Using eq. (2.7) the value of the potential at this stationary point may be written as The second normal stationary point N s (denoted sI in [25]) corresponds to a solution of the minimisation conditions where both the doublets and the singlet Φ S acquire non-zero VEVs. This additionally breaks the singlet Z 2 symmetry -thus the singlet Φ S will mix with the remaining neutral scalars. Starting from the following VEV configuration we define where v 2 = v 1 2 + v 2 2 and m 2 H ± N s is the squared charged scalar mass at the N s stationary point, given by As before, the specific entries of V N s are a consequence of the minimisation conditions, and the eigenvalues of the scalar mass matrices at an N s stationary point. As for the value of the potential, we have As mentioned earlier, another charge and CP conserving vacuum may arise in the model -one for which the singlet field acquires a VEV but the doublets do not. This type of vacuum -dubbed S in [25] -would lead to massless electroweak gauge bosons and fermions, and as such it is unphysical. This stationary point exists if m 2 S < 0, and the singlet VEV is found to be The value of the potential at this stationary point is equal to Both the N and N s phases can accommodate SM-like physics (provided that v 2 1 + v 2 2 ∼ (246 GeV) 2 and v 1 2 + v 2 2 ∼ (246 GeV) 2 , respectively), although each of these phases has a different phenomenology (for N dark matter candidates exist, for N s three CP-even states mix with each other). We will now analyse the stability of both N and N s against the possible existence of deeper minima of different nature. For a large part of the parameter space of the model the minimisation conditions yield a single minimum, and its stability is ensured (at least at tree level). However, for many combinations of the parameters of the potential, multiple minima can coexist. If the tunnelling time from a minimum of type N (or N s) to a deeper minimum is smaller than the age of the universe then the corresponding set of parameters should be excluded.

Stability of normal minima against charge breaking
Since charge breaking minima have to be avoided, it is important to know under what circumstances a normal minimum is safe against eventual tunnelling to a deeper charge breaking minimum. In the 2HDM that question was answered [8,9,14,15] in a conclusive manner: whenever a normal minimum exists, any charge breaking stationary point is necessarily a saddle point lying above the normal minimum. In the N2HDM, as we will now show, the situation is changed. Let us first define both of the possible charge breaking stationary points and introduce some notation concerning them.
• In the first charge breaking stationary point CB (denoted IIb in [25]) the singlet field has no VEV, and the doublet VEVs are Consider also the vectors X and V evaluated at a CB stationary point, given by )/2 is one of the squared scalar masses at the CB stationary point. The entries of V CB are dictated by the CB minimisation conditions.
• In the second charge breaking stationary point CBs (denoted sIIb in [25]) the singlet also acquires a VEV, the VEV configuration being given by Analogously to what we have done for the previous stationary points, we define the following vectors: And again, the entries of V CBs are dictated by the CBs minimisation conditions.
The manipulation of the X and V vectors will allow us to establish analytical formulae relating the value of the potential at two coexisting stationary points. This technique was first used in ref. [14], and it essentially consists in following four basic steps: (1) perform the internal product of X evaluated at one of the stationary points with V evaluated at the second one; (2) repeat, with X evaluated at the second stationary point and V at the first one; (3) use the explicit formulas for V to relate the previous internal products with the value of the potential at each stationary point; (4) the two internal products will have a common term, through which they can be related to one another, thus obtaining a relation between the potentials. The technique is best understood going through some explicit examples of its application, which we will now provide. Note that all of the following conclusions are derived at the tree-level and may be affected by higher order corrections. Let us assume that the parameters of the N2HDM are such that the potential has two stationary points 2 , one of type N and another of type CB. These extrema may or may not be minima, at this time we do not need to specify it. Let us then consider the vectors defined above, containing information about the VEVs and the minimisation conditions in each extremum, for N in eqs. (2.10), for CB in eqs. (2.20).
The internal product of the vectors X CB and V N yields which may also be written as From eq. (2.7), we know that the quantity X T CB A is twice the value of the potential at the extremum CB, and therefore, combining eqs. (2.23) and (2.24), We now perform similar operations on the vectors X N and V CB , yielding The quantity X T N A is twice the value of the potential at the extremum N , hence Since the matrix B (defined in eq. (2.4)) is symmetric, the left-hand sides of eqs. (2.26) and (2.28) are identical. It is then trivial to obtain the following expression comparing the depth of the potential at both extrema, Therefore, if N is a minimum one will have m 2 H ± > 0, and since the terms in square brackets above are surely positive, one will have V CB − V N > 0. Thus we may conclude that: If the potential has a minimum of type N , any CB stationary point, if it exists, lies above N .
As such, no tunnelling to a deeper CB minimum can occur.
Similar conclusions are reached when one compares N and CBs stationary points. Again, the starting point is to analyse the internal products of the vectors X and V for each stationary point. Using eqs. (2.22) and (2.10), we obtain and therefore, subtracting both equations one easily obtains If N is a minimum all of the squared scalar masses computed therein must be positive, and thus V CBs − V N > 0.
If the potential has a minimum of type N , any CBs stationary point, if it exists, lies above N .
As such, no tunnelling to a deeper CBs minimum can occur. We can therefore conclude that minima of type N are completely stable against the possibility of charge breaking.

Extrema N s vs. CB and CBs
The analysis of the previous section can now be extended to the stability of N s minima -but as we will shortly see, the conclusions are different. Let us begin by comparing N s and CBs stationary points. As before, and using eqs. (2.14) and (2.22), we have where we use the subscript "N s " to emphasise that both the squared charged mass and the sum of the square of the VEVs concern the N s stationary point. From these equations, it is trivial to obtain Therefore, as before, if N s is a minimum, any CBs stationary point, if it exists, lies above it, and N s is stable against tunnelling to CBs. However, when one follows these steps whilst comparing N s and CB stationary points, one finds: where, recall, s is the singlet VEV at vacuum N s and m 2 S1 one of the squared scalar masses at CB. From this one obtains There is now no mandatory relationship between the depths of these stationary points -a priori, both of them can be minima, and none is privileged with respect to the other. As such -and numerical analyses prove this -there are situations in which a minimum N s coexists with a deeper CB minimum (or vice-versa). Thus we conclude: Minima of type N s are stable against charge breaking for vacua of type CBs, but not necessarily for those of type CB.
The addition of a real singlet to the 2HDM qualitatively changes the vacuum stability behaviour of the scalar potential. Whereas in the 2HDM a normal minimum is guaranteed to be stable against any possible deeper charge breaking minimum, this is no longer the case in the N2HDM. The addition of the singlet field leads to possible instabilities, where a normal minimum which breaks the Z 2 symmetry of the singlet might coexist with a charge breaking minimum (deeper or not) which does not break that same symmetry. However, any N2HDM normal minimum which preserves the Z 2 symmetry of the singlet is perfectly stable against charge breaking.

Stability of normal minima against CP breaking
Just like in the 2HDM, spontaneous CP breaking is possible in the N2HDM. In fact, CP in the N2HDM can be broken by two different minima, whereas in the 2HDM only one such vacuum can occur. Before discussing the stability of such vacua, however, some general considerations are in order: it only makes sense to study spontaneous CP breaking in models were CP is a well-defined symmetry, i.e. models that are invariant under a given CP symmetry, as is the case with the potential written in eq. (2.1) 3 . Also, care must be taken when discussing CP breaking, as it is not sufficient to have a complex valued VEV to be able to affirm that CP violation is occurring. In fact, there are situations for which CP may be preserved even if complex VEVs arise, and therefore one ought to look for other signs of CP violation, such as the couplings of scalar mass eigenstates to Z bosons. In the N2HDM, however, with the field basis we chose, no such problems arise: if the vacuum state contains a complex VEV, CP breaking occurs and produces scalar states of indefinite CP properties. Finally, as was shown in ref. [25], in the N2HDM it is not possible to have coexisting CP breaking and charge breaking stationary points -if the minimisation conditions can be solved for one type (CP breaking or charge breaking) of vacua, then the other type (charge breaking or CP breaking) admits no solution. Thus the possibility of tunnelling between CP breaking and charge breaking minima is excluded a priori.
As before, the question under which conditions a given normal minimum is stable against tunnelling to a deeper CP breaking vacuum has been previously answered in the 2HDM [8,9,14,15], and the conclusion is analogous to the charge breaking case: whenever a normal minimum exists, any CP breaking stationary point is necessarily a saddle point lying above the normal minimum. In the N2HDM the situation of the CP breaking vacua will differ, as it did for the charge breaking case. The vacua where CP can be spontaneously broken are: • The first CP-breaking stationary point CP (denoted IIa in [25]) preserves the Z 2 symmetry of the singlet but one of the doublets has a complex VEV. We parametrise the VEVs as Let us define (2.38) The entries of V CP are determined by the stationarity conditions and the form of the mass matrices at this vacuum.
• The second CP-breaking stationary point CPs (denoted sIIa in [25]) also breaks the Z 2 symmetry of the singlet and gives a complex VEV to one of the doublets. The VEVs are therefore and we define The entries of V CPs are determined by the stationarity conditions.

Extrema N vs. CP and CPs
Following the strategy employed for comparing normal and charge breaking vacua, we now assume that the potential has two stationary points, one of type N and the other of type CP, each of which may, or may not, be a minimum. With the vector definitions outlined above, we see that the internal product of the vectors X CP and V N yields and thus Eq. (2.7) tell us that X T CP A = 2 V CP and therefore With similar manipulations on X N and V CP , we obtain Since X T N A = 2 V N , and with the definition ofB in eq. (2.38), it is seen that Now, since the pseudoscalar squared mass for an N stationary point is given by we finally obtain Thus, if N is a minimum one will have m 2 A > 0, and therefore inevitably V CP − V N > 0. Following analogous steps for the N and CPs stationary points, one arrives easily at the following formula comparing the depths of the potential at each stationary point, Therefore, one reaches the same conclusions for CP and CPs stationary points, when they coexist with N : If N is a minimum, it is deeper than any CP or CPs stationary points. The conclusions of the previous subsection do not extend unchanged to coexisting N s and CP or CPs stationary points. Starting with CPs, we have Also, we derive that and after similar calculations as before it may be seen that and therefore if N s is a minimum it is certainly deeper than CPs-the same type of result we obtained when comparing N minima and CP ones. On the other hand, if we compare N s and CP stationary points, we obtain Also, it is easy to obtain and hence, after trivial manipulations, This expression shows -as for the pair N s, CB -that N s is not necessarily stable against tunnelling to a deeper CP minimum.
Minima of type N s are stable against CP-breaking minima of type CPs, but not against those of type CP.

Other coexisting neutral minima
Another possibility for vacuum instability is the existence of multiple minima of types N , N s or even S. If for instance two N and N s stationary points coexist, we can follow similar steps to those outlined in the previous sections and arrive at the following formula relating the depths of the potential: Therefore, we see that since either one of N or N s can be minima, none of them is guaranteed to be deeper than the other. Therefore, though N is stable against tunnelling to a deeper charge breaking or CP breaking minimum, it is not guaranteed to be stable against a deeper N s vacuum. Likewise, an N s minimum, which is safe against tunnelling to possible charge breaking or CP breaking minima, may be unstable against a deeper N minimum. Nonetheless, we can derive another conclusion considering this formula in tandem with the results of previous sections: If the parameters of the potential are such that N and N s minima coexist in the potential, then the global minimum of the potential preserves charge and CP.
The demonstration is simple: though from eq. (2.55) we cannot be certain whether N or N s is the global minimum, the existence of an N minimum places it certainly below any charge/CP breaking stationary points that might exist. Therefore, the conclusion becomes that either N or N s is the global minimum.
The other possibility still in play would be the coexistence of an N (or N s) minimum with an S minimum of the type described in eqs. (2.17) and (2.18), where only the singlet acquires a VEV. This is the simplest possibility of vacuum instability to verify: provided we find a solution of the N type, it will be safe against tunnelling to an S minimum provided we verify the following three conditions: • Since the S vacuum only exists if m 2 S < 0, we need not worry about tunnelling from N to S if m 2 S > 0.
• If however m 2 S < 0, then the N vacuum is deeper than S if V N is smaller than V S , with V S having a very simple form given by eq. (2.18).
• If m 2 S < 0 and V N > V S then the tunnelling time between both vacua must be computed.
Likewise for an N s vacuum, the analogous conditions for stability of N s would hold. Finally, there is still another possibility for instability of vacua of types N (or N s): that the minimisation conditions of the N2HDM may yield more than one solution for a given type of vacuum. This means that a solution of the type GeV 2 , as well as another, N ≡ {w 1 , w 2 , 0}/ √ 2 exists, with w 2 1 + w 2 2 = 246 2 GeV 2 . This possibility already arises in the 2HDM [8-12] -therein dubbed "panic vacua" -and it remains in the N2HDM as an avenue for instability of the N vacuum (and also of the N s one, since the minimisation equations of the potential may well yield more than one solution of type N s). We do not study this possibility analytically, but it is included in the numerical analysis presented in section 3.
We end this section with a very interesting scenario for the limit m 2 12 = 0, when all symmetries are exact. The N and N s stationary points are related by eq. (2.55). This equation can re-written as It we set m 2 12 = 0, and N is a minimum it is a global minimum because not only V N s − V N > 0, but also because we proved before that it is stable with respect to other charge breaking or CP breaking minima. However, this conclusion is only valid provided both doublet VEVs are non-zero, that is, the only dark matter candidate has origin in the singlet.

Vacuum stability
The results of the previous sections show that, unlike what happened for the 2HDM, when normal minima occur in the N2HDM they are not necessarily the global minima of the model. We summarise the results we obtained in table I, where we illustrate the relation between the various types of possible minima. If a minimum of type N exists (i.e. a minimum where the singlet has no VEV and its discrete symmetry is preserved even after spontaneous symmetry breaking) then N is certainly deeper than any charge or CP breaking stationary points that the potential might have -the stability of N against CP or charge breaking is perfectly guaranteed in the model. In fact, it is even possible to demonstrate (see appendix A) that in this situation any charge breaking stationary points are necessarily saddle points: an N minimum implies that at least one, but not all, of the squared masses of a CB(s) stationary point is negative. Presumably the same applies to CP(s) stationary points as well, assuming the 2HDM analysis generalizes. Of course, for considerations of stability, the nature (minimum, maximum, saddle point) of extrema that lie above N is of no consequence.
The stability found for N minima does not hold, however, for minima of type N s: for these -the discrete symmetry of the singlet is spontaneously broken in addition to EW symmetry -coexistence with minima of certain types is indeed possible. An N s minimum will certainly be deeper than any stationary points of types CBs or CPs -which break, respectively, charge conservation and CP symmetry, and also break the discrete symmetry of the singlet. But  it is possible to have coexisting N s and CB or CP minima -which break, respectively, charge conservation and CP symmetry, but do not break the discrete symmetry of the singlet. These results underline the curiously unique nature of the vacuum structure in the 2HDM, where the existence of a minimum of a given nature automatically implies that no minima of different types may exist. That property is not shared by models with a different scalar content -even in models with a simpler scalar content, such as the doublet + singlet (real or complex) model, the vacuum structure is much more complex, and no general, 2HDM-like conclusions may be drawn [37]. In models with more than two doublets the 2HDM stability also breaks down, at least concerning charge breaking [38]. What the analysis above has also shown is that the mere addition of just a real singlet to the 2HDM is enough to qualitatively change the vacuum structure of the model. The N2HDM preserves some of the nice vacuum properties of the 2HDM -wherein the N minimum mimics the stability behaviour of the normal minima of the 2HDM -but when N s minima are considered, the possibility of tunnelling to deeper minima of different types arises.

NUMERICAL ANALYSIS
In order to illustrate the impact of the N2HDM vacuum structure on the phenomenologically relevant regions of the parameter space we perform a numerical study. We study combinations of parameters that are allowed by all available theoretical and experimental constraints and analyse their vacuum structure. We first outline our method for scanning the parameter space and present the constraints we apply. In order to judge whether deeper minima are indeed excluded it is necessary to calculate the tunnelling time from the EW vacuum. We use the method developed in [22] to numerically study the vacuum structure of these parameter points and estimate the lifetime of their EW vacua.

Parameter Scan
We performed a scan of the N2HDM parameter space using an improved private version of ScannerS [25,[39][40][41]. We generated parameter points where the EW vacuum is of type N s since -following the analytical analysis -this is the most interesting case for vacuum stability. All of the resulting parameter points fulfil the applied theoretical constraints and are compatible with the applied current experimental constraints at the 2σ level.
The included theoretical constraints are tree-level perturbative unitarity [25] as well as boundedness from below [42]. A global minimum of the scalar potential only exists at finite field values if eq. (2.1) is bounded from below. This is a prerequisite for any study of vacuum stability. The allowed region is given by with Ω 1 = λ 1,2,6 > 0; λ 1 λ 6 + λ 7 > 0; λ 2 λ 6 + λ 8 > 0; λ 1 λ 2 + λ 3 + D > 0;  and and depends on the discriminant In contrast to earlier works [25,43] we do not impose absolute stability of the EW vacuum as a theoretical constraint since we want to study the vacuum structure in detail and take into account that metastable regions of the parameter space are allowed. The experimental constraints include bounds from flavour physics in the m H ± -tan β plane [44] -the B d → µµ constraint being the strongest in type I. We also require compatibility with the oblique parameters S, T and U [45,46] including the full correlation between these quantities [44]. We check for agreement with the collider Higgs data using HiggsBounds (v5.3.2beta) [47][48][49][50][51] and HiggsSignals (v2.2.3beta) [49,[52][53][54]. With HiggsBounds we check for 2σ compatibility with all searches for additional scalars, and with HiggsSignals we employ a cut on ∆χ 2 = χ 2 N2HDM − χ 2 SM < 6.18 (corresponding approximately to a 2σ region). This cut ensures that the N2HDM predictions yield a χ 2 in the fit to the LHC Higgs data that is at most 2σ worse than the one of the SM. The required model predictions for branching ratios and total widths are obtained from N2HDECAY [25,55] and the hadron collider production cross sections from SusHi [56,57].
We use this setup to generate a sample of valid parameter points on which to study the vacuum structure and vacuum stability. One of the CP-even, neutral Higgs masses is fixed to m Hx = m h125 = 125.09 GeV . (3.5) The remaining input parameters are independently drawn from uniform distributions with the ranges given in table II. The three mixing angles in the CP-even scalar sector are scanned through their whole allowed range. In this work we only consider the N2HDM of type I, i.e. where all fermions couple to Φ 2 , just mentioning briefly the results for a type II model as the vacuum structure and vacuum stability behaviour is unaffected by the choice of Yukawa type. Note that we do not specify a mass ordering for m Hx,y,z -the h 125 can be the lightest or heaviest state as well as the one in between.

Numerical Vacuum Stability
We use the approach presented in [22] to numerically study the vacuum structure and vacuum stability of the obtained parameter points. This approach is a highly efficient and numerically reliable method to study vacuum stability at the tree level in BSM models with extended scalar sectors. We will now give a short review of our approach and refer to [22] for more details.
Our code uses polynomial homotopy continuation (PHC) (see e.g. [58] or [59]) to find all stationary points of the scalar potential eq. (2.1). This method reliably finds all solutions of a system of polynomial equations -in our case given by for the real component fields ϕ i of the doublets and singlet. The value of the scalar potential eq. (2.1) at each of these stationary points is then compared to the depth of the EW vacuum. If there is no stationary point deeper than the EW vacuum we consider the EW vacuum at this parameter point as absolutely stable. If stationary points deeper than the EW vacuum exist we calculate the tunnelling time to each of these deeper extrema. The decay width per (space-)volume V S to tunnel to a deeper point in field space is [60,61] We approximate the tunnelling path by a straight path connecting the two minima in field space and use the semianalytic solution given in [62] along this path to obtain the bounce action B. The prefactor K is a subdominant contribution requiring an involved calculation and is therefore estimated on dimensional grounds. We consider the vacuum of the potential for a given parameter point to be short lived and the corresponding deeper minimum dangerous if This is a conservative estimate where only vacua with a survival probability through the age of the universe are considered short lived.

Discussion
In this section we present a numerical and phenomenological analysis of the N2HDM vacuum structure and vacuum stability. The analysis is based on the sample of 10 6 phenomenologically viable parameter points generated according to section 3.1. We aim to investigate whether the possible coexistence of minima discussed analytically in section 2 • is found in a substantial region of the N2HDM parameter space that is compatible with current theoretical and experimental constraints, • can be directly related to phenomenological observations at colliders.
Since we assume the EW vacuum to be of type N s the potentially dangerous minima are CB, CP, N , and a second different minimum of type N s (see below for a discussion of minima of type S). Unless otherwise stated, in the following we will distinguish three possibilities for these potentially dangerous vacua: • they coexist with the EW vacuum (shown in green in the following plots), • they are also deeper than the EW vacuum (shown in blue in the following plots), • they are additionally dangerous, i.e. tunnelling from the EW vacuum is fast (as defined in eq. (3.8)) (shown in red in the following plots). Table III shows the prevalence of these cases for the different possible secondary minima in our sample. While the precise numbers in table III have no physical significance as they depend on the applied method for sampling the parameter space, the displayed results clearly show that the possibilities discussed in section 2 remain relevant even after all other applicable constraints are considered. Especially, dangerous minima of type N (the dark matter phase with wrong EW symmetry breaking pattern) occur frequently in our sample. Table III also shows that the requirement of absolute stability would correspond to a substantially stronger constraint on the parameter space compared to the requirement that the EW vacuum should be sufficiently long-lived. As a consequence, important parts of the parameter space that are actually viable would be discarded if the requirement of absolute stability was imposed.   The only case missing in table III that is allowed by the analytical analysis are secondary minima of type S. However, we have not found a single parameter point in our sample where a stationary point of type S is a minimum. This could mean that minima of type S cannot coexist with an N s vacuum, that all points where this is possible are ruled out by current constraints, or that these minima are exceedingly rare. Either way, since secondary minima of type S do not occur in our sample they are of limited phenomenological interest, and we will not discuss them further here. Figure 1, left, shows the distribution of charge and CP breaking secondary minima in the plane of the pseudoscalar Higgs mass m A and charged Higgs mass m H ± . The overall distribution of the phenomenologically viable parameter points is primarily driven by the EW precision measurements which force the neutral Higgs bosons to be relatively close in mass to the charged Higgs boson. Note that parameter points without any secondary minima as well as parameter points with secondary N minima exist throughout the allowed region. In contrast, secondary CB minima only exist as long as m A > m H ± while CP minima only exist when m H ± > m A .
The origin of this strict separation -making m H ± = m A the boundary between regions where only one of these types of minima exists -can be understood analytically. The pseudoscalar and charged masses in an N minimum are such that (see eq. (2.46)) of the scalar mass matrix are given by Thus a CB minimum will imply λ 4 −λ 5 > 0 and therefore, according to (3.10), m A > m H ± . Similarly, a CP minimum requires λ 4 − λ 5 < 0 which then implies m A < m H ± . The same behaviour can be seen in fig. 1, right, showing the plane of λ 4 and λ 5 . The N minima are again scattered throughout the allowed parameter space while the CP and CB minima can only occur in sharply defined regions. Therefore, λ 4 = λ 5 would be the expected border between the regions where CP and CB can exist. However, fig. 1, right, shows that there is an additional region where neither CP nor CB minima can exist (see appendix A for an explanation).
In fig. 2 we compare the analytical result for the relative depth of N and N s vacua to the numerical results. The relative depth of an N s and N vacuum, as given by eq. (2.55), is shown as a function of tan β at the N s EW vacuum. The plot only includes parameter points where a secondary N minimum exists and shows its depth relative to the depth of the N s EW vacuum. As expected, in all parameter points where V N s − V N > 0 the N minimum is classified as either deep (blue points) or dangerous (red points). The parameter points with dangerous N only begin to appear if V N s − V N 10 7 , and their distribution shows some dependence on tan β. For small tan β 2 the N vacuum is only unstable if the depth difference is 10 9 while for large tan β 12 the majority of deep N vacua in our sample is dangerous. 4 So far we have illustrated how the analytical results of section 2 are reflected in the phenomenologically viable parameter space. We will now discuss the vacuum stability constraints arising from these secondary vacua. In imposing vacuum stability constraints we distinguish the following cases: • parameter points where the EW vacuum is the only vacuum, • absolutely stable parameter points where secondary minima exist but are never deep, • long-lived parameter points where secondary vacua are deep but never dangerous, • short-lived parameter points that have dangerous secondary minima.
as a function of the charged Higgs mass. The short-lived (different shades of red) parameter points are plotted below the grey points, for which no secondary minima exist. This means that any region where only the red parameter points are visible is excluded by vacuum stability. One can see that significant parts of the parameter space corresponding to an enhanced signal strength, µ γγ > 1, are excluded because they have a dangerous N , CP or CB minimum below the EW vacuum. If for instance a charged Higgs is found with a mass of 500 GeV, a bound of about µ γγ 1.03 in the N2HDM of type I can be derived from fig. 3. If on the other hand the charged Higgs mass could be constrained to be larger than 250 GeV (e.g. by a 500 GeV e + e − -collider) enhancements of µ γγ above 1.1 would be excluded in the N2HDM of type I by the vacuum stability constraint. One can also see from fig. 3 that if the constraint of an absolutely stable EW vacuum were imposed, the blue points in fig. 3, which indicate a long-lived EW vacuum, would be excluded, implying possibly misleading conclusions.
The reason for the behaviour observed in fig. 3, i.e. the impact of vacuum stability on the allowed µ γγ values, is the h 125 coupling to a pair of charged Higgs bosons (defined in the appendix of [25]) as shown in fig. 4  stability on the h 125 couplings to gauge bosons and to fermions. Therefore, µ γγ is the observable where the vacuum stability constraint is expected to have the largest impact since, among the currently measured observables, it has the highest sensitivity to the possible effects of a triple scalar coupling. The large impact of the vacuum stability constraint on µ γγ is specific to the N2HDM of type I. This is due to the fact that in type I all Yukawa couplings are rescaled by the same factor c(h 125 ff ). The cancellation of this factor which occurs in µ γγ in the approximation Γ tot (h 125 ) ≈ Γ(h 125 → bb), leads to an increased sensitivity to Γ(h 125 → γγ) and thus to g h125H + H − . In contrast, for Yukawa types where c(h 125 tt) = c(h 125 bb) (e.g. type II) the effect of vacuum stability constraints on µ γγ is no longer visible as the ratio of Yukawa couplings has a much stronger impact on the signal rate than the charged Higgs contribution to Γ(h 125 → γγ).
It is interesting to note that although the allowed range for µ γγ is very similar in the type I 2HDM [25] and in the type I N2HDM, a measurement of µ γγ above 1 for certain charged Higgs masses could exclude the N2HDM but be compatible with the 2HDM due to the different vacuum stability constraints. Figure 5 shows vacuum stability constraints in the plane of the mass m H2 of the second lightest Higgs boson H 2 , with a mass above H 1 = h 125 , and the signal strength µ τ τ of h 125 (defined analogously to eq. (3.13)). In this case, there are hardly any regions where points can be clearly excluded due to the existence of a secondary dangerous vacuum. There are regions where only points with a non-stable vacuum exist, which can be either dangerous or long-lived, but no direct bounds can be derived from the experimental measurements of µ τ τ . This is due to the fact that in contrast to fig. 3 these regions are always populated by long-lived metastable vacua, so that allowed parameter points exist in these regions. Therefore, fig. 5 clearly shows the phenomenological difference between requiring an absolutely stable EW vacuum (keeping only the grey and green parameter points) and a long-lived EW vacuum (additionally keeping the blue parameter points). As discussed above, enforcing absolute stability could lead to misleading phenomenological conclusions. The signal strength µττ of h125 → τ τ as a function of the second lightest neutral scalar mass mH 2 . The parameter points without any secondary minima (grey) are plotted on top, followed by the absolutely stable (green), and long-lived (blue) parameter points. Below these, the points with dangerous secondary minima are shown in different shades of red denoting the type of dangerous minimum present (N -light red, CB -red, CP -dark red).

CONCLUSIONS
We have performed a detailed analysis of the vacuum structure of the N2HDM, an extension of the SM by an extra doublet and an extra real singlet. We have shown that it is possible to derive analytical expressions to compare minima of different nature. In the case where the singlet has no VEV the conclusions are the same as for the 2HDM [14], that is, minima of different nature, N , CB and CP, never coexist. We have also shown analytically that when the singlet acquires a VEV, if a normal N s minimum exists, it is stable against tunnelling to a corresponding charge breaking CBs or CP-breaking CPs extremum. However, that conclusion no longer holds when comparing minima with and without singlet VEV. In fact, minima of different natures can coexist and potentially tunnel into each other. Moreover, it is known that in the 2HDM minima of type N are not unique [8][9][10][11][12] and the existence of a second, normal minimum (panic vacuum) can exist below the one with the correct EW symmetry breaking. In the N2HDM panic vacua of types N and N s can appear for EW vacua of either type. Additionally, minima of type S with only a singlet VEV could also appear as panic vacua. However, we have not found a single parameter point in our sample where a stationary point of type S is a minimum.
Based on this analytical analysis we have conducted a numerical study to investigate the impact of the intricate N2HDM vacuum structure on the phenomenology of the model. We have generated a large sample of parameter points with an EW vacuum of type N s that fulfil all applicable theoretical and experimental constraints (without enforcing that the EW vacuum be a global minimum). This way, we were able to compare minima of different nature and identify regions of parameter space where the EW vacuum is the global minimum, where deeper minima exist but tunnelling is so slow that the EW vacuum is long-lived, and regions that are excluded because the tunnelling time is short compared to the age of the universe.
The first important conclusion of our study was that panic vacua of type N , as well as charge breaking CB, and CP breaking CP minima deeper than the EW vacuum appear in a significant portion of the (otherwise) phenomenologically viable parameter space. We have also shown the distribution of secondary CB and CP minima and established the boundaries of the disjunct parameter regions where these minima can exist.
Studying the impact of vacuum stability on collider observables we have found that a precise measurement of µ γγ above 1 could exclude the model on the grounds of vacuum stability alone, unless the charged Higgs is very light. This is due to the sensitivity of µ γγ to the triple Higgs coupling g h125H + H − , which is constrained by vacuum stability. If the Yukawa sector is of type I this effect is clearly visible in µ γγ because of an approximate cancellation between the modifications of the Yukawa couplings. In the study of other collider observables, such as µ τ τ , we showed that there are large regions where minima which are absolutely stable do not occur, but a long-lived EW vacuum exists. This illustrates the importance of including parameter regions with a metastable vacuum in phenomenological analyses, as enforcing absolute stability may lead to incorrect conclusions.
where we introduced the matrix B and the vector V which are defined, respectively, in eqs. (2.4) and (2.6). Then, since for a CBs stationary point V l = 0, the mass matrix [M 2 ] is reduced to the second term in the equation above. It is then rather easy to reproduce the calculation in section 5.2 of ref. [14] and deduce that one may simplify the expression of [M 2 ] and obtain where C is a 5 × 5 matrix depending only on the VEVs. Eq. (A.4) demonstrates that the eigenvalues of [M 2 ] at a CBs stationary point will be all positive if and only if the matrix B is positive definite. However, we have shown above that when N is a minimum, the matrix B has at least one negative eigenvalue -and therefore [M 2 ] has also at least one negative eigenvalue. However, since B also has positive eigenvalues, so will [M 2 ]. Therefore, if N is a minimum then any CBs stationary point, it if exists, lies above N and is a saddle point, q.e.d.
We can now also justify the conditions of eq. (3.12) for the non-existence of neither CB or CP minima. The matrix B determines the nature of the CBs stationary point, and one can also show that it does the same for the CB extrema. Checking now eqs. (2.4), we see that the (3, 3) entry of B is λ 4 + λ 5 , and therefore, if λ 4 < −λ 5 one of the diagonal elements of B will be negative -thus B cannot be positive definite, and consequently no CB minima can occur (only saddle points). This justifies the second condition of eq. (3.12). As for the first one -λ 5 < 0 -the nature of CP stationary points will, in analogy with the CB cases (and the 2HDM, see [14]) be determined by a matrix of the quartic couplings. For the CP extrema, however, that matrix is not B but rather the matrixB, eq. (2.38). Observe then that the (3, 3) element ofB is 2λ 5 -and therefore, if λ 5 < 0 no CP minima can occur sinceB cannot be positive definite.