Higgs Portals for Thermal Dark Matter - EFT Perspectives and the NMSSM -

We analyze a low energy effective model of Dark Matter in which the thermal relic density is provided by a singlet Majorana fermion which interacts with the Higgs fields via higher dimensional operators. Direct detection signatures may be reduced if blind spot solutions exist, which naturally appear in models with extended Higgs sectors. Explicit mass terms for the Majorana fermion can be forbidden by a $Z_3$ symmetry, which in addition leads to a reduction of the number of higher dimensional operators. Moreover, a weak scale mass for the Majorana fermion is naturally obtained from the vacuum expectation value of a scalar singlet field. The proper relic density may be obtained by the $s$-channel interchange of Higgs and gauge bosons, with the longitudinal mode of the $Z$ boson (the neutral Goldstone mode) playing a relevant role in the annihilation process. This model shares many properties with the Next-to-Minimal Supersymmetric extension of the Standard Model (NMSSM) with light singlinos and heavy scalar and gauge superpartners. In order to test the validity of the low energy effective field theory, we compare its predictions with those of the ultraviolet complete NMSSM. Extending our framework to include $Z_3$ neutral Majorana fermions, analogous to the bino in the NMSSM, we find the appearance of a new bino-singlino well tempered Dark Matter region.


Introduction
While the Standard Model (SM) is extremely successful in describing the known particle interactions, it fails to explain the large scale structure of the Universe, since it does not provide a good Dark Matter (DM) candidate. The simplest addition to the SM particle content would be in the form of a SM gauge singlet and in this work we shall concentrate on the particular example of a fermion as our DM candidate. In recent years, DM-nucleon scattering experiments such as LUX, XENON1T and PandaX-II have set stringent bounds on the possible couplings of DM to SM particles [1][2][3][4][5][6]. In particular, the coupling of the 125 GeV Higgs boson to DM is significantly constrained. In addition, the vector coupling of DM to the Z gauge boson must be very small (see for instance Ref. [7]), therefore we shall concentrate on singlet Majorana fermions, which couple only axially to the Z boson. Such a fermion has the same gauge quantum numbers as a right-handed neutrino. One can define a matter parity, based on the (B −L) quantum numbers of particles, namely P = (−1) 3(B−L) , and demand interactions to be invariant under such a parity. Assuming that the DM carries no baryon B or lepton L numbers, this forbids all renormalizable interactions of the DM with SM particles, while allowing all SM Yukawa terms and Majorana masses for the right-handed neutrinos.
Since the coupling of DM to SM mediators is strongly constrained, we shall consider extending the SM by additional particles which can mediate interactions between DM and SM particles. Experimental precision tests of the SM strongly constrain extensions of the SM gauge sector, while far less is known about the SM Higgs sector. A well studied extension of the SM Higgs sector are so-called two Higgs doublet models (2HDMs) [8], which consist of adding a second Higgs doublet, as commonly found in models that provide a dynamical origin of the electroweak symmetry breaking (EWSB) mechanism. The interactions of DM may quite generally be described by a set of non-renormalizable operators, including Majorana fermion bilinears and SM gauge invariant operators. The lower dimensional operators involve interactions with the Higgs fields and constitute a simple generalization of the socalled Higgs portal models [9][10][11][12]. As we shall see, the extended Higgs sector allows for the existence of blind spots where the interaction of the Higgs bosons with DM particles may be reduced, satisfying direct detection constraints, while still allowing for the possibility of obtaining the observed (thermal) relic density [13][14][15][16][17][18][19].
We shall require DM to be a weakly interacting massive particle (WIMP). In order to obtain a weak scale mass for the Majorana fermion in a natural way, we demand it to proceed from the vacuum expectation value (vev) of a singlet scalar field, which develops in the process of EWSB. The absence of explicit masses may be the result of the presence of an explicit Z 3 symmetry, which also reduces the allowed number of higher dimensional operators and leads to a redefinition of the blind spot condition.
A possible realization of these class of models is provided by supersymmetric (SUSY) extensions of the SM [20] which also allow for a dynamical explanation for the weak scale [20][21][22]. A particular virtue of the SUSY framework is that the stability of the Higgs mass parameter under quantum corrections can be ensured. In minimal extensions, the SM-like Higgs boson is naturally light [23][24][25], and corrections to electroweak precision and flavor observables tend to be small, leading to good agreement with observations. Additionally, low scale SUSY leads to the unification of couplings at high energies and provides a natural DM candidate, namely the lightest neutralino.
Among the simplest SUSY extensions, the Next-to-Minimal Supersymmetric extension of the SM (NMSSM) [26], fulfills all of the above properties while additionally containing a rich Higgs and neutralino spectrum. This may have an important impact on low energy observables. In particular, if the lighter neutralinos and neutral Higgs bosons are mainly singlets, they would be predominantly produced in association with heavier Higgs bosons or from the cascade decays of other SUSY particles, and therefore can easily avoid current direct experimental constraints [27][28][29][30][31][32][33].
The light neutralino in the NMSSM is naturally mostly singlino-like, but with a nonnegligible Higgsino component. Hence, its spin independent direct DM detection (SIDD) cross section is mediated predominantly via the SM-like Higgs boson. The current bounds on the SIDD cross section lead to relevant constraints on the couplings of DM to the SMlike Higgs boson, and demand the theory to be in the proximity of blind spots, where the contributions from the non-standard Higgs bosons become also relevant. The proper relic density may also be obtained; the thermal annihilation cross section is dominated by either resonant contributions of the Higgs bosons, or non-resonant Z boson exchange contributions, with subdominant contributions from the light CP-even or the CP-odd Higgs bosons, the latter also having a large singlet component. In addition, a bino-like neutralino region with non-negligible Higgsino component may be present. In such a case, a sufficiently large thermal annihilation cross section yielding the proper relic density can be obtained by co-annihilation with the next-to-lightest neutralino, generally the singlino.
In Section 2 we use the language of Effective Field Theories (EFT) to outline the generic requirements of a model with singlet Majorana fermion DM and identify the required extended Higgs sector. In particular, we show the correlations of EFT parameters necessary to simultaneously obtain a thermal relic density, satisfy SIDD constraints, and accommodate a phenomenologically consistent Higgs sector. In Section 3 we discuss the NMSSM as a possible ultraviolet completion of our EFT model, and demonstrate the mapping of EFT parameters to NMSSM parameters utilizing a top-down EFT approach. In Section 4 we use the mature numerical tools available for the NMSSM to study the DM phenomenology, taking into account current collider and astrophysical constraints, as well as projections for the future. We identify two viable regions of parameter space with different DM phenomenology: 1) a new well tempered DM region, where the DM candidate is mostly bino-like and thermal production proceeds via resonant annihilation or co-annihilation with the singlino-like state, and 2) the region where the DM candidate is mostly singlino-like and the thermal relic density is mainly achieved via interactions mediated by the longitudinal mode of the Z boson, the neutral Goldstone mode. Much of the phenomenology in both regions can be understood from the properties of the EFT worked out in Section 2, although some details are only found in complete models such as the NMSSM. We reserve Section 5 for our conclusions.

An EFT for Singlet Dark Matter
As motivated in the introduction, we will consider a model of SM singlet Majorana fermion DM, which has no renormalizable interactions with SM particles. In order to couple DM to the SM, we consider a 2HDM Higgs sector, more specifically, we shall take two Higgs doublets with opposite hypercharges Y , H u with Y = +1/2 and H d with Y = −1/2, which are naturally responsible for generating masses for the up and down-type quarks, respectively, as in type II 2HDMs. Since non-renormalizable interactions are suppressed by the scale associated with the masses of a heavy sector that was integrated out, one expects the dominant interactions to be associated with lower dimensional operators.
Including operators of dimension d ≤ 5, the generic Lagrangian density describing interactions of a Majorana fermion χ with the two Higgs doublets H u , H d is where we have imposed a symmetry H d ↔ H u and used a dot notation for SU (2) products, Assuming, as usual, that both Higgs doublets acquire vevs,

