Mitigating Direct Detection Bounds in Non-minimal Higgs Portal Scalar Dark Matter Models

TL;DR

This paper explores non-minimal Higgs portal scalar dark matter models with extended symmetries, demonstrating they can evade current direct detection bounds while still explaining observed dark matter relic density.

Contribution

It systematically studies how enlarging the stabilizing symmetry and including multipartite features relaxes direct detection constraints in scalar dark matter models.

Findings

Non-minimal models relax direct detection constraints.

Resonant semi-annihilation opens new parameter space.

Multi-component frameworks illustrate complex dark matter interactions.

Abstract

Minimal scalar Higgs portal dark matter model is increasingly in tension with recent results form direct detection experiments like LUX and XENON. In this paper we make a systematic study of minimal extension of the $ \mathbb{Z}_2$ stabilised singlet scalar Higgs portal scenario in terms of their prospects at direct detection experiments. We consider both enlarging the stabilising symmetry to $\mathbb{Z}_3$ and incorporating multipartite features in the dark sector. We demonstrate that in these non-minimal models the interplay of annihilation, co-annihilation and semi-annihilation processes considerably relax constraints from present and proposed direct detection experiments while simultaneously saturating observed dark matter relic density. We explore in particular the resonant semi-annihilation channel within the multipartite $\mathbb{Z}_3$ framework which results in new unexplored regions of parameter space that would be difficult to constrain by direct detection experiments in the near future. The role of dark matter exchange processes within multi-component $\mathbb{Z}_3 \times \mathbb{Z}_3'$ framework is illustrated. We make quantitative estimates to elucidate the role of the various annihilation processes in the different allowed regions of parameter space of these models.

Tables

Table 1. Table 1: Tree level unitarity constraints for discussed models.

Model	Lagrangian	Tree level unitarity constraints
Two particle under ${\mathbb{Z}}_2$	Eq. 2.1	where $x$ is the solution of the equation 27.
One particle under ${\mathbb{Z}}_3$	Eq. 8
Two particle under ${\mathbb{Z}}_3$	Eq. 3.1	where $y$ is the solution of the equation 29.
${\mathbb{Z}}_3\times {\mathbb{Z}}_3^{{}^{\prime }}$	Eq. 15	where $z$ is the solution of the equation 31.

Equations

$$-{\mathcal{L}}_{\mathrm{DM}-\mathrm{Higgs}}$$

$$-{\mathcal{L}}_{DM-Higgs}$$

$$\frac{dn}{dt}+3Hn$$

$$⟨\sigma v{⟩}_{eff}=⟨\sigma v{⟩}_{{\varphi }_1{\varphi }_1\to SM}+⟨\sigma v{⟩}_{{\varphi }_1{\varphi }_2\to SM}(1+\frac{\Delta m}{m_1}{)}^{3/2}{e}^{-\frac{\Delta m}{T}}.$$

$${\sigma }_n^{SI}=\frac{{\lambda }_{1h}^2{f}_n^2}{4\pi }\frac{{\mu }_n^2{m}_n^2}{m_h^4{m}_{{\varphi }_i}^2},$$

$$1\mathrm{GeV}<m_{{\varphi }_1}<1000\mathrm{GeV},2\le \Delta m\equiv m_{{\varphi }_2}-m_{{\varphi }_1}\le 1000\mathrm{GeV},$$

$$\displaystyle 0.001\leq\lambda_{1h}=\lambda_{2h}\leq 1,~{}0\leq\lambda_{12h}\leq 1~{}.\$$

$$-$$

$$m_{{\varphi }_1}=200\mathrm{GeV},\mathrm{\Delta }{\mathrm{m}}_{\mathrm{i}}\equiv {\mathrm{m}}_{{\varphi }_{\mathrm{i}}}-{\mathrm{m}}_{{\varphi }_1}=10\mathrm{GeV},0.001\le {\lambda }_{1h}={\lambda }_{\mathrm{ih}}\le 0.1,{\lambda }_{1ih}\le 0.112$$

$$-{\mathcal{L}}_{DM-Higgs}$$

$$\frac{d{n}_{{\varphi }_1}}{dt}+3H{n}_{{\varphi }_1}=-⟨\sigma v{⟩}_{{\varphi }_1{\varphi }_1\to SM}(n_{{\varphi }_1}^2-{n_{{\varphi }_1}^{eq}}^2)-\frac{1}{2}⟨\sigma v{⟩}_{{\varphi }_1{\varphi }_1\to {\varphi }_{1}SM}(n_{{\varphi }_1}^2-n_{{\varphi }_1}{n}_{{\varphi }_1}^{eq}),$$

$$50\mathrm{GeV}\le m_{{\varphi }_1}\le 1000\mathrm{GeV},0.001\le {\lambda }_{1h}\le 1.0,{\mu }_1\le 7.4m_{{\varphi }_1}$$

$$-$$

$$\frac{d{n}_{{\varphi }_1}}{dt}+3H{n}_{{\varphi }_1}$$

$${\lambda }_{1h}={\lambda }_{2h}=0.001,0.0001\le {\lambda }_{12h}\le 0.01,{\mu }_i\le 2m_{{\varphi }_i},{\mu }_{12}\le 2m_{{\varphi }_1}.$$

$$10\mathrm{GeV}<m_{{\varphi }_1}<1000\mathrm{GeV},2\le \Delta m=m_{{\varphi }_2}-m_{{\varphi }_1}\le 1000\mathrm{GeV},$$

$${\mu }_i\le 2m_i,{\mu }_{12}\le 2m_1,0.001\le {\lambda }_{1h}={\lambda }_{2h}\le 1,0.1\le {\lambda }_{12h}\le 1.$$

$$-{\mathcal{L}}_{DM-Higgs}$$

$$\frac{d{n}_{{\varphi }_1}}{dt}+3H{n}_{{\varphi }_1}$$

$$\frac{d{n}_{{\varphi }_2}}{dt}+3H{n}_{{\varphi }_2}$$

$${\sigma }_{eff}^i=(\frac{{\Omega }_i}{{\Omega }_T})({\sigma }_n^{SI}{)}_i=\frac{{\Omega }_i}{{\Omega }_T}\frac{{\lambda }_{ih}^2{f}_n^2}{4\pi }\frac{{\mu }_n^2{m}_n^2}{m_h^4{m}_{{\varphi }_i}^2}(i=1,2),$$

$$125\le m_1\le 500,500\le m_2\le 1000,5m_i\le {\mu }_1={\mu }_2\le 7.4m_i,$$

$$0.01\le {\lambda }_{1h}={\lambda }_{2h}\le 0.1,0.5\le {\lambda }_e\le 1.5$$

$$V({\varphi }_1,{\varphi }_2,H)={\lambda }_H|H{|}^4+\mathsf{M}({\varphi }_1,{\varphi }_2)|H{|}^2+v({\varphi }_1,{\varphi }_2),$$

$$\mathsf{M}({\varphi }_1,{\varphi }_2)=\frac{1}{2}{\lambda }_{1h}{\varphi }_1^2+\frac{1}{2}{\lambda }_{2h}{\varphi }_2^2+{\lambda }_{12h}{\varphi }_1{\varphi }_2$$

