ScannerS: Constraining the phase diagram of a complex scalar singlet at the LHC

We present the first version of a new tool to scan the parameter space of generic scalar potentials, ScannerS. The main goal of ScannerS is to help distinguish between different patterns of symmetry breaking for each scalar potential. In this work we use it to investigate the possibility of excluding regions of the phase diagram of several versions of a complex singlet extension of the Standard Model, with future LHC results. We find that if another scalar is found, one can exclude a phase with a dark matter candidate in definite regions of the parameter space, while predicting whether a third scalar to be found must be lighter or heavier. The first version of the code is publicly available and contains various generic core routines for tree level vacuum stability analysis, as well as implementations of collider bounds, dark matter constraints, electroweak precision constraints and tree level unitarity.


Introduction
The recent discovery of a scalar particle [2,3] at CERN's Large Hadron Collider (LHC) has boosted the activity in probing extensions of the scalar sector of the Standard Model (SM). So far, the experimental results indicate that this scalar is compatible with the SM Higgs boson. However, many of its extensions are also compatible with the present experimental results. In fact, we know that several of them will never be completely disproved even if the Higgs has all the properties one expects from a SM Higgs. The limit where this occurs is known as the decoupling limit and is characterised by a light Higgs with SM couplings to all other known particles while the remaining scalars are very heavy. It is possible though that other scalars are waiting to be found at the LHC or that a more precise measurement of production cross sections and branching ratios of the newly found 125 GeV scalar reveals meaningful deviations from the SM predictions. As one needs to be prepared to address that scenario we have developed a new code to deal with the vacuum structure of scalar extensions of the SM. The code presented in this work, ScannerS [1], is intended to contribute to the analysis of the different scenarios that can be suggested by the LHC experiments.
In the SM, electroweak symmetry breaking (EWSB) is achieved by one complex SU (2) L scalar doublet. Although the pattern of EWSB can be correctly reproduced by this one scalar doublet, the SM is not able to accommodate a number of experimental results such as the existence of dark matter or the measured baryon asymmetry of the universe. Adding a scalar singlet to the potential could provide a viable dark matter candidate [4][5][6][7][8][9][10][11][12][13][14] as well as means of achieving electroweak baryogenesis by allowing a strong first-order phase transition during the era of EWSB [15][16][17][18][19]. Although minimal, this extension provides a rich collider phenomenology leading to some distinctive signatures that can be tested at the LHC [20][21][22][23][24][25][26][27].
In this work we focus on the conditions on the parameters that lead to the existence of a global minimum in the scalar potential. However, contrary to previous works, we impose these conditions not to single out one particular model (or phase), but rather use that information to distinguish between possible coexisting phases. This way, we expect to identify properties of the models by classifying the possible phases as function of the parameter space point. We have also included in the code the most relevant theoretical and experimental constraints both from dark matter and collider searches. Our goal is therefore to investigate the possibility of using measured experimental quantities (e.g. the mass and branching ratios of a new light scalar at the LHC), to automatically exclude one of the possible phases (e.g. a phase with a dark matter candidate) using the phase diagram of the model. To that purpose we scan over the entire parameter space subject to the most relevant constraints and plot the results in projections that include physical quantities whenever possible.
As the ultimate goal of ScannerS is to be a tool that can be used for general scalar sectors, the core routines of the code can already be used for larger extensions of the scalar sector in terms of its field content. The core includes generic local minimum generation routines (with goldstone/flat direction identification as well as a-priori curved directions/symmetries) and tree level unitarity check routines. However, since the first version of the code has been tested extensively with the complex singlet models we present here, routines for testing electroweak observables can be used only for extensions with n-singlets, and global minimum and boundedness from below routines are defined on a model by model basis by the user. The same naturally applies for some of the analysis of experimental bounds (since such analysis will depend on the model that is being analysed) though some data tables from several experiments are generic and included in the core. We expect in the near future to automatise global minimum and boundedness from below routines to be in the core of the program. Details on the structure of the program and on how to use it can be found in [1].
The structure of the paper is the following. In Sect. 2 we describe the scalar potential of the models we will study. We derive the conditions to be fulfilled for the minimum to be global at tree level for those specific models, based on symmetries, and classify the various phases for each model. In Sect. (3) we address the problem of efficiently performing a scan of the parameters space in a more general perspective. We propose a method to generate a local minimum by generating vacuum expectation values (VEVs), mixings and masses uniformly, and obtaining some of the (dependent) couplings from the linear systems of equations which characterise the minimum. In Sect. 3.2 we also discuss the generic implementation of tree level unitarity bounds. In Sect. 4 we implement constraints from various experimental sources, such as electroweak precision observables, LEP bounds, dark matter bounds, and constraints from Higgs searches at the LHC. Finally, we conclude with a discussion of our main results in Sect. 5.

