How to save the WIMP: global analysis of a dark matter model with two s-channel mediators

A reliable comparison of different dark matter (DM) searches requires models that satisfy certain consistency requirements like gauge invariance and perturbative unitarity. As a well-motivated example, we study two-mediator DM (2MDM). The model is based on a spontaneously broken $U(1)'$ gauge symmetry and contains a Majorana DM particle as well as two $s$-channel mediators, one vector (the $Z'$) and one scalar (the dark Higgs). We perform a global scan over the parameters of the model assuming that the DM relic density is obtained by thermal freeze-out in the early Universe and imposing a large set of constraints: direct and indirect DM searches, monojet, dijet and dilepton searches at colliders, Higgs observables, electroweak precision tests and perturbative unitarity. We conclude that thermal DM is only allowed either close to an $s$-channel resonance or if at least one mediator is lighter than the DM particle. In these cases a thermal DM abundance can be obtained although DM couplings to the Standard Model are tiny. Interestingly, we find that vector-mediated DM-nucleon scattering leads to relevant constraints despite the velocity-suppressed cross section, and that indirect detection can be important if DM annihilations into both mediators are kinematically allowed.


Introduction
The idea that dark matter (DM) communicates with the Standard Model (SM) via the exchange of additional new particles (so-called dark mediators) has recently received large amounts of interest [1][2][3][4][5][6][7][8][9][10]. This framework allows for a consistent interpretation of simplified DM models [11][12][13][14][15][16][17][18] and thus for a reliable comparison of DM searches at the LHC [19][20][21][22][23][24] with other experimental constraints on the interactions of DM [25][26][27][28][29][30], as well as with the interaction strength required to obtain the correct relic density from thermal freezeout [31]. For a model with a single s-channel mediator the typical conclusion is that the couplings of the mediator to SM particles have to be rather small in order to be in agreement with experimental bounds and as a result DM overproduction can only be avoided in rather special corners of parameter space [4,6].
It has been pointed out, however, that for a meaningful comparison of different constraints the dark mediator model needs to be consistent with gauge invariance and perturbative unitarity [16]. As a straight-forward way to satisfy these requirements we consider a DM model containing a new U (1) gauge group and a dark Higgs that breaks the U (1) and generates the mass of the fermionic DM particle as well as the mass of the Z gauge boson. Both the Z and the dark Higgs can also couple to SM particles and thereby mediate the interactions of DM. For the Z such a coupling can arise if some or all of the SM fermions are charged under the U (1) gauge group, while the dark Higgs can couple to SM states via mixing with the SM Higgs.
The model we consider can therefore be thought of as a combination of two simplified models, one with a spin-1 s-channel mediator and one with a spin-0 s-channel mediator. Since we require all coupling structures and masses to result from a perturbative and gaugeinvariant UV completion, however, the model has to satisfy certain relations between the different masses and couplings not usually imposed on simplified DM models. In particular, these relations typically imply that it is not possible to completely decouple one of the two mediators while keeping the remaining masses and couplings fixed. We will therefore refer to this framework as two-mediator DM (2MDM). The aim of the present paper is to classify how the simultaneous presence of two mediators changes the conclusions derived for a single mediator.
In the specific example model for 2MDM that we study in detail in this paper, we assume vector-like and flavour-diagonal couplings of the Z to SM quarks. Hence the U (1) can be identified with gauged baryon number [32]. Recently, realistic gauge theories for baryon number that are in agreement with all experimental constraints have been discussed in [33,34]. These models contain a fermionic DM candidate [35][36][37][38], and in the limit of neglecting the additional fermions relevant for anomaly cancellation they reduce to the example that we consider.
The generic prediction is that both mediators contribute to the DM annihilation cross section that sets the DM relic abundance. On the other hand, certain experimental bounds may constrain only one of the two mediators. For example, the observed properties of the SM Higgs boson constrain the couplings of the spin-0 mediator, whereas LHC searches for monojets and dijets are more sensitive to the properties of the spin-1 mediator. This compelling interplay provides the model with additional freedom to evade experimental constraints and therefore reduces the tension with the requirement to reproduce the observed relic abundance.
A particularly interesting situation occurs if only one of the two mediators has sizeable interactions with SM states, while the other one is very weakly coupled to the SM. In this case, the weakly-coupled mediator can be much lighter than the DM particle and therefore provide a new final state for DM annihilation. This configuration opens up new parameter space where the DM relic abundance can be reproduced. Such a configuration has been called "secluded WIMP" [39] or "indirect Higgs portal" [40] in the context of Higgs portal DM. We refer to such a weakly-coupled light mediator as a dark terminator, because it terminates rather than mediates the interactions of DM.
Models containing only a DM particle and a dark terminator are notoriously difficult to constrain due to the smallness of the interactions between DM and SM particles. The situation is different in our context, because the second mediator can still have sizeable interactions, which can be probed by a range of different DM searches. Hence, the dark terminator can be thought of as a tool to relax the relic density constraints on the interactions of the second mediator in a controlled way and to study the phenomenology of the second mediator in the extended viable parameter space.
The presence of a dark terminator can also lead to novel experimental signatures. Two particularly intriguing examples are indirect detection signals from cascade annihilations, which can avoid the p-wave suppression typically present for Majorana DM [7,41], and production of dark terminators at the LHC from final state radiation of other dark sector states [42][43][44][45][46]. This paper is structured as follows. In section 2 we introduce our model and discuss the relations between the various couplings. In section 3 we then consider the case where the interactions between DM and SM particles are dominantly mediated by the Z and the dark Higgs can be a terminator. Section 4 focuses on the converse case, i.e. a dark Higgs mediator and a Z terminator. Then, in section 5 we study additional effects that may be important if both mediators have sizeable couplings to SM states, or if both act as dark terminators. In section 6 we present the results of a global scan of the parameters of the model, with the goal to identify the regions of parameter space which are consistent with a thermal relic abundance, all experimental constraints and perturbative unitarity. We summarise in section 7. Supplementary material is provided in four appendices.

