Multilepton Signatures from Dark Matter at the LHC

Leptonic signatures of Dark Matter (DM) are one of the cleanest ways to discover such a secluded form of matter at high energy colliders. We explore the full parameter space relevant to multi-lepton (2- and 3-lepton) signatures at the Large Hadron Collider (LHC) from representative minimal consistent models with scalar and fermion DM. In our analysis, we suggest a new parametrisation of the model parameter spaces in terms of the DM mass and mass differences between DM and its multiplet partners. This parametrisation allows us to explore properties of DM models in their whole parameter space. This approach is generic and quite model-independent since the mass differences are related to the couplings of the DM to the Standard Model (SM) sector. We establish the most up-to-date LHC limits on the inert 2-Higgs Doublet Model (i2HDM) and Minimal Fermion DM (MFDM) model parameter spaces, by using the complementary information stemming from 2- and 3-lepton signatures. We provide a map of LHC efficiencies and cross-section limits for such 2- and 3-lepton signatures allowing one to easily make model-independent reinterpretations of LHC results for analogous classes of models. We also present combined constraints from the LHC, DM relic density and direct search experiments indicating the current status of the i2HDM and MFDM model.


I. INTRODUCTION
The nature of Dark Matter (DM) is one of the greatest puzzles of modern particle physics and cosmology. While overwhelming observational evidences from galactic to cosmological scales point to the existence of DM [1][2][3], after decades of experimental effort, only its gravitational interaction has been confirmed. At the moment, we do not have information about DM properties, such as its mass, spin, interactions (other than gravitational), symmetry responsible for its stability, number of states associated to it and possible particles that would mediate the interactions between DM and the Standard Model (SM) particles.
If DM is light enough and interacts with SM particles directly, or via some mediators with a strength beyond the gravitational one, its elusive nature can be detected or constrained in direct production at colliders, complementing direct and indirect DM searches in non-collider experiments. Therefore, the search for DM is one of the top priorities of the Large Hadron Collider (LHC) [4] programme and that of future collider experiments.
The most general DM signature at colliders is the mono-X one, where X stands for a SM object, such as jet, Higgs, Z, W , photon, top-quark, etc. that recoils against the missing energy from the DM pair. This signature has limitations, though, especially if there are no light mediators decaying into a DM pair which could enhance the corresponding signal to compete with the large SM background from Z → νν. If such mediators are absent or very heavy even the High-Luminosity LHC (HL-LHC) [5] could probe the DM mass only up to about 250 GeV, as shown, e.g., for the case of higgsino DM in the Minimal Supersymmetric Standard Model (MSSM) [6][7][8][9]. This limitation motivates us to look for signatures beyond the mono-X one.
This study is devoted to a generic class of DM models, where the DM is a part of an Electro-Weak (EW) multiplet and the mass splitting between DM (D 1 ) and its charged partner(s) (D + , ...) or next-to-lightest neutral partners (D 2 , ...) is large enough so as to give rise to multi-lepton signatures at the LHC from processes such as pp → D 1 D ± , pp → D + D − , pp → D 1 D 2 , pp → D ± D 2 and pp → D 2 D 2 , which are then followed by the D ± and D 2 decays to gauge or Higgs bosons which in turn decay into leptons.
These signatures with visible 1-4 lepton(s) and Missing Transverse Energy (MET) can be observed for ∆m + = m D ± − m D 1 and/or ∆m 0 = m D 2 − m D ± above a few GeV. The case when ∆m + is of the order of the pion mass -when D + and D 1 are degenerate at tree-level and their mass split is generated radiatively, due to the quantum corrections involving a loop of the photon and the charged DM multiplet partner(s) -leads to a compelling complementary signature with disappearing tracks which was recently studied in detail in [10].
The 2-and 3-lepton signatures from DM models have been previously explored in the context of the inert 2-Higgs-Doublet Model (i2HDM) [11], Supersymmetric scenarios [12,13] (including those with sneutrino DM [14]) or models with vector-like leptonic DM [15]. These earlier studies have mainly explored the one or two-dimensional parameter space of specific DM models, or just selected benchmarks points from their multi-dimensional parameter space. Our study presents the following new results on multi-lepton signatures from DM at the LHC: • we explore the full three-dimensional parameter space relevant to the LHC for two representative minimal consistent DM models with DM of spin-0 and -1/2, respectively: the i2HDM and Minimal Fermion DM (MFDM) model, respectively, and present the sensitivity of recent LHC data to these; • we suggest a new parametrisation for both models which allows one to better understand and interpret their properties and visualise the no-loose theorem in their full parameter space; • we find that the 3-lepton signature becomes relevant for large D + − D mass gaps and adds important and complementary sensitivity to the LHC searches; • we implement and validate the 8 TeV ATLAS [16] multi-lepton analysis of Ref. [17] in Check-MATE [18] relevant to our study and made available those analyses to the community; • we create a map of the LHC efficiencies and the cross-section limits for 2-and 3-lepton signatures in the simple parameter space expressed through the DM mass and the D + − D and D 2 − D + mass differences which can be used for generic reinterpretation of spin-0 and spin-1/2 DM models.
The rest of the paper is organized as follows: in section II we present the DM models under study, in section III we discuss details of the signal processes leading to 2-and 3-lepton signatures, in section IV we present details of the analysis, including the implementation into the CheckMATE package, and the new LHC sensitivity to the parameter space of the two models under study while in section V we draw our conclusions. Furthermore, details of our study are given in Appendices A-E, including the respective plots on the exclusion of the parameter space (Appendix D) and tables with the efficiencies and cross-section limits expressed in terms of universal parameters such as the DM mass and the D + − D and D 2 − D + mass differences which can be applied for generic reinterpretation of spin-0 and spin-1/2 DM models (Appendices F-G).

II. MODELS
In the following subsections, II A and II B, we present the models under study in this work. Both of these are phenomenologically well-motivated and can be mapped onto well-known Supersymmetric scenarios. They are minimal in nature, due to the small number of parameters that impact the phenomenology, in terms of both collider and DM considerations.
For the purposes of this paper, where we are interested primarily in the collider phenomenology of the model, we do not initially impose strict requirements on relic density and (in)direct detection of DM. Rather, we first set out the collider constraints independently of DM considerations. Only eventually we show how the latter impact the parameter space, in order to illustrate that these are complementary to collider constraints.

