Minimal semi-annihilating $\mathbb{Z}_N$ scalar dark matter

We study the dark matter from an inert doublet and a complex scalar singlet stabilized by $\mathbb{Z}_N$ symmetries. This field content is the minimal one that allows dimensionless semi-annihilation couplings for $N>2$. We consider explicitly the $\mathbb{Z}_3$ and $\mathbb{Z}_4$ cases and take into account constraints from perturbativity, unitarity, vacuum stability, necessity for the electroweak $\mathbb{Z}_N$ preserving vacuum to be the global minimum, electroweak precision tests, upper limits from direct detection and properties of the Higgs boson. Co-annihilation and semi-annihilation of dark sector particles as well as dark matter conversion significantly modify the cosmic abundance and direct detection phenomenology.

presence of the singlet is as essential, since it is impossible to allow only the λ 6 |H 1 | 2 H † 1 H 2 term of the two Higgs doublet model (2HDM) for semi-annihilation without also allowing the λ 7 |H 2 | 2 H † 1 H 2 term which mixes the two doublets. 1

Constraints on charge assignments
To allow the SM Yukawa terms of the Higgs H 1 , and to keep the DM stable, the discrete charges must satisfy certain requirements. On one hand, the H 1 Yukawa terms with u-and d-type quarks and charged leptons must separately have zero discrete charge modulo N . From a low energy point of view, however, we can simply set the charges of all standard model fields to zero. On the other hand, we want to forbid the |H 1 | 2 S and other terms that lead to mixing, 2 together with Yukawa couplings for H 2 . Therefore, the discrete charges must satisfy All possible scalar potentials contain a common piece V 0 since under any Z N and for any charge assignment each field can be paired with its Hermitian conjugate to form an invariant: (2.2)

Doublet-like DM
In case of mixing angle, θ, between the neutral components of the doublet and the singlet, dark matter can be either doublet-like or singlet-like. In the first case, if DM is degenerate with the next-to-lightest stable particle (NLSP) it can only contribute to a small fraction of the relic density Ω DM because of the large coannihilation contribution due to the ZH 0 A 0 coupling, where H 0 and A 0 are the real and imaginary neutral components of H 2 . If dark matter is singlet-like, the coannihilation cross section is suppressed by an additional factor of sin 2 θ . In this case even full degeneracy of DM and NLSP does not affect the relic density significantly.
In the case of the doublet-like DM, there are three ways to lift the degeneracy between the DM and the NLSP particle (we keep X 1 for generality): 1. A nonzero µ 2 S S 2 term (together with a µ SH S † H † 1 H 2 or µ SH SH † 1 H 2 mixing term to convey the mass gap from S to H 2 ). If the N in the Z N is odd, the term is not allowed due to X S > 0. For an even N , X S = N/2 is possible and, in addition, X 2 = N/2 + X 1 (or 1 ↔ 2) is needed.

A λ 5 (H †
1 H 2 ) 2 term. A λ 5 term can only be allowed for even N , because it transforms as ω 2(X 1 −X 2 ) and as X 1 = X 2 , the exponent 2(X 1 − X 2 ) cannot be zero modulo an odd number such as 3.
3. Both the µ SH S † H † 1 H 2 and µ SH SH † 1 H 2 terms. It is only possible to have both these terms for an even N with X S = N/2 and X 2 = N/2 + X 1 (or vice versa). In this case the µ 2 S and λ 5 terms are allowed as well.
None of these work for odd N , in which case the only possibility is singlet-like DM. The semi-annihilation terms λ S12 H † 1 H 2 S 2 and λ S21 H † 2 H 1 S 2 are forbidden if both µ SH and µ SH terms are allowed, nor can they coexist with the µ 2 S term. Thus, for doublet-like DM with semi-annihilation, the only option to lift degeneracy between the DM and NLSP particle is a λ 5 term.
The SM matter fields and a heavy neutrino singlet can be put in a 16 of SO (10), and the SM Higgs field in the 10 of SO (10). The minimal choice of representation to embed the complex singlet and the new doublet in is a scalar 16. Under SU (5) × U (1) X the representations decompose as 16 = 10 16 1 +5 16 −3 + 1 16 5 , S is the scalar analogue of the neutrino singlet in 1 16 5 and H † 2 is the scalar analogue of the lepton doublet in5 16 −3 . Therefore the U (1)-charges of the scalar sector are X S = 5, X 1 = 2 and X 2 = 3 (these are equal, modulo N , to the discrete Z N charges of the fields).
In flipped SU (5), the relation to the SM fields is 2 . The SM Higgs is H 1 ∈ 5 10 −2 . Note that the doublet must be flipped too with respect to standard SU (5).
S is the scalar analogue of the neutrino singlet in 10 16 1 and H † 2 is the scalar analogue of the lepton doublet in5 16 −3 . The U (1)-charges of the scalar sector are X S = 1, X 1 = −2 and X 2 = 3, which are equal, modulo N , to the discrete Z N charges of the fields.
For the sake of completeness, we include the unique scalar potential symmetric under Z 2 : For the given dark sector of H 2 and S, scalar potentials for higher Z N that do not contain semi-annihilation terms will be equivalent to the potential (3.1) with some terms set to zero. The Z 2 potential (3.1) in the case of SO(10) GUT in which some interactions are suppressed was studied in detail in [61][62][63][64][65]. Both H 2 and S are odd under Z 2 but only one of them is dark matter.