Model set-up
The aim of this paper is to study the interaction of a Majorana DM particle χ with two different mediators, namely a massive vector boson Z and a real scalar s. We consider a set-up where the SM gauge group is extended by an additional U (1) that is broken by the vacuum expectation value (vev) w of a complex scalar field. In such a set-up the real scalar s corresponds to the Higgs boson from the spontaneous breaking of the U (1) and the Z is the additional gauge boson that acquires a mass proportional to w. This set-up offers a natural framework for the simultaneous presence of both mediators and is motivated by basic requirements such as gauge invariance and perturbative unitarity [16].
We consider the following interactions of the two mediators with DM and SM quarks: χs , (2.1) L q ⊃ − q g qq γ µ qZ µ + sin θ m q vq qs , (2.2) where v = 246 GeV is the electroweak vev and θ is the mixing angle between the dark and the SM Higgs. Here we only describe the most important features of our example model for 2MDM, more details are given in appendix A. Note that we restrict our discussion to vector couplings of the Z to quarks. In principle it is possible that the Z has additional couplings to the SM, such as axial-vector couplings to quarks, vector and/or axial-vector couplings to leptons, or kinetic mixing with the SM hypercharge gauge bosons. These alternative possibilities are discussed in some detail in appendix C. In general they are strongly constrained by electroweak precision tests (EWPT) and dilepton searches at the LHC [16]. We show in appendix C that these constraints typically imply that neither axial couplings nor kinetic mixing can give large enough DM annihilation cross sections to avoid DM overproduction. To avoid constraints from flavour-changing processes, we assume the vector couplings to quarks to be flavour-independent. The flavour-universal vector couplings to quarks that we consider correspond to the case of gauged baryon number, and a detailed discussion of anomaly cancellation and full viable solutions can be found in [33,34,47,48]. Note that a Majorana DM can easily be realised in these models by an appropriate choice of quantum numbers, with different possibilities for the additional fermion content [37,38]. Let us emphasise that our model is still simplified in the sense that we do not specify the additional particles needed for anomaly cancellation, hence our example model corresponds to a subsector of these models. The coupling structure imposed by the underlying symmetries of the model ensure that there is no colour anomaly [16]. As a result, no new coloured states are required in order to make the U (1) anomaly-free and the additional states are expected not to affect the phenomenology of the model significantly. 1 At first sight, the interactions in eqs. (2.1) and (2.2) resemble the ones usually considered for s-channel simplified models. However, in the present context, there are a number of important differences. Most importantly, the two couplings of DM to the Z and the dark Higgs are not independent, because the DM Yukawa coupling can be re-expressed in terms of the DM mass and the vev of the dark Higgs, y χ = √ 2 m χ /w, and w can be reexpressed in terms of the Z mass and the Z -DM coupling, w = m Z /(2g χ ). We therefore obtain the important equality This relation is essential to ensure that processes like χχ → Z L Z L , where Z L denotes a Z with longitudinal polarisation, do not violate unitarity at large energies [16]. From the particle masses coupling constants DM mass m χ dark-sector coupling g χ or y χ Z mass m Z quark-Z coupling g q dark Higgs mass m s Higgs mixing angle θ Table 1. Summary of the 6 independent parameters of our model. For given masses, there is only one independent dark-sector coupling, since g χ and y χ are related via eq. (2.3).
practical point of view, having only one independent coupling in the dark sector significantly simplifies the analysis. The interactions between DM and SM states are therefore characterised by 6 independent parameters: the three masses m χ , m Z and m s , the two couplings between the mediators and the SM, g q and sin θ, and the single dark sector coupling, either g χ or y χ . This number should be compared to the four parameters conventionally required in simplified models with a single s-channel mediator. The additional two parameters provide the freedom to interpolate between regimes where one of the mediators dominates the phenomenology and to describe the case that both mediators give a relevant contribution. The parameters are summarised in table 1.
In principle, it would be conceivable to simply implement the interactions given above, without specifying the interactions of the new particles with any other SM state. For the purpose of this paper, however, we prefer a different approach, where we include the additional contributions predicted by the simple UV completion in terms of a spontaneously broken U (1) , see appendix A for details. The additional effects compared to the simplified model introduced in eqs. (2.1) and (2.2) can be understood as follows: 1. An essential feature of our model is that the two mediators are linked by their common origin from the broken U (1) . This implies in particular that the two mediators can interact with each other, leading to processes like for example χχ → Z * → Z s or χχ → s ( * ) → Z Z .
2. The interactions between the dark Higgs s and SM quarks arise from mixing between the two Higgs bosons. Such a mixing necessarily introduces couplings of all SM particles to the dark Higgs proportional to their masses, leading in particular to the possibility that DM particles can annihilate into SM gauge bosons. 2 At the same time, this mixing changes the couplings of the SM Higgs to other SM particles and furthermore couples the SM Higgs to both the Z and the DM particle. This leads to additional constraints on sin θ from the observed properties of the SM Higgs boson.
3. In order to ensure that the model remains perturbative, one needs to consider not only g q and g χ , but also y χ and sin θ. In particular, requiring that all interactions between the two Higgs bosons remain perturbative leads to non-trivial bounds on g q sin θ g q ∼ sin θ sin θ g q m s m Z Spin-1 mediator Spin-0 mediator with simplified model spin-1 terminator Spin-1 mediator with Spin-0 mediator spin-0 terminator simplified model Table 2. Simple illustration of the different regimes of the two-mediator model considered in this paper. Note that this table ignores the mass of the DM particle, which ultimately determines whether one of the two mediators can act as a dark terminator.
the coefficients in the scalar potential, which can be translated into bounds on the mixing angle θ in terms of the various masses, see section 2.3 below.
4. In this model the stability of the DM particle is a consequence of the U (1) gauge symmetry. Even after symmetry breaking a Z 2 symmetry remains unbroken, which protects the DM particle from decay.
5. Finally, a general expectation for such a model is that the new U (1) can mix with the U (1) Y hypercharge gauge group of the SM. In principle, the kinetic mixing parameter could be taken as an additional free parameter of our model. However, the magnitude of this mixing is very tightly constrained, so that it typically cannot have a large effect on the DM phenomenology of the model. We therefore assume in the main text that kinetic mixing is absent at some high scale and only introduced radiatively by quark loops. Tree-level kinetic mixing is discussed in appendix C.
We structure the discussion by considering the couplings between the two mediators and SM states. For g q sin θ, we expect the Z to play the dominant role in the phenomenology of the model (discussed in section 3 below). Conversely, for sin θ g q , the dark Higgs will be responsible for mediating the interactions of DM (section 4). Finally, if both couplings are comparable, we can expect an interesting interplay between the two mediators, depending on the ratio of the different masses (section 5). This set-up is illustrated in table 2.

The relic density constraint
Throughout this work we will adopt the hypothesis that the total observed abundance of DM is produced by thermal freeze-out ("WIMP hypothesis"). This implies that for each point in the parameter space all available annihilation channels are completely determined by the structure of the model and included self-consistently. We provide analytic approximations for the most relevant annihilation cross sections in appendix D. For the numerical study we have implemented the model in MicrOMEGAs v4.2.5 [49] in order to obtain an accurate calculation of the predicted relic density. 3 Then we vary the parameters of the model to match the observed relic density Ω DM h 2 = 0.1188 ± 0.0010 [50].

Perturbative unitarity
Partial wave perturbative unitarity can be used to derive a number of conditions that the couplings and masses in a given model have to satisfy in order to obtain a consistent description [16]. Considering the process χχ → χχ in the limit of large centre-of-mass energy √ s m χ and neglecting terms proportional to log(s/m 2 Z ) one gets the following unitarity bounds: Since g χ and y χ are related via eq. (2.3), those inequalities imply also a constraint on the masses. For example the second inequality can be rewritten as g χ m χ /m Z < √ π. Perturbative unitarity also constrains the couplings λ s , λ h and λ hs appearing in the scalar potential, which is fully determined once the dark Higgs mass m s , the mixing parameter sin θ and the dark vev w have been specified (see appendix A). By considering the scattering processes ss → ss and hh → hh, it is possible to derive a combined constraint on the three couplings appearing in the scalar potential: 4 Expressions for the scalar couplings in terms of the free parameters considered in this section are given in eqs. (A.25)-(A.27). In the absence of Higgs mixing (λ hs = 0), the inequality (2.5) simply gives an upper bound on the mass of the dark Higgs: m s < 16π/3 w = 4π/3 m Z /g χ . We note that this expression can also be rewritten as y χ < 32π/3 m χ /m s , thus giving an upper bound on y χ for fixed ratio m χ /m s . For m s /m χ > 4π/3 ≈ 2.05 this bound is stronger than the bound on the Yukawa coupling given in eq. (2.4).
For the numerical analysis below we are going to impose the conditions specified in eqs. (2.4) and (2.5) to determine the region of perturbative unitarity. Note that we do not consider the running of the dark gauge coupling, which may become relevant for g χ close to the perturbative boundary. For such large values of g χ it is expected that the dark gauge coupling develops a Landau pole at some high energy. Requiring that g χ remains perturbative up to high scales would therefore lead to slightly stronger perturbativity bounds. For the largest part of the parameter space of the model, however, the running is sufficiently slow so that these effects are completely negligible. 3 MicrOMEGAs does not take non-perturbative effects due to the multiple exchange of light mediators into account. A detailed study of these effects is beyond the scope of this paper. We note that while such effects might be relevant in some parts of the parameter space, we do not expect these effects to significantly change our conclusions. 4 Our final result differs from the one in [51] by a factor of 2, because we find a factor 1/2 from phase space for two identical final-state particles (see also [52]). Note also that stronger bounds could be obtained by also considering the scattering of W + W − , ZZ and Z Z [53].
We begin by focusing on spin-1 mediation, i.e. we consider the case that the Higgs mixing angle θ is sufficiently small to be negligible for the relic density calculation. If the dark Higgs is significantly heavier than the DM particle, it plays a negligible role in the thermal freeze-out. This scenario has been considered to some extent in [16]. Here we extend this framework by considering also the situation of m s m χ . We start by summarising the most important constraints on the model in sections 3.1 and 3.2, and discuss our results in section 3.3.