A. MFDM Model
The Lagrangian for the MFDM model, which augments the SM by an EW fermion (ψ) doublet of χ 0 1 , χ 0 2 and χ + as well as a Majorana singlet fermion χ 0 s , is [19] L MFDM = L SM + ψ(i / D − m ψ )ψ + 1 2 The fermion SU (2) doublet is while the Higgs doublet Φ, after EW Symmetry Breaking (EWSB), is where h is the SM-like Higgs boson with a mass of 125 GeV, while v = 246 GeV is Vacuum Expectation Value (VEV) of the SM Higgs potential. We impose a Z 2 odd parity on the EW fermion doublet as well as the Majorana singlet fermion, such that the lightest particle after diagonalisation of the fermion mass matrix is stable and a DM candidate. Note that, in the model as written above, we assume that χ 0 1 , χ 0 2 are Majorana fermions. The DM relic density requirements are satisfied by the usual freeze-out mechanism by scatterings mediated by the Z or Higgs boson. The minimality of this model is manifest in the fact there are only 3 new parameters, m ψ , Y DM and m S .
Note that χ 0 1 and χ 0 s mix via Yukawa couplings while χ 0 2 and χ + are mass degenerate. The physical masses are then obtained by diagonalisation of the mass matrix from the (χ 0 s , χ 0 1 , χ 0 2 ) interaction basis to the (D, D 2 , D ) mass basis as The mass matrix is thus diagonalised by one rotation angle θ of the upper 2 × 2 block, while we identify m ψ = m D + . The rotation angle θ that diagonalises the mass matrix can be expressed in terms of two useful relations, (6) and the connection between gauge eigenstates χ 0 s and χ 0 1 and mass eigenstates D and D 2 is given by χ 0 s = D sin θ + D 2 cos θ, χ 0 1 = D cos θ − D 2 sin θ.
The connection between physical masses and mass parameters in the Lagrangian is given by The mass splitting between m D 1 and m D 3 is related to the Yukawa coupling Y DM via where The requirement that Y DM must be real leads to the following mass ordering: The parameter space of the model can be expressed in terms of three physical masses: However, since the Yukawa coupling controls the mass splitting as well as DM mixing and therefore the strength of the interactions between DM and SM, it is also convenient to trade parameters and parametrise the model space using the following set which consists of the DM mass and two mass splittings between D + and DM as well as between D 2 and D + . As we will see, the latter parametrisation allows to better visualise and interpret collider search results for the whole parameter space. The model is also subject to various theoretical constraints, like perturbativity and radiative stability, which are satisfied for the couplings and masses adopted in this paper. As far as DM constraints are concerned, we first assess direct detection constraints. Since there is no tree level Z boson interaction with the DM candidate, the spin-independent direct detection constraint only arises from DM-nucleon scattering controlled by the Higgs coupling to DM. By controlling the mass splitting, this coupling can be made small, leaving a viable parameter space. In order to obtain the correct (not over-abundant) relic abundance, one needs an efficient annihilation mechanism via the Higgs funnel m D m h/2 . These considerations were previously taken into account in a previous paper [19].

B. i2HDM
This scenario [20][21][22] introduces a second scalar doublet φ 2 to the SM Higgs sector, with the same quantum numbers as the SM Higgs doublet φ 1 . However, a discrete Z 2 odd symmetry is imposed on the second doublet such that there are no direct coupling to fermions, rendering it 'inert'.
In the unitary gauge, the two scalar doublets can be written as where D 1 , D 2 and D + are the lightest, next-to-lightest and charged inert scalars, respectively. Within this parametrisation, the Z 2 symmetry is preserved due to the fact that v 2 = 0 and the absence of Yukawa couplings of the φ 2 to fermions. This 'inert minimum', for < φ 0 i >= v i √ 2 , corresponds to v 1 = v and v 2 = 0. This symmetry provides the stability of the lightest inert boson against decaying into SM particles, which makes D 1 a DM candidate.
The scalar Lagrangian is then given by where the most general scalar potential is written as and contains all of the scalar interactions allowed by the Z 2 symmetry. Here, all mass and coupling parameters are defined to be positive. Under the condition [23], where m D 1 is the mass of D 1 and λ 345 = λ 3 + λ 4 + λ 5 ,

6
The physical masses are given by (λ 5 can always be chosen positive due to the φ 2 → iφ 2 , λ 5 → −λ 5 symmetry) which, together with λ 345 and λ 2 , parametrise completely the DM sector of the model, described by five parameters: These can be further reduced for this analysis, by noting that λ 2 does not affect the LHC phenomenology under study at all since it controls DM the self-interaction only. Furthermore, we also exclude λ 345 from our analysis for the following reason. The relevant coupling for the gg → H * → D + D − and gg → H * → D 2 D 2 processes are λ 3 and λ 3 + λ 4 − λ 5 , respectively, which are limited by EW precision tests plus perturbativity [23,24] and related to the mass splitting between D + and D 2 . We have checked that, even if we maximise the cross section of either of the two processes by increasing the respective coupling up to the maximal allowed value, the respective cross section 1 will still be below the signal production cross section via the weak coupling. Moreover, the value of the λ 345 coupling is also strongly limited by DM direct detection constraints as well as by limits on the invisible Higgs boson decay branching ratio in the case when m H > 2m D 1 (see, e.g., [23,25]). Therefore, in order to establish a conservative and generic limit on the i2HDM parameter space, we exclude λ 345 from the current study and set its value to zero. Therefore, the parameter space we focus on is three-parametric: which can easily be visualised in 2D planes, while fixing one parameter. In analogy with the MFDM model we also use the parametrisation, where which allows an even better visualisation and interpretation of the model parameter space. We note one important caveat here. Note that we deliberately choose the mass hierarchy m D 2 > m D ± > m D 1 , such that D 1 is the lightest state. This choice is motivated by phenomenological reasons in order to bring out the intricacies of the collider constraints set by 2-and 3-lepton searches. The other hierarchy m D ± > m D 2 > m D 1 is entirely feasible and will lead to constraints comparable to the ones described in this paper. We leave this possibility for a future work.