Z 3 potential with semi-annihilation
A Z 3 potential that induces semi-annihilation processes is invariant under e.g. the Z 3 charges X 1 = 0, X 2 = X S = 1. Another such potential is obtained from eq. (3.2) by substituting S → S † (with µ SH → µ SH and λ S12 → λ S21 ). From a low energy point of view, the two potentials are indistinguishable. Both the standard SU (5), in which case the Z 3 charges of fields are X S = 2, X 1 = 2, X 2 = 0, and flipped SU (5), with X S = 1, X 1 = 1, X 2 = 0, yield the potential obtained by S → S † in (3.2). Our Z 3 and Z 4 lagrangians are invariant under hypercharge symmetry H 1 → e iφ Y H 1 and H 2 → e iφ Y H 2 . We can use a hypercharge rotation of the doublets [74] to satisfy X 1 = 0 and X 2 > 0. For example, for standard SU (5), we can rotate the charges to X S = 1, X 1 = 0, X 2 = −1 ≡ 2 mod 3, which upon replacing H 1 and H 2 with their conjugates gives X S = 1, X 1 = 0 and X 2 = 1, which gives the Z 3 scalar potential (3.2).
The µ SH term in (3.2) induces mixing between the neutral components of H 2 and S. In terms of the mass eigenstates x 1 , x 2 , we have Considering the masses M h , M H ± , M x 1 , M x 2 and the mixing angle θ as free parameters of the model, we have

Z 4 potential with semi-annihilation
The only potential for Z 4 that contains semi-annihilation terms 3 is invariant under e.g. the assignment of Z 4 charges X 1 = 0, X 2 = 2, X S = 1. The dark sector particles do not mix with each other, because S and H 2 have different Z 4 charges. As a result this model has two dark sectors with the complex scalar S in the first one, and the second one comprising the charged Higgs boson H ± and the real scalars H 0 and A 0 . Any of the neutral particles with a non-zero Z 4 charge can be a dark matter candidate. Considering the masses of the scalars M 2 h , M H ± , M S , M H 0 and M A 0 as independent parameters, we have 14) (3. 16) 4 Experimental and theoretical constraints

Perturbativity
There are several possible definitions of perturbativity constraints on scalar quartic couplings. In [36], for example, it was required that the contribution of each coupling λ i to its own βfunction is less than unity so that the couplings do not run too fast. We follow [75], comparing the couplings λ i in the potential to the vertices in the Feynman rules for mass eigenstates. Barring accidental cancellations, the vertex factors have to be smaller than 4π to ensure that the one-loop level quantum corrections are smaller than the tree level contributions. For example, the quartic singlet self-interaction term λ S |S| 4 yields the vertex factor i 4λ S . Demanding that 4λ S < 4π gives λ S < π. If a quartic coupling λ i in the potential occurs in several vertices, we choose the strongest bound.

Perturbative unitarity
At high energy, the tree-level scalar-scalar scattering matrix is dominated by the quartic contact interaction terms. The s-wave scattering amplitudes should not exceed the perturbative unitarity limit for this partial wave, requiring that the eigenvalues of the S-matrix M must be smaller than the unitarity bound given by The unitarity bounds of the 2HDM were first studied in [76,77]. We will extend the formalism of [78,79] for the 2HDM to states containing the singlet S. The initial states are classified according to their total hypercharge Y (0, 1 or 2), weak isospin σ (0, 1 2 or 1) and discrete Z N charge X. The Z 3 unitarity bounds are reducible to those of the Z 4 case, and therefore we present only the latter.
For the sake of brevity, we list only the two sets of initial states which differ from the 2HDM initial states given in [78,79]. The full set of possible initial states with hypercharge Y = 1 and σ = 1 2 is The full set of possible initial states with hypercharge Y = 0 and σ = 0 is where the first three states have discrete Z 4 charge X = 0 and the last four states have X = 2. We do present all the scattering matrices and bounds on their eigenvalues for the Z 4 model. They reduce to the Z 3 case with λ S = 0, λ 5 = 0, λ S21 = 0. The scattering matrices The eigenvalues Λ X Y σi of the above scattering matrices (where i = ± or 1, 2, 3) can be written as Λ 0,1 |Λ 0 00 1,2,3 | Since Λ 0 00 1,2,3 are too cumbersome to be presented in full, in (4.14) we have given an upper bound on their absolute values by applying Samuelson's inequality [80] to the characteristic equation (in our numerical calculations we use the exact eigenvalues). The inequality arises from the observation that a collection of n points is within √ n − 1 standard deviations of their mean. For the polynomial a n x n + a n−1 x n−1 + . . . + a 1 x + a 0 with only real roots, the roots lie in the interval bounded by a n a n−2 . (4.17)

