aainstitutetext: Department of Physics, Indian Institute of Technology
Guwahati, North Guwahati, Assam - 781039, Indiabbinstitutetext: School Of
Physical Sciences, Indian Association for the Cultivation of Science, 2A and
2B, Raja S.C. Mullick Road, Kolkata 700032, Indiaccinstitutetext: Department
of Physical Sciences, Indian Institute of Science Education and Research
Kolkata, Mohanpur - 741246, India

# Distinguishing two dark matter component particles at $e^{+}e^{-}$ colliders

Subhaditya Bhattacharya, b Purusottam Ghosh a Jayita Lahiri c Biswarup
Mukhopadhyaya<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this work we demonstrate that a distorted double hump like missing energy
($\not{E}$) or missing transverse momentum ($\not{E}_{T}$) or missing mass
($\not{M}$) distribution at $e^{+}e^{-}$ colliders may hint towards the
presence of multipartite dark sector. We illustrate the phenomena using a two
component dark matter (DM) model involving an inert scalar doublet stabilised
under a $\mathcal{Z}_{2}$ symmetry providing a scalar DM, one vector like
fermion doublet and a right handed fermion singlet both stabilised under a
different $\mathcal{Z}^{{}^{\prime}}_{2}$ providing a fermion DM. We indicate
the region of parameter space where the production of the heavy charged
particles and their subsequent decay to DM yield double peak behaviour in
$\not{E}$ spectrum after satisfying DM constraints. Importantly, we illustrate
why and how $\not{E}$ serves as a better variable than $\not{E}_{T}$ in
distinguishing two component DM frameworks and therefore how International
Linear Collider (ILC) does better than the ongoing Large Hadron Collider
(LHC). We also chalk out a set of criteria to identify and segregate the
second peak in $\not{E}$ spectrum, after a careful analysis of the
corresponding Standard Model (SM) background contribution, which plays a
crucial role.

###### Keywords:

Multipartite Dark Matter, International Linear Collider (ILC)

##  1 Introduction

Evidences of dark matter (DM) have been accumulated from different
astrophysical observations like rotation curves of galaxies in galaxy cluster
Rubin:1970zza ; Zwicky:1937zza , collisions of galaxy clusters as in bullet
cluster Hayashi:2006kw , cosmological observations like anisotropy of cosmic
microwave background radiation (CMBR)Hu:2001bc  in WMAP Hinshaw:2012aka ;
Spergel:2006hy  or PLANCK Planck:2018vyg  satellite etc. Observations further
suggest that DM constitutes a large portion $(\sim 24\%)$ of the energy budget
of the universe, often expressed in terms of relic density $\Omega h^{2}\simeq
0.12$ Planck:2018vyg , where $\Omega={\rho}/{\rho_{c}}$ represents
cosmological density with $\rho$ being the DM density, $\rho_{c}$ the critical
density and $h$ represents Hubble expansion rate in units of 100 km/s/Mpc. In
spite of the overwhelming evidences, DM is yet to be discovered! A growing
belief that DM consists of fundamental particles largely explains the known
properties of DM like electromagnetic charge neutrality, stability over the
life time of the universe, collisionlessness or having weak to feeble
interaction with the visible sector. As none of the particles within the known
Standard Model (SM) of particle physics satisfy above, physics beyond the SM
(BSM) becomes inevitable.

Amongst DM genesis, two major ideas have been realised to account for correct
relic density; one where the DM is assumed to be in thermal and chemical
equilibrium in early universe to freeze-out when the rate of interaction falls
below the Hubble expansion rate and gives rise to the so-called Weakly
Interaction Massive particle (WIMP)Bertone:2004pz ; Roszkowski:2017nbc ;
Kolb:1990vq , while in the other case, DM is assumed to be out of equilibrium
in early universe owing to its feeble interaction with SM, and are produced
from the decay or annihilation of bath particles giving rise to feebly
interacting massive particle (FIMP)Hall:2009bx . Due to DM-SM interactions of
weak interaction strength to account for correct relic, WIMPs can show up in
different experiments like direct search via nuclear recoil or in collider
search experiments via missing energy and is the focus of this paper.

Further, given the limited information that we have from the dark sector, it
has been rather impossible to predict the nature of DM, whether WIMP or FIMP,
whether a spin zero scalar, a fermion or a vector boson etc. Bulk of the
analysis so far assumes the presence of one DM candidate, while the
possibility of having more than one DM component is quite possible,
particularly in analogy to the known visible sector having so many fundamental
particles and interactions. There lies the possibility of multiparticle dark
sector, which has been studied in many contexts, particularly to explain the
absence of DM signal in direct search experiments XENON:2018voc ;
PandaX-4T:2021bab . Two or more DM components are often thought of as WIMPs
Aoki:2012ub ; Liu:2011aa ; Cao:2007fy ; Bhattacharya:2013hva ; Esch:2014jpa ;
Karam:2016rsz ; Ahmed:2017dbb ; Poulin:2018kap ; Aoki:2018gjf ;
YaserAyazi:2018lrv ; Aoki:2017eqn ; Biswas:2013nn ; Bhattacharya:2016ysw ;
Bhattacharya:2017fid ; Barman:2018esi ; Bhattacharya:2018cgx ;
Bhattacharya:2019fgs ; Borah:2019aeq ; Chakraborti:2018lso ;
Chakraborti:2018aae ; Bhattacharya:2018cqu ; Yaguna:2021rds ; Belanger:2021lwd
; VanLoi:2021dzv ; Yaguna:2021vhb ; DiazSaez:2021pfw ; Chakrabarty:2021kmr ;
Nam:2020twn ; Betancur:2020fdl ; Nanda:2019nqy ; Bhattacharya:2019tqq ;
Elahi:2019jeo ; Herrero-Garcia:2018lga ; Das:2022oyx , while a combination of
WIMP-FIMP and SIMP-SIMP have also been studied Bhattacharya:2021rwh ;
DuttaBanik:2016jzv ; Choi:2021yps . We focus on a two component WIMP set up in
this paper.

Multicomponent DM although have been studied widely, but what can possibly
indicate the existence of multicomponent dark sector hasn’t been elaborated
greatly. Observation of a kink in the recoil energy spectrum have been
identified as a possible signal of two component WIMPs in direct
searchHerrero-Garcia:2017vrl ; Herrero-Garcia:2018qnz , while not many efforts
in context of collider searches is known to us excepting Hernandez-
Sanchez:2020aop 111Counting number of DMs from the end point of the spectrum
has been studied in Giudice:2011ib .. In this paper, we elaborate such a
possibility. Specifically, we are interested in a situation where two
different dark sector particles provide a similar signature. Now, this is
specifically because often a dark sector provides multiple possibilities of
signal, so two different signals can not be compared to draw an inference for
multipartite dark sector.

Amongst the major classifications of WIMP signal, we mainly focus on
multilepton plus missing energy channel. We elaborate how a distorted missing
energy spectrum, in particular, the presence of a second peak in the
distribution, may hint towards the presence of a second DM component. In fact,
identifying two peaks in measured distributions is a requirement in other
branches of science as well, and numerous studies in this direction have been
taken place. We exemplify this by choosing a sample two component model
framework involving a scalar doublet providing one DM component and a singlet
and a vector like doublet fermion to produce the second DM. To the best of our
knowledge, this hasn’t been addressed before in the literature.

It should be noted further that, we choose electron-positron collider to
elucidate such a possibility. The reason mainly lies in the fact that at
electron positron machine, the full kinematic information can be utilised
together with the fact that SM background contamination can be tackled in a
much better way utilising the beam polarisation. We will illustrate the
possibility in context of Large Hadron Collider (LHC) in a subsequent
analysis.

The paper is organised as follows: we first discuss WIMP signal at collider in
section 2 including the aspects of kinematics to single component and two
component DM models, in section 3 we discuss the model chosen for illustration
and DM constraints, followed by the selection of benchmark points in section
4, in section 5 we discuss the signal from two component dark sector with
corresponding SM background contribution, while in section 6 we analyse them
in detail, in section 7 we set up a few criteria for distinguishing the peaks
and exemplify in context of the chosen benchmark points, and finally we
summarise in section 8. In appendices A, B, C we provide some details omitted
in the main text.

##  2 WIMP signal at colliders

Owing to its weak interaction strength with the visible sector to account for
correct relic density, WIMPs can be produced at colliders 222It may be
contrasted to that of FIMPs, where only the particle in bath responsible for
DM production can be produced at collider leading to a disappearing charge
track Belanger:2018sti  as an indirect signal of DM.. But due to
electromagnetic charge neutrality together with stability, they will miss the
detectors. This will end up providing energy or momentum imbalance, which are
best characterised by the following kinematic variables:

  * •

Missing Transverse Energy or MET ($\not{E}_{T}$): The transverse momentum of
all electromagnetically charge neutral particles (DM or neutrinos) that miss
the detector is somehow called Missing Transverse Energy or MET and it can be
estimated from the momentum imbalance in the transverse direction associated
to the visible particles. Thus MET is defined as:

$\displaystyle\not{E}_{T}=-\sqrt{(\sum_{\ell,j}p_{x})^{2}+(\sum_{\ell,j}p_{y})^{2}};$
(1)

where the sum runs over all visible objects that include leptons, jets and
also the unclustered components.

  * •

Missing Energy or ME ($\not{E}$): The energy which is carried away by the
missing particles can simply be identified from energy of the visible
particles ($E_{\rm vis}$) given the CM energy ($\sqrt{s}$) of the reaction as

$\displaystyle\not{E}=\sqrt{s}-\sum_{\ell,j,\gamma}E_{\rm vis}.$ (2)

  * •

Missing Mass or MM ($\not{M}$): Lorentz invariant “missing mass” of the system
is defined as,

$\displaystyle\not{M}^{2}=(\sum
p_{\text{miss}})^{2}=\left(\sum_{i}p_{i}-\sum_{f}p_{f}\right)^{2}\,,$ (3)

which requires the knowledge of initial state four momentum ($p_{i}$) and
final state ones ($p_{f}$), where $f$ runs over all the visible particles
produced at the collider. For mono photon process
$e^{+}e^{-}\to\chi\bar{\chi}\gamma$, where $\chi(\bar{\chi})$ are DM, $m_{\rm
miss}^{2}=s-2\sqrt{s}E_{\gamma}$. Here, $E_{\gamma}$ is the energy of the
outgoing photon.

We note further that whenever DM (WIMP) is produced at collider, they will
produce all of the three variables MET ($\not{E}_{T}$), ME ($\not{E}$) and MM
($\not{M}$); however, $\not{E}$ and $\not{M}$ are measurable only in
$e^{+}e^{-}$ machines (where CM energy and initial four momentum of the
particles are known). In contrast, hadron colliders can only measure
$\not{E_{T}}$, absent the knowledge of subprocess CM energy ($\sqrt{\hat{s}}$)
of the collision. Note here that scalar sum of transverse momentum, sometimes
referred as Effective mass ($H_{T}$) can also be considered in hadron
colliders, but the predictions are similar to that of $\not{E_{T}}$. Before
entering into the merits of $\not{E}$ over $\not{E_{T}}$ we remind that either
$\not{E}$ or $\not{E_{T}}$ can tell us something meaningful in association
with some visible particles, against which the DM particle(s) recoil and can
be classified into two major signal categories:

$\bullet$ mono-X + $\not{E}~{}(\not{E_{T}}$), where X is a jet, a photon, a
weak boson or a Higgs.

$\bullet$ $n$-leptons + $m$-jets + $\not{E}~{}(\not{E_{T}}$), where $n,m$ can
be any integer.

The mono-X signature usually arises when two DM particles are produced
directly via a portal to dark sector Liew:2016oon ; Kahlhoefer:2017dnp ;
Boveia:2018yeb  or via DM effective operators Abercrombie:2015wmb ;
Abdallah:2015ter ; Barman:2021hhg , where $X$ comes out of initial state
radiation (ISR). Any dark sector consisting of WIMPs, necessarily provide
mono-X signal at collider. However, the kinematic observables at our disposal
in such a case are the four-momentum of particle X and the missing
momentum/energy/mass, which are mildly dependent on DM mass. We show that only
Missing Energy, which is measurable at ILC can have the seed of different DM
masses, that too in the end point of the spectrum not in the peak. Together
with the fact that the SM background is also large for such signals,
extraction of WIMP model/WIMP mass becomes harder.

The second category of final state, namely, multi-jet plus multi-lepton
channel usually arises when a pair of heavy particles (usually members of the
dark sector, we call it ‘Heavier Dark Sector Particles’ (HDSP)) are produced
at the collider, each of which further decays to a DM together with jets and
leptons in the final state. Obviously such an event topology, with a greater
multitude of visible particles will have richer kinematics. For example, the
kinematics in this case will depend on both HDSP and DM mass, a phenomenon
which we discuss in details, are better suited for differentiating the
contributions coming from different DM components. We may further note that
such signal may also arise in situations where two DM particles are produced
via separate decay chains involving HDSPs, as happens in SUSY. Also note that
leptonic channels with zero jets are often better signals to be probed at
Hadron collider for less background, while that with jets are suited at lepton
collider. It is such details that are invaluable in discriminating multi-
component DM particles.

Therefore $\not{E_{T}}~{}\rm{or}~{}\not{E}~{}\rm{or}~{}\not{M}$ distribution
will have direct correlation with the DM mass as well as the mass of the HDSP
(specifically in multilepton, multijet signal) as we elaborate here. In that
case, the distributions can in principle have multiple peaks in the presence
of more than one DM constituent. We propose this to be a tell-tale sign of
multi-component DM framework. However, distinguishability of those peaks is a
challenging task. We systematically address all these issues and suggest some
distinguishing criteria as we proceed. Before entering into all this, let us
discuss the various features of relevant kinematic observables for both one
and two component WIMP scenarios.

###  2.1 Aspects of kinematics

As mentioned above, we focus on multiplepton, multijet signal, a scenario in
which DM is produced at the end of a decay chain involving multiple
intermediate massive unstable particles. A look at the relevant kinematical
observables can give us information about the decay chain, such as, the mass
of DM particle ($m_{\rm DM}$) and mass difference with HDSP ($\Delta m$). Let
us examine the missing momenta distribution and its analytical expression in
terms of the parameters $\Delta m$ and $m_{\rm DM}$.

####  2.1.1 Single-component DM

Let us study the case of one DM candidate with mass $m_{\rm DM}$ and a HDSP
with mass $m_{\rm HD}=m_{\rm DM}+\Delta m$, which decays to DM with one/or
more massless SM fermions. For simplicity, we assume the two HDSPs are pair-
produced at rest 333Such a situation arises when HDSP mass is very close to
the CM energy of the collider ($\sqrt{s}\simeq 2m_{\rm HD}$), as we assume
mostly in the rest of the analysis.. In that case, the lab frame will be the
same as rest frame of the HDSP. At hadron colliders, one can only rely upon
MET ($\not{E_{T}}$) to study DM production. The maximum MET, which evidently
correspond to the end point of $\not{E_{T}}$ distribution, is obtained when
both DMs are moving in the same direction (as it is a vector addition, see
definition in Eq. 1), in particular when $z$-component of missing momentum is
minimal. With these conditions using four-momentum conservation, we arrive at
the following expression for maximum $\not{E_{T}}$ (for details, see Appendix
A):

$\not{E_{T}}^{max}\simeq\Delta m\left(1+r\right)\,;$ (4)

where $r=\frac{m_{\rm DM}}{m_{\rm HD}}$ is the ratio of the DM and HDSP
masses. It is clear from Equation. 4 that, $\not{E_{T}}^{max}$ depends on both
$\Delta m$ and $m_{\rm DM}$. In order to verify, we plot $\not{E_{T}}$
distribution in Fig. 1 for fixed $r=5$ and different $\Delta m=\\{48,98\\}$
GeV in the left, and for fixed $\Delta m=140$ GeV with different $m_{\rm
DM}=\\{100,350\\}$ GeV on the right plot. We use inert scalar doublet model
(detailed discussion is done in the next section)444The shape and end point of
the distribution depends on the model only when the number of intermediate
states differ, otherwise the kinematics remain similar excepting for the
cross-section of the process. to pair produce the singly charged scalar
components (HDSPs here) which further decays to DM with on/off shell $W^{\pm}$
bosons. The CM energy of the reaction is kept at $\sqrt{s}=500(1000)$ GeV, as
necessary to pair produce the HDSPs on-shell. From the left plot, we see that
the end-point shifts significantly towards a higher value of MET with the
increase of $\Delta m$ as suggested from Eqn. 4555To be specific, the end
point of the blue MET curve (with $\Delta m=192$ GeV, $m_{\rm DM}=48$ GeV) in
left hand side of Fig. 1 is at $\sim$200 GeV, whereas the end point is
estimated around $\sim$192 GeV using Eqn. 4; similarly the end point of the
purple plot ( with $\Delta m=392$ GeV, $m_{\rm DM}=98$ GeV) is seen at
$\sim$380 GeV, while the estimated value is at 392 GeV.; thus validating the
estimate to be in the correct ballpark. More importantly, we see that the peak
of the MET distribution shifts to larger values with larger $\Delta m$. From
the right plot we see a relatively weaker dependence on the DM mass $m_{\rm
DM}$. This further validates our previous claim that MET distribution is more
likely to distinguish $\Delta m$ and has much weaker dependence on $m_{\rm
DM}$, thus making multilepton, multijet channel more suitable for
distinguishing multicomponent DM signals.

