The waning of the WIMP? A review of models, searches, and constraints

Weakly Interacting Massive Particles (WIMPs) are among the best-motivated dark matter candidates. No conclusive signal, despite an extensive search program that combines, often in a complementary way, direct, indirect, and collider probes, has been detected so far. This situation might change in near future due to the advent of one/multi-TON Direct Detection experiments. We thus, find it timely to provide a review of the WIMP paradigm with focus on a few models which can be probed at best by these facilities. Collider and Indirect Detection, nevertheless, will not be neglected when they represent a complementary probe.


Introduction
A combination of cosmological observations including (between others) studies of the cosmic microwave background (CMB), distant supernovae, large samples of galaxy clusters, baryon acoustic oscillation measurements has firmly established a standard cosmological model where the Dark Matter (DM), a new yet-to-be discovered form of matter, accounts for about 85% of the matter content of the Universe, and about 27% of the global energy budget [1]. Cold Dark Matter (CDM) is a key ingredient to successfully explain a e-mail: arcadi@mpi-hd.mpg.de b e-mail: maira.dutra@th.u-psud.fr c e-mail: pradipta.ghosh@th.u-psud.fr d e-mail: lindner@mpi-hd.mpg.de e e-mail: yann.mambrini@th.u-psud.fr f e-mail: mathias.pierre@th.u-psud.fr g e-mail: profumo@ucsc.edu h e-mail: queiroz@mpi-hd.mpg.de the formation of large-scale structure, producing theoretical predictions in striking agreement with observations [2,3].
Little is, however, as of yet known about DM as a particle; any candidate for (most of) the DM must nevertheless be consistent with the following five observationally-motivated constraints: (i) The relic abundance of DM needs to account for the observed CDM abundance; (ii) the DM particle should be non-relativistic at matterradiation equality to form structures in the early Universe in agreement with the observation. As a result, if the DM was produced as a thermal relic in the early Universe, its mass cannot be arbitrarily light. Specifically, cosmological simulations rule out DM masses below a few keV [4][5][6]. (iii) The DM should be electromagnetically neutral, as a result of null searches for stable charged particles [7,8] as well as Direct Detection (DD) experiments, which we will review subsequently. (iv) The DM particle must be cosmologically stable since its presence is ascertained today, implying that its lifetime is larger than the age of the Universe. Under certain assumptions, much stronger limits are applicable conservatively requiring a lifetime order of magnitude larger can be derived [9][10][11][12][13][14][15][16]. (v) Cluster collisions, such as the Bullet Cluster [17], constrain the level of self-interactions that DM particles can have (see however, Refs. [18,19] for alternative scenarios).
Within the generous parameter space outlined by the observational requirements listed above, we will argue below

(c) Collider Detection
Jets that the paradigm of WIMPs [20] is one of the most compelling options for DM as a particle. As such, it has undergone very intricate and effective experimental scrutiny. A very schematic representation of the main DM detection strategies is shown in Fig. 1.
In this work we will attempt to give an up-to-date state of the art of ongoing WIMP searches and possible future prospects.
The extreme broadness of the topic makes, however, very difficult to satisfactory cover all the different DM search strategies. Consequently we will mainly focus on DM DD, motivated by the advent of highly sensitive one-and multi-Ton detectors. We will then investigate, a selection of simple, but well motivated, WIMP models whether and to which extent the WIMP paradigm can be tested through this kind of search strategy.
Even within the WIMP framework, other DM search strategies like collider searches and Indirect Detection (ID) provide a complementary and essential contribution. This complementarity will be highlighted in some relevant cases of study.
Before discussing the main topic, we will anyway provide in the next sections a general and pedagogical introduction to the WIMP paradigm for the generation of the cosmological abundance of the DM relic density and to the three main categories of DM searches: DD, ID and collider searches.

The WIMP paradigm
The paradigm of thermal decoupling, based upon applications to cosmology of statistical mechanics and particle and nuclear physics, is enormously successful at making detailed predictions for observables in the early Universe, including the abundances of light elements and the CMB [21]. It is somewhat natural to invoke a similar paradigm to infer the abundance of DM as a thermal relic from the early Universe uniquely from the underlying DM particle properties.
Assuming there exist interactions between a cosmologically stable particle χ -the (generic) DM -with Standard Model (SM) particles, sizable enough so that for a high enough temperature T the DM is in thermal equilibrium with the primordial thermal bath, the cosmological evolution of the DM particle can be traced through the following Boltzmann equation: describing the DM number density n χ , in turn defined as: with f χ being the DM distribution function. g χ is the number of internal degrees of freedom of the DM particle. The quantity σ v , dependent on temperature T , is the thermally averaged pair annihilation cross-section associated to the process χχ → pair of SM particles while H (T ) is the Hubble rate. Assuming a standard cosmological evolution for the Early Universe, the whole process of DM production occurs while the Universe is dominated by the radiation energy density so that the Hubble expansion parameter is given by: where g eff (T ) is the effective number of relativistic degrees of freedom at the temperature T and M Pl = 1.22×10 19 GeV is the Planck mass. Moreover, n χ,eq is the equilibrium number density obtained from Eq. (2) by replacing f χ with the equilibrium distribution function (by convenience one typically adopts the Maxwell-Boltzmann distribution): where m χ is the DM mass while K i is the modified Bessel function of order 'i'. Equation (1) is more easily handled by adopting as dependent variable the 'yield' or comoving number density: being the entropy density (h eff (T ) is the effective number of entropy degrees of freedom at the temperature T ), so that it is possible to get rid of the term dependent on the Hubble expansion rate on the left-hand side of Eq. (1), giving: To obtain the last equation we have used the entropy conservation relation ds dt = −3Hs. Qualitatively, Eq. (7), describes the following picture: If DM interactions are enough efficient, as in the case of WIMPs, at early times the annihilation rate Γ ann = σ v Y χ s exceeds the Hubble expansion rate and Eq. (7) is solved for Y χ = Y χ,eq , meaning that the DM is in thermal equilibrium with the primordial thermal bath. At later times, when the temperature eventually drops below the DM mass, the DM yield becomes Boltzmann suppressed, Y χ,eq ∝ exp(−m χ /T ), so that the annihilation rate falls below the Hubble expansion rate leading to the thermal freeze-out of this "cold" relic, i.e., thereafter Y χ is approximately constant with time. 1 Equation (7) can be solved by adopting the temperature T 2 of the thermal bath or x = m χ /T as independent variable. A good approximate solution is represented by the following semi-analytical expression [23]: where: with T 0 as the present time temperature while T f represents the freeze-out temperature which can be determined by solving the equation: where δ = (Y χ − Y χ,eq )/Y χ,eq is conventionally set to 1.5 while x = m χ /T . The DM relic abundance is usually expressed in terms of the parameter Ω DM h 2 where h ∼ 0.7 is the value Hubble expansion rate at present times in units of 100 (km/s)/Mpc 1 See Ref. [22] for an exception ("relentless" DM) for modified expansion histories. 2 We make also use of: ds dt = 3s while Ω DM represents the ratio between the DM energy density ρ DM and the so called critical energy density ρ cr , namely: where s 0 = s(T 0 ) is the entropy density at present times. By combining the expressions above, the DM relic density can be numerically estimated as: The behavior of the solution of the Boltzmann equation is illustrated in Fig. 2. As expected, the DM relic density is basically set by the inverse value of the thermally averaged cross-section (calculated at the freeze-out temperature), with a logarithmic dependence on m χ . It can be straightforwardly verified that the experimental determination of Ω DM h 2 ≈ 0.12 [1] is matched by a value of the cross-section of the order of 10 −9 GeV −2 corresponding to σ v ∼ 10 −26 cm 3 s −1 .
The WIMP paradigm hence reduces, under the hypothesis of standard cosmological evolution of the Universe, the solution of the DM problem to the determination of a single particle physics input, i.e., the thermally averaged pair annihilation cross-section of the DM.
Its formal definition reads [23]: 3 where σ (s) is the DM annihilation cross-section, as computed using the conventional field theory techniques. Since WIMPs freeze out in the non-relativistic regime, and thus, v c (where v is the relative velocity of the two annihilating WIMPs), a useful approximation consists of a velocity expansion (given in the Appendix) σ v a + bv 2 . The velocity expansion is, however, not valid in some relevant cases, like for example annihilations through the resonant exchange of an s-channel mediator [25]. For this reason, all the numerical results presented in this work will rely on the full numerical determination of σ v , as given in Eq. (14) and on the solution of the DM Boltzmann equation, as provided by the numerical package micrOMEGAs [26][27][28].
The WIMP "miracle" is the observation that the cosmologically favored value of the DM pair annihilation crosssection is met by a DM featuring electroweak (EW) interactions, i.e., with cross-section scaling as σ v ∼ g 4 /m 2 χ , with g being a coupling of the order of the EW gauge couplings and mass in the O(100 − 1000) GeV range and with typical freeze-out temperature T f ∼ m χ 20 . It is important to realize that this coincidence is a statement about cross-sections (and, weakly, masses), and thus, is not unique to the weak scale and weak interactions. In large fraction of this work we will for example consider the scenario in which the DM is a SM singlet interacting (mostly) with the SM fermions through a mediator field. In this case a typically expected scaling, in the case when the DM is lighter than the mediator, for the cross-section would be σ v ∼ λ 2 χ λ 2 This shows that also in this setup, the cosmological relic density can be obtained for EW scale masses of the DM and of the mediator field, provided a suitable assignation of the couplings λ χ , λ f . It is remarkable to notice that the WIMP paradigm provides, independent of theoretical reasons such as naturalness and the hierarchy problem, a plausible argument to expect new physics at and around the EW scale.
As a result, concrete realizations of WIMP models have been developed in different Beyond the Standard Model (BSM) frameworks, accessible to several different search strategies, as reviewed in the next sections.

Direct detection
Astrophysical and gravitational evidences indicates the existence of DM halos, surrounding all the visible structures like galaxies and the cluster of galaxies. The halos are formed by DM particles described by a time dependent velocity distribution. DD experiments aim at detecting DM, through scattering off nuclei, belonging to the halo surrounding our galaxy, flowing through the Earth. Several experiments have played an important role in this direction . In this section we focus on DD experiments looking for WIMPs scattering, but there are important searches stemming from Neutrino telescopes by measuring the neutrino flux from the Sun [60][61][62][63].
Direct DM detection seeks to measure the nuclear recoil imparted by the scattering of a WIMP particle. The WIMPnuclei differential scattering rate can be written as, where E is the recoil energy associated to the scattering events, N T is the number of target nuclei per kilogram of the detector, m χ is the DM mass, ρ χ is the local DM density (ρ χ = 0.3 GeV/cm 3 ) [64][65][66][67][68], v is the velocity of the DM particle relative to the Earth, f E (v, t) is the distribution of velocities of the WIMP in the frame of the Earth, 4 v min = m N E/(2μ 2 ) is the minimum WIMP speed required to produce a detectable event at energy E, v esc is the escape velocity i.e., the velocity for which the WIMP are no longer gravitationally bounded to the Milky Way. μ = m χ m N /(m χ + m N ) is the WIMP-nucleus reduced mass (m N is the nucleus mass), dσ /d E(v, E) is the differential cross-section for the WIMP-nucleus scattering as follows, with q being the momentum transfer. F 2 (q) and S(q) are the Spin-Independent (SI) and Spin-Dependent (SD) form factors, as described e.g., in Refs. [69][70][71][72]. After measuring the scattering rate, the next and fundamental task is to discriminate the signal from the background. This is done by using the detector response to electron and nuclear scattering, which might vary from one experiment to another. For instance, in Germanium detectors ionization yield is used to discriminate signal from background, whereas in experiments that use Xenon, the ionization/scintillation ratio is the discriminating variable. What determines an experiment's sensitivity to a WIMP signal is a combination of: 3 Left Illustrative impact of the energy threshold, exposure and target nucleus. Right Impact of background and exposure on the sensitivity. These plots are taken from Ref. [73] 1. Energy threshold: drives the sensitivity to low WIMP masses, and consequently the sharpening of the DD limits on the scattering cross-section at low masses as shown in Fig. 3; 2. Control over the background and exposure: determine the overall sensitivity of the experiment pushing the limits to lower scattering cross-sections assuming that they are statistically dominated; 3. Target: has an impact on the experimental sensitivity to low and heavy WIMP masses, as well as on the capability to probe SD scatterings.
All these facts are illustrated in Fig. 3 where we plotted the impact of exposure, energy threshold and mass of the nucleus target on DD experiments sensitivity for WIMPnucleon scatterings. Figure 3 is described as follows: Left panel: (i) Comparing the solid black and blue curves at low WIMP mass one can see that the energy threshold determines the smallest WIMP mass accessible to a given DD experiment. The blue line refers to a DD experiment with lower energy threshold; (ii) Notice that stronger bounds on the scattering cross-section are possible with larger exposure as represented in the solid green line; (iii) The target nucleus can influence the WIMP mass where the strongest limit on the scattering cross-section lies at and also the sensitivity to lower and larger WIMP masses. This is visible by comparing the red and black lines. Right panel: (i) Comparing the black and blue curves we can see the importance of increasing exposure; (ii) Red and green curves exhibit the impact of background discriminations.
It is important to highlight that from going to the measured scattering rate in Eq. (16) to the derivation of a limit on the WIMP-nucleon scattering cross-section as a function of the WIMP mass, there are some assumptions that have to be made about the velocity distribution, nuclear form factor, type of WIMP-nucleon scattering, and local DM density that suffers from large uncertainties [74][75][76]. In particular, the common assumptions are that there is a smooth halo of DM particles in our galaxy well described by a Maxwellian velocity distribution [77][78][79], that the nucleus can be treated as a hard sphere as indicated by the Helm form factor [70], and that the WIMP-nucleon scattering is elastic. Our results rely on the same set of assumptions throughout this manuscript (see Refs. [80][81][82][83][84] for further discussions on these topics). Interestingly, if the uncertainties present in the astrophysical input are under control and precise measurements on the scattering cross-section can be realized, then one might even determine the nature of DM using DD experiments alone [85].
In summary, the present measurement of the scattering rate has not yet observed any excess over the background, which after some assumptions, translates into limits on the WIMP-nucleon scattering cross-section as a function of the DM mass. In this work we will be using the following limits and projections: • Current spin-independent limit: We adopt the latest result from XENON1T [86], based on 34.2 live days of exposure and a fiducial mass of 1042 ± 12 kg. The maximal sensitivity corresponds to a DM mass of 35 GeV for which SI scattering cross-sections above 7.7×10 −47 cm 2 are excluded at 90% confidence level. Despite of the limited exposure time, this result improves the previous limits provided by LUX, with strongest exclusion (90% C.L) of 2.2 × 10 −46 cm 2 at the value 50 GeV [57] of the WIMP mass based on 3.35 × 10 4 kg day exposure and, PandaX [87], whose maximal sensitivity 2.5 × 10 −46 cm 2 based on 3.3 × 10 4 kg day exposure is reached for WIMP mass of 40 GeV, hence showing the potential capability of the detector in probing the WIMP parameter space.
• Current spin-dependent limit: We adopt the latest results by LUX considering the full exposure time [88]. For the SD cross-section on neutrons, a limit as strong as 1.6 × 10 −41 cm 2 for a WIMP of mass 35 GeV is obtained. Also in this case similar limits have been reported by the PandaX collaboration, based on 3.3 × 10 4 kg day of exposure. Hopefully PandaX will continue to run and improve their sensitivity and possibly unearth a DM signal.
• Projected spin-independent limit: We adopt the projected XENON1T SI limits outlined in Ref. [89] assuming a 2 years exposure and the LZ collaboration referred as baseline in Ref. [90].
• Projected spin-dependent limit: In the lack of published projections for SD scattering we simply rescaled the current SD limits taking into account the planned exposure. In practice, we derived the scaling factor between latest LUX SI limit and the LZ projection, and then applied this same scaling factor to derive the LZ/XENON1T projection for SD scattering. In the light of no large background, the limits will roughly be improved simply by exposure, justifying our method.
The limits (current and projected) which will be adopted throughout all this work are illustrated in Fig. 4. The two panels report, in a bi-dimensional plane of DM mass and scattering cross-section, the region which are/will be excluded, under the null result, by the concerned experiments. As evident that the experimental sensitivity appears to be very different between SI and SD interactions. This is due to the fact that the SI cross-section of the DM on a nucleus originates from a coherent sum of the contributions from the interactions of the DM with the single nucleons. Targets, like Xenon, with atomic number of the order of 100 feature then an enhanced sensitivity to SI interactions. The case of SD interactions is, instead very different, since the contribution to the different nucleons tend to interfere destructively, so that a sizable cross-section is obtained only for targets with an unpaired nucleon.
As customary, we have expressed the limits on SI interactions in terms of the scattering cross-section of the DM on protons. For SD interactions we have instead reported the scattering cross-section on neutron since Xenon targets, characterizing all the three considered experiments, are mostly sensitive to this process having Xenon isotopes with unpaired neutrons.
In most of the models considered in this work the DM candidate will have unsuppressed (i.e., independent from the DM Fig. 4 Limits (by LUX) and projected sensitivities (by XENON1T and LZ) adopted in this study for SI (top panel) and SD (bottom panel) interactions in the bi-dimensional plane (DM mass, cross-section). For SI interactions we have reported, as usual, the scattering cross-section of the DM on protons; for SD we have, instead, considered scattering on neutrons since we are considering "Xenon based" detectors. The blue region are, respectively, currently excluded by LUX while the magenta/purple regions will be excluded in case of absence of signals at XENON1T/LZ velocity or the momentum transfer) SI interactions. Given the stronger sensitivity, we will consider just the limits from SI interactions in presenting our results while report the ones from SD interactions only when relevant.