$$v({\varphi }_1,{\varphi }_2)=\frac{{\lambda }_{e_1}}{4}{\varphi }_1^2{\varphi }_2^2+\frac{{\lambda }_{e_2}}{3!}{\varphi }_1^3{\varphi }_2+\frac{{\lambda }_{e_3}}{3!}{\varphi }_1{\varphi }_2^3+\frac{{\lambda }_{1s}}{4!}{\varphi }_1^4+\frac{{\lambda }_{2s}}{4!}{\varphi }_2^4$$

$${\lambda }_{1s}>0,{\lambda }_{2s}>0,$$

$$6{\lambda }_{12h}{\lambda }_{1h}({\lambda }_{e_2}+{\lambda }_{e_3})+3{\lambda }_{e_1}{\lambda }_H{\lambda }_s-3{\lambda }_{12h}^2{\lambda }_s$$

$$\displaystyle\Gamma(\phi_{2}\to\phi_{1}\,X)=\frac{1}{16m_{\phi_{2}}(2\pi)^{3}}\frac{\lambda_{12h}^{2}m_{f}^{2}[(m_{\phi_{2}}-m_{\phi_{1}})^{2}-4m_{f}^{2}]}{m_{h}^{4}}\Big{[}\frac{m_{\phi_{1}}^{2}(m_{\phi_{2}}-m_{\phi_{1}})}{m_{\phi_{1}}}$$

$$\displaystyle+\frac{(m_{\phi_{2}}-m_{\phi_{1}})^{2}}{2}-m_{\phi_{1}}^{2}\ln(\frac{m_{\phi_{2}}}{m_{\phi_{1}}})\Big{]}.$$

$$\mathcal{M}(\theta )=16\pi \sum _{l=0}^{l=\infty }{a}_l(2l+1)P_l(\cos \theta )$$

$$|\mathcal{M}|<8\pi$$

$$w^{+}{w}^{-},\frac{hh}{\sqrt{2}},\frac{zz}{\sqrt{2}},\frac{{\varphi }_1{\varphi }_1}{\sqrt{2}},\frac{{\varphi }_2{\varphi }_2}{\sqrt{2}},{\varphi }_1{\varphi }_2,hz,h{\varphi }_1,h{\varphi }_2,z{\varphi }_1,z{\varphi }_2$$

Full text

1 Introduction

Existence of dark matter (DM) is supported from many astrophysical evidences like rotation curve of galaxies Rubin:1967msa ; Rubin:1970zza , anisotropies in CMBR Hu:2001bc and observations in bullet cluster Clowe:2006eq . It is a possibility that DM is particulate and may even have some non-gravitational interaction with the Standard Model (SM) sector. This gives the well discussed possibility of DM composed of thermal relic of cosmologically stable particles Jungman:1995df ; Munoz:2003gx ; Bertone:2004pz ; Bergstrom:2012fi ; Arcadi:2017kky .

The key feature of this paradigm is having a cosmologically stable DM candidate. This is usually ensured in particle physics models by invoking some symmetry arguments. The $\mathbb{Z}_{N}$ discrete symmetries provide the simplest realization of this stabilizing symmetry and is commonly employed in extensions of the SM Batell:2010bp . These symmetries can also arise as subgroup of broken continuous symmetry groups. While the manifest discrete symmetry in the Lagrangian typically prevent any decay it still allows number changing processes between the dark sector and the SM. These number changing processes are crucial in maintaining the thermal equilibrium between the two sectors in the early stages of evolution of the Universe. Once the Universe starts expanding these processes becomes less effective before finally stalling, leading to the standard framework for DM freeze out, leaving behind a relic density observable till the present epoch. The relic abundance has been precisely estimated from CMBR studies at WMAP Komatsu:2010fb and then at PLANCK Ade:2013zuv experiments, to be in the range $0.1133\leq\Omega_{\rm DM}h^{2}\leq 0.1189$.

The simplest number changing process that may be operative between the DM and SM sectors, allowed by the stabilizing symmetry, is the so called DM annihilation. Where typically two DM particles annihilate to produce two or more SM states111For other non-standard annihilation possibilities see Dolgov:1980uu ; Hochberg:2014dra ; Dey:2016qgf . This is effective in reducing the number density of the DM and lead to the freeze out. For weak scale mass and annihilation cross sections, these processes leads to the Weakly Interacting Massive Particle (WIMP) framework Lee:1977ua . However these same processes can be probed by the direct detection experiments Goodman:1984dc ; Lewin .

Non observation of DM in direct detection experiments provides some of the most stringent bound on the DM models constraining the annihilation processes as an effective mechanism to drive freeze out. This correlation provides the motivation for the large number of direct detection experiments that are in operation or have been proposed. The Large Underground Xenon (LUX) experiment is a dual-phase Xenon detector operating at the Sanford Underground Research Facility. First results of LUX Akerib:2013tjd set forth a minimum upper limit on WIMP-nucleon spin independent (SI) cross section of $7.6\times 10^{-46}~{}\rm cm^{2}$ at a WIMP mass of $33$ GeV. XENON is another experiment in operation at the Laboratori Nazionali del Gran Sasso, using ultra pure liquid Xenon as WIMP target. XENON100 experiment Aprile:2012nq gathered data for 13 months between 2011 and 2012, reaching a minimum sensitivity of $2\times 10^{-45}~{}\rm cm^{2}$ at a DM mass of $55$ GeV with $90\%$ confidence level. The upgraded XENON1T which acquired data for $34.2$ days have recently published first results Aprile:2017iyp already reaching sensitivity comparable to LUX and is expected to have increased sensitivity in near future with more data. Future proposals include the next generation XENONnT projected to achieve minimum spin-independent WIMP nucleon cross section $1.6\times 10^{-48}~{}\rm cm^{2}$ at WIMP mass of $50$ GeV Aprile:2015uzo . Dark matter WIMP search with liquid xenon (DARWIN) Aalbers:2016jon will be an experiment for the direct detection of DM using multi-ton liquid xenon. This experiments can be sensitive to Spin-independent DM-nucleon cross section of $2.5\times 10^{-49}~{}\rm cm^{2}$ at a WIMP mass of $40$ GeV Aalbers:2016jon . Note that collider searches at Large Hadron Collider (LHC) puts a considerably weaker bound on the intermediate mass (100-1000 GeV) Higgs portal dark matter candidates compared to the direct detection constraints Arcadi:2017kky .

These ongoing and proposed experiments mandates a closer look at various WIMP DM models and their prospects at these direct detection experiments. In this paper we will confine ourselves to a class of models where the SM singlet scalar dark sector communicates with the SM through a coupling with the Higgs. These, so called, Higgs portal models provides a simple framework for WIMP DM and are subject to extensive discussion in the literature Bertone:2004pz . Null results at direct detection experiments Akerib:2013tjd ; Aprile:2017iyp has already put a strong bound He:2016mls ; He:2008qm on the minimal Higgs portal models where a single scalar DM is stabilized by a discrete $\mathbb{Z}_{2}$ symmetry. Other than the extremely tuned Higgs mass pole region, the well motivated Beniwal:2017eik ; Artymowski:2016tme ; Vaskonen:2016yiu ; Cline:2017qpe minimal Higgs portal model has been pushed to a heavy DM mass region by XENON1T data Aprile:2017iyp , which can be further excluded by continued non observation in the immediate next generation experimental results. In this article we make a systematic study of the simple extension of this framework that can evade these constraints while remaining a viable DM candidate.

