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

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.


Introduction
Existence of dark matter (DM) is supported from many astrophysical evidences like rotation curve of galaxies [1,2], anisotropies in CMBR [3] and observations in bullet cluster [4]. 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 [5][6][7][8][9].
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 Z N discrete symmetries provide the simplest realization of this stabilizing symmetry and is commonly employed in extensions of the SM [10]. 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 JHEP10(2017)088 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 [11] and then at PLANCK [12] experiments, to be in the range 0.1133 ≤ Ω DM h 2 ≤ 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 states. 1 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 [16]. However these same processes can be probed by the direct detection experiments [17,18].
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 [19] set forth a minimum upper limit on WIMP-nucleon spin independent (SI) cross section of 7.6 × 10 −46 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 [20] gathered data for 13 months between 2011 and 2012, reaching a minimum sensitivity of 2 × 10 −45 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 [21] 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×10 −48 cm 2 at WIMP mass of 50 GeV [22]. Dark matter WIMP search with liquid xenon (DARWIN) [23] 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×10 −49 cm 2 at a WIMP mass of 40 GeV [23]. 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 [9].
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 [7]. Null results at direct detection experiments [19,21] has already put a strong JHEP10(2017)088 bound [24,25] on the minimal Higgs portal models where a single scalar DM is stabilized by a discrete Z 2 symmetry. Other than the extremely tuned Higgs mass pole region, the minimal Higgs portal model has been pushed to a heavy DM mass region by XENON1T data [21], 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 Z 2 to Z N or by introducing multi-particle dark sector. These non-minimal models facilitate non-standard number changing channels like semiannihilation [26] and co-annihilation [27], which cannot be explored at direct detection experiments. We find that within multipartite Z 2 framework co-annihilation and mediated annihilation processes ameliorate some of the direct search bounds [28,29]. Further we investigate multipartite 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 Z 3 × Z 3 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 Z 2 framework. In section 3 we briefly review the impact of semi-annihilation in scalar DM models stabilized by a Z 3 symmetry, before we present a rigorous study of multiparticle Z 3 framework including a detailed discussion of resonant semi-annihilation. And in section 3.2, we focus two component DM under Z 3 × Z 3 . In section 4 we briefly discuss tree level vaccum stability and unitarity constraints. Finally we conclude in section 5.
2 Minimal Z 2 model and extensions The operators facilitating the thermal freeze-out involve DM and SM fields which can be decomposed into O SM −DM ∼ O DM 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 (φ) interacting to SM through Higgs portal interactions as φ 2 H † H and has been studied exhaustively in the literature [30]. The Lagrangian can be written as,

JHEP10(2017)088
This setup implies the presence of an unbroken Z 2 symmetry under which φ → −φ while all SM particles are even, making the field φ stable. The relevant annihilation processes that drive the freeze out are depicted in figure 1a. However, this minimal setup has been pushed to uncomfortable corner by recent results from direct search data in LUX 2016 [31] and XENON1T [21]. The remaining unconstrained region are now confined to the tuned Higgs resonance region or for heavy DM mass m φ 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 φ ≈ m h /2 ∼ 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 σv ≈ 24 λ 2 h GeV −2 assuming a Higgs width Γ h = 4 MeV [32]. This saturates the relic density bound when σv is equal to ∼ 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 × 10 −52 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 Z 2 × Z 2 · · · 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 Z 2 × Z 2 has been studied in [33,34]. 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 [35].
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 nonstandard 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.

Two particles under Z 2 : co-annihilation & mediated annihilation
We will consider the case where the SM is augmented by two real scalar particles φ 1 and φ 2 odd under the same discrete Z 2 symmetry and singlet under SM. The relevant part of JHEP10(2017)088 the Lagrangian is given by, Note that the coupling λ 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 figure 1b and figure 1c respectively. Assuming without any loss of generality, m φ 1 < m φ 2 the lightest mass state φ 1 can be identified as the potential DM candidate. 2 The φ 2 state can promptly decay to φ 1 , therefore do not effect the freeze out except through its contribution to the φ 1 number changing process discussed above. This decay occurs through an off-shell Higgs and a schematic calculation can be found in appendix B.