Indirect detection
ID of the DM relies on the detection of the byproducts of WIMPs annihilations over the expected background at galactic or extra-galactic scales, using Earth based telescopes such as HESS and CTA, or satellites such as AMS and Fermi-LAT [91][92][93][94][95][96][97][98][99][100][101][102][103][104][105] 5 In this regard, the search for gamma-rays and cosmic-rays and neutrinos offer an exciting possibility of DM detection. Here we will focus on gamma-rays. The gamma-ray flux from WIMP annihilation is proportional to: -The number density squared of particles, i.e., n 2 χ = ρ 2 /m 2 χ ; -The WIMP annihilation cross-section today, σ ; -The mean WIMP velocity v; -Volume of the sky observed within a solid angle Ω; -Number of gamma-rays produced per annihilation at a given energy and for a given annihilation final state, also known as the energy spectrum (d N/d E).
In summary, it is found to be: Therefore, the differential gamma-ray flux in Eq. (18) is determined with the three different inputs, DM annihilation cross-section, energy spectrum, the integral over the line of sight (l.o.s) of the DM density distribution, as represented above. The DM annihilation cross-section along with the energy spectrum capture the particle physics input of the model, whereas the DM distribution accounts for the astrophysical input. At the end of the day, the most important information needed on the particle physics side is to know whether the DM annihilation cross-section depends on the relative velocity of the WIMP. If it does, then the annihilation cross-section today is bound to be suppressed, rendering ID probes sub-dominant or irrelevant. The energy spectrum plays also an important role, because different final states produce distinct gamma-ray yields. Some final states produce hard gamma-rays at high energies, while others at lower energies. For this reason, the flux of gamma-ray due to DM annihilation directly depends on the final state annihilation. The DM density profile which is not so well-known serves as a normalization quantity, since it does not feature any energy dependence, but gives an overall constant for a given target region.
In Eq. (18) the DM density is integrated over the l.o.s from the observer to the source. The DM density is not tightly constrained, and several DM density profiles have been considered in the literature leading to either spike or core DM densities toward the center of galaxies [115][116][117][118][119][120][121]. In this work we adopt the Navarro-Frenk-White (NFW) profile [117] which reads, where r s = 24.42 kpc is the scale radius of the halo, as used by Fermi-LAT collaboration in Ref. [101], and ρ s = 0.184 is a normalization constant to guarantee that the DM density at the location of the Sun is 0.3 GeV/cm 3 . It is important to emphasize that the NFW profile is known as steep profile, since it leads to a large DM density toward the inner regions of the galactic center. Therefore, it is not a conservative choice [122]. Different profiles, core-like, can significantly weaken limits from ID probes. From Eq. (18) it is clear that ID probes complementary properties of the DM particles. It is sensitive to how the DM is distributed, to the annihilation cross-section today, which might be different than the annihilation cross-section relevant for the relic density, and to the WIMP mass. Therefore, after measuring the flux of gamma-rays from a given source, we compare it with the background expectations. If no excess is observed, we can choose a DM density profile and select an annihilation final state needed for d N/d E, and then derive a limit on the ratio σ v/m 2 χ according to Eq. (18). This is the basic idea behind experimental limits. Although, more sophisticated statistical methods have been conducted such as likelihood analysis.
An interesting aspect of indirect DM detection, when it comes to probe WIMP models, is the fact that if the annihilation cross-section, σ v, is not velocity dependent, bounds on σ v today are directly connected to the DM relic density. In particular, the observation of gamma-rays from dwarf spheroidal galaxies (dSphs) results in stringent limits on the plane of annihilation cross-section vs WIMP mass. If for a given channel the annihilation cross-section of 10 −26 cm 3 s −1 is excluded for DM masses below 100 GeV, it also means that one cannot reproduce the right relic density for WIMP masses below 100 GeV. 6 In other words, in this particular case, ID limits will trace the relic density curve. This effect will be clearly visible in many instances.

Collider searches
LHC proton-proton collisions might result in the production of WIMPs in association with one or more QCD jets, photons as well as other detectable SM debris. Since WIMPs are electrically neutral and cosmologically stable massive particles, they manifest at colliders as missing transverse momentum. For this reason searches for DM are based on the observation of the visible counterpart of the event such as charged leptons, jets or a photon, generally referred to as mono-X searches. By selecting events with large missing transverse momentum/energy one can reduce the SM background and potentially disentangle a DM signal. However, as mentioned above, what colliders identify is missing energy, and therefore they cannot uniquely ascertain the presence of DM in a signal event. They can simply confirm the presence of a neutral and "stable" particle, that might have even decayed outside the detector.
Nevertheless, colliders offer an exciting and complementary search strategy to identify WIMPs. Indeed, assuming that the production of WIMPs at colliders is uniquely connected to the WIMP-nucleon scatterings at underground laboratories, one can use the non-observation signals with large missing transverse momentum to derive limits on the WIMPnucleon scattering cross-section [123][124][125][126][127][128][129][130].
A large part of this work will be devoted to the study of the cases in which the interactions between DM and SM particles are due to a neutral spin-0 or spin-1 s-channel mediator (see next sections). In this kind of scenario a compelling complementary collider probe is represented by the searches of its "visible" (i.e., into SM final states) decay channels.
We now review in some detail the specifics of the aforementioned search channels for WIMPs at colliders.

Mono-X searches
Mono-X searches stand for the search of WIMPs produced in association with one or more QCD jets or potentially other SM particles, such as γ , h, Z etc. The idea is to search for events with a jet of high transverse momentum p T within an event with large missing transverse momentum. In particular, the most recent studies performed at the LHC include up to four jets and require the leading jet to have p T > 250 GeV [131,132], while others do not limit the number of jets while selecting events with at least one jet with p T > 100 GeV [133]. While being more inclusive, these recent searches have become more challenging due to the number of jets analyzed, requiring a substantial improvement on the background coming from Z + jet and W + jet channels.
There are important detector effects, such as fake jets, and QCD backgrounds that weaken the LHC sensitivity to WIMPs, and for these reasons mono-jet searches are subject to large systematics. Nevertheless, fortunately an enormous effort has been put forth in this direction with data driven background and optimized event selections, which combined with the increase in luminosity has led to an overall improvement on the LHC sensitivity to WIMPs.
That said, in the review, we will be using the latest results from CMS and ATLAS collaborations in the search for DM based on mono-X searches [133,134]. Now we discuss the WIMP production at colliders and also address another collider constraint relevant for the DM purposes which has to do with the invisible decay widths of the SM-Higgs as well as Z bosons.

Invisible Higgs decays
If WIMPs are lighter than 62.5 GeV, the Higgs boson might invisibly decay into WIMP pairs. In this case, one can use bounds from LHC on the invisible branching ratio (Br) of the Higgs, Br(h → inv) ≤ 0.25 at 95% C.L. [135,136], to set constraints on WIMP models . Throughout the manuscript whenever applicable we compute the invisible decay rate of the Higgs into WIMPs and impose the aforesaid upper limits to obtain the limits displayed in the figures.

Invisible Z decays
The decay width of the Z boson has been precisely measured and therefore stringent limits can be derived on any extra possible decay mode of the Z boson. In some of the models as we will discuss further, the DM particle does couple to the Z boson, thus, when mass of the DM is smaller than half of the Z mass stringent limits are applicable. In particular, one can use only direct measurements of the invisible partial width using the single photon channel to obtain an average bound which is derived by computing the difference between the total and the observed partial widths assuming lepton universality. The current limit is Γ (Z → inv) ≤ 499 ± 1.5 MeV [137].