5)
where H 0 d and H 0 u denote the neutral components of the respective Higgs doublets. These may be related to the usual type II 2HDM Higgs bosons by the relations The H SM interaction eigenstate has the same couplings to SM particles as a SM Higgs boson, G 0 is the (neutral) Goldstone mode making up the longitudinal polarization of the Z boson after EWSB, and H NSM and A NSM are the non SM-like CP-even and CP-odd states, respectively. In particular, note that the Higgs basis fields are defined such that all the SM vev is acquired by the field corresponding to the neutral component of H SM , hence H SM = √ 2v and H NSM = 0. Since the observed 125 GeV Higgs state h appears to be close to SM-like in nature [41,42], the interactions of χ with h may be obtained from the above, approximating h ∼ H SM to first order. Ignoring the charged Higgs fluctuations, we obtain, at linear order in the fields, Hence, the SM-like Higgs coupling to DM becomes The interaction of a Majorana fermion with the SM-like Higgs boson listed above may be suppressed in three scenarios: 1) suppression of the couplings δ and γ; 2) large values of µ v; and 3) a particular correlation of the two couplings δ and γ resulting in g χχh ∼ 0. The last scenario, the so-called blind spot solution, is given by It is interesting to consider a model in which there are no explicit mass terms or scales and hence the Lagrangian is scale invariant. In such a situation, a natural way to generate the mass m χ and the scale µ is via the vev of a singlet S = S + 1 √ 2 H S + iA S . Hence, without loss of generality we can define m χ = 2κ S and µ = λ S , where κ and λ are dimensionless parameters.
The absence of explicit scale dependence could be understood as originating from a Z 3 symmetry, under which all scalar and fermion fields transform like Ψ → exp[2πi/3] Ψ (therefore also µ → exp[2πi/3] µ). Besides forbidding explicit fermion mass terms, imposing such a Z 3 symmetry also forbids certain interactions. The remaining d ≤ 5 terms are resulting in the following DM-Higgs sector interactions: (2.12) Imposing the Z 3 symmetry removes the possibility of a blind spot as defined in Eq. (2.10). The contributions from the χχS coupling to either the thermal annihilation cross section relevant for the relic density or the SIDD cross section is further suppressed by singlet-doublet mixing since the singlet S does not couple to SM particles beyond the Higgs sector. Hence, the dominant contributions to the SIDD and the thermal annihilation cross section will be proportional to δ 2 . Barring accidental cancellations between contributions from different Higgs bosons, the coupling δ must be suppressed in order to satisfy the stringent bounds from direct detection experiments. Hence, since current data implies that the dimension d = 5 operators must be suppressed, we will include d = 6 operators in the following. As we shall demonstrate, this will again allow for blind spot solutions to appear, enabling the suppression of the SIDD cross section. In addition, we find relevant contributions to the annihilation cross section from d = 6 operators which will allow us to obtain sufficiently large annihilation cross sections to avoid over-closure of the Universe. The most relevant d = 6 operators are suppressed by powers of m χ /µ with respect to the d = 5 ones, and thus become most relevant if the ratio m χ /µ is not very small. One could inquire about the impact of higher dimensional d > 6 operators in such a case. We shall address this question later by considering an ultraviolet completion of the EFT. Although the qualitative features found in the EFT remain valid in the complete theory, the precise quantitative predictions will indeed be affected to some degree by higher dimensional terms.
Assuming that the d > 4 terms in Eq. (2.11) originate from a theory where a heavier SU (2)-doublet Dirac fermion with mass µ has been integrated out, we can write all the allowed d = 6 operators which would arise from integrating out such a field. Ignoring the charged gauge boson interactions, we get where S = µ/λ +Ŝ, Q and T 3 are the charge and weak isospin operators, s W ≡ sin θ W with the weak mixing angle θ W , and g 1 = e/ cos θ W is the hypercharge coupling. Note, that the term proportional toŜ (the fluctuations of S) in the δ-term arises because this originally d = 5 term was actually suppressed by 1/λS, which we have expanded around the vev of S, On the other hand, all the d = 6 terms arising from integrating out a Dirac fermion are suppressed by 1/λ 2 . Moreover, we have not included terms involving higher powers of the singlet field, since they are not expected to arise from integrating out a Dirac doublet fermion. The DM interactions with singlets are dominated by the tree level coupling κ, and the only modification from such terms would be a redefinition of the Sχχ coupling κ → κ(1 + O(m χ /µ). Observe, that if dealing with on-shell χ fields, there is a redundancy in the above terms, since the application of the equation of motion on the terms proportional to the derivative of χ will lead to terms proportional to the χ mass, which also appear from the κ-term when inserting the vev of the field S.
where the dependence on αm χ results from the application of the equations of motion.
In general, we calculate the on-shell relationships by using the fact that, ignoring total derivatives, where Φ is a real scalar field. Note, that the direct expansion of the derivative terms proportional to α in Eq. (2.13) leads to interactions with the CP-even Higgs bosons when the derivative is acting on the Majorana fermion fields, and to derivative interactions with the CP-odd Higgs states when the derivative is acting on the Higgs doublets, as required by hermiticity. We see that the blind spot for the cancellation of the coupling of H SM to pairs of DM now occurs for The χ interactions with H NSM are Note, that there are no terms proportional to γ (or m χ ) and therefore there is no blindspot such as the one for H SM in Eq. (2.16); instead g χχH NSM → 0 for tan β → 1.
On the other hand, the interactions with the CP-even singlet state are given by Here, the dependence on α comes from a field renormalization of χ necessary to retain a canonical kinetic term for χ when including dimension d ≤ 6 operators. In principle, this field renormalization introduces corrections to all couplings of χ. However, we are only considering operators of d ≤ 6. The modification from the field renormalization is suppressed by |µ| −2 , hence, this correction is only relevant for the renormalizable χχŜ interactions.
The interactions of χ with the CP-odd scalars are easy to read from the above as well. For instance, although the Goldstone interactions involve derivatives of the Goldstone fields, for on-shell χ's one can use Eq. (2.15) to obtain the interaction with the (neutral) Goldstone mode The orthogonal state, A NSM , also has relevant interactions with DM, namely Finally, the interactions of the CP-odd singlet state A S are analogous to its CP-even counterpart, (2.22)

Higgs Sector
In the previous section we have motivated a structure for the scalar sector consisting of two Higgs doublets and one singlet, all three of which acquire a vev. We can define the extended Higgs Basis, {H SM , H NSM , H S } for the CP-even states and {A NSM , A S } for the CP-odd states, where the doublet components are as defined in Eqs. (2.3)-(2.6), and the singlet S = S +Ŝ = µ/λ+ 1 √ 2 H S + iA S does not get rotated [29]. These interaction eigenstates mix into mass eigenstates. We denote the CP-even mass eigenstates as h i = {h, H, h S }, and the CP-odd states as a i = {A, a S }, The mixing angles S j i and P j i are obtained from the diagonalization of the corresponding mass matrices. We can write the (symmetric) squared mass matrix of the CP-even Higgs for the doublet-like eigenstate, and for the singlet-like mass eigenstate. If we use the approximate eigenmasses, we find for the SM-like mass eigenstate and for the other mass eigenstates

EFT: Relic Density
In the absence of co-annihilation, the thermally averaged annihilation cross section for a pair of DM particles at temperature T can be expanded as where x = m χ /T . After integrating over the thermal history of the Universe until the freeze-out temperature T F = m χ /x F , the thermal relic density is obtained The interactions of the singlet fermion χ with SM particles depicted in Fig. 1 arise via the couplings to the extended Higgs basis states given in Eqs.(2.14)-(2.22) and the mixing of extended Higgs basis states into mass eigenstates, (2.40) The singlet states H S and A S do not couple to SM particles, thus, assuming a type II 2HDM Yukawa structure, the couplings of the mass eigenstates to up-type quarks are given by and to down-type quarks by For completeness, we record the couplings to pairs of vector bosons The contribution to σv x F from annihilation into pairs of quarks (χχ → qq) from the s-channel exchange of the CP-even Higgs bosons is given by and from the exchange of CP-odd Higgs boson by the width of the Higgs mass eigenstate Φ i . Note, that there is no interference between the contributions listed in Eqs. (2.44) and (2.45) since the scalar Higgs bosons exchange contribution is p-wave suppressed while the annihilation cross section via pseudoscalar Higgs bosons is s-wave. For typical freeze-out temperatures m χ /T F 25, the contribution from CP-even Higgs bosons to σv x F is suppressed by 3T F /4m χ ∼ 1/30 compared to the contribution from CP-odd Higgs bosons, as long as m q /m χ 1 such that the kinematic correction from the quark mass is irrelevant. Besides via Higgs bosons, (χχ → qq) annihilations can also be mediated by the schannel exchange of Z bosons. This is accounted for by extending the sum in Eq. (2.45) to include the s-wave amplitudes mediated by both the longitudinal polarization of the Z boson, i.e. the (neutral) Goldstone mode G 0 , as well as the transversal polarizations of the Z boson The couplings of the Goldstone mode to up-type and down-type quarks, respectively, are 48) and the relevant axial-vector coupling of the transversal polarizations of the Z boson to quarks are The coupling of the Z boson to the Majorana fermion can be read off from Eqs. (2.13) or (A.1) Note, that the s-wave contribution to the annihilation cross section from the transversal polarization of the Z boson is suppressed with respect to that of its longitudinal polarization (the neutral Goldstone mode) by All the amplitudes appearing in Eqs. (2.44)-(2.47) are proportional to the Yukawa couplings. Due to the hierarchy of the Yukawa couplings, the contribution to the thermal cross section from (χχ → qq) annihilations will be dominated by the heaviest accessible quarks, i.e. top-quarks for m χ > m t ∼ 173 GeV and bottom quarks for lighter m χ . In the latter case the p-wave contribution from the transversal polarization of the Z may become relevant, since in contrast to its s-wave contribution listed above it is not chirality suppressed (hence, not proportional to the Yukawa couplings).
It is interesting to consider the size of the thermally averaged cross section obtainable via (χχ → qq) annihilation. For example, if we assume m χ > m t , such that (χχ → tt) is kinematically allowed, and assume the dominant annihilation channel to be via the (neutral) Goldstone mode, for m 2 (2.51) Hence, the correct relic density Ωh 2 ∼ 0.12 [43] is obtained from (χχ → qq) annihilations for couplings |g χχG 0 | ∼ 0.1. In addition to the (χχ → qq) annihilation discussed above, there may be relevant contributions to σv x F from (χχ → Φ i Φ j ) annihilations, where Φ i denotes a scalar or pseudoscalar Higgs mass eigenstate. Such processes can be mediated either via diagrams with a Higgs or a Z boson in the s-channel, or via t/u-channel exchange of the Majorana fermion χ or the (heavy) SU (2)-doublet Dirac fermion we integrated out. In our EFT, the last possibility proceeds via the χχΦ i Φ j contact interaction terms in Eq. (2.13). Regardless of the type of diagram, the annihilation into a pair of CP-even (χχ → h i h j ) or CP-odd Higgs bosons (χχ → a i a j ) is p-wave suppressed, while (χχ → h i a j ) annihilations are swave. The corresponding s-wave contribution to the thermally averaged annihilation cross section is given by where the sum includes the s-wave amplitudes mediated by CP-odd scalars the amplitude mediated by transversally polarized Z bosons in the s-channel the amplitudes mediated by the Dirac fermion in the t/u-channel proceeding via contact interaction terms after integrating out the Dirac fermion Ψ 55) and the amplitude from the t/u-channel exchange of the Majorana fermion χ are the dimensionful trilinear couplings between different Higgs mass eigenstates (between the neutral Goldstone mode and two Higgs mass eigenstates), which are not related to the parameters of our EFT, but arise from the Higgs potential (see Appendices of Ref. [29] for details). The g χχΦ i Φ j are the (χχΦ j Φ k ) couplings of dimension (mass) −1 which can be read off from the Lagrangian Eqs. (2.13) or (A.1) taking into account the mixing of the Higgs mass eigenstates, Eqs. (2.23), (2.24).
After accounting for suppression of couplings arising from the requirement of an m h = 125 GeV SM-like Higgs mass eigenstate and from approximately satisfying the blindspot condition, the most relevant final states for the (χχ → Φ i Φ j ) processes are (χχ → a S h) and (χχ → a S h S ). If kinematically accessible, they can compete with (χχ → tt) annihilation. For both these channels, the amplitudes mediated by the singlet-like CP-odd a S in the s-channel may play an important role. However, their relevance to the total cross section is dictated by the coupling strengths g ha S a S and g h S a S a S respectively, which as mentioned above are not related to our EFT parameters. Hence for simplicity, we will assume that these couplings are small, rendering these processes irrelevant for the relic density.
Ignoring such Higgs exchange diagrams, the amplitude A h i a j Ψ mediated by a t/u-channel Dirac fermion, which we integrated out, is most relevant for the final state (χχ → a S h). Ignoring the kinematic correction in Eq. (2.52) which is relevant only very close to threshold (m h i + m a i = 2m χ ), canonical values of the thermally averaged annihilation cross section σv x F ∼ 2 × 10 −26 cm 3 s −1 can be achieved for where we have assumed the mixing of the CP-odd Higgs bosons to be small. This corresponds to couplings (2.58) For the channel (χχ → a S h S ) the processes associated with a Majorana fermion χ in the t/u-channel can be relevant. Assuming again for simplicity that these processes dominate compared to the one associated to the interchange of a singlet pseudoscalar (i.e. g h S a S a S is small), neglecting the threshold corrections for (m h S + m a S ≈ 2m χ ), and corrections from singlet-doublet mixing (the latter potentially leading to O(1) suppression), the thermally averaged annihilation cross section σv x F ∼ 2 × 10 −26 cm 3 s −1 can be achieved for a DM coupling to the singlets This implies where we assumed v µ for the estimate on κ such that g χχh S ∼ √ 2κ. Annihilations into pairs of vector bosons [χχ → ZZ(W + W − )] do not play an important role for obtaining the thermal relic density as long as m χ > m t . Final states consisting of two vector bosons with longitudinal polarizations are p-wave suppressed since they correspond to annihilations into a pair of CP-odd scalars (i.e. the neutral and charged Goldstone modes for the Z and W bosons, respectively). Annihilations into a pair of transversally polarized vector bosons or one transversally polarized and one longitudinally polarized vector boson are s-wave. However, such annihilations proceeding via t/u-channel exchange of the neutral (charged) components of the SU (2)-doublet fermion we integrated out correspond to χχZZ (χχW + W − ) contact interaction terms in our EFT which would only appear at dimension d ≥ 7 and hence are strongly suppressed. For m χ > m t , annihilations mediated by an s-channel Higgs or Z boson are also suppressed: the coupling of the s-channel mediator to one transversally and one longitudinally (a pair of transversally) polarized vector bosons is proportional to the gauge couplings (squared). The couplings of the Higgs and Goldstone bosons to quarks, instead, are proportional to the top Yukawa coupling. The Higgs mediated channel is furthermore also p-wave suppressed.
To summarize, the proper value of the thermally averaged annihilation cross section σv x F ∼ 2×10 −26 cm 3 s −1 leading to the observed relic density can be easily obtained when (χχ → tt) annihilations are kinematically allowed, i.e. when m χ > m t . The annihilation cross section will then typically be dominated by annihilations into top quarks mediated by the neutral Goldstone mode for which the proper value of σv x F is achieved for g χχG 0 ∼ 0.1, cf. Eq. (2.51). Annihilation into pairs of vector bosons is suppressed because it proceeds through smaller couplings. If kinematically allowed, the annihilation cross section into pairs of Higgs mass eigenstates may become large enough, cf. Eqs. (2.58), (2.59), although the annihilation into pairs of top quarks tends to be competitive or dominant unless the EFT parameters conspire to suppress the g χχG 0 coupling.
For lighter DM candidates, m χ < m t , achieving a sufficiently large annihilation cross section σv x F ∼ 2 × 10 −26 cm 3 s −1 is more difficult. In this case, (χχ → bb) annihilation is dominated by the Z boson mediated p-wave contribution, and couplings are usually not sufficiently large to obtain the proper annihilation cross section. Annihilation into pairs of Higgs bosons requires large couplings between the different Higgs mass eigenstates in order to be sufficiently effective; in addition it is not easy to obtain a light enough Higgs mass spectrum (m h i + m a i < 2m χ ) while simultaneously satisfying collider constraints. Annihilation into pairs of vector bosons usually does not achieve sufficiently large cross sections either since they are controlled by gauge couplings. Hence, for m χ < m t avoiding over-closure of the Universe requires large couplings in the Higgs sector unless the annihilation cross sections via a Higgs boson Φ (a Z boson) in the s-channel is boosted via resonant annihilation, 2m χ ≈ m Φ (2m χ ≈ m Z ).