Figure 1: Normalized $\not{E_{T}}$ distribution for different $\Delta m$ and
same $r$ (left) and different $m_{\rm DM}$ and same $\Delta m$ (right) [see
figure inset for details] to pair produce the singly charged scalar components
(HDSP) of an inert scalar doublet at $e^{+}e^{-}$ collider which further
decays to DM in association with on/off shell $W^{\pm}$ bosons at
$\sqrt{s}=500(1000)$ GeV, as necessary to pair produce the HDSPs on-shell.

We examine next the missing energy (ME: $\not{E}$) distribution, which is a
crucial observable in context of the lepton colliders, and compare it with the
$\not{E_{T}}$ distribution. ME for each DM particle is

$\not{E}=\sqrt{{|\vec{p}|^{2}_{\text{ DM}}}+m_{\rm DM}^{2}}\,;$ (5)

where $|\vec{p}|_{\text{D}M}$ is the magnitude of its three-momentum
associated to DM. Naturally, for an event with two DM (of same mass) produced
in a decay chain, ME of the event will be:

$\not{E}=\sqrt{{|\vec{p}_{1}|^{2}_{\text{ DM}}}+m_{\rm
DM}^{2}}+\sqrt{{|\vec{p}_{2}|^{2}_{\text{ DM}}}+m_{\rm DM}^{2}}\;.$ (6)

Importantly, due to the presence of DM mass $m_{\rm DM}$ in Eq. 5 and in Eq.
6, ME distribution will have a seed of DM mass unlike MET. To illustrate the
aforementioned feature, we plot missing energy $\not{E}$ for the different
$\Delta m$ and same $m_{\rm DM}$ on the left panel and for same $\Delta m$ and
different choices of $m_{\rm DM}$ in the right panel of Fig. 2. In the left
panel, we see that shift of ME spectrum peak to the left side with larger
$\Delta m$, i.e. different $\Delta m$ can be distinguished from ME spectrum
similar to MET. Importantly, in the right panel the $\not{E}$ distribution
shows a clear distinction with the distribution getting sharper as well as the
peaks also shifting to larger $\not{E}$ with larger $m_{\rm DM}$. Therefore,
ME variable, as expected from Eq. 5, depends on DM mass $m_{\rm DM}$, a
feature which is not observed in MET distribution (compare the right panel
plot in Fig. 1). Also, the end point of the ME distribution on the right side,
really shows the CM energy of the reaction, which is kept $\sqrt{s}=1$ TeV in
this case. Similar to Fig. 1, the plots are generated to pair produce singly
charged scalar components (HDSP) of an inert scalar doublet at $e^{+}e^{-}$
collider which further decays to DM in association with on/off shell $W^{\pm}$
bosons. We also note while MET distribution peaks at lower values of
$\not{E_{T}}$, ME distribution peaks at higher values of $\not{E}$, which is
evident from the definitions of these variables. Therefore, $\not{E}$
distribution is likely to be a better kinematic variable to elucidate the
underlying DM model than $\not{E_{T}}$ as ME is sensitive to both $m_{\rm DM}$
and $\Delta m$, while MET is sensitive to $\Delta m$ primarily.

Figure 2: Normalised $\not{E}$ distribution for same DM mass and different
$\Delta m$ (left) and for same $\Delta m$ and different DM mass (on right, see
figure insets for details) to pair produce singly charged scalar components
(HDSP) of an inert scalar doublet at $e^{+}e^{-}$ collider which further
decays to DM in association with on/off shell $W^{\pm}$ bosons at
$\sqrt{s}=1000$ GeV. Figure 3: Normalised $\not{M}$ distribution for two cases
with same DM mass and different $\Delta m$ (see figure insets for details) to
pair produce singly charged scalar components (HDSP) of an inert scalar
doublet at $e^{+}e^{-}$ collider which further decays to DM in association
with on/off shell $W^{\pm}$ bosons at $\sqrt{s}=1000$ GeV.

Finally we plot normalised missing mass ($\not{M}$) distribution for same DM
mass and different $\Delta m$ in Fig. 3 for the same model. One can see that
the distribution remains identical to that of the $\not{E}$. The definition of
$\not{M}$ reduces to

$\displaystyle\not{M}=\sqrt{\not{E}^{2}-|{\sum_{i}}\vec{p}_{\ell_{i}}|^{2}}\,;$
(7)

in context of hadronically quiet leptonic signal. This is similar to
$\not{E}$, unless $|{\sum_{i}}\vec{p}_{\ell_{i}}|$ is very large. However, we
have checked that such a scenario can rarely occur in the entire allowed
limited kinematic range of an $e^{+}e^{-}$ collider. Therefore, we won’t use
$\not{M}$ in our analysis and focus mainly on $\not{E}$ as our primary
observable. We discuss $\not{M}$ and compare it with $\not{E}$ for various
regions of our parameter space together with different choices of CM energy in
Appendix B.

####  2.1.2 Two-component DM

Having discussed the relevant kinematic features of a single-component DM, we
proceed now to the two-component scenario which is the focus of our study.
While the details of a two component DM model will be taken up in the next
section together with its cosmological possibilities, here we focus on a
simple situation having two scalar DMs stemming from two inert scalar doublets
with masses $m_{\rm DM_{1}},m_{\rm DM_{2}}$ and mass splitting with the
corresponding HDSPs as ${\Delta m}_{1}~{}{\rm and}~{}{\Delta m}_{2}$
respectively. Singly charged pairs (HDSP) are produced at $e^{+}e^{-}$
collider and decay to the respective DM components together with leptonically
decaying on/off-shell $W^{\pm}$. Subsequent ME and MET distributions for
different choices of $m_{\rm DM}$ and $\Delta m$ are studied here and some
broad features emerge out of the discussion, which we take up further in the
analysis. Let us consider a few possibilities as below:

  * •

${\Delta m}_{1}\approx{\Delta m}_{2}$ and $m_{\rm DM_{1}}<m_{\rm DM_{2}}$

In this case, the $\not{E_{T}}$ distribution for the two cases is expected to
merge since ${\Delta m}_{1}\approx{\Delta m}_{2}$, whereas the $\not{E}$
distributions are expected to show a double peak, owing to the dependence of
$\not{E}$ on the DM mass with $m_{\rm DM_{1}}<m_{\rm DM_{2}}$. We can see this
feature in Fig. 4, where we have considered, ${\Delta m}_{1}={\Delta
m}_{2}=100$ GeV, $m_{\rm DM_{1}}=50$ GeV and $m_{\rm DM_{2}}=390$ GeV. The
difference in DM masses for the two candidates is showing up in the form of
two peaks in the right hand side of Fig. 4 for $\not{E}$ distributions, while
$\not{E_{T}}$ distribution (on the left panel) doesn’t show any distinction
from that of a single component case 666For simplicity, we assume here the
production cross-sections for the HDSP pairs in both cases to be equal. Cross-
sections indeed play a vital role, which will be taken up in context of a
specific model.. Thus the difference in $m_{\rm DM}$ in two-component DM case
leads to a distinguishability in the $\not{E}$ distribution, even for
comparable $\Delta m$.

Figure 4: Normalized (a) $\not{E_{T}}$ on the left and (b) and $\not{E}$
distribution on the right for ${\Delta m}_{1}={\Delta m}_{2}=100$ GeV, $m_{\rm
DM_{1}}=50$ GeV and $m_{\rm DM_{2}}=390$ GeV. The production cross-sections
for the two HDSP pairs are assumed to be equal here.

  * •

${\Delta m}_{1}<{\Delta m}_{2}$ and $m_{\rm DM_{1}}\approx m_{\rm DM_{2}}$

We consider next a scenario, where the two DM sectors are characterized by
different mass splittings but similar DM masses. As an example, we chose a
specific case where $m_{\rm DM_{1}}=m_{\rm DM_{2}}=100$ GeV and ${\Delta
m}_{1}=10$ GeV and ${\Delta m}_{2}=380$ GeV. This is shown in Fig. 5. Here,
unlike the previous case, $\not{E_{T}}$ distribution shows a double peak
nature (see the left plot of Fig. 5). However, $\not{E}$ distribution produces
a better distinction where the second peak is even more prominent, while being
well separated, see right panel of Fig. 5.

Figure 5: Normalized (a) $\not{E_{T}}$ on the left and (b) and $\not{E}$
distribution on the right for ${\Delta m}_{1}=10$ GeV and ${\Delta m}_{2}=380$
GeV, $m_{\rm DM_{1}}=m_{\rm DM_{2}}=100$ GeV. The production cross-sections
for the lighter HDSP pairs are assumed to be 50% of that of the heavier HDSP
here.

  * •

${\Delta m}_{1}<{\Delta m}_{2}$ and $m_{\rm DM_{1}}<m_{\rm DM_{2}}$

We further consider the case where both DM mass and splitting with the
corresponding HDSPs are different, following a hierarchy ${\Delta
m}_{1}<{\Delta m}_{2}$ and $m_{\rm DM_{1}}<m_{\rm DM_{2}}$. As an example, we
consider ${\Delta m}_{1}=100$ GeV, $m_{\rm DM_{1}}=50$ GeV and ${\Delta
m}_{2}=$ 200 GeV, $m_{\rm DM_{2}}=390$ GeV and corresponding
$\not{E}$($\not{E_{T}}$) distributions are shown in Fig. 6. Here also we see
that the $\not{E_{T}}$ distribution (on the left) is almost overlapping with
each other for each DM and therefore show a single peak. On the other hand,
the $\not{E}$ (on the right) distribution shows two peaks coming from two
different DM candidates with different masses and splittings. Again for
simplicity, the production cross-sections for the heavier HDSP candidate is
assumed to be $1/4$th of that of the lighter one. The difference in $\Delta m$
thus accentuates the distinction between two types of DM particles when
$\not{E}$ is studied.

Figure 6: Normalized (a) $\not{E_{T}}$ and (b) and $\not{E}$ distribution for
${\Delta m}_{1}=100$ GeV, $m_{\rm DM_{1}}=50$ GeV and $m_{\rm DM_{2}}=700$
GeV, $\Delta m_{2}=200$ GeV. The production cross-sections for the two heavier
DM candidates is assumed to be $1/4$th of that of the lighter DM candidate.

  * •

${\Delta m}_{1}<{\Delta m}_{2}$ and $m_{\rm DM_{1}}>m_{\rm DM_{2}}$

Finally we explore another possibility, where the hierarchy in the DM masses
is opposite of that of the mass splittings with corresponding HDSPs. For
illustration, we have chosen $m_{\rm DM_{1}}=390$ GeV, $m_{\rm DM_{2}}=50$
GeV, ${\Delta m}_{1}=100$ GeV and ${\Delta m}_{2}=400$ GeV to demonstrate our
case. A comparison between $\not{E_{T}}$ and $\not{E}$ distribution is shown
in Fig. 7. The peak on the right corresponds to the higher DM mass. Here also,
one can see that $\not{E}$ distribution shows a clear two-peak behaviour while
the $\not{E_{T}}$ distribution shows a mere distortion.

Figure 7: Normalized (a) $\not{E_{T}}$ and (b) and $\not{E}$ distribution for
${\Delta m}_{1}=100$ GeV, $\Delta m_{2}=400$ GeV, $m_{\rm DM_{1}}=390$ GeV and
$m_{\rm DM_{2}}=50$ GeV. The production cross-sections for the two heavier DM
is assumed to be $1/4$th of that of the lighter DM.

It is clear from the preceding discussion that multipartite DM scenario can
show up in terms of double peak in ME or MET distribution. However, energy
observable is more useful in distinguishing multiple peak behaviour as
$\not{E}$ is more sensitive to $m_{\rm DM}$ than $\not{E_{T}}$ is, in both the
limits, $m_{\rm DM}>>\Delta m$ and $\Delta m>>m_{\rm DM}$. At $e^{+}e^{-}$
colliders, the complete knowledge of total centre-of-mass energy of each
collision and the very possibility of having $\not{E}$ distribution can
therefore be advantageous over the LHC where one has to resort to
$\not{E_{T}}$, so long as DM particle productions are sufficiently copious and
kinematically allowed.

So far we have been concerned mainly with the relative position of the peaks.
Larger the splitting between $m_{\rm DM_{1}}$ and $m_{\rm DM_{2}}$, or
${\Delta m}_{1}$ and ${\Delta m}_{2}$, the better is the separation between
two peaks in $\not{E}$ distribution and can hint for two component DM.
However, another important factor is the relative heights of the peaks, which
are determined by the cross-section of the process. For clear distinction, the
two peaks better be of comparable size. Assuming for simplicity that the HDSPs
decay primarily via one intermediate state with almost 100% branching ratio,
we can estimate the size of the peaks from the pair production cross-section
of the respective HDSPs. In the discussion above, we choose the cross-section
by hand for simplicity, however the HDSP cross-section crucially depends on
the parameters $m_{\rm DM}$ and $\Delta m$, both of which are crucially
involved in segregating the peaks in resultant ME (or MET) distribution. The
production moreover also depends on the model crucially, whether the particles
produced are scalar, fermion or vector boson and the mediators which
contributes to the production process. We will refer to these features
quantitatively when we encounter the model. However, we must note that HDSP
production cross-section diminishes with larger $m_{\rm DM}$ and $\Delta m$,
so whichever component is heavier contributes with a smaller cross-section at
least kinematically. So, the suitable scenario for discriminating DM
components should be chosen in a way, that we do not lose cross-section
significantly due to kinematics. This is better addressed if one dark sector
is fermionic, while the other is scalar, as the fermion production cross-
section is larger than that of the scalar for a similar mass and same
mediators777Having both DM components scalar or fermion can also provide such
distinctive ME (MET) distribution, but with a more fine tuned choices of
$\Delta m$ and $m_{\rm DM}$.. We therefore study the pair production of HDSPs
in the form of singly charged scalars ($\phi^{\pm}$) and singly charged
fermions ($\psi^{\pm}$) and their subsequent decays to respective DM
components at $e^{+}e^{-}$ collider. We shall examine in specific how they can
yield distinguishable peaks in ME (MET) distributions with comparable cross-
sections after satisfying cosmological constraints, in context of an
illustrative model, which we elaborate next.

##  3 A model with two DM components

We illustrate a model example here which consists of two dark sectors having
two stable DM candidates in the form of: $(i)$ SDM: lightest neutral state of
an inert scalar doublet
$\Phi^{T}_{i}=\left(\begin{matrix}\phi^{+}&&\frac{\phi^{0}+iA^{0}}{\sqrt{2}}\end{matrix}\right)$
and $(ii)$ FDM: lightest neutral state that arises from an admixture of a
vector like fermion doublet
$\Psi^{T}=\left(\begin{matrix}\psi&&\psi^{-}\end{matrix}\right)$ and a right
handed (RH) fermion singlet $\chi_{R}$. Stability of both DM components is
ensured by additional ${\mathcal{Z}_{2}\otimes\mathcal{Z}^{\prime}_{2}}$
symmetry, where the inert scalar doublet transforms under $\mathcal{Z}_{2}$,
while the fermion fields transform under a different
$\mathcal{Z}^{\prime}_{2}$. The quantum numbers of additional fields are shown
in the Table 1. The minimal renormalizable Lagrangian for this model then
reads,

$\displaystyle\mathcal{L}\supset\mathcal{L}^{\rm SDM}+\mathcal{L}^{\rm FDM}.$
(8)

Fields | $\underbrace{SU(3)_{C}\otimes SU(2)_{L}\otimes U(1)_{Y}}$ $\otimes\underbrace{\mathcal{Z}_{2}\otimes\mathcal{Z}^{\prime}_{2}}$
---|---
SDM | $\Phi=\left(\begin{matrix}\phi^{+}\\\ \frac{1}{\sqrt{2}}(\phi^{0}+iA^{0})\end{matrix}\right)$ |  1 2 1 - +
FDM | $\Psi_{L,R}=\left(\begin{matrix}\psi\\\ \psi^{-}\end{matrix}\right)_{L,R}$ |  1 2 -1 + -
| $\chi_{{}_{R}}$ |  1 1 0 + -