III. SIGNAL PROCESSES
In this section, we study production and decay processes inclusively, the combination of which gives an idea about the LHC event rate of the signatures under study. Specifically, we study lepton signatures from DM multiplet partners generically appearing when the mass splitting between DM and those partners is large enough to give rise to sufficiently energetic leptons, which can originate from the following processes: where D n = D 2 or D n = D , D 2 depending on the scenario, i2HDM or MFDM model, respectively. Since D ± and D n decay via Z or W ± either hadronically or leptonically in association with the DM candidate, they can provide signatures with several charged leptons plus MET, which are the subject of this study. The common diagrams for both models under study for the signal production processes that provide the multi-lepton signatures we study here are presented in Fig. 1. In addition to these common Feynman diagrams, the ZD 2 D 1 and ZZD 1 D 1 vertices, specific only to the i2HDM, provides additional Feynman diagrams and corresponding new kinematic topologies, shown in Fig. 2 providing the required multi-lepton signatures through the topologies shown in Fig. 2(c), (d) and (e), respectively. The cross sections for the above production processes presented here are calculated using CalcHEP [26], a package designed for the evaluation of the tree-level processes and their respective For the latter, we use here the combination PYTHIA [29] and DELPHES [30], respectively. The relevant cross sections are presented in Fig. 3 as a function of the ∆m + parameter for a 100 GeV DM mass and two values of ∆m 0 =1 and 100 GeV (left and right panels, respectively) for the i2HDM (top panels) and MFDM model (bottom panels). Here, one can see that for ∆m 0 =1, D 1 D 2 production has the highest rate for the i2HDM while, in the case of the MDFM model, the analogous D 1 D channel is highly suppressed. This can be explained by the fact that the ZD 1 D 2 coupling controlling this process in the i2HDM model is just a pure (weak) gauge coupling while in the MFDM model the ZD 1 D coupling is the product of a (weak) gauge coupling and the cosine of the χ 0 1 − χ 0 s mixing angle, which is suppressed when ∆m + ∆m 0 = 1, as shown by the blue line in Fig. 4 In contrast, for ∆m 0 =1 GeV, the production cross sections for D + D − , D 2 D ± and D D ± (represented by red, green and brown lines, respectively) are close to each other in each model. This can be explained by the fact that the D + and D 2 masses are about the same as well as the couplings controlling these processes (which are purely (weak) gauge ones). Furthermore, one should note that the ZZD 1 D 1 coupling (unique to the i2HDM) contributes to ZD 1 D 1 production, the 2 → 3 process with a subdominant cross section in comparison with 2 → 2 production. However, when the D 2 → ZD 1 decay is open, the ZD 1 D 1 process will include the corresponding resonant 2 → 2 production and decay. Therefore, to avoid double counting, we do not present the cross section for this 2 → 3 process in Fig. 3(a).
For the MFDM model, the additional production process D D 2 (pink line), characteristic of this scenario, has a similar cross section to D + D − and D D ± for the very same reason. One should also note that the cross sections for scalar DM production are smaller than those for fermion DM production by the spin factor   Let us now consider the ∆m 0 =100 GeV case presented in the right panels of Fig. 3. In the i2HDM, the ∆m 0 = 1 → 100 GeV change (which increases the D 2 mass) equally suppresses D 1 D 2 and D 2 D ± production by about a factor of 4. For the MFDM model, the ∆m 0 = 1 → 100 GeV change affects both the D 2 mass (see Eq. (10)) and the χ 0 1 − χ 0 s mixing angle θ (the ZD 2 D and W + D 2 D − couplings are both proportional to sin θ), leading to a suppression of D D 2 and D 2 D ± production from both causes (see pink and green lines, respectively, in Fig. 3(d)). In contrast, D 1 D production is enhanced, since it is proportional to cos 2 θ, which increases with ∆m 0 (see Fig. 4(b)). Finally, the D + D − processes are not affected by a ∆m 0 variation in both models. One should also note that in the MFDM model D D ± is also not affected by a ∆m 0 variation since such a variation does change neither the D mass nor the W + D D − coupling (which is purely weak and is not affected by the χ 0 1 − χ 0 s mixing). In order to get an idea of the rate at which leptons are produced in the final state we also need to discuss Branching Ratios (BRs) of decay chains leading to DM particles and connect these with their pair production cross sections. In Fig. 5, we present the decay patterns common to both   the i2HDM and MFDM model. In Fig. 6, we then present the unique decays of each model: the D 2 → ZD 1 decay permitted in the i2HDM while the D 2 → HD 1 and D decays are so in the MFDM model. The BRs for the decays leading to our DM signatures (see Figs. 5 and 6) are presented in Fig. 7 as a function of the ∆m + parameter, again, for a 100 GeV DM mass and the two values of ∆m 0 =1 and 100 GeV. The D ± decay to W ± dominates in all cases. This is because, for the i2HDM, only m D 2 > m D ± is considered and so D ± can only decay to W ± D 1 while, for the MFDM model, D and D ± are mass degenerate and their decay to one another is not permitted. For ∆m 0 =1, in the i2HDM ( Fig. 7(a)), the D 2 decay to leptonically decaying Z and D 1 is favoured over the D 2 decay to leptonically decaying W + and D − which is suppressed due to a small D 2 − D ± mass splitting. This BR falls rapidly with increasing ∆m + , due to D + in the final state increasing in mass, while the mass of the D 1 in the D 2 → ZD 1 final state is fixed.
The behaviour of the D 2 BRs changes when considering ∆m 0 = 1 → 100 GeV ( Fig. 7(b)) since      the D 2 − D ± mass gap is increased, initially suppressing the D 2 decay to ZD 1 . However, the trend over increasing ∆m + remains the same and the two D 2 BRs cross for ∆m + around 60 GeV. Meanwhile, for the MFDM model, the D → W + D − decay is not permitted and so the BR of the D decay to a leptonically decaying Z and D 1 remains unchanged between Fig. 7(c) and (d).
The MFDM model decays from D 2 include a unique Higgs interaction, D 2 → D 1 H ( * ) , where the D 2 D 1 H coupling is a product of the Yukawa Y DM (which is proportional to √ ∆m 0 ∆m + ) and cos(2θ) of the χ 0 1 − χ 0 s mixing angle θ. This interplay is seen in Fig. 4(a) (blue line), where Y DM suppresses this coupling at small ∆m 0 and ∆m + , but the latter rises with increasing ∆m + . For ∆m 0 =1 the shape of the D 2 → D 1 H * (brown) line in Fig. 7(c) follows the trend just described, then it levels out as √ ∆m + . As the Higgs boson becomes on-shell with increasing ∆m + , a jump in the BR of leptonic decays via D 2 → D 1 H is observed in Fig. 7(c), due to the additional contribution H → W + W − . Meanwhile, the D 2 decays to W + D − and ZD are both controlled by the sine of the mixing angle, i.e., the amount by which the D 2 is χ 0 1 . Therefore, their shapes are similar to one another but the blue and black lines fall as the D 2 → D 1 H brown line increases with ∆m + (note that the D 2 → W + D − BR is a factor ∼ 8 larger than the D 2 → D Z one due to combinatorics and the ratio m Z /m W ). Now, considering ∆m 0 = 100 for the MFDM model, D 2 decays are strongly affected by any ∆m 0 variation. In Fig. 7(d), D 2 → ZD and D 2 → W + D − are now on-shell and these decays are then enhanced. When ∆m + = 30 GeV and m D 2 − m D 1 = 130 GeV, the D 2 → HD 1 decay (brown line) becomes on-shell and boosted, reducing the other D 2 decays proportionally. However, D 2 changes from mostly χ 0 s to mostly χ 0 1 as ∆m + reaches ∆m 0 from below, where the coupling reaches zero at ∆m + = ∆m 0 = 100 GeV (see Fig. 4(a) red line). Beyond this value of ∆m + , the coupling then increases with √ ∆m + as ∆m + > ∆m 0 which in turn reduces the D 2 → ZD and D 2 → W + D − leptonic BRs above ∆m + > 100 GeV as seen in Fig. 7(d).
The discussed combinations of production and decay rates provide the expected rates for the 2and 3-lepton signatures which we study below at the (fast) detector level and compare to published LHC data.

IV. RESULTS
In order to understand the constraints on the models we study in this paper, we reinterpret existing multi-lepton searches at the LHC. We first provide a brief summary of the models, then describe briefly the reinterpretation tools for this work, finally followed by an assessment of the impact of these searches on our two benchmark models.