The models
We consider an extension of the scalar sector of the SM that consists on the addition of a complex singlet field to the SM field content. Our starting point is the most general renormalisable model, invariant under a global U (1) symmetry. We then consider two models with explicit breaking of this symmetry, but where Z 2 -type symmetries are preserved for one or two of the components of the complex scalar singlet. Each model is classified according to phases associated with the spontaneous symmetry breaking pattern of the vacuum. Considering the allowed parameter space, each model can exhibit more than one global minimum with the correct pattern of EWSB for fixed couplings. In these models, the mixing matrix is usually lower dimensional, allowing for one or two unstable scalar bosons that could be detected at the LHC and at future colliders, together with two or one dark matter candidates, respectively. However, if all Z 2 symmetries are broken we can end up with three unstable scalars (and no dark matter candidate) complicating the signatures and making a clear identification of the Higgs boson a much more difficult task. Our aim is to investigate the parameter space of these extensions as to identify the properties of the models which are not yet excluded by theoretical and experimental constraints. Furthermore, we want to investigate the possibility of, given a set of experimental measurements related to the scalar sector at the LHC, that we can automatically use the phase diagram of the model to exclude one of the phases, e.g. a phase with dark matter or a phase where all the scalars are mixed.
The scalar field content of the model is as follows. There is the SM Higgs doublet H which is a singlet under SU (3) C and is in the fundamental representations of SU (2) L × U (1) Y , i.e. (1, 2, 1/2). We add a complex scalar field S = S +iA which is a singlet under the SM gauge group. This is equivalent to adding two real singlet fields. Several models have been studied in the literature by imposing special symmetries on this scalar sector. Our starting point are the models discussed in [25], where the Higgs potential has a global U (1) symmetry in S. Besides spontaneous breaking of this symmetry by the vacuum, we can break it explicitly at various levels, by soft linear and quadratic terms which are technically natural and do not generate other soft-breaking terms through renormalisation [25]. The scalar potential is then (soft breaking terms in parenthesis) We must ensure that the potential is bounded from below and in this simple model an analytic condition can be found easily. Noting that the quartic operators of V are quadratic forms in the real positive quantities H † H and |S| 2 , positive definiteness of such forms is equivalent to boundedness from below and implies Before specifying the models, we note that the stationarity conditions are: which have no closed form solution in general. We will expand the Higgs doublet and the complex scalar field around the vacuum according to where the Higgs Vacuum Expectation Value (VEV) is v = 246 GeV. As previously stated, simpler versions of this model were already discussed in the literature. In most cases a specific model is singled out by imposing additional symmetries on the model. Because one of the main motivations for adding a scalar singlet to the SM is to provide a dark matter candidate, models with one or two dark matter candidates are then analysed and confronted with experimental results. Here, we are interested in applying our new tool ScannerS to various versions of such models, while applying the latest experimental bounds, together with theoretical constraints. The points that pass all constraints during the scan are then classified according to the phase they are in (see classification below) and plotted, whenever possible, in a physical projection of the parameter space, such as a mass of a new particle or a measurable rotation angle between group and mass eigenstates. A physical measurement allows us in some cases to discriminate a phase with dark matter candidates from one with no dark matter candidates. Consequently, the phase which is realised for a given model, can in such cases be decided by experiment.
Theoretically, the independent models and their phases are classified according to their symmetry group and spontaneous breaking, respectively, as follows: • Model 0, U (1) symmetry with up to two dark matter candidates. This is obtained by imposing a U (1) symmetry on the complex singlet field, which eliminates the soft breaking terms, thus a 1 = b 1 = 0. There are two possible phases: 1. S = 0 at the global minimum (symmetric phase). Then we have two degenerate dark matter particles. In this scenario the model is equivalent to two independent real singlets of the same mass and quantum numbers.
2. S = 0 at the global minimum (broken phase). The U (1) symmetry is spontaneously broken, then there is an extra scalar state mixing with the Higgs and a (massless) Goldstone boson associated with the breaking of the symmetry. Since the phase of the complex singlet is unobservable, without loss of generality we can take S = 0 and A = 0 for such phase and A is the dark Goldstone particle. This phase is however strongly disfavoured by observations of the Bullet Cluster [7,25,28,29]. These observations can be used to constrain the mass of the dark matter particle as a function of the value of δ 2 . Hence, unless δ 2 is vanishingly small, a zero mass dark matter particle is ruled out.
As one of the phases is not allowed, we will discard this model from our discussion.
• Model 1, Z 2 × Z 2 symmetry with up to two dark matter candidates. This model is obtained by imposing a separate Z 2 symmetry for each of the real components of the complex singlet. The Z 2 symmetries imply that the soft breaking couplings are a 1 = 0 and b 1 ∈ R (6 real couplings in the scalar potential & no other couplings are generated through renormalisation). Specialising the minimum conditions (2.3) we obtain the following qualitatively different possibilities for minima with v 2 = 0: 1. S = 0, no mixing and two dark matter candidates (symmetric phase).
2. S = 0 or A = 0, one of the singlet components mixed with the Higgs doublet and one dark matter candidate (spontaneously broken phase). One can show, by noting that swapping S ↔ A only changes the sign of b 1 , that without loss of generality we can take A = 0, while still covering the full parameter space (this is so because b 1 ∈ R, and the potential only depends on squares of the VEVs). This is true in our scans only because we will adopt the strategy of first generating a locally viable minimum and, only after, to check all possibilities for minima below the one generated.
• Model 2, One Z 2 symmetry with up to one dark matter candidates. This is obtained by imposing a Z 2 symmetry on the imaginary component A. Then the soft breaking couplings must be both real, i.e. a 1 ∈ R and b 1 ∈ R. Looking at the minimum conditions we find the following cases: mixing between h (SM Higgs doublet fluctuation) and S only (symmetric phase). In this case we can take S ∈ R + as long as a 1 runs through positive and negative values.
To summarise this classification, we show in table 1 a list of the three possible models (labelled by their symmetry group), with a description of the particle content of the two possible phases (symmetric or broken) as well as the VEV/symmetry breaking pattern. Model Higgs+2 degenerate dark 2 mixed + 1 Goldstone Higgs + 2 dark 2 mixed + 1 dark 2 mixed + 1 dark 3 mixed VEVs at global minimum Table 1. Phase classification for the three possible models.
Once we have picked one of the cases above, we have to check for all possible minima by evaluating the potential at all possible stationary points. A list of the possible cases is provided in appendix A.
We have chosen to start the studies with our new tool with these models, because they already have a great diversity of physically different phases (which provide some interesting results -see Sect. 5), while allowing us to test the routines of the code we propose to develop. We have reproduced several results presented in the literature with a very good agreement. In the next sections, we describe the ScannerS scanning strategy and the main tree level theoretical constraints to generate a stable vacuum with the correct symmetry breaking pattern.