Vacuum stability
In order to have a finite minimum of the potential energy, the scalar potential has to be bounded below, especially in the limit of large field values. The quadratic and cubic terms are negligible in this limit and therefore it suffices to consider only the quartic terms to find the constraints of vacuum stability. To ensure that the quartic potential is bounded below, we write the matrix of quartic interactions in a basis of non-negative field variables and demand this matrix to be copositive [81]. The Higgs doublet bilinears can be parameterised as [82] The parameter |ρ| ∈ [0, 1] as implied by the Cauchy inequality 0 |H † 1 H 2 | |H 1 ||H 2 |. The singlet can be written in the polar form as S = se iφ S .
Then the quartic part of the Z 4 -symmetric potential (3.10) takes the form We discuss only the Z 4 case as its vacuum stability conditions reduce to the Z 3 case with λ S = 0, λ 5 = 0, λ S21 = 0 and further cos(φ + 2φ S ) = −1 so that the λ S12 term can always be chosen negative in (3.2). The parameters r 2 1 , r 2 2 and s 2 are non-negative and can be used as a basis for the matrix of quartic couplings. If the potential has semi-annihilation terms, then it contains terms with r 1 r 2 besides r 2 1 , r 2 2 . This leads to an ambiguity in the matrix because of r 2 1 r 2 2 = (r 1 r 2 ) 2 . In that case we define r 1 = r cos γ, r 2 = r sin γ, where 0 γ π 2 is a free parameter. In the vacuum stability conditions, the potential must be minimised with respect to all free parameters such as φ, φ S or γ.

Globality of the Z N -symmetric vacuum
The Z N -symmetric, EW-breaking vacuum has to be the global minimum of the scalar potential. 4 To check whether the correct vacuum is a global minimum, the solutions or a given point in the parameter space will be substituted back in the scalar potential and compared to each other. Below we give a classification of the stationary points of the scalar potentials with the field content H 1 , H 2 , S. In the general 2HDM there are three possible forms of vacua due to SU (2) invariance. Via SU (2) L × U (1) Y transformations, the vacuum expectation values can always be reduced to the form where v 1,2,+ are real and θ = 0 when v + = 0 (see [83]). In the classification of the extrema, we extend the notation used for the inert doublet model in [84] by adding the index S to the symbol when a singlet VEV is present (if only the singlet VEV breaks Z N , we prefix the name of the vacuum by 'non-').
In the fully symmetric vacuum (EW) all VEVs are zero, while in the EW-symmetric EW S , Z N is broken by the singlet VEV. In the inert vacuum I 1 , only the SM Higgs H 1 gets a VEV and the neutral component of H 2 can be DM candidate -this is the vacuum which we require to be the global minimum of the potential. Z N is broken by v S in the non-inert I 1S . Table 1. Classification of the stationary points of the model.

CP-violating CP S
In the inert-like I 2 , the rôles of H 1 and H 2 are reversed [84], except for the Yukawa couplings of H 1 that break Z N at loop level; in the non-inert-like I 2S , Z N is already broken at tree level by v S . In the mixed vacuum M both doublets get a VEV, and in M S the singlet gets a VEV too. To have a charge-breaking vacuum, all the v 1 , v 2 and v + have to be nonzero and v S can be zero (CB) or not (CB S ). In a CP-violating extremum CP S there is a phase of θ between the VEVs of H 2 and H 1 .
In the 2HDM, the existence of a normal extremum I 1 , I 2 or M that is a minimum implies that it has lower potential energy than the charge-breaking or CP-violating vacua [85,86]. It is not necessarily true when the model is extended with the singlet, since in a vacuum with a singlet VEV the effective doublet mass terms can be different. Also note that the CP-violating vacuum CP S does not exist in the pure inert doublet model and is enabled only by v S = 0: the singlet VEV generates the necessary mixing term for H 1 and H 2 .
Actual calculations are much simplified if in the scalar potential the doublet bilinears H † i H j are expressed in terms of gauge orbit variables [87] (see also [88,89]). All the bilinears can be arranged into the Hermitian 2 × 2 matrix which is decomposed as where σ a are the Pauli matrices. The four real gauge orbit variables are Positive semidefiniteness of the matrix K implies the 'future light cone' conditions Inverting (4.33), the potential can be written in terms of K µ . Each term in the potential is at most quadratic in K µ which reduces the degree of the minimisation equations. Of the solutions to the equations, only the stationary points that satisfy (4.34) are physical. In terms of the doublet VEVs (4.30), the VEVs of the gauge orbit variables are given by One can see that the condition that the doublet VEVs preserve Z N is The vacua EW and EW S are in the tip K µ = 0 of the doublet light cone. If we choose µ 2 1 = −M 2 h /2, then the SM Higgs mass is always negative at the tip and this point is by construction never a minimum, but a saddle point.
The charge-breaking vacua are inside the forward light cone: This point is a minimum if in the basis (K 0 , K 1 , K 2 , K 3 , ReS, ImS), the leading principal minors of the Hessian are all positive. The vacua I 1 , I 1S , I 2 , I 2S , M, M S and CP S , where the full electroweak gauge group SU (2) L × U (1) Y is broken into U (1) EM , are on the null surface of the future light cone: The latter condition is enforced by adding to the potential a Lagrange multiplier term The inequality K 0 > 0 in (4.37) has to be checked separately. This point is a minimum if u > 0 and the last five leading principal minors of the bordered Hessian in the basis (u, K 0 , K 1 , K 2 , K 3 , ReS, ImS) are all negative. The vacuum expectation values can be calculated in analytical form in most of the stationary points. For the extrema where the VEVs of both the doublets and the singlet are all non-zero, we solve the minimisation equations numerically with the PHCpack equation solver [90].
The solutions can then be substituted back in the scalar potential, and the global minimum be found by the smallest value. We require the inert vacuum I 1 to be the global minimum. In addition, in order for the stationary point to be a (local) minimum, the scalar masses or eigenvalues of the Hessian matrix at the stationary point have to be positive.
This requirement limits the size of µ S in the Z 3 model as in [24] and requires µ 2 To ensure that we are in the inert and not in the inert-like vacuum, we must have (4.40)

