A stable Higgs portal with vector dark matter

We explore an extension of the Standard Model by an additional U(1) gauge group and a complex scalar Higgs portal. As the scalar is charged under this gauge factor this simple model supplies a vector dark matter candidate satisfying the observed relic abundance and limits from direct dark matter searches. An additional Higgs-like state, that may be heavier or lighter than the observed Higgs, is present and satisfies LEP and LHC bounds whilst allowing for absolute stability of the electroweak vacuum in a range of the parameter space.


Introduction
The discovery of a 125 GeV scalar [1,2] completes the Standard Model (SM) of particle physics and appears to confirm some rather intriguing features of nature. Firstly, the Higgs potential of the SM develops an instability at high field values, [3][4][5][6][7][8][9][10]. Secondly, the theory on its own does not account for nearly 95% of the matter in the universe, known as dark matter (DM) [11] and dark energy.
In this paper we explore a simple extension of the SM by an extra U(1) X gauge symmetry factor, supplemented with an additional complex scalar charged solely under this U(1) X . It was studied by [12][13][14][15] and by [16][17][18] in its non-Abelian version (for other realizations of vector fields in the context of dark matter see [19,20] the real component of the complex scalar obtains a vev) generate a viable massive vector dark matter candidate. It also provides an additional Higgs mass eigenstate and introduces additional freedom in the theory to completely alleviate the issue of vacuum stability.
In this paper we use a collection of HEP-tools and software to give an accurate account of the currently viable parameter space of such a model. In particular we combine SARAH [21][22][23][24], Pyrate [25] and SPheno [26,27] to explore stability at 2-loops and micrOMEGAs [28][29][30][31] for a detailed study of vector dark matter for this model. We also compare our results to various experiments including LUX [32] and XENON100 [33].
The paper is organized as follows. In section 2 the model is defined and a detailed exposition of the parameters of the model is presented. In section 3 the renormalisation group equations of the theory are used to constrain the parameter space from physical consideration such as positivity, stability and the absence of Landau poles in the theory. In section 4 collider constraints are applied and the dark matter abundance constraint is imposed on the parameter space that remained. A detailed study of the DM-nucleon cross section is also performed. Finally in section 5 the findings are summarized.
In appendix A we give a detailed account of the tree-level calculation of the relic abundance. Appendix B supplies the SARAH model files which fixes the conventions and allows for the generation of code for RGES and for the code used in micrOMEGAs [28][29][30][31].

A model of vector boson dark matter
The vector dark matter (VDM) model [12][13][14][15][16] is an extension of the SM by an additional U(1) X gauge symmetry together with a complex scalar field S, whose vev generates a mass for this U(1)'s vector field. The quantum numbers of the scalar field are S = (0, 1, 1, 1) under U(1) Y × SU(2) L × SU(3) c × U(1) X . (2.1) None of the SM fields are charged under the extra gauge group. In order to ensure stability of the new vector boson a Z 2 symmetry is assumed to forbid U(1)-kinetic mixing between U(1) X and U(1) Y . The extra gauge boson A µ and the scalar S field transform under Z 2 as follows A µ X → −A µ X , S → S * , where S = φe iσ , so φ → φ, σ → −σ. cos θ W = g g 2 + g 2 , sin θ W = g g 2 + g 2 . (2.4)