EFT: Direct Detection
Elastic (χq − χq) scattering relevant for direct detection proceeds via the same diagrams as annihilation in the early Universe depicted in Fig. 1, but interpreting them as t-channel exchanges of Higgs bosons. The exchange of CP-odd Higgs bosons leads to spin-dependent scattering, and the contribution of the Goldstone mode G 0 is furthermore suppressed by q 2 /m 2 Z , where q is the momentum transfer. In contrast, the exchange of CP-even Higgs bosons leads to spin-independent scattering. Since bounds on the spin-independent DMnucleon scattering (SIDD) cross section [1][2][3][4] are much stronger than those on the spindependent scattering (SDDD) cross section [5,6] we focus on SIDD in the remainder of this section.
Summing coherently over the contribution from all CP-even Higgs mass eigenstates, the χ-proton SIDD cross section can be written as where m p is the mass of the proton, and the form factors (at zero momentum transfer) are 2 13. The (a q /m q ) i parametrize the contribution to the scattering amplitude from one Higgs mass eigenstate, where the g χχh i and g qh i are given in Eqs. (2.40) and (2.41). As discussed in Section 2.1, the observation of the 125 GeV Higgs boson at the Large Hadron Collider (LHC) with couplings close to that of a SM Higgs implies that our model must contain a CP-even Higgs eigenstate h with m h ≈ 125 GeV and composition S SM In addition, to avoid bounds on additional Higgs bosons from the LHC and the Large Electron-Positron Collider (LEP), the remaining CP-even mass eigenstates H and h S must be either heavy m H m h or dominantly composed of H S . The coupling of the H NSM Higgs boson to down-type quarks is enhanced by tan β, and therefore at large values of tan β the suppression induced by its large mass may be compensated by an enhancement of the coupling. This case allows for effective destructive interference between the H SM and H NSM contributions to the SIDD cross section [15]. For values of tan β = O(1), as we shall use in our work, the contribution of the non-standard Higgs bosons to the SIDD cross section will either be suppressed by (a q /m q ) 2 ∝ 1/m 4 Under such conditions, the SIDD cross section will be dominated by the contribution from the SM-like state h. Hence, the current stringent bounds on the SIDD cross section lead to relevant constraints on g χχh .
The bounds on g χχh may be estimated by considering the SIDD cross section, taking into account only the diagrams with an h in the t-channel. We find from Eq. (2.61) while the experimental limit is σ SI p (m χ = 300 GeV) 3.3 × 10 −10 pb [4]. Hence, values of g χχh 0.025 are necessary to fulfill these constraints. This range of values of g χχh may be compared with the values of g χχG 0 ∼ 0.1 necessary to obtain the observed relic density as discussed in the previous section, cf. Eq. (2.51). In general, the couplings of χ to the Higgs mass eigenstates are expected to be of similar magnitude, or larger, than g χχG 0 , since they arise at the same order, or lower, in our EFT expansion. In particular, the coupling g χχh arises from dimension d = 5 operators with the leading contribution suppressed by v/µ, while g χχG 0 arises via d = 6 operators and is suppressed by m χ v/|µ| 2 . We therefore conclude that under the requirement of obtaining an acceptable relic density, the values of