Electroweak precision tests
The measurements of electroweak precision data put strong constraints on physics beyond the SM. The latest electroweak fit by the Gfitter group [91] gives for the oblique parameters S and T the central values with a correlation coefficient of +0.89.
To calculate electroweak precision parameters S and T , we use the results for general models with doublets and singlets [92,93]. The usual loop functions are defined as with F (I, I) = 0 in the limit of J → I, and In terms of these functions, the non standard model contributions to the S and T parameters for the Z 3 invariant potential (3.2) are and where in the limit of vanishing θ only the middle term survives. We see that the T parameter is in general positive and that M H ± cannot be too different from M x 2 .
For the Z 4 invariant potential (3.10), the complex singlet does not mix with the neutral components of the doublet and does not contribute to the EWPT at 1-loop level. The S parameter is For the T parameter, we reproduce the result of [22]:

LEP limits
The results of precision measurements at the Large Electron-Positron Collider (LEP) exclude decays of the SM Z and W ± bosons into invisible particles, requiring [94,95]

Higgs diphoton signal and invisible decays
In these models the Higgs couplings are SM-like, except for the radiatively generated diphoton coupling which can receive a contribution from the charged Higgs. Modifications to the h → γγ rate are expected to be large only for a light charged Higgs, as in the inert doublet model [95,[98][99][100][101]. The fit to the latest experimental data from TeVatron [102], ATLAS [103][104][105][106][107][108] and CMS [109][110][111][112][113][114][115][116][117][118][119][120][121] gives for the diphoton rate R γγ = 1.06 ± 0.10, (4.51) if all the other rates are fixed to their SM values [122]. Furthermore there is the possibility of invisible Higgs decays with the 125 GeV Higgs decaying into the singlet or scalar/pseudoscalar components of the doublet: the invisible branching ratio is BR inv < 0.24 at 95% C.L. [122,123]. In the Z 3 model this range is most likely ruled out by direct detection as seen previously in the Z 3 singlet dark matter model [24]), while a large invisible rate can be generated in the Z 4 model as will be discussed below.

Cosmic density of dark matter
The PLANCK collaboration has recently released results for the cosmological parameters, in particular for the DM relic density [1]. When averaged with the WMAP-9 year data [124], it leads to the very precise value Ωh 2 = 0.1199 ± 0.0027. (4.52) We will use the 3σ range below.
To compute the relic density in the Z 3 model, we use micrOMEGAs_3.5 which takes into account all annihilation, coannihilation and semi-annihilation channels [125]. Final states with virtual gauge bosons that can be present in (co-)annihilation and semi-annihilation processes are also included.
In the Z 4 model there can be two dark matter candidates: the singlet with X = 1 and the lightest component of the doublet, H 0 , A 0 , with X = 2. To compute the relic density, we use the generalized equations for the abundance Y i = n i /s: where σ abcd v stands for the thermally averaged cross section σv for the reactions ab → cd, a, b, c, d = 0, 1, 2 represent any particle with a given X (SM particles have X = 0), and Y a are the equilibrium abundances. In σ abcd v all annihilation and coannihilation processes are taken into account as well as annihilation into virtual gauge bosons. Semi-annihilation processes include all those where 2 DM particles annihilate into one DM and one standard particle, specifically σ 1110 v , σ 2220 v , σ 1120 v and σ 1210 v . The method of solution for these equations as implemented in micrOMEGAs is described in [25] and [126].
The abundances Y 1 and Y 2 will be modified by the interactions between the two dark sectors. After the light DM freezes-out interactions such as hh → ll lead to a decrease of the abundance of the heavy component h and to an increase in the light component l . Such is also the case for semi-annihilation processes of the type hh → l 0 (or its reverse h0 → ll ). The semi-annihilation process 12 → 10 has no influence on Y 1 while leading to a decrease of Y 2 . Note that this process is always kinematically open unless the only SM particle in the final state is heavier than H 0 , A 0 . Finally, semi-annihilations involving only particles of a given sector always lead to a decrease of the abundance of the corresponding dark matter.

Dark matter direct detection
The best upper limit on the spin independent (SI) scattering cross section on nuclei has been obtained by the LUX experiment [9]: σ SI < 7.6×10 −46 cm 2 for M DM = 33 GeV. We also show the results from XENON100 (2012) [8]. Future detectors such as SuperCDMS(SNOLAB) [127], XENON1T [128] or LZD [129] will increase the sensitivity by one to four orders of magnitude.
To compute the model predictions for the SI cross section we use micrOMEGAs_3.5, and we assume that both DM and anti-DM have the same local density. In the Z 3 model, the DM candidate x 1 a mixture of the complex doublet and singlet scalar and there can be large differences in the scattering rate on protons and neutrons. To compare directly with the experimental limits, which are obtained assuming isospin conservation, we compute the normalised cross section of DM on a point-like nucleus (that we take to be xenon) where x denotes the DM candidate, µ x the reduced mass, f p , f n the amplitudes for protons and neutrons, and the average over x and x * is assumed implicitly.
In the Z 4 model with two dark matter candidates, we also compute the normalised cross section on xenon for each dark matter candidate after rescaling by the relative density of each component.