Possible augmentation of the minimal model can be done by enlarging the stabilizing symmetry group from $\mathbb{Z}_{2}$ to $\mathbb{Z}_{N}$ or by introducing multi-particle dark sector. These non-minimal models facilitate non-standard number changing channels like semi-annihilation D'Eramo:2010ep and co-annihilation Edsjo:1997bg , which cannot be explored at direct detection experiments. We find that within multipartite $\mathbb{Z}_{2}$ framework co-annihilation and mediated annihilation processes ameliorate some of the direct search bounds Casas:2017jjg ; Ghorbani:2014gka . Further we investigate multipartite $\mathbb{Z}_{3}$ model which have a significantly enriched DM phenomenology. Here the interplay of semi-annihilation and resonant semi-annihilation together with co-annihilation and annihilation processes uncover a large region of unexplored parameter space which satisfy both relic density and direct search bounds. Finally we explore the DM exchange processes in two component DM scenario with a $\mathbb{Z}_{3}\times\mathbb{Z}_{3}^{\prime}$ stabilizing symmetry. Interestingly, in certain regions of parameter space of this scenario both the DM states can be detected at the next generation experiments. We make an organized study of these frameworks taking each scenario one by one in increasing order of complexity. A detailed numerical scan of the parameter space is performed to explore the intricate interface of the various number changing process that are operative in a given framework.

The rest of the paper is organized as follows. In section 2 we present multi-particle $\mathbb{Z}_{2}$ framework. In section 3 we briefly review the impact of semi-annihilation in scalar DM models stabilized by a $\mathbb{Z}_{3}$ symmetry, before we present a rigorous study of multi-particle $\mathbb{Z}_{3}$ framework including a detailed discussion of resonant semi-annihilation. And in section 3.2, we focus two component DM under $\mathbb{Z}_{3}\times\mathbb{Z}_{3}^{\prime}$. In section 4 we briefly discuss tree level vaccum stability and unitarity constraints. Finally we conclude in section 5.

2 Minimal $\mathbb{Z}_{2}$ Model and Extensions

The operators facilitating the thermal freeze-out involve DM and SM fields which can be decomposed into $\mathcal{O}_{SM-DM}\sim\mathcal{O}_{DM}\mathcal{O}_{SM},$ assuming SM sector do not transform under stabilizing symmetry of the dark sector. The simplest renormalizable operator of such kind can be written by ideating the existence of a real scalar singlet DM ($\phi$) interacting to SM through Higgs portal interactions as $\phi^{2}H^{\dagger}H$ and has been studied exhaustively in the literature Lerner:2009xg . The Lagrangian can be written as,

[TABLE]

This setup implies the presence of an unbroken $\mathbb{Z}_{2}$ symmetry under which $\phi\to-\phi$ while all SM particles are even, making the field $\phi$ stable. The relevant annihilation processes that drive the freeze out are depicted in Fig. 1a. However, this minimal setup has been pushed to uncomfortable corner by recent results from direct search data in LUX 2016 Akerib:2016lao and XENON1T Aprile:2017iyp . The remaining unconstrained region are now confined to the tuned Higgs resonance region or for heavy DM mass $m_{\phi}\gtrsim 500$ GeV. The heavy DM region lies just below the present direct detection constraints and most of the allowed region will be explored in the next generation experiments. Admitting a tuning of the DM mass to be $m_{\phi}\approx m_{h}/2\sim 125/2$ GeV allows considerable relaxation of the direct search bounds while still satisfying the relic density calculations. Back of the envelope calculation show that at the resonance $\langle\sigma v\rangle\approx 24~{}\lambda^{2}_{h}~{}\rm GeV^{-2}$ assuming a Higgs width $\Gamma_{h}=4$ MeV Khachatryan:2016vau . This saturates the relic density bound when $\langle\sigma v\rangle$ is equal to $\sim 0.1$ pb. Estimating the Higgs portal coupling from here we find that the direct detection cross section can be as low as $1.69\times 10^{-52}~{}\rm cm^{2}$. Thus is expected to remain mostly unconstrained by direct detection experiments in near future.

A permissible Higgs portal DM with intermediate mass necessarily imply extensions of the this framework. The simplest possible generalization to this framework is to add more singlet scalars states to the dark sector. If they transform under different symmetries, say $\mathbb{Z}_{2}^{{}^{\prime}}\times\mathbb{Z}_{2}^{{}^{\prime\prime}}\cdots$ and so on, then all of them are stable leading to a multi-component DM framework. The impact of the DM-DM interactions in the two component framework under $\mathbb{Z}_{2}\times\mathbb{Z}_{2}^{{}^{\prime}}$ has been studied in Bhattacharya:2016ysw ; Karam:2016rsz . The dark sector exchange interaction play the role of a see-saw between the two components. While the lighter component behave like the single component framework with early detection possibility, the heavier component may have suppressed direct detection cross section. This suppression for the heavier component arises when its contribution to the total DM relic abundance is relatively small. When the DM masses and couplings to SM are same the symmetry is enlarged to $O(N)$ for $N$ component DM scenario. The exchange process now vanishes and requires all the DM components to have larger annihilation cross-section and larger DM-SM couplings making them more constrained by direct detection experiments Drozd:2011aa .

When the dark sector particles are charged under the same stabilizing symmetry, the heavier one can decay to the lighter component yielding a single component DM framework while the dark sector remains multipartite. The presences of extra particles open new non-standard number changing processes in the dark sector like co-annihilation and mediated annihilation. As will be detailed in the rest of this section this leads to considerable easing in many tensions of the minimal model.

2.1 Two particles under $\mathbb{Z}_{2}$ : Co-annihilation & Mediated annihilation

We will consider the case where the SM is augmented by two real scalar particles $\phi_{1}$ and $\phi_{2}$ odd under the same discrete $\mathbb{Z}_{2}$ symmetry and singlet under SM. The relevant part of the Lagrangian is given by,

[TABLE]

Note that the coupling $\lambda_{12h}$ arises as both of the components transform under same symmetry and give rise to the novel co-annihilation and mediated annihilation channels depicted in Fig. 1b and Fig. 1c respectively. Assuming without any loss of generality, $m_{\phi_{1}}<m_{\phi_{2}}$ the lightest mass state $\phi_{1}$ can be identified as the potential DM candidate222In principle we can have a term like $m^{2}\phi_{1}\phi_{2}$ in the Lagrangian which is allowed by all symmetries of the theory. This can lead to a mass mixing between the states. However, note that the Lagrangian in Eq. 2.1 is written in the mass basis of $\phi_{1}$ and $\phi_{2},$ assuming the mass matrix has been diagonalized. . The $\phi_{2}$ state can promptly decay to $\phi_{1},$ therefore do not effect the freeze out except through its contribution to the $\phi_{1}$ number changing process discussed above. This decay occurs through an off-shell Higgs and a schematic calculation can be found in Appendix B.

2.1.1 Relic density

The new set of Feynman graphs corresponding to the number changing processes of DM are shown in Fig. 1b and 1c. With a small annihilation cross section, as mandated by direct detection results, relic density bound can be saturated using the co-annihilation process depicted in Fig. 1b if the masses of the two sates are relatively degenerate. However, for large mass gap the $t$-channel mediated annihilation process in Fig. 1c takes up a major role in controlling DM relic density. These novel process are of interest because while they change the number density of the DM and thus aid freeze out they do not contribute to direct detection cross section. The co-annihilation channels do not contribute as the mass gap between the co-annihilating states kinematically forbid the corresponding direct detection process. The t-channel mediated annihilation with Higgs in the final state couple to the nucleon at 1-loop. So remains relatively unconstrained by direct detection.