Bounds on Couplings and Parameters of the EFT
In the previous section we argued that in order to suppress the SIDD cross section below experimental limits, the model must sit close to the blind spot, Eq. (2.16). It is also possible that the amplitude mediated by the h interferes destructively with those mediated by the other CP-even mass eigenstates H and h S , or both mechanisms may be at work simultaneously. As argued above, the amplitudes of the diagrams mediated by H and h S are in general suppressed with respect to the amplitude of the h-mediated diagram, such that destructive interference will only be effective when the contribution from the SMlike Higgs h is already suppressed by proximity to the blind spot condition. In order to demonstrate these properties, we shall briefly consider the simple case in which the singlet fields are heavy and play no role in the low energy effective theory. Fig. 2 shows a representative example, illustrating the bounds on the parameters of our EFT and the couplings of DM to the CP-even and CP-odd Higgs bosons necessary to satisfy the SIDD experimental constraints and obtain the proper relic density concurrently.
We used characteristic values of the DM mass, m χ = 300 GeV, and the non-standard Higgs boson masses m H NSM = m A NSM = 500 GeV, and moderate values of tan β = 2. The singlets are decoupled {m H S , m A S } 500 GeV and perfect alignment is assumed in the doublet sector such that the Higgs basis states coincide with mass eigenstates. In the left panel we show the couplings of the DM to the non-standard Higgs bosons as a function of the DM coupling to the SM-like Higgs boson. As explained in the last section, ignoring the other CP-even Higgs bosons contributions, the coupling to the SMlike Higgs boson needs to be suppressed to satisfy the SIDD constraints, g χχh 0.025. As can be seen in Fig 2, this bound may be slightly relaxed in the presence of destructive interference with other Higgs bosons. The orange shaded region denotes the values of the couplings of DM to the CP-even Higgs bosons consistent with the SIDD bounds, with the boundaries denoting the values of these couplings saturating the current SIDD bound σ SI p = 3 × 10 −10 pb (dashed and solid lines). The couplings of the CP-odd Higgs bosons are constrained only by the demand of obtaining a proper relic density. 3 The solid and dashed ellipses show the values of g χχA NSM (black) and g χχG 0 (blue) needed to obtain Ωh 2 ∼ 0.12, consistent with the SM-like Higgs couplings and the boundary values of the heavy CPeven Higgs couplings denoted by the solid and dashed lines, respectively. Thick and thin lines represent two different solutions for the CP-odd/Goldstone couplings for each set of values of the CP-even Higgs bosons couplings, with points in each of the ellipses in one-to-one correspondence with similar points in the other ellipses. The clear correlation between the couplings to the Goldstone mode and to the heavy CP-odd Higgs boson may be understood from the fact that for the given Higgs and DM masses, the suppression associated with the propagator contributing to the annihilation cross section is about a factor 3 weaker for the heavy CP-odd Higgs boson than for the Goldstone mode of mass For couplings of DM to the CP-even Higgs bosons in the shaded area, for which the SIDD cross section would be smaller than the experimental limit, the couplings of the CP-odd Higgs boson and the Goldstone mode would take intermediate values between the ones represented by the two ellipses.
The right panel of Fig. 2 shows the three independent combinations of the parameters of the EFT that play a role in the determination of the relic density and the direct detection constraints. The shaded regions, solid lines, and dashed lines are in one-to-one correspondence with those in the left panel, and the shown correlations may be easily understood from the dependence of the couplings on these parameters, Eqs. (2.14), (2.18), (2.20) and (2.21).