Dark matter indirect detection
The annihilation of DM in the Milky Way halo can lead to excesses in the cosmic ray fluxes (photons, positrons, antiprotons) that provide indirect evidence for DM. The measurements of the gamma-ray flux from Dwarf Spheroidal Galaxies by the Fermi-LAT satellite provide the most stringent constraint for light DM annihilating into bb or τ τ [130]. The canonical cross-section, σv ≈ 3 × 10 −26 cm 3 /s is ruled out for DM masses below 30 GeV [131]. For heavier DM, the measurements of the antiproton flux from PAMELA [132], for the MED set of propagation parameters [133] have roughly the same sensitivity as Fermi-LAT, with a limit of σv ≈ 10 −25 cm 3 /s for a DM mass of 200 GeV in the bb, τ τ or W W channels [134]. Since in our models the light masses are severely constrained by the Higgs invisible width we will in the numerical analysis consider only the PAMELA limit from antiprotons and we will assume the MED set of propagation parameters.
The measurements of the positron spectra by PAMELA and AMS are very powerful to constrain models which favour DM annihilation into e + e − or τ + τ − pairs, this is not the case in our models, we will therefore not consider these signatures.

Results for the Z 3 model
To explore the phenomenology of the model, we perform a random scan over the parameter space. The masses and the cubic term are generated with uniform distribution in the ranges 1 GeV < M x 1 < 1000 GeV, 1 GeV < M x 2 , M H ± < 2000 GeV, 0 GeV < µ S < 3500 GeV, and 124 GeV < M h < 127 GeV. The range of the mixing angle θ between the neutral components of H 2 and S is 0 θ π/2 and we choose M x 1 < M x 2 without loss of generality. In practice, we take the mixing angle in the range 0 θ 0.06 with uniform distribution. This guarantees that the DM candidate x 1 is dominantly singlet-like and so does not lead to a too large direct detection rate.
The quartic couplings are generated with triangular distribution (with the mode at zero) in the ranges allowed by perturbativity: 5 except for λ 4 from (3.9) and λ 1 from (3.8), the perturbativity of which is subsequently checked. We then apply the unitarity, vacuum stability and globality bounds, the upper limit on the Higgs invisible decays, the 3σ range for the DM relic density and the 3σ range for the electroweak precision parameters S and T . We present results for points that satisfy this set of constraints. We thten compute the Higgs diphoton signal strength as well as the DM direct detection rate and self-annihilation cross section relevant for indirect detection.
To emphasize the role of semi-annihilation on the DM properties in the model, we will present our results with a colour code characterising the fraction of semi-annihilation defined 5 Since the mixing is very small, we ignore it in the derivation of the perturbativity bounds. vσ xx→x * X vσ xx * →XX + 1 2 vσ xx→x * X , where x stands for x 1 , x 2 , H + and X for any SM particle. The dominant DM annihilation processes lead to gauge boson and Higgs pairs while the dominant semi-annihilation process is generally Other semi-annihilation processes such as When the singlet-doublet mixing is small these processes depend on λ S12 and/or µ S . For small mass differences between x 1 and x 2 (and/or H ± ) coannihilation occurs and new semiannihilation processes become possible, for example x 1 x 2 → Zx 1 , hx 1 or even In figure 1 we show the regions allowed by our set of constraints in the planes λ S1 , µ S and λ S12 vs. M x 1 . Basically it is possible to satisfy the Planck constraint for any mass of DM, while the upper bound on the Higgs invisible decay rules out all DM masses below ≈ 50 GeV. The same range of masses are also incompatible with the upper limit on the direct detection rate as will be discussed in Section 5.2. Similarly to the case of the pure Z 3 singlet DM [24], the relic density constraint determines the range of allowed values for λ S1 /M x1 , the combination of parameters that control DM annihilation, while smaller values of λ S1 are possible when semi-annihilation is important. Significant contributions from semiannihilation is associated with large values for µ S and/or λ S12 , the former contributing to the semi-annihilation process x 1 x 1 → x * 1 h (which is more important for relatively light DM masses), while the latter to processes with the dominantly doublet Higgs in the final state, these occur only if M x 2 ,H ± + M SM < 2M x 1 where SM refers to the scalar or gauge boson produced in the semi-annihilation process. The requirement that the SM vacuum be the global one constrains the possible maximum value of µ S . A few benchmarks satisfying these constraints are described in Appendix A. Figure 2 shows the results of the electroweak precision parameters S and T . Due to the fact that the allowed mixing angle is very small, the parameters virtually do not depend on the mass of the singlet-like DM. The constraint on the T parameter basically imposes |M H ± − M x 2 | 120 GeV on the mass splitting of H ± with the doublet-like neutral scalar x 2 . Therefore although it restricts the parameter space there is no direct correlation with other observables, in particular the spin-independent direct detection cross section σ SI . In the following figures we impose the 3σ constraint from the S and T parameters. Figure 3 shows the h → γγ signal strength (normalised to the SM) vs. the DM mass. The rate is systematically below the current average, and generally within the 2σ error band. Nevertheless there can be larger deviations from the SM for low DM masses. Note that here the constraint from invisible Higgs width has been applied.