Table 1: Charge assignment of the BSM fields under the gauge group
$\mathcal{G}\equiv\mathcal{G}_{\rm SM}\otimes\mathcal{G}_{\rm DM}$ in the toy
scenario. Here $\mathcal{G}_{\rm SM}\equiv SU(3)_{C}\otimes SU(2)_{L}\otimes
U(1)_{Y}$ and $\mathcal{G}_{\rm
DM}\equiv\mathcal{Z}_{2}\otimes\mathcal{Z}^{\prime}_{2}$.

The Lagrangian for the scalar dark sector (SDM), having inert scalar doublet
$\Phi$ can be written as :

$\displaystyle\mathcal{L}^{\rm SDM}$ $\displaystyle=$
$\displaystyle\Big{|}\Big{(}\partial^{\mu}-ig_{2}\frac{\sigma^{a}}{2}{W^{a}}^{\mu}-ig_{1}\frac{Y}{2}B^{\mu}\Big{)}{\Phi}\Big{|}^{2}-V(\Phi,H)\,;$
$\displaystyle V(\Phi,H)$ $\displaystyle=$
$\displaystyle\mu_{\Phi}^{2}(\Phi^{\dagger}\Phi)+\lambda_{\Phi}(\Phi^{\dagger}\Phi)^{2}+\lambda_{1}(H^{\dagger}H)(\Phi^{\dagger}\Phi)+\lambda_{2}(H^{\dagger}\Phi)(\Phi^{\dagger}H)+\frac{\lambda_{3}}{2}[(H^{\dagger}\Phi)^{2}+h.c.]~{}\,.$

Note here, we assume $\mu_{\Phi}^{2}>0$ which keep vacuum extension value
(vev) $\langle\Phi\rangle=0$ and ensure $\mathcal{Z}_{2}$ remains unbroken.
$H$ represents SM Higgs isodoublet. After Electroweak Symmetry Breaking
(EWSB), the SM Higgs acquires non zero vev, $\langle H\rangle=v~{}(246~{}{\rm
GeV})$ and generates physical states in the form of CP even state $\phi^{0}$
(mass $m_{\phi^{0}}$), CP-odd state $A^{0}$ (mass $m_{A^{0}}$) and charged
scalar states $\phi^{\pm}$ (mass $m_{\phi^{\pm}}$) with adequate mass
splitting between them. Using minimization conditions of the scalar potential
$V(H,\Phi)$ along different fields directions and the definitions of physical
states, one can obtain the following relations between the mass eigenvalues
and parameters involved in the scalar potential (for details, see Dolle:2009fn
; Bhattacharya:2019fgs ):

$\displaystyle\mu_{H}^{2}~{}=~{}\frac{m_{h}^{2}}{2},~{}~{}~{}~{}\mu_{\Phi}^{2}~{}=~{}m_{\phi^{0}}^{2}-\lambda_{L}v^{2},~{}~{}\lambda_{1}~{}=~{}2\lambda_{L}-\frac{2}{v^{2}}(m_{\phi^{0}}^{2}-m_{\phi^{\pm}}^{2})~{},$
$\displaystyle\lambda_{2}~{}=~{}\frac{1}{v^{2}}(m_{\phi^{0}}^{2}+m_{A^{0}}^{2}-2m_{\phi^{\pm}}^{2})~{}~{}{\rm{and}}~{}\lambda_{3}~{}=~{}\frac{1}{v^{2}}(m_{\phi^{0}}^{2}-m_{A^{0}}^{2})~{}~{};$
(10)

where $\lambda_{L}=\frac{1}{2}(\lambda_{1}+\lambda_{2}+\lambda_{3})$ and
$m_{h}$ represents the mass of SM-like neutral scalar found at LHC with
$m_{h}=125.09$ GeV. In our analysis we consider the mass hierarchy:
$m_{\phi^{0}}<m_{A^{0}},m_{\phi^{\pm}}$ and hence the lightest neutral state
$\phi^{0}$ with mass $m_{\phi^{0}}$ is a viable DM candidate888Note here that
other neutral state, $A^{0}$ can also be a viable DM candidate with
$m_{A^{0}}<m_{\phi^{0}},m_{\phi^{\pm}}$, which is equally possible by
adjusting model parameters.. The parameters relevant for DM phenomenology from
SDM sector are given by:

$\displaystyle\\{m_{\rm DM_{1}},~{}\Delta m_{1},~{}\lambda_{L}\\}.$ (11)

Here $m_{\rm DM_{1}}=m_{\phi^{0}},~{}\Delta
m_{1}=m_{\phi^{\pm}}-m_{\phi^{0}}\simeq m_{A^{0}}-m_{\phi^{0}}$, where we
assumed $m_{\phi^{\pm}}\simeq m_{A^{0}}$ for the sake of simplicity, which
does not alter DM phenomenology or subsequent collider analysis significantly.

Let us now turn to fermion DM component. The minimal renormalizable Lagrangian
for the fermion dark sector (FDM) having one vector like doublet ($\Psi$) and
one right handed singlet ($\chi_{R}$) reads Dutta:2020xwn :

$\displaystyle\mathcal{L}^{\rm FDM}$ $\displaystyle=$
$\displaystyle\overline{\Psi}_{L(R)}~{}[i\gamma^{\mu}(\partial_{\mu}-ig_{2}\frac{\sigma^{a}}{2}W_{\mu}^{a}-ig_{1}\frac{Y^{\prime}}{2}B_{\mu})]~{}\Psi_{L(R)}+\overline{\chi_{{}_{R}}}~{}(i\gamma^{\mu}\partial_{\mu})~{}\chi_{{}_{R}}$
(12) $\displaystyle-$ $\displaystyle
m_{\psi}\overline{\Psi}\Psi-\Big{(}\frac{1}{2}m_{\chi}\overline{\chi_{{}_{R}}}(\chi_{{}_{R}})^{c}+h.c\Big{)}-\frac{Y}{\sqrt{2}}\Big{(}~{}\overline{\Psi_{L}}~{}\widetilde{H}{\chi_{R}}+~{}\overline{\Psi_{R}}~{}\widetilde{H}{\chi_{R}}^{c}+h.c\Big{)}.$

In the above
$\Psi_{L(R)}=P_{L(R)}\Psi;~{}P_{L/R}=\frac{1}{2}(1\mp\gamma_{5})$. After EWSB,
a mixing is generated amongst neutral component of doublet $\Psi$ and singlet
$\chi_{R}$ thanks to the Yukawa coupling $Y$. The resulting mass matrix in the
basis of $({\psi_{R}}^{c},\psi_{L},{\chi_{{}_{R}}}^{c})^{T}$ is given by:

$\mathcal{M}=\left(\begin{array}[]{ccc}0&m_{\psi}&\frac{Yv}{2}\\\
m_{\psi}&0&\frac{Yv}{2}\\\ \frac{Yv}{2}&\frac{Yv}{2}&m_{\chi}\\\
\end{array}\right)\,,$ (13)

where $v=\langle H\rangle~{}(246~{}{\rm GeV})$ is the vev of SM Higgs. The
above symmetric mass matrix, $\mathcal{M}$ can be diagonalised by an unitary
matrix $\mathcal{U}$ to yield $\mathcal{M}_{\rm
diag.}=\mathcal{U}.\mathcal{M}.\mathcal{U}^{T}$ with

$\mathcal{U}=\left(\begin{array}[]{ccc}1&0&0\\\ 0&e^{i\pi/2}&0\\\ 0&0&1\\\
\end{array}\right)\left(\begin{array}[]{ccc}\frac{1}{\sqrt{2}}\cos\theta&\frac{1}{\sqrt{2}}\cos\theta&\sin\theta\\\
-\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0\\\
-\frac{1}{\sqrt{2}}\sin\theta&-\frac{1}{\sqrt{2}}\sin\theta&\cos\theta\\\
\end{array}\right)\\\ ,$ (14)

where $\theta$ is the singlet-doublet mixing angle, which in terms of the
model parameters is given as:

$\tan 2\theta=\frac{Yv}{\sqrt{2}(m_{\psi}-m_{\chi})}\,.$ (15)

The mass eigenvalues of the physical states are then:

$\displaystyle m_{\psi_{{}_{1}}}$ $\displaystyle=$ $\displaystyle
m_{\chi}\cos^{2}\theta+m_{\psi}\sin^{2}\theta-\frac{Yv}{\sqrt{2}}\sin
2\theta,$ $\displaystyle m_{\psi_{{}_{2}}}$ $\displaystyle=$ $\displaystyle
m_{\psi},$ (16) $\displaystyle m_{\psi_{{}_{3}}}$ $\displaystyle=$
$\displaystyle
m_{\psi}\cos^{2}\theta+m_{\chi}\sin^{2}\theta+\frac{Yv}{\sqrt{2}}\sin
2\theta,$ $\displaystyle m_{\psi^{\pm}}$ $\displaystyle=$ $\displaystyle
m_{\psi}~{}.$

For $Yv/\sqrt{2}<m_{\chi}\ll m_{\psi}$, we have
$m_{\psi_{{}_{3}}}>m_{\psi_{{}_{2}}}(\equiv m_{\psi^{\pm}})>m_{\psi_{{}_{1}}}$
making $\psi_{{}_{1}}$ the DM. Note here all the neutral states,
$\psi_{1,2,3}$ are Majorana fermions. The DM phenomenology of fermion sector
(FDM) mainly then depend on the following independent parameters:

$\displaystyle\\{m_{\rm DM_{2}},~{}\Delta m_{2},~{}\sin\theta\\}\,;$ (17)

where $m_{\rm DM_{2}}=m_{\psi_{1}},\Delta
m_{2}=m_{\psi^{\pm}}-m_{\psi_{1}}=m_{\psi_{2}}-m_{\psi_{1}}$. The parameters
$Y$, $m_{\psi}$ and $m_{\chi}$ are expressed in terms of above independent
parameters as:

$\displaystyle Y$ $\displaystyle=$ $\displaystyle\frac{\sqrt{2}~{}\Delta
m_{2}~{}\sin 2\theta}{v}\,,$ $\displaystyle m_{\psi}$ $\displaystyle=$
$\displaystyle m_{\psi_{1}}\sin^{2}\theta+m_{\psi_{3}}\cos^{2}\theta\,,$ (18)
$\displaystyle m_{\chi}$ $\displaystyle=$ $\displaystyle
m_{\psi_{1}}\cos^{2}\theta+m_{\psi_{3}}\sin^{2}\theta.$

Note here that the Yukawa coupling ($Y$) is proportional to mass splitting as
well as mixing angle, which plays an important role in the consequent DM
phenomenology.

Before moving on, we note further that the two dark sectors have no
renormalizable interactions in between999The interaction terms appear in
dimension five like
$(\overline{\Psi_{L}}\Psi_{R})(\Phi^{\dagger}\Phi),~{}(\overline{\chi_{{}_{R}}}(\chi_{{}_{R}})^{c})(\Phi^{\dagger}\Phi)$
with hermitian conjugates or other higher order terms., so that DM-DM
interactions mainly proceed via Higgs portal interaction after EWSB.

###  3.1 Constraints on the model parameters

The parameters of the model are subject to the following constraints:

  * •

The requirement of a stable vacuum at the electroweak scale.

  * •

Perturbativity at the electroweak scale.

  * •

Non-observation of the BSM particles at collider search, the most significant
of those here are $m_{\chi^{\pm}}>102.7$ GeV LEP  and $m_{\phi^{\pm}}>70$ GeV
Pierce:2007ut .

  * •

The requirement that the two DM components $\phi^{0}$ and $\psi_{1}$ saturate
the relic-density according to the Planck data Planck:2018vyg  at $2\sigma$
level to adhere to a two-component DM scenario, where

$\displaystyle\Omega_{\rm DM}h^{2}=\Omega_{\rm DM_{1}}h^{2}+\Omega_{\rm
DM_{2}}h^{2}=\Omega_{\phi^{0}}h^{2}+\Omega_{\psi_{1}}h^{2}\equiv 0.12\pm
0.001.$ (19)

  * •

Non-observation of a WIMP like DM in direct search experiments, the latest
being XENON-1T, where the spin independent DM-nucleon cross-section limits at
the $2\sigma$ level XENON:2018voc ; PandaX-4T:2021bab  goes as down as

$\displaystyle\sigma^{SI}_{\rm{min}}\lesssim 10^{-47}~{}{\rm{cm}^{2}}.$ (20)

Some details of these constraints are presented in Appendix C. In this section
we focus mainly on the DM constraints and show the allowed parameter space of
the model, which in turn provides some characteristic benchmark points that is
taken up for subsequent collider analysis.

Dark matter constraints

The most important constraint for the DM comes from the observed relic
density. In a two component model, the total DM relic density as observed in
PLANCK data Planck:2018vyg  will be the sum of the contributions from
individual components as

$\displaystyle\Omega_{\rm DM}h^{2}=\Omega_{\rm DM_{1}}h^{2}+\Omega_{\rm
DM_{2}}h^{2}\,,$ (21)

as mentioned before. In the equation above, $\Omega_{\rm DM_{i}}(i=1,2)$
represents the share of a particular DM component in the total relic density.
We may further recall that relic density of WIMP is inversely proportional to
its annihilation cross-section to visible sector particles or to the other DM,
which governs the relic under abundance or over abundance of the particular DM
component:

$\displaystyle\Omega_{\rm DM_{i}}\sim\frac{1}{\langle\sigma v\rangle}\,.$ (22)

Now, one may assume that the DM number density of any component around the
earth is proportional to the share of that component in the total relic
density. In such a situation, the upper limit on the scattering cross-section
of the $i$th DM component $\sigma^{i}\rm{(max)}$ with the detector nucleon for
direct search is related to $\sigma^{i}_{\rm{tot}}$, the cross-section if the
corresponding component were the sole WIMP DM constituent, by the fraction of
relic density with which it is present in a two component set up, following:

$\sigma^{i}\rm{(max)}={\frac{\Omega^{i}_{\rm{DM}}}{\Omega_{\rm{DM}}}}\sigma^{i}_{\rm{tot}}\,.$
(23)

It is the above direct search limit that should be applicable to the elastic
scattering cross-section for a given mass of a particular candidate in a
multi-DM scenario. Thus the direct search limit on this component is different
from the case where it is the only WIMP candidate.

Since any two-component WIMP DM is subject to both the limits from relic
density and direct search, it is useful to recall two more points. First, the
tension between the limit on the cross-section in direct searches and the
relic density limit is relaxed in scenarios where the dark sector permits co-
annihilation. Secondly, if we look at the fermionic DM component in our
illustrative model, the right-chiral fermion $\chi_{R}$ is introduced there.
Since $\chi_{R}$ has no Z-coupling (having $Q=Y=0$), its mixing with
$\psi_{0}$ suppressed the Z-induced contribution to the elastic scattering and
relieves the pressure from direct search limits. Annihilation, co-annihilation
and conversion processes for both DM components are mentioned in Appendix C.

In the following, we study the impact of relic density and direct detection
constraints on this two component set up and then choose suitable benchmarks
for our collider analysis. In order to find the consistent DM parameter space,
we implement the model in Feynrules Degrande:2011ua  and in MicrOmegas
Belanger:2018ccd  to perform a numerical scan on the relevant parameters in
the following range:

$\displaystyle 30~{}\leq m_{\phi^{0}}\leq 600~{}{\rm GeV},~{}1~{}\leq\Delta
m_{1}\leq 500~{}{\rm GeV},~{}0.01~{}\leq\lambda_{L}\leq 0.1,$ $\displaystyle
30\leq~{}m_{\psi_{1}}\leq 600~{}\rm GeV,~{}1~{}\leq\Delta m_{2}\leq 500~{}{\rm
GeV},~{}0.01~{}\leq\sin\theta\leq 0.3~{}.$ (24)