NMSSM Ultraviolet Completion
The EFT model discussed in the previous section extends the SM particle content by an additional SU (2)-doublet Higgs field, as well as a complex scalar field and a Majorana fermion both of which transform as singlets under the SM gauge group. There is also an additional SU (2)-doublet Dirac fermion which is assumed to be heavy and integrated out, yielding interactions of the Majorana fermion with the Higgs doublets. This particle content is very similar to the Higgs and neutralino sector of the NMSSM. Hence, decoupling the gluinos as well as the SUSY partners of the SM fermions from the theory, the NMSSM can serve as an ultraviolet completion for the EFT, preserving all the essential features discussed above. In addition, a range of mature numerical tools are available for the NMSSM which allows us to compute particle spectra and couplings and subsequently study the DM and collider phenomenology of the model. As we shall see, the NMSSM also contains a second region of parameter space where the DM candidate does not transform under the Z 3 symmetry, however, this region can nonetheless quite simply be mapped onto our EFT. In the SUSY case such a Z 3 symmetry is defined so that all chiral superfields transform by e 2πi/3 , while gauge superfields transform trivially.
The Higgs sector is analogous to that of our EFT model, and can be rotated to the extended Higgs basis using Eqs. (2.3)-(2.6). Including the usual soft SUSY breaking and F -and D-terms [26], and ignoring radiative corrections for presentational purposes, see e.g. Ref. [29], the symmetric squared mass matrix for the neutral CP-even Higgs bosons in the basis {H SM , H NSM , H S } is 4 3) where c β ≡ cos β, s β ≡ sin β, and we have introduced In the basis {A NSM , A S }, the symmetric tree-level squared mass matrix for the CP-odd neutral Higgs bosons is Radiative corrections may be relevant and are incorporated at the two-loop level in the numerical results obtained with NMSSMTools [58]. As in the case of our EFT model, the NMSSM must contain a CP-even 125 GeV mass eigenstate mostly composed out of H SM to accommodate the SM-like Higgs boson observed at the LHC [41,42]. This can be achieved either by making the remaining CPeven states H and h S heavy, 5 M S,22 , M S, 33 M S,11 , the decoupling limit, or by setting M 2 S,12 ≈ M 2 S,13 ≈ 0, the alignment-without-decoupling limit [29,32]. Perfect alignment is achieved for [29] and in the alignment limit the mass of the SM-like Higgs mass eigenstate is given by where ∆(mt h ) are radiative corrections that are common to the MSSM. Note, that with respect to the MSSM, one obtains an additional contribution (λvs 2β ) 2 to m 2 h which allows for a 125 GeV SM-like Higgs boson mass without the need for large radiative corrections for moderate values of tan β 3 if λ ∼ 0.7. Intriguingly, the first alignment condition in Eq. (3.6), which suppresses the H SM − H NSM mixing, is satisfied for the same values of λ, namely 0.6 λ 0.7. Alignment with the singlet [the second condition in Eq. (3.6)] is also easily achieved by judicious choices of M A and µ.
The neutralino sector of the NMSSM consists of the superpartners of the neutral electroweak gauge bosons, the bino B and neutral wino W 3 , the neutral Higgsinos H 0 d and 4 Note, that Ref. [29] uses the parameter M Z ≡ m 2 Z − λ 2 v 2 and the v = 246 GeV convention, while we use v = 174 GeV. 5 Note, that we use the same notation for the Higgs states as in Section 2.
H 0 u belonging to the respective Higgs doublets, and the singlino S, the fermionic component of S. In the basis { B, W 3 , H 0 d , H 0 u , S}, the symmetric tree level neutralino mass matrix is where M 1 and M 2 are the bino and wino soft SUSY breaking masses. The neutralino mass eigenstates are given in terms of the interaction eigenstates by   .
In terms of the mixing angles, the couplings of the lightest neutralino to the Higgs basis states are 14) From these, the couplings to the Higgs mass eigenstates can be obtained using Eq. (2.40) and the mixing angles of the Higgs mass eigenstates. The singlino will play the role of the Majorana singlet in our EFT model. However, the NMSSM contains a second SU (2)-singlet neutralino, the bino. Unlike the singlino, the bino does not transform under the Z 3 symmetry since it stems from a gauge superfield. Besides allowing for an explicit mass term ( 1 2 M 1 B B + h.c.), this also results in different couplings to the Higgs doublets and the singlet than those found for the singlino. The region where the bino is the DM candidate can nonetheless be connected to our EFT in a straightforward way, as we will see in the following section.

Top-down EFT
In order to connect the NMSSM to our EFT from Section 2, we can construct a top-down EFT for the NMSSM Higgs and neutralino sector by integrating out the Higgsinos. This approach is valid as long as the Higgsino mass is large compared to the mass of the lightest neutralino. We will also assume the winos to be heavy, M 2 {M 1 , µ}. Neglecting the Yukawa terms and ignoring the charged current interactions, the terms in the Lagrangian involving the neutral components of the Higgsinos are The corresponding equation of motion for H 0 u is and for H 0  20) and (3.21) Substituting these into Eq. (3.17) and keeping terms of dimension d ≤ 6 (the same order used for our generic EFT in section 2), we find where in the last line we have included the standard bino mass term and the singletsinglino-singlino interaction term. Note, that because the singlino transforms under the Z 3 symmetry, it only gets a mass when the singlet S acquires a vev. In contrast, the bino does not transform under the Z 3 and hence a soft SUSY breaking mass term M 1 is allowed. Expanding S around its vev S → µ/λ +Ŝ and correspondingly we can appreciate the similarities with the Lagrangian of our EFT model. In particular, the singlino has the same structure for the couplings as the Majorana fermion χ in Section 2, and we can map the couplings in Eq. (2.13) directly to those in the NMSSM via of (H 0 d )(H 0 u )χχ, which can be compensated for by changing the sign of the coupling of the binos to the CP-odd states coming from the Higgs doublets, i.e. A NSM and G 0 . Keeping this in mind, we can map the couplings of the bino in Eq. (3.22) to those in the EFT, Eq. (2.13), via The blind spot condition for the bino region is then sin 2β = −m χ /µ . (3.27) Note, that the bino couples with characteristic strength g 2 1 /2 ≈ 0.06 to Higgs doublet states, whereas the singlino couples to Higgs doublet states with characteristic strength λ 2 ∼ 0.4, recalling that the presence of the 125 GeV SM-like Higgs implies λ ∼ 0.6.
The couplings of pairs of (on-shell) singlinos to the extended Higgs basis states and the (neutral) Goldstone mode can be directly read off from Eq.
Similarly, the couplings of pairs of (on-shell) binos are given by 37) in agreement with the mapping of parameters in Eq. (3.26), keeping in mind the switch in the sign of the (χχA NSM ) and (χχG 0 ) couplings and the additional mass term m χ = M 1 . Note that Eq. (3.22) also induces bino-singlino-Higgs couplings which might be relevant for thermal production via co-annihilation when the bino and singlino are approximately mass degenerate, M 1 ≈ 2κµ/λ. For on-shell binos and singlinos, these couplings are