Dark matter observables
In figure 4 we show the results of the DM spin-independent cross section σ SI vs. the DM mass M x 1 . The colour variation from black to pink (black to light grey) shows the fraction of semi-annihilation. The solid grey lines are the XENON100 (2012) [8] and LUX [9] exclusion limits at 90% C.L. and the dashed grey lines are the projected 90% C.L. exclusion limits for SuperCDMS(SNOLAB) [127], XENON1T [128] or LZD [129] (which will be at the limit of liquid xenon based experiments).  The current LUX upper limit severely constrains M x 1 < 120 GeV as well as points with a large doublet singlet mixing. As expected scenarios dominated by semi-annihilation lead  Figure 4. Spin-independent direct detection cross section on xenon, σ SI vs. M x1 for the Z 3 model. The solid grey lines are the XENON100 (2012) [8] and LUX [9] exclusion limits at 90% C.L. and the dashed grey lines are the projected 90% C.L. exclusion limits for SuperCDMS(SNOLAB) [127], XENON1T [128] or LZD [129]. The colour code shows the fraction of semi-annihilation α.
to suppressed cross sections. Note also the dark points with small semi-annihilation but low σ SI at high DM masses. They correspond to co-annihilation -that is, the relic density is dominated by self-annihilation of either x 2 or H ± and the cross section for annihilation of x 1 is small, resulting in a small cross section with nucleons as well. Some of the points will remain out of reach of Xenon1T. We also investigate the indirect DM signatures in this model. The annihilation channels for DM in the Galaxy are, as in the early Universe, often dominated by the bb, τ τ, W W, ZZ channels, with some contribution from the hh channel. The main new feature is the possibility of semi-annihilation with final states such as x 1 Z, x 1 h or even x 2 Z, x 2 h or H ∓ W ± .
In figure 5 we show the results of the DM annihilation cross section, σv for v ≈ 10 −3 c, the quantity relevant for indirect detection. We find that generally σv ≈ 3 (5)×10 −26 cm 3 /s when α → 0 (1), although the predictions can span several orders of magnitude. In particular σv is suppressed when coannihilation dominates (indeed coannihilation is relevant only for the computation of the relic density), while kinematic effects can lead to either a large enhancement or suppression of the (semi-)annihilation cross section. This occurs when M x 1 ≈ M x 2 /2. The cross section can be enhanced by more than one order of magnitude because of the near resonant x 2 exchange in the s-channel. Large suppression of σv can also be found when the thermal annihilation relevant for the relic density benefits from a resonance effect while the cross section at small velocities in the galaxy does not for kinematical reason.
The largest values of σv can potentially lead to a strong enhancement of the antimatter flux, in particular of the anti-proton flux. To ascertain the viability of all our scenarios in a quantitative way, we have performed a χ 2 fit to the data measured by PAMELA for energies between 10 − 100 GeV. To avoid large solar modulation effects we ignore the data at lower energies. For the background we assume the analytical parametrisation in Ref. [133] and determine the 95%C.L. allowed region by imposing that χ 2 = χ 2 background + 4. Despite the enhancement in the annihilation cross-section, we find that only a handful of points are constrained for the MED propagation parameters [133]. This is because semi-annihilation channels have only one SM particle in the final state decaying into anti-matter, hence a flux reduced by a factor 2 and a softer anti-matter spectra with a shape similar to the one for DM annihilation into a pair of SM particles (here typically Z or h). Furthermore the largest σv are found for DM masses above a few hundred GeV where indirect detection have a reduced sensitivity. Note that for MIN propagation parameters our results are always compatible with the background only hypothesis.

Renormalisation group running
The interaction couplings depend on the energy scale via renormalisation group equations (RGEs). Due to the RGE running, above some scale Λ the model may become non-perturbative or the scalar potential may not be bounded from below. At energies greater than Λ, new physics, for example a Grand Unified Theory, is expected to appear to ensure that the full theory is perturbative and stable up to the Planck scale.
Since we are interested in the influence of semi-annihilation on the relic density and direct detection, we show in figure 6 the bounds on the λ S1 vs. λ S , λ S12 vs. λ S and λ S12 vs. λ S1 planes at the EW scale. 6 All the other couplings are set to zero with the exception of, obviously, λ 1 , and 0 λ 3 , λ S2 0.6 which are chosen to roughly maximise the scale of vacuum stability of the model in each point (e.g. prevent the Higgs quartic coupling from running to negative values). We use the SM two-loop β-functions for the gauge couplings and the top Yukawa coupling and the one-loop β-functions for the scalar quartic couplings given in Appendix B. The white contour lines show the logarithm of the scale log 10 (Λ/1 GeV) where perturbativity is lost, while the black contour lines show the combined bound from loss of perturbativity and vacuum stability.
A large fraction of semi-annihilation is associated with small values of λ S1 , figure 6 shows that the model can then be valid up to the GUT scale provided λ S and λ S12 are not too large. In fact when semi-annihilation into doublet final states is dominant (recall that this requires x 2 not to be too heavy compared to x 1 ), the value of |λ S12 | that produces maximal semi-annihilation is about 0.5, therefore the model can be valid up to the GUT or Planck scale. For light enough DM, large semi-annihilation rather arises from the cubic µ S term in which case a large value of λ S is needed in order to for the SM vacuum to be the global minimum of the potential (see [24]) and perturbativity will be lost close to the TeV scale.