The scanning method
To implement the scan of the parameter space of the models, we have developed a dedicated tool ScannerS which can be used for more generic potentials. Our method is based on a strategy of reducing as many steps as possible to linear algebra, since these are computationally less expensive. Before describing the method in detail, let us describe the general idea.
A possible way of performing the scan for a generic scalar potential (strategy 1), and determining the spectrum of scalars is: i) scan the couplings λ a uniformly in chosen ranges, ii) determine all stationary points by solving the (non-linear) stationarity conditions, iii) check if any are minima and choose the global one iv) if yes, accept the point, compute the mass matrix, diagonalise and check if the masses of the scalars (m 2 i ) and mixing matrix (M ij ) are consistent with the symmetry breaking imposed. This method contains two steps which are quite expensive, computationally, which are executed before we know if the minimum has the desired properties. The first is the determination of the stationary points, which in general is a problem of finding the solutions of a polynomial system of equations in the VEVs. The second, not as expensive, is the diagonalisation of the mass matrix at the global minimum.
Our alternative strategy (strategy 2) relies on two observations. First, a generic scalar potential V (φ i ) (where the φ i are the real fields used to decompose the scalar multi-plets/singlets into real components), is a linear form in the couplings λ a where (for renormalisable models) V a (φ i ) are monomials in the fields, of degree up to four. Second, in strategy 1, the couplings λ a are taken as independent parameters, and the VEVs, masses and mixing parameters are determined. However, any independent set of parameters is equally good to label a certain point in parameter space, so if instead we use VEVs, masses and mixing matrix elements as the parameters to be scanned over, then the determination of dependent couplings λ a becomes a linear algebra problem due to the linearity property of (3.1). Furthermore, this set of parameters (VEVs, masses and mixings) is more directly related to physical properties of the scalar states, so it is a more natural choice of parametrisation.