A. LHC Searches and Tools
We begin by identifying the LHC searches that can potentially be useful to constrain the latter. For the i2HDM, constraints were derived previously in [11] by reinterpreting 8 TeV 2-and 3-lepton searches using MadAnalysis5 [31]. In preparation for this publication, the previous work is verified using the public recast tool CheckMATE [18] (see Appendix A). The corresponding recast software is publicly available at [32,33]. In this paper, we extend this result to 13 TeV for both the i2HDM and MFDM model based on searches available in CheckMATE.
In Tab. I, we identify the 13 TeV searches that are relevant for the reinterpretations in this paper. The most stringent constraints for these models are expected to emerge from LHC Supersymmetry searches targeting EW gauginos. Since we are investigating models with a variety of mass gaps in this paper, we look for searches that constrain both small and large mass gaps. We therefore note the following.  • For small mass gaps the CMS [41] soft lepton searches cms_sus_16_025 [39] and cms_sus_16_048 [40] can potentially constrain the parameter space under study here. The search cms_sus_16_025 constrains χ + 1 χ 0 2 pair production followed by decay to leptons and missing energy via off-shell W ± and Z bosons for mass gaps 5 − 50 GeV. It also constrains direct stop pair production followed by a decay to leptons via off-shell W ± 's for mass gaps up to 70 GeV. The search cms_sus_16_048 targets the same EWino pair production as above and constrains mass gaps up to about 50 GeV at an integrated luminosity of 35.9 fb −1 . Both of the above searches require Opposite Sign (OS) ee/µµ/eµ with leading leptons p T < 20 GeV, E miss T < 200 GeV and at least one jet.
• For large mass gaps the ATLAS [16] search atlas_conf_2016_096 [36] as well as the CMS search cms_sus_16_039 [38] are the most constraining ones. These publications target EWino pair production with the same decay pattern as the soft lepton searches. The searches look for OS leptons and constrain mass gaps above 50 GeV.
The i2HDM and MFDM model input LHE files were produced with different mass parameters, as described in section III, with CalcHEP [26]. This is followed by showering and hadronisation using PYTHIA8 [29]. Jets, with final-state hadrons are constructed using FASTJET [42], while detector simulation is performed using DELPHES [30]. The entire process (barring parton level event generation) is performed within CheckMATE. The built-in AnalysisHandler processes the detector-level events with the user selected analyses. The signal size is determined based on the efficiency, acceptance, signal cross section and integrated luminosity [43] of the analysis.

B. Constraints on the i2HDM and MFDM Model
In this section, we present the results of our study for the i2HDM and MFDM model. For each of these scenarios, we constrain the parameter space in the (m D 1 , ∆m + ) plane using the searches quoted above to calculate the r-value: where S DM is the number of DM events expected to pass the signal selection and S 95 is the 95% Confidence Level (CL) upper limit on the number of selected events. Any point with r ≥ 1 is excluded by current LHC limits. In general we observe the following.
• Low Mass Gap: 2l + p miss T channel -For low mass gaps the soft lepton analyses cms_sus_16_025 and cms_sus_16_048 constrain a large part of the parameter space in the 14 2l + p miss T channel. For the MFDM model, they exclude phase space within 30 < m D 1 < 160 GeV and 10 < ∆m + < 50 GeV for ∆m 0 = 1, 10, 100 GeV. For the i2HDM, the two soft lepton analyses exclude phase space within 1 < m D 1 < 60 GeV and 1 < ∆m + < 60 GeV for ∆m 0 = 1 GeV as well as phase space within 1 < m D 1 < 40 GeV and 1 < ∆m + < 40 GeV for ∆m 0 = 10, 100 GeV.
• 3l + p miss They all require three leptons with at least two e or µ forming an opposite-sign same-flavour (OSSF) pair with m ll < 75 GeV and including a τ -lepton. In the MFDM model, these signal regions exclude a phase space within 30 < m D 1 < 110 GeV and 20 < ∆m + < 300 GeV for ∆m 0 = 1, 10, 100 GeV. For the i2HDM, substantial r-value contributions are provided by the cms_sus_16_039 analysis with signal regions for three leptons, SR_A30 (three leptons with one OSSF pair 75 < m ll < 105 GeV, 250 < p miss T < 400 GeV), SR_K02 (four leptons, including two τ -leptons, 50 < p miss T < 100 GeV), SR_F10 (three leptons including two τ -leptons) and SR_I03 (four or more leptons, including one τ ), within the phase space 1 < m D 1 < 20 GeV and 5 < ∆m + < 70 GeV for ∆m 0 = 1, 10 GeV. For ∆m 0 = 1 GeV, the dominant DM pair production process is D 2 D 1 (see Fig. 3(a)) and the dominant D 2 decay contribution is D 2 → ZD 1 (see Fig. 7(a)), providing two leptons. The 3-lepton final states can be provided by the second largest DM production process, D 2 D ± , with decays as D 2 → ZD 1 and D ± → W ± D 1 . This process and its decays can also fulfil the 2-lepton criteria if the W decays hadronically, or a lepton is misidentified.
The horizontal wedge of large r-value in Fig. 8(a) within ∆m + < 60 GeV, m D1 < 70 GeV is excluded by analyses cms_sus_16_025 [39] and cms_sus_16_048 [40], both with signal region SR1_weakino_1low_mll_2, requiring two leptons with m ll < 20 GeV and at least one jet. In the i2HDM, this signature would be mostly provided by the leptonic Z decay and the hadronic W decay in the D 2 D ± pair production. Most of this phase space is already excluded by LEP-II observations [44] (green line) and LEP-I limits from Z width measurements (grey line), excluding on-shell Z → D 2 D 1 decays. However, at ∆m + < 8 GeV, there is a small wedge of allowed phase space from soft leptons, when m D 1 > 50 GeV and ∆m + < 8 GeV which is not covered by the LEP-I or LEP-II limits. The second, broader region of large r-value around 40 < ∆m + < 100 GeV is excluded by analysis atlas_conf_2016_096 [36], signal region 2LASF. The dominant D 2 D 1 process with leptonic Z decay provides a signal that strongly contributes to the r-value in this region. Note  [44], [45]. This includes m D 1 + m D 2 < m Z (grey line) and 2m D + < m Z (blue), from LEP-I Z boson width measurements excluding on-shell Z → D 1 D 2 and Z → D + D − decays. The same applies to the W width measurement excluding on-shell W ± → D ± D 1 decays. The region m D 1 < 80 GeV and m D 2 < 100 GeV and m D 2 − m D 1 > 8 GeV (green lines) is excluded by LEP-II observations. Where these lines are absent, they are overlapped completely by the other LEP limits.
that a significant portion in Fig. 8(a) within the region 60 < ∆m + < 95 GeV and m D 1 < 50 GeV not excluded by LEP, is excluded by these LHC limits.
Similar excluded regions are found for ∆m 0 = 10 GeV in Fig. 8(b), but without the region of allowed phase space at small ∆m + , for which the LEP-II limit now overlaps with this region (the grey line corresponding to the LEP limit m D 1 + m D 2 < m Z in Fig. 8(a) would be overlapped completely in Fig. 8(b) by the green line, so is not plotted here). Again, we see a significant contribution from LHC limits in 60 < ∆m + < 95 GeV and m D 1 < 50 GeV where LEP does not already exclude.
As ∆m 0 is increased further to 100 GeV in Fig. 8(c), the D 2 D 1 and D 2 D ± cross sections are now suppressed at small ∆m + (see Fig. 3(b)) compared to the dominant production D + D − for ∆m + < 130 GeV, with decays D ± → W ± D 1 (see Fig. 7(b)). The M ll > 100 GeV cut applied in the atlas_conf_2016_096 [36] analysis now removes all events for larger ∆m + in Fig. 8(c), as heavier D 2 production cross section with harder lepton decays has decreased. In addition, the combination of LEP limits covers the totality of LHC limits for ∆m 0 = 100 GeV in Fig. 8(c).