We present in Fig. 8, the allowed region of the model in the planes of $m_{\rm
DM_{1}}-\Delta m_{1}$ (left) and $m_{\rm DM_{2}}-\Delta m_{2}$ (right). Red
points satisfy relic density and the blue points satisfy additionally the
direct search bound from XENON1T. It is seen that for the scalar DM, relic
(under)abundance and direct detection bounds are both satisfied for almost the
entire range of $m_{\rm DM_{1}},~{}\Delta m_{1}$ scanned here. This is
achieved mainly due to two reasons: one, the gauge mediated depletion
processes (see Appendix C for details) ensure significant annihilation cross-
section, while suitable choice of the parameter $\lambda_{L}$, which governs
the interaction between the DM and SM Higgs, ensures the consistency with
direct detection bound. For FDM case (see right panel figure), both relic
density and direct search constraints are satisfied consistently for $\Delta
m_{2}\lesssim 20$ GeV. This is due to co-annihilation processes contributing
significantly to relic density with smaller $\Delta m_{2}$, which do not
contribute to DM elastic scattering for direct search. We also note that for
$m_{\psi_{1}}\simeq\frac{m_{h}}{2}$, in the Higgs resonance region, one can
only simultaneously satisfy both relic (under) abundance and direct detection
bound even for larger $\Delta m_{2}$, where the co annihilation contribution
is compensated by the resonance enhancement. We further note that there is a
region in the right panel of Fig. 8 marked in white, where none of the points
satisfy relic under-abundance, sandwiched between red and blue points. This
region actually provides relic over abundance with smaller annihilation cross-
section. Here with $\Delta m_{2}sim20-100$ GeV, the co annihilation
contribution turns drastically low, while the annihilation cross-section isn’t
large enough given the range of $\sin\theta$ chosen for the scan. For $\Delta
m_{2}\gtrsim 100$ GeV, the Yukawa coupling ($Y\sim\Delta m_{2}$) turns large
enough to provide required under abundance. Also note that the allowed
parameter space in this two component model most often has a dominant fermion
DM contribution over the scalar one, simply because satisfying relic under
abundance together with direct search constraint for fermion DM relies heavily
on co-annihilation contribution given a minimal $\sin\theta$, limiting the
minimum relic density that FDM can achieve.

Figure 8: Relic and direct search outcome for scalar DM(left) and fermion
DM(right). The red points satisfy relic density requirement and the blue
points satisfy direct detection bound as well as relic density constraints.
Figure 9: Regions allowed by indirect detection constraints from $W^{+}W^{-}$
mode for scalar and fermionic DM. Their individual contributions have been
scaled by their respective number density.

In Figure 9, we present the region allowed by indirect detection constraints.
Due to gauge interaction, the major annihilation channel for both scalar and
fermionic DM occurs into $W^{+}W^{-}$ final state. Therefore, the most
stringent upper limit comes from this particular channel in MAGIC
MAGIC:2016xys  \+ FERMI satellite data Fermi-LAT:2015att . Similar to direct
search, while calculating the annihilation cross-section of each DM into
$W^{+}W^{-}$ final state, we have scaled their individual contribution with
the fraction of respective DM densities to the total relic density. We can see
that, the entire range of DM parameters scanned here satisfy the upper limit
from indirect searches. So, the relic density and direct search allowed
parameter space of the model automatically addresses indirect search
constraints.

##  4 Selection of benchmark points

Having discussed all the constraints on the model, and carefully examining
various criteria for distinguishing two peaks in kinematic distributions for
two-component dark matter, we are in a position to choose a few benchmark
points, suited for the likely centre-of-mass energies of $e^{+}e^{-}$
colliders planned for the near future, including the International Linear
Collider (ILC). We note the relevant kinematic regions for the study where the
peaks are maximally distinguishable are:

  * •

Region I: $m_{\phi^{0}}>m_{\psi_{1}}$ and $\Delta m_{1}>\Delta m_{2}$ ,

  * •

Region II: $m_{\phi^{0}}<m_{\psi_{1}}$ and $\Delta m_{1}<\Delta m_{2}$ ,

  * •

Region III: $m_{\phi^{0}}<m_{\psi_{1}}$ and $\Delta m_{1}>\Delta m_{2}$ ,

  * •

Region IV: $m_{\phi^{0}}>m_{\psi_{1}}$ and $\Delta m_{1}<\Delta m_{2}$ ,

Allusions to these two kinematical regions have have already been made in the
previous section. However, choosing the benchmark points that show
distinguishable peaks at collider crucially depends on the relative cross-
sections of the two dark sector particles can contribute with. Here we make a
quantitative comparison between the scalar and fermionic HDSP’s in the context
of our illustrative model. The cross section of the process
$e^{+}e^{-}\rightarrow\phi^{+}\phi^{-}/\psi^{+}\psi^{-}$ as a function of
centre-of-mass (CM) energy ($\sqrt{s}$) is shown in Fig. 10. ILC is supposed
to start running with $\sqrt{s}=$ 500 GeV with a possible upgrade to 1 TeV.

Figure 10: Production cross-section of HDSP pair for scalar $\phi^{+}\phi^{-}$
(left) and fermionic $\psi^{+}\psi^{-}$ (right) case as a function of
$\sqrt{s}$ for different choices of HDSP masses for unpolarised beam at ILC.

It is clear from Fig. 10, the scalar DM has smaller production cross-section
compared to the fermionic DM at a given HDSP mass, because of the absence of
polarization sum in the former. Therefore, we would ideally choose a lighter
$\phi^{\pm}$ particle and heavier $\psi^{\pm}$, in order for them to have
comparable production rates. From this viewpoint, Region I proves to be
underwhelming. Also, we have seen earlier that small $\Delta m$ would imply a
narrow distribution. A narrow distribution is preferred for the scalar case,
as the production cross-section is in general less than the fermionic case.
Therefore, for the scalar DM to have a detectable/distinguished peak,
$m_{\phi^{\pm}}<m_{\psi^{\pm}}$ and small $\Delta m_{1}$ is preferred, which
means $\phi^{\pm}$ and $\phi^{0}$ almost degenerate. One has to also keep in
mind, the collider bound from LEP, which forces us to choose the mass of the
charged scalar $\phi^{\pm}>70$ GeV. In case of Region II and IV where $\Delta
m_{1}<\Delta m_{2}$, we restrict ourselves in the Higgs-resonance region for
the fermionic DM, ie $m_{\psi_{1}}\approx\frac{m_{h}}{2}$ to satisfy the
observed relic density as well as direct detection bound. Lastly, for the sake
of completeness, we also explore the Region III, where
$m_{\psi_{1}}>m_{\phi^{0}}$, and $\Delta m_{1}>\Delta m_{2}$, i.e. the
fermionic DM is degenerate with its HDSP and the scalar DM and its
corresponding HDSP have larger mass gap. In Tables 3 and 2, we identify two
pairs of benchmark points, each of which can be probed with $\sqrt{s}=500$ GeV
and 1TeV respectively. We emphasize further that we have chosen the benchmarks
in such a way, that we can probe different regions of the phase-space after
satisfying DM constraints, and at the same time explore their detectability at
different CM energies.

Benchmark | $\\{m_{\phi^{0}},~{}\Delta m_{1},~{}\lambda_{L}\\}$ | $\\{m_{\psi_{1}},~{}\Delta m_{2},~{}\sin\theta\\}$ |  $\sigma$(fb) |  $\sigma$(scalar) |  $\sigma$(fermion)
---|---|---|---|---|---
BP1 | 100, 10, 0.01 | 60.5, 370, 0.022 | 2.75 fb | 0.74 fb | 2.0 fb
BP2 | 100, 10, 0.01 | 58.91, 285, 0.032 | 3.3 fb | 0.74 fb | 2.6 fb
Table 2: Signal benchmark points and HDSP production cross-sections at ILC
with $P_{e^{-}}=0.8,P_{e^{+}}=-0.3$ polarization at $\sqrt{s}=1$ TeV.
$m_{A}=m_{\phi^{\pm}}$ is assumed for all the benchmarks. $WW$ background
production cross-section is 18 fb.

In Table 2, we show benchmark points BP1 and BP2 chosen from Region IV, where
$\Delta m$ for scalar DM is smaller than that of fermionic DM and
$m_{\phi^{0}}>m_{\psi_{1}}$. For both benchmarks, the fermionic DM is around
the Higgs resonance to satisfy relic density and direct search constraints.
BP1 and BP2 differ from each other in terms of different $\Delta m$ and
therefore their cross-sections are different. The HDSP masses for BP1 and BP2
are such that they can only be produced at $\sqrt{s}=1$ TeV.

Benchmark | $\\{m_{\phi^{0}},~{}\Delta m_{1},~{}\lambda_{L}\\}$ | $\\{m_{\psi_{1}},~{}\Delta m_{2},~{}\sin\theta\\}$ |  $\sigma$(fb) |  $\sigma$(scalar) |  $\sigma$(fermion)
---|---|---|---|---|---
BP3 | 100, 10, 0.01 | 58.87, 176, 0.04 | 7.8 fb | 2.3 fb | 5.5 fb
BP4 | 100, 10, 0.01 | 58.87, 190, 0.042 | 3.9 fb | 2.3 fb | 1.6 fb
Table 3: Signal benchmarks and HDSP production cross-sections at ILC with
$P_{e^{-}}=0.8,P_{e^{+}}=-0.3$ polarization at $\sqrt{s}=500$ GeV.
$m_{A}=m_{\phi^{\pm}}$ is assumed for all the benchmarks. $WW$ background
cross-section is 50 fb. All masses are in GeV.

In Table 3, we choose benchmark points which can be probed at 500 GeV run of
the proposed ILC. Here, $\Delta m$ for scalar DM is smaller than fermionic DM
and $m_{\phi^{0}}>m_{\psi_{1}}$. BP3 and BP4 have different $m_{{\psi}^{\pm}}$
and therefore yield quite different production cross-sections.

##  5 Signal and background

Figure 11: Dominant process contributing to $\ell^{+}\ell^{-}+\not{E}$ from
the two component dark sector signal (top) and $WW$ and $ZZ$ background
(bottom).

The final state that we are interested in is $\ell^{+}\ell^{-}+\not{E}$. The
following production and subsequent decay chains involving new physics
contribute dominantly to this final state.

$e^{+}e^{-}\rightarrow\phi^{+}\phi^{-},\phi^{+}\rightarrow
W^{+(*)}\phi^{0},W^{+(*)}\rightarrow\ell^{+}\nu_{\ell},\phi^{-}\rightarrow
W^{-(*)}\phi^{0},W^{-(*)}\rightarrow\ell^{-}\bar{\nu}_{\ell}$

$e^{+}e^{-}\rightarrow\psi^{+}\psi^{-},\psi^{+}\rightarrow
W^{+(*)}\psi_{1},W^{+(*)}\rightarrow\ell^{+}\nu_{\ell},\psi^{-}\rightarrow
W^{-(*)}\psi_{1},W^{-(*)}\rightarrow\ell^{-}\bar{\nu}_{\ell}$

The SM background processes that will lead to same final state are as follows.

$e^{+}e^{-}\rightarrow
W^{+}W^{-},W^{+}\rightarrow\ell^{+}\ell_{\nu},W^{-}\rightarrow\ell^{-}\bar{\nu_{\ell}}$

$e^{+}e^{-}\rightarrow
W^{+}W^{-}Z,W^{+}\rightarrow\ell^{+}{\nu_{\ell}},W^{-}\rightarrow\ell^{-}\bar{\nu_{\ell}},Z\rightarrow\nu_{\ell}\bar{\nu}_{\ell}$

$e^{+}e^{-}\rightarrow
ZZ,Z\rightarrow\ell^{+}\ell^{-},Z\rightarrow\nu_{\ell}\bar{\nu}_{\ell}$

The first background in the aforementioned list is $e^{+}e^{-}\rightarrow
W^{+}W^{-}$, which has the largest production cross-section. The most
efficient way to reduce this background is the use of beam polarization, which
we will discuss in the next subsection. $e^{+}e^{-}\rightarrow ZZ$ constitutes
the second most copious background, while $e^{+}e^{-}\rightarrow W^{+}W^{-}Z$
has much lesser production cross-section and therefore will not be considered
in the rest of the analysis. In Fig. 11, we show all the major Feynman
diagrams contributing to the final state $\ell^{+}\ell^{-}+\not{E}$ from dark
sector signal and the dominant SM background processes.

###  5.1 Effect of beam polarization

From Fig. 11, we see that $e^{+}e^{-}\rightarrow W^{+}W^{-}$ receives
contribution mostly via $t$-channel neutrino exchange, which happens due to
interactions with left-handed electrons and right-handed positrons. Such
background contribution can thus be suppressed in a linear collider with the
$e^{-}$ beam being maximally right polarized and $e^{+}$ beam to be maximally
left polarised. At the beginning of ILC run, the electron polarization is
proposed to be maximum 80% Behnke:2013lya  whereas, the positron polarization
is expected to be 30% Behnke:2013lya  with a possible upgradation to 60%
Behnke:2013lya . In Fig. 12 and 13, we will compare $\not{E}$ distributions
for BP1 (for 1TeV run) and BP3 (for 500 GeV run) for unpolarized beams and
various degrees of beam polarization and see the effect of polarization in
reducing the dominant $WW$ background. We will therefore use polarized beams
with an 80% right-polarized electron beam ($P_{e^{-}}=+0.8$) and a 30% left-
polarized positron beam ($P_{e^{+}}=-0.3$) for our analysis. We note further
that the sub-dominant background $e^{+}e^{-}\rightarrow ZZ$, gets major
contribution from lepton mediated $t$-channel diagrams (see Fig. 11). It
receives little suppression for the beam polarization of the above type. The
cross-sections from DM signal are mentioned in Tables 2 and 3, while those
from SM background is mentioned in captions of Tables 2, 3, Fig. 12 and Fig.
13. Importantly, the signal cross-section from dark sector stemming from
vector like fermion and scalar do not diminish with beam polarisation.

Figure 12: Unnormalized $\ell^{+}\ell^{-}$ event distribution of $\not{E}$ at
ILC with $\sqrt{s}=$1 TeV and 1000 $fb^{-1}$ luminosity for BP1 with (a)
$P_{e^{-}}=-0.8,P_{e^{+}}=0.3$, (b) $P_{e^{-}}=0.0,P_{e^{+}}=0.0$ and (c)
$P_{e^{-}}=0.8,P_{e^{+}}=-0.3$. Note here that $\ell^{+}\ell^{-}$ cross-
section from $W^{+}W^{-}$ production for case (a) left-handed electron and
right-handed positron corresponds to $\sigma=296$ fb, for (b) Unpolarized case
corresponds to $\sigma=128$ fb and for case (c) right-handed electron and
left-handed positron corresponds to $\sigma=18.3$ fb. $\ell^{+}\ell^{-}$
cross-section from $Z(\ell^{+}\ell^{-})Z(\nu\bar{\nu})$ for (a) left-handed
electron and right-handed positron corresponds to $\sigma=7.5$ fb, for (b)
Unpolarized case corresponds to $\sigma=4.4$ fb and for (c) right-handed
electron and left-handed positron corresponds to $\sigma=3.5$ fb.

Figure 13: Unnormalized $\ell^{+}\ell^{-}$ event distribution of $\not{E}$ at
ILC with $\sqrt{s}=$ 500 GeV and 1000 $fb^{-1}$ luminosity for BP3 with (a)
$P_{e^{-}}=-0.8,P_{e^{+}}=0.3$, (b) $P_{e^{-}}=0.0,P_{e^{+}}=0.0$ and (c)
$P_{e^{-}}=0.8,P_{e^{+}}=-0.3$. $\ell^{+}\ell^{-}$ cross-section from
$W^{+}W^{-}$ production for (a) left-handed electron and right-handed positron
corresponds to $\sigma=802$ fb, for case (b) Unpolarized beam corresponds to
$\sigma=342$ fb and for (c) right-handed electron and left-handed positron
corresponds to $\sigma=51fb$. $\ell^{+}\ell^{-}$ cross-section from
$Z(\ell^{+}\ell^{-})Z(\nu\bar{\nu})$ for (a) left-handed electron and right-
handed positron corresponds to $\sigma=21$ fb, (b) Unpolarized case
corresponds to $\sigma=12.0$ fb and for case (c) right-handed electron and
left-handed positron corresponds to $\sigma=9.6$ fb.

##  6 Analysis with selected benchmark points

###  6.1 Analysis at $\sqrt{s}=$ 1 TeV

Figure 14: Unnormalized distribution of (a) missing energy, (b) energy of the
leading lepton, (c) missing energy after applying the cut on the energy of
leading lepton $<$ 150 GeV for DM signal and SM backgrounds, (d) missing
energy distribution of signal from BP1 plus background after applying the cut
at ILC with $\sqrt{s}=$ 1 TeV, and $P_{e^{-}}=0.8,P_{e^{+}}=-0.3$.

Figure 15: Same as Fig. 14 but for BP2.