Generation of a local minimum
The details of the method we use are as follows. We first want to express the Lagrangian in terms of physical propagating scalar degrees of freedom H i with masses m i , which are the fluctuations around a minimum with VEVs v i . The most general expansion of the original fields such that the kinetic terms remain canonical is then where M ij is a generic rotation matrix, and we have assumed that the original kinetic terms for the φ i fields are canonically normalised. The (linear) stationarity conditions (vacuum conditions) are where we have defined the derivative of the a-th operator with respect to the field φ i evaluated at the VEVs by ∂ i V a . Such object is as matrix with indices {i, a} acting on a vector of couplings λ a . These conditions are independent of the mixing matrix. If we choose values for the VEVs (by scanning over an interval or fixing some of them such as the SM Higgs VEV), this system of equations can be analysed as to identify a sub-set of couplings to be eliminated in favour of the remaining couplings (let them be grouped in a sub-vector λ a 1 where the sub-index a 1 runs only over dependent couplings). Otherwise, if the system is over-determined we must reject such choice of VEVs. The dependent couplings λ a 1 can be eliminated in favour of the remaining ones, whose subset we denote λ a 2 1 . This elimination is represented in a matrix form as follows where the matrix Λ a 1 a 2 can be found numerically through Gaussian elimination using the system (3.3). This procedure amounts to replacing the λ a 1 couplings by the VEVs that have been scanned over. The vacuum conditions can be implemented in any VEV of some operator O acting on the potential, using the matrix Λ a 1 a 2 : where in the last step we have defined the "VEV reduced" action of O on the potential, OV a 2 , which contracts (linearly) with the sub-vector of couplings λ a 2 . The second step is to impose that the stationary point is a minimum, which is done by assuming non-negative masses squared. Thus we write the quadratic derivative conditions which involve the masses and mixing matrix elements, i.e.
where we have used Eqs. (3.2) and (3.5), and in the last line we have replaced latin indices by a bold face matrix notation. The hatted indexî denotes no summation over i. For a given model, the mass matrix before diagonalisation may have eigen-directions which are independent of the values of the couplings λ a 2 . In such case we say that there are "a priori" eigen-directions at the minimum. This can be tested by finding the eigen-vectors of a matrix ∂ 2 V a 2 for fixed a 2 , and then check which eigen-vectors remaining eigen-vectors of the other matrices with different a 2 . We may then find "a priori" flat eigen-directions (which would be Goldstone bosons) or curved eigen-directions (which would be massive particles) 2 . Once those directions are identified, the form of the mixing matrix is restricted to a block diagonal form. The flat directions fix the mass of the corresponding particles to zero regardless of λ a 2 , so effectively they eliminate a set of conditions from system (3.6). For each curved eigen-directions, there is a condition with the mass squared of the particle on the right hand side. For the block where mixing occurs, since the mixing matrix is symmetric, we have n(n + 1)/2 conditions (n is the dimension of the block). Along the diagonal of the block, n conditions contain a mass squared on the right hand side. The n(n − 1)/2 conditions off the diagonal conditions contain a zero on the right hand side. At this stage the undetermined parameters are: the λ a 2 couplings; the mixing matrix elements M ij ; and the physical masses m 2 i . Because we want to avoid solving non-linear equations, we choose to generate first the mixing matrix uniformly 3 . Then we are left with a set of conditions relating λ a 2 with the masses m 2 i , which can be re-arranged as a homogeneous system. If we define the vector of parameters which are still undetermined by v T = (λ a 2 , m 2 i ) such homogeneous system is where the matrix D is read out from (3.6). Using again Gaussian elimination we can solve for a sub-set of parameters of v, as a function of another subset which we are left to scan over (if the system is over-determined we must reject the vacuum).
With this procedure, we end up generating a point in parameter space, and all the properties of the physical states are determined (such as masses and mixing matrices), having avoided the problem of solving a system of non-linear equations. The price to pay is that we have delayed checking if the minimum is global. However, the advantage of this procedure is that we do not spend any time in points of parameter space where there is no stationary point with the correct properties. Furthermore, we can add to this procedure more constraints on the local minimum (if they are computationally quick to check), before checking if it is global. If we use strategy 1, the first steps involve a computationally intensive non-linear problem, which will be wasted each time a point in parameter space is rejected, whereas with strategy 2 we generate a local minimum with the desired properties quickly and only then do we have to perform the computationally intensive steps.