Presence of $\lambda_{e_{i}}$ in Eq. 2.1 essentially influences the evolution of the DM density in two ways: (i) From the decay of $\phi_{2}$, DM is produced and (ii) DM scattering processes: $\phi_{2}\phi_{2}\to\phi_{1}\phi_{1},~{}\phi_{2}\phi_{2}\to\phi_{2}\phi_{1},$ and $\phi_{2}\phi_{1}\to\phi_{1}\phi_{1}.$ However, due to the prompt decay $\phi_{2}\rightarrow\phi_{1}X$ , these exchange processes have no role in setting the relic density for $\phi_{1}.$ Assuming that all the $\phi_{2}$ will ultimately transform themselves through decay processes to $\phi_{1}$ one can write down the Boltzmann equation for this case in terms of the total DM relic density $n=\sum{n_{\phi_{i}}}=n_{\phi_{1}}+n_{\phi_{2}}$ as Edsjo:1997bg ,

[TABLE]

where

[TABLE]

The relic density can easily be obtained approximately from the above equation as $\Omega h^{2}\approx(0.1~{}{\rm pb})/\langle\sigma v\rangle_{eff}$ Kolb:1990vq ; Bhattacharya:2016ysw . Note here that co-annihilation effect reduces with larger mass differences $\Delta m$ due to the Boltzmann suppression of $\exp(-\Delta m/T)$. In our numerical scans we will consider DM relic density to lie within: $0.1133\leq\Omega_{\rm DM}h^{2}\leq 0.1189$ Ade:2013zuv .

2.1.2 Direct Detection

In this section we will discuss about direct search constraints on multipartite $\mathbb{Z}_{2}$ model. In these experiments incoming DM flux scatter with the nuclei in the target crystals and the recoil can be searched for as a signal of the DM. Within the Higgs portal framework the direct search of the DM goes via $t-$ channel exchange of Higgs as depicted in Fig. 2. The spin independent direct search cross section of DM-Nucleon scattering reads Cline:2013gha ,

[TABLE]

where $\mu_{n}=m_{n}m_{\phi_{i}}/(m_{n}+m_{\phi_{i}})$, $m_{n}$ is the mass of the nucleon and nucleon form factor, $f_{n}\approx 0.28$ Alarcon:2011zs ; Alarcon:2012nr . As we are dealing with scalar DM and also we do not have any axial interaction term, therefore relevant bound comes from spin independent interaction cross section only. We will consider limits from the recent LUX 2016 Akerib:2016lao and XENON1T Aprile:2017iyp data from non observation of DM in direct detection experiments and compare projected sensitivity in XENONnTAprile:2015uzo and DARWIN Aalbers:2016jon experiments to validate the model. Note that unless otherwise stated, throughout this paper we use of micrOmegas Belanger:2014vza to study the spin independent direct detection cross section.

2.1.3 Numerical Scans and Analysis

The parameters of this model which govern DM phenomenology are essentially DM mass, mass of the co annihilating particle, their couplings to SM i.e $\{m_{\phi_{1}},m_{\phi_{2}},\lambda_{1h},\lambda_{2h},\lambda_{12h}\}.$ We numerically scan the parameters of the model to find the relic density allowed parameter space and then show the compatibility of the model with direct search experiments. We utilize micrOmegas Belanger:2014vza to estimate both the relic density and the direct detection spin independent cross sections as summarized in Appendix A. In the scans presented here the parameter ranges are chosen as follows,

[TABLE]

In Fig. 3 we have shown spin-independent direct detection cross section as function of DM mass ($m_{\phi_{1}}$) for the parameter space scanned. All the points in the plot satisfy the relic density constraint. The relic density allowed points are further categorized in terms of the dominant underlying number changing process that drives the freeze-out as,

• •

Co-annihilation dominant (points in red)

• •

Annihilation dominant (Cyan points)

• •

Mixed (Orange points)

If the contribution of a particular process, co-annihilation or annihilation $\geq 80\%$ to the relic density we assume that as dominant channel. For mixed cases, we choose all those points which are neither co-annihilation nor annihilation is dominant. For the sake of comparison, we also depict relic density allowed parameter space points in minimal model, i.e. with one particle under $\mathbb{Z}_{2}$ in purple. The scan also indicates the background limit from solar, atmospheric and diffuse supernovae neutrinos in gray shaded region called neutrino floor, where detection of DM signal through direct search will be difficult Billard:2013qya .

From the plot in Fig. 3, we observe that points that satisfy relic density in this model easily survives the LUX 2016 Akerib:2016lao , XENON1T Aprile:2017iyp bound and can go beyond sensitivity of DARWIN, hitting the neutrino floor. Thus, co-annihilation and mediated annihilation not only resuscitate the intermediate mass scale of $\mathbb{Z}_{2}$ Higgs portal scenario but in some regions of parameter space it remains unconstrained upto the projected limit of direct detection experiments. More interesting features arise when one investigates the underlying channels that contribute dominantly to the relic density calculations. For DM mass greater than $125~{}\rm GeV$ and below $400~{}\rm GeV$, dominant contribution to relic density and direct search allowed points come from mediated annihilation Fig. 1c. However as we increase the mass of the DM consequently the NLSP mass increases implying progressively enhanced propagator suppression and the effect fades away at larger masses. With larger DM mass, above $400~{}\rm GeV$, for surviving points the dominant contribution comes from co-annihilation. Admittedly this requires the NLSP to be relatively degenerate with $\Delta m/m_{\phi_{1}}\lesssim 20\%.$ The relic density and direct search allowed points are plotted in $m_{\phi_{1}}-m_{\phi_{2}}$ plane in Fig. 4a. The co-annihilation dominated points depicted in red predictably populate the region near $m_{\phi_{2}}\sim m_{\phi_{1}}$ due the Boltzmann factor. The first hump in the allowed parameter space for $m_{\phi_{1}}$ between $125-400$ GeV (cyan points in Fig. 4a) corresponds to mediated annihilation which fades out beyond $m_{\phi_{1}}\sim 400$ GeV. The second allowed region for heavier masses $>600$ GeV arises due to the traditional annihilation through the Higgs portal coupling of $\phi_{1}.$ The Fig. 4b in the $m_{\phi_{1}}-\lambda_{1h}$ plane clearly depicts that with large contribution from the co-annihilation and mediated annihilation the usual Higgs portal couplings can be suppressed. And as can be seen from Eq. 4 this effectively reduces the direct detection cross section.

Note that a complex scalar in $\mathbb{Z}_{2}$ framework essentially inherits two degenerate degrees of freedom, here they will amount to two degenerate DMs. Operationally this will effectively yield : $\Omega_{\rm complex}=2~{}\Omega_{\rm real}$ and the allowed parameter space can easily be scaled from above analysis. It has already been pointed out, that such two component model are relatively more constrained from direct search bounds. The allowed parameter space of this complex scalar scenario shown in Fig 5. However, in presence of co-annihilation and mediated annihilation, some of these tensions eased.

2.2 N-Scalar $\mathbb{Z}_{2}$ Model