We now analyse multi-component DM signal. We take up first the benchmark
points BP1 and BP2, which, due to their higher masses, can only be probed at
$\sqrt{s}=1$ TeV. In Fig. 14 (top left), we show the corresponding $\not{E}$
distribution for signal (BP1) and dominant backgrounds at 1000$fb^{-1}$ of
integrated luminosity. As one can clearly see, the missing energy distribution
for the signal peaks at a higher value compared to that for the backgrounds.
It is observed that, in case of $W^{+}W^{-}$ background, the leptons carry
most of the energy and therefore, the missing energy peaks at a lower value.
In case of $ZZ$ background, the energy is shared comparably between the two
$Z$ bosons. Therefore $\not{E}$ peaks at $\frac{\sqrt{s}}{2}$ i.e. at 500 GeV
in this case. The signal distribution on the other hand, has two peaks. The
peak at lower missing energy corresponds to the fermionic DM and the higher
one arises due to the scalar DM, as the mass of the fermionic DM is smaller
than the scalar DM mass at this benchmark. Its worthy to remind here that the
choice of masses also take care DM constraints, for example, $m_{\psi_{1}}\sim
m_{h}/2$ provides adequate depletion via Higgs resonance in satisfying the
relic density limits, as discussed in the previous section. The scalar sector
being almost degenerate shows much narrower peak compared to the fermionic
case, where $\Delta m_{2}$ is comparatively large and t therefore the
distribution is flatter. In Fig. 14 (top right), we plot the energy of the
leading lepton. As at ILC, the total energy of collision is fixed, the energy
distribution of the leading lepton shows a sharp complementarity with the
$\not{E}$. It peaks at a higher value for the backgrounds and at a lower value
for signal. Lepton energy distribution also retains the double-peak behaviour
like $\not{E}$ distribution. Next we apply a cut on the energy of the leading
lepton and the $\not{E}$ distribution of the $W^{+}W^{-}$ background becomes
significantly softer (see bottom left panel of Fig. 14). The first peak of the
signal also shifts to the right because of this cut, whereas the second peak
remains almost unaffected. However, the resulting position and size of the
peak of the $W^{+}W^{-}$ background distribution is subject to the lepton
energy cut applied. In this case we apply $E_{\ell_{1}}<150$ GeV, so that the
peak of the background distribution coincides with the first peak of the
signal, and the distinction between the resulting first and second peak become
even more prominent. The extent at which the value of this cut affects our
analysis will be discussed shortly. It may be noted that the contribution from
$ZZ$ background becomes negligible after applying the leading lepton energy
cut and therefore is not shown explicitly in Fig. 14 (bottom left). We further
show in bottom right panel of Fig. 14 the final signal plus background
distribution after applying the aforementioned cut. We can see that the two
peaks can be distinguished rather well, even in the presence of background
events.

In Fig. 15, we perform the same analysis for BP2. We can see that the
separation between two signal peaks is slightly worse here, see top left and
top right panel figures of Fig. 15. This is because $\Delta m_{2}$ in this
case is smaller than the previous case. However, after applying the leading
lepton energy cut, the first peak, as well as $WW$ background gets diminished
and the heights of the two peaks in the signal plus background $\not{E}$
distribution become comparable. One can note that, in this case after applying
the cut $E_{\ell_{1}}<150$ GeV, the background distribution does not exactly
coincide with the first peak of the signal (see bottom left panel figure) and
therefore, the first peak of the signal plus background distribution (bottom
right) is flatter compared to the previous case BP1 as in Fig. 14.

###  6.2 Analysis at $\sqrt{s}=$ 500 GeV

In BP3 and BP4, it is possible to produce HDSPs on-shell at $\sqrt{s}=500$
GeV. The case of BP3 is shown in Fig. 16. It is evident that the production
cross-section of the charged fermion is much larger here due to low mass.
Therefore, the size of the first peak is also bigger compared to the second
peak (see top left and top right panel figures). However, the lepton energy
cut somewhat reduces the area under the first peak. The additional $WW$
background events (see bottom left figure) nevertheless reintroduces a
difference in height of the two peaks in the resulting distribution (see
bottom right).

The case for BP4 is shown in Fig. 17. The charged fermion mass is larger and
consequently the size of the fermionic peak is smaller (see top left and top
right figures). Therefore, after applying the lepton energy cut and finally
summing up the signal and background events, the heights of the two peaks are
comparable. This is clear from bottom left and bottom right panel figures of
Fig. 17. Based on such qualitative conclusions, we attempt a set of
distinction criteria in the following section.

Figure 16: Same as Fig. 14 but for BP3 at $\sqrt{s}=$ 500 GeV and cut on
energy of leading lepton $<$ 75 GeV applied.

Figure 17: Same as Fig. 16, but for BP4.

###  6.3 A fifth benchmark

Benchmark |  $m_{\phi^{0}}$ and $\Delta m_{1}$ |  $m_{\psi_{1}}$ and $\Delta m_{2}$ |  $\sigma$(fb) |  $\sigma$(scalar) |  $\sigma$(fermion)
---|---|---|---|---|---
BP5 | 10 GeV, 120 GeV | 455 GeV, 40.0 GeV | 1.3 fb | 0.66 fb | 0.66 fb
Table 4: Signal benchmark point and cross section for
$P_{e^{-}}=0.8,P_{e^{+}}=-0.3$ polarization at $\sqrt{s}=$ 1 TeV centre-of-
mass energy.

Figure 18: Same as Fig. 14 but for benchmark point BP5. Figure 19: Signal from
a single component DM together with SM background distribution at
$\sqrt{s}=$500 GeV and 1000 $fb^{-1}$ luminosity.

Before proceeding further let us now examine the situation when the fermionic
dark sector is degenerate instead of the scalar dark sector (Region III
mentioned above). We choose another benchmark point BP5 keeping this criterion
in mind. As the production cross-section of charged scalar is usually lower
than that of the charged fermions, we chose $m_{\phi}^{\pm}<m_{\psi^{\pm}}$.
Non-degenerate scalar and degenerate fermion sector in that case immediately
imply $m_{\phi^{0}}<m_{\psi_{1}}$. Therefore the missing energy peak
pertaining to the scalar DM will appear on the left of the fermionic peak,
will be smaller and broader due to small cross-section and large $\Delta m$.
On the contrary, the peak from fermionic dark sector will be larger and
sharper due to larger cross-section and smaller $\Delta m$. In such a case,
the peak from scalar sector will be mostly buried under $WW$ tail. We show the
resulting distributions in Fig. 18. The missing energy distribution in this
case not only fails to produce clearly separated peaks of comparable sizes,
but yields very similar distribution to the scenario when there exists only a
single DM component and the background distribution contributes to a second
peak-like behaviour. If we compare Fig. 18 (bottom right) with Fig. 19, where
only inert doublet has been used to plot the signal $\not{E}$ distribution, a
qualitative similarity is obvious. This is a rather general consequence of the
lower peak from DM signal having a much flatter distribution.

From the discussion above, one may thus conclude that Region III is by and
large more difficult from the angle of peak distinction, as compared to the
other regions in our classification. We reiterate that the distinguishability
of the two peaks in Regions I, II and IV is facilitated by our capacity to use
$\not{E}$ as a kinematic variable, something that is possible in an
$e^{+}e^{-}$ machine only.

###  6.4 Impact of lepton energy cuts

As pointed out earlier, at the $e^{+}e^{-}$ collider, the lepton energy
distribution shows a complementarity with the missing energy distribution. We
have applied a cut on the energy of the leading lepton to suppress the
background and we also see that putting such cut also helps us retain the
Gaussian nature of the resulting distribution. However, the choice of the cut
is crucial in this case. As the total energy of collision in each event is
fixed, the higher lepton energy naturally correspond to smaller $\not{E}$.
Therefore, an upper cut on the lepton energy immediately indicate a lower cut-
off in $\not{E}$ distribution. In case of $WW$ background $\not{E}$
distribution peaks at a much lower value compared to the DM signal. Hence an
upper cut on lepton energy distribution eliminates the low $\not{E}$ region. A
harder lepton energy cut pushes the $\not{E}$ distribution of the $WW$
background towards larger $\not{E}$. On the other hand, the first peak of the
signal is also somewhat affected by this cut. But, the second peak remains
unaffected as long as the lepton energy cut is not too hard. Our goal
therefore is to shift the $WW$ distribution to the right so that it coincides
with the first peak of the signal which elucidates the double hump behaviour.

Figure 20: Signal (BP1) and SM background distribution after applying cuts
$E_{\ell}<150$ GeV (left) and $E_{\ell}<200$ GeV (right) with $\sqrt{s}=$ 1
TeV for $P_{e^{-}}=0.8,P_{e^{+}}=-0.3$ at ILC.

In Fig. 20, we compare $\not{E}$ distribution after applying a lepton energy
cut of 150 GeV (left), and 200 GeV (right). We can see that the background
distribution is pushed towards right with harder lepton energy cut. Thus there
is a need to optimise the lepton energy cut depending on the DM particle
masses that are being probed.

##  7 Distinction criteria for two peaks

As we argue here that two different dark sectors having mass differences
$\Delta m_{1}$ and $\Delta m_{2}$ between the HDSP and DM giving same final
state, will in general have missing energy peaks at different $\not{E}$, which
will partially overlap when the full signal produced from both dark sector is
analysed. It is understood that when the difference between $\Delta m_{1}$ and
$\Delta m_{2}$ is significant, it will show a distorted missing energy
distribution where two peaks can be distinctly identified. The missing energy
distribution (of the signal with reduced SM background) for this two component
scenario looks like a two-peak asymmetric Gaussian distribution as below, to
which the simulated $\not{E}$ distributions can be fitted.

$\displaystyle G(x)$ $\displaystyle=$ $\displaystyle G_{1}(x)+G_{2}(x)+B$ (25)
$\displaystyle=$ $\displaystyle
A_{1}~{}e^{-\frac{(x-\mu_{1})^{2}}{2\sigma_{1}^{2}}}+A_{2}~{}e^{-\frac{(x-\mu_{2})^{2}}{2\sigma_{2}^{2}}}+B~{}~{}.$

Here $G_{1}(x)$ is the Gaussian function with amplitude $A_{1}$, mean
$\mu_{1}$, and standard deviation $\sigma_{1}$ which correspond to the first
peak, while $G_{2}(x)$ is the function containing the corresponding quantities
for the second DM peak. The constant (or slowly varying) parameter $B$ is
further introduced, in order to account for various theoretical uncertainties
as well as those due to $\not{E}$ mismeasurement etc 101010In the results
presented below, we have treated $B$ as a constant function of $\not{E}$ for
simplicity. However, we have checked that the nature of the peaks do not
change appreciably on fitting $B$ as a polynomial function upto second degree.
Differences, if any, are noticeable mainly in regions away from both peaks.. A
schematic description of such fitting is shown in Fig. 21.

Figure 21:  A cartoon of a missing energy distribution with two peaks.

In Fig. 21, we can identify two peaks at $\not{E}$ values $t_{1}(=\mu_{1})$
and $t_{2}(=\mu_{2})$. The number of events at those peaks are roughly denoted
by $y_{1}=y(t_{1})$ and $y_{2}(t_{2})$ respectively; actually this will be the
area under the curve in a small interval around the peak as indicated in the
figure. The minima between two peaks is identified as $t_{\rm min}$ and the
corresponding value along y-axis is denoted by $y_{\rm min}$.

We are now ready to set up the criteria for distinguishing the two peaks in
missing energy spectra. Each such criterion must address either of both of the
following questions:

(a) How to elicit the prominence of the second (read smaller) peak relative to
the first (bigger) peak?
(b) How to resolve best the separation between the two peaks?

  * •

C1: The first condition examines how much the second peak distorts the
symmetry of the distribution about the first peak. This will require us to
compare the number of events within n$\sigma$ ($n\geq 1$) range of the first
peak on both sides. Let us assume, $t$ (on left) and $t^{\prime}$ (on right)
are the two positions which are n$\sigma$ ($n\geq 1$) away from the first peak
$t_{1}$, respectively, as in Figure 21. Let us also suppose,

$\displaystyle\Delta N_{1}=\int_{t}^{t_{1}}ydt,~{}~{}\Delta
N_{2}=\int_{t_{1}}^{t^{{}^{\prime}}}ydt$ (26)

Then, if

$\displaystyle\frac{|\Delta N_{2}-\Delta N_{1}|}{\sqrt{\Delta N_{1}}}>2\,,$
(27)

we can stipulate that the peaks are resolved to 2$\sigma$ significance. Note
that the denominator in Equation 27 represents the fluctuation of the
distribution corresponding to the first DM particle. The criterion C1, defined
in terms of integrated quantities, is useful when the two peaks are not
visually prominent but a distortion to the $\not{E}$ spectrum occurs. It is
also worth registering that this criterion has its limitation when the second
peak is much smaller than the first.

  * •

C2: The second condition addresses how high does the second peak go. Let us
assume $y_{3}=y(t^{{}^{\prime\prime}})$ is the number of events at
$t^{{}^{\prime\prime}}$ which is assumed to be $k\sigma$ ($k\geq 1$) away from
the first peak $t_{1}$ on the right side and
$y_{3}^{\prime}=y^{\prime}(t^{{}^{\prime\prime}})$ is the corresponding number
of events in the absence of second peak upon a Gaussian fit. Then if the
difference between these two numbers follows the relation given below:

$\displaystyle y(t^{{}^{\prime\prime}})-y^{\prime}(t^{{}^{\prime\prime}})$
$\displaystyle>$ $\displaystyle 2\sqrt{y^{\prime}(t^{{}^{\prime\prime}})},$
$\displaystyle{\rm or~{}}~{}y_{3}-y_{3}^{{}^{\prime}}$ $\displaystyle>$
$\displaystyle 2\sqrt{y_{3}^{{}^{\prime}}}\,.$ (28)

If C2 is satisfied, then we will be able to say that the fluctuation is more
than 2$\sigma$ and can be termed as a second peak. Note that C2 offers a
criterion at the differential level, in contrast with C1, which comes at an
integral level.

  * •

C3: Another possibility to check how significant the second peak is to compare
the number of events within a close window ($\Delta t\sim$ 5 GeV) around two
peaks. One may take a ratio between the difference and sum of the following
quantities as:

$\displaystyle\epsilon$ $\displaystyle=$
$\displaystyle\frac{\int_{t_{1}-\Delta t}^{t_{1}+\Delta
t}ydt-\int_{t_{2}-\Delta t}^{t_{2}+\Delta t}ydt}{\int_{t_{1}-\Delta
t}^{t_{1}+\Delta t}ydt+\int_{t_{2}-\Delta t}^{t_{2}+\Delta
t}ydt}~{}~{}\underrightarrow{\\{\Delta t\to
0\\}}~{}~{}\frac{y_{1}-y_{2}}{y_{1}+y_{2}}.$ (29)

The smaller is $\epsilon$, the more significant is the second peak. This
essentially compares the number of events about the two maxima, when the peaks
are comparably high.

  * •

C4: This condition addresses how significant is the second peak with respect
to the minimum in between. Obviously if there exists a second peak (maxima)
along with the first, there should be a local minimum between them
($t_{1}<t_{\rm min}<t_{2}$). Then the significance of the second peak with
respect to the local minimum can be obtained as follows:

$\displaystyle\frac{y(t_{2})-y(t_{\rm min})}{\sqrt{y(t_{\rm
min})}}>2\,;~{}~{}~{}~{}~{}~{}~{}{\rm
or~{}equivalently}~{}~{}~{}~{}~{}~{}\frac{y_{2}-y_{\rm min}}{\sqrt{y(t_{\rm
min})}}>2\,.$ (30)

This says that the second peak rises more than 2$\sigma$ with respect to the
fluctuation of the intermediate minimum. This too is a criterion at the
differential level.

It is obvious from the discussion above that the criteria to distinguish two
peaks significantly depend on statistics/luminosity and the parameters $n(k)$.
We will explore the effect of these factors explicitly for our chosen BP’s, in
the context of collider signature of two component DM scenario.

BP1

Figure 22: Gaussian Fitting of $\not{E}$ distribution for signal (BP1) +
background with different integrated luminosities. The polarization of initial
beams $P_{e^{-}}=0.8,P_{e^{+}}=-0.3$ has been used.

Figure 23: Dependence of selection criteria (top left) C1 on $\int{\cal L}dt$
for different $n$ (defined in Equation 27), (top right) C2 on $\int{\cal L}dt$
for different $k$ (defined in Equation 28), (bottom left) C3 on $\Delta t$
(defined in Equation 29) and (bottom right) C4 on $\int{\cal L}dt$ for BP1.

BP2

Figure 24: Fitting of $\not{E}$ distribution for signal (BP2) + background
with different integrated luminosities. The polarization of initial beams
$P_{e^{-}}=0.8,P_{e^{+}}=-0.3$ has been used.

Figure 25: Dependence of selection criteria (top left) C1 on $\int{\cal L}dt$
for different $n$ (defined in Equation 27), (top right) C2 on $\int{\cal L}dt$
for different $k$ (defined in Equation 28), (bottom left) C3 on $\Delta t$
(defined in Equation 29) and (bottom right) C4 on $\int{\cal L}dt$ for BP2.