Results for the Z 4 model
We perform a random scan over the parameter space. The masses are generated with uniform distribution in the ranges 1 GeV < M S , M H 0 , M A 0 < 1000 GeV, 1 GeV < M H ± < 2000 GeV, and the Higgs mass 124 GeV < M h < 127 GeV. We consider only the cases when M H 0 ,A 0 < 2M S , since otherwise the neutral doublet would decay before the freeze-out of S and the situation would therefore be analogous to the one particle DM model. Here and below, M H 0 ,A 0 ≡ min(M H 0 , M A 0 ) stands for the mass of the lighter neutral component of the doublet H 2 . The quartic couplings are generated with triangular distribution (with the mode at zero) in the ranges allowed by perturbativity: |λ S1 | < 4π, |λ S2 | < 4π, |λ S12 | < 4π, |λ S21 | < 4π, |λ S12 − λ S21 | < 2π, (6.1) except for λ 1 , λ 4 and λ 5 which are derived from the free parameters, (3.14), (3.15) and (3.16), for these the perturbativity condition is subsequently checked. We then apply the unitarity, vacuum stability and globality bounds, the upper limit on the Higgs invisible decays, the 3σ range for the DM relic density and the 3σ range for the electroweak precision parameters S and T . We present results for points that satisfy this set of constraints. We compute the Higgs diphoton signal strength as well as the DM direct detection rate and self-annihilation cross section relevant for indirect detection.  Figure 7. Electroweak precision parameters ∆T vs. ∆S for the Z 4 model. We show the 1, 2, and 3 σ contours [91]. The colour code shows Ω 2 /(Ω 1 + Ω 2 ). Figure 7 shows the results of the electroweak precision parameters S and T . Since there is no mixing between the singlet and doublet components, the S and T parameters depend only on the latter. They tend to restrict the three mass splittings between H 0 , A 0 and H ± , but cancellations in the parameters can occur for specific combinations of masses. Restricting S and T to the 3σ range, however, does not directly exclude any specific region on e.g. the direct detection plots. Figure 8 shows the h → γγ signal strength (normalised to the SM) vs. M H 0 ,A 0 . For higher masses the rate is between 0.9 − 1.0, but higher or lower rates are possible for doublet masses below 150 GeV. The allowed rates are similar to those of the inert doublet model (see [101]).

Dark matter observables
We find that most of the time the relic density is dominated by one or the other of the dark sector. First note that for the doublet DM -as in the inert doublet model -annihilation channels into gauge bosons are efficient and the relic density is typically too small unless the DM mass is below the mass of W or above 500 GeV. For the singlet component the relic density constraint can be satisfied over the full mass range [25] when ignoring the interactions between the two dark sectors. The interactions between the two dark sectors influence this picture. Figure 9 shows the allowed region in the M H 0 ,A 0 − M S plane, the color code indicate which is the dominant component from blue (singlet dominated) to orange (equal contribution) to green (doublet dominated). In the different regions we have indicated some benchmark points, numbered 1 to 6, which are also shown in plots below and described in Appendix A. Typically the lighter component is dominant.
When S is the heavier DM component, both DM conversion and semi-annihilation, notably SS → Hh, lead to a decrease in Y 1 (see (4.53)), the DM is therefore usually dominated by H 0 , A 0 , and falls into two regions one with M H 0 ,A 0 < 100 GeV the other with M H 0 ,A 0 > 500 GeV. These regions correspond roughly to the medium and high mass ranges of the inert doublet model. However in some cases when the relic abundance of S falls in the measured range, then the subdominant DM component, H 0 , A 0 can have any mass, see the points below the diagonal in figure 9.
When H 0 , A 0 are heavier than S, both self-interactions and semi-annihilations tend to reduce Y 2 , thus the relic density is typically dominated by the singlet. The lion's share of the points are quite similar to the model with a single complex singlet. As a result the doublet  mass can span the whole range up to M H 0 ,A 0 < 2M S . Note that for both components to contribute significantly (orange points in figure 9), their mass must be of the same order, therefore they are restricted to the region M S , M H 0 ,A 0 < 80 GeV and to the region between 300 GeV and 1000 GeV, otherwise the doublet component would be sub-dominant. Annihilation of S depends on the λ S1 coupling and annihilation of H 0 , A 0 respectively on λ L,R = λ 3 + λ 4 ± λ 5 . Dark matter conversion is mediated by the λ S2 coupling and the semi-annihilation couplings for H 0 , A 0 are λ ShH 0 ,ShA 0 = λ S12 ± λ S21 . In general, if the lighter component dominates, the absolute value of its annihilation coupling is less than 0.5, its (semi-)annihilation couplings and λ S2 have a large range, while for the heavier component to dominate, the latter must also be within 0.5.
In figure 10 we show the predictions for the SI cross section for each DM component independently. The cross section has been rescaled by a factor Ω i /(Ω 1 + Ω 2 ) to take into account the relative abundance of each DM component. On top, the SI cross section on protons for the doublet component is displayed. As in the inert doublet model most of the points where the DM mass is near the electroweak scale have a large cross section -exceeding the current bounds -this is so even when the doublet is a subdominant DM component (blue points). Only a few points around M Z /2 or M h /2 remain below the current bound. For heavier DM the current limits are generally satisfied. It is intriguing that we also find scenarios where the doublet leads to a cross section detectable by the next generation of detectors even if it forms a subdominant DM component. This means that direct detection experiments could detect two DM components of different masses in this model.
In figure 10 (bottom) the SI cross section on xenon for the singlet component is displayed. As in the singlet DM model, the cross section exceeds the current limit for masses below ≈ 100 GeV except for a few cases where M S ≈ M h /2. For higher masses the predictions are below the current limit but mostly in a range accessible by future detectors such as Super-  Figure 10. Normalised spin independent cross section on xenon for H 0 , A 0 (top) and S (bottom). The solid grey lines are the XENON100 (2012) [8] and LUX [9] exclusion limits at 90% C.L. and the dashed grey lines are the projected 90% C.L. exclusion limits for SuperCDMS(SNOLAB) [127], XENON1T [128] or LZD [129]. The colour code shows Ω 2 /(Ω 1 + Ω 2 ). CDMS or XENON1T. The lowest cross sections can be related to strong semi-annihilation in the singlet sector or simply to the fact that the abundance of the singlet component is small (recall that points in green are doublet-dominated). Comparing both figures one can see that points where both component have a comparable abundance typically lead to similar cross sections, see e.g. points 1 and 5.  Figure 11. Dominant spin independent cross section on xenon for either S or H 0 , A 0 . The solid grey lines are the XENON100 (2012) [8] and LUX [9] exclusion limits at 90% C.L. and the dashed grey lines are the projected 90% C.L. exclusion limits for SuperCDMS(SNOLAB) [127], XENON1T [128] or LZD [129]. The colour code shows Ω 2 /(Ω 1 + Ω 2 ).
These two figures might give the impression that for many input parameters, DM would escape direct detection. However, figure 11 which displays for each point in parameter space the value of the largest σ S SI or σ H 0 ,A 0 SI , clearly shows that most of the model parameter space leads to a signal that could be detected in the near future. Only some points remain out of reach of the projected sensitivity of Xenon1T or even LZD.
In this model the σv relevant for indirect detection can be very large for either DM candidate, however after rescaling by the fraction of DM density for each component, we obtain σv ≈ 3 × 10 −26 cm 3 /s. Therefore the limits from PAMELA on antiprotons are easily satisfied.