Dark Matter phenomenology
Besides serving as an ultraviolet completion of our EFT, the NMSSM also provides a convenient computational basis since mature numerical tools are available for the analysis of both collider and DM phenomenology. In this section, we study the phenomenological properties of the NMSSM, going beyond the EFT validity conditions. Doing so, we can identify those properties that are shared with the EFT approach, while also determining the differences between the EFT and the full theory predictions. We use NMSSMTools 5.1.2 [58-62] to compute NMSSM spectra and couplings and to subsequently test parameter points against a subset of the constraints implemented in NMSSMTools (see Ref. [58] for details). In particular, points are excluded if they feature an unphysical global minimum, soft Higgs masses much larger than the SUSY scale, and if the lightest neutralino is not the lightest SUSY particle. Furthermore, we require the spectrum to contain an m h ≈ 125 GeV Higgs boson with couplings to SM particles compatible with those of the SM-like Higgs observed at the LHC. We also require compatibility with constraints from LEP, Tevatron, and the LHC on additional Higgs bosons and sparticles as implemented in NMSSMTools. For points passing these constraints we compute the relic density and the direct detection cross section with micrOMEGAs 4.3.5 [44-47, 63, 64].
We perform a random scan over 10 9 parameter points, drawing the parameters from linear-flat distributions over the ranges listed in Table 1. The choice of parameter ranges is motivated by the phenomenology of a SM-like Higgs: for tan β ≤ 5, a 125 GeV SMlike Higgs boson is obtained without the need for large radiative corrections to its mass, cf. Eq. (3.7). The range for A λ is chosen to be larger than those for µ and A κ because approximate alignment implies, from Eqs. (3.4) and (3.6), The range of the bino mass M 1 is chosen such that we obtain both the case where the lightest neutralino χ 1 is mostly bino-like (i.e. when M 1 < {µ, 2κµ/λ}) and the case where the bino component of χ 1 is negligible (i.e. when min(µ, 2κµ/λ) M 1 ). In addition to the parameters listed in Table 1, we scan over the stop mass M 3 Q = M 3 U in the range [0.75; 2.5] TeV and, for definiteness, we set the stop and sbottom mixing parameters X t ≡ (A t −µ cot β) = 0 and X b ≡ (A b −µ tan β) = 0, 6 allowing for moderate radiative corrections to the Higgs masses. We decouple the remaining SUSY particles by setting all remaining sfermion mass parameters to 3 TeV, the gluino mass parameter to M 3 = 2 TeV, and, in order to minimize the wino component of the lightest neutralino, the wino mass parameter to M 2 = 10 TeV. Note, that due to this choice of parameters we also satisfy all LHC bounds on sfermions and gluinos.
Due to our choice of parameters M 2 {M 1 , |µ|, |2κµ/λ|}, the heaviest neutralino χ 5 will be wino like, while the wino component of the lighter mass eigenstates will be negligible. Furthermore, since we chose |κ| ≤ 0.3, the singlino mass parameter |2κµ/λ| will practically always be smaller than the Higgsino mass parameter |µ|, while we chose the range of the bino mass parameter M 1 such that the bino can be both lighter and heavier than the singlino and the Higgsinos. This also ensures that the Higgsinos with mass parameter µ are always heavier than the lightest neutralino χ 1 , such that we omit the phenomenologically disfavored Higgsino DM region and we can map results onto our EFT, where m χ 1 /µ appears as an expansion parameter. Hence, our DM candidate, the lightest neutralino, will be either mostly singlino-like or mostly bino-like with small Higgsino admixture facilitating Higgs mediated processes. In all figures presented in this section, we compute the couplings of the neutralinos to Higgs mass eigenstates and of the SM particles to the Higgs mass eigenstates from the mixing angles N ij for the neutralinos and S ij (P ij ) for the CP-even (CP-odd) Higgs bosons as output from NMSSMTools. Computing the couplings from mixing angles, Eqs. (3.12)-(3.16), takes into account sub-leading effects from the (small) admixture of additional neutralino components. These would only appear in the EFT through higher dimensional operators, or for wino effects from including operators arising from integrating out an SU (2)-triplet fermion in addition to the operators from integrating out an SU (2)-doublet fermion. Furthermore, using the output from NMSSMTools for the mixing angles also captures effects induced by the running of the parameters from the SUSY scale to the electroweak scale.
In the NMSSM, apart from the processes discussed in detail in section 2 , the annihilation cross section can be enhanced either by resonant annihilation or by co-annihilation [65]. For the non-resonant annihilation cross sections, we checked numerically that for points with an acceptable relic density Ωh 2 ∼ 0.12, annihilations into top quarks typically make up for O(80 %) of the thermally averaged cross section σv x F , while annihilations into pairs of Higgs bosons usually account for O(20 %) of σv x F . This can be quiet generically understood to be due to the difficulty in obtaining a Higgs spectrum light enough to allow for the second final state, while evading Higgs phenomenology and collider constraints.
Concentrating on (χχ → tt), Fig. 3 shows the SIDD cross section vs. the contribution to the thermal annihilation cross section from the diagrams mediated by the Goldstone mode obtained from Eq. (2.45) when taking into account only the Goldstone amplitude in Eq. (2.46). The left panel shows all points from our scan passing the Higgs and LHC constraints described above but before requiring the correct relic density or compatibility with direct detection limits. The dominant composition of the DM candidate is color coded and denoted as bino B if N 2 11 ≥ 0.95 (red), singlino S if N 2 15 ≥ 0.95 (blue). Similarly, points are denoted as mixed if the sum of the square of corresponding mixing angles is N 2 1i ≥ 0.95 but none of the individual contributions is sufficiently large to put them in one of the previous categories. The behavior shown can be understood quite intuitively from our findings in sections 2 and 3.1: first of all, we note that both the SIDD and the annihilation cross section for mostly bino DM are smaller than those for mostly singlino DM because the couplings to Higgs bosons are proportional to g 2 1 /2 ≈ 0.06 for binos while for singlinos the couplings are proportional to λ 2 ∼ 0.4. Secondly, for all the different DM compositions, a striking feature of this plot is the relative independence of the SIDD cross section and σv G , which can be suppressed independently from each other. The SIDD cross section is suppressed close to the blind spot conditions m χ /µ → ± sin 2β. However, in this case the coupling to the Goldstone mode g χχG 0 ∝ m χ cos 2β/|µ| 2 remains sizable and thus the corresponding contribution to the annihilation cross section is not suppressed. Note, that the contributions to the thermal annihilation cross section from the Higgs mass eigenstates are usually smaller than those from the Goldstone mode: while the couplings are typically of the same order, the contributions from the CP-odd eigenstates are suppressed by and the contributions from CP-even Higgs bosons are p-wave suppressed. As argued in section 2.2, the contribution from the transverse polarizations of the Z boson is generally also much smaller than the contribution from the Goldstone mode. The right panel of Fig. 3 shows the points from our scan having approximately the correct relic density Ωh 2 = 0.12 ± 50 % and satisfying bounds from SIDD and SDDD experiments [3][4][5][6]. For mostly singlino DM, the contribution from the Goldstone mode mediated amplitudes to the thermally averaged annihilation cross section is of the order of the canonical value σv ∼ 2 × 10 −26 cm 2 /s yielding the observed relic density Ωh 2 ∼ 0.12 [43] 7 . In contrast, we find that for mostly bino DM the Goldstone mode does not mediate sufficiently large amplitudes to avoid over-closure of the Universe. Hence in the bino DM scenario, both resonant annihilation and co-annihilation play a very important role. In Fig. 4 we show points from our numerical scan which have an acceptable relic density, Ωh 2 = 0.12 ± 50 %, with the color coding indicating the possibility of resonant annihilation. In the right panel of Fig. 4 we show points with mostly singlino DM candidate in the (singlino mass)-(Higgsino mass) plane, demonstrating that neither co-annihilation nor resonant annihilation is relevant for the singlino region (M 1 is always large in this region). In the left panel of Fig. 4 we show bino-like points in the (bino mass)-(singlino mass) plane. We find that points either resonantly annihilate via the Z boson or one of the Higgs mass eigenstates, or, feature binos approximately mass degenerate with the singlino such that co-annihilation yields the correct relic density. In the latter case, it is in fact the annihilations of the mostly singlino like m χ 2 which set the relic density. We have thus found a new well tempered bino region: m χ 1 is mostly bino-like with very small couplings, evading direct detection constraints easily. However due to the presence of an almost mass degenerate singlino-like m χ 2 (which does not play a role in direct detection) which has significantly larger couplings, an observationally consistent relic density is easily obtained. The value of the µ parameter, and consequently the Higgsinos, tends to be about the same order as shown for the singlino-like DM in the right panel.
In the left panel of Fig. 5 we show the SIDD cross section vs. m χ /(µ sin 2β) for points from our scan passing the Higgs and collider constraints described above but before requiring the correct relic density or compatibility with direct detection limits. For bino DM, the SIDD cross section is suppressed when the blind spot condition m χ /µ = − sin β is approximately satisfied, while for all other compositions of the DM candidate, usually singlino dominated, the SIDD cross section is suppressed for m χ /µ ≈ sin β. Note that for bino DM candidates we also find suppression of the SIDD cross section for m χ /µ ≈ 0, which corresponds to M 1 µ for which we do not find any parameter points with an acceptable relic density in our scan. This is because in this region neither co-annihilation with the singlino nor resonant annihilation with the Higgs bosons is possible, one of which would be required to boost the thermal annihilation cross section to avoid over-closure of the Universe.
In the right panel of Fig. 5 we show the SIDD cross section vs. the DM coupling to the SM-like 125 GeV Higgs mass eigenstate. Besides the coupling to the H SM Higgs basis state, which vanishes at the respective blind spots, this takes into account the contributions to the coupling from the (small) admixtures of the H NSM and H S interaction eigenstates to the 125 GeV mass eigenstate. We find that for mostly bino points the SIDD cross section is very tightly correlated with the coupling to the 125 GeV SM-like Higgs mass eigenstate, and SIDD cross sections satisfying the current experimental bounds can be achieved by suppression of the g B Bh coupling. In the case of mostly singlino DM we find this correlation to be looser, indicating that the SIDD cross sections must be suppressed by additional mechanisms.
In Fig. 6 we show the contributions to the SIDD cross section when taking into account only one Higgs mass eigenstate at a time as obtained from Eq. (2.61) ignoring the sum over the CP-even Higgs mass eigenstates. We show these contributions plotted against the full SIDD cross section in the left (right) panel for mostly bino (singlino) points from our dataset satisfying all Higgs/collider constraints described above and featuring an acceptable relic density. For mostly bino DM, we find that SIDD cross sections as small as σ SI p ∼ 10 −13 pb can be obtained by suppression of the coupling to the SM-like Higgs mass eigenstate alone. Destructive interference between different Higgs mass eigenstates is needed only for even For mostly singlino DM, as shown in the right panel of Fig. 6, destructive interference between different Higgs mass eigenstates is almost always required to satisfy the experimental bounds on the SIDD cross section. This can be understood from the typical strength of singlino couplings σ SI p ∝ g 2