Constraints on the MFDM Model Channels
Based on the re-parametrisation of Eq. (10), Fig. 9 (overlaid in Fig. 20 of Appendix D with cross-section yields) shows the r-value exclusion contours of the MFDM model 2-and 3-lepton signatures discussed in section III as a function of ∆m + and m D1 for various choices of ∆m 0 . The LEP-II limit in the MFDM model case corresponds to bounds on fermionic DM [45], which covers m D + < 100 GeV, a smaller ∆m + range than the LEP limits for the i2HDM. Since we allow DM-Higgs coupling, the Higgs-to-invisible limit [46] of ∼ 0.15 BR (magenta region) is also plotted.
Contributions to the 2-lepton r-value are provided by the dominant D + D − production (see Fig. 3(c)-(d)) with its leptonic decays D ± → W ± D 1 (see Fig. 7(c)-(d)), by D D ± and, specifically for ∆m 0 = 1 GeV, D 2 D ± productions, where the latter two require D and D 2 to decay leptonically. Cascades from D D 2 production can also contribute to 2-lepton final states for ∆m 0 = 1 GeV. However, D 2 D ± and D D 2 production become suppressed with ∆m 0 = 100 GeV (see Fig. 3(d)). Meanwhile, D D 1 production, which is suppressed for ∆m 0 = 1 GeV, becomes enhanced with ∆m 0 = 100 GeV as detailed in section III, contributing to 2-lepton final states. Contributions to the 3-lepton r-value are provided by the D D ± production with fully leptonic decays, although these are less likely to satisfy the signal regions that require τ -leptons than D D 2 and D 2 D ± production with D 2 → D 1 H * → D 1 τ τ decay. However, as stated, D 2 productions become suppressed with ∆m 0 = 100 GeV (see Fig. 3(d)).
For ∆m 0 = 1 GeV, the sharp excluded region of large r-value in Fig. 9(a) within ∆m + < 30 GeV, m D 1 < 150 GeV was seen similarly for the i2HDM case in Fig. 8. As with the i2HDM case, this signal of soft leptons for the MFDM model in Fig. 9(a) is also excluded by the 2-lepton analyses cms_sus_16_048 and cms_sus_16_025 both with signal regions stop_1low_pt, weakino_1low_mll. These require two leptons and at least one jet, a signal also provided by the D D ± pair production, dominant here in the MFDM model (see Fig. 3(c)).
As ∆m + increases, the 3-lepton r-value becomes larger than the 2-lepton-only r-value in the region 30 < ∆m + < 100 GeV and m D 1 < 80 GeV, excluded by the cms_sus_16_039 [38] analysis signal regions requiring three leptons including τ . This is where the most dominant decay of D 3 changes from 1-lepton final states to 2-lepton final states including τ -leptons through virtual Higgs decays, as detailed in section III and Fig. 7(c). This gives more 3-lepton final states, with τ -leptons, from D + D 2 production. The additional third excluded area of smaller r-value than the previous two areas in Fig. 9(a), but still within an excluded region 120 < ∆m + < 260 GeV, m D 1 < 100 GeV has contributions to the signal from on-shell Higgs decay, where D 2 changes from decaying through H * , into decaying to real Higgs which further decays to two W bosons or τ -leptons. This contributes more to the r-value as ∆m + > M H for real Higgs production as detailed in section III. This region is also excluded by analysis cms_sus_16_039 [38], with its 3-lepton (including one τ ) signal regions.
As ∆m 0 is further increased to 10 GeV and 100 GeV in Fig. 9(b) and (c) respectively, the same r-value exclusion patterns are consistently observed. This no-lose theorem appears with the 2-lepton r-value due to the following scenarios: for small ∆m 0 , the cross section is large for light D 2 production, but with suppressed coupling between D 1 − D (see Fig. 7(c)). Then, for large ∆m 0 , this coupling is increased and (see Fig. 7(d)), while the heavy D 2 leads to suppressed production cross section.
For the 3-lepton r-value, while D 2 production becomes suppressed for increasing ∆m 0 as shown in Fig. 3, its decays D 2 → W + D − and D 2 → ZD are enhanced, which can easily provide three leptons, including τ -leptons required by the relevant signal regions. In addition, D D ± production cross section is not affected by the ∆m 0 variation.
As ∆m 0 is increased, the Higgs-to-invisible limit [46] of ∼ 0.15 BR excludes larger regions in Fig. 9(c)-(f), where H → D 1 D 1 becomes the Higgs' dominant decay channel (since this coupling is proportional to Y DM ), until D 1 is too heavy to be produced on-shell. Other than this, the similarities in r-value plots for these three ∆m 0 scenarios means we only need to consider the 2D plane (∆m + , m D1 ). Fig. 10 visualises the relation of r-values and BR by superimposing the decay   The DM mass is fixed to m D 1 = 100 GeV, corresponding to a vertical slice of Fig. 9(a) and (e). BR from Fig. 7, rotated to match the contour plots in Fig. 9, with the 2-(dot-dashed line) and 3-(dashed line) lepton r-values from a vertical slice of Fig. 9 at m D 1 = 100 GeV as a function of ∆m + .
For ∆m 0 = 1 GeV (Fig. 10(a)), the dominant contribution to the r-value switches from 2lepton final states to 3-lepton final states around ∆m + > 45 GeV, due to the change in dominant branching of D 2 → D ∓ W * ± (→ lν) (blue line) to D 2 → D 1 H * (→ τ + τ − ) (brown line) as ∆m + (and m D2 ) increases. This corresponds to a phase space excluded specifically by the analysis cms_sus_16_039 [38], requiring three leptons with at least one τ -lepton, which can be provided in significant quantities by the H * decay. Since the dominant productions at ∆m 0 = 1 include D 2 D + and D 2 D (see Fig. 3(c)), the dominating contribution to the r-value changes from D 2 D + → lνD + lνD 1 + X (total of two charged leptons) to D 2 D + → llD 1 lνD 1 + X and D 2 D → τ + τ − D 1 llD 1 (total of at least three charged leptons, including τ -leptons). As the Higgs boson becomes on-shell at ∆m + > 130 GeV, the H → W + W − channel opens, in addition to the real H → τ + τ − , noticeably contributing further to the 3-lepton r-value line in Fig. 10(a).
For ∆m 0 = 100 GeV in Fig. 10(b), while D 2 production is suppressed (see Fig. 3(d)), the co-dominant pair production D D ± also provides three leptons via D → D 1 Z (green line) and D ± → D 1 W ± (pink line). In addition, although D 2 production being suppressed, the D 2 decays D 2 → D ± W ± (→ lν) and D 2 → D Z(→ ll) become strongly enhanced throughout (since W and Z are on-shell). Combined, they provide enough leptons, including τ -leptons, to compensate for the suppression of D 2 production and maintain a significant contribution to the 3-lepton r-value with increasing ∆m 0 . As ∆m + increases, an increase in the D 2 → D 1 H branching, which also decreases the D 2 → D ± W ± (→ lν) branching, contributes with τ -leptons to the 3-lepton r-value around 40 < ∆m + < 60 GeV. As detailed in section III, the D 2 D 1 H coupling falls to zero again as ∆m + reaches ∆m 0 = 100 GeV where the mixing angle θ is such that cos(2θ) and therefore the D 2 D 1 H coupling falls to zero (see Fig. 4(a) red line). This decrease does not occur for the ∆m 0 = 1 GeV case (see Fig. 4(a) blue line), since ∆m + > ∆m 0 . However, this reduced contribution in Fig. 10(b) is compensated by the increase in D 2 → W ± D ∓ and D 2 → ZD . Finally, there is a boost in the 3-lepton r-value for ∆m + > 100 GeV where the D 2 D 1 H coupling becomes more enhanced.