Searches of visible decays of the mediator
As already mentioned a relevant case of study of this review is represented by the case in which a new neutral field mediates interactions between SM fermion pairs and DM pairs. Provided that its mass is within the reach of a collider, it can be singly produced thanks to its coupling with the SM fermion pairs. The mono-X searches discussed before essentially probe the decay channels of the mediator into DM pairs (the mono-X is radiated by one of the initial state fermions). A potentially even more powerful probe is, however, provided by the visible decay channels of the BSM mediator. Its on-shell production, and subsequent decay into SM states, may lead to spectacular signals represented by dijet and/or dilepton resonances 7 peaked at the mass of the mediator field. Among the models which will be discussed here these kinds of signal are particularly prominent in the case of spin-1 mediators, since their gauge-like couplings with the SM fermions allow a high production cross-section for these resonances. The considered models with spin-0 mediators feature instead poorer prospects since a Yukawa-type structure will be assumed for their couplings with the SM fermions.
We will consequently specialize our discussion to the case of spin-1 resonances.
Among dijet and dilepton resonance searches, the latter have typically the potential to exclude larger portions of the parameters space. This can be understood by noticing the different backgrounds, QCD jets and Drell-Yan production for dijet and dileptons respectively. Dilepton resonances typically represent a much clearer signal, hence provide the strongest constraints unless the couplings of the spin-1 mediator with SM leptons are very suppressed, with respect to the ones with the SM quarks, or even null.
In this work we will mostly consider models with spin-1 mediators with couplings with the SM quarks and leptons of similar size. We have then applied to them the LHC limits from searches of narrow resonances decaying into dileptons, considering an integrated luminosity of 37 fb −1 and 13 TeV of center-of-mass energy [138]. Under the assumption, considered in [138], that the new spin-1 neutral resonance is coupled only with the SM states, the present absence of signals is translated, according the assignations of the couplings of the spin-1 mediator, into very stringent lower bound on its mass. Masses below 4 TeV are, for example, excluded in the case when the couplings of the mediator with quarks and leptons are of the same size of the ones of the Z -boson. This picture might be potentially altered by the presence of the DM. In case the decay of the spin one mediator into a DM pair is kinematically allowed, the actual limit follows the simple rescaling relation given by [139]: where Br DM is the decay branching fraction of the neutral spin-1 resonance into DM pairs and σ ex p ll is the limit obtained from experimental searches. It is then evident that the limit on the mass of the spin-1 mediator from the searches of its visible decay channels is actually also dependent on the DM mass and couplings.
Keeping a similar spirit we will discuss a selection of DM setups, with increasing degree of refinement, which we will globally refer to as "Dark Portals". These are characterized by spin-0, spin-1/2 and spin-1 DM candidates (we will distinguish whenever relevant the cases of self-and not selfconjugated DM) interacting with the SM fields, typically the electrically neutral bosons, and/or possibly a spin-0 or spin-1 mediator field. 8 A schematic representation of the Dark Portal models considered in this work is provided by the following Lagrangians (for simplicity we will report only the interaction terms): We will describe in detail the different Lagrangians in dedicated subsections. A brief overview of some relevant details will nevertheless be provided at the end of this section.
We have generically labeled a (real) scalar 9 and vector mediators as S and Z , respectively. We have indicated a scalar, fermionic and vectorial DM as χ , ψ and V , respectively while f generically refers to a SM fermion. Throughout this work we refer to models described by Lagrangians of the type as schematically shown in Eqs. (21) and (22) as "s-channel portals" since the DM annihilation processes into the SM states are described by diagrams in which the mediator is exchanged in the s-channel. We will also assume that the s-channel mediators interact with all the SM fermions (sums over all the SM fermion species are, hence, implicitly 8 We will briefly consider a single case of spin-1/2 mediator. We won't consider higher spin assignations. For these cases the interested reader can look for example at Refs. [164,167]. 9 We will also consider the case of a pseudoscalar mediator. As will be clarified later, this is a more specific scenario since, under the main assumptions considered in this work, it can occur only for a fermionic DM. We refer to the dedicated section later for the corresponding Lagrangian. intended in Eqs. (21) and (22)). In addition to these models we will consider "t-channel portals", schematically described by: with f and Ψ f being, respectively, a scalar and a fermionic field with the same quantum numbers as of the SM fermion f . Here annihilations into the SM states are associated with the exchange of mediator field in the t/u-channel. Contrary to the case of s-channel mediators, which can be coupled with all the SM fermions, t-channel mediators, according their quantum number, are coupled only to one type of the SM fermions, i.e., left/right-handed leptons, left/right-handed uptype quarks or left/right-handed down-type quarks. Some further comments about Lagrangians of Eqs. (21) and (22) are discussed subsequently. The trilinear coupling between a spin-0 mediator and spin-0 (spin-1) DM is dimension-full. This has been re-expressed as the product of a dimensionless quantity λ S χ (η S V ) and a mass scale μ S χ (μ S V ). We will assign specific value, according to the considered theoretical framework, for the latter while regards λ S χ (η S V ) as free parameter. The coupling of the scalar mediator with the SM fermions has been written as the SM Yukawa coupling, i.e., m f / √ 2v h , times a free scaling factor c S . The first two Lagrangians in Eq. (21) and the second Lagrangian of Eq. (22) are substantially equivalent for self-conjugate (i.e., real scalar or Majorana fermion) and non self-conjugate (i.e., complex scalar and/or Dirac fermion) DM with ξ = 1/2 and 1, respectively. For this reason we are explicitly reporting Lagrangians only for the case of a real scalar and a Dirac fermionic DM. Note, however, that V Z ψ is identically zero for a Majorana fermionic DM. The first Lagrangian of Eq. (22) can only be written for a complex scalar DM. In the case of a spin-1 DM with spin-1 mediator, two sensitively different Lorentz structures are allowed for non-self conjugate (third Lagrangian of Eq. (22)) and self-conjugate DM (fourth Lagrangian of Eq. (22)) which we have reported individually. In the case of a non self-conjugate vectorial DM the interactions are substantially analogous to the one of two SM W bosons with the Z boson. Thus, and for this reason we will often refer to this scenario as non-Abelian DM. We, however, will use Abelian DM for the case of a self-conjugate vectorial DM coupled through the Levi-Civita operator to the spin-1 mediator. We finally remark that in all the cases when the DM is interacting with a spin-1 mediator, we have normalized the couplings with the SM gauge coupling g. The Lagrangians of Eqs. (21)-(23) describe interactions of the DM candidates which are singlet under the SM gauge group. In the last part of this work we will, instead, relax this assumption and consider a DM (at least partially) charged under EW interactions. These kind of models have a richer phenomenology and are better motivated from the theoretical point of view though they still feature a limited number of free parameters. In addition they are still somehow correlated to the Dark Portal models and Eq. (21), Eq. (22) represent a possible UV completion for some of them.
In summary we will discuss a collection of models mostly built according to a bottom-up approach and maintains a substantial degree of simplicity. Nevertheless, contrary to the aforementioned "simplified models", they also aim at a more faithful description of theoretically motivated scenarios. For example, in most cases we will adopt specific choices of the parameters inspired by possible theoretical completions of the models under consideration. In the same fashion we will account, in our analysis, besides experimental constraints, the theoretical limitations of these frameworks (see e.g., Refs. [151,152,161,162,168] for more extensive discussions.) Besides these considerations the choice of the models under study is dictated by the main purpose, already stated in the introduction, of assessing the capability of present and, more importantly, near future DD facilities of probing the WIMP paradigm. For each of the models under scrutiny we will determine the portion of the parameter space excluded by current DD limits, as set by the XENON1T experiment, and the parameter space which would be ruled out in case of absence of signals, after 2 years of exposure time at the XENON1T experiment itself, and at the multi-TON detector LZ. Most of the considered models will then be characterized by efficient interactions of the DM with nucleons. At the same time we will discuss possible scenarios capable of surviving even eventual limits from high sensitivity future experiments like LZ. Limits from other kind of DM searches, i.e., ID and LHC, will not be neglected when competitive or complementary with the ones from DD.
To make simpler the orientation among the many different models which will be individually discusses in the following sections we have adopted the following classification: 1. SM portals: These represent a special case of s-channel portals, i.e., Eqs. (21) and (22), in which S = h and Z = Z with h and Z being the SM Higgs and Z bosons, respectively. 10 These are the most predictive models since the couplings of the mediators with the SM states are fixed. Moreover, μ h χ = μ h V = v h . These leave us with just two free-parameters, the DM mass and an adimensional coupling. On the other hand, these models are also the most strongly constrained ones so that they are already strongly disfavored by the present limits; 2. BSM s-channel portals: These are the models fully described by Eqs. (21) and (22). Here the DM is coupled with the SM fermions either by a new spin-0 (real) or by a spin-1 electrically neutral state. These scenarios, besides the SM portals, are the most sensitive ones to DD of the DM; 3. Portals evading DD: Here we will instead consider the case of a pseudoscalar s-channel mediator, in which constraints from DD are particularly weak so that the complementarity with other search strategies becomes crucial. We will also considered a more theoretically refined scenario in which the pseudoscalar mediator is part of a complex scalar field. The presence of an additional scalar component will reintroduce DD bounds. A broad region of the parameter space for thermal DM is nevertheless reopened by considering a very light pseudoscalar, which can be interpreted as a pseudo-Goldostone boson of a global symmetry carried by the original complex scalar field; 4. t-channel portals: In this alternative version of the Dark Portals the mediator field has non-trivial quantum numbers with respect to the SM gauge group. In these kind of setups DM annihilation arise from t-channel interactions while s-channel interactions are responsible for DD; 5. Portals to secluded sectors: We assume here that the mediator field cannot be directly coupled with the SM fermions. Even in these kind of constructions a portal can be originated by mass mixing with the SM Higgs, in the case of spin-0 mediators, and by kinetic mixing with the Z-boson in the case of a spin-1 field, so that the DM actually interacts through a double s-channel mediator, represented by a SM and a BSM state; 6. Portals with DM (partially) charged under SU (2)×U (1): In these models the DM is charged under the EW symmetry group, being either the lightest electrically neutral component of a (or more) SU (2) multiplet (we will discuss in detail only the case of SU (2) doublets) or a mixed state between the latter and a SM singlet. Note that this category of models also contain renormalizable completions of the Higgs and Z-portal models although the actual phenomenology features sensitive differences.
We also remark that we have assumed, throughout all our study, real couplings and hence, no CP-violation.
The SM portals are the simplest models since they are characterized by only two free parameters, the DM mass and couplings. As a consequence the comparison between DM relic density and DD limits/prospects (and eventual additional constraints) can be performed in full generality in a bi-dimensional plane of the two quantities. s-and t-channel portals, as appear for the aforementioned model categories 2, 3 and 4, feature two mass scales, the DM mass and the mediator, and one or at most two relevant couplings (in some cases this actually implies additional assumptions, clarifications will be provided in the dedicated sections). For these scenarios we have adopted the two masses as freely varying parameters and assigned O(1) values to the couplings unless this option is precluded by theoretical considerations.
There is very modest loss of generality in this choice. The DM annihilation and scattering cross-section feature a very similar dependence on the free couplings, lower values of the couplings correspond to weaker limit/prospect from DD but, at the same time, progressively disfavor the achievement of the correct relic density so that the overall picture is only marginally affected. For the last category of models the interplay between DD and relic density is less trivial. In this case we have considered, besides representations in bidimensional planes of two relevant parameters, also scans in which all the free parameters of the models are varied.
Concerning, in particular, the DM mass, we have typically assumed a range of variation from 10 GeV to a few TeV. The choice of the lower bound is motivated by our focus on DD as DM search strategy. For masses below 10 GeV, experimental sensitivity is strongly reduced because of the approaching of the energy threshold of the experiments under consideration. We also remark that achieving the correct relic density becomes increasingly more difficult as the DM mass decreases, because of the reduced number of accessible final states (for this reason in the case of DM interactions mediated by BSM fields the lowest value of the DM mass shown in the figures is 100 GeV). Finally for several realizations, namely the cases in which the annihilation cross-section is velocity independent (s-wave), thermal DM below 10 GeV is already ruled out by ID as well as CMB measurements [1].
Concerning the maximal value of the DM mass, i.e., a few TeV, this should be regarded mostly as an assumption. We remark, nevertheless, that in bottom-up setups it might result difficult to accommodate viable thermal DM without encountering issues of theoretical consistency. In other cases, like the ones discussed in Sect. 12, for a high value of the DM mass the correct relic density cannot be achieved because of a systematic suppression of the annihilation cross-section. Relevant viable WIMP scenarios, featuring DM heavier than a few TeV, like the so-called "minimal DM", have been proposed in Refs. [169][170][171]. These kind of models will not be explicitly revised here; we will nevertheless discussed them in slightly more detail in Sect. 12.

SM portals
The first class of models which will be the object of study are the SM Dark Portals, 11 i.e., models in which the DM interacts with the SM state through the Higgs or the Z-boson. In the case when the DM is a pure SM singlet, gauge invariant renormalizable operator connecting the DM with the Z or the SM-Higgs boson can be build only in the latter case and only for scalar and vectorial DM. In the other cases one should rely either on higher dimensional operators, or on the case that the coupling with the Higgs and/or the Z is originated by their mixing with new neutral mediators. The latter case can imply the presence of additional states relevant for the DM phenomenology and will be then discussed later on in the text. We will instead quote below some example of higher dimensional operator but we will not refer to any specific construction for our analysis. Alternatively one could assume that the DM has some small charge under SU (2) L or U (1) Y , see e.g., Refs. [173][174][175][176][177][178][179]. We will not review these scenarios here.

Higgs portal
The most economical way to connect a SM singlet DM candidate with the SM Higgs doublet H is through four field operators built to connect the Higgs bilinear H † H , which is a Lorentz and gauge invariant quantity, with a DM bilinear. Assuming CP conservation, the possible 12 operators connecting the SM-Higgs doublet with scalar, fermion and vector DM are given by [181][182][183][184][185][186][187][188][189][190][191]: where, in the unitary gauge, with h, v h denoting the physical SM Higgs boson, Vacuum Expectation Value (VEV) and ξ = 1/2(1) in case the DM is (not) its own antiparticle. From Eq. (25) note that stability of the DM is protected either by a discrete Z 2 (for ψ, V μ and when χ = χ * ) or by a U (1) (for χ = χ * ) symmetry.
As already pointed out, in the case of a scalar and vectorial DM it is possible to rely on a dimension-4 renormalizable operator; on the contrary a fermionic DM requires at least a dimension-5 operator which depends on an unknown Ultra-Violet (UV) scale Λ.
After EW symmetry breaking (EWSB), trilinear couplings between the Higgs field h and DM pairs are induced. In the case of fermionic DM it is possible to absorb the explicit Λ dependence by a redefinition of the associated coupling, i.e., λ H ψ v h Λ as λ H ψ , so that it does not appear explicitly in computations.
The models defined by Lagrangians of Eq. (25) have only two free parameters, the DM masses m χ,ψ,V and couplings λ H χ,ψ,V with the SM-Higgs. The constraints on these models can be then easily summarized in bi-dimensional (i.e., DM mass vs its coupling with the SH-Higgs) planes.
In Fig. 5 we summarize our results for scalar, fermionic and vectorial DM, respectively. All the plots report basically three set of constraints. 13 The first one (red contours) is represented by the achievement of the correct DM relic density. 14 The DM annihilates into SM fermions and gauge bosons, through s-channel exchange of the SM-Higgs boson, and, for higher masses, also into Higgs pairs through both sand t-channel diagrams (in this last case a DM particle is exchanged). Since the coupling of the SM-Higgs with SM fermions and gauge bosons depends on the masses of the particles themselves, the DM annihilation cross-section is suppressed, at the exception of the pole region m χ ∼ m h /2, until the W + W − , Z Z and tt final states are kinematically accessible. Even in this last case, the cosmologically allowed values for the couplings are in strong tension with the constraints from DD of the DM, which for all the considered spin assignation of the DM, arise from SI interactions of the DM with the SM quarks originated by t-channel exchange of the SM-Higgs boson. As can be easily seen that the entire parameter space corresponding to thermal DM is already ruled out, at the exception, possibly, of the pole region, for DM masses at least below 1 TeV. Eventual surviving resonance regions will be ruled-out in case of absence of signals at the forthcoming XENON1T, assuming a 2 years of exposure time. As expected, the most constrained scenario is the fermionic DM one because of the further suppression of the p-wave suppression of its annihilation cross-section.

Z-portal
An interaction between the Z -boson and a SM singlet DM candidate is not gauge invariant for any dimension-4 operators. In the case of scalar and fermionic DM models, the simplest option is to consider a dimension 6 operator 15 [178,[206][207][208]. In the case of scalar DM it is of the form: which give rise to a trilinear interaction between the Z -boson and a DM pair once the SM-Higgs field in the Lagrangian is replaced by its VEV, so that H Λ is again the relevant cutoff scale of the effective theory. Similar to the case of a fermionic DM in Higgs portal we can absorb it in the definition of a dimensionless coupling as In addition, after EWSB, an effective dimension-4 inter- For simplicity we maintain a rescaling with powers of the SU (2) L gauge coupling g.
The interaction Lagrangian for the DM, along with the relevant SM parts, can thus be written as: 15 Similar to the case of the Higgs portal we just quote, as an example, the lowest dimensional operator. This is however, not the only possible option.
where c W = cos θ W and θ W is Weinberg angle [137,209]. Note that we have used a normalization of g/4 cos θ W throughout in analogy to the SM f f Z couplings. The interaction Lagrangian for fermion DM is built in a similar fashion as the scalar case. In the case of Dirac DM the starting operator is: which, after the EWSB, together with the apposite SM part leads to: In the case of Majorana DM V Z ψ = 0 and we rescale the remaining DM coupling by a factor of 1/2.
In the case of spin-1 DM we will consider two possible kind of interactions, namely, self-(Abelian) and not selfconjugated (non-Abelian) DM, respectively. For the latter, along with the necessary SM parts, we can write the following Lorentz invariant interaction: coupling is normalized as g/4 cos θ W while the model specific information are parametrized as η Z V . In the case of self-conjugate spin-1 DM, an interaction with the gauge boson can be built through the Levi-Civita symbol as Ref. [210]: Similar to the previous cases the coupling η Z V in Eq. (31) encodes a cut-off scale (see e.g., Refs. [211][212][213][214] for the construction of effective theories). The theoretical derivation of Eq. (30) is however, more contrived.
Further, similar to the Higgs portal, the Z-portal models are fully defined by two parameters so that one can repeat the same kind of analysis performed in the previous subsection. The results are summarized in Figs. 6, 7 and 8. 16 As evident, in all but the Majorana Z -portal case, thermal DM is already excluded, even for masses above the TeV scale, by current constraints from XENON1T. These constraints are even stronger with respect to the case of the SM-Higgs portal. This is because, apart from the lighter mediator, the scattering cross-section on Xenon nuclei is enhanced by the isospin violating interactions of the Z with light quarks. Low DM masses, possibly out of the reach of DD experiments, are instead excluded by the limit on the invisible decay width of the Z -boson. As already pointed out, the only exception to this picture is represented by the case of Majorana DM where the SI component of the DM scattering cross-section is largely suppressed due to the absence of a vectorial coupling of the DM with the Z . This scenario is nevertheless already (partially) within the reach of current searches for a SD component of the scattering cross-section. The increased sensitivity of XENON1T will allow to exclude DM masses below 300 GeV, except the "pole" region. The latter, however, will meet the same fate from a projected future LZ sensitivity.

BSM s-channel portals
The results presented in the previous cases for the SM-Higgs and Z -boson portals will be generalized and discussed in more details in the case of generic, BSM spin-0 and spin-1 mediators interacting with a pair of scalar, fermion or vector DM fields. Contrary to the case of SM portals, interactions of the mediators with the gauge bosons are not mandatory. We will thus stick, in this section to the case, analogous to the socalled simplified models, in which the DM is coupled only to the SM fermions. The case of interactions with gauge bosons will be discussed separately later in the text, where unlike the aforementioned models, we will also assume interactions with both quarks and leptons.

Scalar dark matter
We will consider the following Lagrangian: where S is a real scalar field and ξ denotes the normalization factor, accounting, similar to the previous section, for the case when the DM coincides (or not) with its own antiparticle. In the case of SM fermions we have assumed a Yukawa-like structure of the couplings with the mediator while for the scalar DM (χ ) we have parametrized all the information, including possible normalization factors (e.g., factor of 1/2 in the second term of Eq. (32)), in the respective couplings. Note that μ S χ parameter has the dimension of mass. Unless differently stated we will assume μ S χ = λ S χ m S with λ S χ being a dimensionless coupling and m S as the mass of S. We will also add self interaction term for the scalar field given by: The assignation for the dimensional couplings, as well as the introduction of the Lagrangian term in Eq. (33), are inspired to scenarios in which the scalar field S acquires a VEV. In this setup the Lagrangian of Eq. (33) originates from the quartic term in the scalar potential whose presence cannot be forbidden by any symmetry argument. In the same fashion a quartic interaction term S 2 H † H with the SM-Higgs doublet, responsible for a mixing of the S and h states, should also be included. For simplicity we will assume here that the coupling of this last operator is negligible and postpone the discussion of the most general case to a dedicated section.
Contrary to the case of the SM portals, which have only the DM mass and its coupling as free parameters, we have expressed, as reported on Fig. 9, our main results in the bidimensional plane (m χ , m S ) for the three free coupling assig- Figure 9 hence shows the comparison between current DD limits, as well as the projected sensitivities from XENON1T and LZ, and the requirement of the correct DM relic density.
The results reported in Fig. 9 can be explained as follows. A t-channel exchange of the scalar mediator induces SI interactions of the DM, which are written, in the case of the proton as: Here A, Z represent, respectively, the atomic and proton number of the material constituting the detector, μ χ p = m χ m p /(m χ + m p ) denotes reduced mass of the WIMPproton system with m p representing the mass of the latter while f p and f n represent the effective couplings of the DM with protons and neutrons. In the case of a scalar mediator we have: with f N ( f N q ) being the form factor whose physical meaning is associated to the contribution from all the six quark flavours (up, down and strange quark flavours) to the mass of the proton and the neutron. The contribution of the heavy quarks in f N is described by a unique form factor (35) we have implicitly summed over the three heavy quark flavours). For their numerical values we have adopted the default assignations by micrOMEGAs. Notice that the factor f p Z A + f n 1 − Z A is actually a rescaling factor which is introduced for a consistent comparison with the experimental limits which customarily assume f p = f n [215]. This assumption is justified in the case of the spin-0 mediator since f p and f n differ only by the contributions of up and down quarks which are sub-dominant with respect to the contribution from the strange quark (and then from the heavy quarks, see Eq. (35)), which is the same for proton and neutron. In the following numerical estimates we will then automatically set f p As will be shown in the next section, for spin-1 mediators one expects in general f p = f n . This often translates into an enhancement of the cross-section and, hence, stronger limits on the model parameters. This can be already noticed by comparing the limits in the case of the SM-Higgs and Z -boson portals.
Current limits exclude then low values for both the mass of the DM and the one of the mediator. These limits will become, of course, progressively stronger, in case of absence of signals at XENON1T and/or LZ.
Concerning the DM relic density for m χ < m S , m t the DM annihilation cross-section is suppressed by Yukawa structure of the couplings so that the correct relic density is obtained only around the resonance region m χ ∼ m S /2. This corresponds to the wide region between the two lines in Fig. 9. This region is of course determined by λ S χ × c S and is then the same as shown in the left and right plots of Fig. 9. The difference between these two figures is the disappearance of the region corresponding to the SS final state (see right plot) where, annihilation cross-section is proportional to (λ S χ ) 4 (see Eq. (37)). This channel is then not sufficient to avoid an over-density of the Universe for λ S χ = 0.25. In the middle plot, the pole region is enlarged, to the point of covering almost all the parameter space, even joining the SS final state at m S 1 TeV. For higher DM masses, instead, the correct relic density is achieved also far from s-channel resonances through either the tt channel or the SS channel, whether kinematically open. In such a case, the following analytical estimates, through the conventional velocity expansion, of the DM annihilation cross-section, can be obtained: and: As evident that both the tt and SS cross-sections are swave dominated and thus, velocity independent. As a con- Fig. 10 The same as Fig. 9 but for a Dirac fermion DM, i.e., replacing λ S χ , m χ by g ψ , m ψ , respectively sequence residual annihilation would occur at present times which can be probed by DM ID strategies. Similar to the case of the SM-Higgs portal, DD limits are much more competitive with respect to the ones from ID, hence the latter have not been explicitly exhibited on Fig. 9. We also notice that the dominant contribution of the annihilation cross-section into SS depends only on the λ S χ coupling; as a consequence the scalar self-coupling λ S does not play a relevant role for DM phenomenology. 17

Fermionic dark matter
The interaction of a fermionic DM and a scalar s-channel mediator can be described by the following phenomenological Lagrangian: where L S is defined in Eq. (33). Contrary to the case of a scalar DM, the operator ψψ S is of dimension 4, so that g ψ is already a dimensionless parameter. Note that similar to Eq. (32) we have parametrized g ψ to contain all the information of the ψψ S vertex including a normalization factor. One could think that an eventual VEV of the scalar mediator S can be the origin of the DM mass, so that g ψ ∼ m ψ /v S , with v S being the VEV of S [216][217][218][219]. We won't make this assumption in this work and regard g ψ as a generic dimensionless constant.
The main results of our analysis have been summarized in Fig. 10. We have once again considered the DM and scalar 17 This statement is strictly valid in the case, considered here, of λ S χ = 1. In the case λ S χ λ S the contribution of the trilinear couplings λ S to the annihilation cross-section into SS might be sizable and event dominant. However, we expect in such a case, since the annihilation cross-section would scale at least as (λ S χ ) 2 , the DM to be in general overabundant. masses as free parameters and an analogous assignation of the couplings as in the previous subsection.
The results shown in the figure can be described analytically as follows: the DD of the DM is again principally determined by SI interactions whose cross-section is given by: where μ ψ p = m ψ m p /(m ψ + m p ) denotes reduced mass of the associated WIMP-proton system. As evidenced from Fig. 10, DM masses even above the TeV scale, are excluded by current DD limits for m S 400− 500 GeV. Values below the TeV scale for both the DM and mediator masses will be excluded in the absence of signals from the next generation experiments.
The correct DM relic density can be achieved, without relying on s-channel resonances, only when at least one between the tt and SS final states is kinematically accessible. In such a case the DM pair annihilation cross-section can be approximated as: Fig. 11 The same as Fig. 9 for a vector DM with scalar mediator, i.e., trading λ S Here v 2 ∼ 0.23. We notice again that in the limit m ψ m S , the scalar self-coupling λ S does not influence the DM relic density. The dependence on the couplings between the three plots of Fig. 10 is the same as of the scalar case. However contrary to the case of scalar DM, all the annihilation channels are now velocity suppressed, hence cannot account for a sizable ID signals.

Vector dark matter
For the description of the vectorial DM case we consider the following Lagrangian: which is inspired by the construction proposed in Refs. [220,221] (see also [222]). Note that the first three terms of Eq. (42) appear after the spontaneous symmetry breaking, once the portal field is expanded as (S + v S )/ √ 2 with v S as the concerned VEV. The quantity m V is expressed as η S V v S /2. A similar construction is also possible from a gauge invariant D μ S(D μ S) * operator for a complex scalar field S with This scenario has been analyzed with the same procedure as the scalar and fermionic DM cases. The results, reported in Fig. 11 appear not to be very different from what obtained in the case of a scalar DM in Fig. 9. This can be explained by the fact that a massive vectorial DM can be viewed as three scalar degrees of freedom. The DM scattering rate on protons and its most relevant annihilation channels are described by the following analytical expressions: The parameter μ V p = m V m p /(m V + m p ) as usual represents reduced mass of the relevant WIMP-proton system and In this subsection we will analyze, in analogous fashion as the previous subsection, the case of a s-channel spin-1 mediator. Unlike the scalar case, this kind of scenario offers a much richer collider phenomenology since one could assume gauge-like interactions (contrary to Yukawa-like interactions for spin-0 mediators) of the mediator with the SM light quarks and leptons, leading to visible signals, implying stronger constraints [223][224][225][226][227][228][229][230][231][232][233][234][235][236][237][238][239][240][241]. As a consequence we will consider a wider mass range, for both DM and the spin-1 mediator, in our analysis. New spin-1 s-channel mediators can be straightforwardly associated to gauge bosons of extra U (1) groups. Extra U (1) symmetries are particularly common in extensions of the SM and, in particular, in Grand Unified Theories (GUT). The new Z particles, i.e., the gauge bosons of these U (1) groups, can be coupled to the SM fermions either indirectly, through kinetic mixing [242][243][244][245][246][247][248][249] with the Z -boson, or directly in case the SM fermions have non-trivial charges under the new symmetry group (see e.g., Refs. [250][251][252][253][254][255][256]). We will focus, in this section, on this last case while the case of kinetic mixing will be reviewed later in the text. Scalar or fermionic DM can be easily embedded in this kind of construction since they can be assumed to be new states charged under the new U (1) but singlet with respect to the SM group. In such a case portal interactions simply arise from their covariant derivatives. 18 In the case of spin-1 DM we re-propose the two constructions: the self-conjugate (Abelian) and not self-conjugate (non-Abelian) DM as already proposed in the SM Z-portal setup.

Scalar dark matter
Following the discussion above the interaction between a scalar DM and spin-1 (Z ) mediator, together with a piece connecting Z to the SM fermions, is described by the following Lagrangian: 18 Interestingly, couplings of the Z with SM particles as presented here are similar to the ones in several existing EW extensions of the SM which can be embedded in different GUT models [257][258][259][260].
Notice that trilinear interaction between DM pairs and the Z , of the form reported above, is possible only in the case of a complex scalar DM.
Similar to the case of scalar mediator, our main parameters will be represented by the DM and Z masses. For what regards the couplings of the Z with the SM fermions we will consider some definite assignations, as dictated by the Sequential Standard Model (SSM), i.e., same couplings as the Z -boson, and some GUT-inspired realizations. According to this the coupling g will be set to g ≈ 0.65 in the case of SSM and to g GUT = √ 5/3 g tan θ W ≈ 0.46 for the GUT realizations. Finally, unless differently stated, we will (46)) considered in our analysis for the three cases of SSM, E 6 χ and E 6 ψ realizations are exhibited in Table 1. For the same three realizations the effect of different constraints are summarized in Fig. 12.
As first thing we notice, in the case of SSM and E 6χ models, a much stronger impact of the limits from direct DM searches with respect to the case of scalar mediator. The reason lies on the fact that SI interactions, with cross-section given by (as usual for the case of SI we will refer to scattering on protons): are particularly efficient since, as evident from the fact that, for spin-1 mediators, the effective couplings of the DM with the proton and the neutron, f p and f n , are just linear combination of couplings of the Z with up and down quarks. 19 A further enhancement comes, in general, as already remarked, from the fact that f p = f n . As a consequence, an absence of signal from XENON1T would exclude values of the masses of the DM and of the Z even above 5 TeV.  Table 1), i.e., SSM (left), E 6χ (middle) and E 6ψ (right) in the relevant bi-dimensional plane m Z , m χ for the choice of λ Z χ = 1. In these plots the red contours represent the correct DM relic density. The blue region is already excluded by DM DD while the magenta and purple regions are currently allowed but within the sensitivity of XENON1T (assuming 2 years of exposure) and LZ, respectively. Finally, the green and orange regions represent the exclusions from dilepton searches by the LEP/Tevatron and LHC experiments Stringent limits from DD, although weaker with respect to the previous two cases, are remarkably present also for the E 6ψ realization, despite the assignations of the charges of the quarks under the new U (1) imply a null vectorial combination. Indeed non-null vectorial couplings, at the typical energy scale of DM scattering with nucleons, are radiatively generated by the axial couplings of the Z , in particular with the top quark [261][262][263]. An approximate expression, mostly valid for m Z > m Z , for this Renormalization Group induced couplings are given by: Here with y f as the SM Yukawa couplings and μ N is the characteristic energy scale of scattering interaction, here taken to be 1 GeV.
The DM annihilation cross-section into SM fermions is instead velocity suppressed: where the sum runs over the kinematically accessible final states. The parameter n f c represents color factor for the final state fermions. The correct relic density is thus achieved only around the pole region, namely m χ ∼ m Z /2, unless the annihilation into Z Z is kinematically accessible. This crosssection is s-wave dominated and can be simply approximated, for m χ m Z as: where, for definiteness, we have considered the SSM for the numerical estimates. As already pointed out that strong collider limits complement the ones from DD. The plots of Fig. 12, compared to the previous figures, report two new exclusion regions, green and orange. Both of these regions are related to limits associated to the couplings of the Z with SM leptons. The first ones, associated to the green regions, come from the LEP and Tevatron [264,265] and are based  on possible modifications of the dilepton production crosssection. Since they do not necessarily rely on on-shell production of the Z , once its couplings with the SM fermions are fixed, like in our case, they are straightforwardly translated into lower bounds on m Z . More specifically, these lower bounds are 1789, 853 and 804 GeV for the SSM, E 6 χ and E 6 ψ configurations, respectively. These constraints are combined with the limits (orange regions in the plots) from LHC searches of dilepton resonances [138,266,267]. Contrary to the previous case, these limits are in principle sensitive to modification of the decay branching fraction of the Z as a consequence, for example, of couplings with the DM [139]. For the chosen assignation of the couplings, the decay branching fraction of the Z into DM pairs is small so that the limits substantially coincide with the ones reported by experimental collaborations. Other limits stemming from flavour and muon (g −2) are weaker compared to the collider bounds [268][269][270][271].