One can easily extend this framework by populating the dark sector with more than two real scalar singlet particles transforming under same $\mathbb{Z}_{2}$ symmetry. The Lagrangian with $N$ such particles $\phi_{i}$ is given by,

[TABLE]

The lightest $\mathbb{Z}_{2}$ odd particle will be cosmologically stable and thus a DM candidate. Due to additional states in the odd sector we now have multiple copies of the co-annihilation and mediated annihilation channels assisting the freeze out of the DM. It is then easy to appreciate that the limit on each parameter gets much more relaxed compared to the two component case to survive the direct search bound after satisfying relic density. A complete analysis of the model is computationally expensive and does not introduce any novel feature in the general discussion of non-standard number changing mechanism of the DM. As an illustration of the impact of additional states we perform a simplified scan:

[TABLE]

where $\phi_{1}$ is assumed to be the DM. All other interactions between the DM states have been put to zero. Mass difference is deliberately kept low to have appreciable contribution from co-annihilation. The direct search cross-section for relic density allowed points are shown in Fig. 6 as a function of $N$ that denote the number of scalar sates in the $\mathbb{Z}_{2}$ odd sector including the DM candidate. We can see that, given the choice of parameters with relatively smaller values of $\lambda_{12h}$ in Eq. 7, the $N=2$ case can not satisfy LUX direct search bounds, whereas $N=4$ scenario can go upto the DARWIN limit 333 This should be contrasted with a multi-component DM scenario with multiple particles that are odd under different $\mathbb{Z}_{2}\times\mathbb{Z}_{2}^{\prime}\times\mathbb{Z}_{2}^{\prime\prime}\cdots$. In this case the scenario is expected to more constrained from direct detection bounds, as the number of such DM states are increased. .

3 Minimal $\mathbb{Z}_{3}$ Model and Extensions: Semi-Annihilation

A complementary approach that can accommodate processes which contribute to the DM decoupling but are not bounded by direct detection experiment is to enlarge the stabilizing symmetry $\mathbb{Z}_{2}\rightarrow\mathbb{Z}_{N}.$ In this section we will consider the scenario where a discrete $\mathbb{Z}_{3}$ symmetry stabilizes the dark sector 444For a study of $\mathbb{Z}_{4}$ models including comparison with the minimal $\mathbb{Z}_{3},$ see Belanger:2014bga .. A novel feature of this framework is the existence of the semi-annihilation processes that can potentially lead to relaxation in the direct detection bounds within Higgs portal models Belanger:2012zr ; Karam:2015jta .

A Higgs portal model with a single SM singlet complex scalar stabilized by a $\mathbb{Z}_{3}$ has been discussed in Belanger:2012zr . In this case, one requires necessarily a complex scalar field $\phi_{1}$ which transforms non-trivially under $\mathbb{Z}_{3}$ as $\phi_{1}\rightarrow\omega^{n}\phi_{1}$ where $\omega=\exp(i2\pi/3)$ and $n=1,2$ . Then invariant Lagrangian is given by,

[TABLE]

The novel features in this framework arises from the $\phi^{3}_{1}$ term proportional to the dimensionful coupling $\mu_{1}.$ This leads to DM to semi-annihilate by $\phi_{1}\phi_{1}\to\phi_{1}h$ through a s-channel and a t-channel processes shown in Fig. 7. These are the new number changing channels that are now available in addition to the usual processes shown in Fig. 1a which are common to all Higgs portal models. The semi-annihilation process is operative when the DM mass becomes heavier than Higgs mass $(m_{\phi_{1}}>m_{h})$. Note that its contribution reduces with increasing DM mass because of propagator suppression.

The Boltzmann equation for relic density in presence of semi-annihilation is given by Belanger:2012zr ,

[TABLE]

where semi-annihilation is present in addition to the annihilation cross-sections. Semi-annihilation cross-section crucially depends on the dimensionful coupling $\mu_{1}$ and also the SM-DM coupling $\lambda_{1h}$. The pronounced effect from semi-annihilation is around lower DM mass ($m_{h}\lesssim m_{\phi_{1}}\lesssim 400$) GeV where the propagator suppression is minimal. The right relic density can now be achieved for smaller values of the Higgs portal coupling ($\lambda_{1h}$) because of the assistance from the new channels. This improves the direct detection prognosis for the $\mathbb{Z}_{3}$ model.

To illustrate the essential feature of this framework we perform a three dimensional scan of the parameters as follows,

[TABLE]

where we have set the upper limit on $\mu_{1}$ from vacuum stability considerations as detaile in section 4. We follow the methodology outlined in Appendix A. The results of the scan is shown in DM mass versus direct search cross-section plane in Fig. 8. We show different choices of semi-annihilation parameter $\mu_{1}$ in different colors and it is obvious that the larger the $\mu_{1}$ is, for example, when we choose $5<\mu_{1}/m_{\phi_{1}}\leq 7.4$, we can reach the maximum sensitivity of the model by lowering spin independent cross section to nucleons while still saturating the relic limits. However, this is still not sufficient to survive the XENONnT Aprile:2015uzo except in the Higgs resonance region where the cross section can be as low as $10^{-52}~{}\rm{cm}^{2}$ as sketched previously. The apparent lower limit along the Higgs resonance branch as depicted in Fig. 8 is a result of our choices of the scan parameters.

3.1 Two particles under $\mathbb{Z}_{3}$ : Resonant semi-annihilation

We now turn to multipartite model by introducing two complex scalar $\phi_{1}$ and $\phi_{2}$ charged under the same $\mathbb{Z}_{3}$ symmetry. The Lagrangian is an generalization of Eq. 9 and can be written as,

[TABLE]

The $\mathbb{Z}_{3}$ symmetry group elements are $\{1,~{}\omega,~{}\omega^{2}\}$ and the Lagrangian above is invariant if the dark sector particles $(\phi_{1},~{}\phi_{2})$ have identical charges. If the charges are different the Lagrangian is modified, however, the essential features remains identical in both cases. Without loosing any essential physics we will assume that $m_{\phi_{1}}<m_{\phi_{2}}$ making $\phi_{1}$ as the DM candidate.

3.1.1 Relic Density

The lightest particle transforming under $\mathbb{Z}_{3}$ will be stable and serve as the DM candidate of the model, while the heavier one will have prompt decay to the DM. In this analysis, we will assume $\phi_{1}$ as the lightest stable particle and DM while $\phi_{2}$ is the next to lightest stable particle (NLSP). For further simplification we will assume $\lambda_{1h}=\lambda_{2h}$. The Boltzmann equation governing the freeze-out of this DM will be given by:

[TABLE]

where contributions from both semi-annihilation and co-annihilation dictates the thermal freeze-out of the DM on top of the annihilation cross-section of the DM ($\phi_{1}$). The late time relic density for $\phi_{1}$ depend on the number densities $n_{\phi_{1}}$ and $n_{\phi_{2}}$ at the instance of freeze out of $\phi_{1},$ since any residual $\phi_{2}$ eventually decay to $\phi_{1}.$ An estimation of this requires the simultaneous solution of the coupled Boltzmann equations for the two concerned species, however in our numerical results we consistently utilize micrOmegas. Note that the bound on the values of $\mu$’s comes from vacuum stability considerations and we adopt a conservative choice of $\mu\leq 2m_{\phi_{1}}$. Additional processes that contribute to the freeze out of the DM are shown in Fig. 9: (i) additional semi-annihilation channels $\phi_{1},\phi_{1}\to\phi_{1,2},h$ are shown in Fig. 9a, (ii) mediated annihilation shown in Fig. 9b, (iii) Co-annihilation channels shown in Fig. 9c and (iv) some novel co-annihilation channels that have numerically suppressed contribution, are shown in Fig. 9d. All these processes while crucially assisting freeze out do not contribute to the tree level DM-nucleon coupling, hence remains unconstrained from direct detection experiments.