JHEP09(2015)162
At leading order the vector bosons masses are given by: gv, M Z = 1 2 g 2 + g 2 v and M Z = g x v x , (2.5) where v and v x are H and S vacuum expectation values (vev's): ( H , S ) = 1 √ 2 (v, v x ). The scalar potential for this theory is given by It will also be useful to define, for future reference, the parameter λ SM ≡ M 2 h /(2v 2 ) = 0.13, where M h = 125.7 GeV.
The requirement of positivity for the potential implies the following constraints that we impose in all further discussions: (2.7) Hereafter the above conditions will be referred to as the positivity or stability conditions. It is easy to find the minimization conditions for scalar fields (without losing generality one can assume v, v x > 0): For the VDM model only the latter case is relevant, since both vevs need to be non-zero to give rise to the masses of the Standard Model fields and dark vector boson. Both scalar fields can be expanded around corresponding vev's as follows The mass squared matrix M 2 for the fluctuations (φ H , φ S ) and their eigenvalues read (2.11) The matrix M 2 could be diagonalized by the orthogonal rotation R, such that M 2 diag = R −1 M 2 R. The convention adopted for the ordering of the eigenvalues and for mixing angle α is the following Note that since vev of H, if fixed at 246.22 GeV, with κ = 0 (no mass mixing) and λ H = λ SM it is only h 2 which can have the observed Higgs mass of 125.7 GeV. Even though the mass matrix is diagonal in this case, however in order to satisfy our convention that M h 1 = 125.7 GeV a rotation by α = ±π/2 is required in such a case. There are 5 real parameters in the potential: µ H , µ S , λ H , λ S and κ. Adopting the minimization conditions (2.8) µ H , µ S could be replaced by v and v x . The SM vev will be fixed at v = 246.22 GeV. Using the condition M h 1 = 125.7 GeV, v 2 x could be eliminated via (2.11) Therefore eventually there are 4 independent unknown parameters in the model x ) together with (2.14) one finds the following universal formula for the mass of the non-standard Higgs: It is easy to see that positivity of M 2 h 2 is guaranteed if the following conditions are satisfied: • for λ H > λ SM , λ H > κ 2 4λ S + λ SM (same as (2.15)).
From (2.11) one can easily derive the following useful inequalities Therefore for λ H > λ SM we have M + > M h 1 , so the lighter scalar is the SM-like Higgs particle, while when λ H < λ SM the observed Higgs particle is the heavier one. Consequently λ H is the sole parameter that distinguishes between these two scenarios. To illustrate that behavior we show in figure 1 contours of non-standard Higgs masses in the plane (κ, λ H ). Note the presence of M 2 h 2 = 0 contour (the most external blue contour in the lower part of the figure) that corresponds to vanishing determinant of the matrix (2.11): In fact it is interesting to consider the special case of κ 2 − 4λ H λ S = 0, JHEP09(2015)162 then if one requires that the potential has a minima for v, v x = 0 it is necessary to assume that µ 2 H /λ 1/2 S , which implies that the potential could be written as H . This potential has equipotential contours (in the unitary gauge) on ellipses such that λ 1/2 H |H| 2 + λ 1/2 S |S| 2 = constant. Fluctuation along the ellipses that corresponds to the minimum is the massless mode. Note that this is a mode that exists even if v x = 0, so this parameter choice is different from the case discussed below at the end of this section where v x is approaching 0 (so det(M 2 ) → 0 as well).
The behavior of v x is presented in figure 2. For λ H < λ SM and λ H > κ 2 4λ S (as required by the scalar mass positivity) one finds Similarly for λ H > λ SM and  and it is diverging at the parabola λ H = κ 2 /(4λ S ) + λ SM . For λ H < λ SM only the physical region corresponding to M 2 In order to understand the behaviour of the mixing angle α one can directly adopt the formulae (2.13). Since α varies in the range [−π/2, π/2], the absolute value of |α| and its sign can be read from the inverse of cos(2α) and sin(2α), respectively. The coupling of h 1 eigenstate (the one that has 125.7 GeV mass) to V V is proportional to cos α therefore the LHC data favours regions of α ∼ 0. As can be found from (2.13) small α corresponds to either to v x → ∞ (for λ H > λ SM ) or to v x → 0 (for λ H < λ SM ). One can see the same behaviour from figure 2.
It is worth investigating the SM limit of the VDM model. Figure 2 is a good starting point as it is easy to recognize regions in the parameter space that imply vanishing corrections (relative to the SM) to h 1 couplings (this is what we define by the SM limit). So, for λ H > λ SM this is the parabola λ H = κ 2 4λ S + λ SM where v x → ∞ and α vanishes, while for λ H < λ SM it is the vicinity of λ H = λ SM . The case λ H > λ SM is less interesting as it is just the decoupling limit with M h 2 ∼ 2λ S v 2 x + · · · → ∞, see figure 1, when the effective low energy theory is just the SM with h 2 being integrated out. The mixing angle α behaves in that region as α ∼ κ/λ S (v/v x ) 2 + · · · . More interesting is the region with λ H approaching λ SM from below, as there the mass of h 2 goes to zero (also v x → 0), as seen in figure 2. In this region the model contains a SM-like scalar with the mass of 125.7 and almost massless state (h 2 ) that decouples from V V , h 1 h 1 and fermions, however its cubic (κh 1 h 2 2 ) and quartic (κh 2 1 h 2 2 ) scalar couplings remain. This limit of the VDM model turns out to be phenomenologically unattractive since the Higgs boson h 1 would decay invisibly into pairs h 2 h 2 .

JHEP09(2015)162
3 The renormalisation group equations As will be illustrated shortly, in order to investigate vacuum stability in the VDM model it is necessary to use 2-loop RGEs. We have adopted SARAH [21][22][23][24] to obtain the full set of 2-loop beta functions. However, in order to facilitate further discussion we show below the beta functions of the gauge couplings at 1-loop level, even though whenever we impose any constraints on the parameter space that rely on the RGE running we always use 2-loop beta functions.
The various Higgs quartic couplings at 1-loop run as Above Y u,d,e denote the corresponding Yukawa matrices. Since here we are mainly concerned with the case of the masses for the extra scalar h 2 and Z of the order of the electroweak scale and since initial conditions for RGE running will be specified at Q = m t (relatively large scale) therefore we will adopt the above beta functions neglecting decoupling of extra degrees of freedom below the scale Q = M h 2 , M Z . We have verified that, in the SM limit, our 2-loop running of λ H (Q) agrees with known results [34].
In order to explore the stability and positivity of the theory we used the two-loop RGEs and the tree-level potential. To improve the precision of this work further one would likely need the one-loop improved effective potential before extending to 3-loop RGEs.
An example of a representative point in the parameter space at Q = m t is pictured in figure 3, where we show both 1-and 2-loop running of the gauge and scalar quartic couplings. A few comments are here in order. The running of gauge couplings (left panel of 3) is rather stable and similar for 1-and 2-loop beta functions. However the Yukawa coupling Y t already shows (left panel) some sensitivity to the approximation adopted for the beta functions. It is important to note the relevance of 2-loop running of λ H (Q); as seen in the right panel of figure  nearly the same for 1-and 2-loop beta functions. It is worth noticing that β λ S is always positive, so the stability condition λ S > 0 can never be violated radiatively. The evolution of κ is rather mild as β κ ∝ κ therefore at least for small κ the evolution is quite suppressed. In the right panel of figure 3 we also show the running of the positivity condition κ(Q) Here the 2-loop effects are again important. The 1-loop curve terminates already around Q ∼ 10 11 GeV, the scale at which 1-loop λ H (Q) becomes negative, so that the positivity condition κ(Q) + 2 λ H (Q)λ S (Q) > 0 can not be verified. On the other hand the 2-loop running of this condition shows that it is satisfied up to Q ∼ 10 12 GeV. The choice of initial conditions for this plot illustrates the fact that there exist initial conditions (i.e. points in the parameter space (λ H (m t ), κ(m t ), λ S (m t ), g x (m t ))) such that even though the condition κ(Q) + 2 λ H (Q)λ S (Q) > 0 is satisfied at low scale it fails at high energies.

Stability
The constraints (2.7) can be used to determine areas of parameter space (λ H (m t ), κ(m t ), λ S (m t ), g x (m t )) in which the conditions for stability/positivity of the potential are satisfied at all renormalization scales. In the SM, the absolute stability is ensured just by the positivity of the quartic coupling at all energy scales Q: λ H (Q) > 0. Since the mass of the Higgs boson is known experimentally the initial condition for running of 13 and for this initial value λ H (Q) becomes negative at some scale causing the instability. However here, in the presence of the extra scalar S this is not necessarily the case; the LHC Higgs mass measurement fixes the following combination of couplings and vev's: It is easy to see that the VDM model has the freedom to increase the value of λ H at low scales; a freedom which the SM does not possess. Larger initial values of λ H such that λ H (m t ) > λ SM are allowed delaying the instability (by shifting up the scale at which λ H (Q) < 0). There is also another remedy for the instability within the VDM model; even JHEP09(2015)162 if the initial λ H is smaller than its SM value, λ H (m t ) < λ SM , still there is a chance to lift the instability scale if appropriate initial value of the portal coupling κ(m t ) is chosen. This effect is caused by the positive κ 2 contribution to the beta function β λ H that partially compensates the negative top-quark effect, see (3.3). Figure 4 illustrates the way the SM stability problem encoded by λ H (Q) < 0 for Q > Q could be relaxed within the VDM model. The white region above the horizontal line λ H (m t ) = λ SM shows the region of λ H (m t ) > λ SM so the positivity of λ H up to the Planck scale could be easily guaranteed. On the other hand the white region below the line λ H (m t ) = λ SM shows those pairs of (λ H , κ) for which even though the starting point for λ H evolution is lower than for the SM, nevertheless the extra positive contribution to β λ H makes λ H positive up to the Planck scale. Clearly for large κ the stability region increases (for negative κ the other stability condition gives tighter constraint). The colorful regions show the scale at which λ H becomes negative. The three panels shown in figure 4 correspond to three different pairs of initial values for (g x (m t ), λ S (m t )). As seen, the sensitivity to those choices is very weak even though g x (m t ) and λ S (m t ) vary in a wide range, in fact this is understandable since the evolution of λ H (Q) is influenced by g x (m t ) and λ S (m t ) only indirectly through the presence of κ 2 in the beta function β λ H , see (3.3).
In the VDM model the SM stability problem (positivity of λ H ) is easily solved as was illustrated above. However in this case positivity requires two extra constraints: λ S > 0 and κ + 2 √ λ H λ S > 0. Since β λ S > 0 therefore whenever λ S (m t ) > 0 the positivity is preserved during the evolution. However the second extra condition is non-trivial, as illustrated in figure 3 it is possible that κ(Q) + 2 λ H (Q)λ S (Q) changes sign while running from low energies up. Figure 5 shows the scale at which κ(Q) + 2 λ H (Q)λ S (Q) becomes negative as a function of (λ(m t ), κ(m t )) for three fixed sets of (g x (m t ), λ S (m t )).
Stability of the U(1) X dark matter model was also discussed in [35]. The vacuum stability induced by the dark matter has also been considered in the context of complex [36] and real [37,38] extra scalars serving as dark matter candidates. JHEP09(2015)162 Figure 5. The "in between" stability frontier: these plots identify the scale t * = Log 10 (Q * ) at which the positivity condition κ(Q) + 2 λ H (Q)λ S (Q) fails and the vacuum becomes unstable, as a function of (λ(m t ), κ(m t )) for fixed choices of (g x (m t ), λ S (m t )) specified above each panel. The horizontal solid black line corresponds to λ H (m t ) = λ SM 0.13. The gray area is excluded by the requirement that there is no Landau poles up to the Planck mass.

The Landau poles
As we have discussed above, the additional freedom in the Higgs sector, that is due to the presence of the Higgs portal κ|S| 2 |H| 2 allows one to increase the low scale value of λ H sufficiently to avoid its negative value (instability in the H direction) at high scales. However, this possibility is bounded from above by the requirement that there are no Landau poles in the evolution of λ H (or any other parameter -a pole in the evolution of any coupling implies divergence of all of them at the same energy) up to a chosen high scale, e.g. the Planck mass. In figure 6 we show contour plots of λ H (M P l ) in (λ H (m t ), κ(m t )) space for fixed g x (m t ) and λ S (m t ). It is clear that too large κ(m t ) or λ H (m t ) implies early divergence of λ H (Q). Also when g x (m t ) and/or λ S (m t ) grow (from left to the right panel) the safe region shrinks in agreement with expectations.  Left and right panels corresponds to positive and negative κ, respectively. The allowed area decreases as the magnitude of κ increases. The blue area denotes region where h 1 coupling to vector bosons is away from its SM value by more than 15%. The region over dashed red line is excluded at 95% CL by the analysis of the Peskin-Takeuchi S,T parameters. not be too large since otherwise the expansion has no chance to converge. The maximum adopted by various authors is to some extent subjective and usually varies between 1 and 4π. Here we have the advantage of knowing both 1-and 2-loop beta functions for the RGE evolution of the couplings, therefore the relevance of the 2-loop could be quantitatively estimated for different (large) values of couplings. We found that requiring the 2-loop correction to be smaller than 100% of the 1-loop result for any quartic scalar coupling, i.e. |(λ i | < 1, implies that the coupling should not exceed a value close to 2π. Therefore, in the numerical results, whenever it is stated that a Landau pole appears, it is meant that the corresponding coupling reaches a value of 2π. For larger couplings the 2loop contributions start to dominate, so that one can not trust the perturbative expansion, truncated to this order. The VDM model without the additional U(1) X reduces to a Higgs portal extension of the Standard Model, which has been explored recently in [39][40][41]. These papers note that it can be useful to consider a change of parameterisation from (λ H , κ, λ S ) → (M h 2 , κ, sin α), in which one can re-express the quartic couplings as follows: not have Landau poles at any renormalisation scale before the Planck scale, as pictured in figure 7. We show here those plots in order to compare our results with those of [41], which differ in so far as our study include 2-loop (not 1-loop) running, and where our model includes the additional parameter g x .

Constrained parameter space
In this section we show results of scans over the parameter space (λ H , κ, λ S , g x ). Usually λ S will be fixed and specified, g x is either fixed (in figures 9-11) or varied 0.1 < g x < 1 while λ H and κ will be scanned such that 0 < λ H < 0.25 and −0.5 < κ < 0.5 (with some exceptions when −1 < κ < 1). If RGE running is employed then the scan range should be regarded as for initial values of running quantities at Q = m t .

Collider constraints
In addition to the theoretical requirements discussed above we are going to impose some experimental constraints. First of all there exist limits on branching ratio for invisible Higgs boson decays. Searches for Higgs decaying invisibly has been carried out by both the ATLAS and CMS collaborations at the LHC for various production and decay channels. ATLAS [42] considered an associated Higgs production with a vector boson (V = W ± or Z), assuming SM production they found an upper limit of 29% at 95% confidence level on the branching ratio of Higgs bosons decaying to invisible particles. A search for invisible decays of Higgs bosons was also performed by CMS [43,44]. Assuming Standard Model Higgs boson cross sections and acceptances, the observed upper limit on the invisible branching fraction was found to be 57% at 95% confidence level. It turns out that within the VDM model unless the U(1) X coupling constant g x is tiny or 2M Z > 125.7 GeV decays of the observed Higgs boson h 1 into Z Z would dominate with branching ratio exceeding the experimental limits. Similar comments apply for decays of h 1 into pairs of h 2 (for the scenario with λ H < λ SM ). Therefore in our analysis we simply exclude points in the parameter space such that h 1 decays to Z Z or h 2 h 2 are kinematically allowed. If λ H < λ SM Higgs boson h 2 is light so that it could have been prodeced at LEP. Then the LEP limits for e + e − → Zh 2 apply and should be imposed. Here we adopt the data from [45] where limits on ZZh 1 coupling normalized to the SM ZZh coupling (κ Z = cos α) are tabularised as a function of M h 2 . For fixed λ S the limit on α could be easily translated into allowed region in the (κ, λ H ) plane.
Higgs boson couplings are being measured at the LHC. The ATLAS [46] and CMS [47] collaborations limit e.g. ratios κ V of Higgs boson V V couplings normalized to their SM values. The conclusion is that the observed Higgs has SM-like couplings to the vector bosons. Here we will assume, somehow arbitrarily, that the ratio is limited by 0.85 < κ V < 1. Note that because of the orthogonal mixing κ V = cos α can not exceed 1 within the VDM model. As we have already mentioned above the SM-like nature of h 1 favours regions of large and small v x for λ H > λ SM and λ H < λ SM , respectively. The constraint that originates from 0.85 < κ V < 1 could be expressed as an allowed region in the (κ, λ H ) plane both for λ H < λ SM and λ H > λ SM case.

JHEP09(2015)162
In order to estimate relevance of electroweak precision data we adopt the Peskin-Takeuchi S, T and U parameters [48]. 2 At the 1-loop level, beyond the SM radiative corrections to the SM vector boson progators δΠ V V are not affected by the presence of the Z boson, therefore δΠ V V are the same as those found in the analysis of the plain Higgs portal [41]. At this order shifts in δΠ Zγ and δΠ γγ vanish, therefore the S and T parameters can be expressed as whereas the parameter U is too small to be relevant. Using the fit obtained in [49] we found that at 95% confidence level the S and T parameters do not constrain further the parameter space. These bounds are entirely embedded in region, where the scalar mixing angle is too large or couplings are non-perturbative (see figure 7). It is worth to emphasize that they constrain | sin α| only moderately weaker than the full set of the electroweak precision observables [41]. In particular, in the important mass range M h 2 > 200 GeV, the allowed value of | sin α| in S and T approximation is larger by 25%. For low M h 2 mass differences grow however that region of parameter space is anyway disfavoured by the requirement of absolute stability and excluded by limits on κ V .

Dark Matter abundance
Before we proceed to constrain the parameter space by measurements of DM abundance, in figure 8 we show results for Ω DM h 2 as a function of the DM mass (M Z ) obtained varying coupling constant g x and choosing a few representative values of other parameters. We have calculated Ω DM h 2 adopting standard textbook methods for cold relics, see e.g. [50]. Relevant vertices, Feynman diagrams and corresponding contributions to Z Z annihilation cross section for various final states are collected in appendix A. We have checked our results for Ω DM h 2 against calculations done adopting the micrOMEGAs3 [31]. It turned out that except resonance regions (such that 2M Z ∼ a mass of s-channel resonance) and vicinities of thresholds for new final states, Ω DM h 2 determined via the cold relics technique agrees pretty well with the result provided by the micrOMEGAs. However, in order to have also those special regions under control we have decided to adopt the micrOMEGAs hereafter. For illustration, in 8 we compare Ω DM h 2 obtained by adopting results contained in appendix A with those from micrOMEGAs.
There is a comment here in order. In figure 8 and similar that will follow, one sees that for large M Z it is typical that Ω DM h 2 decreases as M Z grows. In fact it is easy to understand such behaviour. In those plots potential parameters are fixed, so is v x , therefore increasing M Z implies growing g x , so that σv increases and therefore Ω DM h 2 decreases. This fact can be seen in the easiest way by looking at the contribution from direct DM-scalar interaction coming from the vertex V Z ij ∝ g 2 x (appendix A); then σv ∝ |V Z ij | 2 /M 2 Z ∝ g 2 x /v 2 x and therefore Ω DM h 2 ∝ 1/g 2 x .
2 Since the scale of new physics in our case is not always much above the electroweak scale therefore the S, T and U parameters should be used with some extra care. In particular, for M h 2 < M h 1 they do not provide a viable estimation of radiative corrections. Luckily that region of parameter space is not allowed by other constraints.

Allowed parameter space
All constraints that we have considered are collected in figures 9, 10 and 11. Choosing randomly points in the parameter space in the region (white) allowed by perturbativity, stability, LEP and LHC data we require that the DM abundance remains within the 5σ limit of [51] Ω DM h 2 = 0.1199 ± 0.0022 (4.2) Points that fit into the allowed region of (4.2) are shown in figures 9-11 as dots within the white region, for them all the constraints are satisfied.

Direct detection of Dark Matter
Interactions of VDM with nucleons are mediated by the Higgs particles. For the elastic scattering cross section we obtain where µ = M Z M N /(M N +M Z ) is nucleon-DM reduced mass, g hN N = 1/v N | q m qq q|N is the effective Higgs to nucleon coupling [52]. In this subsection we are going to show results for σ Z N as a function of M Z for those points (green points in figures 9-11) in the parameter space which satisfy all the constraints including the DM abundance. In the following subsection we are showing results obtained with micrOMEGAs for λ S fixed at λ S = 0.2.

Light dark matter
We start with the case of λ H < λ SM . Then, since v x is limited from above as in (2.20) and the scanning interval for g x is 0.1 < g x < 1, therefore the DM mass is bounded from above by M Z < v(λ SM /λ S ) 1/2 ∼ 200 GeV. It turns out that in practice there is no consistent points in the parameter space with M Z ≥ 120 GeV. Note also that in order to prevent JHEP09(2015)162  invisible h 1 decays h 1 → Z Z it is required that M Z > 62.9 GeV, so the allowed range for DM mass is 62.9 GeV < M Z < 120 GeV. In figures 12 and 13 we are showing σ Z N coloured with respect g x , M h 2 and sin 2 (2α) in order to learn properties of the points that are plotted. We also show experimental limits for σ Z N from XENON100, LUX (2013) and anticipated results for XENON 1T. As it is seen from the figures, points that are consistent with the present data correspond to medium gauge coupling g x ∼ 0.5, h 2 slightly lighter than the observed Higgs M h 2 ∼ 110 ÷ 125 GeV and small mixing angle 0 α π/8. The DM mass varies between 60 and 120 GeV with heavier states favoured. Figure 14 is helpful in order to understand why there is no points in the parameter space with M Z 120 GeV. As it is seen in the right panel of the figure for large M Z the DM abundance Ω DM h 2 is very small since almost all annihilation channels (except h 1 h 1 and tt) are open for Z of that mass (so that the annihilation cross section is large). It is also instructive to look at correlations between M h 2 and M Z shown in the left panel of figure 14. It turns out that there are two regions consistent with the constraints: i) M Z ∼ 65 GeV and ii) M Z ∼ M h 2 − 5 GeV. The first one corresponds to the vicinity of h 1 resonance in the right plot in the figure. The dip at 62.9 GeV is quite steep such that deviation by about 5 GeV is sufficient to reach Ω DM h 2 ∼ 0.1. The other side of the resonance is excluded by the requirement of no invisible h 1 decays. In fact it also easy to understand the second region that corresponds to the other side of the summit seen to the right of the h 1 resonance in the right panel. There the sudden drop of Ω DM h 2 is caused by the opening of h 2 h 2 final state. In a real case, since annihilating Z pairs are not exactly at rest, therefore the h 2 h 2 channel opens by 5 GeV earlier. Figure 13. Same as in figure 12, but colouring is with respect to M h2 and sin 2 (2α) for the left and the right panel, respectively.