BP3

Figure 26:  Fitting of $\not{E}$ distribution for signal (BP3) + background
with different integrated luminosities. The polarization of initial beams
$P_{e^{-}}=0.8,P_{e^{+}}=-0.3$ has been used.

Figure 27: Dependence of selection criteria (top left) C1 on $\int{\cal L}dt$
for different $n$ (defined in Equation 27), (top right) C2 on $\int{\cal L}dt$
for different $k$ (defined in Equation 28), (bottom left) C3 on $\Delta t$
(defined in Equation 29) and (bottom right) C4 on $\int{\cal L}dt$ for BP3.

BP4

Figure 28:  Fitting of $\not{E}$ distribution for signal (BP4) + background
with different integrated luminosities. The polarization of initial beams
$P_{e^{-}}=0.8,P_{e^{+}}=-0.3$ has been used.

Figure 29: Dependence of selection criteria (top left) C1 on $\int{\cal L}dt$
for different $n$ (defined in Equation 27), (top right) C2 on $\int{\cal L}dt$
for different $k$ (defined in Equation 28), (bottom left) C3 on $\Delta t$
(defined in Equation 29) and (bottom right) C4 on $\int{\cal L}dt$ for BP4.

In Figs. 22, 24, 26 and 28, we show the resulting (signal plus reduced
backgrounds) $\not{E}$ distribution for benchmark points BP1, BP2, BP3 and BP4
respectively for various benchmark luminosities. In each plot we have
mentioned all the required parameters namely, $\mu_{1,2}$, $\sigma_{1,2}$ and
$A_{1,2}$ for the Gaussian fits. One can see that at low luminosities, meagre
statistics leads to large fluctuations in the distributions. Therefore, to
achieve smooth distributions, and thereby pin-down the actual physics effects,
we have to analyse events at high luminosities (upto 10 $ab^{-1}$). Such high
luminosities are in fact unrealistic in the context of the proposed future
$e^{+}e^{-}$ collider. However, the distributions with reasonable design
luminosities (100, 500, 1000 $fb^{-1}$) are also presented for the sake of
realistic analysis.

Let us now describe each cases separately. In Fig. 23, we analyse BP1 in
detail. Fig. 23 (top left), shows the validation of C1 as a function of
$\int{\cal L}dt$ for different values of $n$. One can see from the figure that
C1 increases with $\int{\cal L}dt$ (evident from Equation 27) as well as $n$.
In top right panel of Figure 23, C2 is shown as a function of luminosity for
various $k$ values. C2 increases with $\int{\cal L}dt$ following Equation 28.
For a specific luminosity it increases with $k$ as well. In bottom left of
Fig. 23, we show the variation of C3 as a function $\Delta t$ for three
representative luminosities. One can see that C3 remains almost constant for a
small range $\Delta t$, but starts increasing with $\Delta t$ after a point.
Looking at the definition of C3 in Equation 29, it is clear that C3 should
remain constant with increasing $\int{\cal L}dt$, which is also seen from
bottom left of Fig. 23. The small deviation that takes place is due to
statistical fluctuation. Bottom right panel figure in Fig. 23, presents
dependence of C4 on $\int{\cal L}dt$. One can see that C4 also increases with
increasing luminosity, as was clear from the definition in Equation 30. The
organisation of Fig. 25, 27 and 29, which correspond to BP2, BP3 and BP4
respectively, remains the same as Fig. 23.

A comparison between our chosen benchmarks in the light of the aforementioned
distinction criteria is in order. One can see that in case of BP3, the
amplitude of the second peak is much smaller compared to that of the first.
Therefore, in this case the distortion of the first peak due to the second is
minimal. Consequently, C1 and C2 criteria take the lowest value for a specific
$n$ and ${\cal L}dt$ in this case, among all the benchmarks. C1 takes the
largest value for BP2, since the relative height of the second peak w.r.t the
first peak as well as the width of the second peak is large in this case. On
the other hand, C2 is highest for BP4, since the height of the second peak is
largest there. C3 is largest for BP3, due to significant asymmetry in the
heights of the two peaks. BP1 takes the smallest value under the C3 criterion.
C4 again, will take the maximum value for BP4, due to large height and small
width of the second peak. It is seen that BP3 actually performs poorly under
all the conditions. BP1 performs best under C3, BP2 under C1 and BP4 under C2
as well as C4 criteria.

Figure 30: Fitting of $\not{E}$ distribution for signal (BP1) + background for
$P_{e^{-}}=0.4,P_{e^{+}}=-0.6$ (top left), $P_{e^{-}}=0.8,P_{e^{+}}=-0.3$ (top
right) and $P_{e^{-}}=0.8,P_{e^{+}}=-0.6$ (bottom centre). Figure 31: Left: C1
as a function of $\int{\cal L}dt$ for different $n$ values and Right: contour
plot in the plane spanned by $n$ and $\int{\cal L}dt$ for various fixed C1
values for various polarization configurations.  Figure 32: Left: C2 as a
function of $\int{\cal L}dt$ for different $k$ values and Right: contour plot
in the plane spanned by $k$ and $\int{\cal L}dt$ for various fixed C2 values
for various polarization configurations.  Figure 33: Left: C3 as function of
$\Delta t$ and Right: C4 as a function of $\int{\cal L}dt$ for various
polarization configurations.

After the comparison among our chosen benchmarks is done, we move on to see
the effect of polarization in two-peak identification for a specific benchmark
(BP1). To omit statistical fluctuations and achieve smooth Gaussians which
will encapsulate only the role of polarization, we use an unrealistically
large $\int{\cal L}dt$ value. However, we reiterate that, such luminosity is
chosen with the sole purpose of brining out the polarization effects without
an interference of statistical fluctuations. In Fig. 30, we show the Gaussian
fits of the $\not{E}$ distribution for BP1 for three different polarization
choices. It should be noted that initial state polarization largely affects
the background distribution as discussed earlier. The signal distribution
remains mostly constant. Furthermore, after applying the lepton energy cut the
background distribution alters the first peak significantly, while the second
peak remains unaffected. We have checked that the background cross-section is
largest for polarization configuration P1
$\equiv(P_{e^{-}}=0.4,P_{e^{+}}=-0.6$) followed by P2
$\equiv(P_{e^{-}}=0.8,P_{e^{+}}=-0.3$) and P3
$\equiv(P_{e^{-}}=0.8,P_{e^{+}}=-0.6$). Consequently, the relative size of the
second peak w.r.t the first peak follows the reverse order, which is evident
from the Fig. 30.

In the left panel of Fig. 31, we show the dependence of C1 as function of
$\int{\cal L}dt$, for different $n$ for all the polarization configurations.
One can see from the figure that all three polarization configurations yield
very similar results for all $n$ and $\int{\cal L}dt$ values. In the right
panel of Fig 31, we present contour plot in the plane spanned by $\int{\cal
L}dt$ and $n$ for fixed C1 ($=2$) value. Here we see that with decrease in
$\int{\cal L}dt$ one needs to choose larger $n$ values, in order to achieve C1
$>$ [-0.07cm] $\sim$ 2\. We mention here that we have taken
$n<\frac{\mu_{2}-\mu_{1}}{\sigma_{1}}$, since we are interested in scanning
for $n$ values starting from the first peak upto the second peak.

In the left panel of Fig. 32, we show the variation of the quantity C2 with
$\int{\cal L}dt$ for specified $k$ values, for the aforementioned polarization
configurations (P1, P2 and P3). Here too, one can see that the polarization
affects the criterion only marginally. A contour plot in a plane spanned by
$\int{\cal L}dt$ and $k$ for P1, P2 and P3 is depicted in the right panel of
Fig. 32. Larger $k$ values will be necessary to achieve C2
$\raisebox{-3.69899pt}{~{}\shortstack{$>$ \\\ [-0.07cm] $\sim$}}~{}2$ at lower
luminosities. Similar to the C1 case, here we have taken
$k<\frac{\mu_{2}-\mu_{1}}{\sigma_{1}}$ since we have only considered $k$
values in the range $\mu_{2}-\mu_{1}$.

In the left panel of Fig. 33, we show the dependence of C3 on $\Delta t$ for
various polarizations. We know the height difference between the first and
second peak is largest in case of P1 and smallest in P3. Hence, following
Equation 29, the quantity C3 decreases from P1 to P3. This behaviour is clear
from the figure. In the right panel of Fig. 33, C4 is plotted for all the
aforementioned polarizations as a function $\int{\cal L}dt$. Since the second
peak in our case is unaffected by background addition but the first peak is
altered, only the minima will change from P1 to P3 case (see Equation 30).
Clearly, C4 value for P1 will be smallest followed by P2 and then P3.
Therefore, after the discussion of distinction criteria in the context of
initial beam polarization, we conclude that the polarization effects although
necessary to reduce the background events, are not very significant for
distinction between two peaks.

Figure 34: Fitting of $\not{E}$ distribution for signal (BP1) + background
with different lepton energy cut. The polarization of initial beams
$P_{e^{-}}=0.8,P_{e^{+}}=-0.3$ and $\int{\cal L}dt=20ab^{-1}$ has been used.

Next we analyse the effect of lepton energy cut on the distinction criteria.
In Fig. 34, we show the distribution of $\not{E}$ for signal (BP1) +
background with a lepton energy cut (top left) 100 GeV, (top right) 150 GeV
and (bottom centre) 200 GeV. In case of 100 GeV cut, the background gets
reduced to a large extent. However, significant portion of the first peak of
the signal also gets rejected by the stringent lepton energy upper cut off.
Consequently, the two-peak nature of the distribution disappears and only a
small bump in the distribution is observed. We mention that in this case too,
C1 and C2 criteria can be analysed in principle, to quantify the impact of the
second DM candidate on the distribuiton of the first one. However, we do not
carry out this analysis, since our focus is on the distinctness of two peaks
in the distribution. With 200 GeV cut on the leading lepton energy, the
reduction of background events is less compared to the 150 GeV case. However,
both of these cases show clear two-peak signature. Let us now quantify the
distinguisability of the peaks for various lepton energy cuts.

Figure 35: (top left) C1 as a function of $\int{\cal L}dt$ for different $n$
values, (top right) a contour plot in $\int{\cal L}dt-n$ plane, (centre left)
C2 as a function of $\int{\cal L}dt$ for different $k$ values, (centre right)
a contour plot in $\int{\cal L}dt-k$ plane, (bottom left) C3 as function of
$\Delta t$ and (bottom right) C4 as a function of $\int{\cal L}dt$. All the
plots are given for two specified lepton energy cuts.

In Fig. 35 (top left), we plot C1 as a function of $\int{\cal L}dt$ for
various $n$ values with 150 GeV and 200 GeV lepton energy cuts. In Fig. 35
(top right), we draw the corresponding contours of C1 = 2 in $\int{\cal
L}dt-n$ plane, with the aforementioned cuts. Fig. 35 (centre left and right)
show the analogous distributions for C2. Next, C3 as a function of $\Delta t$
for the specified cuts is presented in Fig. 35 (bottom left). Lastly, C4 as a
function of $\int{\cal L}dt$ for both the cut values is shown in Fig. 35
(bottom right). One can see that at a given luminosity one has to choose
higher $n(k)$ in order to get same C1(C2) values with 200 GeV cut as compared
to 150 GeV cut. The height of the first peak becomes bigger in case of 200 GeV
cut compared to 150 GeV cut. Therefore C3 for the former is larger compared to
the latter, see bottom left of Fig. 35. C4 values remain almost identical in
both cases, as can be seen from Fig. 35 (bottom right).

Let us finally summarize the key findings of this section.

  * •

At low luminosity, large statistical fluctuation becomes a roadblock while one
tries to identify the two-peak signature in the $\not{E}$ distribution in two-
component DM scenario.

  * •

C1, C2 (for specific $n(k)$) and C4 increase with integrated luminosity
$\int{\cal L}dt$.

  * •

With increase in luminosity, significant C1 and C2 values can be obtained at
lower $n(k)$ values.

  * •

C3 increases with $\Delta t$ and remains almost constant with integrated
luminosity.

  * •

The distinction criteria are not too sensitive to the initial state
polarization, given significant background reduction is already achieved.

  * •

The distinction criteria are sensitive to the cut applied on the leptons in
the final state. If the cut is too stringent, the two-peak nature is lost. On
the other hand, if the cut is too relaxed, the second peak will become
insignificant compared to the first. Therefore, one needs to choose the lepton
energy cut optimally.

##  8 Summary and conclusion

In this analysis we illustrate how two component WIMP models can manifest them
in collider signal. In particular, we study the case of extended dark sector
which produces HDSPs at collider with subsequent decay chains to produce DM.
The DM can be identified via kinematic variables like missing energy
($\not{E}$) or missing transverse momentum ($\not{E}_{T}$) or missing mass
($\not{M}$) distribution, where the seed of two dark sectors coexist. We
demonstrate that distributions in above kinematic variables can show a double
peak behaviour where the distinct peaks may arise from different dark sectors.
We show that the key variables that play a role are both the DM masses
($m_{\rm DM_{1}},m_{\rm DM_{2}}$) and their splitting with the corresponding
HDSPs ($\Delta m_{1},\Delta m_{2}$), whenever they are segregated from each
other the double peak may show up in the distribution.

However, not all the kinematic variables are equivalent. For example, DM mass
does not play a role in variables like $\not{E}_{T}$, while it does in
$\not{E}$ or in $\not{M}$. So, $\not{E}_{T}$ has limitations in exploring the
two component dark sector compared to $\not{E}$. However, $\not{E}$ or
$\not{M}$ can’t be measured in context of hadron collider or LHC, absent the
knowledge of partonic CM energy. We therefore choose to illustrate the
phenomena in context of $e^{-}e^{+}$ collider. We will discuss the
possibilities of exploring the same in context of LHC in a subsequent
publication.

When it comes to the kind of model that can be tested under this class of
signals like $\ell^{+}\ell^{-}$ as explained here, we argue that this is
possible when each of the dark sector has some transformation under the SM
weak gauge group, like an isodoublet or an isotriplet etc. This will provide
some HDSPs particles with EM charge that can have larger production cross-
section with subsequent decay to DM plus SM leptons or jets. The model we
choose to illustrate here inherit an inert scalar doublet stabilised under a
$\mathcal{Z}_{2}$ symmetry providing a scalar DM, one vector like fermion
doublet and a right handed fermion singlet both stabilised under a different
$\mathcal{Z}^{{}^{\prime}}_{2}$ providing a fermion DM. This exercise can
easily be done containing both fermion DMs or both scalar ones, just the
parameter space where such observations can be made turns different after
satisfying DM and other constraints. We briefly mention that a distinction of
two DM components in $\not{E}$ distribution for mono-X signals which stems
from non-extended dark sectors turn much harder. We will also elaborate on it
in a coming publication.

When it comes to the presence of distinct peaks from two different dark
sectors, the role of relative signal cross-sections from both dark sectors as
well as that from the SM background turn immensely important. Regarding
signal, we adjudge that unless both the dark sectors provide similar
production cross-section, the peak of one may get buried under the other even
if segregated. In this analysis, we show further that $WW$ background to
$\ell^{+}\ell^{-}$ signal can both spoil or elucidate the double hump
behaviour of the signal $\not{E}$ distribution. A judicious cut on the lepton
energy may be of help to keep $WW$ peak coincide to one of the DM signal.
Apart, we use right handed electron polarisation and left handed positron
polarisation that is maximally accessible at ILC to reduce the $WW$ background
contamination, while retaining the signal as they are not chiral in nature.

Finally distinguishing the peaks from the kinematic distribution is a
challenging task and involves a careful introspection. To this end we provide
a set of criteria, which can be checked for seeing how the second peak is
separated from the first one, how sharp and large it is and at what luminosity
they can be segregated in statistically significant way, again keeping the SM
background in consideration. We carefully verify all such conditions for the
chosen benchmark points and show that the distinction can indeed be achieved
at a large luminosity, while at low luminosity, statistical fluctuation
becomes an issue. Further, the distinguishability criteria is seen to be
insensitive to beam polarisation, after the SM background reduction but they
depend on different lepton energy cuts.

## Acknowledgments

SB and JL would like to acknowledge DST-SERB grant CRG/2019/004078 from Govt.
of India. PG would like to acknowledge the support from DAE, India for the
Regional Centre for Accelerator based Particle Physics (RECAPP), Harish
Chandra Research Institute.

## Appendix

##  Appendix A $\not{E_{T}}$ and $\not{E}$ in single-component DM