3.1.2 Resonant Semi-Annihilation

In this section we highlight the phenomenon of resonant semi-annihilation which is possible within the multipartite $\mathbb{Z}_{3}$ scenarios. Note that for the process $\phi_{1}\phi_{1}\rightarrow\phi_{2}\rightarrow\phi_{1}h$ in Fig. 9a, has a resonance in the vicinity of $m_{\phi_{1}}\sim m_{\phi_{2}}/2$, where the semi-annihilation cross-section shoots up, reducing relic density. This is a novel feature of the non-minimal $\mathbb{Z}_{3}$ model and should be contrasted with the Higgs resonance as this does not put any restriction on the DM mass. One can tune the NLSP mass approximately to achieve this for any value of DM mass $m_{\phi_{1}}$. The couplings involved here are $\lambda_{1h},\lambda_{2h},\lambda_{12h},\mu_{1},\mu_{2},\mu_{12}$ out of which only $\lambda_{1h}$ contributes to direct search. We demonstrate this with a limited scanning in the region of interest in the parameter space as specified below,

[TABLE]

In Fig. 10a we show the variation in relic density due to change in $m_{\phi_{2}}$ for a fixed $m_{\phi_{1}}$ ($m_{\phi_{1}}=300$ GeV). The pronounced resonance effect is evident near $m_{\phi_{2}}\sim 2m_{\phi_{1}}$ where the relic density drops sharply555The features observed in the region $m_{\phi_{2}}>2m_{\phi_{1}}$ arises due to convergence issues in obtaining the thermally averaged cross section in micrOmegas and do not signify any underlying physics. Resonant semi-annihilation effect is also shown in $m_{\phi_{1}}-\Omega h^{2}$ plane for different choices of $m_{\phi_{2}}$ in Fig. 10b which clearly shows that the resonant semi-annihilation can contribute significantly for a large range of DM masses. The parameters for the scan are identical to Eq. 13 except now we vary the NLSP mass in the range $1.5~{}m_{\phi_{1}}\leq m_{\phi_{2}}\leq 2.5~{}m_{\phi_{1}}$ and $m_{\phi_{1}}$ between $125-1000$ GeV. In both cases (Fig. 10a and Fig. 10b) the couplings are deliberately chosen so that the relic density is only satisfied at resonance semi-annihilation. This implies that the Higgs portal couplings are small, easily evading direct search constraints.

3.1.3 Numerical Scans and Analysis

To investigate the generic features of the non-minimal $\mathbb{Z}_{3}$ model we perform a large five parameter scan. The relevant parameters are varied in the following range,

[TABLE]

The tree level DM-nucleon coupling relevant for direct search experiments is still mediated by the usual t-channel exchange of Higgs and is proportional to the coupling $\lambda_{1h}$ as in the minimal framework discussed in section 2.1.2. We plot the relic density allowed parameter space points emerging from the scan in the spin independent direct search cross-section versus DM mass plane in Fig. 11. We note that there exist a lot of relic density allowed points beyond XENON1T Aprile:2017iyp limit and proposed XENONnT Aprile:2015uzo limit. They can go beyond DARWIN sensitivity getting submerged into the neutrino floor for a wide range of DM mass. Once again, we highlight the dominant underlying channels by considering three regions:

• •

Mixed & Annihilation dominant (orange points),

• •

Semi-annihilation dominant (blue points),

• •

Co-annihilation dominant (red points).

By dominant, we mean that more than 80$\%$ contributions to the required annihilation cross-section appears from these channels respectively. Parameter space for which both semi-annihilation and co-annihilation is sub dominant represented by orange color and labeled as Mixed & Annihilation dominant.

To demonstrate the importance of co-annihilation and resonant semi-annihilation in simultaneously addressing the relic density prediction while avoiding stringent direct detection constraints we plot the allowed points on the $m_{\phi_{1}}-m_{\phi_{2}}$ in Fig. 12. On the top left panel Fig. 12a, we show only relic density allowed parameter space. Those satisfying both relic density constraint and XENON1T Aprile:2017iyp limit are shown in Fig. 12b. In the bottom panel 12c, we depict the parameter space that will survive even XENONnT Aprile:2015uzo limit. We again use different colors to identify points with the dominant underlying process that contribute to relic density calculations. The plots clearly shows that there are two distinctive domains in the parameter space where this model is expected to remain relatively unconstrained by the present and proposed direct detection experiments. One of them is bunched in the region $m_{\phi_{2}}\sim m_{\phi_{1}}$ and is dominated by the co-annihilation of DM with the NLSP. This also represent the only region where we can expect the multipartite $\mathbb{Z}_{2}$ models to survive. An entire new wing is obtained around $m_{\phi_{2}}\approx 2m_{\phi_{1}}$, for the resonant semi-annihilation channel. These kinematic boundary lines like $m_{\phi_{2}}=m_{\phi_{1}}$ and $m_{\phi_{2}}=2m_{\phi_{1}}$ shown in the bottom panel of Fig. 12. This is an entire new region of unconstrained parameter space that is obtained in the non-minimal $\mathbb{Z}_{3}$ model. The two regions while remaining unconstrained by direct searches should lead to distinct different consequences in the indirect searches Arcadi:2017vis . The resonant semi-annihilation region can lead to striking indirect detection signals while the co-annihilation processes do not contribute to the this signal strength at all. A detailed study of this will be carried out elsewhere.

The relic density allowed points projected in the plane of $m_{\phi_{1}}-\lambda_{1h}$ is shown in Fig. 13. One can clearly see, that most of mixed & annihilation dominated region of non-minimal $\mathbb{Z}_{3}$ is disfavored by current XENON1T limits. For large contribution from co-annihilations and semi-annihilation (red and blue points which superpose on each other), the required $\lambda_{1h}$ coupling can be brought down significantly and all these regions will be allowed by present direct search constraints.

From the experience with the $\mathbb{Z}_{2}$ case we can conclude that inclusion of more states in the dark sector charged under $\mathbb{Z}_{3},$ will add copies of the channels that are operative in this model. Proliferation of annihilation modes in the $N$-state $\mathbb{Z}_{3}$ model will contribute to further relaxation in the direct detection constraints for the relic density allowed region. It is expected that similar reduction of the direct detection cross section can now be achieved for smaller values of the co-annihilation and semi-annihilation couplings as compared to the two scalar case that has been numerically explored in this section.

3.2 Two Component DM in $\mathbb{Z}_{3}\times\mathbb{Z}_{3}^{{}^{\prime}}$ : DM Exchange

For completion we consider the two component DM that transform under two different $\mathbb{Z}_{3}$. The $\mathbb{Z}_{3}\times\mathbb{Z}_{3}^{{}^{\prime}}$ invariant Lagrangian is given by,

[TABLE]