Heavy dark matter
In this section we are considering the case λ H > λ SM . Then, since v x is limited from below as in (2.21) and the scanning interval for g x starts at 0.1 therefore the DM mass is bounded from below by M Z > 0.1v(λ SM /λ S ) 1/2 ∼ 20 GeV. However since invisible Higgs decays should be prevented therefore Z must be even heavier, so that here M Z > 62.9 GeV. Results for σ Z N are presented in a similar manner as in the case of the light DM, so in figures 15 and 16 we are showing the cross section plotted against the DM mass M Z with colouring corresponding to g x , M h 2 and sin 2 (2α). As it is seen from the figures there exist points that lay below the LUX 2013 upper limit, they correspond to medium gauge coupling g x ∼ 0.2÷0.6, wide range of M h 2 varying from 130 GeV up to 1000 GeV and small mixing angle 0 α π/8. The allowed DM mass varies between 62.9 GeV to 1000 GeV.  It is worth understanding the origin of points that compose the blue and green hills in the vicinity of 100 GeV in figure 15. As seen from the figure they correspond to different strengths of the gauge coupling constant: the left one is made of points with g x ∼ 0.2 while for the right one the coupling is typically larger g x ∼ 0.5. The purpose of right panels in figures 17 and 18 is to illustrate Ω DM h 2 dependance on M Z for representative points (λ H , κ, λ S ) in the parameter space. Note in the left panels of the figures that points in the space (M Z , M h 2 ) are grouped into two rough sets, one for M Z M h 2 /2 and the other for M Z M h 2 /2. Inspecting the right panel of figure 17 it is easy to see that the later group is composed of points sitting on the left slope of a heavy Higgs (h 2 ) resonance (at M h 2 = 400 GeV for this particular example). Since the resonance dip is relatively wide in this case, therefore the shift above the the nominal resonance is substantial, typically of the order of 100 GeV. Then when M Z increases eventually h 2 h 2 channel opens up, the cross section increases and Ω DM h 2 drops leading to higher M Z consistent with the observed DM abundance, that explains the other group of points in e.g. figure 17 located above M Z ∼ M h 2 /2. The two groups are separated by the presence of the summit in Ω DM h 2 between them. The figure 18 shows that the two groups have different g x and therefore they could be identified as the points that compose the two hills in figure 15: the blue hill is made of points that lay below M Z ∼ M h 2 /2 while the green one of those with M Z M h 2 /2. The right panel in figure 18 illustrates the mechanism of the very strong correlation between M Z and M h 2 that is observed for large M Z along the line M Z ∼ M h 2 . As it is seen from the figure the correlation is caused by the steep drop in Ω DM h 2 that corresponds to opening of h 2 h 2 final state.