Relic density
The new set of Feynman graphs corresponding to the number changing processes of DM are shown in figure 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 figure 1b if the masses of the two sates are relatively degenerate. However, for large mass gap the t-channel mediated annihilation process in figure 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 λ e i in eq. (2.2) essentially influences the evolution of the DM density in two ways: (i) From the decay of φ 2 , DM is produced and (ii) DM scattering processes: However, due to the prompt decay φ 2 → φ 1 X , these exchange processes have no role in setting the relic density for φ 1 . Assuming that all the φ 2 will ultimately transform themselves through decay processes to φ 1 one can write down the Boltzmann equation for this case in terms of the total DM relic density n = n φ i = n φ 1 + n φ 2 as [27], In principle we can have a term like m 2 φ1φ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.2) is written in the mass basis of φ1 and φ2, assuming the mass matrix has been diagonalized. The relic density can easily be obtained approximately from the above equation as Ωh 2 ≈ (0.1 pb)/ σv eff [33,36]. Note here that co-annihilation effect reduces with larger mass differences ∆m due to the Boltzmann suppression of exp(−∆m/T ). In our numerical scans we will consider DM relic density to lie within: 0.1133 ≤ Ω DM h 2 ≤ 0.1189 [12].

Direct detection
In this section we will discuss about direct search constraints on multipartite 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 figure 2. The spin independent direct search cross section of DM-Nucleon scattering reads [37], where µ n = m n m φ i /(m n + m φ i ), m n is the mass of the nucleon and nucleon form factor, f n ≈ 0.28 [38,39]. 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 [31] and XENON1T [21] data from non observation of DM in direct detection experiments and compare projected sensitivity in XENONnT [22] and DARWIN [23] experiments to validate the model. Note  that unless otherwise stated, throughout this paper we use of micrOmegas [40] to study the spin independent direct detection cross section.

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 φ 1 , m φ 2 , λ 1h , λ 2h , λ 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 [40] 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, In figure 3 we have shown spin-independent direct detection cross section as function of DM mass (m φ 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 ≥ 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 Z 2 in purple. The scan also indicates the background limit from JHEP10(2017)088     solar, atmospheric and diffuse supernovae neutrinos in gray shaded region called neutrino floor, where detection of DM signal through direct search will be difficult [41].
From the plot in figure 3, we observe that points that satisfy relic density in this model easily survives the LUX 2016 [31], XENON1T [21] 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 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 GeV and below 400 GeV, dominant contribution to relic density and direct search allowed points come from mediated annihilation figure 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 GeV, for surviving points the dominant contribution comes from co-annihilation. Admittedly this requires the NLSP to be relatively degenerate with ∆m/m φ 1 20%. The relic density and direct search allowed points are plotted in m φ 1 − m φ 2 plane in figure 4a. The co-annihilation dominated points depicted in red predictably populate the region near m φ 2 ∼ m φ 1 due the Boltzmann factor. The first hump in the allowed parameter space for m φ 1 between 125−400 GeV (cyan points in figure 4a) corresponds to mediated annihilation which fades out beyond m φ 1 ∼ 400 GeV. The second allowed region for heavier masses > 600 GeV arises due to the traditional annihilation through the Higgs portal coupling of φ 1 . The figure 4b in the m φ 1 −λ 1h plane clearly depicts that with large contribution from the coannihilation and mediated annihilation the usual Higgs portal couplings can be suppressed. And as can be seen from eq. (2.4) this effectively reduces the direct detection cross section.
Note that a complex scalar in Z 2 framework essentially inherits two degenerate degrees of freedom, here they will amount to two degenerate DMs. Operationally this will effectively JHEP10(2017)088 yield: Ω complex = 2 Ω 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 figure 5. However, in presence of co-annihilation and mediated annihilation, some of these tensions eased.

N-scalar 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 Z 2 symmetry. The Lagrangian with N such particles φ i is given by, The lightest 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 JHEP10(2017)088 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: where φ 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 figure 6 as a function of N that denote the number of scalar sates in the Z 2 odd sector including the DM candidate. We can see that, given the choice of parameters with relatively smaller values of λ 12h in eq. (2.7), the N = 2 case can not satisfy LUX direct search bounds, whereas N = 4 scenario can go upto the DARWIN limit. 3 3 Minimal 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 Z 2 → Z N . In this section we will consider the scenario where a discrete Z 3 symmetry stabilizes the dark sector. 4 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 [43,44]. A Higgs portal model with a single SM singlet complex scalar stabilized by a Z 3 has been discussed in [43]. In this case, one requires necessarily a complex scalar field φ 1 which transforms non-trivially under Z 3 as φ 1 → ω n φ 1 where ω = exp(i2π/3) and n = 1, 2 . Then invariant Lagrangian is given by, The novel features in this framework arises from the φ 3 1 term proportional to the dimensionful coupling µ 1 . This leads to DM to semi-annihilate by φ 1 φ 1 → φ 1 h through a s-channel and a t-channel processes shown in figure 7. These are the new number changing channels that are now available in addition to the usual processes shown in figure 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 φ 1 > m h ). Note that its contribution reduces with increasing DM mass because of propagator suppression. Figure 7. Feynman diagrams for Semi-annihilation in Z 3 model as in eq. (3.1). The Boltzmann equation for relic density in presence of semi-annihilation is given by [43],