Unlike in the previous section, here we obtain a two component DM. On top of number changing processes like annihilations and semi-annihilations, there will be exchange processes as shown in the Feynman graphs in Fig. 14. Note that co-annihilation and resonant semi-annihilation process are not present in this minimal setup. The existence of two DM candidates requires large reduction of number densities and this is more disfavored by direct search experiments. To maximize the impact of semi-annihilation we use the indulgent limit of $\mu_{i}\lesssim 7.4m_{\phi_{i}}$ Belanger:2012zr .

The coupled Boltzmann equation for this two component DM case can be written as Aoki:2012ub ,

[TABLE]

We have inserted this two component framework into micrOmegas and looked for relic density allowed parameter space using the algorithm given in Appendix. A. However, the direct detection cross section has been calculated manually using following formula Bhattacharya:2013hva ; Cao:2007fy ,

[TABLE]

where the DM-nucleon cross section in detailed in Section 2.1.2.

The key features that emerges out of the analysis is shown in Fig. 15 obtained by scanning over the following parameters,

[TABLE]

In the allowed parameter space we have scenarios where the two components are separated in mass by $>300$ GeV and remains in the range for the next generation XENON experiments. This provides the tantalizing possibility of detecting two different DM particles in direct detection experiment signaling a multi-component dark sector. This should be contrasted with the $\mathbb{Z}_{2}\times\mathbb{Z}_{2}^{\prime}$ model detailed in Bhattacharya:2016ysw , where the allowed DM pushes both components to be beyond $\gtrsim 400$ GeV. This forbids the concurrent discovery of two DM candidates in near future.

A two component DM model with $\mathbb{Z}_{3}\times\mathbb{Z}_{3}^{\prime}$ can further be extended with more dark sector particles transforming either $\mathbb{Z}_{3}$ or $\mathbb{Z}_{3}^{\prime}$ or both to have co-annihilation and resonant semi annihilation to evade direct search constraints to a great extent. However the phenomenology will be simply guided by the analogy in section 3.1 and 3.2 together and not illustrated here.

4 Brief Sketch of Vaccum Stability and Unitarity Constraints

Here we briefly summerize the constraints on the parameter space of the models considered above from tree level unitarity and vacuum stability. We will consider the $\mathbb{Z}_{2}$ and $\mathbb{Z}_{3}$ models in turn.

The stability of the minimal $\mathbb{Z}_{2}$ Higgs portal scalar DM models have been studied in Lerner:2009xg , while the extension to $\mathbb{Z}_{2}\times\mathbb{Z}^{{}^{\prime}}_{2}$ is straight-forward. For the multipartite $\mathbb{Z}_{2}$ stabilized Higgs portal models as described in the Lagrangian given in Eq. 2.1 the condition for global minimum can be obtained following Kannike:2016fmd . The potential given in Eq. 2.1 can be recast in the following form

[TABLE]

where,

[TABLE]

and

[TABLE]

Since the fields have positive nonzero mass we only consider the quartic terms in fields. Condition for the stability of the potential in Eq. 19 implies that the discriminant of it with respect to $\lvert H\rvert^{2}$ is positive. This gives rise to the following conditions,

[TABLE]

where, for simplicity, we have assumed $\lambda_{1h}=\lambda_{2h}$, $\lambda_{1s}=\lambda_{2s}=\lambda_{s}\gg\lambda_{1h}$. These are the condition for which the potential given in Eq. 19 has only one global minimum. As explained in main text $\lambda_{e_{i}},\lambda_{s}$ does not play any role in DM phenomenology, so by tuning them to be sufficiently large, but well within the perturbative limits, one can easily ensure the stability of the potential in all regions of the parameter space that has been scanned.

For tree level constraint from unitarity we utilize the Lee, Quigg and Thacker (LQT) method Lee:1977eg widely used in various BSM context Horejsi:2005da ; Das:2014fea ; Mondal:2015fja ; Chakrabarty:2014aya ; Bhattacharyya:2015nca . The scattering processes which goes through dimensionful couplings have propagator suppression remaining relatively unconstrained. The constrains on the dimension less parameters are obtained as detailed in Appendix C and is given in table 1. It is easily seen that for region of parameter space scanned the constrained are satisfied. For example in a typical choice of parameters in the region of interest: $m_{\phi_{1}}=775~{}\mbox{GeV},~{}m_{\phi_{2}}=778~{}\mbox{GeV},~{}\lambda_{1h}=\lambda_{2h}=0.002,~{}\lambda_{12h}=0.35,~{}\lambda_{e_{i}}=1~{}\mbox{and}~{}\lambda_{is}=1$ the absolute values of the solutions of Eq. 27 are $\sim 0,~{}0.1,~{}1.2~{}\mbox{and}~{}2.$

Vaccum stability for the potential having $\mathbb{Z}_{3}$ symmetry has been studied extensively in Belanger:2012zr . In this case we can have several possible stationary points. The stationary points are $(i)~{}\left\langle H\right\rangle=\left\langle\phi\right\rangle=0$, $(ii)~{}\left\langle H\right\rangle\neq 0$ and $\left\langle\phi\right\rangle=0$, $(iii)~{}\left\langle H\right\rangle=0$ and $\left\langle\phi\right\rangle\neq 0$, $(iv)~{}\left\langle H\right\rangle\neq 0$ and $\left\langle\phi\right\rangle\neq 0$. The desired SM vaccum with a stable DM is defined by $(i)$. Assuming all the quartic couplings are positive in the Lagrangian given by 1 condition for this to be the global minima sets an upper bound on the trilinear coupling $\mu_{1}\lesssim 2m_{\phi_{1}}$. A metastable vaccum with the desired property and lifetime greater than the age of the universe, relaxes the bound to $\mu_{1}\lesssim 7.4m_{\phi_{1}}$. In the $\mathbb{Z}_{3}\times\mathbb{Z}^{{}^{\prime}}_{3}$ model discussed in section 3.2 these limits gets extended on every trilinear coupling $\mu_{i}$ in the Lagrangian defined in Eq. 15. The situation for multipartite $\mathbb{Z}_{3}$ model is more involved due to the additional trilinear couplings involving multiple dark sector state. A detailed study entail the estimation of lifetime of the desired metastable vacuum which essentially constraints the trilinear couplings to be smaller and relatively large quartics. Therefore we adhere to a conservative limit of $\mu\lesssim 2m_{\phi}.$

Following the algorithm portrayed in Appendix C unitarity constraint on the parameter space for the Lagrangian stated in Eq. 3.1 are given in table 1. It is clear that the given constraints are satisfied for the region scanned. For a representative parameter point with $\lambda_{1h}=\lambda_{2h}=0.001,~{}\lambda_{12h}=0.5,~{}\lambda_{e}=1~{}\mbox{and}~{}\lambda_{is}=1,$ the solutions of Eq. 29 are $1,\,3,\,5,\,5.2$. An identical discussion also ensure unitary behavior of the $\mathbb{Z}_{3}\times\mathbb{Z}^{{}^{\prime}}_{3}$ model for the constraints given in table 1.

5 Conclusions

Ever increasing precision of direct detection experiments set considerable constraints on the otherwise well motivated scalar Higgs portal DM framework. A significant region of the parameter space within the minimal model has been ruled out by the current bounds on DM-nucleon cross section from LUX 2016 and XENON1T. Except the tuned Higgs pole region, this model is on the verge of being precluded by the next generation experiments. In this paper we perform a systematic study of simple extensions of the minimal Higgs portal framework in terms of their viability of surviving next generation direct detection experiments while saturating the DM relic abundance estimates. Here we have considered both enlarging the stabilizing symmetry from $\mathbb{Z}_{2}\rightarrow\mathbb{Z}_{3}$ and introducing multipartite dark sectors.