B Bh
∝ g 4 1 /4 ∼ 0.004. In addition, compared to bino DM, singlino DM has a much larger coupling to the scalar singlet state H S due to the presence of the tree-level coupling κ, cf. Eqs. (3.30) and (3.35). Hence we see the necessity of destructive interference between the contributions from the singlet like mass eigenstate h S and the SM-like mass eigenstate h to suppress the SIDD cross section below the experimental limits.
Although blindspot cancellation or destructive interference arguably require some fine tuning, we stress that we readily find points in our dataset with SIDD cross sections below σ SI p 10 −13 pb, out of reach of direct detection experiments for the foreseeable future. Such small cross sections are challenging to probe with current direct detection strategies due to the presence of the neutrino floor.
In Fig. 7 we show the constraints from direct detection experiments for points from our scan with an acceptable relic density and satisfying collider constraints. The left panel shows the SIDD cross section vs. the DM mass. Note, that almost all parameter points with DM masses below the top mass m χ < m t ≈ 173 GeV are ruled out by current SIDD Figure 7. Left: The SIDD cross section vs. the mass of the DM candidate m χ for points from our parameter scan passing the required experimental collider constraints with acceptable relic density Ωh 2 = 0.12±50 %. The solid black and dashed red lines indicate the most constraining experimental upper limits on the SIDD cross section from XENON1T [3] and PandaX-II [4], respectively. Right: SIDD vs. SDDD cross section in units of the respective observed limit for the same points. For SIDD cross section, at each respective DM mass we use the stronger of the two limits from XENON1T and PandaX-II. For SDDD scattering we use the more constraining of the current bounds for either SDDD scattering of neutrons from LUX [5], or SDDD scattering of protons from PICO-60 [6]. To guide the eye we indicate the current bounds with thin dashed lines; points lying in the lower left quadrant satisfy all current direct detection bounds. The color coding in both panels is the same as in Fig. 3. constraints. The right panel shows a comparison of the constraining power of current SIDD and SDDD experiments [3][4][5][6] for the same points. As noted above, current experimental SIDD limits are much more constraining for our model than SDDD limits: as we can see in Fig. 7, parameter points satisfying SIDD limits almost always satisfy current constraints on the SDDD scattering cross section from direct detection experiments, while the reverse is not true. Fig. 7 also shows that improvements in SDDD experiments will probe the remaining parameter space of our model more efficiently than a correspondingly large improvement in the bounds from SIDD. In particular, improving the sensitivity of SDDD scattering by approximately two orders of magnitude with respect to current bounds will probe most of the singlino region. A naïve estimate of the sensitivity can be obtained by rescaling current bounds from LUX on the spin-dependent WIMP-neutron scattering cross section which were obtained with an exposure of = 129.5 kg yr [5]. The XENON1T experiment will have a sensitivity roughly one order of magnitude better assuming an exposure of = 1 t yr, while XENONnT could probe most of the singlino region assuming an exposure of = 20 t yr.
Finally, let us mention that out of the points passing the main phenomenological constraints described at the beginning of this section, which in particular lead to a SM-like Higgs boson compatible with LHC measurements and evade collider constraints as implemented in NMSSMTools, ∼ 15 % have an acceptable relic density Ωh 2 = 0.12 ± 50 %, and of the points with the correct relic density ∼ 35 % satisfy SDDD and SIDD constraints.

Indirect Detection
In Section 4 we have focused on direct detection constraints on our model. Any DM candidate is also subject to constraints from indirect detection experiments. Due to large astrophysical uncertainties on indirect detection constraints arising from charged particle production, we focus on constraints from photon emission here. The currently strongest indirect detection constraints on WIMP DM come from Fermi-LAT and MAGIC analyses of Milky Way satellite galaxies, ruling out canonical values of the thermally averaged annihilation cross sections σv ∼ 2 × 10 −26 cm 3 s −1 for DM masses m χ 100 GeV if the annihilation is dominantly into pairs of b quarks or τ leptons [69,70]. For DM masses m χ m t indirect detection bounds are much weaker, constraining thermal annihilation cross sections much larger than the typical values σv x F ∼ 2×10 −26 cm 3 s −1 consistent with an acceptable relic density. Relevant constraints might further arise from the annihilation of DM matter particles captured via elastic scattering in the Sun, potentially giving rise to a flux of neutrinos observable at Earth. Under certain assumptions [71], most relevant for our case if DM annihilates dominantly into pairs of W bosons, current bounds from DM capture and annihilation yield the most constraining upper limits on the spin dependent DM-proton scattering cross section [72][73][74][75][76].
The parameter region in our model yielding an acceptable relic density and compatible with current SIDD constraints typically features DM masses heavier than the top quark m χ m t ≈ 173 GeV. Furthermore, SIDD constraints require the coupling of DM to the SMlike Higgs boson to be suppressed which in turn implies that (χχ → W + W − ) annihilation is suppressed. This is because of all the Higgs basis states, only H SM couples to pairs of W bosons. Thus, current indirect detection limits do not pose relevant constraints on our model since the annihilation modes most constraining for indirect detection bounds (χχ →bb/τ + τ − /W + W − ) are suppressed for our preferred region of parameter space. Note however, that DM with mass m χ > m t decaying dominantly into pairs of top quarks may explain the excess of photons emitted from the Galactic Center observed by Fermi-LAT [77][78][79][80].