JHEP10(2017)088
where semi-annihilation is present in addition to the annihilation cross-sections. Semiannihilation cross-section crucially depends on the dimensionful coupling µ 1 and also the SM-DM coupling λ 1h . The pronounced effect from semi-annihilation is around lower DM mass (m h m φ 1 400) GeV where the propagator suppression is minimal. The right relic density can now be achieved for smaller values of the Higgs portal coupling (λ 1h ) because of the assistance from the new channels. This improves the direct detection prognosis for the Z 3 model.

JHEP10(2017)088
To illustrate the essential feature of this framework we perform a three dimensional scan of the parameters as follows, where we have set the upper limit on µ 1 from vacuum stability considerations as detailed 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 figure 8. We show different choices of semi-annihilation parameter µ 1 in different colors and it is obvious that the larger the µ 1 is, for example, when we choose 5 < µ 1 /m φ 1 ≤ 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 [22] except in the Higgs resonance region where the cross section can be as low as 10 −52 cm 2 as sketched previously. The apparent lower limit along the Higgs resonance branch as depicted in figure 8 is a result of our choices of the scan parameters.

Two particles under Z 3 : resonant semi-annihilation
We now turn to multipartite model by introducing two complex scalar φ 1 and φ 2 charged under the same Z 3 symmetry. The Lagrangian is an generalization of eq. (3.2) and can be written as, The Z 3 symmetry group elements are {1, ω, ω 2 } and the Lagrangian above is invariant if the dark sector particles (φ 1 , φ 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 φ 1 < m φ 2 making φ 1 as the DM candidate.

Relic density
The lightest particle transforming under 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 φ 1 as the lightest stable particle and DM while φ 2 is the next to lightest stable particle (NLSP). For further simplification we will assume λ 1h = λ 2h . The Boltzmann JHEP10(2017)088 equation governing the freeze-out of this DM will be given by: 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 (φ 1 ). The late time relic density for φ 1 depend on the number densities n φ 1 and n φ 2 at the instance of freeze out of φ 1 , since any residual φ 2 eventually decay to φ 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 µ's comes from vacuum stability considerations and we adopt a conservative  choice of µ ≤ 2m φ 1 . Additional processes that contribute to the freeze out of the DM are shown in figure 9: (i) additional semi-annihilation channels φ 1 , φ 1 → φ 1,2 , h are shown in figure 9a, (ii) mediated annihilation shown in figure 9b, (iii) Co-annihilation channels shown in figure 9c and (iv) some novel co-annihilation channels that have numerically suppressed contribution, are shown in figure 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.

Resonant semi-annihilation
In this section we highlight the phenomenon of resonant semi-annihilation which is possible within the multipartite Z 3 scenarios. Note that for the process φ 1 φ 1 → φ 2 → φ 1 h in figure 9a, has a resonance in the vicinity of m φ 1 ∼ m φ 2 /2, where the semi-annihilation cross-section shoots up, reducing relic density. This is a novel feature of the non-minimal 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 φ 1 . The couplings involved here are λ 1h , λ 2h , λ 12h , µ 1 , µ 2 , µ 12 out of which only λ 1h contributes to direct search. We demonstrate this with a limited scanning in the region of interest in the parameter space as specified below, In figure 10a we show the variation in relic density due to change in m φ 2 for a fixed m φ 1 (m φ 1 = 300 GeV). The pronounced resonance effect is evident near m φ 2 ∼ 2m φ 1 where the relic density drops sharply. 5 Resonant semi-annihilation effect is also shown in m φ 1 − Ωh 2 5 The features observed in the region m φ 2 > 2m φ 1 arises due to convergence issues in obtaining the thermally averaged cross section in micrOmegas and do not signify any underlying physics.  figure 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.

Numerical scans and analysis
To investigate the generic features of the non-minimal Z 3 model we perform a large five parameter scan. The relevant parameters are varied in the following range, 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 λ 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 crosssection versus DM mass plane in figure 11. We note that there exist a lot of relic density allowed points beyond XENON1T [21] limit and proposed XENONnT [22] limit. They can JHEP10(2017)088 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 φ 1 − m φ 2 in figure 12. On the top JHEP10(2017)088 left panel figure 12a, we show only relic density allowed parameter space. Those satisfying both relic density constraint and XENON1T [21] limit are shown in figure 12b. In the bottom panel 12c, we depict the parameter space that will survive even XENONnT [22] 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 φ 2 ∼ m φ 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 Z 2 models to survive. An entire new wing is obtained around m φ 2 ≈ 2m φ 1 , for the resonant semi-annihilation channel. These kinematic boundary lines like m φ 2 = m φ 1 and m φ 2 = 2m φ 1 shown in the bottom panel of figure 12. This is an entire new region of unconstrained parameter space that is obtained in the non-minimal Z 3 model. The two regions while remaining unconstrained by direct searches should lead to distinct different consequences in the indirect searches [45]. 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 φ 1 −λ 1h is shown in figure 13. One can clearly see, that most of mixed & annihilation dominated region of non-minimal 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 λ 1h coupling can be brought down significantly and all these regions will be allowed by present direct search constraints. From the experience with the Z 2 case we can conclude that inclusion of more states in the dark sector charged under Z 3 , will add copies of the channels that are operative in this model. Proliferation of annihilation modes in the N -state 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.

Two component DM in Z 3 × Z 3 : DM exchange
For completion we consider the two component DM that transform under two different Z 3 . The Z 3 × Z 3 invariant Lagrangian is given by, 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 figure 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 µ i 7.4m φ i [43].

JHEP10(2017)088
The coupled Boltzmann equation for this two component DM case can be written as [46], 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 [47,48], 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 figure 15 obtained by scanning over the following parameters, 125 ≤ m 1 ≤ 500, 500 ≤ m 2 ≤ 1000, 5m i ≤ µ 1 = µ 2 ≤ 7.4m i , 0.01 ≤ λ 1h = λ 2h ≤ 0.1, 0.5 ≤ λ e ≤ 1.5 (3.11) 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 Z 2 × Z 2 model detailed in [33], where the allowed DM pushes both components to be beyond 400 GeV. This forbids the concurrent discovery of two DM candidates in near future.
A two component DM model with Z 3 × Z 3 can further be extended with more dark sector particles transforming either Z 3 or Z 3 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.

Brief sketch of vaccum stability and unitarity constraints
Here we briefly summarize the constraints on the parameter space of the models considered above from tree level unitarity and vacuum stability. We will consider the Z 2 and Z 3 models in turn. Figure 15. Spin-independent DM-nucleon effective cross-section vs DM mass for relic density allowed points in Z 3 × Z 3 model indicated in eq. (3.8) where both φ 1 and φ 2 are allowed by XENON1T limit for the scan over parameters specified in eq. (3.11).

JHEP10(2017)088
The stability of the minimal Z 2 Higgs portal scalar DM models have been studied in [30], while the extension to Z 2 ×Z 2 is straight-forward. For the multipartite Z 2 stabilized Higgs portal models as described in the Lagrangian given in eq. (2.2) the condition for global minimum can be obtained following [49]. The potential given in eq. (2.2) can be recast in the following form where, Since the fields have positive nonzero mass we only consider the quartic terms in fields. Condition for the stability of the potential in eq. (4.1) implies that the discriminant of it with respect to |H| 2 is positive. This gives rise following conditions on the parameter where, for simplicity, we have assumed λ 1h = λ 2h , λ 1s = λ 2s = λ s λ 1h . These are the condition for which potential given in eq. (4.1) has only one global minimum. As explained in main text λ e i , λ 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 [50] widely used in various BSM context [51][52][53][54][55]. 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 JHEP10(2017)088

Model
Lagrangian Tree level unitarity constraints x < 8π where x is the solution of the equation (C.6).
z < 8π where z is the solution of the equation (C.10).  Vaccum stability for the potential having Z 3 symmetry has been studied extensively in [43]. In this case we can have several possible stationary points. The stationary points are (i) H = φ = 0, (ii) H = 0 and φ = 0, (iii) H = 0 and φ = 0, (iv) H = 0 and φ = 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 (2.1) condition for this to be the global minima sets an upper bound on the trilinear coupling µ 1 2m φ 1 . A metastable vaccum with the desired property and lifetime greater than the age of the universe, relaxes the bound to µ 1 7.4m φ 1 . In the Z 3 × Z 3 model discussed in section 3.2 these limits gets extended on every trilinear coupling µ i in the Lagrangian defined in eq. (3.8). The situation for multipartite 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 µ 2m φ .
Following the algorithm portrayed in appendix C unitarity constraint on the parameter space for the Lagrangian stated in eq. (3.4) are given in table 1. It is clear that the given constraints are satisfied for the region scanned. For a representative parameter point with λ 1h = λ 2h = 0.001, λ 12h = 0.5, λ e = 1 and λ is = 1, the solutions of eq. (C.8) are 1, 3, 5, 5.2. An identical discussion also ensure unitary behavior of the Z 3 × Z 3 model for the constraints given in table 1.

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 JHEP10(2017)088 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 Z 2 → Z 3 and introducing multipartite dark sectors.
In the two particle 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 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 ∼ 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 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 DMnucleon 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 N LSP ∼ 2m DM . Interestingly these two possibil- ities while being insensitive to probing through direct detection experiments will provide strikingly different signatures for indirect detection experiments.

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 (φ 2 ) to DM (φ 1 ) is given by And decay time,τ (φ 2 → φ 1 X) = 1/Γ(φ 2 → φ 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 × 10 42 GeV −1 ). A simple estimation shows that within DM mass m φ 1 , 10 − 10 4 GeV, the limit on the coupling, λ 12h ∼ 10 −13 . In our analysis, we have consistently used the co-annihilation coupling λ 12h ∼ 0.1 − 1.0. For this values the heavier component φ 2 decays promptly to φ 1 .

C Unitarity constraints
The amplitude of a scattering processes can be written in terms of Legendre polynomial as M(θ) = 16π l=∞ l=0 a l (2l + 1)P l (cos θ) 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 |Re a 0 | < 1/2. This imply the following bound on scattering amplitude [51] |M| < 8π (C.1) In our analysis we will confine ourself to two particle scattering processes. So the assignment is to calculate amplitude of all possible 2 → 2 scattering processes and employ eq. (C.1) on the eigenvalues of amplitude matrix.
Two particle under Z 2 . The model defined by the Lagrangian given in eq. (2.2) have 11 neutral two particle states i.e.

JHEP10(2017)088
Minimal Z 3 model. Identically, for the Lagrangian given in eq. (3.1) there will be a 11 × 11 amplitude matrix for neutral particle states and 4 × 4 for singly charged particle states. Corresponding distinct eigenvalues are mentioned in table 1.