In the two particle $\mathbb{Z}_{2}$ framework we can have novel DM number changing processes like, co-annihilation and mediated annihilation. Phase space barrier of co-annihilation and final state Higgs in the mediated annihilation channels make these topologies insensitive to direct search experiments though they contribute to freeze out. This provides a handle to disentangle these two phenomenon which are usually correlated within the WIMP paradigm. A release of imminent tension with direct detection experiments is obtained, allowing them to survive upto and even beyond the proposed sensitivity of DARWIN. We perform an extensive numerical simulation of the model to study the interplay of these processes in alleviating the direct detection bounds. We observe that mediated annihilation processes are effective for the DM mass below $400$ GeV, while co-annihilation plays a similar role for larger masses, provided that the DM and the NLSP are relatively degenerate. Expectedly increasing the number of states contributing to these processes by populating the dark sector with more species will further alleviate the tension.

Within the minimal $\mathbb{Z}_{3}$ scalar singlet Higgs portal models the semi-annihilation processes are brought into play. This can help sustain the model beyond XENON1T for DM mass above the Higgs mass and below $\sim 400$ GeV. However, the semi-annihilation processes are constrained from above by vacuum stability considerations and will be unable to survive an absence at XENONnT. Whereas a non-minimal $\mathbb{Z}_{3}$ facilitate the possibility of co-annihilation, mediated annihilation, resonant semi-annihilation etc. which are efficient in allowing the DM direct detection cross section to as low as the sensitivity of direct detection experiments set by the ambient neutrino flux at these class of experiments. We perform a multidimensional numerical scan to present the allowed parameter space and highlight the role of various underlying processes that contribute to the relic density and direct detection calculations. We briefly discuss the possible constraints on the allowed parameter space from tree level unitarity and vacuum stability conditions.

We conclude that a no show at next generation experiments will push the intermediate scale Higgs portal DM paradigm into the multipartite era. Having multi-component DM stabilized with individual symmetries admittedly worsen the situation. However, the prognosis is much better for models where there is a multipartite dark sector all charged under the same stabilizing symmetry, where the lightest charged state is the cosmologically stable DM candidate. In this case we find that there are two interesting possibilities of co-annihilation and resonant semi-annihilation, that enable effective reduction of the DM-nucleon cross section to remain virtually unconstrained by direct detection while driving the relic density to the right ball park. Co-annihilation require a relative degeneracy in the DM and the NLSP while the possibility of the resonant semi-annihilation is effective provided we admit a tuning of the form $m_{NLSP}\sim 2m_{DM}.$ Interestingly these two possibilities while being insensitive to probing through direct detection experiments will provide strikingly different signatures for indirect detection experiments.

Acknowledgments :

We thank Ujjal Kumar Dey and Sayan Das Gupta for discussions. SB would like to acknowledge the DST-INSPIRE research grant IFA13-PH-57 at IIT Guwahati. SB and PG acknowledges the hospitality at IIT Kharagpur where the work was planned. PG and TNM would like to thank MHRD, Government of India for research fellowship. TSR acknowledges the hospitality provided by ICTP, Italy, under the Associate program, during the completion of the project. TSR is partially supported by the Department of Science and Technology, Government of India, under the Grant Agreement number IFA13-PH-74 (INSPIRE Faculty Award).

Appendix A Numerical Procedure

In this appendix we briefly summarize the numerical method we have followed to scan the parameter space through out this paper. Each model has been implemented utilizing the LanHEP Semenov:2014rea platform. The extract model files from LanHEP have been sanity checked using CalCHEP Belyaev:2012qa . Then we implement model files into micrOmegas Belanger:2014vza which have been used for making extensive scan of the models presented in the text. micrOmegas has been used to calculate both the relic density $(\Omega h^{2})$ and spin independent direct detection cross section $(\sigma_{n}^{\rm SI})$. We also extract the information for the contribution from different channel to relic density for every parameter point simulated. In the case of the two component DM presented in section 3.2, the direct detection cross section is calculated utilizing the analytical expressions given in Eq. 17, while micrOmegas has been used to obtain relic density values.

Appendix B NLSP decay

Here we point out lifetime of Next to Lightest Stable Particles which undergoes a three body decay. The decay width for the heavier particle ($\phi_{2}$) to DM ($\phi_{1}$) is given by

[TABLE]

And decay time,$\tau(\phi_{2}\to\phi_{1}\,X)=1/\Gamma(\phi_{2}\to\phi_{1}\,X)$. We would like to evaluate the limit on the parameters of the model so that the decay width is not larger than the age of the universe ($0.66\times 10^{42}$ $~{}{\rm GeV}^{-1}~{}$). A simple estimation shows that within DM mass $m_{\phi_{1}}$, $10-10^{4}$ GeV, the limit on the coupling, $\lambda_{12h}\sim 10^{-13}$ . In our analysis, we have consistently used the co-annihilation coupling $\lambda_{12h}\sim 0.1-1.0$. For this values the heavier component $\phi_{2}$ decays promptly to $\phi_{1}$.

Appendix C Unitarity Constraints

The amplitude of a scattering processes can be written in terms of Legendre polynomial as

[TABLE]

In high energy limit the $a_{0}$ partial wave will determine the leading energy dependence of the scattering processes. Therefore the unitarity constraint on the process translates to $\lvert\rm Re~{}a_{0}\rvert<1/2$. This imply the following bound on scattering amplitude Horejsi:2005da

[TABLE]

In our analysis we will confine ourself to two particle scattering processes. So the assignment is to calculate amplitude of all possible $2\to 2$ scattering processes and employ Eq. 22 on the eigenvalues of amplitude matrix.

**Two particle under $\mathbb{Z}_{2}$:

**The model defined by the Lagrangian given in Eq. 2.1 have $11$ neutral two particle states i.e.

[TABLE]

and $4$ singly charged particle states i.e.

[TABLE]

Therefore for neutral states the $11\times 11$ amplitude matrix is given by

[TABLE]

Where

[TABLE]

and

[TABLE]

For singly charged particle state the $4\times 4$ amplitude matrix is given by

[TABLE]

The distinct eigenvalues of matrix given in Eq. 23 and 24 are

[TABLE]

and the solutions of the equation given below

[TABLE]

Where,

[TABLE]

**Minimal $\mathbb{Z}_{3}$ model:

**Identically, for the Lagrangian given in Eq. 8 there will be a $11\times 11$ amplitude matrix for neutral particle states and $4\times 4$ for singly charged particle states. Corresponding distinct eigenvalues are mentioned in table 1.

**Two particle under $\mathbb{Z}_{3}$ :

**Similarly for multiparticle $\mathbb{Z}_{3}$ model we will have a $22\times 22$ amplitude matrix for neutral particle states and a $6\times 6$ matrix for the singly charged particle states. Constraints on the parameter space are given in table 1. Where $y$ is the solution of the equation given by

[TABLE]

Where

[TABLE]

$\mathbb{Z}_{3}\times\mathbb{Z}^{{}^{\prime}}_{3}$** Model :

**Constraints on dimensionless parameters for this model has been given in table 1. Note that $z$ is the solution of the following equation

[TABLE]

Where

[TABLE]