We consider pair production of two heavy intermediate particles(HDSPs), which
will be produced at rest. Under this assumption, the lab frame will be same as
rest frame of the intermediate heavy particles. Let us also assume, for
simplicity, each of these pair-produced HDSPs with mass $m_{HDSP}$, decays to
massless SM fermion and DM. After using four-momenta conservation the
following relations hold.

When the HDSP is at rest, 3-momentum conservation on both ends of the decay
chain implies

$|p_{DM}|=|p_{f}|$ (31)

and energy conservation at both ends of the decay chain implies

$E_{DM}+E_{f}=m_{HDSP}$ (32)

Now using equation of motion for DM (true for both DMs at either end of the
decay chain), we get

$\sqrt{|p_{DM}|^{2}+m_{DM}^{2}}+|p_{f}|=m_{HDSP}$ (33)

Substituting $|p_{f}|$ in terms of $|p_{DM}|$, from Equation 31, we get

$\displaystyle\sqrt{|p_{DM}|^{2}+m_{DM}^{2}}+|p_{DM}|$ $\displaystyle=$
$\displaystyle m_{HDSP}$ $\displaystyle|p_{DM}|^{2}+m_{DM}^{2}$
$\displaystyle=$ $\displaystyle(m_{HDSP}-|p_{DM}|)^{2}$ $\displaystyle
m_{HDSP}^{2}-m_{DM}^{2}$ $\displaystyle=$ $\displaystyle 2m_{HDSP}|p_{DM}|$
$\displaystyle|p_{DM}|$ $\displaystyle=$
$\displaystyle\frac{(m_{HDSP}^{2}-m_{DM}^{2})}{2m_{HDSP}}$
$\displaystyle|p_{DM}|$ $\displaystyle=$ $\displaystyle\frac{1}{2}\Delta
m\left(1+\frac{m_{DM}}{m_{HDSP}}\right)$ $\displaystyle|p_{DM}|$
$\displaystyle=$ $\displaystyle\frac{1}{2}\Delta
m\left(1+\frac{m_{DM}}{m_{DM}+\Delta m}\right)$ (34)

Now this equation is true for both DMs at the end of the decay chain. However,
let us now denote the magnitude of their momenta as $|p^{1}_{DM}|$ and
$|p^{2}_{DM}|$ and the angle between their momentum three-vector as $\theta$.
Therefore, the missing momenta of the total system ($\not{E_{T}}$) will be
given as

$\not{E_{T}}=\sqrt{|p^{1}_{DM}|^{2}+|p^{2}_{DM}|^{2}+|p^{1}_{DM}||p^{2}_{DM}|\cos\theta}$
(35)

It is evident that the maximum of this quantity will occur when $\theta=0$ ie.
the two DM particles are moving along the same direction. Therefore,

$\not{E_{T}}_{max}=|p^{1}_{DM}|+|p^{2}_{DM}|=\Delta
m\left(1+\frac{m_{DM}}{m_{DM}+\Delta m}\right)$ (36)

On the other hand the energies of the two DM particles will be

$\displaystyle E^{1}_{DM}=\sqrt{|p^{1}_{DM}|^{2}+m_{DM}^{2}}$ (37)
$\displaystyle E^{2}_{DM}=\sqrt{|p^{2}_{DM}|^{2}+m_{DM}^{2}}$ (38)

The total contribution two the missing energy therefore will be

$\not{E}=E^{1}_{DM}+E^{2}_{DM}=\sqrt{|p^{1}_{DM}|^{2}+m_{DM}^{2}}+\sqrt{|p^{2}_{DM}|^{2}+m_{DM}^{2}}$
(39)

##  Appendix B Comparison between $\not{E}$ and $\not{M}$

We present a detailed comparison between the two observables $\not{E}$ and
$\not{M}$(missing mass). For the sake of simplicity we restrict our analysis
in the context of single-component DM. The definition of $\not{M}$ is given as
follows.

$\not{M}=\sqrt{\not{E}^{2}-|{\sum_{i}}\vec{p}_{\ell_{i}}|^{2}};$ (40)

It is clear from the definition of missing mass that it is it can be similar
or different compared to $\not{E}$ depending on the momenta of the leptons and
the relative orientaion of them. Keeping this in mind, we consider a few
scenarios below and explore various possibilities.

  * •

${\Delta m}=10$ GeV and $m_{DM}=50$ GeV