Renormalisation group running
Without loss of generality, we consider here only H 0 as the doublet component of DM, because the RGEs are symmetric under the exchange of λ L with λ R and λ ShH 0 with λ ShA 0 . There are four relevant parameters: λ S1 , λ L , λ S2 and λ ShH 0 . We show the bounds on various parameter planes at the EW scale in figure 12. The values of other couplings are set to zero, except λ 1 , and λ 2 , λ 3 and λ S which are chosen to roughly maximise the scale of vacuum stability of the model.
Each subplot can be thought of as relevant to a limiting case where only two of the processes of annihilation of S and H 0 , DM conversion and semi-annihilation are relevant.
The top left panel corresponds to the bounds in the λ L vs. λ S1 plane in the special case where interaction between the two sectors is negligible, but independent annihilation of S and H 0 ensures the correct relic density. The bottom left panel shows the bounds in the λ ShH 0 vs. λ S1 plane for the case where S is lighter and H 0 decays into it via semi-annihilation; the bottom right panel shows the bounds in the λ ShH 0 vs. λ S1 plane for the similar case where H 0 is lighter and S decays into it via semi-annihilation.
For these special cases, except if S is lighter and H 0 decays into it via semi-annihilation (bottom left), one can find points with realistic relic density where the model is valid up to the GUT scale. However, for generic points where all the quartic couplings are of O(1), for most of the points the model loses perturbativity at about the TeV scale.

Conclusions
We have explored the phenomenology of an inert doublet and complex scalar dark matter model stabilized by Z N symmetries, with explicit investigation of the Z 3 and Z 4 cases. The new feature of these models as compared to the Z 2 case is the possibility of semi-annihilation and dark matter conversion. This has important consequences for all dark matter observables.
In the Z 3 model, semi-annihilation processes, e.g. x 1 x 1 → x 1 h, can give the dominant contribution to the relic abundance through the cubic (µ S S 3 ) or quartic (λ S12 S 2 H † 1 H 2 ) terms in the scalar potential. This means that the λ S1 parameter which sets the coupling of DM to the Higgs and thus the direct detection cross section is not uniquely determined by the relic density constraint as occurs in the Z 2 model. Large semi-annihilation is therefore associated with suppressed direct detection rate. While the bulk of the points will be testable by ton-scale detectors, it is possible to satisfy the constraints from vacuum stability and globality of the minimum of the potential with very small values of λ S1 -hence to escape all future searches, in particular when the DM is near the TeV range. The direct detection limits from LUX almost completely rule out the region where dark matter masses are below 120 GeV since for kinematic reasons the semi-annihilation does not play an important rôle (the Higgs cannot be produced in the final state). In this model, because there can be resonance enhancement of the annihilation cross section in the Galaxy when the mass of the DM is tuned to be half the mass of the doublet Higgs, the indirect detection cross section can be much enhanced as compared to the canonical value. The semi-annihilation processes will however lead to softer spectra since the DM particle in the final state drains part of the energy of the reaction. Furthermore we have shown that the model can be perturbative up to the GUT scale even with a large fraction of semi-annihilation.
Enlarging the symmetry to Z 4 entails two dark sectors, hence two dark matter candidates: a singlet and a doublet. In this case both semi-annihilation and dark matter conversion significantly affects the dark matter phenomenology of the model. While this model shares many characteristics of the inert doublet model especially when interactions between the two dark sectors can be ignored, the presence of the singlet dark matter candidate means that the doublet DM could only contribute to a fraction of the relic density (and vice versa). This means in particular that the doublet DM can have any mass instead of being confined to be at the electroweak scale or heavier than 500 GeV as in the inert doublet model.
We found that for the sub-dominant dark matter component, it is possible to have a detectable signal in future direct detection experiments even after taking into account the fraction of each component in the DM density. This occurs in particular when the sub-dominant component is the doublet since it typically has a large direct detection rate. Furthermore in some cases a detectable signal in future ton-scale experiments is predicted for each DM component, opening up the exciting possibility of discovering two DM particles.