Tree level unitarity
An important constraint on the region of parameter space to be scanned over is given by imposing tree level unitarity at high energies. This was pioneered by Lee, Quigg and Thacker in [31]. It was shown that the Goldstone high energy theorem applies and all we need is to compute all the scalar quartic interaction amplitudes, describing 2 → 2 processes between any (suitably normalised) scalar two particle states. The basis of such states for N real scalars is  where we have emphasised the 1 √ 2! term which arises from the normalisation of harmonic oscillator (indistinguishable) 2-particle states, and the second block contains all the pairings of different field states. Then, one constructs an s-wave amplitude matrix where T (0) is the J = 0 contribution to the usual transition matrix. In the second step we have used the Feynman rules for the scalar sector to express each matrix element as a sum over quartic vertices labelled by the index a 4 . If the |Φ i state contains identical particles, n i = 2 otherwise it is one. P a 4 is defined as the product of the factorials of the powers of each field in the corresponding monomial V a 4 (φ). For a generic scalar potential, (3.10) yields real coefficients. Tree level s-wave unitarity, is then implemented by requiring that the absolute value of each eigenvalue of a ij is smaller than 1/2. This matrix can be computed efficiently at each point of parameter space, and diagonalised, as to check if unitarity is preserved. For the complex singlet model we are considering, one can write down exactly the most restrictive conditions which are which correctly reduces to the SM tree level unitarity bound when d 2 = δ 2 = 0.

Electroweak precision observables
Another source of constraints comes from electroweak precision observables. Here we focus on the S, T, U variables [32,33]. For an extended sector with scalar singlets coupling only to the Higgs doublet, the only extra contributions compared with the SM are in the self energies Π ZZ (q 2 ) and Π W W (q 2 ) for the Z and W particles respectively. Using the general expressions in [24], one can show that the variations with respect to the SM contributions from a single Higgs doublet are where the terms inside the brackets are to be evaluated with the values of the mixing matrix elements M hj that represent the mixing between the SM Higgs and the new singlet field components. The masses m 2 i correspond to the new scalar physical eigenstates. The functions f (Y ) and g(Y ) are defined as where the covariance matrix is defined in terms of the correlation matrix, ρ ij , and the standard deviation, σ i , of each parameter through σ 2 ij ≡ σ i ρ ij σ j . To test these observables, we have used results for the SM global fit of the Gfitter collaboration [34]

LEP and LHC bounds
An important set of experimental constraints comes from collider searches for the Standard Model Higgs boson. The standard strategy is to compute, for each search channel, the predicted signal strength defined as where σ New (H i ) and σ SM (h SM ) are the production cross sections of H i and the SM Higgs respectively, both evaluated at the mass of H i ; Br New (H i → X SM ) is the H i branching ratio (BR) to SM particles and Br SM (h SM → X SM ) is the SM Higgs BR to SM particles evaluated at the mass of H i . In the models we are considering, the scalars couple to the SM particles always through the same combination h = M T hi H i = M ih H i . Therefore, both the production cross sections and the decay widths are just rescaled by the factor M 2 ih . We can then write µ i as . (4.10) However, because there are new particles involved, the ratio of BRs is now where the term Γ(H i → newscalars) is only present when the channels for which the SM Higgs decays to the new scalars are open, and once again Γ(h SM → X SM ) denotes the SM Higgs width evaluated at the mass of H i . We only consider two-body final states for the Higgs boson decaying to other scalars. Then, when kinematically allowed, the decay widths for a process of the type H i → H j H j and H i → H j H k are given respectively by where g ijj , g ijk are the coupling strengths between the corresponding scalars i, j, k and m j is the mass of the scalar state H j . The combined LEP data [35] from the four LEP collaborations sets a 95 % confidence level upper bound on the H i ZZ coupling in non-standard models relative to the same coupling in the SM, through the quantity where g BSM H i ZZ designates the non-standard H i ZZ coupling while g SM HZZ stands for the corresponding SM coupling. In the last line we have related this quantity to the signal strength µ i for this particular channel, as defined in Eq. (4.10). We apply the LEP limits for each scalar mass eigenstate H i in the bb and τ + τ decay channels separately.
With the recent discovery of a Higgs like boson by the ATLAS and CMS collaborations, we now have a very strong constraint on BSM models. In any extension of the SM one of the scalars is bound to have a mass of approximately 125 GeV. At the same time both LHC experiments have also constrained any new scalar couplings to the SM particles, that is, they have provided us an exclusion region for µ i as function of the mass of the scalar H i . These bounds are also included in the program and the details are as follows. We assume a Higgs boson with mass m h = 125 GeV, and allow for a signal strength µ i in the interval 1.1 ± 0.4 [36]. In the mass regions where a scalar particle is excluded, we apply the 95 % CL combined ATLAS upper limits on µ i as a function of m i [36] for all other non 125 GeV scalars, as long as their production is allowed at the LHC.