In this case, HDSP mass low ($\approx 60$ GeV), therefore the HDSP will be
produced with significant boost. Therefore, the lepton and DM will be almost
collinear. It will also imply that the two leptons from the two HDSP’s will be
produced almost back-to-back for conservation of four-momenta. The effective
visible momenta (the second term in Equation 40 will be negligibly small. In
this case, $\not{M}$ and $\not{E}$ will be almost overlapping on each other,
which is clear from Figure 36.

Figure 36: Comaparison between $\not{E}$ and $\not{M}$ distribution for
${\Delta m}=10$ GeV and $m_{DM}=50$ GeV

  * •

${\Delta m}=100$ GeV and $m_{DM}=390$ GeV

In this case, HDSP mass is heavy ($\approx 490$ GeV $\approx\sqrt{S}/2$),
therefore will be produced at rest. Hence, its decay products will necessarily
cancel their net momenta among each other. Consequently, the two leptons from
the two HDSP’s produced need not be back to back. However, in this case the
magnitude of lepton momenta itself will be negligible(since the DM is really
massive). Hence, here too, the two observables, namely $\not{E}$ and $\not{M}$
will overlap with each other (see Figure 37).

Figure 37: Comaparison between $\not{E}$ and $\not{M}$ distribution for
${\Delta m}=100$ GeV and $m_{DM}=390$ GeV

  * •

${\Delta m}=200$ GeV and $m_{DM}=50$ GeV

Let us further consider a possibility where the HDSP produced is of
intermediate mass (250 GeV here). We also take a DM mass not very large.
Therefore, we expect the HDSP will have some boost and consequently the
leptons will have a partial cancellation of four-momenta among themselves.
Also, the mass of the DM, being not so large enables the leptons to have non-
negligible absolute momenta. As a result the second term in Equation 40, will
be somewhat important and the two distributions will show some distinction
especially in the tail(left-end) of the distribution, since there the lepton
momenta is substantially larger compared to the rest of the distribution. The
distributions are plotted in Figure 38.

Figure 38: Comaparison between $\not{E}$ and $\not{M}$ distribution for
${\Delta m}=200$ GeV and $m_{DM}=50$ GeV

  * •

${\Delta m}=450$ GeV and $m_{DM}=20$ GeV

Lastly, we consider a scenario where, HDSP is massive and therefore produced
almost at rest. In such cases, as discussed above, momenta of the lepton and
DM produced at each end will almost fully cancel each other. Therefore, the
leptons from the decay of both HDSP’s will not necessarily be produced back to
back. However, here the DM being very light ensures the magnitudes of lepton
momenta are substantial. Therefore, in this case $\not{M}$ distribution the
maximum noticeable departure from the $\not{E}$ distributoon. The deviation is
largest at the left tail-end of the distribution due to considerable lepton
energy in that region. In Figure 39, we show the corresponding distributions.

Figure 39: Comaparison between $\not{E}$ and $\not{M}$ distribution for
${\Delta m}=450$ GeV and $m_{DM}=20$ GeV

We mention here that even in cases where the two distributions are different,
the difference between $\not{E}$ and $\not{M}$ are found to be  $<$ [-0.07cm]
$\sim$ 10%. However, one should remember that in principle it is possible to a
have discernible difference between the two distributions, especially at the
tail-end, when leptons have considerable energy and are not back to back. At
limited centre of energy of ILC such situations are rare to occur. But with
high energy muon-colliders there can be a large kinematic region where such a
situation can arise and in that case the two distributions can be
significantly different.

##  Appendix C Constraints on the model

Here we will breifly discuss possible theoretical and experimental constraints
on model parameters.

Stability: In order to ensure stability of the scalar potential, $V(H,\Phi)$
(including SM Higgs and inert-scalar doublet) along any field direction, the
quartic couplings of the scalar potential must obey following co-positivity
conditionsKannike:2012pe ; Chakrabortty:2013mha :

$\displaystyle\lambda_{H},\,\lambda_{\Phi}>0\,;$
$\displaystyle\lambda_{1}+2\sqrt{\lambda_{H}\lambda_{\Phi}}>0\,;$
$\displaystyle\lambda_{1}+\lambda_{2}-|\lambda_{3}|+2\sqrt{\lambda_{H}\lambda_{\Phi}}>0\,.$
(41)

Perturbativity: The theory remains perturbative at any given scale if the
scalar quartic couplings and the Yukawa coupling satisfy the following limits:

$\displaystyle|\lambda_{H}|<{4\pi}~{};~{}~{}|\lambda_{\Phi}|<{4\pi}~{};~{}|\lambda_{1,2,3}|<{4\pi}~{}~{}{\rm
and}~{}~{}|Y|<\sqrt{4\pi}~{}.$ (42)

Collider constraints: Higgs decay width measurement at LHC constrained the
Higgs invisible decay branching as: $Br(h\to{\rm inv})\leq 19\%$ Aaboud_2019 ;
CMS:2018yfx  which further translates into our model parameters when Higgs,
$h$ decay into a pair of dark particles is kinematically allowed.
The observed LEP data Pierce:2007ut  restrict the decay of gauge bosons ($W$
and $Z$) to BSM scalars and provide bounds on the dark scalar fields:
$m_{Z}<m_{\phi^{0}}+m_{A^{0}}$, $m_{Z}<2m_{\phi^{\pm}}$ and
$m_{W}<m_{\phi^{\pm}}+m_{\phi^{0}}(m_{A^{0}})$. The direct search charged
scalar particles at LEP-II has excluded the mass region $m_{\phi^{\pm}}<70$
GeV Pierce:2007ut . LEP data also exclude charged fermion mass
$m_{\psi^{\pm}}<102.7$ GeV LEP .

Details of DM constraints

$\phi^{0}$$\phi^{0}$$h$$h$$\phi^{0}$$\phi^{0}$$\phi^{0}$$h$$h$$H^{0}$$H^{0}$$h$SMSM

$\phi^{0}$$\phi^{0}$$W^{+}(Z)$$W^{-}(Z)$$\phi^{0}$$\phi^{0}$$\phi^{\pm}(A^{0})$$W^{\pm}(Z)$$W^{\mp}(Z)$$\phi^{0}$$\phi^{0}$$h$$W^{+}(Z)$$W^{-}(Z)$

Figure 40: Feynman diagrams for DM annihilation to SM particles for the Scalar
DM ($\phi^{0}$).

$H^{0}$$A^{0}(H^{\pm})$$Z(W^{\pm})$SMSM$H^{0}$$H^{\pm}$$Z,A$$W^{\pm}$

$\phi^{0}$$A^{0}$$\phi^{0},A^{0}(\phi^{\pm})$$h,Z(W^{\pm})$$Z,h(W^{\pm})$$\phi^{0}$$H^{\pm}$$\phi^{\pm}(A^{0})$$W^{\pm}(Z,h,A)$$h,Z,A(W^{\pm})$

Figure 41: Feynman diagrams for DM co-annihilation to SM particles for the
Scalar DM ($\phi^{0}$) associated with heavy states ($A^{0},\phi^{\pm}$).
$\psi_{1}$$\psi_{1}$$h$f /W /Z$\overline{f}$ /W
/Z$\psi_{1}$$\psi_{1}$$\psi^{\pm}$$W$$W$$\psi_{1}$$\psi_{1}$$\psi_{2}/\psi_{1,3}$$Z/h$$Z/h$
Figure 42: Feynman diagrams for DM annihilation to SM particles for the
Fermion DM ($\psi_{1}$).

$\psi_{2}$$\psi_{1}$$Z$f/W
/Z$\overline{f}$/W/h$\psi_{2}$$\psi_{1}$$\psi^{\pm}$$W$$W$$\psi_{2}$$\psi_{1}$$\psi_{1,3}$$Z$$h$

$\psi_{3}$$\psi_{1}$$h$f/ W
/Z/h$\overline{f}$/W/Z/h$\psi_{3}$$\psi_{1}$$\psi^{\pm}$$W$$W$$\psi_{3}$$\psi_{1}$$\psi_{1,3}$$h$$h$

$\psi^{-}$$\psi_{1}$$W^{+}$f/W/W/W
$f^{\prime}$/h/Z/A$\psi^{-}$$\psi_{1}$$\psi_{2}$$W$$Z$$\psi^{-}$$\psi_{1}$$\psi_{1,3}$$W^{-}$$h$

Figure 43: Feynman diagrams for DM co-annihilation to SM particles for the
Fermion DM ($\psi_{1}$) associate with the heavy states, $\psi_{2}$,
$\psi_{3}$ and $\psi^{\pm}$.  SDMSDM$h/W/Z$FDMFDM Figure 44: Feynman diagrams
for DM conversion between scalar and fermion DM.

First, we would like to discuss the annihilation and co-annihilation processes
in detail, which contribute to the observed relic density. The scalar DM
($\phi^{0}$) can annihilate and co-annihilate with other heavy states, $A^{0}$
and $\phi^{\pm}$ to SM particles via Higgs and gauge mediated interactions as
shown in Figure 40 and Figure 41 respectively. Similarly, the fermion DM has
Higgs and gauged mediated annihilation as well as co-annihilation type number-
changing processes to SM as shown in Figure 42 and Figure 43 respectively.
Along with the standard annihilation and co-annihilation, the model have DM-DM
conversion type number changing processes where heavier DM can annihilate into
lighter ones as shown in Figure 44.

Due to gauge and Higgs mediated interactions with SM thermal bath, both scalar
and fermionic dark sector particles can be assumed to be in equilibrium in the
early universe, so that both can attain relic density in today’s universe via
thermal freeze out and therefore can assume WIMP-like nature. This naturally
enables us to keep the masses of the DM in the range of $\mathcal{O}(100){\rm
GeV}$ and provides relic density of correct order. The scalar DM component,
$\phi^{0}$ will have large annihilations and co-anihilation via SM gauge
coupling and that naturally providing under-abundance as soon as the gauge
final states start contributing, which can be easily understood by looking
into the gauged mediated Feynman diagrams shown in Figure 40 and Figure 41. As
a result, the scalar DM mass in the range $m_{W}\lesssim m_{\phi^{0}}\lesssim
525$ GeV becomes underabundant ($\Omega_{\rm DM}h^{2}<0.12-0.001$). This, in
fact, motivated us to think about two-component DM scenario where the presence
of second DM candidate can produce rest of the relic.

On the other hand, the fermionic dark sector where the lightest neutral
fermion $\psi_{1}$ is considered to be the second DM-component can also have
both Higgs and gauge mediated annihilation and co-annihilation interactions
with SM. However, its gauged-mediated annihilation processes
($\psi_{1}~{}\psi_{1}\to{\rm SM~{}SM}$) are suppressed by the doublet-singlet
mixing angle $\sin\theta$. Hence relic-underabundance through annihilation is
possible in general at $Z$ and Higgs resonances ($m_{\psi_{1}}\sim
m_{Z}/2,~{}m_{h}/2$). However, relic underabundace can also be acheived at
very low $\Delta m_{2}$, where DM co-annihilation becomes dominant and at
large $\Delta m_{2}$ ($\propto Y$) where Higgs-mediated annihilation to SM
fermions provide required annihilation cross-section. The presence of DM
conversion processes are also important in order to modify freezout of heavier
DM component, and subsequent relic density. We consider only the scenario
where the combined contribution from the two DMs actually satisfy the observed
relic density within 2$\sigma$.

$\phi^{0}$$n$$h$$\phi^{0}$$n$$\psi_{1}$$n$$h$$\psi_{1}$$n$ Figure 45: Feynman
diagrams for spin independent DM-neucleon scattering process for scalar DM in
left and for fermion DM in right.

On the other hand, one has to also satisfy the limits coming from direct
detection experiments, namely XENON1T XENON:2018voc  for the two-component DM
scenario. Both scalar ($\phi^{0}$) and fermion ($\psi_{1}$) DM can be detected
through Higgs mediated $t-$ channel spin-independent (SI) DM-neucleon
scattering events: $\phi^{0}~{}n\to\phi^{0}~{}n$ and
$\psi_{1}~{}n\to\psi_{1}~{}n$ as depicted in Figure 45. At the tree-level, the
DM-nucleon scattering cross-section for scalar DM, $\sigma_{\rm
SI}(\phi^{0})\propto\lambda_{L}^{2}/m_{\phi^{0}}^{2}$, while for fermion DM,
$\sigma_{\rm SI}(\psi_{1})\propto Y^{2}/m_{\psi_{1}}^{2}$ ($Y\propto\Delta
m_{2}\sin\theta$). It is clear that in case of fermion DM, large $\Delta
m_{2}$ will result in large direct-detection cross-section, and therefore be
disfavored from the data, unless, the mixing is suppressed, ie. $\sin\theta$
is small. Such small values of $\sin\theta$ will of course lead to small
annihilation and relic over-abundance. Therefore, $Z$ or Higgs resonant region
will be the only allowed range for fermion DM in such cases.

Since both scalar and fermion DM contribute to today’s DM density, the direct
detection cross-section for each individual component should be modified by
the fraction ($\Omega_{i}/\Omega_{\rm DM}$) with which it is present in the
present universe. Hence the effective direct search cross-section can be
modified as $\sigma_{\rm SI}^{\rm
eff}(\phi^{0})=\big{(}\Omega_{\phi^{0}}/\Omega_{\rm DM}\big{)}\sigma_{\rm
SI}(\phi^{0})$ and $\sigma_{\rm SI}^{\rm
eff}(\psi_{1})=\big{(}\Omega_{\psi_{1}}/\Omega_{\rm DM}\big{)}\sigma_{\rm
SI}(\psi_{1})$.

## References

  * (1) V. C. Rubin and W. K. Ford, Jr., _Rotation of the Andromeda Nebula from a Spectroscopic Survey of Emission Regions_ , _Astrophys. J._ 159 (1970) 379–403.
  * (2) F. Zwicky, _On the Masses of Nebulae and of Clusters of Nebulae_ , _Astrophys. J._ 86 (1937) 217–246.
  * (3) E. Hayashi and S. D. M. White, _How Rare is the Bullet Cluster?_ , _Mon. Not. Roy. Astron. Soc._ 370 (2006) L38–L41, [astro-ph/0604443].
  * (4) W. Hu and S. Dodelson, _Cosmic microwave background anisotropies_ , _Ann. Rev. Astron. Astrophys._ 40 (2002) 171–216, [astro-ph/0110414].
  * (5) WMAP collaboration, G. Hinshaw et al., _Nine-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Parameter Results_ , _Astrophys. J. Suppl._ 208 (2013) 19, [1212.5226].
  * (6) WMAP collaboration, D. N. Spergel et al., _Wilkinson Microwave Anisotropy Probe (WMAP) three year results: implications for cosmology_ , _Astrophys. J. Suppl._ 170 (2007) 377, [astro-ph/0603449].
  * (7) Planck collaboration, N. Aghanim et al., _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6, [1807.06209].
  * (8) G. Bertone, D. Hooper and J. Silk, _Particle dark matter: Evidence, candidates and constraints_ , _Phys. Rept._ 405 (2005) 279–390, [hep-ph/0404175].
  * (9) L. Roszkowski, E. M. Sessolo and S. Trojanowski, _WIMP dark matter candidates and searches—current status and future prospects_ , _Rept. Prog. Phys._ 81 (2018) 066201, [1707.06277].
  * (10) E. W. Kolb and M. S. Turner, _The Early Universe_ , _Front. Phys._ 69 (1990) 1–547.
  * (11) L. J. Hall, K. Jedamzik, J. March-Russell and S. M. West, _Freeze-In Production of FIMP Dark Matter_ , _JHEP_ 03 (2010) 080, [0911.1120].
  * (12) XENON collaboration, E. Aprile et al., _Dark Matter Search Results from a One Ton-Year Exposure of XENON1T_ , _Phys. Rev. Lett._ 121 (2018) 111302, [1805.12562].
  * (13) PandaX-4T collaboration, Y. Meng et al., _Dark Matter Search Results from the PandaX-4T Commissioning Run_ , _Phys. Rev. Lett._ 127 (2021) 261802, [2107.13438].
  * (14) M. Aoki, M. Duerr, J. Kubo and H. Takano, _Multi-Component Dark Matter Systems and Their Observation Prospects_ , _Phys. Rev._ D86 (2012) 076015, [1207.3318].
  * (15) Z.-P. Liu, Y.-L. Wu and Y.-F. Zhou, _Enhancement of dark matter relic density from the late time dark matter conversions_ , _Eur. Phys. J._ C71 (2011) 1749, [1101.4148].
  * (16) Q.-H. Cao, E. Ma, J. Wudka and C. P. Yuan, _Multipartite dark matter_ , 0711.3881.
  * (17) S. Bhattacharya, A. Drozd, B. Grzadkowski and J. Wudka, _Two-Component Dark Matter_ , _JHEP_ 10 (2013) 158, [1309.2986].
  * (18) S. Esch, M. Klasen and C. E. Yaguna, _A minimal model for two-component dark matter_ , _JHEP_ 09 (2014) 108, [1406.0617].
  * (19) A. Karam and K. Tamvakis, _Dark Matter from a Classically Scale-Invariant $SU(3)_{X}$_, _Phys. Rev._ D94 (2016) 055004, [1607.01001].
  * (20) A. Ahmed, M. Duch, B. Grzadkowski and M. Iglicki, _Multi-Component Dark Matter: the vector and fermion case_ , _Eur. Phys. J._ C78 (2018) 905, [1710.01853].
  * (21) A. Poulin and S. Godfrey, _Multicomponent dark matter from a hidden gauged SU(3)_ , _Phys. Rev._ D99 (2019) 076008, [1808.04901].
  * (22) M. Aoki and T. Toma, _Boosted Self-interacting Dark Matter in a Multi-component Dark Matter Model_ , _JCAP_ 1810 (2018) 020, [1806.09154].
  * (23) S. Yaser Ayazi and A. Mohamadnejad, _Scale-Invariant Two Component Dark Matter_ , _Eur. Phys. J._ C79 (2019) 140, [1808.08706].
  * (24) M. Aoki, D. Kaneko and J. Kubo, _Multicomponent Dark Matter in Radiative Seesaw Models_ , _Front.in Phys._ 5 (2017) 53, [1711.03765].
  * (25) A. Biswas, D. Majumdar, A. Sil and P. Bhattacharjee, _Two Component Dark Matter : A Possible Explanation of 130 GeV $\gamma-$ Ray Line from the Galactic Centre_, _JCAP_ 1312 (2013) 049, [1301.3668].
  * (26) S. Bhattacharya, P. Poulose and P. Ghosh, _Multipartite Interacting Scalar Dark Matter in the light of updated LUX data_ , _JCAP_ 1704 (2017) 043, [1607.08461].
  * (27) S. Bhattacharya, P. Ghosh, T. N. Maity and T. S. Ray, _Mitigating Direct Detection Bounds in Non-minimal Higgs Portal Scalar Dark Matter Models_ , _JHEP_ 10 (2017) 088, [1706.04699].
  * (28) B. Barman, S. Bhattacharya and M. Zakeri, _Multipartite Dark Matter in $SU(2)_{N}$ extension of Standard Model and signatures at the LHC_, _JCAP_ 1809 (2018) 023, [1806.01129].
  * (29) S. Bhattacharya, P. Ghosh and N. Sahu, _Multipartite Dark Matter with Scalars, Fermions and signatures at LHC_ , _JHEP_ 02 (2019) 059, [1809.07474].
  * (30) S. Bhattacharya, P. Ghosh, A. K. Saha and A. Sil, _Two component dark matter with inert Higgs doublet: neutrino mass, high scale validity and collider searches_ , _JHEP_ 03 (2020) 090, [1905.12583].
  * (31) D. Borah, R. Roshan and A. Sil, _Minimal Two-component Scalar Doublet Dark Matter with Radiative Neutrino Mass_ , 1904.04837.
  * (32) S. Chakraborti and P. Poulose, _Interplay of Scalar and Fermionic Components in a Multi-component Dark Matter Scenario_ , 1808.01979.
  * (33) S. Chakraborti, A. Dutta Banik and R. Islam, _Probing Multicomponent Extension of Inert Doublet Model with a Vector Dark Matter_ , 1810.05595.
  * (34) S. Bhattacharya, A. K. Saha, A. Sil and J. Wudka, _Dark Matter as a remnant of SQCD Inflation_ , _JHEP_ 10 (2018) 124, [1805.03621].
  * (35) C. E. Yaguna and O. Zapata, _Fermion and scalar two-component dark matter from a $Z_{4}$ symmetry_, 2112.07020.
  * (36) G. Belanger, A. Mjallal and A. Pukhov, _Two dark matter candidates: The case of inert doublet and singlet scalars_ , _Phys. Rev. D_ 105 (2022) 035018, [2108.08061].
  * (37) D. Van Loi, N. M. Duc and P. V. Dong, _Dequantization of electric charge: Probing scenarios of cosmological multi-component dark matter_ , 2106.12278.
  * (38) C. E. Yaguna and O. Zapata, _Two-component scalar dark matter in Z 2n scenarios_, _JHEP_ 10 (2021) 185, [2106.11889].
  * (39) B. Díaz Sáez, K. Möhling and D. Stöckinger, _Two real scalar WIMP model in the assisted freeze-out scenario_ , _JCAP_ 10 (2021) 027, [2103.17064].
  * (40) N. Chakrabarty, R. Roshan and A. Sil, _Two Component Doublet-Triplet Scalar Dark Matter stabilising the Electroweak vacuum_ , 2102.06032.
  * (41) C. H. Nam, D. Van Loi, L. X. Thuy and P. Van Dong, _Multicomponent dark matter in noncommutative $B-L$ gauge theory_, _JHEP_ 12 (2020) 029, [2006.00845].
  * (42) A. Betancur, G. Palacio and A. Rivera, _Inert doublet as multicomponent dark matter_ , _Nucl. Phys. B_ 962 (2021) 115276, [2002.02036].
  * (43) D. Nanda and D. Borah, _Connecting Light Dirac Neutrinos to a Multi-component Dark Matter Scenario in Gauged $B-L$ Model_, _Eur. Phys. J. C_ 80 (2020) 557, [1911.04703].
  * (44) S. Bhattacharya, N. Chakrabarty, R. Roshan and A. Sil, _Multicomponent dark matter in extended $U(1)_{B-L}$: neutrino mass and high scale validity_, _JCAP_ 04 (2020) 013, [1910.00612].
  * (45) F. Elahi and S. Khatibi, _Multi-Component Dark Matter in a Non-Abelian Dark Sector_ , _Phys. Rev. D_ 100 (2019) 015019, [1902.04384].
  * (46) J. Herrero-Garcia, A. Scaffidi, M. White and A. G. Williams, _Time-dependent rate of multicomponent dark matter: Reproducing the DAMA/LIBRA phase-2 results_ , _Phys. Rev. D_ 98 (2018) 123007, [1804.08437].
  * (47) A. Das, S. Gola, S. Mandal and N. Sinha, _Two-component scalar and fermionic dark matter candidates in a generic U $(1)_{X}$ model_, 2202.01443.
  * (48) S. Bhattacharya, S. Chakraborti and D. Pradhan, _Electroweak Symmetry Breaking and WIMP-FIMP Dark Matter_ , 2110.06985.
  * (49) A. Dutta Banik, M. Pandey, D. Majumdar and A. Biswas, _Two component WIMP–FImP dark matter model with singlet fermion, scalar and pseudo scalar_ , _Eur. Phys. J. C_ 77 (2017) 657, [1612.08621].
  * (50) S.-M. Choi, J. Kim, P. Ko and J. Li, _A multi-component SIMP model with U(1) X→ Z2 × Z3_, _JHEP_ 09 (2021) 028, [2103.05956].
  * (51) J. Herrero-Garcia, A. Scaffidi, M. White and A. G. Williams, _On the direct detection of multi-component dark matter: sensitivity studies and parameter estimation_ , _JCAP_ 1711 (2017) 021, [1709.01945].
  * (52) J. Herrero-Garcia, A. Scaffidi, M. White and A. G. Williams, _On the direct detection of multi-component dark matter: implications of the relic abundance_ , _JCAP_ 1901 (2019) 008, [1809.06881].
  * (53) J. Hernandez-Sanchez, V. Keus, S. Moretti, D. Rojas-Ciofalo and D. Sokolowska, _Complementary Probes of Two-component Dark Matter_ , 2012.11621.
  * (54) G. F. Giudice, B. Gripaios and R. Mahbubani, _Counting dark matter particles in LHC events_ , _Phys. Rev. D_ 85 (2012) 075019, [1108.1800].
  * (55) G. Bélanger et al., _LHC-friendly minimal freeze-in models_ , _JHEP_ 02 (2019) 186, [1811.05478].
  * (56) S. P. Liew, M. Papucci, A. Vichi and K. M. Zurek, _Mono-X Versus Direct Searches: Simplified Models for Dark Matter at the LHC_ , _JHEP_ 06 (2017) 082, [1612.00219].
  * (57) F. Kahlhoefer, _Review of LHC Dark Matter Searches_ , _Int. J. Mod. Phys. A_ 32 (2017) 1730006, [1702.02430].
  * (58) A. Boveia and C. Doglioni, _Dark Matter Searches at Colliders_ , _Ann. Rev. Nucl. Part. Sci._ 68 (2018) 429–459, [1810.12238].
  * (59) D. Abercrombie et al., _Dark Matter benchmark models for early LHC Run-2 Searches: Report of the ATLAS/CMS Dark Matter Forum_ , _Phys. Dark Univ._ 27 (2020) 100371, [1507.00966].
  * (60) J. Abdallah et al., _Simplified Models for Dark Matter Searches at the LHC_ , _Phys. Dark Univ._ 9-10 (2015) 8–23, [1506.03116].
  * (61) B. Barman, S. Bhattacharya, S. Girmohanta and S. Jahedi, _Catch ’em all: Effective Leptophilic WIMPs at the $e^{+}\,e^{-}$ Collider_, 2109.10936.
  * (62) E. M. Dolle and S. Su, _The Inert Dark Matter_ , _Phys. Rev. D_ 80 (2009) 055012, [0906.1609].
  * (63) M. Dutta, S. Bhattacharya, P. Ghosh and N. Sahu, _Singlet-Doublet Majorana Dark Matter and Neutrino Mass in a minimal Type-I Seesaw Scenario_ , _JCAP_ 03 (2021) 008, [2009.00885].
  * (64) _Searches for supersymmetric particles in e + e- collisions up to 208 gev and interpretation of the results within the mssm_ , _The European Physical Journal C_ 31 (Dec, 2003) 421–479.
  * (65) A. Pierce and J. Thaler, _Natural Dark Matter from an Unnatural Higgs Boson and New Colored Particles at the TeV Scale_ , _JHEP_ 08 (2007) 026, [hep-ph/0703056].
  * (66) C. Degrande, C. Duhr, B. Fuks, D. Grellscheid, O. Mattelaer and T. Reiter, _UFO - The Universal FeynRules Output_ , _Comput. Phys. Commun._ 183 (2012) 1201–1214, [1108.2040].
  * (67) G. Bélanger, F. Boudjema, A. Goudelis, A. Pukhov and B. Zaldivar, _micrOMEGAs5.0 : Freeze-in_ , _Comput. Phys. Commun._ 231 (2018) 173–186, [1801.03509].
  * (68) MAGIC, Fermi-LAT collaboration, M. L. Ahnen et al., _Limits to Dark Matter Annihilation Cross-Section from a Combined Analysis of MAGIC and Fermi-LAT Observations of Dwarf Satellite Galaxies_ , _JCAP_ 02 (2016) 039, [1601.06590].
  * (69) Fermi-LAT collaboration, M. Ackermann et al., _Searching for Dark Matter Annihilation from Milky Way Dwarf Spheroidal Galaxies with Six Years of Fermi Large Area Telescope Data_ , _Phys. Rev. Lett._ 115 (2015) 231301, [1503.02641].
  * (70) H. Abramowicz et al., _The International Linear Collider Technical Design Report - Volume 4: Detectors_ , 1306.6329.
  * (71) K. Kannike, _Vacuum Stability Conditions From Copositivity Criteria_ , _Eur. Phys. J._ C72 (2012) 2093, [1205.3781].
  * (72) J. Chakrabortty, P. Konar and T. Mondal, _Copositive Criteria and Boundedness of the Scalar Potential_ , _Phys. Rev._ D89 (2014) 095008, [1311.5666].
  * (73) M. Aaboud, G. Aad, B. Abbott, O. Abdinov, B. Abeloos, D. Abhayasinghe et al., _Search for invisible higgs boson decays in vector boson fusion at s=13tev with the atlas detector_ , _Physics Letters B_ 793 (Jun, 2019) 499–519.
  * (74) CMS collaboration, A. M. Sirunyan et al., _Search for invisible decays of a Higgs boson produced through vector boson fusion in proton-proton collisions at $\sqrt{s}=$ 13 TeV_, _Phys. Lett. B_ 793 (2019) 520–551, [1809.05937].