Fermionic dark matter
As already mentioned we will describe the interactions of a fermionic DM ψ and a Z , mediating its interactions with the SM fermions, through a Lagrangian of the form: where ξ = 1(1/2) for Dirac (Majorana) fermions. We remind that in the case of Majorana fermions V Z ψ = 0. The combination of constraints is reported in Fig. 13 for a Dirac fermion DM, following the color coding of Fig. 12  (Fig. 14).
The SI cross-section, from t-channel exchange of a Z , in the case of Dirac fermion DM, exactly coincides with the similar one for a complex scalar DM. As a consequence the excluded regions in the (m ψ , m Z ) plane are the same as shown in the previous subsection.
In the case of Majorana DM, DD principally relies on SD interactions, to which Xenon based detectors are also sensitive. We have then reported, together with the most recent constraints [53,272], an estimation of the XENON1T and LZ sensitivities. As evident, even in the case of LZ, we have much weaker limits, not competitive with bounds from dilepton searches.
On the contrary the regions corresponding to the correct DM relic density are sensitively different with respect to the case of scalar DM. Indeed the pair annihilation cross-section is not velocity suppressed and can be schematically expressed as: In addition, the t-channel mediated annihilation process ψψ → Z Z is particularly efficient, where the corresponding cross-section is given by: As evident that a strong enhancement is originated by the velocity dependent term being proportional to is well-known that this kind of behavior leads to a violation of perturbative unitarity unless new degrees of freedom, like a dark Higgs [168], are added to cure the pathological behavior of the theory. In absence of a UV completion, we have imposed, in our simplified framework, a unitarity constraint on the axial coupling A Z ψ of the form:

Vector dark matter
As stated already we will discuss separately the cases of Abelian (real vector) and non-Abelian (complex vector) DM, in order to exploit different scenarios regarding DD. Similar to the case of SM Z -portal, we will consider the following two constructions for non-Abelian and Abelian DM, respectively: where the second terms represent interactions among Z and the SM fermions. As a convention we have normalized, in both cases, the DM coupling to the new gauge coupling g . As already discussed in the case of Abelian DM we have considered a Chern-Simons type interaction [210]. The interaction term of the complex vector DM can instead arise at the renormalizable level by considering the DM as vector boson of an additional non-Abelian group, the minimal option would be SU (2), and the exact mimic of the SM EW group SU (2) L × U (1) Y (this would require the presence of an additional Z which we assume to be heavy enough to have a negligible impact in the low-energy phenomenology). The parameter η Z V contains the model specific information for The most important difference among the two scenarios relies in the DM-nucleon scattering cross-section. In the non-Abelian case, interaction with the vectorial current qγ μ q are possible, thus leading to the SI cross-section: In the Abelian case, on the contrary, the only (momentum) unsuppressed interaction, is with the axial-vector quark current qγ μ γ 5 q, so that the interaction with nucleons is SD with cross-section given by: Here Δ p i denotes spin content of the 'i'-th quark flavour inside proton and S A p , S A n represent proton and neutron contribution to the spin of nucleus, respectively. For Xenon based detectors S A p S A n so that the reference cross-section is the one of DM on neutrons.
The Abelian and non-Abelian DM have very different properties also concerning annihilations. In the first case we see that the annihilation cross-section into SM fermions is strongly suppressed, being, in fact, given by: Its first non-zero contribution, the p-wave, is suppressed by the final state fermion mass. Ad exception of the case in which the mass of the Z is not too far from the one of the top quark, the DM annihilation cross-section is actually dominated by the d-wave contribution and hence, suppressed as v 4 . On the contrary, the DM features a very efficient annihilation into Z pairs, when allowed kinematically: Its contribution to the DM relic density is nevertheless limited by the unitarity constraint [273,274].
In the case of non-Abelian vector DM the annihilation cross-section into SM fermions is only p-wave suppressed: while the annihilation cross-section into Z Z , instead, features a similar behavior with respect to the Abelian case: The impact of different theoretical and experimental limits on the two scenarios are shown, as customary, in the plane (m Z , m V ), in Fig. 15. In the case of non-Abelian DM the weaker limits from DD do not coincide with a larger viable region for thermal DM since the contemporary suppression of the DM annihilation cross-section into fermions (the annihilation into Z Z is strongly limited by unitarity) allows for the correct relic density only above the limit from LHC dilepton searches ad exception of a tiny region corresponding to the s-channel resonance. Much worse is the situation in the case of a non-Abelian vectorial DM. Indeed, already current limits from DD overpower the ones from LHC and exclude thermal DM for both m Z and m V below 5 TeV.  Table 1). Colour scheme is the same as Fig. 13