C. Complementarity between the LHC and Non-collider Experiments
Finally we comment on the cosmological aspects of these models, first presenting relic density constraints on the i2HDM, then followed by relic density and direct detection constraints on the MFDM model.
The velocity-averaged annihilation cross section is dominated by DM scatterings to SM particles via the Z boson or the Higgs or via the quartic interaction D 1 D 1 V V . For the i2HDM, the relevant coupling that controls the annihilation cross section here is λ 345 , described in section II. For the MFDM, instead the annihilation cross section is controlled by the Yukawa coupling Y DM , Eq. (9) and the D 2 −D 1 mixing angle Eq. (6). Ignoring direct detection constraints, the i2HDM can account for the relic density in three distinct regions of parameter space [24]. The first is the so-called Higgs funnel region m D 1 m h 2 , the second around m D 1 72 GeV and finally the heavy mass region m D 1 ≥ 500 GeV. The Higgs funnel region is essentially a resonance annihilation regulated by the width and the coupling λ 345 and constrains the coupling to 10 −4 < λ 345 < 10 −2 if relic density has to be exactly satisfied. For the region around m D1 72 GeV, the annihilation proceeds primarily through the D 1 D 1 V V quartic interaction, via a pure gauge coupling. As m D 1 approaches the WW mass threshold, this annihilation is over efficient, and results generically in under abundant dark matter. Finally, there is a region of parameter space with m D 1 > 500 GeV with destructive interference between the quartic interaction and a t-channel diagram involving heavier i2HDM states. Moving on to direct detection considerations, we realise that since DM-nucleon scattering cross section scales as λ 2 345 , almost all of the parameter space in the low mass region is highly constrained by Xenon100 [47] and XENON1T [48]. Note that the constraints on the invisible Higgs force the coupling λ 345 ≤ 6 × 10 −3 . Therefore the only viable region of parameter space consistent with relic density and direct detection constraints remaining corresponds to the Higgs funnel region.
In view of the above discussions, we set λ 345 0, such that the only region where relic density is satisfied is the high mass region via pure quartic gauge couplings. In principle the Higgs funnel region exists for non zero λ 345 , and this region should be treated with caution. Meanwhile, for the MFDM model, DM scattering through the Z boson is suppressed due to a the D 1 − D mass split. Therefore in MFDM, direct detection constraints rely instead on DM scattering via Higgs exchange, which is proportional to the Yukawa coupling Y DM (see Eq. (9)) and D 2 − D 1 mixing angle (see Eq. (6)).
We now present the results from Fig. 8 combined with the non-collider constraints on the i2HDM DM masses and mass splitting. Fig. 11 shows the excluded regions from relic density (RD: orange) LEP (grey) and LHC limits (pink). The allowed regions include predicted DM that explains <95% of observed DM (blue) and DM that explains 95%-100% of observed DM (green). Since we assumed that λ 345 ∼ 0 in section II, the HD 1 D 1 vertex necessary for DM scattering is not permitted and direct detection results therefore do not constrain this phase space for our scenario. As discussed in the context of Fig. 11, the LEP limits cover the majority of the LHC exclusion region, with the exception of a small region with m D 1 > 40 GeV, 2 < ∆m + < 8 GeV for ∆m 0 = 1 GeV, and a region around ∆m + > 60 GeV, m D 1 < 50 GeV for ∆m 0 = 1, 10 GeV. The relic density limit excludes regions around m D 1 < 70 GeV before co-annihilation of D 2 − D 1 through Z and D + D 1 through W + channels (proportional to the gauge weak coupling) open up, and m D 1 > 600 GeV, ∆m + < 9 GeV for ∆m 0 = 1 GeV. This second region then reduces to m D 1 > 1000 GeV, ∆m + < 3 GeV for ∆m 0 = 10 GeV and even further for ∆m 0 = 100 GeV. For large DM mass but small mass split, the relic density remains too large for relic density constraints. However, for large DM masses and mass split, direct annihilation through W W and ZZ opens, dropping the relic density.
In Fig. 12, we now present the results from Fig. 9 combined with the non-collider constraints on the MFDM model DM masses and mass splitting. Again, the allowed regions include predicted DM that explains <95% of observed DM (blue) and DM that explains 95%-100% of observed DM (green). Relic density (orange), LEP limits (grey) and LHC limits (pink) are also displayed presented in Fig. 12. For the MFDM model we now include regions excluded from direct detection (DD: red) since we allow D 1 D 1 H interactions for this model. As discussed in the context of Fig. 9, the new LHC limits extend the phase-space coverage significantly beyond the LEP limits, including a substantial portion of the Higgs-funnel region around m D 1 = 60 GeV.
The direct detection exclusion region from XENON1T [48] limits an increasing range of ∆m + with increasing ∆m 0 as Yukawa coupling Y DM (see Eq. (9)) increases, since direct detection in the MFDM model comes from DM scattering with Higgs. This is approximately within the region ∆m + > 10 GeV and m D 1 < 350 GeV for ∆m 0 = 1 GeV in Fig. 12(a) and ∆m + > 6 GeV, for ∆m 0 = 10 GeV in Fig. 12(b). However, when ∆m 0 is sufficiently large (Fig. 12(c)) there is little change in this exclusion region, above ∆m + > 6 GeV, as the Y DM coupling relevant for direct detection converges towards √ ∆m 0 . In the region of m D 1 > 800 GeV the relic density excludes a larger range of m D 1 with increasing ∆m + as the DM co-annihilation of DD + through W and of D 1 D through Z becomes suppressed, enhancing the relic density above the exclusion limit for smaller values of m D 1 . However, starting with m D 1 = 100 GeV, as the DM mass increases, the region of excluded ∆m + decreases due to an enhancement of the D 1 D − W + and D 1 D Z couplings proportional to cos θ. The same couplings are relevant for DM self-annihilation through D ± to W + W − or through D to ZZ respectively. This trend continues until the required relic density is too high as m D 1 > 300 GeV, for which less DM annihilation suppression is needed. This happens roughly within the regions m D 1 > 600 GeV, ∆m + < 15 GeV and m D 1 > 80 GeV ∆m + > 15 GeV for ∆m 0 = 1 GeV in Fig. 12(a). Then, for ∆m 0 = 10 GeV, the lower limit for m D 1 < 600 GeV increases to ∆m + > 20 GeV, and the limit for smaller ∆m + < 10 GeV changes to m D 1 > 900 GeV in Fig. 12(b). This excluded space reduces further for ∆m 0 = 100 GeV in Fig. 12(c) to m D 1 > 1050 GeV within ∆m + < 30 GeV. This effect is due to increasing ∆m 0 and is attributed to an enhancement of the DM co-annihilation of DD + through W and of D 1 D through Z, with couplings proportional to cos θ (see Fig. 4(b)).