Dark matter experimental bounds
We have applied the limits on the total dark matter relic density from WMAP (Wilkinson Microwave Anisotropy Probe) 7-year measurements of cosmic microwave background (CMB) anisotropies [37] Ω cdm h 2 = 0.112 ± 0.006 (4.15) where h is the Hubble constant (h = 0.704 ± 0.025 in units of 100 km/Mpc/s). We have calculated numerically the relic density for the DM candidate with the software package micrOMEGAS [38,39]. We have allowed the contribution from the dark matter candidate in each model to be below the limit for the relic density, taking into account that there could be another DM contributor.
Another important constraint on these models is the elastic cross section from DM scattering off nuclear targets. The most recent direct detection DM experiments have placed limits in the spin-independent (SI) scattering cross section (σ SI ) of weakly interacting massive particles (WIMPS) on nucleons. The most restrictive upper bound on the SI elastic scattering of a DM particle is the one from XENON100 [40]. In ScannerS we have included the limits of σ SI as a function of the dark matter mass as presented in [40].
In these singlet models, the scattering cross section of the dark matter candidate with a proton target is given by [41] where m p is the proton mass, v = 246 GeV is the SM-Higgs VEV and f pi are the proton matrix elements with central values [42]: The dominant contribution for SI scattering comes from t-channel scalar exchange. We note that the cross-sections for scattering off protons is very similar to the one for neutrons. Therefore we present only the results for the scattering off proton targets. Finally we should mention that we have used micrOMEGAS [38] to perform an independent calculation of the cross sections and that we have found a very good agreement with our result. We have written the cross section expressions for the models considered in Appendix B.