Direct detection
For our assumptions of Majorana DM and vector couplings between the Z and quarks we obtain an effective four-fermion interaction relevant for DM-nucleus scattering of the axial-vector type: where we have integrated out the Z . For this operator, the DM-nucleon scattering cross section is suppressed in the non-relativistic limit. In most studies the resulting constraints from direct detection experiments have therefore simply been neglected. Here we show, however, that these constraints can in fact be relevant. The leading non-relativistic DM-nucleus operators have been classified systematically [55]. It was shown that in the non-relativistic limit the operator (3.1) can be decomposed into a piece proportional to the DM velocity v and a piece proportional to the momentum transfer q. In the notation of [56] one finds where v ⊥ = v + q 2µ χN with µ χN being the DM-nucleon reduced mass, and S N and S χ denote the spin of the nucleus and the DM particle, respectively. This corresponds to the operators O 8 and O 9 of [56]. Crucially, the first term is independent of the nucleus spin and therefore receives a coherent enhancement proportional to the mass of the target nucleus squared. This enhancement factor, which can be as large as 2×10 4 for xenon-based experiments, partially compensates for the velocity suppression and allows direct detection experiments to retain some sensitivity to the interaction in eq. (3.1).
Since these operators lead to a recoil spectrum that differs substantially from the standard spin-(in)dependent interactions, the limits reported by direct detection experiments do not apply directly. Currently LUX [57] has the best sensitivity, and we translate their results into an upper bound on the interactions considered here.
As we do not have access to detailed information on the observed events or the background expectation, we cannot repeat the full likelihood analysis performed by the collaboration. In our simplified analysis we exploit the fact that most of the events observed by LUX fall above the median in the S1-S2 plane (see Fig. 2 of [57]) and limit the region of interest to the lower half of the LUX search region. We extract the detection efficiency as a function of the nuclear recoil energy from Fig. 1 of [57]. Given that the efficiency is close to unity for a substantial range of nuclear recoils, the smaller signal region can be accounted for by reducing the total efficiency by a factor 1 2 . We use the implementation of the nuclear response functions provided by [56] to calculate the differential recoil rate in a xenon detector and, after convolving the differential rate with the efficiency of LUX, we obtain the number of expected events in the signal region. Taking one observed event below the median in the S1-S2 plane and assuming no background (which gives the weakest upper limit on the signal) we find a 90% CL upper limit on the expected number of events of 3.89. In order to validate our results we repeat our analysis with the LUX 2013 data [58] and recover the bound on the effective operator in eq. (3.1) reported in [59].

Collider constraints
Important collider constraints on the model arise from a number of searches at the LHC and at LEP. Some of the signatures are directly related to the production of DM, such as monojet searches, while other observables, for example searches for dijet and dilepton resonances, constrain new physics interactions between SM particles. Finally, indirect precision measurements of SM relations, in particular EWPT, can be relevant.
Monojets. Searches for the production of DM in association with jets have been conducted by ATLAS and CMS both at 8 and 13 TeV [60][61][62][63]. Since the 13 TeV results are not yet competitive with the final limits from the 8 TeV run of the LHC we focus on the 8 TeV CMS results [60] in our analysis. We simulate monojet events with CalcHEP v3 [64] and pass them to Pythia v8 [65] for showering and hadronisation. Detector effects are taken into account with the help of DELPHES v3 [66].