V. CONCLUSIONS
In this paper, we have explored the full parameter space relevant to 2-and 3-lepton signatures at the LHC emerging from two representative minimal consistent scenarios with scalar and fermion DM: the i2HDM and MFDM model, respectively. In our analysis, we have suggested a new parametrisation for both frameworks in terms of the DM mass, m D 1 , the mass difference between it and its charged multiplet partner, ∆m + = m D + − m D 1 , as well as the mass difference between D + and D 2 (the next-to-lightest neutral Z 2 -odd particle), ∆m 0 = m D 2 − m D + . This parametrisation allowed us to understand and interpret better the properties of the models under study and visualise a no-lose theorem covering their full parameter spaces. This approach is generic and quite model-independent since the mentioned mass differences are related to the couplings of the DM particle to the SM sector.
Our numerical results have led to the newest and most up-to-date LHC limits on the i2HDM and MFDM model parameter spaces, coming from the complementary combination of the afore-mentioned 2-and 3-lepton signatures. In particular, in the case of the i2HDM, we have found new regions which the LHC can cover above the LEP limits, for example: a) for a small ∆m 0 1 GeV, the LHC excludes regions with 45 GeV m D1 55 GeV and ∆m + 8 GeV as well as with 10 GeV m D1 50 GeV and 60 GeV ∆m + 95 GeV; b) for a larger ∆m 0 10 GeV, the LHC excludes regions with 10 GeV m D1 50 GeV and 50 GeV ∆m + 85 GeV. Furthermore, in the case of the MFDM model, the LHC probes a notably larger DM parameter space, in comparison, for example: a) for ∆m + 10 GeV, the LHC excludes DM masses up to 150-160 GeV, hence, going well beyond not only the LEP limits but also the LHC ones in mono-jet; b) for a much larger ∆m + 150 − 200 GeV, the LHC excludes DM masses up to 100 GeV, which is also well beyond such LEP and LHC constraints. Therefore, the parametrisation that we have suggested thus allows one to clearly establish a no-lose theorem for the MFDM model parameter space: the increase of ∆m 0 leads to the increase of the W ∓ D 1 D ± coupling and respective increase of the D 1 D ± production rate while kinematically suppressing the D 2 D ± production rate. These two trends roughly compensate each other making the exclusion picture in the (m D 1 , ∆m + ) plane quasi-independent of ∆m 0 .
We have also found another important complementarity between the 2-and 3-lepton signatures. The first one mainly covers low values of ∆m + while the second one is enhanced in the region of large ∆m + values. To corroborate our quantitative results, we have implemented and validated a 8 TeV ATLAS multi-lepton analysis into the CheckMATE package and made available this analysis publicly to the community. We have then created a map of the LHC efficiencies and cross-section limits for such 2-and 3-lepton signatures for the simple parametrisation of the parameter space that we have suggested, which would then allow for a quick and easy model-independent reinterpretation of the LHC limits for analogous classes of models.
Finally, we have produced combined constraints from the LHC, DM relic density and DM direct search experiments delineating the current status of the i2HDM and MFDM model parameter spaces. These combined limits indicate yet another important complementarity, the one between non-collider experiments and the LHC, both of which will continue to probe unknown territory of DM parameter space, hopefully resulting in a DM signal discovery. In such a quest, we deem the multi-lepton signatures studied here to be of the utmost importance at the LHC.

Appendix A: 8 TeV Validation: i2HDM
Appendix A details the validation of our CheckMATE recast for 8 TeV LHC exclusion limits for 2-lepton final states by comparing with the existing MadAnalysis implementation [49]. The CheckMATE analysis code was written based on and validated with the original experimental results, which searched for direct production of charginos, neutralinos and sleptons in final states with two leptons and missing transverse momentum [50] at the LHC. This was implemented using the CheckMATE's AnalysisManager in the current public build, and the SUSY analysis is available at [32]. The cutflows are given in tables II and III.
Global Cut   [51] was also written for the CheckMATE analysis performed here, as it also looks for final states with two leptons and missing energy. The public code of the Higgs analysis is available at [33]. The cutflow is given in table IV.
The events used for the validations were generated with CalcHEP, with 100000 events produced for each 9 benchmark points, using the SLHA files provided from HepData[52] and 1 benchmark point for HZ → invisible, with M H = 125.5 GeV. Leptonic decays of Z in the HZ production, χ ± χ2 production and W in the χ + χ − , were also specified in CalcHEP to improve efficiency, which were then showered with CheckMATE's built-in PYTHIA8 module. Detector effects are also applied via CheckMATE with a DELPHES module. Validation for the SUSY analysis is available at [53] and Higgs analysis available at [54].
The motivation behind fixing m D+ is because it is mostly only important for D + D − production, which give an additional EW coupling factor compared to D 2 D 1 production, only providing  significant contributions to r-value at very light m D+ . The ZD 1 D 1 production is also less dominant as the ZZD 1 D 1 coupling is quadratic and therefore weak compared to other couplings. The lowest allowed LEP limit of m D+ = 85 GeV is used and the higher value of m D+ = 150 GeV for comparison. The r-value contour plots for m D+ = 85 GeV in Fig. 13(a) for the W W b signal region shows that the Run 1 ATLAS 2-lepton analysis excludes for the lightest DM mass of m D1 ≤ 35 GeV for m D2 = 100 GeV, reaching a maximum of ∼ 40 GeV as m D2 is increased. In the case of Fig. 13(b), the W W c signal region reports stronger limits, excluding the lightest DM mass of m D1 ≤ 40 GeV for m D2 = 100 GeV, reaching a maximum of ∼ 45 GeV as m D2 is increased. For m D1 → 0, both cases exclude up to a mass of m D2 = 130 GeV.
By increasing m D + to 150 GeV in Fig. 14(a) increased.
Constraints increase with larger m D+ in both W W b and W W c signal regions, due to increased contributions from D + D − production and from D 2 D + production.
For larger m D2 , the signal events coming from D 2 D 1 production is smaller as the Z decay from D 2 is replaced in favour of W decay to D + and its decay to an additional W + and D 1 , which produces much softer leptons than required by the signal region cuts.
As m D2 is increased, m D1 is further constrained when considering the phase space above m D2 − m D1 = m Z , due to harder lepton production Z decaying from D 2 . Beyond this line, for a large enough m D2 − m D1 mass splitting, the Z-veto required by the SUSY analysis is no longer fulfilled by the signal, as real Z decay emerges in the production. Instead, the Z-window required by the Higgs analysis accepts these events, and becomes the dominant signal region as the mass splitting is further increased, independent of m D+ . From Fig. 13(a) and 14(a) the W W b+Higgs analyses agree with the general shape in [17] Fig.  1, for both m D + = 85 GeV and 150 GeV, but with lower r-value overall. However, results with larger r-value are obtained when considering W W c+Higgs signal regions in Fig. 13(b) and 14(b). While [17] considered both W W b and the HZ →invisible, they did not consider the W W c signal region. This is because, although the W W c signal region gives a larger observed r value than W W b, the expected r value is lower than W W b (deeming it the better channel by analysis tools). It is worth noting that the MadAnalysis validation for this signal region is overestimated compared to the experimental paper by a small amount, while the CheckMATE analysis implemented is closer to the experimental findings for survived number of MC events. contour exclusion limits for W W c, where a higher m T 2 cut on the leptons is implemented, shows higher observed r-value than the original paper's W W b contour exclusion limits.