JHEP09(2015)162
Constraints coming from the direct detection can be compared with results of [17]. However, a detailed comparison is quite complicated, therefore we limit ourself to a conclusion that qualitatively results obtained in [17] for the Abelian case agree with those found here. The results presented in the figure 2 therein present similar behaviour to those of figure 17 in this work. It can be seen that for a given M h 2 , when M Z approaches M h 2 /2 the nucleon scattering cross section diminishes and LUX bounds can be easily satisfied, whereas for M Z between (M h 1 +M h 2 )/2 and M h 2 the cross section is substantially larger (the points for M h 2 < 600 GeV). Similarly the LUX bounds constrain the vicinity of M Z = 80 GeV. We have compared the sensitivity of the running of the scalar quartic couplings and the stability of the vacuum on the loop order of the RGEs and found that the 2-loop RGEs make a significant effect on the 1-loop result, which could be improved further with higher order corrections. We found that an increase in λ H (m t ) by even a modest amount easily eliminates the electroweak vacuum stability problem up to the Planck scale, which is possible due to the additional degrees of freedom in this theory. When the Higgs portal κ|S| 2 |H| 2 is present (so κ = 0) then the initial value λ H (m t ) could be chosen larger than its Standard Model value (where it is fixed by the Higgs mass measurement) in this way escaping the dangerous possibility of λ H (Q) becoming negative. In addition the model discussed here can have a stable vacuum even if λ H (m t ) is chosen below its Standard Model value, since the RGE beta function for λ H contains an extra positive contribution proportional to κ 2 which can lift λ H efficiently enough while running up. For this mechanism to work κ(m t ) must be large enough.