Discussion
In this section we analyse the results of full scans over all parameter space, for the various phases of each model. We focus on models 1 and 2, which we denote by cSM 1 and cSM 2 respectively. We have performed two main scans for each model. One with a smaller hypercubic box in parameter space to allow for a better resolution of the region being scanned over, and another wider scan to check which boundaries did not change significantly. Unless stated otherwise, we always use the wider scan 5 . We will also present a scan with some of the constraints removed to clarify the appearance of certain boundaries in the allowed regions of parameter space. We start by presenting some results for model 1 (cSM 1) where we label the phase with two dark matter candidates as "2DM" and the phase with only one dark matter candidate by "1DM". The key of all figures contains an extra label for each colour of the points, indicating whether all the bounds/constraints discussed in previous sections have been included, and when not, the corresponding constraints that have been removed, are indicated.
In Fig. 1 we show two particular projections in parameter space which show regions which are exclusive of a given physical phase. On the left panel we display m heavy as a function of m D,light where m D,light is the mass of the lightest dark matter particle and m heavy is the mass of the the other non-SM scalar; on the right we show m heavy − m D,light as a function of m D,light . Note that the heavy scalar is in fact a dark matter candidate in the phase "2DM". There are regions exclusive to the phase with two dark matter candidates, and also very small regions where just one dark matter state is possible. A measurement of the mass of a dark matter candidate could give us a hint on the phase the model is in. However, even if it is possible in some cases to distinguish between phases with actual physical measurements, it is rather hard to discriminate between phases in model 1, even more so because in some regions it requires independent measurements of the properties of a dark matter candidate and an unstable scalar mixing with the Higgs.
Model 2 (cSM 2 ) on the other hand displays more interesting possibilities. Here we label the phase with one dark matter candidate and a new scalar by "1DM", whereas the phase with no dark matter candidates and all three scalars mixed, by "mix". Fig. 2 (left) shows the phase diagram projection in the plane (λ, |M ih |), where λ is the doublet quartic coupling and |M ih | is the mixing matrix element with the SM Higgs doublet component of (any of) the new mixed scalar(s) which is not the 125 GeV scalar. Three sets of points are displayed in this plot (where all bounds have been taken into account): the phase "mix" with a wide scan and the phase "1DM" both with a standard and a wide scan. The wide and standard scans allow us to discriminate fixed boundaries from moving boundaries between the two phases. It is clear that the left most boundary is moving towards the vertical axis and therefore cannot be considered a physical boundary between the two phases. This is On the left we plot λ (the doublet quartic coupling) as a function of |M ih | which is the mixing matrix element with the SM Higgs doublet component of (any of) the new mixed scalar(s) which is not the 125 GeV scalar. In this case, we show a standard and a wide scan for the "1DM" phase. On the right we display the effects of applying the different bounds to the "1DM" phase obtained for a wide scan. For such case there is only one extra mixed scalar so |M 2h | is the mixing matrix element of such non-dark state with the SM Higgs doublet component. The points have been overlaid following the order in the key (first in the key list is the bottom layer in the plot).
observed by comparing the boundary between the navy and blue points. Concentrating on the two wide scans (pink and navy colour) it is clearly possible to separate the two phases for some regions of the parameter space; for instance, for small values of |M ih |, λ has to be close and above its SM value in phase "1DM" while in the "mix" phase it can be either above or below that value in a much wider interval. Thus, in this projection, one can clearly identify two pink regions which are exclusive of the "mix" phase (excluding the region close to the moving boundary along the vertical axis).
In the right panel of Fig. 2 we display the effects of applying the different bounds to the "1DM" phase obtained for a wide scan. In this case the horizontal axis has the variable |M 2h | which is the mixing matrix element with the SM Higgs doublet component of the non-dark matter state. The points have been overlaid following the order in the key -first in the key list is the bottom layer in the plot. It can be seen that EWPO constraints cut out a considerable portion of the parameters space and are responsible for the upper right boundary between the allowed region for the one dark matter phase as compared to the mixed phase in the left panel. Nevertheless, the theoretical constraints alone are responsible for the bottom boundary as seen from comparing the left and the right plots. Note that the apparent top left boundary is not meaningful, and should shrink to the vertical axis if we perform increasingly wider scans in parameter space, as discussed in the previous paragraph.
In Fig. 3 we plot the projection on the (λ, m light ) plane, where m light is the mass of the lightest non-dark matter state. If some other scalar is detected at the LHC there are values of λ that are only allowed in the "mix" phase. The SM value for λ (fixed by the Higgs and W-boson masses) is the horizontal line slightly above 0.5. If the new scalar is lighter than the SM Higgs the "1DM" phase only exists below that value while if it heavier than the SM Higgs that phase is only present for values above the SM λ .
More interesting are the cases presented in Fig. 4 that show scatter plots of various projections of the parameter space points obtained for a wide scan of the two phases of model 2. In the horizontal axis we have either m light or m heavy which are the masses of the lightest or heaviest of the new mixed scalars which are neither the 125 GeV Higgs nor the dark matter candidate. In the vertical axis we have either |M ih | (mixing matrix element with the SM Higgs doublet component of the non-dark matter scalar corresponding to the mass in the horizontal axis), or |M 1h | which is the 125 GeV scalar mixing matrix element. Note that if the 125 GeV scalar is not allowed to decay to any of the other two scalars, |M 1h | is just √ µ, that is, the square root of the signal strength parameter. As before, the points have been overlaid following the order in the key (first in the key list is the bottom layer in the plot). In this case we are dealing with directly measurable quantities only. Finding a new particle and measuring its decay rates will gives us access to both the masses and |M ih | while there are already results for |M 1h | from the LHC measurements. In fact note that current LHC bound of the signal strength already impose a strong constraint on the mixing yielding |M 1h | 0.84 [36]. We should note however that this is just the 1σ bound for √ µ -had we taken the 2σ result the limit on |M 1h | would be relaxed.  bounds also force |M ih | 0.55. The plots in Fig. 4 clearly show that a combination of measurements of a new scalar at the LHC may decide for a given phase or exclude the scenario altogether. The top plots show that a scalar with a mass close the SM Higgs mass is allowed to exist in the "1DM" phase. Furthermore, there are regions which are clearly exclusive of a no dark matter phase. The same trend appears in the bottom plots even if the distinction between phases is not so striking. A very interesting property, that is observed in both bottom and top plots, is that given a measurement of |M 1h | or |M ih | and the mass of the new scalar in a region which is exclusive of the "mix" phase, one can infer whether a heavier or lighter scalar is expected to be observed (within this model). For concreteness consider for example the pink (light grey) regions in the left plots of the figure, there is clearly a big portion which does not exist in the corresponding region of the right plots, so if a measurement falls in that region, we can immediately say that if we are in the cSM 2 we are looking at the light scalar of the "mix" phase. Similarly, there is a pink region in the right plots which do not exist in the left plot, so in such case one can say we are observing the heavy scalar in the "mix" phase of the cSM 2.