Appendix B: i2HDM 13 TeV, 2-or 3-lepton Final states
In appendixB the analysis in [17] is then extended to higher centre of mass energies, using all available ATLAS and CMS 13 TeV analyses in CheckMATE. There is no equivalent to the Higgs ZH → + − + invisible recasted code as of writing, which is why the r-value beyond the m D2 − m D1 = m Z line is negligible. Extended to 13 TeV, Fig. 15 shows that this does not necessarily improve results, due to scaling of vector boson backgrounds in this phase space. Also shown between Fig. 15(a) and (b) is the improvement to r-value due to the inclusion of 3-lepton final states. This introduces 6 additional contributing diagrams where D 2 ,D + production (a) plays an important role. D 2 provides two leptons via D 2 → Z( + − ), D 1 and D + gives one extra lepton via D + → W + ( + ν), D 1 which contributes in heavier D 1 (and thus softer leptons) phase spaces than 2-lepton exclusive searches.
If D + is heavy enough, it can decay via D 2 , but this only occurs in Fig 16. The m D1 extends from 40 GeV to 45 GeV limit, while m D2 can extend from 120 GeV to 125 GeV.  Fig. 16 is the improvement from increasing lepton multiplicity for the larger m D+ value. Although less apparent than in Fig. 15, there are extensions to exclusion limits from under m D1 = 40 GeV to above the 40 GeV line and similarly small extensions to the m D2 limit.

Shown in
The diagrams with three leptons in the final state can now contribute to the relevant phase space, but as this chain of decays contains more steps than other contributions, the soft leptons it produces mostly do not pass the cuts here. Contrary to the 8 TeV case, there is not much increase in limits when moving to heavier D 2 . In fact, for 2-and 3-lepton searches, Fig. 15 In principle one would combine 13 TeV result with the 8 TeV exclusion limits for a comprehensive picture of LHC limits. Although 8 TeV limits are stronger than 13 TeV limits for the phase space in [17] and in this appendix, 8 TeV results do not improve within the phase pace of our 13 TeV results in new parametrisation, shown in section IV.

29
Appendix C: MFDM 8 TeV In appendix C 8 TeV r-value contour plots for MFDM are next discussed, starting with the same 8 TeV analyses that were used in the i2HDM case. Plots in Fig. 17 show fixed m D+ results, with scans in the m D 1 -m D 2 plane. The shaded region shows the Higgs to invisible excluded region, covering a large phase space of the CheckMATE excluded region. In figure 17(a), where m D + = 150 GeV, the region above m D 1 = 60 GeV has limits reaching to m D 1 < 85 GeV excluded while m D 2 > 180 excluded. As m D 2 increases, this increases the split between D 2 and D + /D , increasing the Yukawa coupling in Eq. (9). This facilitates more decays that would lead to leptonic final states, thus r-value is increased in the positive x-axis direction. In the y-axis, as DM mass is increased, the exclusion changes from atlas_higgs_2013_03 analysis at low DM mass, to the atlas_1403_5294 analysis at higher DM mass. Looking at Tables II, IV, this is due to the Z-window changing to a Z-veto with softer leptons being produced in association with harder DM for the same input energy. Fig. 17(b) shows m D + = 200 GeV, where no exclusion outside of the Higgs to invisible limit of m D 1 < 60 GeV excluded for m D 2 > 220 GeV. This is because m D + and m D become too heavy to produce at these masses, so the sources of leptonic final states are suppressed. The plots in Fig. 18 for MFDM at 8 TeV with fixed m D 2 , showing the m D 1 -m D + plane, are more analogous to those displayed for the i2HDM results shown previously. Figure 18(a) with m D 2 = 150 GeV is mostly excluded by Higgs to invisible limits, of m D 1 < 60 GeV excluded.
On the other hand, Fig. 18(b) limits, for m D 2 = 500 GeV, extend much further. It excludes a peak at m D 1 = 100 GeV, m D + = 180 GeV, and follows along M Z in terms of the mas split between m D + and m D 1 . This is when D can decay via real Z boson to two leptons, and D 1 in the final state. It is also close to when this mass split equals M W , facilitating two D ± decays via real W ± to a charged lepton and neutrino, along with D 1 in the final state.
In Fig. 18(b) there is an additional region of large r-value at 160 GeV < m D + < 225 GeV, with a limit of m D 1 = 65 GeV excluded, coming from the atlas_higgs_2013_03 analysis with harder leptonic decays, and lighter D 1 .

Appendix D: Numerical Overlaid Plots
In appendix D we present the exclusion plots overlaid with cross sections in f b for 2 or 3-leptonic final states. We first present the results for i2HDM in Fig.19 and MFDM in Fig.20. The magenta region and grey region indicate the current Higgs-to-invisible limit [46] of 0.15 branching fraction and LEP bounds on charginos for the fermionic DM case [45] respectively. 33

Appendix E: Sample Exclusion Formulae
In appendix E we describe the samples and formulae used to understand the tables presented in Appendix F,G. The input samples are separated between A,B and C, in both i2HDM and MFDM cases to improve efficiency.
For i2HDM, set A produces 50,000 events of D + D − and D 1 D 2 pairs, while specifying the decays D ± → ± , ν, D 1 (via W )and D 2 → + , − , D 1 (via Z). Set B produces 150,000 events of D ± D 2 , while specifying D 2 → + , − , D 1 (via Z) to obtain at least two leptons. Set C produces 100,000 events of the D 1 D 1 Z production, while specifying leptonic Z decay.
For MFDM, set A produces 50,000 total events, of D + D − and D 1 D pairs, while specifying decays of D ± → ± , ν, D 1 (via W ) and decays D → + , − , D 1 (via Z). Set B produces 150,000 D D 2 events, without specifying decays. Set C produces 100,000 events of D ± D and D ± D 2 pairs, while requiring D → + , − , D 1 decays, D 2 → ± , ν, D ± (via W ) or D 2 → − , + , D decays (via Z). This last set means D or D ± coming from D 2 can also decay hadronically, to fulfil more 2-lepton and 3-lepton thresholds of the analyses.
The exclusion limits between input samples, for the same signal region of a given analysis, can be related by the equation allowing for exclusion r values to be calculated, also noting that 34