Collider Constraints
Our EFT DM model extends the neutral particle content of the SM by four Higgs mass eigenstates, two of them CP-even and two CP-odd, and a SM singlet Majorana fermion. The Dirac fermion we integrated out may be accessible at the LHC if its mass µ is below the TeV scale [53]. Note, that in our NMSSM ultraviolet completion we set the masses of all additional particles such as the SUSY partners of the SM fermions and the charged superpartners of the SM gauge bosons to values ≥ 2 TeV such that they are not directly accessible at the LHC. 8 Furthermore, DM is coupled to the SM via the Dirac fermion (the Higgsinos in the case of the NMSSM) which in turn couples to the Higgs and the Z bosons. Thus, this model can be tested at the LHC either by DM pair production via Z or Higgs bosons, or, via production of the additional Higgs bosons.
As discussed previously in Sections 2.1 and 3, the observation of a 125 GeV Higgs boson with couplings compatible with those of a SM Higgs [41,42] requires the presence of an m h ∼ 125 GeV Higgs mass eigenstate approximately aligned with the H SM Higgs basis interaction eigenstate. Direct detection limits furthermore strongly constrain the coupling of DM to both the SM-like Higgs and the Z boson, such that the additional decay width of the Z boson or the SM-like Higgs into DM particles must be small, even before taking into account the fact that m χ m h /2 > m Z /2 in our model. The non-observation of signals from additional Higgs bosons at the LHC furthermore implies that additional Higgs bosons must either have masses heavier than 2m t ∼ 350 GeV, or, be dominantly composed of the additional singlet interaction states H S or A S such that their production cross section is suppressed. A comparison of the constraining power of conventional searches for additional Higgs bosons decaying into pairs of SM particles can for example be found in the Appendix of Ref. [32]. Note, that most of the conventional Higgs searches such as (gg → H/A → τ + τ − /γγ/Zh) will only have small increases in sensitivity at future LHC runs because they are increasingly dominated by systematic errors. On the other hand, models with extended Higgs sectors and potentially large couplings between different Higgs mass eigenstates such as ours can be probed effectively by decays of heavy Higgs bosons into lighter Higgs bosons or a light Higgs and a Z boson. Due to the presence of the SM-like Higgs, decays into pairs of SM-like Higgs bosons (H → hh) or into a SM-like Higgs boson and a Z (A → Zh) are suppressed since the corresponding couplings are proportional to the mixing of the H SM Higgs basis state with the H NSM and H S states. If kinematically allowed, the branching ratios of (H → hh S /Za S ) and (A → Zh S /ha S ) can however be sizeable. The collider signatures arising through such decays can be effective probes of our model [29,32,33].

Conclusions
Current bounds from direct detection experiments [1][2][3][4][5][6] have set stringent limits on WIMP DM scenarios where the relic density proceeds from thermal production mediated by Higgs and gauge bosons. In particular, direct detection experiments strongly constrain the coupling of the SM-like Higgs boson to DM particles and the vector-like coupling of DM to Z bosons [7]. The latter constraint is naturally satisfied for Majorana DM, since Majorana fermions couple to Z bosons only via axial couplings. In this article, we explore an EFT describing the interactions of Majorana fermion WIMPs with an extended Higgs sector, comprised of a type II 2HDM and a SM gauge singlet. This model can be interpreted as an extension of Higgs portal models. In particular, we study the case where the EFT mass scale is identified with the mass of a heavy SU (2)-doublet Dirac fermion integrated out from the theory. Furthermore, we assume that all explicit mass terms are forbidden, which from the LHC allow for much smaller masses of the SUSY particles, probeable at ongoing and future LHC runs. may be realized by imposing a Z 3 symmetry. In such a case the new fermions acquire mass through the vev of the singlet field. If this vev is generated by the same mechanism that induces electroweak symmetry breaking, we naturally obtain a weak scale mass for the DM candidate. In addition, the Z 3 symmetry reduces the number of allowed operators involving the singlet Majorana fermion and the doublet and singlet scalar fields, such that we can easily extend the EFT to dimension 6, including derivative operators. We derive the low energy Higgs spectrum and the couplings of DM to the neutral Higgs and gauge bosons in this model, and use these results to compute the relic density and spin independent direct detection cross section. We find blind spot conditions that allow for the evasion of spin independent direct DM detection constraints, while yielding the observed relic density. These blind spot conditions can be easily characterized in terms of the EFT parameters.
In order to test the validity of the results derived in our EFT model, we compare them to those obtained in the supersymmetric NMSSM, which features an analogous Higgs sector. We demonstrate that in the case of heavy gauge and fermion superpartners the NMSSM can be reduced to our EFT model at low energies, such that we can use it as an explicit computational basis. The Majorana singlet in our EFT is identified with the singlino in the NMSSM, and the EFT scale with the mass of the heavier Higgsinos. We show that the qualitative features obtained in the EFT are preserved in the complete theory, while the precise values of the couplings, relic density, and cross sections are modified by the presence of small corrections which could be included by dimension d > 6 operators in the EFT description. We also discuss the case of a Z 3 invariant Majorana fermion DM candidate, e.g. the bino in the NMSSM, and show that our EFT model can be generalized in a straightforward manner to include this scenario.
Both in the EFT and in the NMSSM we show that the coupling of (singlino-like) DM to the SM-like Higgs is constrained to values below g χχh 0.1 by direct detection experiments, while simultaneously consistency with the DM relic density can be obtained through thermal annihilation via couplings of DM to the remaining Higgs and gauge bosons. The neutral Goldstone mode, comprising the longitudinal mode of the Z boson, plays a prominent role in obtaining the thermal annihilation cross section. Moreover, in order to evade direct detection bounds, not only the coupling to the SM-like Higgs boson must be reduced, but in addition destructive interference of the SM-like Higgs boson with the singlet and/or non SM-like CP-even doublet Higgs boson is required to further suppress the cross section.
In the NMSSM, when considering the case of light binos, we find a new well tempered DM region. For bino-like DM the couplings to the Higgs bosons tend to be smaller than for singlino-like one and direct detection bounds are hence evaded by a mild proximity to the blind spot conditions for the DM coupling to the SM-like Higgs state. The relic density in the well tempered region is obtained via co-annihilation of the bino with the singlino. Beyond the well tempered region, the correct relic density for bino-like DM may also be obtained via resonant annihilation with either the Z boson or a Higgs mass eigenstate.
In both the cases of bino-like and singlino-like DM, we find that the current constraints coming from SDDD are subdominant compared to those coming from SIDD. However, we also find that an improvement of two orders of magnitude in the bounds on the SDDD cross sections could efficiently probe most models consistent with a blind spot in SIDD.
Finally, collider searches for the Majorana DM particles are mostly restricted to the usual channels of single SM-particles plus missing energy, in particular those associated with Higgs bosons. Going beyond the effective theory, these searches may be complemented by missing energy searches in the production and decay of the Dirac doublet fermion and, in the NMSSM, by related searches for heavy superpartners.

A EFT Lagrangian
Rescaling the Majorana fermion field χ in order to retain a canonical kinetic term and keeping terms to order O(µ −2 ), the EFT Lagrangian (2.13) reads in terms of the Higgs basis states, Eqs. (2.3)-(2.6), (A.1) The couplings of the Majorana fermions to the Higgs basis states and the Z boson can be read off from here. Note, that couplings g χχΦ i Φ j ... to Higgs basis states written explicitely in this work are normalized as