t-channel portals
We consider, in this section, the case in which the DM, instead of coupling in pairs, is coupled with one mediator state and a SM quark. Keeping the assumption that the DM is a SM singlet, the mediator field should carry at least non-trivial charges under the SU (3) C and U (1) Y of the SM.
We will then consider the following Lagrangians: 20 These Lagrangians describe the couplings of a fermionic (scalar) DM ψ (χ ) with scalar (fermion) mediator fields u , d , Q , (Ψ u , Ψ d , Ψ Q ) having quantum numbers: Hence, these mediators can interact with the left-handed quark doublet, up-type and down-type right-handed quarks, 21 respectively. To ensure the cosmological stability of the DM some suitable additional quantum numbers for both the DM and mediator should be assumed, so that couplings between any of these two states and only SM states are forbidden. 22 It is also evident that the DM can be stable only if it is lighter than the t-channel mediator. In order to avoid possible occurrence of flavour violating effects, with consequent strong constraints on the couplings λ Ψ u,d,Q (λ u,d,Q ), we assume that the mediator carries also a flavour quantum number (a "flavored DM" [278][279][280] would be equally feasible), for example u ≡ (σ u , σ c , σ t ) (Ψ u ≡ (ψ u , ψ c , ψ t )), and that the interactions of Eq. (63) are flavour conserving. This is achieved by assuming the components of the mediator field to be degenerate in masses. We call these masses m u,d,Q and m Ψ u,d,Q and the couplings λ Ψ u,d,Q , λ u,d,Q (being actually matrices) and considered them to be diagonal in the flavour space. For further simplification, we will assume all of these couplings to be equal and thus, drop the flavour indices.
Concerning DM phenomenology the relic density is determined, for all the three types of candidates, by annihilation processes into fermion pairs induced by t-channel exchange of the mediator field. For close values of the DM and mediator masses co-annihilation processes, like ψ u,d,Q (Ψ u,d,Q χ) → (u, d, Q)g and mediator pair annihilation processes like, 20 Given the non-trivial gauge charges, the t-channel mediators can also couple with the gluons, the photon and the Z -boson. We do not explicitly write such interactions for simplicity. In the case of the scalar mediator a renormalizable 4-field coupling with the SM-Higgs doublet, i.e., H † H † q q , can arise at the tree-level. We will assume that the corresponding coupling is negligible. 21 In Eq. (63) the labels u and d refer globally to up-and down-type quarks. The couplings and the masses of the mediator fields carry also a generation index which is not explicitly reported (see main text for further clarification). 22 By relaxing this hypothesis it is possible to have a viable decaying DM candidate model with characteristic phenomenology [275][276][277]. by gauge interactions, might also become important. In analogous fashion, as the other models reviewed in this work, we have focused on the assignations λ Ψ u,d,Q , λ u,d,Q = 1. In such a case the dominant contribution to the DM relic density comes from the DM pair annihilations into SM fermions. We can then achieve an analytical description through some simple approximations of the corresponding cross-sections, based on the velocity expansion (more complete expressions are presented in the Appendix): In the cases of complex scalar and Majorana fermion DM the s-wave term of the annihilation cross-section is helicity suppressed so that, ad exception of values of the DM mass close to the top-quark mass, the dominant contributions comes from the p-wave term, leading to the following estimate: On the contrary, the annihilation cross-section of a Dirac fermion DM is s-wave dominated and can be estimated as: For simplicity we have referred, for analytical expression and numerical estimates, to the case of u (Ψ u ) mediator. The extension to the other cases is straightforward; as can be easily noticed that the main difference is for the case in which the DM and the mediator are coupled with right-handed downtype quarks and the terms dependent on the top-quark mass are absent.
Concerning DD, it relies on the scattering of the DM with up-quarks through s-channel exchange of the mediator. In the cases of complex scalar and Dirac fermion DM these interactions lead to SI cross-sections (in the case of a Dirac fermion DM also SD scattering is present but its impact is negligible given the much weaker experimental limits) whose expressions are [61,69,140]: (Dirac fermion), (68) On the contrary, in the case of a Majorana fermion DM, one should consider SD interactions described by the following cross-section: here μ ψn represents reduced mass of the WIMP-neutron system and Δ n u , Δ n d , Δ n s denote spin content of the up, down and strange quark inside neutron.
The results of our analysis are presented, in the usual fashion, in Figs. 16 and 17 for scalar and fermionic (both Dirac and Majorana) DM, respectively.
For simplicity we will refer only to the case in which the DM interacts only with one kind of mediator, chosen to be u (Ψ u ). As evidenced by the expressions above, differences between the various cases mostly arise in the scattering crosssection because of the different interactions with the lightquarks, leading in turn to different interactions with proton and neutrons. This however does not lead to dramatic variations in the picture presented in Figs. 16 and 17. Similarly, there is not loss of generality in restricting to the hypothesis  16 Combined constraints for a complex scalar DM χ coupled with the right-chiral up-type quarks through a Dirac fermionic t-channel mediator Ψ u . The results are in the bi-dimensional plane (m Ψu , m χ ) and the coupling λ Ψu has been set to 1. The red curve corresponds to the contour of correct DM relic density. The blue region is excluded by the current constraints while the magenta region will be excluded in the absence of signals from XENON1T after 2 years of exposure. The grey region corresponds to m χ > m Ψu where the DM is not cosmologically stable of a single mediator field. Adding different mediators would lead to simultaneous enhancement of the DM annihilation and scattering cross-section causing a global shift of the isocontours in the bi-dimensional planes of the figures. There is a potential exception though. In the case when both Q and at least one between u and d are considered and a mass mixing is allowed between these fields (this occurs for example in Supersymmetric theories), Majorana fermion DM would acquire a sizable SI scattering cross-section [281][282][283]. We will not explicitly discuss this kind of scenario. Figures 16 and 17 clearly show that the cases of complex scalar and Dirac fermion DM, in which unsuppressed SI interactions are present, are already substantially ruled out. 23 More interesting is the case of a Majorana DM. A sizable portion of the parameter space corresponding to the correct DM relic density survives from the present constraints but is within the reach, for masses of both the DM and the mediator up to a few TeV, of the near future DD experiments. Possible complementarity with LHC searches is also worth of being explored.  [140,284]. Similar limits for the cases of Dirac fermion and complex scalar DM have been omitted since they are irrelevant with respect to DD exclusions Collider searches can probe the DM pair production through mono-X events and pair production of the scalar or fermionic mediator with subsequent decay into a DM particle and a SM quark through dijet + missing energy events. The latter signature is also considered in supersymmetry (SUSY) searches since it would originate from squark pair production, with subsequent decay into a neutralino and a quark, under the assumption of decoupled gluino. Limits from SUSY searches cannot be straightforwardly applied to the scenario under consideration since the pair production cross-section of the scalar mediator field is enhanced by the presence of diagrams with t/u channel exchange of the DM (the corresponding diagrams in SUSY are suppressed since the concerned couplings of the neutralino with quarks and squarks cannot be of O(1)). We have then applied the limits from the dedicated studies [140,284] (see also [285,286]), in which searches of both monojet and dijet + missing energy events have been considered, and reported them in the bottom panel of Fig. 17 (orange region with dashed orange border). Similar constraints have been also obtained in the case of Dirac and complex scalar DM. These are, however, overwhelmed by the constraints from DM DD and thus, omitted in the figures.
Our discussion can be extended further by considering mediators charged only under EW interactions while being colour singlet. Even if no interactions between the DM and the quarks are present at the tree-level, DD might still represent a more effective probe with respect to other DM searches.
A compelling case is represented by Dirac fermion DM where SI interaction can be induced at the loop level. The DM indeed interacts with the quarks through an effective vertex generated by a triangle loop formed by one scalar field and two SM leptons. The mediator of the interaction is mostly the photon. The corresponding scattering rate is given by [278,287,288]: where v E is the DM velocity in the Earth frame, E R is the recoil energy and F em is the electromagnetic form factor [289]. For definiteness we have considered a scalar field e with the same quantum numbers as of the right-handed charged leptons (the assumptions on the flavor structure of the couplings with the DM are the same as stated at the beginning of the section).
The impact of the corresponding current/projected constraints is shown in the usual fashion in Fig. 18. Although sensitively weaker with respect to the case of charged scalar mediator, DD constraints are still capable of ruling out the region corresponding to the correct DM relic density, for order one assignation of the coupling. As further test of this Fig. 18 The same as the top panel of Fig. 17 but for a scalar mediator with the same quantum numbers as of the right-handed charged leptons scenario one could apply the limits from searches of direct production of sleptons in SUSY models [290]. These limits, however, are weaker with respect to the ones from DD and hence, for simplicity, have not been reported in Fig. 18.