Dijets.
A combination of dijet searches from ATLAS and CMS at 8 TeV and 13 TeV [67][68][69][70][71] was recently performed in [5]. The resulting model-independent bounds on the Z coupling as a function of its mass and width can be directly applied to our case. These bounds are valid as long as Γ Z /m Z ≤ 0.3, which is sufficient for all coupling combinations considered in this paper (with the exception of the global scans discussed in section 6). While larger widths can in principle be probed using dijet angular correlations (see [70]), we do not consider these constraints here. We also note that there are presently no LHC bounds on dijet resonances with an invariant mass below 500 GeV. In principle, this mass range can be constrained using dijet searches from previous hadron colliders and LHC searches for dijet resonances produced in association with SM gauge bosons [4], but we find that stronger constraints are obtained from the searches for dilepton resonances discussed next.
Bounds due to gauge boson kinetic mixing. It has been known for a long time that fermions charged under hypercharge and a U (1) induce kinetic mixing between the corresponding gauge bosons [72] which can be parametrized by where B µν and F µν are the field strength tensors of the SM hypercharge U (1) Y and the U (1) , respectively. This interaction respects the full gauge symmetry of the theory and could therefore be present at tree level so that would be a free parameter. This option is discussed in appendix C, where we show that it is strongly constrained. Here we adopt the natural assumption that kinetic mixing is absent at some high scale (see e.g. [73]) and consider only the contribution expected from loops. Under the assumption that vanishes at a scale Λ the induced kinetic mixing at a scale µ is given by [74] ( Hence, in this approach there is a natural size for that depends only logarithmically on Λ. To estimate the magnitude of this effect, we set Λ = 10 TeV and identify µ with m Z in the following. Effectively, is then no longer a free parameter. As detailed in appendix B, the mixing between the Z and the Z generates a coupling between leptons and the Z which can be searched for using dilepton resonances. In our numerical analysis we include the 8 TeV ATLAS search for dilepton resonances [75] and a CDF search [76], which is more sensitive in the low-mass regime. Furthermore, kinetic mixing also modifies the properties of the Z, see appendix B for details. In the gauge boson sector any deviation from the SM expectation is tightly constrained by precision measurements at LEP. The impact of heavy new physics is conveniently parametrized by the S and T parameters. For m 2 Z > 2m 2 Z we confront the model predictions with the limits from the combined fit of S and T given in [77]. For lighter m Z this procedure becomes unreliable and we use the limit on the ρ parameter ρ ≡ m 2 W /(m 2 Z c 2 W ) = 1.0004 ± 0.00024 [77] instead. Our bound on is in good agreement with the limit obtained from a general fit to precision observables presented in [78]. We emphasise that the bounds from dilepton searches as well as EWPT could be modified due to the UV sensitivity of the kinetic mixing. The resulting constraints are discussed in more detail in appendix C.

Results
All DM annihilation processes involve the single dark sector coupling g χ . As a result, it is always possible to fix this coupling (for given masses m Z , m χ and m s and fixed coupling g q ) by the requirement that the observed relic abundance is reproduced (up to unitarity constraints). This reduces the dimensionality of the parameter space and thus we can get a better qualitative understanding of the relevance of the various constraints. Figure 1 shows the value of g χ obtained from the relic density requirement as a function of m Z (top row) and of m χ (bottom row) for different values of the remaining masses and couplings. For the dark Higgs boson mass we consider two representative cases: the solid black line corresponds to the case that the dark Higgs is significantly heavier than the DM particle (m s = 3 m χ ), whereas the dashed black line corresponds to the case that it is much lighter (m s = 0.1 m χ ). It can be seen that away from the resonance corresponding to m χ ≈ m Z /2 smaller values of g χ are required in case of a light dark Higgs. This is because  DM can very effectively annihilate into final states with dark terminators, e.g. χχ → ss, as will be discussed in more detail below.
We also show in figure 1 the constraints discussed above both in the m Z -g χ and in the m χ -g χ parameter plane. Here and in the following, all constraints are shown at 95% confidence level, except for direct and indirect detection constraints, which are shown at 90% confidence level as customary. As we consider negligible mixing of the dark Higgs with the SM Higgs, these constraints are independent of the value of m s , with the only exception of the unitarity constraint, which is indicated by a solid grey line for m s = 3 m χ and by a dashed grey line for m s = 0.1 m χ . For m s = 3 m χ and θ ≈ 0 the unitarity bound from the scalar potential given in eq. (2.5) is always stronger than the direct bound on the Yukawa coupling in eq. (2.4) and leads to the strongest constraint for m Z < 500 GeV (top row) and m χ > 40 GeV (bottom row). For larger values of m Z or smaller values of m χ this bound becomes weaker, so that the unitarity bound results simply from the requirement g χ < √ 4π. For m s = 0.1 m χ the requirement of perturbative unitarity of the scalar potential is less constraining and the relevant bound for small m Z and large m χ comes from the condition y χ < √ 8π. In the lower panels of figure 1 the required thermal coupling g χ , as well as the unitarity bound on g χ , are seen to decrease with increasing m χ in a seemingly counter-intuitive manner. This behaviour can be understood once we consider the dominant annihilation channels. For m Z m χ the relic abundance is either set by χχ → Z L Z L or by χχ → sZ L and both these processes are controlled by the Yukawa coupling y χ ∝ g χ mχ m Z which grows with m χ just as expected. We note that for m s m χ the DM mass can be as large as m χ ∼ 40-50 TeV before the requirement of perturbative unitarity becomes incompatible with thermal freeze-out [79].
We now fix the dark sector coupling via the relic density requirement as shown in figure 1. We can then study the phenomenology of the model in the m χ -m Z parameter plane. The resulting constraints are shown in figure 2. The two rows correspond to different values of g q , the two columns to different ratios of m s /m χ . The constraints can be understood qualitatively by referring to figure 1: as expected, the weakest constraints are found for the resonance region m Z ∼ 2 m χ , where g χ can be very small. For increasing m Z larger values of g χ are required to reproduce the relic abundance up to the point where g χ is in conflict with the requirement of perturbative unitarity. For small Z masses, on the other hand, constraints from direct detection and EWPT are very strong. Dijet and dilepton constraints are most relevant for 100 GeV ≤ m Z ≤ 3 TeV. Finally, monojet constraints dominate for m χ < m Z /2, where DM production at the LHC receives a resonant enhancement.
In the left column, we have taken m s = 3 m χ . For such a large mass (and the assumed small mixing), the dark Higgs remains completely irrelevant for the phenomenology of the model. Its presence becomes nevertheless visible in the fact that large DM masses are excluded by the requirement of perturbative unitarity of the scalar potential. In other words, if we insist on m s m χ , the dark Higgs self-coupling necessarily becomes nonperturbative for DM masses at the TeV scale. This constraint can be significantly relaxed if we instead consider the case m s = 0.1 m χ , as shown in the right column of figure 2. In this case, the bound on y χ from perturbative unitarity only requires m χ 40-50 TeV (see figure 1).
Taking the Higgs to be light compared to the DM particle leads to additional changes, which are visible even for DM masses well below the unitarity bound. The reason is that, once the dark Higgs mass becomes comparable to or lighter than the DM mass, DM particles can annihilate into it. These processes can happen even for very small Higgs mixing, as long as the dark Higgs ultimately decays into SM particles. The two processes of interest are χχ → ss and χχ → sZ . The corresponding annihilation cross sections It is important to note that the dominant contribution to the annihilation into two dark Higgs bosons results from the diagram with a t-channel DM particle, so that the resulting annihilation cross section is proportional to y 4 Consequently, the requirement that g χ remains perturbative still places an upper bound on m Z even if the dark Higgs is the dominant annihilation channel.

Spin-0 mediation
Let us now consider the case g q 1, such that the dark Higgs is the only state that can mediate sizeable interactions between DM and SM particles. If the Z couples only very weakly to SM states, it also makes sense to assume that kinetic mixing is negligible. The independent parameters of importance for this case are therefore the three masses m Z , m χ and m s , the coupling between DM and the dark Higgs y χ , and the Higgs mixing parameter sin θ. Note that in contrast to the previous section, we now take y χ as the independent dark sector coupling, which is related to g χ via eq. (2.3). We now review the relevant constraints for this case.

Higgs mixing and signal strength
The magnitude of the mixing angle θ between the SM Higgs and the dark Higgs is constrained by the observed Higgs signal strength µ, which is required to satisfy µ > 0.89 at 95% confidence level [80]. In our model, there are three separate effects that potentially lead to a reduction of the signal strength: • Mixing between the two Higgs bosons reduces the production cross section for the SM Higgs.
• For m χ < m h /2, the SM Higgs can have a sizeable invisible branching fraction. In principle, such invisible decays can be tested directly by looking for invisibly decaying Higgs bosons produced either in vector boson fusion or in association with a Z boson. The resulting constraints are however generally weaker than the ones obtained from the Higgs signal strength.
• For m s < m h /2 or m Z < m h /2, the SM Higgs can decay into two singlets or two Z bosons, which depletes the cross section in the remaining channels. Such decays would give rise to decays of the form h → 4f (with f being SM quarks or leptons), which should in principle be observable. A detailed study of such decay processes is beyond the scope of this work, so we only consider the impact on the Higgs signal strength here.
The partial decay width for invisible decays is given by (assuming m χ < m h /2) while the partial decay width for decays into a pair of dark Higgs bosons with m s < m h /2 is Note that Γ ss is sensitive to the sign of θ due to an interference between contributions proportional to 1/w and 1/v. A negative sign of θ leads to a partial cancellation between these contributions and thus suppresses the resulting bounds. The partial decay width for decays into a pair of Z bosons with m Z < m h /2 is given by Finally, if Γ 0 denotes the total width of the SM Higgs for cos θ = 1, the partial width into all SM final states will be given by Furthermore, for cos θ < 1, the Higgs production cross section will be reduced proportional to cos 2 θ. Putting all of these contributions together, we find that the predicted Higgs signal strength in our model is given by Even for Γ ss = Γ Z Z = Γ inv = 0, the observed lower bound on µ [80] implies θ < 0.34. If additional decay channels are kinematically allowed, even stronger bounds will be obtained. The mixing angle θ is also constrained by resonance searches for additional Higgs bosons. These bounds depend on the mass of the dark Higgs and can be somewhat more stringent than the mass-independent bound of θ < 0.34 for 200 GeV m s 500 GeV [81,82]. However, for the values of θ considered in the remainder of this paper these bounds play no role. A similar comment applies to bounds on θ from EWPT, see e.g. [83].

Direct detection
The scalar mediator induces unsuppressed spin-independent DM-nucleus interactions, which are strongly constrained by direct detection experiments. Hence these experiments provide a stringent bound on the Higgs mixing angle θ. The spin-independent DM-nucleon scattering cross section is given by where µ χN = m N m χ /(m N + m χ ) is the reduced mass of the DM-nucleon system and f N = 0.3 is the effective coupling of DM to nucleons [18]. Recent results from the LUX experiment [57] place strong bounds on σ SI , which lead to constraints on θ that can be comparable or even more constraining than the ones from the Higgs signal strength. We implement this constraint by comparing our predicted value of σ SI to the 90% CL upper bound from LUX [57].

Results
Let us first consider the case that the Z is significantly heavier than the DM particle. Since we consider very small quark couplings, the Z then plays no role for the phenomenology of the model. In the presence of Higgs mixing, however, both the dark Higgs and the SM Higgs can mediate DM annihilation processes (see appendix D). Moreover, for m χ > m s , the DM particles can also annihilate into pairs of dark Higgs bosons, which subsequently decay into other SM particles (see section 3). Crucially, all these annihilation channels depend on y χ , so it is always possible (for given masses and mixing angle) to fix the DM Yukawa coupling in such a way that the observed DM relic abundance is reproduced (up to perturbativity). The required value of y χ can then be used to compare to the various constraints discussed above.
The results of such an analysis are shown in the left column of figure 3 as a function of m χ and m s for different values of θ (and m Z large enough that the precise value is irrelevant). The following features can be identified: For m χ ∼ m h /2, annihilation via the SM Higgs receives a resonant enhancement, whereas for m χ ∼ m s /2 annihilation via the dark Higgs is resonantly enhanced. For m χ > m s direct annihilation into two dark Higgs bosons dominates. In this parameter region, the relic density is independent of sin θ and therefore the constraints from direct detection and the Higgs signal strength can be weakened arbitrarily by making sin θ very small. Nevertheless, for m Z = 3m χ the required value of g χ may violate perturbative unitarity, limiting the allowed parameter region to m χ < 3.5 TeV. For small values of m χ , on the other hand, large regions of parameter space violate the perturbativity condition on y χ . The requirement of perturbative unitarity in the scalar potential only becomes relevant if both m s and m χ are large.
Let us now turn to the case where the Z is significantly lighter than the DM particle, so that DM can directly annihilate into a pair of Z bosons. Since we assume that the quark coupling of the Z is very small, these Z bosons can have a very long lifetime, but it is easily possible to ensure that they decay before Big Bang nucleosynthesis. As shown in appendix D, the relic density becomes independent of m Z in the limit m Z → 0, so we simply fix m Z = 0.1 m χ for convenience. In this case, the bound g χ < √ 4π can easily be satisfied up to very large DM masses, so only the bound on y χ and the requirement of a perturbative scalar potential remain visible.
In the case that 2 m χ > m Z + m s , one also needs to consider to the process χχ → sZ . In contrast to all other processes considered so far, this annihilation process can proceed via s-wave and is enhanced for m χ m Z due to the DM Yukawa coupling becoming large. In analogy, also the process χχ → hZ can become large in the presence of Higgs mixing. It is then not even necessary for the dark Higgs to be light.
As shown in the right column of figure 3, having a light spin-1 terminator significantly relaxes the constraints. In particular, since smaller values of y χ are sufficient to reproduce the observed relic abundance, the requirement of perturbative unitarity is less constraining. Nevertheless, it is still not possible to raise m s arbitrarily above m χ . If the channel χχ → sZ is open, indirect detection experiments give relevant constraints for m χ 100 GeV, which are independent of both θ and g q (see section 5.1). The excluded parameter region is shown in yellow in figure 3.
To conclude this section, we return to the choice of sign for θ. As discussed above, this sign is relevant for the trilinear vertices coupling one SM Higgs to two dark Higgs and vice versa. Considering θ < 0 thus modifies the prediction for the decay h → ss as well as for the annihilation process χχ → hh via an s-channel dark Higgs. The latter process, however, is typically subdominant compared to χχ → W + W − , so that the relic density calculation is not significantly affected by the sign of θ (this is also true for the process χχ → hs). To make the impact of the sign choice for θ explicit, we show in figure 4 the case θ = −0.2 for m Z = 3 m χ and m Z = 0.1 m χ , which can be directly compared to the first row of figure 3. We observe that the only visible change is that the bound from the Higgs signal strength for m s < m h /2 is significantly relaxed for the negative sign. Nevertheless, this parameter region is independently excluded by direct detection experiments (which are not sensitive to the sign of θ). For smaller values of |θ| the differences between positive and negative sign are even smaller. Since considering θ < 0 does not open up any additional parameter space, we will focus on the case θ > 0 for the remainder of this work.

Two mediators
Having discussed the case of a Z mediator with g q sin θ in section 3 and the case of a dark Higgs mediator with sin θ g q in section 4, we now turn to the case where the couplings of the two mediators are comparable. To visualise this situation we consider the m s -m Z plane, i.e. we vary the masses of both mediators, while keeping the DM mass fixed.
We will be particularly interested in the case that both g q and sin θ are very small (typically of order 0.01 or smaller), so that both mediators are secluded from the SM. In this case, most experimental constraints can be suppressed. However, as recently pointed out in [7], it may still be possible to probe this set-up in indirect detection experiments, which we now discuss in more detail.

Indirect detection
For almost all annihilation channels discussed above the relic density is dominantly set via p-wave processes. An exception is the s-wave contribution from χχ → Z Z which is comparable to the p-wave one for 0.5 m χ < m Z < 0.9 m χ . For smaller m Z the p-wave contribution is enhanced by m 4 χ /m 4 Z and dominates the cross section during freeze-out.
As a result, the annihilation rate in the present Universe is suppressed compared to the thermal cross section. 6 However, one interesting possibility is the case that both the Z and the dark Higgs are sufficiently light that the annihilation process χχ → Z s is allowed. This process proceeds dominantly via s-wave, with an annihilation cross section given by (for θ → 0 and m Z , m s m χ ) When kinematically allowed, this process is typically the dominant one for thermal freezeout. In this case, observable indirect detection signals may be obtained from cascade annihilations [41], i.e. the subsequent decays of the Z and the dark Higgs into SM particles.
Fermi Large Area Telescope (LAT) observations of Milky Way dwarf spheroidals [84] place tight constraints on the γ-ray flux from DM annihilations. To study these constraints, we calculate this γ-ray spectrum as a function of g χ and the three masses m χ , m Z and m s using MicrOMEGAs v4.2.5 [49]. For each set of parameters we then calculate the likelihood using the publicly available likelihood functions from Fermi-LAT [84]. Finally, we combine the likelihoods from 15 different dwarf spheroidals using the J-factors from [84].
We find that for m Z , m s < m χ 100 GeV, the resulting constraints are typically sensitive to the thermal cross section, so that the possibility that the DM relic abundance is set by the process χχ → Z s in this mass region can be excluded.

Results for fixed couplings
For the remainder of this section, we proceed in the same way as above by fixing the couplings θ and g q and then calculating g χ as a function of the three masses by imposing the relic density requirement. As expected, the various constraints are found to depend sensitively on the chosen values of the couplings. In section 6 we therefore perform a global scan, where for each combination of masses we vary all three couplings simultaneously in order to determine whether there is a combination that reproduces the observed relic density while evading all experimental constraints.
We first consider the case where both g q and sin θ are sizeable, i.e. of order 0.1. As discussed above, significantly larger values of g q would be in strong tension with searches for dijet resonances. Significantly larger values of θ, on the other hand, would be in conflict with the observed Higgs signal strength and electroweak precision observables. For g q and sin θ being sizeable, both the dark Higgs and the Z can mediate the interactions of DM with SM particles and all the constraints discussed in the two previous sections apply.  Figure 5. Constraints for the case of two mediators, with the (common) DM-mediator coupling determined by the requirement to reproduce the observed relic abundance. Plots in the same row correspond to constant θ, plots in the same column correspond to constant g q . In all panels, we have fixed m χ = 100 GeV. In the grey shaded regions (solid lines) at least one of the couplings leads to violation of perturbative unitarity. The yellow shaded regions (dotted) are excluded by dijet searches, the green shaded regions (short dashed) by monojet searches, the red shaded regions (long dashed) by direct detection, the purple shaded regions (double dash-dotted) by the observed Higgs signal strength and the bound on invisible Higgs decays. The dark blue regions (short dashdotted) are excluded by EWPT and the light blue regions (long dash-dotted) by dilepton searches, both for loop-induced kinetic mixing.
Figs. 5 and 6 show the resulting constraints in the m s -m Z plane for different combinations of g q and sin θ as well as for m χ = 100 GeV and m χ = 500 GeV, respectively. We make the following general observations: • As before, the weakest constraints are found if DM annihilations in the early universe receive a resonant enhancement, i.e. if either m χ ≈ m Z /2 or m χ ≈ m s /2.
• For m χ < m s , m Z the relic density is controlled by direct annihilation of DM into SM final states. Depending on the exact values of the couplings and masses, the  • For either m s < m χ or m Z < m χ , the relic density can be easily reproduced via annihilation into two dark terminators. However, if g q and sin θ are fixed to relatively large values, these regions are typically tightly constrained by direct detection and the Higgs signal strength.
• Moreover, for m s < m χ the requirement that y χ remains perturbative places an upper bound on m Z , while for m Z < m χ we obtain an upper bound on m s above which the Higgs potential becomes non-perturbative.
For m χ = 100 GeV we see that -for the coupling combinations considered in figure 5 -only small regions of parameter space close to the two resonances remain viable, whereas much larger regions are still allowed for m χ = 500 GeV. The reason is that in the latter case, it is possible for m s or m Z (or both) to be smaller than m χ without immediately being strongly constrained by direct detection or Higgs physics.  The fact that direct detection bounds are so strong in figures 5 and 6 is a direct consequence of the assumption that g q and sin θ are sizeable. For m Z , m s < m χ , however, direct annihilation into two terminators typically gives the dominant contribution to DM annihilation. Consequently, it should be possible to reproduce the required DM relic abundance for much smaller values of g q and sin θ, provided g χ is sufficiently large. Figure 7 summarises the constraints for g q = θ = 0.01. The four different panels are constructed in analogy to figures 2-6, i.e. they show the m χ -m Z plane, the m χ -m s plane and the m s -m Z plane for m χ = 100 GeV and m χ = 500 GeV, respectively. We find that the parameter region with m Z , m s > m χ is very tightly constrained, essentially because it is impossible to reproduce the relic abundance via annihilation into SM final states using perturbative couplings. For smaller masses of either m Z or m s , large allowed regions open up, which are very difficult to probe experimentally. If both m Z and m s are small, however, the indirect detection constraints discussed above become important, provided the DM mass is sufficiently small so that Fermi LAT is sensitive to the thermal cross section.
To conclude this section, we note that for m χ ∼ 30-50 GeV it may be possible within this framework to provide a viable explanation for the Galactic centre excess [85][86][87][88][89]. For example if 0.5 m χ < m Z < 0.9 m χ and m s + m Z > 2 m χ the s-wave cross section of χχ → Z Z is naturally close to the thermal value. Alternatively, for m s m Z and m s +m Z ≈ 2 m χ , it is possible that both χχ → Z s and χχ → ss contribute at comparable level to thermal freeze-out. The annihilation cross section for the first processes can then be a factor of a few below the thermal cross section, while the second process is velocity suppressed and therefore becomes negligibly small in the present Universe. This way, it may be possible to evade constraints from dwarf spheroidals while still obtaining an acceptable fit to the Galactic centre excess within astrophysical uncertainties.

Global scan of couplings
While the various panels shown in figures 5 to 7 are helpful in order to develop an intuitive understanding of how the different bounds depend on the couplings g q and sin θ, it is difficult to derive general conclusions based on these plots. For example, for m χ = 100 GeV and m Z = m s = 400 GeV, the couplings considered in figure 5 are excluded by LHC searches, while for the couplings considered in figure 7, these masses already violate the perturbativity constraints. This immediately leads to the question whether there might be an intermediate value for g q and θ for which LHC constraints can be avoided and the perturbativity requirement satisfied.
To answer this question, we have performed the following global analysis. For fixed masses we scan over g q and θ and determine for each choice of those parameters the dark sector coupling by the relic abundance. Each combination of masses then falls into one of three categories: 1. All combinations of g q and θ are excluded by at least one experimental constraint or perturbative unitarity.
2. There is at least one combination of g q and θ which is consistent with all experimental constraints and perturbative unitarity.
3. There is at least one combination of g q and θ, for which current experimental constraints do not apply and which therefore cannot presently be excluded.
The third situation occurs whenever a combination of couplings is most strongly constrained by LHC searches for dijet resonances. These searches typically assume Γ Z /m Z < 0.3, corresponding to g q 0.8, so that the bounds do not apply for larger widths/couplings. While it is to be expected that larger quark couplings lead to stronger rather than weaker constraints, we prefer not to perform an extrapolation and instead treat this special case separately.
We have performed such scans for a large number of different values of m Z and m s . The results of these scans for m χ = 100 GeV are summarised in the top-left panel of figure 8. The three categories discussed above are shown in red (excluded), white (allowed) and orange (inconclusive bound for broad mediator width). To illustrate these different situations, we consider in the remaining panels of figure 8 five different benchmark cases, each of which corresponds to a set of m χ , m Z and m s (indicated by black dots in the top-left panel). For each benchmark case we show the different constraints as a function of the couplings g q and θ (with g χ fixed by the relic density requirement) using the same colour coding as in figures 5-7. We make the following observations: • General observations: As expected, the constraints from dijets, dileptons and EWPT exclude large values of g q , while the constraints from the Higgs signal strength exclude large values of θ. Monojet searches are sensitive to an intermediate range of g q , as very large values of g q suppress the invisible branching ratio of the mediator and hence the monojet signal, while very small values of g q suppress the total production cross section of the mediator below the observable level. The behavior of the direct detection bound, on the other hand, is less intuitive: if the relic density is dominantly set by annihilation into quarks via the Z , then for small values of θ both the annihilation and the direct detection rate scale proportionally to g 2 q g 2 χ . In other words, if g χ is fixed by the relic density requirement, direct detection will either exclude all values of g q or no values of g q . For larger values of θ the direct detection bound on Higgs exchange becomes relevant for small g q (corresponding to large g χ and hence large y χ ). Thus for m χ < m s , m Z , direct detection typically gives the strongest constraint on small g q and large θ.
• Benchmark point 1: The combination m χ = 100 GeV, m Z = 75 GeV and m s = 75 GeV is allowed due to the presence of two light terminators. Since DM annihilation proceeds dominantly via χχ → Z s, the two couplings g q and θ can be made almost arbitrarily small, thus evading all experimental constraints. We note that in principle it is sufficient to have one light terminator, with the second mediator being much heavier than the DM particle. However, perturbativity limits place an upper bound on the mass of the second mediator. For m χ = 100 GeV we find approximately m s 800 GeV and m Z 2 TeV. Conversely, if both m Z and m s are very small, the indirect detection constraints discussed above become relevant.
• Benchmark point 2: For m χ = 100 GeV, m Z = 1000 GeV and m s = 210 GeV, it is impossible to evade all constraints for θ ≈ 0. While small values of g q (and correspondingly large values of g χ ) would be excluded by perturbativity, larger values are constrained by direct detection as well as searches for dijet and dilepton resonances. However, for θ > 0 there is plenty of parameter space to evade all experimental constraints due to the resonant enhancement of the annihilation processes χχ → s → SM.
• Benchmark point 3: For m χ = 100 GeV, m Z = 300 GeV, m s = 300 GeV, one can observe the strong complementarity between different experimental probes. Direct  . Global scan for m χ = 100 GeV. In the first panel (top-left) a global scan has been performed to determine those values of the two SM-mediator couplings that give the weakest constraints (see text for details). The red shaded region is excluded for all possible combinations, while in the white region all constraints can be evaded. In the orange shaded region it is not possible to exclude large values of g q , corresponding to Γ Z /m Z > 0.3. The remaining five panels show the various constraints as a function of the couplings between the two mediators and SM states for different benchmark cases (indicated by black dots in the top-left panel) with the common DM-mediator coupling fixed by the relic density requirement. The colour coding is the same as in figures 5-7. detection experiments and bounds on the Higgs signal strength exclude large values of θ, while searches for monojets and for dilepton resonances exclude wide ranges of g q . Nevertheless, all constraints can be evaded for g q ≈ 0.04, so that this parameter point is allowed.
• Benchmark point 4: For m χ = 100 GeV, m Z = 350 GeV, m s = 800 GeV all constraints become somewhat stronger compared to benchmark 3. In particular, the perturbativity constraint now excludes g q < 0.1 (for θ 1), closing the previously unconstrained gap around g q ≈ 0.04. Large values of g q are excluded by searches for dilepton resonances. However, these searches cannot be applied to g q 0.85 due to the broadening of the resonance. We can therefore only reliably exclude smaller couplings of g q .
• Benchmark point 5: The case m χ = 100 GeV, m Z = 500 GeV, m s = 500 GeV is again more strongly constrained than benchmark 4. Crucially, direct detection experiments are now also sensitive to large values of g q , thus covering the parameter region where dilepton searches lose sensitivity. The combination of all these constraints therefore allows to exclude this combination of masses for all coupling values that reproduce the relic abundance.
The fact that in many cases a combination of different constraints is necessary to conclusively rule out a given combination of masses (see e.g. benchmark point 5) illustrates the necessity for comprehensively scanning over the two couplings. We observe that indeed large parts of the parameter space are excluded for all combinations of couplings that reproduce the relic abundance. For m Z > 2 m χ , LHC constraints are typically very important and push the model to very broad widths. In this parameter region even stronger constraints can be expected in the near future due to improved direct detection and collider constraints. Benchmark point 3 for example may soon be excluded for all combinations of g q and sin θ. Indeed, for a traditional dark mediator with mass larger than twice the DM mass, the simple thermal freeze-out picture appears under considerable pressure. Nevertheless, sizeable allowed regions remain for m Z , m s < m χ , unless both masses are small enough for indirect detection constraints to become important.
In figure 9 we present the results of the global scans for different values of m χ . We make qualitatively the same observations as in the case m χ = 100 GeV, noting that small DM masses are already very tightly constrained by the data. We clearly observe that consistent points can only be obtained at the s-channel resonances or with at least one terminator. For m χ < 100 GeV and m s + m Z < 2 m χ indirect detection constraints from χχ → Z s give a constraint that cannot be avoided by varying g q and sin θ. For m χ = 30 GeV there is an additional constraint from χχ → Z Z , which is relevant for m s > m χ and 0.5 m χ m Z 0.9 m χ . For larger DM masses the inconclusive large-mediator-width regions become more important, but heavy mediators are still tightly constrained. On the other hand, the size of the parameter region with m Z , m s < m χ grows and indirect detection constraints are absent. The configuration with one or two dark terminators thus remains viable up to   1. We consider a Majorana DM particle which is a singlet under the SM gauge group but charged under a new U (1) gauge group.

The dark Higgs mediates DM-SM interactions via mixing with the SM Higgs doublet.
To obtain a Z mediated interaction we assume that quarks are charged under the U (1) as well. 4. We assume that left-handed and right-handed quarks have the same U (1) charge (vector-like coupling to Z ), that the quark-coupling is flavour universal, and that the charge of leptons is zero. Hence the U (1) corresponds to gauged baryon number.
5. We take into account only loop-induced kinetic mixing between the SM hypercharge gauge boson and the Z , which implies that the kinetic mixing parameter is not an independent parameter of the model. 6. The DM-related phenomenology is completely described by the model. In particular, we require that the total observed abundance of DM is obtained by thermal freezeout within this model, where all available annihilation channels are self-consistently determined within the adopted choices of interactions.
We denote the model by 2MDM. It is characterised by 6 parameters (3 masses and 3 couplings) as given in table 1, and it can be considered as a joint framework for simplified DM models with the simultaneous presence of a scalar and a vector s-channel mediator (see table 2). The 2MDM model is gauge invariant at the Lagrangian level and renormalisable, and we consider only regions of the parameter space where perturbative unitarity is guaranteed. The stability of the DM particle is a consequence of the U (1) gauge symmetry. In principle the model requires additional fields in order to cancel the gauge anomalies. However, the vectorial coupling of the Z to quarks ensures that the new fermions are not colour-charged and, therefore, the new particles are not expected to have a substantial impact on the phenomenology that we consider. Including additional states to cancel anomalies would lead us to recover the models discussed in the context of Baryonic DM [33,34].

Main results
We have analysed the 2MDM model under the assumptions specified above, taking into account a large variety of constraints: direct searches (LUX) and indirect searches (Fermi-LAT γ-ray observations of dwarf spheroidals), various observables at colliders (monojets, dijets, dileptons), electroweak precision tests, and Higgs observables at the LHC.
We find that generally the WIMP hypothesis is under severe pressure. Within the considered model there are basically only two possibilities to obtain the correct relic abundance without running into conflict with experimental constraints and/or perturbative unitarity: • either the DM and mediator masses are tuned close to an s-channel resonance, or • one or both mediators are lighter than the DM, such that DM annihilations into the dark sector control the relic abundance, while the coupling to the SM can be quite small. In this case, one of the mediators (or both) actually becomes a terminator.
The presence of a terminator makes it difficult to test this region of the parameter space. In our model, typically the annihilation cross section is p-wave and/or helicity suppressed, which makes indirect detection signals from the decay of the terminator into SM particles unobservable. However, there is one exception, namely if the annihilation channels χχ → Z s or χχ → Z h become kinematically allowed. Those processes are a generic feature of the two-mediator set-up, they appear at s-wave, and typically dominate the DM annihilation in the present Universe. Fermi-LAT observations of dwarf spheroidals constrain this configuration for DM masses m χ 100 GeV.
It is well known that direct detection strongly constrains a thermal WIMP with an s-channel scalar mediator, e.g., [40]. The vector mediator considered in our model has axial couplings to DM and vectorial couplings to quarks, which leads to a velocity suppressed scattering cross section with nuclei, and therefore bounds from direct detection are usually expected to be much weaker. However, we have shown that current limits from LUX are strong enough that despite the velocity suppressed interaction, direct detection provides a relevant constraint on the parameter space of the vector-mediator model, under the assumption of the correct thermal freeze-out abundance.

Outlook
To conclude, let us discuss possible future directions. First of all, significant progress is expected for the near future in our knowledge about the SM-like Higgs boson. More precise measurements of the total signal strength may further constrain or in fact provide evidence for mixing between the observed Higgs boson and additional Higgs bosons. Similarly, relevant constraints are expected from searches for invisibly decaying Higgs bosons and from searches for non-standard Higgs decays such as h → 4b or h → 4τ . In the presence of a light Z terminator, it is also conceivable in our model to have a significant branching ratio for h → 4j.
We furthermore believe that the results presented in the present paper motivate a new class of LHC searches, namely searches for a light dark Higgs. If the DM relic abundance is set via χχ → ss, the dark Higgs would by definition decay visibly (because m s < m χ ), with bb being the dominant decay mode in large regions of parameter space. The difference to conventional searches for additional light Higgs bosons (as conducted for example in the context of the NMSSM) is that the dark Higgs would generically be produced in association with large amounts of missing transverse energy, because the dominant production modes are dark Higgs Strahlung (pp → Z * → Z s → χχs) and final-state radiation of a dark Higgs (pp → Z → χχ → χχs). In contrast to existing searches for heavy quarks and missing energy, on the other hand, the fact that both b-quarks result from the decay of an on-shell resonance means that the distinctive invariant mass distribution can be used to distinguish the signal from potential backgrounds. 7 A detailed study of the LHC sensitivity for this signature will be left for future work.
Another potential constraint not taken into account in the present work arises from DM capture in the Sun and the resulting neutrino signal in IceCube (see e.g. [6,29]). However, since in the set-up we consider DM-nucleus scattering is dominantly spin-independent, we expect solar capture not to be competitive with direct detection experiments. Nevertheless, future progress in indirect detection experiments is crucial in order to further constrain the case of one (or several) dark terminators. In particular, CTA has the potential to substantially extend the sensitivity of indirect searches up to DM masses above the TeV range [90]. Furthermore, it will be interesting to see whether the Galactic centre excess can be attributed to unresolved point sources as recently suggested [91,92]. If the possibility of a DM interpretation persists, 2MDM can easily accommodate such a signal via cascade annihilations. In particular, the observed relic abundance can be obtained with comparable contributions from s-wave and p-wave annihilation, so that the present-day annihilation cross section would be a factor of a few below the thermal cross section.
To conclude, taking the 2MDM model as a benchmark scenario, we have seen that the WIMP hypothesis of a thermal relic DM abundance from weak-scale physics is under severe pressure. While "classic" WIMP scenarios with mediators comparable or heavier than the DM particle are largely excluded, the 2MDM model is flexible enough to save the WIMP due to resonances or the presence of a terminator. Although our study is performed within the specific 2MDM model, we expect our conclusions to apply to a larger class of models where the DM particle is a SM gauge singlet. Examples for alternatives not directly covered by our results are DM models where DM is charged under the SM gauge group [93,94] or where co-annihilations are relevant [95]. Generically we can conclude that if DM is a WIMP then the dark sector is most likely more complicated than just a weak-scale DM particle with effective interactions with the SM.
where L SM is the usual SM Lagrangian, and Here, S denotes the complex dark Higgs, H is the SM Higgs doublet, and the gauge kinetic terms are defined as where B µ is the SM U (1) Y gauge field. After spontaneously breaking the U (1) and the electroweak symmetries we go to the unitary gauge, where Here v = 246 GeV is the SM Higgs vev. Then, we obtain Using the phase freedom for S and χ we can choose both w and y χ real without loss of generality, which ensures a real mass for χ and a pure scalar coupling of s; since w is the only source of U (1) symmetry breaking we cannot obtain a pseudo-scalar coupling of s to χ [40]. The model has the following independent new parameters in the Lagrangian: since the dark Higgs can only couple to the DM (and eventually give it a mass) if q S = −2q χ , and the minimisation conditions of the scalar potential enforce the relation Additionally, note that the g gauge coupling only appears in connection with the charges of the fermion fields, which is expected since the normalisation of an Abelian force is not well defined. Thus, freedom in g can be hidden in the fermion charges, and g should not be taken as an independent parameter. We can define 8 and describe the mixing in the Higgs sector as where θ is required to lie between −π/4 and π/4, such that by definition H 1 is the mostly SM-like state. Correspondingly, we denote the mass of H 1 as m h and the mass of H 2 as m s . These are given by Note that these expressions allow for both m h < m s (if v 2 λ h − w 2 λ s < 0) and m h > m s (if v 2 λ h − w 2 λ s > 0). Then, the free parameters of the model can be taken as g χ , g f , m χ , m Z , m s , θ, , (A. 21) and the relevant parts of the Lagrangian are given in terms of these parameters by where the parameters from the Higgs potential are given in terms of the free parameters in eq. (A.21) and m h and v as To derive these equations, we have used eq. (A.12) with w replaced appropriately. If we allow for a U (1) charge for the SM Higgs, the term appears in the Lagrangian of the model. Once the SM Higgs develops a vev, this interaction leads to mass mixing between the SM Z and the Z gauge bosons: where s W = sin θ W and c W = cos θ W with the weak mixing angle θ W . Furthermore, the invariance of the SM Yukawa terms under SU (2) L ×U (1) Y ×U (1) implies axial interactions of the Z with the SM fermions, It should be noted that the interaction strength g A f is not a free parameter here and can be related to the Higgs charge by g A f = 1 2 g q H . This implies in particular that g A f has to be independent of flavour and equal for quarks and leptons. Therefore, searches for dilepton resonances strongly constrain this option (see appendix C).

B Gauge boson mixing
Our treatment of Z-Z mixing follows closely the discussion outlined in [16] and, therefore, we limit ourselves to a schematic description in the following. In its most general form the gauge part of the Lagrangian presented in the previous section can be written as where we use hats to denote non-normalised fields and defineẐ ≡ĉ WŴ 3 −ŝ WB and F µν ≡ ∂ µX ν − ∂ νX µ . First, the kinetic term can be brought to the canonical form by the transformation which generates an additional contribution to the off-diagonal mass term. In a second step the mass terms are diagonalised by the rotation with the rotation angle ξ given by As a consequence of the mixing, the physical mixing angle s W differs from the fundamental mixing angleŝ W . This difference, however, is only relevant at higher orders in the mixing parameters and will therefore be neglected in the following. Both kinetic and mass mixing can have a profound impact on the phenomenology of the model. The field redefinition induces corrections to the properties of the SM gauge bosons. These deviations from the SM expectation can be parametrized by the electroweak precision observables S and T and, expanding in ξ and , we find to leading order [96] where α = e 2 /4π. For small values of m Z , it is more appropriate to use the ρ parameter, ρ ≡ m 2 W /(m 2 Z c 2 W ), which in the presence of mixing will deviate from unity: (B.6) Furthermore, due to the mixing with the Z, the Z acquires couplings to all SM particles. From a phenomenological point of view the most interesting of these are the couplings to leptons which can lead to observable signals in searches for dilepton resonances at colliders. If the tree level coupling to leptons is zero, the induced coupling between the new gauge boson and leptons is given by where g and g are the gauge couplings of SU (2) L and U (1) Y .

C Relic density from mixing
In this appendix we consider two additional ways in which the Z can couple to SM states: tree-level kinetic mixing and mass mixing. The latter case is most easily realised if the SM Higgs carries a charge under the new U (1) . This also implies that the Z has couplings to SM fermions and, as a result, the axial coupling is directly related to the Higgs charge by  Figure 10. Bounds on scenarios with kinetic mixing (left) and axial couplings (right). We fix the dark-sector coupling g χ = 1 and vary (left) or g A q (right) in order to obtain the correct relic density. The grey shaded region (solid line) corresponds to > 1 (g A q > 1). The light blue (long dash-dotted) and dark blue (short dash-dotted) shaded regions are excluded by searches for dilepton resonances and EWPT, respectively. The green shaded region (dashed) corresponds to m χ > m Z . In this parameter region, the relic abundance depends only on g χ and not on the mixing parameter, which can therefore be arbitrarily small. g A f = −2g q H due to gauge invariance. The Z vector couplings to fermions, on the other hand, are not uniquely fixed in this model. In the following we will set them to the minimal value in agreement with gauge invariance, g V f = g A f , see [16] for a detailed discussion of this choice.
In principle, either of the two mixing scenarios could conceivably generate large enough couplings of the Z to the SM to provide the necessary annihilation channels in order to reproduce the observed relic abundance. However, gauge boson mixing is strongly constrained by electroweak precision observables. The contribution of mixing to EWPT is given by eq. (B.5). Moreover, both mixing scenarios imply that the Z couples to leptons either through mixing or directly and, therefore, it can be constrained by searches for dilepton resonances. In the following we use the same experimental constraints as described in section 3.2.
The bounds discussed above can be compared to the coupling strength required by the relic density. The results of this analysis are shown in figure 10. For both panels in the figure we have conservatively set g χ = 1 in order to reduce the necessary interaction strength with the SM and fix m s = 3 m χ to remove the contribution of the dark Higgs to the relic density. The two panels correspond to varying kinetic mixing and varying Higgs charge, respectively, such that the correct relic abundance is obtained. Since we fix the dark sector coupling to a relatively large value, we generally predict DM underproduction in the parameter region where m χ > m Z and the process χχ → Z Z is allowed. We shade this region in green, since mixing plays a sub-dominant role here.
As can be seen in the left panel of Fig. 10 the values of which are necessary to account for Ω DM are excluded by EWPT and dilepton searches in most of the remaining parameter space. Nevertheless, kinetic mixing could explain thermal DM if the annihilation cross section receives a resonant enhancement from the Z or Z . Similarly, we find that processes induced by mass mixing (i.e. axial couplings) cannot account for the observed relic density once EWPT and dilepton searches are taken into account unless there is a resonant enhancement of the annihilation rate.

D.1 Standard model final states
Annihilation into quarks via the Z is given by whereas annihilation into quarks via the dark Higgs is given by σv χ (χχ → s → qq) = 3 g 2 χ sin 2 θ cos 2 θ 2π The same expression (with m s replaced by m h ) is obtained for annihilation via the SM Higgs. For m s comparable to m h , the interference between the two diagrams becomes important. Due to the different sign in the mixing matrix, the interference is destructive. The total annihilation cross section is therefore proportional to (m 2 The SM Higgs (as well as the dark Higgs) also mediates additional annihilation processes, most notably annihilation into SM gauge bosons and the SM Higgs: σv χ (χχ → s → W + W − ) = g 2 χ sin 2 θ cos 2 θ 4π These processes are enhanced relative to annihilation into quarks by a factor m 2 χ /m 2 q and therefore become dominant for large DM mass.

D.2 Dark terminators
For the process χχ → ss, we are only interested in the limit m s m χ . We then obtain σv χ (χχ → ss) where for the term proportional to v 2 χ we only retain the leading p-wave contribution; additional terms are relevant if m Z and m χ are comparable. As can be seen the p-wave contribution is enhance by a factor m 4 χ /m 4 Z relative to the s-wave part. The reason for this behaviour is that the dominant contribution to the p-wave channel results from annihilation into pairs of longitudinal Z bosons. These are the Goldstone modes and couple in the same way as the dark Higgs, so that the annihilation cross section is proportional to y 4 χ or, equivalently, to g 4 χ m 4 χ /m 4 Z . Finally the annihilation cross section for χχ → Z s is given by (D. 6) In the limit that both the Z and the dark Higgs are light compared to the DM particle, this expression reduces to eq. (5.1). When kinematically accessible this processes typically dominates over χχ → ss and χχ → Z Z due to their velocity suppression and smaller prefactors.