Summary
The presence of an extra neutral scalar makes the phenomenology of the scalar sector more attractive and richer but still testable at the LHC and future colliders. We have focused here on dark matter aspects of the model. The parameter space of the models has been investigated in detail and regions consistent with theoretical, collider and cosmological constraints have been determined. In particular we have shown that the DM-nucleon cross section could be consistent with the LUX and XENON100 limits, and also looked at why this is so. It has also been shown that the anticipated XENON 1T limits on σ Z N may be satisfied within this model.

JHEP09(2015)162
Acknowledgments BG and MD acknowledge partial support by the National Science Centre (Poland) research project, decision no DEC-2014/13/B/ST2/03969. This work was partially supported by the Foundation for Polish Science International PhD Projects Programme co-financed by the EU European Regional Development Fund. This work has also been partially supported by National Science Centre (Poland) under research grant DEC-2012/04/A/ST2/00099.

A Relic abundance
In this appendix we present a detailed account of the tree-level calculation of the relic abundance in this model. This has been used to compare with the result generated with micrOMEGAs.

A.1 Thermally averaged cross-section
Thermally averaged cross-section is defined as the annihilation cross seciton σ times Møller velocity v averaged over the Boltzmann distribution (A.1) In the nonrelativistic limit m T it can be approximated using the formula T andσ(s) = 2 s(s − 4M 2 Z )σ(s) can be written as the sum of contributions coming from all possible final stateŝ Cross sectionsσ for fermions and vector bosons (presented in A.4) can be easily expressed as a functions of s only, but for annihilations into scalars the amplitudes squared depend also on the Mandelstam variable t. Therefore respectiveσ(s) needs to be written aŝ where f (s, t) is the respective sum of all amplitudes squared for given i, j and t can be expressed by s and θ using