Pseudoscalar portal
We will investigate in this subsection the phenomenology of a pseudoscalar s-channel portal. Under the assumption of CP-invariant interactions, as performed throughout this paper, only a fermionic DM can be considered in this case. For simplicity we will describe only a Dirac fermion DM since the Majorana case features no substantial differences. We will then consider the following Lagrangian: where we have assumed, similar to the case of a scalar mediator, Yukawa-like couplings of the mediator with the SM fermions, ensuring a SU (2) L invariant construction (more complete realizations of this scenario have been considered e.g. in [165,[291][292][293]). The DM relic density is determined by annihilation into the SM fermions pairs and, when kinematically accessible, into aa pairs. At the leading order in the velocity expansion the corresponding cross-section can be analytically approximated as follows: This cross-section is s-wave dominated (hence capable of indirect DM signals). Given the Yukawa structure of the couplings with the SM fermions this cross-section is sizable, away from s-channel resonances, only for m ψ > m t . In such a case a numerical estimate is given by: The annihilation cross-section into aa pairs is, instead, p-wave suppressed: The most peculiar feature of the pseudoscalar portal scenario is, nevertheless, the weakening of the interactions possibly responsible of DD. 24 Indeed, tree-level interactions between the DM and the SM quarks (and gluons), mediated by a pseudoscalar, are momentum suppressed. This can be described by the following differential interaction crosssection for a target nucleus of mass m T [295]: where v E represents the DM speed in the Earth frame, Δ N q denotes quark spin content of the nucleon, E R is the nuclear recoil energy with recoil velocity q, and F N N are (squared) form factors whose (approximate) analytical expressions are given in Ref. [296]. The cross-section of Eq. (75) corresponds neither to SI nor to SD (although the latter is a good approximation) interactions; for this reason we have expressed it directly in terms of a differential cross-section . The red lines are the iso-contours of the correct DM relic density. In the yellow regions the DM annihilation crosssection into SM fermions, computed at the present time, exceeds the limit given by FERMI from searches of signals in dSphs [105]. In the orange regions the loop-induced annihilation cross-section into photon lines exceeds limits set by searches of gamma-ray lines [102]. For (λ a ψ , c a ) = (1, 1) an exclusion region (grey) from the LHC searches of monojet events [133], as well as projected exclusion regions from XENON1T (magenta) and LZ (purple) are also reported on a nucleus. As can be easily argued that the q 4 dependence of Eq. (75) implies a very suppressed scattering rate so that a potentially detectable signal is produced only for m a ∼ O(10 MeV) [295]. These low values are, however, subject to bounds from low energy observables and rare flavor processes (see e.g., Ref. [297] for an extensive analysis). For these reason we will consider here (and also in this section) values of m a above 1 GeV. 25 We also remark that Xenon based detectors, like the ones considered in our study, would be in any case not suited for probing the interaction crosssection of Eq. (75) since it originates from a coupling of the DM mostly with unpaired protons, not present in Xenon, which instead has isotopes with an odd number of neutrons.
In the considered setup the most relevant interactions originate at the loop level, from a box diagram in which a pair of pseudoscalars, one DM and one quark state are exchanged [291], and are of SI-type. The corresponding cross-section is given by: 25 Values of m a below 10 GeV are still partially in tension with low energy constraints for c a = 1. We will ignore this issue since it has marginal impact on our discussion.
where the sum in the first row involves all quarks. The coefficients f p q are the same as the ones defined in the case of real scalar portals (we remind f p c,b,t = 2 27 f p T G , as mentioned in the context of Eq. (35)).
As evidenced in Fig. 19, this cross-section is very suppressed so that no constraints come from present experiments. On the contrary, for O(1) values of both the λ a ψ and c a couplings, next generation detectors can partially probe the parameter space corresponding to the thermal DM.
More stringent constraints come from ID of the DM. Contrary to the scalar mediator case the DM annihilation crosssection into SM fermions is now s-wave dominated. These processes lead to potential signals in the gamma-ray continuum which can be probed by the FERMI satellite [105]. In addition, the DM can also pair annihilate into two photons, through an effective coupling between the pseudoscalar mediator and photons, generated by triangle loops of the SM fermions [105]. This process, responsible for the generation of gamma-ray lines, is strongly constrained by the negative results in present searches [102].
As already anticipated that the most stringent constraints comes from indirect DM searches in dSphs. 26 Indeed, the absence of signals excludes thermal DM for mass below approximately 50 GeV. The most disfavored scenario turns to be the one corresponding to the assignation λ a ψ = c a = 1 (middle plot of Fig. 19). Indeed a complementary constraint comes from searches of gamma-ray lines so that the excluded ranges of the DM masses reaches the order of 100 GeV. This specific assignation of the couplings has also been investigated at the LHC through searches of events with monojet and missing transverse momentum. Due to the absence of any signal, the regions of the parameter space corresponding to 100 m a 500 GeV as well as m a > 2m ψ are currently excluded.

Scalar + light pseudoscalar portal
In this subsection we will consider the case in which the pseudoscalar mediator is actually a component of a complex scalar field Φ → (S + ia)/ √ 2, described by the following Lagrangian: After the EWSB, the scalar component of Φ gets non-zero VEV (v Φ ), generating a mass term of its scalar component, m S = √ 2λv Φ while leaving the pseudoscalar component massless. A mass, m a = √ 2 Φ , for this second field is originated by an explicit mass term 2 Φ Φ 2 /2 in Eq. (77). In this kind of setup it is rather natural to identify the pseudoscalar component with a pseudo-Goldstone boson associated to a U (1) global symmetry carried by the complex scalar Φ and then assumes that m a m S [298][299][300]. We further assume that after the EWSB it is possible to write interaction terms of the fields S and a both with the SM fermions and a fermionic DM candidate as follows: 26 For the mass range assumed for the DM, the ID signal probed by dSph relies on the bb final state for m ψ < m t and on the tt final state for m ψ > m t . For simplicity we have adopted everywhere the limit from the bb final state. This is a good enough approximation for the purposes of our study.
Here we have again assumed Yukawa-like interactions among S, a and the SM fermions where the concerned couplings, including a normalization factor of 1/ √ 2, are parametrized as c S . The other free parameters are the DM mass and coupling, m ψ and g ψ . 27 The main feature of this double s-channel portal scenario is the lower amount of correlation between the direct DM detection and the relic density, due to the presence of a light mediator state (see also Ref. [301] for a similar idea).
Indeed, DD of the DM relies only on the coupling of the DM with the scalar component of the Φ field. The DM relic density, instead, is determined by the different annihilation channels, including f f , aa, Sa and SS involving interactions of both S and a fields. The expressions for the annihilation cross-sections into fermion and S pairs can be derived straightforwardly from the cases of the scalar and pseudoscalar portals and thus, won't be re-discussed in detail. The annihilation cross-section into a pairs, despite the presence of an additional contribution from s-channel exchange of S is only moderately altered with respect to the case of the pseudoscalar portal (see Sect. 10.1) and we just then refer to its detailed expression presented in the Appendix. The most prominent feature, relevant for the DM relic density, is the presence of Sa as final state for annihilation processes, when allowed kinematically. To this corresponds, in fact, an efficient s-wave annihilation cross-section which can be analytically approximated as: The comparison between the DM relic density and the experimental (present and future) constraints, is performed, in the usual fashion, in Fig. 20. Here we have chosen the DM mass m ψ and the scalar mass m S as the two free parameters while m a has been set to 5 GeV, in order to avoid dangerous constraints from low energy physics. We have then considered three assignations for (g ψ , c S ) couplings, i.e., (1, 0.25), (1, 1) and (0.25, 1) while the value of coupling λ of the scalar potential has been set to 1.
As evident that the presence of annihilation channels, like aa and Sa, involving non-SM light states allows to achieve the correct relic density, compatibly with constraints from The presence of s-wave unsuppressed annihilation channels like f f (contribution from s-channel exchange of the pseudoscalar) and Sa 28 makes mandatory to consider, besides DD, also limits from ID. However, for the chosen parameter assignation, these last constraints have no impact in the region corresponding to the correct DM relic density.

Portals to secluded sectors
In this section we will consider cases in which the mediator of the DM interactions has no direct coupling with the SM fermions. Dark portals can nevertheless be realized, at the renormalizable level. Indeed the SM features two Lorentz and gauge invariant bilinears, i.e., H † H and B μν , which together with the similar appropriate structures can connect the SM and BSM sectors. The first one, H † H , can be coupled to another scalar bilinear. In case this second scalar field also has non-zero VEV, a mass mixing with the SM-Higgs is generated so that the portal interactions with the SM fermions, as well as the gauge and the SM-Higgs bosons itself are produced. The field strength B μν can be coupled with the field strength of another U (1) gauge bosons. Such kinetic mixing 28 For simplicity, we have considered in our analysis only limits from searches in dSphs. As pointed out in Refs. [299,302,303] the annihilation into Sa can lead to box shaped gamma-ray signals which could be probed in the near future by CTA. However, our choice for m a , implies a too suppressed decay branching ratio for a → γ γ process. For this reason we have not considered explicitly this possible signal in our analysis.
term after the EWSB, is the origin of a mixing between the Z and the new Z boson.
In both these two scenarios, the DM interacts with the SM fields (both fermions and bosons) through a double schannel mediator. The relevant phenomenological processes are substantially the same as already investigated in the cases of the SM and BSM s-channel dark portals. Contrary to these scenarios, we cannot consider O(1) couplings between the SM states and the BSM mediators, as well as between the DM and the SM ones, since the mixing between the SM-Higgs and an additional scalar or the Z and the Z' are required to be small by several experimental and theoretical constraints.

SM-Higgs + spin-0 portal
In this subsection we will revisit the phenomenology of a real spin-0 mediator Φ considering the more realistic case in which it also features interaction with the SM-Higgs doublet H . We will thus assume the following scalar potential: where V H,SM is the SM scalar potential ( here μ h , μ Φ are parameters with the dimension of mass and μ 2 h , μ 2 Φ < 0 for spontaneous symmetry breaking. One should note that in V H,SM + V (H, Φ), the condition for getting a positive definite mass spectrum requires 4λ h λ Φ > λ 2 hS while λ h , λ Φ > 0 are necessary to get V H,SM + V (H, Φ) bounded from below. Combining these two conditions we see that λ hS can take both the positive and negative values. We denote non-zero VEV of the scalar field Φ as v Φ , so that it can be expanded as The coupling of Φ with the SM-Higgs doublet H induces mass mixing such that, after the EWSB (assuming unitary gauge), one can define the following two mass eigenstates [304]: where H (0) represents the electrically neutral scalar part of the SM Higgs doublet H and the mixing angle θ is defined by: The phenomenology of the mediator sector, thus, can be expressed as functions of the five free parameters λ h , λ Φ , v h , v Φ , λ hS or equivalently in terms of m 2 h , m 2 S , v h , sin θ, λ hS using the following relations [304][305][306]: − sin 2 θ , Note that in reality we have only three free parameters, i.e., λ hS , sin θ and m S as m h , v h have known measured values [137].
The mixing between H (0) and φ (see Eq. (82)) indicates that both the mass eigenstates (h, S) will couple to the SM states as well as to the DM and thus, it represents a two-portal scenario.
The couplings of the two s-channel mediators with the W ± , Z -bosons and the SM fermions are described by: while their cubic self-couplings are given by: with: The parameters λ hS and sin θ are subject to several experimental and theoretical constraints (see e.g., Refs. [304,[307][308][309] for an extensive discussion). For example, a non-zero value of θ modifies the couplings of the SM-Higgs with SM particles and thus, is constrained by the measurement of the Higgs signal strengths. The coupling λ hS is instead constrained by the stability of the scalar potential. Most of these constraints become increasingly stringent as m S decreases; for this reason we will focus our analysis on the case m S > m h .
Similar to the other spin-0 mediator scenarios, this extended Higgs sector will be coupled to a scalar DM χ , a fermionic (we will restrict to the Dirac case) DM ψ and a spin-1 DM V μ . We will consider the following Lagrangians for the corresponding interactions.
In the case of a scalar DM 29 we consider the following simplified (i.e., assuming vanishing coupling for |χ | 4 term) interaction: which, after the EWSB, leads to the following effective Lagrangian: where: 29 Stability of a scalar DM is protected either by a U (1) (complex scalar) or by a Z 2 symmetry (real scalar) and thus, the associated Lagrangian contains |χ| 2 or χ 2 term, respectively. In this analysis we have used the notation |χ| 2 for generality although the DM is considered to be a real scalar.
using Eq. (84). From Eq. (90) it can be seen that without the loss of generality one of the couplings λ χ H and λ χ Φ can be set to zero. We will thus, pursue this minimal choice setting λ χ H = 0. With this choice one can consider the following five parameters λ χ Φ , λ hS , sin θ, m S and m χ as the free inputs. In the case of a fermionic DM it is natural, in this setup, to assume a Yukawa-type coupling with the field Φ of the form: L ψ = −y ψ ψψΦ; (91) here the DM mass is dynamically generated by the VEV of Φ so that the DM mass and its coupling with Φ are not independent but are related as y ψ ∝ m ψ /v Φ . This choice allows to reduce the number of free parameters compared to the scalar DM scenario and hence, λ hS , m S , sin θ and m ψ remain the four free input parameters. The effective couplings of the DM with the S and h fields can be straightforwardly derived by using Eqs. (82) and (84). A dynamical generation of the DM mass has also been considered in the case of a vectorial DM. We indeed identify the DM as the stable gauge boson of a U (1) dark gauge group, 30 spontaneously broken by the VEV of a complex scalar field Φ. The interactions between the latter and the DM are then embedded in the covariant derivative 31 After spontaneous symmetry 30 The DM as a gauge boson of the dark sector connected to the SM through the SM-Higgs sector has also been proposed, for different DM production mechanisms, in Refs. [310,311]. 31 We have adopted the same definition of the covariant derivative as of Ref. [312].
breaking the DM Lagrangian reads: where we have used the expression of Φ in the unitary gauge It is now apparent that one can consider the following four parameters, i.e., m S , sin θ , η S V and m V as the free inputs for a vectorial DM in the SM-Higgs + spin-0 portal scenario. Hence, one should replace the new derived parameter λ hS judiciously in Eqs. (83), (84), (87) and (90) for numerical analysis.
The processes responsible for the DM relic density and DD have been already discussed in detail for the cases of the SM-Higgs and scalar portal individually; we thus just illustrate the results of our analysis, as reported in Fig. 21. We just remind here that all the considered types of DM candidates feature SI interactions with nuclei.
The results are reported in the usual bi-dimensional planes, i.e., (m S , m χ ), (m S , m ψ ) and (m S , m V ), respectively. The angle θ has been conservatively set to sin θ = 0.1 in order to comply with the constraints from the SM-Higgs signal strengths. Similarly, the coupling λ hS has been set to a value of − 0.1 in order to satisfy various constraints on the SM-Higgs sector (see e.g., Ref. [304] for a detailed discussion). The couplings λ χ Φ , η S V have been set to 0.1 and 1, respectively while for the fermionic DM the associated Yukawa coupling y ψ is not a free parameter in our construction.
It is evident that the outcome of our analysis presents some sensitive differences with respect to the case of the spin-0 mediator discussed in the previous sections. In the case of a fermionic and a vectorial DM the limits from DD are rather effective, due to the presence of an additional light mediator. In the case of a vectorial DM the only viable region, for masses of the DM and the mediator below the TeV scale, corresponds to the case m V > m S , thanks to the enhancement of the DM annihilation cross-section due to V V → SS process. In the case of a fermionic DM the region surviving the current constraints corresponds to the "pole" m ψ ∼ m S /2. In both cases, the next generation of DD detectors will probe the WIMP paradigm for masses of the DM and the mediator up to a few TeV. More particular is the case of a scalar DM where the shape of the DD contours is rather different compared to the other spin-0 mediators. This is due to the different choice of the energy scale, the VEV v Φ rather than the mass of the mediator, which implies a larger cross-section at higher values of m S . On the contrary, as can be noticed from Eq. (90), for λ χ H = 0 the DM couplings become smaller as m h and m S get close in their values. 32 For the chosen parameter assignations, a large region of the viable thermal DM is present in the regime m χ > m S for a rather light m S . A sizable part of this region will be excluded in the absence of signals at XENON1T and LZ.

Kinetic mixing
We will reconsider in this subsection the scenario in which the DM is coupled to the gauge boson of a new U (1) group. In this case, however, SM fermions won't be charged under the new gauge group so that no direct couplings with the new gauge boson are induced. The "dark" and visible sectors can nevertheless be connected by a kinetic mixing operator B μν X μν which can already exist at the tree level, being both Lorentz and gauge invariant, or generated radiatively (for example, if the new gauge sector features new fermions having non-trivial quantum numbers also under the SM gauge group [315]). 32 In the case of a scalar DM it would even be possible to generate a "blind spot" in the scattering cross-section for a rather specific choice of the couplings λ χ H,Φ [221,313] which would induce a destructive interference between the amplitudes associated to the h and S exchange. A different solution for relaxing DD constraints in the presence of multiple scalar singlets coupled to the SM-Higgs has recently been proposed in Ref. [314].
We will then consider the following Lagrangian with δ as the kinetic mixing parameter: where B μν and X μν are the fields strength of the U (1) Y and the new U (1) group (associated with X μ ), respectively. L SM is the SM Lagrangian besides the already written kinetic term of B μ , while L DM is the DM Lagrangian including the kinetic and mass terms as well as a coupling with the X μ boson. These terms will depend on the spin assignation of the DM. We will consider the following cases: As already discussed, a natural option to couple a scalar or a fermionic DM to a new gauge boson is to assume it charged under the new symmetry group so that its interactions originate from the covariant derivative D μ = ∂ μ − ig X X μ . 33 Slightly more complicated is the case of a vectorial DM. Here we have assumed an analogous coupling, as the one of the SM, between the Z and two W bosons. 34 The coupling η X V encodes the gauge coupling g X and other eventual extra factors arising from a specific model construction. m χ, ψ, V denote the respective DM masses of different spin assignments.
The kinetic term of Eq. (94) should be diagonalized and canonically normalized. After the EWSB, it is possible to define three mass eigenstates, i.e., A μ , Z μ and Z μ , for the electrically neutral gauge bosons through the following two transformations [246][247][248]316]: 33 Note that in the definition of g X we have also encoded the value of U (1) charge of the DM. 34 The case of an Abelian vectorial DM in the kinetic mixing scenario will be addressed in a dedicated future publication.
where t δ , c δ = tan δ, cos δ, cŴ , sŴ = cos θŴ , sin θŴ , c ξ , s ξ = cos ξ, sin ξ and the angle ξ is defined by: Note that the transformation given in Eq. (96) leads to physical solutions only if one of these two conditions is met [247]: In the expressions above sŴ , mẐ do not represent the experimental measures of Weinberg angle and the Z -boson mass but are related to the latter, i.e., s W ≡ sin θ W , m Z , as function of δ. Indeed, since the photon coupling does not change once passing to the "physical/mass" basis, one can write: In an analogous fashion, the invariance of the W -boson mass under the transformations of Eq. (96) allows to relate the kinetic mixing parameter to the ρ parameter as: which can be reformulated as: where Δ = ρ − 1 and t 2 W = tan 2 θ W . Hence, ρ − 1 = 4 +8 −4 × 10 −4 measurement [209] can be used to constrain the parameter δ.
The kinetic mixing parameter is further constrained by EW Precision Tests (EWPTs) [317,318] which reads as: The spectrum of the neutral gauge bosons features one massless eigenstate, coinciding with the SM photon, and two massive states. Their masses, in the experimentally favoured limit, i.e., sin δ, sin ξ ∼ δ, ξ 1, are given as: where m Z must coincide with the experimentally measured value of the Z boson mass. The interactions (relevant for the DM phenomenology) of the Z , Z with the SM states are described by [247]: where: , , Here T 3 , Q are the isospin quantum number and electric charge of the associated SM fermions. The coupling of the DM of various spin allocations, i.e., scalar, fermion and vector, with the two mass eigenstates are respectively given by: with g Z DM = cos ξ cos δ and g Z DM = − sin ξ cos δ . As evident, in the physical basis, the DM is connected to the SM sector by two s-channel mediators, the Z and the Z .
The DM relic density is determined by annihilation processes into the SM fermion pair final states, W W and Z (Z )h, induced by s-channel exchange of the mediators, and Z Z, Z Z and Z Z , induced by t-channel exchange of a DM state. The corresponding rates can be straightforwardly derived from the cases of Z /Z portals so won't be re-discussed in detail here. Similar to the scenarios already described, the DM DD relies on the SI interactions, which induces a scattering cross-section written, for the case of proton, as: Here μ χ p/ψ p/V p , as already defined, denotes reduced mass of the concerned WIMP-proton system. c X encodes the DM couplings g 2 X (for complex scalar and fermionic DM) or η X V (for vectorial DM) and other eventual overall factors, depending on the kind of DM candidate.
The interplay between the DM relic density and DD is shown, for the various spin assignations of the DM, in Fig. 22. We have considered two assignations of the kinetic mixing parameter, the first one is derived from the present limit of the EWPT (i.e., taking the equal sign in Eq. (102)) while the second one corresponds to a constant value δ = 0.01. This last choice is inspired by models in which the kinetic mixing parameter is radiatively generated upon integrating out heavy degrees of freedom charged under both the new U (1) gauge group and the U (1) Y of the SM [315,319].
By comparing the outcome of Fig. 22 with the scenarios where the Z is directly coupled to the SM fermions, we notice that DD probes a more limited region of the parameter space. This is because the scattering cross-section depends on coupling suppressed by the smallness of the parameter δ. However, this kind of suppression affects also the DM annihilation processes, ad exception of the Z Z final state so that a strong tension with the experimental constraints, for example similar to the SSM (see Sect. 8.2), still persists. We, indeed notice that the thermal DM is excluded for masses below the TeV scale unless small values, O ∼ (0.01), of the kinetic mixing parameter are taken. Even in such a case, the correct relic density is achieved only at the pole m Z /2 or for m χ,ψ,V > m Z , where annihilation into Z Z is accounted without relying on the kinetic mixing parameter. These setups will nevertheless be excluded in the absence of signals in the next generation multi-TON experiments.
As explicitly indicated in the plots of Fig. 22 that we have considered only the case m Z > m Z . The opposite regime (often dubbed as dark photon) would also be feasible, although the constraints on δ would be even stronger because of Eq. (102) and additional constraints from muon (g − 2) and parity violating effects in the atomic physics [247] as well as from low-energy colliders [320]. We have checked that the low mass Z regime is substantially excluded by the current DD limits, unless DM masses below the sensitivity reach of DD experiments are considered. Thus, we have not reported it explicitly.

DM (partially) charged under SU(2) × U(1)
All the models considered until now have assumed the DM to be a SM singlet. In this section we will, instead, consider the case in which the DM particle has non-trivial quantum numbers under the EW component of the SM gauge group. A rather straightforward realization consists into considering the DM as the lightest neutral component of a SU (2) multiplet. Extensive discussions can be found in Ref. [321] and Refs. [169,170] for scalar and fermionic DM, respectively. These models cannot strictly be classified as Dark portals since the DM relic density is mostly set by annihilations processes into Z and W boson which are determined by the gauge couplings between the latter and the components of the DM multiplet. These are, nevertheless, very interesting and predictive models, having the DM mass (for this reason they are dubbed as "Minimal DM") as the only free parameter. For these models the DM relic density is typically strongly suppressed by efficient annihilation processes into W + W − and Z Z final states unless the DM is sufficiently heavy, possibly above the TeV scale. For these high values of the DM masses the treatment presented in Sect. 2 should be refined since, besides co-annihilation processes with the other states belonging to the DM multiplet, one should also account for Sommerfeld enhancement effects [170,171,321] and even effects from bound states formation [322]. These effects are typically associated with an enhancement, with respect to the conventional treatment, of the DM annihilation cross-section so that the correct DM relic density is reached, for some realizations of the Minimal DM, for masses even of the order of 10 TeV. Concerning detection strategies the best probe, at the moment, is represented by ID [323]. On  (right column) DM, respectively, interacting with a Z , kinetically coupled with the SM Z boson. In the top-row plots the kinetic mixing parameter δ has been set to the maximal value, as a function of m Z , consistent with the EWPT constraints while for the bottom-row plots δ has been set to a constant value of 0.01. We set g X = 1 for all these plots. In this figure the red curve represents the contour of correct DM relic density. The blue region is excluded by the current constraints from XENON1T while the magenta and purple regions would appear excluded in the absence of signals from XENON1T (after 2 years of exposure) and LZ, respectively the contrary, scattering cross-section on nucleons are typically rather suppressed [324][325][326] and hence, DD probes are not very efficient. Further, given the high values of the DM masses corresponding to the correct DM relic density, the collider searches could, similarly, hardly probe these scenarios in the near future [326].
As already mentioned that the minimal DM scenarios, despite being phenomenologically very interesting, do not strictly belong the category of Dark Portals. Moreover, they are not the best target for direct DM searches. For these reasons they will not be discussed here in further detail.
Models more similar to Dark portals, achieving the correct DM relic density at lower values of the DM mass, consist of the case in which the DM is the mixture of the neutral component of a SU (2) multiplet and a SM SU (2) singlet. In such a case the coupling of the DM with the gauge bosons is suppressed through the mixing between the SU (2) singlet and the multiplet(s). Similar to Dark portals, annihilation into the SM fermions, through s-channel mediation of the SM Higgs (or the Z boson), can play a relevant role for the DM relic density. Viable relic density is potentially achieved also for the DM masses of the order of a few hundreds GeV. The most popular example of these kind of scenarios is represented by Supersymmetric theories where the DM candidate, the lightest neutralino, is indeed a mixed state comprised of the bino (SM singlet), the Wino (electrically neutral component of a SU (2) triplet) and two higgsinos (electrically neutral components of SU (2) doublets).
We also remark that these kind of constructions represent renormalizable completions of the SM Higgs and Z portal models discussed in the Sect. 7. As will be shown subsequently, the phenomenology of these more realistic constructions differs from the simple SM Higgs and Z-portals in particular because of the presence of additional states belong-ing to the DM sector which can affect the DM relic density through co-annihilations.

Inert doublet model
The IDM enlarges the SM Higgs sector with an additional doublet, giving rise to the following scalar potential: with H 1 and H 2 being, respectively, the SM and the new scalar doublet. μ i and λ i are real parameters. The additional Higgs doublet is odd under a discrete Z 2 symmetry which forbids its direct couplings with the SM fermions. In addition, the new doublet does not acquire VEV when the SU (2) L × U (1) Y symmetry is spontaneously broken. After the EWSB, four stable mass eigenstates are obtained: two electrically charged, labeled as H ± , and two electrically neutral, labeled as H 0 and A 0 . Their masses are given by: A viable DM candidate, thus, is achieved as one of the states between H 0 and A 0 by setting the quartic couplings λ 3 , λ 4 and λ 5 appropriately. For simplicity, we will restrict to the case when H 0 is the DM candidate. The quartic couplings are also constrained by the vacuum stability, which imposes (at tree level): and, λ i < 4π from perturbativity. As will be shown subsequently that the most relevant DM observables depend on the combination of couplings λ L = 1 2 (λ 3 +λ 4 +λ 5 ). We will thus adopt, as free parameters for our analysis, λ L as well as the masses of the new Higgs states. The constraints from vacuum stability and perturbativity are then translated into bounds on m H 0 , m A 0 , m H ± through Eq. (112). The latter are also constrained by the EWPTs. Finally, a lower bound of 79.3 GeV on the mass of the charged Higgs is applied from LEP [137].
The direct DM detection relies on a SI interaction whose cross-section is substantially analogous to the one used for the Higgs/spin-0 portal: 35 (according to the usual convention we quote the DM scattering cross-section on protons): where m p is the mass of the proton, μ H 0 is the reduced mass of the DM-proton system and the coefficients f p and f n are the ones defined in Eq. (35), and A (Z ), is the atomic mass (number) of the target nucleus.
Concerning the DM relic density, there are similarities with the SM Higgs portal only at the low DM masses, namely below the mass m W of the SM W boson. Here the DM relic density is mostly determined by annihilation into SM fermions mediated by the h boson. As the mass of the DM approaches towards m W , the annihilation into W W finals states enhances, with respect to the SM Higgs portal model discussed in Sect. 7, by t-channel diagrams involving the exchange of H ± and contact four field interactions between a H 0 and a W ± pairs. The contributions of these diagrams to the cross-section depend on the SU (2) L gauge coupling g rather than the quartic couplings of the scalar potential. A strong enhancement of the cross-section is already produced for m H 0 m W by the process H 0 H 0 → W W * → W f f (this process is illustrated in detail, including an analytical expression of its amplitude, in Ref. [328]).
For positive values of λ L the correct DM relic density is achieved for DM masses around 70 GeV while for higher values of the DM masses, the annihilation rate into W W is so strong that the DM turns out to be underabundant. 36 For λ L < 0, destructive interference among the different diagrams contributing to the annihilation rate into W W final state might occur so that the region of viable DM can be extended, according to the values of λ L , m A 0 and m H ± , towards the DM masses slightly above 100 GeV [329].
An overview of the DM phenomenology of the IDM in the low DM mass regime is presented in Fig. 23. The two 35 In the case of extreme mass degeneracy between the H 0 and A 0 states, one should consider an additional SI interaction, namely H 0 p → A 0 p, mediated by the Z boson [327]. We will neglect this possibility implicitly by assuming a mass splitting of at least 1 GeV between these two states. This, in turn, implies that the coupling λ 5 is always not negligible. 36 There are elegant ways to overcome this underabundant DM regime in the two Higgs Doublet Model using axions [343].

Fig. 23
Summary of constraints for the IDM model in the light DM mass regime, in the bi-dimensional plane (m H 0 , |λ L |) for λ L >0 (upper panel) and λ L < 0 (lower panel). The color coding is the same as used for the SM Higgs portal model in Sect. 7: the red curve represents the iso-contour of the correct DM relic density. The blue, magenta and purple regions are, respectively, the XENON1T excluded region and the projected exclusions by XENON1T (assuming 2 years of exposure time) and LZ. The brown region is excluded by the current experimental determination of the SM Higgs invisible branching fraction panels of the figure, corresponding to the cases of λ L > 0 and λ L < 0 report, similar to the case of SM Higgs portal model, the combined constraints from the DM relic density, DD including the projected sensitivities of XENON1T and LZ and, invisible branching fraction of Higgs boson, in the bi-dimensional plane m H 0 , |λ L |. Moreover, we assumed a sizable mass splitting between the DM and the extra scalar states in order to forbid co-annihilation effects. As already anticipated, for m H 0 m h /2 the outcome is substantially analogous to the SM Higgs portal model; the viable DM region is excluded by the combination of DD and Higgs invisible branching fraction constraints, except for the pole m H 0 ∼ m h /2 which will be fully probed by the next generation DD experiments. For m H 0 > m h /2, an additional viable region is present corresponding to a scenario where the DM relic density is mostly determined by the annihilation process H 0 H 0 → W f f . Given the dependence of the corresponding rate on the gauge coupling, the correct DM relic density is achieved for very low values of λ L , passing the current DD constraints.
As the DM mass increases further the annihilation into two on-shell W (as well as into Z pair) bosons becomes accessible. An efficient annihilation channel is represented by H 0 H 0 → hh process as well. These rates are in general very efficient to account for the correct DM relic density. It can nevertheless be achieved for relatively high DM masses provided that some specific requirements on the model parameters are met. The annihilation cross-section into hh is proportional to the coupling λ L , a sufficient suppression is then achieved by going to the limit λ L → 0. Interestingly, this would also imply a suppression of the DM scattering rate on nucleons which is also proportional to the same coupling. More contrived is the case of the annihilation into W + W − and Z Z final states. The corresponding cross-sections normally receive contributions depending on the coupling λ 4 +λ 5 (in the case of W + W − ) and λ 5 (for Z Z) which are responsible of enhancement factors of respectively. It is thus possible to suppress the annihilation crosssections into gauge bosons by taking almost degenerate m H 0 , m A 0 and m H ± states. In the limit of m H 0 ∼ m A 0 ∼ m H ± , the DM annihilation cross-section into gauge bosons scales as α 2 W /m 2 H 0 and the thermally favored values are matched for masses of the order of 500 GeV [327]. In this nearly degenerate mass scenario an interesting additional phenomenon occurs. If only self annihilations of the DM are taken into account then the expected value of its abundance is typically small (i.e., much below the thermally favored value). However, since the pseudoscalar and the charged Higgs belong to the same doublet and share similar cross-sections, they get decoupled after the DM freeze-out. Therefore, despite of a small DM abundance shortly after the freeze-out, an eventual enhancement of the DM abundance appears when both the pseudoscalar and charged Higgs fields decay back to the DM pair. This is the mechanism for getting the correct DM relic density as explained in detail in Ref. [344].
In order to properly investigate the effects described above, we have conducted a scan on the parameters m H 0 , m A 0 , m H ± and λ L over the following ranges: Notice that now we have extended the scan also for the low DM mass regime including the possibility, neglected in the case represented by Fig. 24, of mass degeneracy between m H 0 , m A 0 , m H ± , so that co-annihilation processes can take place.
The model points featuring the correct DM relic density have been reported in Fig. 24 in the bi-dimensional plane (m H 0 , |λ L |) together with the existing limits and projected sensitivities from direct DM detection and Higgs invisible branching fraction. As evident that the presence of coannihilations does not allow to evade the strong constraints in the low DM mass regime. On the contrary, a new viable DM region is present for m H 0 500 GeV where although the present and expected future limits are effective but still there will be a sizable amount of model points that would survive even in case of absence of signals at LZ.
We emphasize that in this high DM mass regime, coannihilations are responsible for producing the correct DM relic density. The DM annihilations into gauge bosons are rather strong to get the correct DM relic density. That said, today these co-annihilations are absent since the coannihilating particles first got decoupled from the thermal bath and then decayed back into DM. Therefore, we are left with the DM self annihilations into gauge bosons which are of the order of 10 −25 − 5 × 10 −26 cm 3 s −1 . This falls perfectly within the sensitivity reach of Cherenkov Telescope Array (CTA). Indeed, it has been shown that such telescope array provides the most effective probe for this model, being possibly able to rule out the model till DM masses of 2.8 TeV [344,345] (see also [346][347][348] for possible indirect signals in the Inert Doublet Model).

Singlet-doublet DM
We will discuss in the following texts two scenarios of fermionic DM originating from the mixing between a SM singlet and the electrically neutral component of a SU (2) L doublet. As already pointed out that this kind of setup allows for renormalizable couplings between the DM and Higgs and gauge bosons, in particular the Z , and thus, represents a theoretical framework where the SM portal models discussed in Sect. 7 can be embedded. We also remark that the new fermions which couple with the Higgs should belong to a real representation (i.e., Majorana fermions) or form vector-like pairs, since chiral fermions, with masses originating from the EWSB, are experimentally strongly disfavored [349].

Dirac DM
The Dirac fermionic Singlet-Doublet model is realized by adding to the SM particle spectrum vector like fermions L L ,R and N L ,R , transforming as SU (2) L doublets (with hypercharge 1/2) and singlets (with hypercharge 0), respectively. The relevant phenomenology is described by the following Lagrangian [177]: where, as usual for simplicity we have omitted the kinetic terms. H is the conjugated Higgs field. A discrete Z 2 symmetry, under which the new fermions are odd, is assumed to avoid couplings of the latter (through the Higgs) with the SM fermions which in turn guarantees the cosmological stability of the DM. The DM candidate is the lightest of the two neutral Dirac fermionic states arising from the bi-diagonalization procedure, of the mass matrix: where: with: so that the two physical masses, with cθ L , sθ L ≡ cos θ L , sin θ L are given by: The physical spectrum of the theory features also a charged state ψ ± with mass m ψ ± M L . The model features four free parameters, the two masses M N , M L and the couplings y 1 , y 2 . As can be easily realized that the couplings of the DM with the SM fields essentially depends on the angles θ L and θ R which determine the amount of its component charged under SU (2) The DM relic density is determined by a large variety of annihilation processes, including fermion pair final states (from s-channel exchange of the Z and Higgs boson), W + W − final state (from s-channel exchange of the Z / h boson and t-channel exchange of ψ ± ), Z Z final state (from s-channel exchange of Higgs boson and t-channel exchange of ψ 1,2 ) and finally Zh final state (from s-channel exchange of the Z boson and t-channel exchange of ψ 1,2 ). To these, one should possibly add co-annihilation processes in analogous combination of the final states in case the additional neutral and charged fermions have masses close enough to the same of the DM.
The most relevant contribution to the DM scattering processes on nucleons comes from SI interaction mediated by the Z boson. The corresponding cross-section is very similar to the one written for the SM Z -portal model 37 where V Z u and V Z d are the vectorial couplings of the SM Z -boson with up and down-type quarks.
Similar to all the models discussed in this work, we will provide an overview of the interplay between the requirement of the correct DM relic density and the experimental constraints, mostly from DD, by considering a bi-dimensional plane in the mass parameters (M L , M N ) in this case and, a 37 SI interactions are also induced by interactions of the DM with the SM Higgs boson. Their contributions are sub-dominant and for simplicity we will omit them in our discussion. It has nevertheless been included in our numerical analysis. The correct DM relic density is achieved, through thermal freeze-out, along the red contour. The blue region is excluded by limits on SI interactions from the XENON1T experiment while the magenta and purple regions will be excluded in the case of absence of signals at XENON1T, assuming 2 years of exposure time, and LZ, respectively fixed assignation of the couplings. This kind of study is presented in Fig. 25. Contrary to most of the previously studied cases we have not chosen O(1) values of the couplings here but a sensitively lower value, i.e., y 1 = y 2 = 0.1. This is because the possibility of vector-like fermions having O(1) couplings with Higgs is theoretically rather contrived since it would result potentially dangerous effects on the stability of the scalar potential [179,350,351].
The behavior of the contour of the correct DM relic density in the bi-dimensional plane (M L , M N ) is explained as follows. In order to match the DM annihilation cross-section with the thermally favored value, a sizable mixing between the singlet and the doublet components is needed, implying M N ∼ M L (in this regime co-annihilations are also important). The correct DM relic density can be obtained for a maximal value of the DM mass of approximately 733 GeV [177]. At this value the correct relic density can be achieved when the DM mostly coincides with the neutral component of a SU (2) doublet (hence co-annihilating with its almost degenerate charged partner) and, is guided by gauge interactions of the DM with the Z and the W bosons. Note that a similar effect concerning the issue of mass degeneracy of the DM multiplet, as the one discussed in the previous subsection, occurs also in this case. Here, the presence of a state nearly mass degenerate with the DM is associated with an increase of the DM relic density, rather than a decrease as happens in the conventional co-annihilation paradigm.
Because of the Dirac nature of the DM, the associated sizable interactions with the Z boson lead to an extremely high SI cross-section, already excluded by far by the present experiments. This is confirmed by the outcome of  which shows that the region of viable DM relic density is completely excluded by the XENON1T constraint alone.
In order to evade DD constraints one should instead approach the pure singlet limit, i.e., sin 2 θ L + sin 2 θ R 1. This would, however, imply a very suppressed pair annihilation cross-section for the DM. The correct DM relic density can nevertheless be achieved through the pair annihilations of the ψ 2 and ψ ± states, provided that they are almost degenerate in masses with the DM. This would require M N ∼ M L and to avoid an enhancement of sin 2 θ L + sin 2 θ R , y 1,2 1.
The aforesaid features cannot be represented through a simple bi-dimensional plot. Hence, to investigate those characteristics, similar to the previous subsection, we have conducted a parameter scan of the four free parameters M N , M L , y 1 , y 2 over the following ranges: The results of our analysis are reported in Fig. 26 in the bi-dimensional plane (m ψ 1 , σ SI ψ 1 ). We have reported there the model points reproducing the correct DM relic density and compared them with the current and anticipated exclusions by XENON1T and LZ, respectively.
As can be seen from the figure that values of the DM scattering cross-section span over a wide range of values. The highest ones, which can be even above 10 −42 GeV 2 , correspond to the case of a highly mixed DM. In such cases the interactions of the DM with the Z boson are of similar size as expected when they are coupled via ordinary gauge couplings and are, by far, excluded by DD. At the same time, model points with the correct DM relic density are obtained also for the case of a almost pure singlet-like DM with y 1 , y 2 ∼ 10 −6 , implying negligible Z -mediated scattering rate, even surviving upon the absence of signals at LZ. For these points the DM relic density is entirely determined by the (co)annihilation processes of the quasi mass degenerate doublet partner. Notice, that the viable points are present only for m ψ 1 500 GeV. Above this value, all the points with the correct DM relic density are characterized by M L < M N , corresponding to a doublet-like DM. These points cluster at very high values of the DM scattering crosssection, beyond the range reported in the y-axis, and thus, are largely ruled-out by DD.
The presence of extra fermions besides the DM would in principle offer collider tests for the scenario under consideration. The heavier neutral fermion ψ 2 and ψ ± could indeed be produced through EW processes at colliders, followed by subsequent decay into the DM and an on-or off-shell gauge boson. However, the strong bounds form direct DM detection can be evaded only through a very compressed spectrum for the new fermions. Furthermore, the decay rate of the heavier fermions into the DM is suppressed when the DM is mostly singlet-like. This, as already discussed, is a way to evade DD constrains. The only possible detection prospects, compatible with constraints from the DM phenomenology, would be represented by the detection of the charged fermion ψ ± , stable at the collider scales, giving large missing transverse momentum with charge tracks.

Majorana DM
The Majorana version of the Singlet-Doublet DM model is described by the following Lagrangian: where we have assumed the definition: As will be pointed out subsequently that the Singlet-Doublet model with a Majorana fermion DM can be interpreted as simplified limit of some specific SUSY realizations. To make this comparison easier, we have assigned to L L and L R fields the same set of quantum numbers as the higgsinos H d and H u (notice that this choice of quantum assignments is different from the case studied in the previous subsection). Further note that, given the transformation properties of the fields N , L L , L R under the discrete Z 2 symmetry, only one between the signs of the mass term M L and of the couplings y 1 and y 2 is physical [340]. We will assume, in our analysis, positive sign for M L and assign freely the signs of y 1, 2 .
Contrary to the previous case, the mass matrix of the electrically neutral new sector is now a 3 × 3 matrix, which is diagonalized through a unitary transformation leading to three (Majorana) mass eigenstates: the lightest of which, ψ 1 , is the DM candidate. The mass spectrum of the new states is again completed by a electrically charged Dirac fermion ψ ± with mass m ψ ± ≈ M L . The DM and the other neutral states interact with the SM states through s-channel mediation of the Z and Higgs bosons. The associated couplings are given by: where the labels V and A refer to vectorial and axial couplings, respectively. As expected, for Majorana fermions the coupling of the DM with the Z is only axial. The couplings with the Higgs boson instead, is given by: For the DM phenomenology, in particular for the DM annihilation, the couplings between the neutral states ψ i , their charged partner and the W ± -boson are relevant. These can be written as: We will trade in our numerical study the parameters y 1 , y 2 with y, θ, similar to Ref. [176], so that: These relations allow us to identify our setup as the limiting case of a SUSY theory in which the DM candidate neutralino is a bino-higgsino mixture once these substitutions are considered: with M 1 , μ and tan β being, respectively, the bino mass parameter, the mass parameter associated with a term bilinear in the up and down-type Higgs superfields and the ratio of the up and down-type Higgs VEVs. Interestingly, the couplings of the DM with the Z and Higgs bosons can manifest "blind spots", i.e., vanish for specific assignations of the parameters of the theory. The blind spot in the coupling of the DM with the Higgs can be found by rewriting the coupling as [173,340]: which, as evident, vanishes when: Given the chosen sign convention, the blind spot appears for sin 2θ < 0. A blind spot in the coupling of the DM with the Z boson is, instead, achieved, when |U 12 | 2 = |U 13 | 2 which corresponds to: tan θ = ±1, and/or M L = m ψ 1 .
The presence of these blind spots is particularly relevant for direct DM detection. Indeed, the two components, SI and SD, of the DM scattering cross-section on nuclei are sensitive to the couplings of the DM with Higgs boson and the Z boson, respectively. The corresponding cross-sections are given by: in the case of SI interactions, while SD cross-section, originated by the (axial) interactions with the Z boson, reads: As already discussed in the case of the Singlet-Doublet model with a Dirac fermion DM, the phenomenology related to the DM relic density features many similarities with the SM Z and Higgs portal models studied in Sect. 7. The DM annihilation rates into W + W − and Z Z/Zh final states are enhanced by t-channel exchange of the charged particle ψ ± , for the first, and by the additional neutral fermions, for the latter. A further enhancement can be provided by coannihilation processes. The annihilation cross-sections into W W and Z Z final states have sizable s-wave contributions, making them sensitive also to indirect signals (this would occur also in the case of a Dirac fermion DM. This possibility has not been explicitly discussed in the previous subsection because of the overwhelming constraints from DD). The annihilation cross-section into SM fermion pairs features also a s-wave contribution which is, however, helicity suppressed. Consequently, the latter is relevant only for the DM masses not very far from the one of the top-quark.
A first illustration of the DM constraints, for Singlet-Doublet model with a Majorana fermion DM, is provided in Fig. 27. We have reported there the constraints from the DM relic density and detection strategies in the bi-dimensional plane (M L , M N ). The coupling y has been set to 0.2, i.e., the value corresponding to the SUSY limit, while we have chosen two assignations for tan θ , namely ±2. Similar to the case of a Dirac fermion DM the correct DM relic density can be achieved only up to a maximal value of the DM mass, of approximately 1.1 TeV. This corresponds to the case of pure higgsino-like neutralino DM in supersymmetric models with the DM relic density mostly determined by annihilations into W W final state controlled by the gauge interactions. 38 The current limits from DD mostly affect the region M L ∼ M N which would correspond to maximal mixing between the singlet and the doublet components of the DM and then to maximal couplings with Higgs boson. 39 The anticipated sensitivities of XENON1T and LZ will, instead, allow to sensibly increase the coverage of the parameter space. Limits from indirect DM detection by Fermi show potential complementarity in the low DM mass region (see also [355]). We notice in particular the different shapes of DD contours according to the different signs of tan θ ; in the case of tan θ = −2, the presence of the blind spot corresponding to Eq. (135) is evident.
Similar to the case of a Dirac fermion DM, we have conducted a more extensive analysis, focused on the relation between the DM relic density and DD, by performing a scan over the four free parameters M N , M L , y, θ. The ranges chosen this time, are the following: The points featuring the correct DM relic density have been reported in Fig. 28 in the bi-dimensional plane (m ψ 1 , σ SI ψ 1 p ). The plot shows also the current and expected limits by 38 The maximal value of the DM mass is actually higher because of Sommerfeld enhancement [352][353][354]. For simplicity we have neglected this effect since it would marginally influence the discussion. 39 Figure 27 reports only limits from SI interactions. For the chosen assignation of the parameters, limits from SD interactions are always sub-dominant. For simplicity, thus, these have not been reported in the figure.  as the eventual effect of co-annihilations, a sizable number of points with viable DM relic density would survive even in the absence of signals at LZ.
We conclude by commenting briefly on the possible collider tests of the model under consideration. As already pointed out in the previous subsection, because of the sizable couplings with the SM gauge bosons, the charged fermion ψ ± and the neutral fermions ψ 2,3 can be produced through Drell-Yann processes, followed by successive decay into the DM and a (either on-shell or off-shell) W ± /Z boson, respectively, leading to events with 2 -3 leptons and missing energy. Given the similarity with SUSY setups, one could recast the results from the corresponding searches of chargino/neutralino production [290,356] in this framework. The present relevant collider constraints are, however, not competitive with the ones from other DM searches [176,337] and thus, for simplicity, have not been reported in Fig. 27.

Summary and discussion
We have discussed impact of the current and possible future DD limits, possibly complemented by the ones from ID and/or collider searches, in several simplified realizations of the WIMP DM.
The first and simplest classes of models considered are the ones in which the interactions of a pair of SM singlet scalar, fermionic and vectorial DMs and a pair of SM fermions, are mediated by electrically neutral s-channel (portal) mediators. In the minimal most case the particle spectrum of the SM should be complemented by just a new state, i.e., the DM candidate, since portal interactions can be mediated either by the SM-Higgs or by the Z-boson, although in the last case a theoretically consistent construction is more contrived. In the case of SM-Higgs portal, for all the DM spin assignations, SI interactions with nucleons are induced. The consequent very strong limits, due to the light mediator, are incompatible with the thermal relic density ad exception of DM masses above the TeV scale or the "pole", i.e., m DM m h /2, region. This last scenario would nevertheless be ruled out in the absence of signals at XENON1T (assuming a 2 years of exposure time) and LZ. In the Z -portal scenario current limits on the SI cross-section already exclude the pole region. These strong limits can nevertheless be partially overcome in two setups: (i) a fermionic DM with only axial couplings with the Z , as naturally realized in the case of a Majorana fermion DM and, (ii) a vectorial DM coupled through Chern-Simons term. In these two cases the DM features SD interactions with nuclei, whose constraints are sensitively weaker. In particular, in the case of a Majorana fermion DM, the thermal DM with mass of a few hundreds GeV would remain viable even in the absence of signals at the next generation detectors.
The SM-Higgs and Z-portal setups are easily extended to the cases of BSM spin-0 and spin-1 mediators, respectively. In the case of scalar mediators we have imposed, in order to preserve SU (2) L invariance, a Yukawa structure for the couplings of the mediator with the SM fermions. This, on one side, implies a suppression of the DM annihilation crosssection for masses below the one of the top-quark (unless the SS final state is kinematically accessible). At the same time, possible collider signals are also strongly suppressed so that the corresponding limits are not competitive with respect to the ones from DD and have been neglected for simplicity. Despite of the different velocity dependencies of the annihilation cross-sections, the regions of the correct DM relic density are then mostly determined by Yukawa structure of the couplings between the mediator and the SM fermions. The correct DM relic density is indeed obtained, far from the resonance regions, only when the tt and/or SS annihilation channels are kinematically open. Regarding DD, the limits are associated to SI component of the DM scattering crosssection for all the different assignations of the DM spin. The shape of the DD iso-contours are, however, different for the various DM scenarios. This is due to the different assignations of the couplings with the dimension of mass for the scalar and vectorial DM. Theoretical considerations suggest, indeed, to parametrized these couplings in terms of a fundamental mass scale, the mass of the mediator and the DM mass in the cases of scalar and vectorial DM, respectively, and an unknown dimensionless coupling. The current limits still allow masses of a few hundreds GeV for both the DM and the mediator while XENON1T, in the absence of signals after 2 years of exposure, will exclude mediator masses up to approximately 1 TeV and DM masses up to a few TeV. Given the several free parameters, for clarity of the picture, we have focused our investigation on the masses of the new particle states and fixed the couplings to be close to O (1) (see for alternative, e.g., Ref. [357]). We notice on the other hand that lowering the couplings would simultaneously suppress both the DD rate and the DM annihilation cross-section, in particular the SS channel becomes negligible as soon as the DM couplings λ S χ , g ψ , η S V deviate sensitively from O (1) values. As a consequence, in this setup, the thermal DM is achieved only in the pole region which requires particular fine tuning because of the typical small decay width of the scalar mediator.
The scenario of spin-1 BSM s-channel mediator is even more constrained than the spin-0 case. Indeed, the constraints from SI cross-section are typically much stronger, because of an effective enhancement of the cross-section due to the isospin violating interactions of the Z with nucleons, so that masses of the DM and the mediator approximately below 5 TeV are already excluded. In the case of no signals at the next generation DD experiments, the exclusion regions will extend up to masses of O ∼ 10 TeV, beyond the reach of LHC. In addition, the (reasonable) assumption of a Z coupled with both the SM quarks and leptons implies a strong complementarity with the LHC searches of dilepton resonances. The corresponding limits, exclude, for the models considered here, masses of the Z between 2 and 3 TeV (the exclusion can be even above 4 TeV in other realizations [358]), even in setups in which the SI component of the DD cross-section is suppressed or absent. We remark again that although in our analysis we have limited to some fixed assignations of the couplings, our results, nevertheless, have general validity because of the strong correlation between the DM relic density and its scattering rate of nucleons. For example, reducing the sizes of the couplings would actually reduce the viable parameter regions since the correct DM relic density would then be achieved only in correspondence of the s-channel resonances.
Despite of the fact that our work is focused on scenarios already probed by the current and will be tested by the near future DD experiments, we have nevertheless also discussed a setup in which DD is, in general, evaded: the pseudoscalar portal. Under the assumption of CP conservation only a fermionic DM is considered in this case. Most of the parameter space is substantially insensitive to DD (we remind here that we have, conservatively, considered values of the pseudoscalar mass above 1 GeV in order to avoid flavour constraints) since tree level interactions with nucleons are momentum suppressed and, furthermore, are not subject to coherent enhancement. A rather limited region of the parameter space might still be probed by 1-and multi-TON detectors because of an one-loop induced SI cross-section. The thermal DM is nevertheless sensitively constrained from ID. In addition, there is again a strong complementarity from the collider constraints, dominating for this scenario from monojet searches. A light pseudoscalar mediator can be interpreted as the pseudo-Goldstone boson of a spontaneously broken global U (1) symmetry. We have then considered the case of a complex scalar mediator which can be decomposed into a scalar and a light pseudoscalar components. Although in this case sizable DD limits are reintroduced, the thermal DM still remains viable in the large portions of the parameter space due to the presence of efficient annihilation processes in the aa and Sa final states.
Relaxing the hypothesis of a SM singlet BSM mediator and assigning it non-trivial quantum numbers under the SM gauge groups, one can construct a simple and predictive class of model which we have labeled as t-channel portal. In this case a single DM state is coupled with the mediator and a SM fermion, according to the gauge charge assignation of the mediator (for simplicity we have restricted most of our analysis to couplings with the right-handed up-type quarks). Contrary to the other scenarios considered in this work, here the DM pair annihilation occurs through t-channel exchange of the mediator while DD scattering is induced by its schannel exchange. Focusing, for simplicity, mostly to the case in which the mediator field has the same quantum numbers with respect to the SM gauge group as the right-handed up-type quarks, the scenarios, i.e., a complex scalar and a Dirac fermion DM, in which SI interactions are present are excluded for O(1) values of the couplings and for mediator masses up to O ∼ 10 TeV. On the contrary, thermal Majorana fermion DM is still viable for masses below a TeV and will be extensively probed by the next generation of DD experiments.
We have then performed some steps towards more theoretically motivated realizations of dark portals. As well known, the bilinears H † H and B μν (together with a new BSM field strength) are Lorentz and gauge invariant, so naturally lead to portal interactions with a dark sector, even if this is completely secluded. We have thus, considered the cases of (i) a scalar mediator coupled both to the DM and the SM-Higgs boson and mixes with the latter because of a non-zero VEV and, (ii) a Z coupled to the Z boson through a kinetic mixing term, also responsible for a mixing between the two spin-1 states. These two mixing in turn allow the BSM states to interact with the other SM states. Also concerning the coupling of the DM with the mediators we have considered less generic assignations with respect to the cases considered in the previous frameworks. In the case of SM-Higgs + spin-0 portal we have explicitly considered a dynamical origin for the DM mass. Indeed, a fermionic DM has been assumed to have Yukawa interaction so that its mass is originated by the VEV of the new scalar field. Similarly, a vectorial DM has been assumed to be the vector boson of a spontaneously broken dark U (1) gauge symmetry and its mass is again related to the VEV of the new scalar field, also charged with respect to the same U (1) group.
As the last case of study we have considered models in which the DM features SM gauge interactions, either being the lightest neutral component of a SU (2) L multiplet (in our studied examples we have focused on the case of SU (2) L doublets) or a mixing between the latter and a SM singlet. These models are characterized by the fact that the DM sector is composed of multiple states. Particularly relevant is the case in which some of these states are very close in mass to the DM. In this setup the DM annihilation cross-section features a twofold enhancement with respect to the other Dark portals. First of all the annihilation into W + W − is enhanced by the t-channel exchange of electroweakly charged states belonging to the DM sector (this, at the same time, requires that the DM has a non-negligible SU (2) component). A further enhancement comes from co-annihilation effects. Despite of a limited number of free parameters, these models are capable of encompassing many theoretically motivated scenarios like for example SUSY models. Interestingly, the Singlet-Doublet models considered at the end of Sect. 12 allows renormalizable couplings between the DM with the SM Higgs and the Z bosons. From the phenomenological point of view, the Singlet-Doublet model with a Dirac fermion DM features the greatest similarities with the SM Z-portal model, sharing with the latter the extremely strong constraints from DD. This model appears viable only in a very fine tuned coannihilation configuration. On the contrary, only the response of the future experimental facilities can fully probe the case of a Majorana fermion DM. More particular is the case of the Inert Doublet Model. Its main limitation is represented by the very efficient annihilation cross-section into gauge bosons. The viable DM relic density can be obtained when this annihilation channel is kinematically allowed and only very specific assignations of the masses of the new particle sector are considered.

Conclusions
We have reviewed the theoretical foundations of the WIMP paradigm and discussed the limits, prospect and challenges of direct DM detection in a multitude of models encompassing scalar, vectorial and fermionic DM setups. In the light of extensive programme of direct DM searches and including a broad variety of complementary probes from indirect and collider searches, we assessed the status of the WIMP paradigm in the context of simplified models, accounting for the current and projected limits.
In particular, we have reviewed well known portals such as the SM-Higgs portal and the Z-portal. We have also addressed the popular dark Z portal and many others models that possess in their spectrum more than one mediator. Moreover, we have also investigated new models dictated by the kinetic mixing, often used in dark photon models.
We concluded that the simplest constructions, i.e., the SM dark portals will be substantially ruled out, ad exception of the case of a fermionic DM with only axial couplings with the Z-boson (e.g., a Majorana fermion DM), in the absence of signals in the next generation of DD experiments.
The most straightforward extension of the SM dark portals, represented by the introduction of BSM s-channel mediators, are, similarly, strongly constrained in the presence of SI interactions of the DM with the nuclei. In particular, the case of spin-1 mediator is strongly disfavored because of the presence of complementary constraints from searches of resonances at the LHC, pushing the DM mass towards the multi-TeV scale.
The tension with DD constraints can be relaxed somehow in the next-to-minimal scenarios, featuring multiple mediators or new states lighter than the DM (we have reviewed the example of a light pseudoscalar).
In summary, we combined a plethora of experimental data set and theoretical models, computed the DM relic density, direct, indirect and collider observables to have a clear picture of where the WIMP paradigm stands and the future prospects. It is clear that most of the WIMP models will be scrutinized in the next decades, highlighting the paramount role of the next generation of experiments.
Moreover, our work shows that some DM constructions will survive the null results from the collider, direct and indirect experiments,suggesting that a further step in sensitivity reach is needed to falsify the WIMP paradigm.

Appendix A: Annihilation cross-sections
The thermally averaged pair annihilation cross-section is defined as: Away from resonances, a manageable analytical expression for the thermally averaged pair annihilation crosssection is obtained by performing the formal velocity expansion, in the non-relativistic limit, as defined in Ref. [23]. The thermally averaged cross-section can be computed as: where: This kind of integral can be computed analytically by considering an expansion in the series of for σ v lab , namely: σ v lab = a 0 + a 1 + a 2 2 + · · · . (A.4) We will report subsequently the full analytical expressions for the thermally averaged cross-section in the various dark portal scenarios investigated in this work, using the velocity expansion. Apart from some exceptions, we will adopt general notations for the couplings so that the various expression can be applied to the different models discussed in the main text as well as being of more general utility. In other words, we will use the generic notations g x xm , g yym , g xxyy for the relevant ones involving a generic DM particle "x", the mediator "m" and the SM particle "y", respectively (notice that couplings g χχ S , g χχ SS in these notations should not be confused with the ones defined in Eq. (90) with similar notations, for the case of SM-Higgs + spin-0 portal).
Further, note that we have used the symbols g Z Zh , g W W Z to parametrized the SM Z Zh, W + W − Z couplings, consistent with the aforesaid notations. We will adopt a slightly different notation for the couplings of a scalar (pseudoscalar) mediator S (a) with the SM fermions, which will be indicated as g S f (g a f ) (in this work we have considered the case g S(a) f = c S(a) m f v h ), and the couplings of the Z mediator with fermions (both the DM and SM), which will be indicated as g Z I i , I = V, A, i = ψ, f (contrary to what discussed in the main text, these couplings encode the gauge couplings in their definition). Notice that the couplings of scalar mediators with the scalar and vector fields are dimensional, rather than being decomposed, as in the main text, into a physical energy scale and a dimensionless parameter. In the case of spin-0 and spin-1 portals, the expressions for hh, Z Z and Zh final states can be derived directly, by suitable substitutions, from the ones corresponding to the SS, Z Z and Z h final states, respectively. The SM-Higgs + spin-0 portal and the kinetic mixing models feature, for the DM, also annihilation processes into hS and Z Z final states, respectively. We won't explicitly report the corresponding expressions since they are particularly complicated. For the same reason we also do not report analytical expressions for the relevant annihilation cross-sections of the models discussed in Sect. 12.
Finally, we remind that the velocity expansion is not valid in the vicinity of s-channel resonances and kinematic threshold of the annihilation channels [23,25]. ψψ → W + W − : ψψ → Z Z: ψψ → SS: (A.12)

Appendix A.5: t-channel portals
As already mentioned in the main text that in the case of O(1) coupling, the DM relic density is mainly accounted by the DM pair annihilations into the SM fermion pairs. We will, thus, explicitly report only the cross-sections corresponding to these processes.