Conclusions
We have presented a new tool, ScannerS, devoted to the search for global minima in multi-Higgs models. In this work we have applied it to some versions of a simple extension of the SM -the addition of a complex singlet to the SM doublet, with some symmetries. The code includes the most relevant theoretical and experimental bounds. Our main focus is in distinguishing the possible phases of each model by using the present experimental data both from the LHC and from dark matter experiments.
Once we have identified our working models based on symmetries, we have excluded all phases that did not display the correct electroweak symmetry breaking pattern. We ended up with two possible phases for each model. In the first model, which we called model 0, one of the phases leads to a massless dark matter candidate which is already excluded by the Bullet Cluster results. The other phase has two dark matter candidates. Because in the allowed phase there is no mixing between scalars, the only way to tell them apart is by actually detecting a dark matter particle. Hence, a study with ScannerS would add no advantageous information to what we know experimentally. The cases of model 1 and (especially) of model 2 are the most interesting. We have shown that by measuring physical quantities like the particle masses, mixing angles, or quartic couplings we are able, in some particular cases, to pinpoint the phase that is realised in Nature if one of these models applies. Most importantly, in model 2 a simultaneous measurement of the mass of a non-dark matter scalar together with its mixing angle could be enough to exclude a dark matter phase, and simultaneously indicate whether we are observing the lightest or the heaviest of the new scalar states expected in the model. Nevertheless, as we move closer and closer to the SM limit the phases become more and more indistinguishable.
In summary, although the differences that we found between phases of these singlet extensions are restricted only to some measurable quantities, we found it possible to fall into regions where measurements will definitely exclude one of the phases and predict properties of the scalar spectrum. An interesting question is whether such differences between phases can be identified for more complicated scalar extensions or if they can even become more impressive and predictive.
The ultimate goal of ScannerS is to provide the community with a tool that can be used to search for global minima in general scalar sectors. Although the core routines can already be used for an arbitrary scalar sector, the present release [1] requires the user to define the boundedness from below and global minimum routines. In the next release, we expect to provide core routines for such tasks, as well as further analysis examples such as the important case of the two-Higgs doublet model. The present publicly available code includes examples with all the analysis used in this article, which we hope will be a useful starting point for users who intend to explore this tool.

C Ranges for the scans in parameter space
In this appendix we present the ranges that were used for the scans over the couplings λ a , masses and VEVs of the scalars. We have performed two different scans, which are indicated in the tables respectively by "standard run" and "wide run". In the latter we allow the hypercubic box in parameter space to be wider to investigate which boundaries in the phase diagram do not changed significantly. Each table contains the two possible phases of each model indicated appropriately. Table 2, Left, shows the ranges of couplings, masses and VEVs of the scalars in the symmetric phase "2DM" of model cSM 1 (two dark matter phase) and for the two different scans (standard and wide). h refers to the Standard Model Higgs and D 1 and D 2 the two dark matter candidates of the model. Table 2, Right, is for the same model (cSM 1 ) but in the broken phase where there is only one dark matter particle (1DM); h is the Standard Model Higgs which mixes with the scalar H 1 and D the dark matter candidate.
The last table is for model 2 (cSM 2 ). Table 3, Left, is for the symmetric phase, which contains a dark matter particle (1DM). Similarly to model 1, the scalar spectrum is given by the Standard Model Higgs, h, which mixes with the scalar H 1 and D is the DM candidate. For the broken phase, Table 3, Right, the scalar spectrum contains the Standard model Higgs, h, which mixes with the two scalars H 1 and H 2 .  Table 3. Parameter ranges for cSM 2 . Left: "1DM" symmetric phase, Right: "mix" broken phase.