JHEP09(2015)162
In the nonrelativistic limit, using the formula (A.2), this tedious integration of f (s, t) can be avoided by changing the order of the limit and the integral [53], then note that this expression does not depend on cos θ, therefore the integration is trivial and gives factor 1. Similarly for the derivative we havê . (A.7) The presence of the last term comes from the fact, that t (s) = a(s) + b(s) cos θ, where a(4M 2 Z ) = −1/2 and b(4M 2 Z ) = ∞. On the other hand integral of cos θ vanishes, therefore one needs to calculate the second derivative to determine the limit. • Propagators

A.2 Relevant vertices
A.4 Z annihilation cross section formulae Squares of amplitudes and interference terms that contribute to the 4 diagrams contributing to the h i h j (i, j = 1, 2) final state shown above are listed below where f k stands for the square of the kth amplitude while f kl for 2 (f k f l ) (k, l = 1, 2, 3, 4).

B Model files for SARAH
In this appendix are included the model files used for SARAH [21][22][23][24], to study the model discussed in this paper. Currently SARAH does not implement a Z 2 symmetry for the imaginary scalar component to give S → S * , nor for the Z 2 for the U(1) X gauge field. Instead the kinetic mixing couplings have been set to vanish, g 1,x = g x,1 ≡ 0, to preserve this symmetry. Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.