# Lepton portal dark matter at muon colliders: Total rates and generic
features for phenomenologically viable scenarios

Adil Jueid<EMAIL_ADDRESS>Center for Theoretical Physics of the
Universe, Institute for Basic Science (IBS), 34126 Daejeon, Republic of Korea
and Salah Nasri<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of physics,
United Arab Emirates University, Al-Ain, UAE  The Abdus Salam International
Centre for Theoretical Physics, Strada Costiera 11, I-34014, Trieste, Italy

###### Abstract

Lepton portal dark matter (DM) models are a class of models where the DM
candidates solely couple to charged leptons through a mediator carrying a
lepton number. These models are very interesting since they avoid constraints
from direct detection experiments even for coupling of order ${\cal O}(1)$,
they have small annihilation cross sections, and can be probed efficiently at
lepton colliders. In this work, we consider a minimal lepton portal DM model
which consists of extending the SM with two $SU(2)_{L}$ singlets: a charged
scalar singlet and an electrically neutral right-handed fermion. We
systematically study the production mechanisms of DM at multi-TeV muon
colliders. After considering all the possible theoretical and experimental
constraints and studying the phenomenology of lepton flavour violation and DM
in the muon-philic scenario, we analyse the production rates of 54 channels
(26 channels for prompt DM production and 28 channels for charged scalar
production) at multi-TeV muon colliders. Finally, we discuss the possible
collider signatures of some channels and the corresponding backgrounds. We
find that at least 9 channels for DM production can be very efficient in
testing DM with masses up to about $1$ TeV.

††preprint: CTPU-PTC-2023-02

##  I Introduction

Supported by various astrophysical and cosmological observations, it is now
widely accepted that dark matter (DM) exists in the universe (see e.g.
Jungman:1995df ; Bergstrom:2000pn ; Bertone:2004pz ; Feng:2010gw  for
comprehensive reviews). On the other hand, the measurements of the
anisotropies in the cosmic microwave background (CMB) implies that DM is the
dominant component of the matter budget in the universe with a density of
$\Omega_{\rm DM}h^{2}=0.1198\pm 0.0015$ Ade:2015xua . The standard theories of
structure formation require that the DM should be non-relativistic at the
matter-radiation equality. In particle physics models, this can be easily
realised by extending the SM with weakly interacting massive particles (WIMPs)
under the standard thermal freeze-out mechanism. The search for WIMPs was one
of the major programmes at the Large Hadron Collider (LHC). A special
characteristic of WIMPs production at the LHC is that one can probe it through
the recoil of a SM particle against a large missing transverse energy
($E_{T}^{\rm miss}$). Examples of these processes are mono-jet Beltran:2010ww
, mono-$Z$ Carpenter:2012rg ; Bell:2012rg  or mono-Higgs Berlin:2014cfa  among
others. Unfortunately various searches for WIMPs at the LHC were unsuccessful
to find such signals and limits were put on the production cross section
versus the DM mass ATLAS:2021kxv ; ATLAS:2021shl ; ATLAS:2021hza ; CMS:2019ykj
; CMS:2020ulv  which were interpreted in various particle physics
realizations. Furthermore, these constraints were even more stringent when the
void bounds from direct detection experiments Aprile:2018dbl ; Cui:2017nnn
are included GAMBIT:2017gge ; GAMBIT:2017zdo . The situation is not very
different in the case the DM production is mediated through colored mediators
or leptoquarks with the main mechanisms for DM density in the early universe
being the co-annihilation or conversion-driven freezeout mechanisms
Baker:2015qna ; Choi:2018stw ; Mohan:2019zrk ; Belanger:2021smw ;
Belfatto:2021ats . The interpretation of these searches exclude DM masses of
about $0.1$–$1$ TeV and mediator masses of about $0.5$–$5$ TeV depending on
the theoretical model.

In the light of this current situation, an important question arises: what if
DM only couples to the lepton sector? From the theoretical standpoint, there
is a priori no fundamental principle that can prevent DM from coupling to
leptons only. This class of models has been proposed some time ago in ref.
Liu:2013gba  and was widely studied in the literature Bai:2014osa ;
Chang:2014tea ; Agrawal:2014ufa ; Garny:2015wea ; Jueid:2020yfj ;
Horigome:2021qof ; Liu:2021mhn . There are many interesting implications for
these models. First, the scattering of the DM off the nucleus is induced at
the one-loop order and therefore these models can evade easily direct
detection constraints even for model parameters of order $\mathcal{O}(1)$.
Second, except for electron-philic scenarios, constraints from positron
indirect detection searches are also not important since their annihilation is
dominated by $p$-wave amplitudes which are suppressed by the square of the DM
velocity. Finally, the DM can be produced at the LHC through the decay of
charged scalars and therefore the corresponding bounds are not as strong as in
the case of mono–X searches especially in the case of $SU(2)_{L}$ gauge
singlet mediators Jueid:2020yfj . Therefore an efficient probe of this
category of models is through leptonic colliders such as the International
Linear Collider (ILC), Chinese Electron Positron Collider (CEPC), and the
future muon colliders. Recently, future muon colliders are attracting high
interest due to their capability to probe new physics beyond the SM at very
high scales Delahaye:2019omf ; Long:2020wfp ; AlAli:2021let  and therefore
competing with the future circular colliders (FCC–hh). On the other hand,
these machines can achieve very high energies thanks to the expected excellent
cooling systems and the weaker synchrotron radiation. Finally at very high
energies, muon colliders are necessarily vector-boson colliders where the
dominant production channels are through vector-boson fusion (VBF)
Costantini:2020stv ; Ruiz:2021tdt . Phenomenology of both the SM and beyond at
muon colliders has been extensively studied in the literature (see e.g.
Capdevilla:2020qel ; Chiesa:2020awd ; Han:2020uid ; Han:2020uak ; Yin:2020afe
; Huang:2021nkl ; Capdevilla:2021rwo ; Capdevilla:2021fmj ; Asadi:2021gah ;
Casarsa:2021rud ; Liu:2021akf ; Han:2021udl ; Han:2021kes ; Han:2021lnp ;
Lv:2022pts ; Liu:2022byu ; Azatov:2022itm ; Yang:2022fhw ; Bao:2022onq ;
Chen:2022msz  and references therein).

In this work, we study the production of DM at muon colliders within the
minimal lepton portal DM model in which we extends the SM with two $SU(2)_{L}$
singlets: a charged scalar that plays the role of the mediator and a neutral
right-handed fermion (or, equivalently, Majorana particle) that plays the role
of the DM candidate. We first comprehensively the impact of the different
theoretical and experimental constraints on the model parameter space in the
muon-philic scenario, i.e. the scenario where the DM couples predominantly to
muons. We then select a few benchmark points that define phenomenologically
viable scenarios that can be probed at high energy muon colliders. We study
the production cross sections and the expected backgrounds for a set of
production channels totaling 26 production channels for DM and 28 production
channels for the charged singlet scalar. A particular feature of this model is
that the DM is a Majorana fermion and therefore does not couple to gauge
bosons directly and therefore the direct production of DM does not receive any
contribution from VBF channels. We select a few production channels that can
have high discovery potential and discuss the possible signatures and the
associated backgrounds. This work is an introduction for future projects where
a complete exploration of the model at muon colliders will be performed.

The remainder of this paper is orgnized as follows. We discuss the model and
its UV completion in section II along with the constraints from LEP searches,
$H_{\rm SM}\to\gamma\gamma$ and theoretical constraints. In section III we
discuss the constraints from charged lepton flavour violation in
$\ell_{\alpha}\to\ell_{\beta}\gamma$, $\ell_{\alpha}\to 3\ell_{\beta}$ and
$H_{\rm SM}\to\ell_{\alpha}\bar{\ell}_{\beta}$. A detailed analysis of DM
phenomenology in this model is presented in section IV where we discuss the DM
relic density, direct detection constraints and Higgs invisible decays. A
study of DM production at muon colliders, the interesting signatures and the
associated backgrounds is performed in section V. In section VI we study the
production of charged scalars at muon colliders. We draw our conclusions in
section VII.

##  II Theoretical setup

###  II.1 The model

We consider a minimal extension of the SM by two gauge singlet fields: a
charged scalar ($S$) and a right-handed fermion ($N_{R}$). We further assume
that the two extra singlets are odd under $Z_{2}$ symmetry while all the SM
particles are even; i.e., $\\{S,N_{R}\\}\to\\{-S,-N_{R}\\}$ and
$\\{\ell,q,\nu,\Phi,V^{\mu}\\}\to\\{\ell,q,\nu,\Phi,V^{\mu}\\}$. To ensure
that the $N_{R}$ state is a suitable DM candidate within our model, we impose
the condition $M_{N_{R}}<M_{S}$. Furthermore, the charged singlet is assumed
to carry a lepton number and therefore couples only to charged leptons.111This
charged singlet is also called a scalar lepton Baum:2020gjj  and the relevant
interaction Lagrangian is similar to the case of interaction of supersymmetric
slepton with a neutralino and a charged lepton. The difference here is that we
assume a single charged scalar to couple to all the leptons instead of three
scalars, usually denoted by $\tilde{e}_{R},\tilde{\mu}_{R},~{}{\rm
and}~{}\tilde{\tau}_{R}$, where each scalar couples to a specific lepton
generation. The full Lagrangian is given by

$\displaystyle\mathcal{L}=\mathcal{L}_{\rm SM}+\mathcal{L}_{S}-V(\Phi,S),$ (1)

where $\Phi$ refers to the SM Higgs doublet, $\mathcal{L}_{S}$ is the
interaction Lagrangian for the singlet scalar (including the kinetic term),
and $V(\Phi,S)$ is the scalar potential. The interaction Lagrangian for the
$S$ field is given by

$\displaystyle\mathcal{L}_{S}=\sum_{\ell=e,\mu,\tau}Y_{\ell
N}\overline{\ell}^{c}_{R}SN_{R}+(\mathcal{D}^{\mu}S)^{\dagger}(\mathcal{D}_{\mu}S)+{\rm
h.c.},$ (2)

with $\mathcal{D}_{\mu}S=(\partial_{\mu}-ig_{2}Y_{S}B_{\mu}/2)S$ being the
covariant derivative, $Y_{S}=2$ is the hypercharge of the scalar singlet and
$g_{2}$ is the $U(1)_{Y}$ gauge coupling. The kinetic term in equation (2)
gives rise to interaction with $A_{\mu}$ and $Z_{\mu}$ which are given, after
field rotations, by

$\displaystyle\mathcal{L}_{S;\rm gauge}$ $\displaystyle=$
$\displaystyle-(eA^{\mu}-e\tan\theta_{W}Z^{\mu})S^{\dagger}\bar{\partial}_{\mu}S+e^{2}A_{\mu}A^{\mu}S^{\dagger}S$
$\displaystyle+$ $\displaystyle
e^{2}\tan^{2}\theta_{W}Z_{\mu}Z^{\mu}S^{\dagger}S-2e^{2}\tan\theta_{W}A_{\mu}Z^{\mu}S^{\dagger}S,$

where $e=\sqrt{4\pi\alpha_{\rm EM}}$ is the electric charge, $\theta_{W}$ is
the Weinberg mixing angle, and $A\bar{\partial_{\mu}}B\equiv
A(\partial_{\mu}B)-(\partial_{\mu}A)B$. The most general CP-conserving,
renormalizable and gauge invariant scalar potential is given by

$V(\Phi,S)=-M_{11}^{2}|\Phi^{\dagger}\Phi|+M_{22}^{2}|S^{\dagger}S|+\lambda_{1}|\Phi^{\dagger}\Phi|^{2}+\lambda_{2}|S^{\dagger}S|^{2}+\lambda_{3}|\Phi^{\dagger}\Phi||S^{\dagger}S|.$
(3)

All the parameters of the scalar potential are assumed to be real valued as a
consequence of CP conservation. The process of electroweak symmetry breaking
leads to three physical scalars: $H_{\rm SM}$ identified with the recently
discovered $125~{}{\rm GeV}$ SM Higgs boson and a pair of charged scalars
denoted by $H^{\pm}$. Their masses are given at the lowest order in
perturbation theory by

$\displaystyle M_{H_{\rm SM}}^{2}=\lambda_{1}\upsilon^{2}=-2M_{11}^{2},\quad
M_{H^{\pm}}^{2}=M_{22}^{2}+\frac{1}{2}\lambda_{3}\upsilon^{2},$ (4)

with $\upsilon$ being the vacuum expectation value (VEV) of the SM Higgs
doublet. This model involves seven additional free parameters which we
parametrise as follows

$\displaystyle\\{M_{H^{\pm}},M_{N_{R}},\lambda_{2},\lambda_{3},Y_{eN},Y_{\mu
N},Y_{\tau N}\\}.$ (5)

For convenience we define the combination of the couplings $Y_{\ell N}$
by222This is equivalent to a definition of a system of spherical coordinates
wherein the new parameters are $Y_{\ell N},\theta$ and $\varphi$ such that
$\theta\in[0,\pi]$ and $\varphi\in[0,2\pi]$. The couplings in equation (2) are
defined here as $Y_{eN}=Y_{\ell N}\cos\varphi\sin\theta$, $Y_{\mu N}=Y_{\ell
N}\sin\varphi\sin\theta$ and $Y_{\tau N}=Y_{\ell N}\cos\theta$.

$Y_{\ell N}=\sqrt{Y_{eN}^{2}+Y_{\mu N}^{2}+Y_{\tau N}^{2}},$

which is a very good parametrisation in case the charged leptons are assumed
to be massless.

###  II.2 Theoretical and experimental constraints

The parameters of the model in equation (5) are subject to various theoretical
and experimental constraints. We start with a brief discussion of the
constraints influencing the scalar potential parameters and $M_{H^{\pm}}$
where more details can be found in Jueid:2020yfj . The width of the SM Higgs
boson is only affected by the rate of its decay to $\gamma\gamma$. In this
model, this process receives new contributions from charged singlet scalar
which give rise to destructive or constructive contributions depending on the
sign of $\lambda_{3}$ Swiezewska:2012eh ; Arhrib:2012ia ; Jueid:2020rek 333We
have found a typo in the analytical expression in ref. Arhrib:2012ia  which
may influences their numerical results.. In the present work, we have used the
most recent ATLAS-CMS combined measurement of $|\kappa_{\gamma}|$
Khachatryan:2016vau

$|\kappa_{\gamma}|\equiv\sqrt{\Gamma(H\to\gamma\gamma)/\Gamma(H\to\gamma\gamma)_{\mathrm{SM}}}=0.87^{+0.14}_{-0.09}.$

We assume the theoretical prediction to be in agreement with the experimental
measurement at the $2\sigma$–level. We found that the enhancement of
$|\kappa_{\gamma}|$ always occur for $\lambda_{3}<0$ which excludes charged
scalars with masses up to $\sim 380~{}{\rm GeV}$ Jueid:2020rek . For
$\lambda_{3}>0$, we get three possible regimes: (i) large and negative
contribution that implies an enhancement of $\kappa_{\gamma}$, (ii) positive
but small contribution which makes $\kappa_{\gamma}$ consistent with the
experimental measurement and (iii) exact or almost exact cancellation between
the $H^{\pm}$ and the $W$-boson contributions which make $\kappa_{\gamma}$
very small. Therefore, for $\lambda_{3}>0$, charged singlet masses up to
$380~{}{\rm GeV}$ are excluded but with small region where the constraints
completely vanish.

In addition to constraints from Higgs decays, the parameters of the scalar
potential are subject to a number of theoretical constraints. We note that the
bounds on the scalar potential of this model can be obtained from those in
e.g. the inert doublet model by setting $\lambda_{4}=\lambda_{5}=0$. In this
study, we impose constraints from vacuum stability conditions (or boundness-
from-below) Branco:2011iw , perturbativity, perturbative unitarity
Kanemura:1993hm ; Akeroyd:2000wc  and false vacuum Ginzburg:2010wa . The false
vacuum condition plays a very important role in constraining the parameters
$\lambda_{2},\lambda_{3}$ and $M_{H^{\pm}}$. We get

$\displaystyle
M_{H^{\pm}}^{2}\geqslant\frac{1}{2}\bigg{(}\lambda_{3}\upsilon^{2}-M_{H_{\rm
SM}}^{2}\sqrt{\frac{\lambda_{2}}{\lambda_{1}}}\bigg{)}.$ (6)

We found that: (i) $\lambda_{3}$ cannot be larger $5$ for all charged scalar
masses and (ii) there is a parabola in the plane defined by $\lambda_{3}$ and
$M_{H^{\pm}}$ which simply tells us that the smaller is the minimum allowed
value of $M_{H^{\pm}}$ the smaller is the maximum allowed value of
$\lambda_{3}$. These conclusions are mildly dependent on the choice of
$\lambda_{2}$ and, therefore, we choose $\lambda_{2}=2$ in the remainder of
this manuscript without loss of generality.

Figure 1: Summary of the collider constraints on the parameter space of the
model displayed on the plane of $(M_{H^{\pm}},M_{N_{R}})$. We show the
constraints from LEP searches of sleptons and charginos (red), LHC searches
for sleptons in the compressed regime (blue) and constraints from LHC searches
of sleptons and charginos for large mass splittings (green). Here, we assume
that the charged singlet scalar decays to $\mu^{\pm}N_{R}$ with a branching
fraction of $100\%$ and assume the Narrow Width Approximation (NWA) by
selecting parameters for which we have $\Gamma_{H^{\pm}}/M_{H^{\pm}}<0.15$.
The gray dashed line corresponds to the kinematical boundary above which the
$N_{R}$ particle is not a suitable dark matter candidate.

The model can be constrained by using the null results of LEP and LHC searches
for supersymmetric particles Abbiendi:2003ji ; ATLAS:2019lff ; ATLAS:2019lng .
The OPAL collaboration of the LEP experiment has searched for charginos
decaying into a charged lepton and the lightest supersymmetric neutralino
using $680~{}{\rm pb}^{-1}$ of integrated luminosity Abbiendi:2003ji .
Assuming that the branching ratio of $H^{\pm}\to\mu^{\pm}N_{R}$ is $100\%$,
the production of charged singlet pairs occurs through gauge interactions
($s$–channel diagrams with the exchange of $\gamma^{*}/Z^{0}$). This search
constrain the mass of the charged singlet to be not heavier than $100~{}{\rm
GeV}$ for any value of $Y_{\mu N}$. This can be seen clearly in the red
contour of figure 1. The ATLAS collaboration at the LHC has also searched for
sleptons and charginos assuming $100\%$ branching fraction to a charged lepton
and neutralino. These searches targeted large mass splitting
$\Delta=m_{\tilde{\ell}}-m_{\chi^{0}}\geq 80~{}{\rm GeV}$ ATLAS:2019lff  and
compressed spectra for a mass splitting as low as $0.55~{}{\rm GeV}$
ATLAS:2019lng . The two searches utilized a total luminosity of $139~{}{\rm
fb}^{-1}$. The first search constrained scalar singlet masses to be lighter
than about $440$ GeV while the search for compressed spectra constrain the
whole compressed region for $M_{H^{\pm}}$ up to $150$ GeV (see the blue and
green contours of figure 1).

###  II.3 Examples of UV completions

In this section, we discuss the UV completions of this minimal framework. In
general, there are two ways to UV complete the first term in
$\mathcal{L}_{S}$: (i) assume it to be a part of a radiative neutrino mass
model or (ii) embed it in a grand-unified theory – $SU(5)$ for example –. We
start by the radiative neutrino mass models. The most economical way to extend
this model is through the so-called Krauss-Nasri-Trodden (KNT) three-loop
radiative neutrino mass model Krauss:2002px . In addition to $S$ and $N_{R}$,
the KNT model extends the SM with an additional scalar singlet that is even
under $Z_{2}$. Another possibility is through the so-called the scotogenic
model which extends the SM with one inert doublet and three right-handed
fermions Ma:2006km . The phenomenology of the scotogenic model has been widely
studied in the literature Ahriche:2017iar ; Baumholzer:2018sfb ; Borah:2018rca
; Kitabayashi:2018bye ; Ahriche:2018ger ; Baumholzer:2019twf ; Borah:2020wut ;
Sarazin:2021nwo . The relevant interaction becomes

$\displaystyle\mathcal{L}\supset
h_{\alpha\beta}\bar{L}_{L\alpha}(i\sigma_{2})\Phi_{\rm IDM}N_{\beta}\supset
h_{\alpha\beta}\bar{\ell}_{L\alpha}SN_{\beta},$ (7)

where $\Phi_{\rm IDM}=(S,(h_{2}+ia_{2})/\sqrt{2})^{T}$, and $\alpha,\beta$ are
generation indices. Identifying (7) with the first term in $\mathcal{L}_{S}$
we have $Y_{eN}=h_{11},Y_{\mu N}=h_{21},Y_{\tau N}=h_{31}$. We must stress out
that the gauge interactions of the singlet scalar in this model are different
from the scotogenic model due to the fact that $S$ is a member of $SU(2)_{L}$
doublet while it is a singlet in the present model.

The first term in equation (2) can be obtained from a grand-unified theory;
For example, by embeding the SM into a $SU(5)$ gauge group with the matter
fields belonging to the ${\bf{10}}_{F}$ and $\bar{{\bf{5}}}_{F}$
representations, the charged singlet belongs to the ${\bf{10}}_{H}$
representation, and the right handed neutrino belongs to the singlet
representation ${\bf{1}}_{\alpha}$, which in this case we can write

$\displaystyle\mathcal{L}_{\text{int}}=g_{\alpha\beta}\overline{{\bf{10}}}_{\alpha}\otimes{\bf{10}}_{H}\otimes{\bf{1}}_{N_{\beta}}\supset
g_{\alpha\beta}\ell_{R\alpha}^{T}CN_{\beta}S^{+}.$ (8)

In addition to the minimal $SU(5)$, we can obtain the first term of equation
(2) from a flipped-$SU(5)\otimes U(1)_{X}$ grand-unified theory. Here, the
right-handed charged lepton field is a singlet under $SU(5)$ while the right-
handed neutral fermion ($N_{R}$) is a member of the ${\bf{10}}_{\alpha}$
representation. In this case, we have

$\displaystyle\mathcal{L}_{\text{int}}$ $\displaystyle=$
$\displaystyle\frac{h_{\alpha\beta}}{\Lambda}\overline{{\bf{10}}}_{\alpha}\otimes\bar{{\bf{1}}}_{\beta}\otimes{\bf{10}}_{H}\otimes{\bf{1}}_{S}+h.c.$
(9) $\displaystyle\supset$
$\displaystyle\frac{h_{\alpha\beta}\langle{\bf{10}}_{H}\rangle}{\Lambda}N^{T}C\ell_{R}S^{-},$

where we integrated out a heavy intermediate state with a scale
$\Lambda\gg\Lambda_{\rm GUT}$.

##  III Charged Lepton Flavour violation

The interaction Lagrangian in equation (2) conserves total lepton number to
all orders in perturbation theory since the charged singlet possesses a lepton
number444To generate a Majorana neutrino mass one has extend the Lagrangian
(2) so that the total lepton number is violated. The minimal realization of
such breaking can be achieved by having, in addition to $S$, a second $SU(2)$
singlet charged scalar with lepton number equals to two units; which is the
KNT model.. However, the charged singlet scalar can give rise to processes
violating flavor lepton numbers $L_{\alpha};\alpha=e,\mu,\tau$ at the one-loop
order. These processes called charged lepton flavor violating (CLFV) processes
are categorised into three categories: (i) $\ell_{\alpha}=\ell_{\beta}\gamma$,
(ii) $\ell_{\alpha}\to\ell_{\beta}\ell_{\beta}\bar{\ell}_{\beta}$ and (iii)
$e$–$\mu$ conversion in nuclei. In this section we discuss the impact of the
CLFV constraints on the model parameter space. The most stringent bounds on
the couplings $Y_{\ell_{\alpha}N}$ come from the branching ratio of $\mu\to
e\gamma$ decay. The analysis of the CLFV decays in this work are heavily based
on the results of refs. Hisano:1995cp ; Arganda:2005ji ; Ilakovac:2012sh ;
Toma:2013zsa . A summary of the current and future bounds on the CLFV decays
is shown in Table 1.

CLFV decay |  Present limit | Future sensitivity
---|---|---
$\mu\to e\gamma$ |  $5.7\times 10^{-13}$ Adam:2013mnn  |  $6\times 10^{-14}$ Baldini:2013ke
$\tau\to e\gamma$ |  $3.3\times 10^{-8}$ Aubert:2009ag  |  $\sim 10^{-8}-10^{-9}$ Hayasaka:2013dsa
$\tau\to\mu\gamma$ |  $4.4\times 10^{-8}$ Aubert:2009ag  |  $\sim 10^{-8}-10^{-9}$ Hayasaka:2013dsa
$\mu\to eee$ |  $1.0\times 10^{-12}$ SINDRUM:1987nra  |  $\sim 10^{-16}$ Blondel:2013ia
$\tau\to eee$ |  $2.7\times 10^{-8}$ Hayasaka:2010np  |  $\sim 10^{-9}-10^{-10}$ Hayasaka:2013dsa
$\tau\to\mu\mu\mu$ |  $2.1\times 10^{-8}$ Hayasaka:2010np  |  $\sim 10^{-9}-10^{-10}$ Hayasaka:2013dsa
$H_{\rm SM}\to\mu\tau$ |  $1.5\times 10^{-3}$ CMS:2021rsq  | $-$
$H_{\rm SM}\to e\tau$ |  $2.2\times 10^{-3}$ CMS:2021rsq  | $-$
$H_{\rm SM}\to e\mu$ |  $3.5\times 10^{-4}$ CMS:2016cvq  | $-$
Table 1: Current experimental bounds and future sensitivities for low-energy
CLFV decays and high-energy Higgs boson LFV decays.

###  III.1 $\ell_{\alpha}\to\ell_{\beta}\gamma$

The radiative decays of charged leptons ($\ell_{\alpha}\to\ell_{\beta}\gamma$)
receive contributions from the exchange of the charged singlet scalar and
Majorana DM. After computing the one-loop integrals we get the effective
magnetic dipole operator
$\mu^{M}_{\beta\alpha}\overline{\ell}_{\beta}\sigma^{\mu\nu}\ell_{\alpha}F_{\mu\nu}/2$
with $\mu_{\beta\alpha}=em_{\alpha}A_{\rm M}/2$ and $A_{\rm M}$ is given by

$A_{\rm
M}=\frac{Y_{\ell_{\beta}N}Y_{\ell_{\alpha}N}}{2(4\pi)^{2}}\frac{1}{M_{H^{\pm}}^{2}}\mathcal{F}(\xi),$

where $\xi=M_{N_{R}}^{2}/M_{H^{\pm}}^{2}$ and
$\mathcal{F}(x)=(1-6x+3x^{2}+2x^{3}-6x^{2}\log x)/(6(1-x)^{4}$ is the one-loop
function which have the following limits $\mathcal{F}(x)\to 1/6~{}(1/12)$ for
$x\to 0~{}(1)$. The resulting decay branching ratio can be computed easily to
give

${\rm
BR}\left(\ell_{\alpha}\to\ell_{\beta}\gamma\right)=\frac{3(4\pi)^{3}\alpha_{\rm
EM}}{4G_{F}^{2}}|A_{\rm M}|^{2}\times{\rm
BR}\left(\ell_{\alpha}\to\ell_{\beta}\nu_{\alpha}.\overline{\nu}_{\beta}\right),\hskip
14.22636pt$ (10)

Here, $G_{F}=1.166\times 10^{-5}~{}{\rm GeV}^{-2}$, $\alpha_{\rm EM}=1/137$,
and ${\rm
BR}\left(\ell_{\alpha}\to\ell_{\beta}\nu_{\alpha}\overline{\nu_{\beta}}\right)$
is the SM decay branching ratios. We choose ${\rm BR}(\mu\to
e\nu\bar{\nu}),{\rm BR}(\tau\to e\nu\bar{\nu}),{\rm
BR}(\tau\to\mu\nu\bar{\nu})\approx 1,0.1783,0.1741$ ParticleDataGroup:2018ovx
.

Using the most recent experimental bounds on ${\rm
BR}(\ell_{\alpha}\to\ell_{\beta}\gamma)$ from the Meg Adam:2013mnn  and BaBar
Aubert:2009ag  experiments, we can use equation (10) to derive the following
bounds on the products of the couplings:

$\displaystyle|Y_{eN}Y_{\mu N}|<\bigg{(}\frac{2.855\times
10^{-5}}{\mathrm{GeV}}\bigg{)}^{2}\frac{M_{H^{\pm}}^{2}}{|\mathcal{F}(\xi)|},$
$\displaystyle|Y_{eN}Y_{\tau N}|<\bigg{(}\frac{4.428\times
10^{-4}}{\mathrm{GeV}}\bigg{)}^{2}\frac{M_{H^{\pm}}^{2}}{|\mathcal{F}(\xi)|},$
(11) $\displaystyle|Y_{\tau N}Y_{\mu N}|<\bigg{(}\frac{4.759\times
10^{-4}}{\mathrm{GeV}}\bigg{)}^{2}\frac{M_{H^{\pm}}^{2}}{|\mathcal{F}(\xi)|}.$

Since the one-loop function varies roughly between $1/12$ and $1/6$, the upper
bound on the coupling $Y_{\ell_{\alpha}N}Y_{\ell_{\beta}N}$ is proportional to
the square of the charged singlet mass with almost no dependence on
$M_{N_{R}}$. Therefore, limits are expected to be strong for light $H^{\pm}$
and become very weak for heavy $H^{\pm}$. This can be clearly seen in figure 2
where the maximum allowed values of $|Y_{\ell_{\alpha}N}Y_{\ell_{\beta}}N|$ by
the CLFV decays ${\rm BR}(\ell_{\alpha}\to\ell_{\beta}\gamma)$ are shown as a
function of $\xi=M_{N_{R}}^{2}/M_{H^{\pm}}^{2}$ for
$M_{H^{\pm}}=500,1000,~{}{\rm and}~{}5000~{}{\rm GeV}$. As expected, the
bounds on $|Y_{eN}Y_{\mu N}|$ are the strongest ones while the bounds on
$|Y_{\tau N}Y_{\mu N}|$ and $|Y_{eN}Y_{\tau N}$ are similar. We must stress
out that the CLFV decays do not constraints the $Y_{\ell_{\alpha}N}$ couplings
per se but only their products. Following this finding, there is some freedom
regarding the choice of the couplings which we call here benchmark scenarios
(see next sections). Given that this study is mainly concerned about the
phenomenology of the leptophilic DM models at muon colliders, we choose a
scenario where the coupling of dark matter to the muon is quite large while
the other couplings are chosen such that they fulfill the experimental bounds
on CLFV decays: $Y_{\mu N}\simeq\mathcal{O}(1)\gtrsim Y_{\tau N}\gg Y_{eN}$.

Figure 2: The maximum value of the products of the
$Y_{\ell_{\alpha}N}Y_{\ell_{\beta}N}$ as a function of $\xi$ for different
values of the charged singlet mass $M_{H^{\pm}}$. The results are shown for
$M_{H^{\pm}}=500~{}{\rm GeV}$ (dashed), $M_{H^{\pm}}=1000~{}{\rm GeV}$
(dotted) and $M_{H^{\pm}}=5000~{}{\rm GeV}$ (dash-dotted).

###  III.2 $\ell_{\alpha}\to\ell_{\beta}\ell_{\beta}\bar{\ell}_{\beta}$

It is noteworthy to discuss the constraints from the CLFV decays
$\ell_{\alpha}\to\ell_{\beta}\ell_{\beta}\bar{\ell}_{\beta}$. These processes
receive four contributions at the one-loop order: penguin diagrams with the
exchange of $\gamma$, $Z$ and $H_{\rm SM}$ and box diagrams. The contribution
of the SM Higgs boson is suppressed due to the smallness of the Higgs-lepton
Yukawa coupling. The corresponding branching ratio is given by Toma:2013zsa

$\displaystyle{\rm
BR}(\ell_{\alpha}\to\ell_{\beta}\ell_{\beta}\bar{\ell_{\beta}})$
$\displaystyle=$ $\displaystyle\frac{3(4\pi)^{2}\alpha_{\rm
EM}}{8G_{F}^{2}}\bigg{[}\overbrace{|A_{\rm ND}|^{2}+|A_{\rm
M}|^{2}\left(\frac{16}{3}\log\left(\frac{m_{\alpha}}{m_{\beta}}\right)-\frac{22}{3}\right)}^{\gamma~{}{\rm
penguin}}+\overbrace{\frac{1}{3}(2|Z_{\rm RR}|^{2}+|Z_{\rm
RL}|^{2})}^{Z~{}{\rm penguin}}$ (12) $\displaystyle+$
$\displaystyle\frac{1}{6}|B_{\rm box}|^{2}+\underbrace{2~{}{\rm
Re}\left(-2A_{\rm ND}A_{\rm M}^{*}+\frac{1}{3}A_{\rm ND}B_{\rm
box}^{*}-\frac{2}{3}A_{\rm M}B_{\rm box}^{*}\right)}_{{\rm
Interference}}\bigg{]}\times\mathcal{B},$

where $\mathcal{B}\equiv{\rm
BR}(\ell_{\alpha}\to\ell_{\beta}\nu_{\alpha}\bar{\nu}_{\beta})$. The
contribution of the $\gamma$-penguins consist of the magnetic or dipole
($A_{\rm M}$) and the non-dipole ($A_{\rm ND}$) contributions. The dipole
contribution is the same as of ${\rm BR}(\ell_{\beta}\to\ell_{\beta}\gamma)$
but enhanced by a factor of $16\times(\log(m_{\alpha}/m_{\beta})-22)/3$ which
varies between $7$ and $36$ for $\tau\to 3\mu$ and $\tau\to 3e$ respectively.
The non-dipole contribution is given by

$A_{\rm
ND}=\frac{Y_{\ell_{\alpha}N}Y_{\ell_{\beta}N}}{6(4\pi)^{2}}\frac{1}{M_{H^{\pm}}^{2}}\mathcal{G}(\xi),$

with $\mathcal{G}(x)=(2-9x+18x^{2}-11x^{3}+6x^{3}\log x)/(6(1-x)^{4})$ being
the one-loop function for the non-dipole $\gamma$–penguin. This function has
the following limits: $\lim_{x\to 0}\mathcal{G}(x)=1/3$ and $\lim_{x\to
1}\mathcal{G}(x)=1/4$. Therefore, the dipole $\gamma$-penguin contribution is
large as compared to the non-dipole contributions; $\lim_{x\to
0}~{}(\lim_{x\to 1})A_{\rm M}/A_{\rm
ND}\times(16/3\log(m_{\alpha}/m_{\beta})-22/3)\approx\\{3.5,11,18\\}~{}(\\{2,7,12\\})$
for $\tau\to 3\mu,\mu\to 3e,~{}{\rm and}~{}\tau\to 3e$ respectively. The
$Z$-penguin contribution is given by

$\displaystyle Z_{\rm RR}=\frac{g_{R}^{\ell}Z_{\rm
ND}}{g_{1}^{2}\sin^{2}\theta_{W}M_{Z}^{2}},~{}Z_{\rm
RL}=\frac{g_{L}^{\ell}Z_{\rm ND}}{g_{1}^{2}\sin^{2}\theta_{W}M_{Z}^{2}},$ (13)

where $g_{R}^{\ell},g_{L}^{\ell}$ are the right and left-handed components of
the $Z$-boson couplings to charged leptons, $g_{1}$ is the $SU(2)_{L}$ gauge
coupling, $\sin\theta_{W}$ is the sine of the Weinberg mixing angle, and
$Z_{\rm ND}$ is the momentum-independent $Z$-boson form factor which is given
by

$Z_{\rm
ND}=\frac{Y_{\ell_{\alpha}N}Y_{\ell_{\beta}N}}{2(4\pi)^{2}}\frac{m_{\alpha}m_{\beta}}{M_{H^{\pm}}^{2}}\frac{g_{1}}{\cos\theta_{W}}\mathcal{F}(\xi).$

We can see that the $Z$-penguin contribution involves an extra suppression by
a factor of $m_{\alpha}m_{\beta}$ as compared to the dipole
$\gamma$-contribution. Finally, the box contribution is given by

$\displaystyle B_{\rm
box}=\frac{Y_{\ell_{\alpha}N}Y_{\ell_{\beta}N}^{3}}{2^{7}\pi^{3}\alpha_{\rm
EM}M_{H^{\pm}}^{2}}\bigg{[}\mathcal{D}_{1}(\xi)+2\xi\mathcal{D}_{2}(\xi)\bigg{]},$
(14)

where $\mathcal{D}_{1,2}(x)$ are the one-loop box functions given by
$\mathcal{D}_{1}(x)=(-1+x^{2}-2x\log x)/(1-x)^{3}$ and
$\mathcal{D}_{2}(x)=(-2+2x-(1+x)\log x)/(1-x)^{3}$. The contribution of the
box diagrams, contrarily to penguins, has an extra factor of
$Y_{\ell_{\beta}N}^{2}$. Therefore, it may dominate for large couplings of the
daughter lepton to DM. In this work, we check that the benchmark scenarios
satisfy the bounds from the $\ell_{\alpha}\to 3\ell_{\beta}$ decays (see Table
2).

###  III.3 $H_{\rm SM}\to\ell_{\alpha}\bar{\ell}_{\beta}$

We close this section by a brief discussion of the CLFV decays of the SM Higgs
boson. These decays have been searched for by the ATLAS and the CMS
collaborations with the most strongest bounds are reported on by CMS
collaboration CMS:2021rsq ; CMS:2016cvq . In this model, the CLFV decays of
the SM Higgs boson are degenerate to the radiative CLFV decays of the charged
leptons. The constraints from CLFV of charged leptons imply that the CLFV
decays of the SM Higgs boson are extremely suppressed and may even be beyond
the future reach of the LHC and future colliders. The SM Higgs boson decay
into $\ell_{\alpha}\ell_{\beta}$ is given by Herrero-Garcia:2016uab

$\displaystyle{\rm BR}(H_{\rm SM}\to\ell_{\alpha}\ell_{\beta})\simeq 1.2\times
10^{3}\times|y_{\ell_{\alpha}}Y_{\ell_{\alpha}N}Y_{\ell_{\beta}N}|^{2}\bigg{(}\frac{\lambda_{3}}{4\pi}\bigg{)}^{2}\bigg{(}\frac{\upsilon}{M_{H^{\pm}}}\bigg{)}^{4},$
(15)

with $y_{\ell_{\alpha}}=m_{\ell_{\alpha}}/(\sqrt{2}\upsilon)$ is the Higgs-
lepton Yukawa coupling of the heavier lepton (chosen here to be
$\ell_{\alpha}$). In this formula, the contribution of the lighter lepton is
neglected. We expect the bounds from $H_{\rm
SM}\to\ell_{\alpha}\bar{\ell}_{\beta}$ searches to be very weak. This can be
clearly seen in from table 2 for the benchmark points we have used in this
study.

##  IV Dark Matter

In this section, we discuss the DM phenomenology within this model. We start
with the calculation of the relic density of the$N_{R}$ particles in section
IV.1 and then move to a detailed analysis of the spin-independent DM-nucleus
scattering cross section in section IV.2. Next, we derive the constraint on
the couplings $Y_{lN}$ by analysing the the Higgs invisible decays and
conclude by a selection of the benchmark points that are compatible with all
the theoretical and the experimental constraints in section IV.4.

###  IV.1 Relic density

The relic density of the $N_{R}$ particles receives contributions from both
the annihilation and the co-annihilation. The co-annihilation becomes active
when the mass splitting $\Delta\equiv M_{H^{\pm}}-M_{N_{R}}<0.1\times
M_{N_{R}}$ while the annihilation contributes for the whole parameter space.
For the annihilation, there are two major contributions: (i)
$N_{R}N_{R}\to\ell_{\alpha}^{+}\ell_{\beta}^{-}$ from the exchange of the
charged scalar singlet in $t$\- and $u$-channels, and (ii)
$N_{R}N_{R}\to\sum_{X\in{\rm SM}}X\overline{X}$ which arises from the exchange
of the SM Higgs boson via $s$-channel diagrams. Note that $s$–channel
contributions to the relic density are negligible in our model if one demands
perturbativity of the couplings. The reason is that the leading order
contribution to the $s$-channel annihilation amplitudes arises at the one-loop
order. To obtain the relic density of the $N_{R}$ particles, one must solve
the Boltzmann equations given by Scherrer:1985zt ; Griest:1990kh ;
Jungman:1995df

$\displaystyle\frac{{\rm d}n_{N_{R}}}{{\rm
d}t}+3Hn_{N_{R}}=-2\langle\sigma_{N_{R}}v\rangle\bigg{[}(n_{N_{R}})^{2}-(n_{N_{R}}^{\rm
eq})^{2}\bigg{]},$ (16)

with $H=\dot{a}/a$, $n_{N_{R}}$ is the number density of the $N_{R}$ particle
and $n_{N_{R}}^{\rm eq}\approx
g_{N_{R}}\left(\frac{M_{N_{R}}T}{2\pi}\right)^{3/2}e^{-M_{N_{R}}/T}$ is its
number density at the thermal equilibrium. Note that in the absence of
interactions that change the number density of $N_{R}$, the right handed side
of equation (16) would be equal to zero and $n_{N_{R}}\propto a^{-3}$. This
equation can be solved to give approximately

$\displaystyle\Omega_{\rm DM}h^{2}\simeq\frac{3\times 10^{-27}~{}{\rm
cm}^{3}{\rm s}^{-1}}{\langle\sigma(x_{f})v\rangle},$ (17)

where $\langle\sigma(x_{f})v\rangle$ is the thermally-averaged annihilation
cross section for the $N_{R}$ particle

$\displaystyle\langle\sigma(x_{f})v\rangle=\frac{1}{8M_{N_{R}}^{4}T_{f}K_{2}^{2}(M_{N_{R}}/T_{f})}\sum_{\alpha,\beta}\int_{4M_{N_{R}}^{2}}^{\infty}{\rm
d}\hat{s}\sqrt{\hat{s}-4M_{N_{R}}^{2}}K_{1}(\sqrt{\hat{s}}/T_{f})\sigma_{N_{R}N_{R}\to\ell_{\alpha}\ell_{\beta}}(\hat{s}),$
(18)

where $K_{1}(x)$ and $K_{2}(x)$ are the modified Bessel functions of the
second kind and $\sigma_{N_{R}N_{R}\to\ell_{\alpha}\ell_{\beta}}(\hat{s})$ is
the annihilation cross section into charged lepton which is given by

$\sigma_{N_{R}N_{R}\to\ell_{\alpha}\ell_{\beta}}(\hat{s})=\frac{1}{2^{3}\pi}\frac{|Y_{\ell_{\alpha}N}Y_{\ell_{\beta}N}|^{2}}{\hat{s}~{}\hat{\kappa}_{1}}\bigg{[}(m_{\ell_{\alpha}}^{2}+m_{\ell_{\beta}}^{2})(\hat{s}-2M_{N_{R}}^{2})+\frac{1}{6}\frac{\hat{\kappa}_{2}}{\hat{\kappa}_{1}}\hat{s}(\hat{s}-4M_{N_{R}}^{2})\bigg{]},$
(19)

where
$\hat{\kappa}_{i}\equiv\hat{\kappa}_{i}(M_{H^{\pm}}^{2},M_{N_{R}}^{2},\hat{s})$,
$\hat{\kappa}_{1}(x,y,z)=(2x+2y-z)^{2}$ and
$\hat{\kappa}_{2}(x,y,z)=(4x-4y+z)^{2}-2z^{2}$. To simplify the discussion
about the relic density, we consider the annihilation cross section in the
limit $\hat{s}\to 4M_{N_{R}}^{2}$

$\displaystyle\sigma_{N_{R}N_{R}\to\ell_{\alpha}\ell_{\beta}}\approx\frac{|Y_{\ell_{\alpha}N}Y_{\ell_{\beta}N}|^{2}}{2^{6}\pi
M_{H^{\pm}}^{4}}(m_{\ell_{\alpha}}^{2}+m_{\ell_{\beta}}^{2})\bigg{(}1+\frac{M_{N_{R}}^{2}}{M_{H^{\pm}}^{2}}\bigg{)}^{-2}.$

This equation simply tells us that the contribution of the annihilation to the
relic density becomes very small for very heavy charged singlet scalar and one
needs to have large $Y_{\ell N}$ to produce the correct relic density. On the
other hand, for large values of the mass splitting and heavy charged singlet
scalar one cannot reproduce the correct relic abundance if one demands
perturbativity of the couplings. The co-annihilations are more involved in
this model as we can have additional contributions that have different
dependence on the model parameters. There are two generic co-annihilation
channels: $N_{R}H^{\pm}\to{\rm SM}$ and $H^{\pm}H^{\mp}\to{\rm SM}$. Below, we
list the individual contributions and the overall dependence of the
corresponding cross section

$\displaystyle N_{R}H^{\pm}\to\ell_{\alpha}^{\pm}H_{\rm SM}$ $\displaystyle:$
$\displaystyle\quad\sigma\propto\lambda_{3}^{2}Y_{\ell_{\alpha}N}^{2},$
$\displaystyle
N_{R}H^{\pm}\to\ell_{\alpha}^{\pm}Z,\ell_{\alpha}^{\pm}\gamma,\nu W^{\pm}$
$\displaystyle:$ $\displaystyle\quad\sigma\propto Y_{\ell_{\alpha}N}^{2},$
$\displaystyle H^{\pm}H^{\mp}\to\ell_{\alpha}^{\pm}\ell_{\beta}^{\mp}$
$\displaystyle:$
$\displaystyle\quad\sigma\propto|Y_{\ell_{\alpha}N}Y_{\ell_{\beta}N}|^{2}\mathcal{A}_{1}+|Y_{\ell_{\alpha}N}Y_{\ell_{\beta}N}|\mathcal{A}_{2}+\mathcal{A}_{3},$
$\displaystyle H^{\pm}H^{\mp}\to q\bar{q}$ $\displaystyle:$
$\displaystyle\quad\sigma\propto\lambda_{3}^{2}\mathcal{B}_{1}+\lambda_{3}\mathcal{B}_{2}+\mathcal{B}_{3},$
$\displaystyle H^{\pm}H^{\mp}\to ZZ,H_{\rm SM}Z,W^{\pm}W^{\mp}$
$\displaystyle:$
$\displaystyle\quad\sigma\propto\lambda_{3}^{2}\mathcal{C}_{1}+\lambda_{3}\mathcal{C}_{2}+\mathcal{C}_{3},$
$\displaystyle H^{\pm}H^{\mp}\to H_{\rm SM}H_{\rm SM}$ $\displaystyle:$
$\displaystyle\quad\sigma\propto\lambda_{3}^{4}\mathcal{D}_{1}+\lambda_{3}^{2}\mathcal{D}_{2},$

with $\mathcal{A}_{i},\mathcal{B}_{i},\mathcal{C}_{i}$ and $\mathcal{D}_{i}$
are real-valued coefficients that depend on the dark matter mass, the charged
singlet scalar mass and the final-state particles. The co-annihilation becomes
very active for quite large $\lambda_{3}$ and $Y_{\ell N}$ and may even drive
the relic density to very small values ($\sim$ orders of magnitudes smaller
than the observed abundance). In general, the co-annihilation is dominated by
contributions of the following two processes $H^{\pm}H^{\mp}\to 2H_{\rm SM}$
and $N_{R}H^{\pm}\to\ell_{\alpha}^{\pm}H_{\rm SM}$. In the presence of co-
annihilations, the Boltzmann equations become

$\displaystyle\frac{{\rm d}n_{N_{R}}}{{\rm d}t}+3Hn_{N_{R}}$ $\displaystyle=$
$\displaystyle-2\langle\sigma_{\rm
eff}v_{r}\rangle\bigg{[}(n_{N_{R}})^{2}-(n_{N_{R}}^{\rm
eq})^{2}\bigg{]}+N\Gamma_{H^{\pm}}n_{H^{\pm}},$ (20) $\displaystyle\frac{{\rm
d}n_{H^{\pm}}}{{\rm d}t}+3Hn_{H^{\pm}}$ $\displaystyle=$
$\displaystyle-\Gamma_{H^{\pm}}n_{H^{\pm}},$ (21)

where $N$ is the mean number of $N_{R}$ particles, $n_{H^{\pm}}$ is the number
density of $H^{\pm}$ and $\Gamma_{H^{\pm}}$ is its total width. Note that here
we have replaced the thermally-averaged annihilation cross section in equation
(16) by the effective cross section

$\displaystyle\langle\sigma_{\rm
eff}v_{r}\rangle=\sum_{i,j\in\\{N_{R},H^{\pm}\\}}\langle\sigma(ij\to{\rm
SM})v_{r}\rangle\frac{n_{i}^{\rm eq}n_{j}^{\rm eq}}{(n_{N_{R}}^{\rm
eq})^{2}}.$ (22)

The relic density of $N_{R}$ is obtained from the numerical solutions of the
coupled Boltzmann equations (21). MadDM version 3.0 is used to solve the
Boltzmann equations and compute the relic density of $N_{R}$ Ambrogi:2018jqj .
In figure 3, we show the values of the coupling $Y_{\ell N}$ consistent with
the measurement of the relic density by the Planck collaboration projected on
the mass of the dark matter and the mass of the charged singlet scalar. We can
see that the relic abundance of the $N_{R}$ is consistent with the Planck
measurement only for very specific regions. If the mass splitting between
$H^{\pm}$ and $N_{R}$ is large, we need large values of the $Y_{\ell N}$.
However, even for $Y_{\ell N}$ near the perturbativity bound the mass
splitting can not be arbitrary large: $\Delta_{\rm max}\approx
600~{}(2000)~{}{\rm GeV}$ for $M_{N_{R}}=10~{}(100)~{}{\rm GeV}$. The relic
density becomes almost independent of $Y_{\ell N}$ for large $M_{N_{R}}$ in
the co-annihilation regions. We conclude this section by noting that the model
can not reproduce the correct relic density with the standard freeze-out
mechanism for the region marked in blue in figure (3) as it breaks the
perturbativity of the coupling $Y_{\ell N}$.

Figure 3: Values of the coupling $Y_{\ell N}$ consistent with the measurement
of the relic density by the Planck collaboration projected on the mass of the
dark matter and the mass of the charged singlet scalar. The isolines
corresponding to $\Omega h^{2}\approx 0.12$ are shown for $Y_{\ell
N}=1,2,3,4,5,7.5$ and $4\pi$. The blue shaded area corresponds to the region
where the perturbativity is broken while the shaded gray region correspond to
the kinematically forbidden region $M_{N_{R}}>M_{H^{\pm}}$ in which $N_{R}$ is
not stable and therefore not a suitable dark matter candidate.

###  IV.2 Direct detection

We turn now into a discussion of the constraints from direct detection
experiments on the model parameter space. In this model, the scattering cross
section of $N_{R}$ off the nucleus with atomic number ($A$) occurs at the one-
loop order where the SM Higgs boson plays the role of a portal. The generic
formula for the spin-independent cross section is given by555The spin-
dependent cross section is very small in our model as the exchanged particle
is the SM Higgs boson which is a scalar particle with $J^{P}=0^{+}$.
Nevertheless, we will compute this observable for some benchmark scenarios and
estimate their consistency with the current bounds from the PICO experiment
PICO:2017tgi .

$\displaystyle\sigma_{\rm
SI}=\frac{4}{\pi}\mu_{A}^{2}\bigg{(}Z\cdot\mathcal{S}_{p}+(A-Z)\cdot\mathcal{S}_{n}\bigg{)}^{2},$
(23)

with $\mathcal{S}_{p,n}$ being the scalar current nucleon ($p/n$) form factors
and $\mu_{A}\equiv M_{N_{R}}m_{A}/(M_{N_{R}}+m_{A})$ is the reduced mass of
the $N_{R}$–$A$ system. The nucleon form factors have two contributions: (i)
from particle physics which is connected to the scattering amplitude of the
$N_{R}$–$(q/g)$ process and (ii) from low-energy nuclear physics that are
computed using chiral perturbation theory DelNobile:2013sia ; Hill:2014yxa ;
Bishara:2017pfq ; Ellis:2018dmb .

Figure 4: The spin-independent cross section as function of the dark matter
mass $M_{N_{R}}$ while the colored scatter points correspond to the charged
singlet scalar mass $M_{H^{\pm}}$. In the same plot, we show the current
bounds from Xenon 1T Aprile:2018dbl  in dashed sienna, and the future
expectations from Xenon nT XENON:2020kmp , LUX LZ LZ:2018qzl  and DarkSide G2
DarkSide-20k:2017zyg . The shaded orange area marked by ‘neutrino floor’
corresponds to the backgrounds from the coherent scattering with solar
neutrinos, atmospheric neutrinos and supernova neutrinos Billard:2013qya . The
spin-independent cross section was scaled by a factor of $\xi_{\rm
Planck}=\Omega_{N_{R}}h^{2}/\Omega_{\rm Planck}h^{2}$ with $\Omega_{\rm
Planck}h^{2}\approx 0.12$. All the calculations were performed for $Y_{\ell
N}=2$ and $\lambda_{3}=4$.

The generic formula of $f_{p,n}$ is given by

$\displaystyle\mathcal{S}_{p,n}=m_{p,n}\sum_{u,d,s}\frac{\mathcal{A}_{q}}{m_{q}}\mathcal{S}^{q}_{p,n}+\frac{2}{27}m_{p,n}\mathcal{S}_{p,n}^{g}\sum_{c,b,t}\frac{\mathcal{A}_{q}}{m_{q}},$
(24)

where $m_{p}=938.27~{}{\rm MeV}$, $m_{n}=939.56~{}{\rm MeV}$ are the proton
and neutron masses respectively. The values of the scalar nucleon low energy
form factors are chosen to be Backovic:2015cra

$\displaystyle\mathcal{S}_{p}^{u}=1.53\times
10^{-2},\quad\mathcal{S}_{p}^{d}=1.91\times
10^{-2},\quad\mathcal{S}_{p}^{s}=4.47\times
10^{-2},\quad\mathcal{S}_{p}^{g}=1-\sum_{q=u,d,s}\mathcal{S}_{p}^{q}=92.09\times
10^{-2},$ $\displaystyle\mathcal{S}_{n}^{u}=1.10\times
10^{-2},\quad\mathcal{S}_{n}^{d}=2.73\times
10^{-2},\quad\mathcal{S}_{n}^{s}=4.47\times
10^{-2},\quad\mathcal{S}_{n}^{g}=1-\sum_{q=u,d,s}\mathcal{S}_{n}^{q}=91.70\times
10^{-2}.$

The parton-level scattering amplitude is

$\displaystyle\mathcal{M}_{qN_{R}\to
qN_{R}}=\mathcal{A}_{q}\bar{\psi}_{q}(p_{\rm out})\psi_{q}(p_{\rm in}),$ (25)

where $\mathcal{A}_{q}$ is connected to the non-hadronic part of the
amplitude. The term $\bar{\psi}_{q}(p_{\rm out})\psi_{q}(p_{\rm in})$ should
be incorporated in a hadronic current
$\langle\mathcal{N}|\cdot|\mathcal{N}\rangle$

$\displaystyle\langle\mathcal{N}|\bar{\psi}_{q}\psi_{q}|\mathcal{N}\rangle=\Bigg{\\{}\begin{array}[]{lr}\frac{m_{\mathcal{N}}}{m_{q}}\cdot\mathcal{S}_{\mathcal{N}}^{q},&\text{for
}q=u,d,s,\\\
\frac{2}{27}\frac{m_{\mathcal{N}}}{m_{q}}\cdot\mathcal{S}_{\mathcal{N}}^{g},&\text{for
}q=c,b,t,\end{array}$ (28)

where $\mathcal{N}=p,n$. The model-dependent non-hadronic form factor is given
by

$\displaystyle\mathcal{A}_{q}=\frac{\tilde{y}(Q^{2}\approx 0)}{M_{H_{\rm
SM}}^{2}}\cdot\frac{m_{q}}{v}\bar{\psi}_{N_{R}}(k_{\rm
out})\psi_{N_{R}}(k_{\rm in}),$ (29)

here $\tilde{y}(Q^{2}\approx 0)$ is the effective $H_{\rm SM}N_{R}N_{R}$
coupling computed in the low energy limit. With the help of the Package X
Patel:2015tea , we can obtain it from equation (32)

$\displaystyle\tilde{y}(Q^{2}\approx 0)\simeq-\frac{\lambda_{3}v|Y_{\ell
N}|^{2}}{16\pi
M_{H^{\pm}}}\frac{1}{\varrho_{N}}\bigg{[}1-\left(1-\varrho_{N}^{-2}\right)\log{\left(1-\varrho_{N}^{2}\right)}\bigg{]}\equiv-\frac{\lambda_{3}v|Y_{\ell
N}|^{2}}{16\pi M_{H^{\pm}}}\mathcal{H}(\varrho_{N}),$ (30)

where $\varrho_{N}=M_{N_{R}}/M_{H^{\pm}}$. $\mathcal{H}(x)$ is monotonous and
increasing function of $x$ in the interval $[0,1]$ and has the following
limits $\lim_{x\to 0}\mathcal{H}(x)=0$ and $\lim_{x\to 1}\mathcal{H}(x)=1$.
Note that the first limit correspond to a small dark matter mass and a heavy
charged scalar for which the model cannot reproduce the correct relic
abundance while the second limit corresponds to the nearly degenerate scenario
where co-annihilation is the most active component in the relic abundance
calculation. In addition the effective coupling involves an extra suppression
by $1/M_{H^{\pm}}$ which simply means that the direct detection spin-
independent cross section is always below the neutrino floor for heavy
$H^{\pm}$. From equation (30) one also expect that the spin-independent cross
section is always proportional to $|Y_{\ell N}|^{4}$. Therefore, large
$Y_{\ell N}$ regions with large $\sigma_{\rm SI}$ would also correspond to
small relic density (which is proportional to $1/|Y_{\ell N}|^{4}$)666This is
consistent with our previous finding in Jueid:2020yfj  where a strong anti-
correlation between $\sigma_{\rm SI}$ and $\Omega_{N_{R}}h^{2}$ was observed.
and for these scenarios $\sigma_{\rm SI}$ needs to be scaled by a factor
$\xi_{\rm Planck}\equiv\Omega_{N_{R}}h^{2}/\Omega_{\rm Planck}h^{2}$. This
means that the spin-independent cross section would always be consistent with
the current Xenon 1T bounds Aprile:2018dbl  for most regions of the parameter
space as we can see clearly in figure 4.

Benchmark point BP1 BP2 BP3 BP4 Parameters $M_{N_{R}}~{}({\rm GeV})$ $50$
$200$ $598$ $1000$ $M_{H^{\pm}}~{}({\rm GeV})$ $500$ $500$ $600$ $1500$
$Y_{Ne}$ $10^{-4}$ $5\times 10^{-4}$ $10^{-3}$ $5\times 10^{-3}$ $Y_{N\mu}$
$2.8$ $1.6$ $1$ $2$ $Y_{N\tau}$ $5\times 10^{-2}$ $5\times 10^{-1}$ $5\times
10^{-1}$ $2$ $\lambda_{3}$ $4$ $5$ $5$ $6$ Decays of $H^{\pm}$ ${\rm
BR}(H^{\pm}\to eN_{R})$ $1.27\times 10^{-9}$ $8.89\times 10^{-8}$ $8.98\times
10^{-7}$ $3.12\times 10^{-6}$ ${\rm BR}(H^{\pm}\to\mu N_{R})$ $99.96\times
10^{-2}$ $91.10\times 10^{-2}$ $89.70\times 10^{-2}$ $50.0\times 10^{-2}$
${\rm BR}(H^{\pm}\to\tau N_{R})$ $3.18\times 10^{-4}$ $8.89\times 10^{-2}$
$10.29\times 10^{-2}$ $49.99\times 10^{-2}$ $\Gamma_{H^{\pm}}~{}({\rm GeV})$
$76.45$ $19.72$ $5.88\times 10^{-4}$ $73.68$ $\Gamma_{H^{\pm}}/M_{H^{\pm}}$
$15.29\times 10^{-2}$ $3.94\times 10^{-2}$ $9.81\times 10^{-7}$ $4.91\times
10^{-2}$ ${\rm BR}(\ell_{\alpha}\to\ell_{\beta}\gamma)$ and ${\rm
BR}(\ell_{\alpha}\to 3\ell_{\beta})$ ${\rm BR}(\mu\to e\gamma)$ $2.68\times
10^{-14}$ $1.51\times 10^{-13}$ $4.31\times 10^{-14}$ $1.89\times 10^{-13}$
${\rm BR}(\tau\to e\gamma)$ $1.52\times 10^{-18}$ $2.64\times 10^{-15}$
$1.92\times 10^{-15}$ $3.38\times 10^{-14}$ ${\rm BR}(\tau\to\mu\gamma)$
$1.17\times 10^{-9}$ $2.64\times 10^{-8}$ $1.87\times 10^{-9}$ $5.28\times
10^{-9}$ ${\rm BR}(\mu\to eee)$ $1.47\times 10^{-16}$ $8.21\times 10^{-16}$
$2.27\times 10^{-16}$ $1.01\times 10^{-15}$ ${\rm BR}(\tau\to eee)$
$1.51\times 10^{-20}$ $2.58\times 10^{-17}$ $1.85\times 10^{-17}$ $3.29\times
10^{-16}$ ${\rm BR}(\tau\to\mu\mu\mu)$ $1.21\times 10^{-8}$ $9.79\times
10^{-9}$ $2.63\times 10^{-12}$ $1.17\times 10^{-9}$ ${\rm BR}(H_{\rm
SM}\to\ell_{\alpha}\ell_{\beta})$ ${\rm BR}(H_{\rm SM}\to\mu\tau)$ $2.31\times
10^{-8}$ $1.18\times 10^{-6}$ $2.22\times 10^{-7}$ $5.24\times 10^{-7}$ ${\rm
BR}(H_{\rm SM}\to e\tau)$ $2.95\times 10^{-17}$ $1.15\times 10^{-13}$
$2.22\times 10^{-13}$ $3.27\times 10^{-12}$ Dark matter observables
$\Omega_{N_{R}}h^{2}$ $9.84\times 10^{-2}$ $9.25\times 10^{-2}$ $2.11\times
10^{-3}$ $8.53\times 10^{-2}$ $\langle\sigma\upsilon\rangle~{}({\rm cm}^{2})$
$2.40\times 10^{-9}$ $2.55\times 10^{-9}$ $7.32\times 10^{-8}$ $2.69\times
10^{-9}$ $\sigma_{\rm SI}^{p}~{}({\rm cm}^{2})$ $1.60\times 10^{-47}$
$3.45\times 10^{-47}$ $2.28\times 10^{-48}$ $1.47\times 10^{-46}$ $\sigma_{\rm
SD}^{p}~{}({\rm cm}^{2})$ $6.51\times 10^{-62}$ $6.29\times 10^{-62}$
$1.98\times 10^{-65}$ $8.29\times 10^{-60}$ XENON1T $\checkmark$ $\checkmark$
$\checkmark$ $\checkmark$ PICO $\checkmark$ $\checkmark$ $\checkmark$
$\checkmark$ DarkSide G2 $\checkmark$ X $\checkmark$ X LZ X X $\checkmark$ X
Neutrino floor X X X X

Table 2: Characteristics of the four benchmark points in our model. Here, we
show the values of the independent parameters, the decay branching ratios and
total width of the charged singlet scalar, the CLFV decay branching ratios and
dark-matter observables. A checkmark ($\checkmark$) indicate that the
parameter point yields a smaller $\sigma_{\rm SI}$ than the experimental bound
(present or expected) while a cross mark (X) indicates that $\sigma_{\rm SI}$
is above the experimental bound.

###  IV.3 Higgs invisible decay

The Higgs invisible decay occurs at the one-loop order with the exchange of
charged scalar and right-handed fermion. The partial decay width is given by

$\displaystyle\Gamma(H_{\rm SM}\to N_{R}N_{R})=\frac{M_{H_{\rm
SM}}|\tilde{y}|^{2}}{8\pi}\bigg{(}1-\frac{4M_{N_{R}}^{2}}{M_{H_{\rm
SM}}^{2}}\bigg{)}^{3/2},\hskip 14.22636pt$ (31)

with $\tilde{y}$ is the one-loop induced effective $H_{\rm
SM}$–$N_{R}$–$N_{R}$ coupling which is given by

$\displaystyle\tilde{y}=\frac{\lambda_{3}\upsilon
M_{N_{R}}}{16\pi^{2}}\sum_{\ell}|Y_{\ell N}|^{2}(C_{0}+C_{2}),$ (32)

with $C_{i}\equiv C_{i}(M_{N_{R}}^{2},M_{H_{\rm
SM}}^{2},M_{N_{R}}^{2},m_{\ell}^{2},M_{H^{\pm}}^{2},M_{H^{\pm}}^{2}),i=0,2$
being the Passarino-Veltman three-point functions Passarino:1978jh . The
computation of the Feynman amplitudes has been performed using FeynArts,
FormCalc, and LoopTools Hahn:1998yk ; Hahn:2000kx . We have used a Python
interface to LoopTools to evaluate numerically the one-loop
integrals777pylooptools is a Python binding to LoopTools and can be found in
this github directory: https://github.com/djukanovic/pylooptools.git.. We
define the Higgs invisible branching ratio as

$\displaystyle B_{\rm inv}\equiv\frac{\Gamma(H_{\rm SM}\to
N_{R}N_{R})}{\Gamma(H_{\rm SM}\to N_{R}N_{R})+\Gamma_{H}^{\rm SM}},$ (33)

where $\Gamma_{H}^{\rm SM}=4.07~{}{\rm MeV}$. Using equations (31), (32) and
(33), we can obtain bounds on the coupling $Y_{\ell N}$. The bound is
analytically defined by

$\displaystyle Y_{\ell
N}<\bigg{(}\frac{2048\pi^{5}\Gamma_{H}^{\mathrm{SM}}}{\beta_{N}^{3/2}M_{H_{\rm
SM}}\lambda_{3}^{2}\upsilon^{2}M_{N_{R}}^{2}|C_{0}+C_{2}|^{2}\left(\frac{1}{B_{\rm
bound}}-1\right)}\bigg{)}^{1/4},$

where $\beta_{N}\equiv(1-4M_{N_{R}}^{2}/M_{H_{\rm SM}}^{2})$ and $B_{\rm
bound}$ is the upper bound on ${\rm BR}_{\rm inv}$.

Figure 5: The present and future exclusions on values of $Y_{\ell N}$ for
$M_{H^{\pm}}=500~{}{\rm GeV}$ and $\lambda_{3}=4$. Here we show the contours
obtained from the LHC (navy), HL-LHC (turquoise), FCC–ee (magenta), ILC
(orange), CEPC (dark red) and FCC–hh (gray). All the bounds were obtained
assuming SM Higgs boson mass of $M_{H_{\rm SM}}=125~{}{\rm GeV}$, and SM Higgs
boson production rates. The Higgs diphoton rate is assumed to be equal to the
SM prediction at LO.

Searches for Higgs invisible decays have been carried by the ATLAS and the CMS
collaborations ATLAS:2015ciy ; CMS:2016dhk ; Sirunyan:2018owy . The strongest
and up-to-date stringent bound on $\textrm{B}_{\mathrm{inv}}$ was reported by
the Cms collaboration using a combination of previous Higgs to invisible decay
searches at $7,8$ and $13$ TeV, where it has been found
$B_{\mathrm{inv}}<B_{\rm bound}=0.19$ at $95\%$ CL Sirunyan:2018owy  assuming
that the rates of the Higgs boson production are equal to the SM predictions.
On the other hand, several groups have carried global analyses using recent
Higgs boson measurements and obtained stringent limits Belanger:2013xza ;
Ellis:2013lra . Finally, several studies have been devoted to the projected
sensitivities of the future collider experiments to Higgs invisible decays
from HL–LHC CMS-PAS-FTR-16-002 , FCC–ee Cerri:2016bew , ILC Asner:2013psa ,
CEPC CEPCStudyGroup:2018ghi  and FCC–hh Selvaggi:2018obt . In figure 5, we
show the excluded values of $Y_{\ell N}$ from present and future bounds on
${\rm BR}_{\rm inv}$ assuming $M_{H^{\pm}}=500~{}{\rm GeV}$ and
$\lambda_{3}=4$888This choice of the charged singlet scalar mass is consistent
with the limits from searches of sleptons at the LHC. We note that increasing
$M_{H^{\pm}}$ would weaken the bounds on $Y_{\ell N}$ from Higgs invisible
decays.. As we can see clearly that the present bounds are extremely weak
which excludes $Y_{\ell N}\sim 6$ for $M_{N_{R}}\sim 49~{}{\rm GeV}$. The
future experiments are expected to exclude smaller values of $Y_{\ell N}$;
e.g. FCC–hh can exclude values up to $0.7$ for $M_{N_{R}}\sim 49~{}{\rm GeV}$.

###  IV.4 Benchmark points

From the discussions in sections II and IV, we can conclude the following:

  * •

The scalar singlet cannot be lighter than $440$ GeV for mass splittings with
the dark matter of order $\leq 80$ GeV.

  * •

CLFV can constrain only the product of the Yukawa-type couplings and not their
individual values. Therefore, benchmark points have to be chosen.

  * •

DM direct detection constraints are not very strong as expected since the
spin-independent cross section is one-loop induced.

  * •

The constraints from the consistency with the measurement of the DM relic
density forbids large mass splittings if the Yukawa-type couplings are of
order $\mathcal{O}(1)$.

The benchmark points used in the discussion of the general features of DM
production at Muon colliders are shown in Table 2. There are four of these
benchmarks and each one has its own phenomenological implications.

##### BP1.

This benchmark point is characterised by a relatively light DM
($M_{N_{R}}=50~{}{\rm GeV}$) and a charged singlet mass in near the exclusion
limit reported on by the LHC (see figure 1). On the other hand, the Yukawa-
type couplings are chosen such that $Y_{\mu N}\gg Y_{\tau N}>Y_{eN}$. This
choice leads to a charged singlet decaying predominantly into $\mu N_{R}$ with
a branching fraction approaching $100\%$. On the other hand, the charged
lepton flavour violating decays of charged leptons are such that ${\rm
BR}(\tau\to e\gamma)$ is well below the sensitivity reach of future
experiments in the foreseeable future. The other branching ratios, are below
the current experimental bounds but can be tested in the near future. For DM
observables, the relic density for this BP is about $90\%$ of the observed
abundance and the spin-independent DM-nucleon cross section is below Xenon1T
bound and the expected DarkSide G2 bound but can be excluded or discovered by
LZ.

##### BP2.

For this point, we choose $M_{H^{\pm}}=500~{}{\rm GeV}$ and $M_{N_{R}}=200$
GeV. The Yukawa-type couplings are chosen using the same hierarchy as BP1 but
with relatively different values, i.e., $Y_{\mu N}=1.6$, $Y_{\tau N}=5\times
10^{-1}$ and $Y_{eN}=5\times 10^{-4}$. This leads to the following branching
ratios ${\rm BR}(H^{\pm}\to\mu^{\pm}N_{R})\simeq 91\%$, ${\rm
BR}(H^{\pm}\to\tau^{\pm}N_{R})\simeq 9\%$ and ${\rm BR}(H^{\pm}\to
e^{\pm}N_{R})\simeq 0\%$. The charged singlet is narrow in this case as
$\Gamma_{H^{\pm}}/M_{H^{\pm}}\simeq 0.04$. The CLFV decays of charged leptons
exhibits similar features as in BP1 with the exception that ${\rm BRs}$ of
$\tau\to\mu N_{R}$ and $\mu\to e\gamma$ can be probed in the future
experiments as they are slightly below the current bounds. The spin-
independent DM-nucleon cross section can be tested by the DarkSide G2
experiment.

##### BP3.

For this point, we choose the following values of the particle masses:
$M_{H^{\pm}}=600$ GeV and $M_{N_{R}}=598$ GeV and therefore a small mass
splitting of $2$ GeV. We choose the following values for the Yukawa-type
couplings: $\\{Y_{\mu N},Y_{\tau N},Y_{eN}\\}=\\{1,0.5,10^{-3}\\}$ which leads
to the following branching fractions: ${\rm
BR}(H^{\pm}\to\mu^{\pm}N_{R})\simeq 90\%$, ${\rm
BR}(H^{\pm}\to\tau^{\pm}N_{R})\simeq 10\%$ and ${\rm BR}(H^{\pm}\to
e^{\pm}N_{R})\simeq 0\%$. On the other hand, the branching ratios of CLFV
decays of $\mu\to e\gamma$ and $\tau\to\mu\gamma$ can be tested in future
experiments. Since the mass splitting is equal to $2$ GeV, the most active
component in the calculation of the relic density comes from coannihilation-
based freezout and therefore the choice of $\lambda_{3}$ is pivotal in this
case. We found that for this BP, the relic density of the $N_{R}$ is below
$2\%$ of the total observed DM relic density. Finally, this BP is not
sensitive to the direct detection experiments and the cross section is above
the neutrino floor.

##### BP4.

Here, we choose relatively heavy DM and charged singlet scalar;
$M_{N_{R}}=1000$ GeV and $M_{H^{\pm}}=1500$ GeV. The Yukawa-type couplings are
chosen such that $Y_{\mu N}=Y_{\tau N}=2\gg Y_{eN}=5\times 10^{-3}$. With this
choice, one gets: ${\rm BR}(H^{\pm}\to\mu^{\pm}N_{R})\simeq{\rm
BR}(H^{\pm}\to\tau^{\pm}N_{R})\simeq 50\%$ while ${\rm BR}(H^{\pm}\to
e^{\pm}N_{R})$. Similar features to BP1 and BP2 are observed for CLFV and DM
phenomenology.

##  V Production of dark matter at muon colliders

###  V.1 Total cross sections

In this section, we discuss the general features of DM production at muon
colliders999The cross sections for both DM and charged scalar production at
muon colliders are computed at Leading order using Madgraph_aMC@NLO
Alwall:2014hca  with a UFO model file Degrande:2011ua  that can be found in
the FeynRules model database
https://feynrules.irmp.ucl.ac.be/wiki/ModelDatabaseMainPage.. In this model,
DM can be produced through a variety of processes:

  * •

DM production in association with one SM particle dubbed as mono-X. Given the
nature of the interaction Lagrangian and the fact that the initial state has a
zero total electric charge, DM can only produced in association with one
neutral boson. Therefore, we have mono-$\gamma$, mono-$Z$ and mono-Higgs (a
full analysis of these channels will be done in future work Jueid:2023txs ).

  * •

DM production in association with two SM particles. For this category, we have
seven different processes. The rates of those processes are slightly smaller
than the mono-X production channels. However, these processes have smaller
backgrounds (a full analysis of these channels will be done in future work
Jueid:2023tzy ).

  * •

DM production in association with three SM particles. The rates of the
$N_{R}N_{R}$ in association with three SM particles are even smaller than the
other two categories. The signal-to-background optimisation for these channels
are even more complicated while the backgrounds, on the other hand, are
extremely small.

Figure 6: Production cross section of $N_{R}N_{R}+X$ as a function of the
center-of-mass energy ($\sqrt{s_{\mu\mu}}$) for the benchmark points BP1 (left
upper panel), BP2 (right upper panel), BP3 (left lower panel) and BP4 (left
lower panel). For each pane, we show the production cross section for
$N_{R}N_{R}$ plus one SM particle, plus two SM particles and in association
with three SM particles.

$\sigma\times{\rm BR}$ [fb] (number of events) Dominant backgrounds
$\sqrt{s_{\mu\mu}}$ [TeV] $3$ $10$ $30$ $N_{R}N_{R}\gamma$ $1.11\times
10^{3}~{}(1.11\times 10^{6})$ $1.80\times 10^{2}~{}(1.80\times 10^{6})$
$2.65\times 10^{1}~{}(2.65\times 10^{6})$
$\nu\bar{\nu}+\gamma,2\nu\bar{\nu}+\gamma$ $1.13\times 10^{2}~{}(1.13\times
10^{5})$ $1.88\times 10^{1}~{}(1.88\times 10^{5})$ $2.83\times
10^{0}~{}(2.83\times 10^{5})$ $1.18\times 10^{1}~{}(1.18\times 10^{3})$
$2.65\times 10^{0}~{}(2.65\times 10^{4})$ $0.41\times 10^{0}~{}(4.10\times
10^{4})$ $3.92\times 10^{1}~{}(3.95\times 10^{4})$ $3.20\times
10^{1}~{}(3.20\times 10^{5})$ $5.94\times 10^{0}~{}(5.94\times 10^{5})$
$N_{R}N_{R}Z(\to\ell\ell)$ $1.68\times 10^{1}~{}(1.68\times 10^{4})$
$4.44\times 10^{0}~{}(4.44\times 10^{4})$ $0.91\times 10^{0}~{}(9.10\times
10^{4})$ $1.62\times 10^{0}~{}(1.62\times 10^{3})$ $0.46\times
10^{0}~{}(4.58\times 10^{3})$ $9.39\times 10^{-2}~{}(9.39\times 10^{3})$
$\gamma/Z(\to\ell\ell)+\nu\bar{\nu}$ $0.13\times 10^{0}~{}(0.13\times 10^{3})$
$0.58\times 10^{-1}~{}(0.58\times 10^{3})$ $1.30\times 10^{-2}~{}(1.30\times
10^{3})$ $W(\to\ell\nu_{\ell})W(\to\ell\nu_{\ell})$ $0.28\times
10^{0}~{}(0.28\times 10^{3})$ $0.61\times 10^{0}~{}(0.61\times 10^{4})$
$0.17\times 10^{0}~{}(1.70\times 10^{4})$ $N_{R}N_{R}Z(\to q\bar{q})$
$1.59\times 10^{2}~{}(1.59\times 10^{5})$ $4.20\times 10^{1}~{}(4.20\times
10^{5})$ $8.61\times 10^{0}~{}(8.61\times 10^{5})$ $1.53\times
10^{1}~{}(1.53\times 10^{4})$ $4.33\times 10^{0}~{}(4.33\times 10^{4})$
$0.89\times 10^{0}~{}(8.89\times 10^{4})$ $\gamma/Z(\to
q\bar{q})+\nu\bar{\nu},H_{\rm SM}(\to b\bar{b})+\nu\bar{\nu}$ $1.26\times
10^{0}~{}(1.26\times 10^{3})$ $0.55\times 10^{0}~{}(5.54\times 10^{3})$
$0.12\times 10^{0}~{}(1.23\times 10^{4})$ $W(\to\ell\nu_{\ell})W(\to
q\bar{q}),t\bar{t}$ $2.67\times 10^{0}~{}(2.67\times 10^{3})$ $5.73\times
10^{0}~{}(5.73\times 10^{4})$ $1.57\times 10^{0}~{}(1.57\times 10^{5})$
$N_{R}N_{R}H_{\rm SM}(\to b\bar{b})$ $2.05\times 10^{1}~{}(2.05\times 10^{4})$
$1.02\times 10^{0}~{}(1.02\times 10^{4})$ $3.67\times 10^{-2}~{}(3.67\times
10^{3})$ $5.83\times 10^{0}~{}(5.83\times 10^{3})$ $0.31\times
10^{0}~{}(0.31\times 10^{4})$ $1.12\times 10^{-2}~{}(1.12\times 10^{3})$
$H_{\rm SM}(\to b\bar{b})Z(\to\nu\bar{\nu}),H_{\rm
SM}\nu_{\mu}\bar{\nu}_{\mu}$ $0.47\times 10^{0}~{}(0.47\times 10^{3})$
$0.47\times 10^{-1}~{}(0.47\times 10^{3})$ $1.81\times 10^{-3}~{}(1.81\times
10^{2})$ $t\bar{t},Z(\to\nu\bar{\nu})Z(\to q\bar{q})$ $0.11\times
10^{0}~{}(0.11\times 10^{3})$ $0.21\times 10^{0}~{}(0.21\times 10^{4})$
$1.47\times 10^{-2}~{}(1.47\times 10^{3})$

Table 3:  The total cross sections times the branching ratio
($\sigma\times{\rm BR}$) and the expected number of signal events for the
$N_{R}N_{R}$ production in association with $\gamma$, $Z(\to\ell\ell)$, $Z(\to
q\bar{q})$, and $H_{\rm SM}(\to b\bar{b})$. We consider three representative
center-of-mass energies of $3$, $10$ and $30$ TeV. For each process, we show
four entries that correspond to the benchmark points considered in this study
along with the associated background contributions. Here $\ell$ refers to
either an electron or a muon.

In figure 6 we show the total cross sections for DM production in $\mu\mu$
collisions as function of the center-of-mass energy ($\sqrt{s_{\mu\mu}}$) for
the four benchmark points defined in table 2. Starting with mono–X processes,
it is clear that the mono-$\gamma$ channel has the highest rate which varies
from $\simeq 1$ pb for $\sqrt{s_{\mu\mu}}=3$ TeV to about $80$ fb for
$\sqrt{s_{\mu\mu}}=30$ TeV in BP1101010Note that for mono–$\gamma$, we have
applied some generator-level cuts by requiring that $p_{T}^{\gamma}>25$ GeV
and $|\eta^{\gamma}|<2.5$.. Mono-$Z$ production has the second highest cross
sections which varies between $200$ fb for $\sqrt{s_{\mu\mu}}=3$ TeV and about
$2$ fb for $\sqrt{s_{\mu\mu}}=30$ TeV. Finally, mono-Higgs production has the
lowest rates among all the mono–X processes with cross section approaching
$63$ fb for $\sqrt{s_{\mu\mu}}=3$. The rates for mono–X decrease by about a
factor of 10 for BP2, by a factor of 100 for BP3 and by a factor of 10 for
BP4. Notice that the decrease in the production cross sections is not due to
the DM mass but also to the change in the value of $Y_{\mu N}$ since the total
rates are proportional to $Y_{\mu N}^{4}$. An exception to this rule is in the
mono-Higgs production cross section which decrease by factors of $6$–$200$
since it scales as $\lambda_{3}^{2}Y_{\mu N}^{4}$.

The rates of the production of DM in association with two SM particles are
shown in figure 6. We can see that, as expected, they are suppressed as
compared to the case of mono–X channels. The process with the highest is
$N_{R}N_{R}+\gamma\gamma$ whose cross section is between $50$ fb and $2$ fb.
This process is followed by $N_{R}N_{R}\gamma Z$ and $N_{R}N_{R}W^{+}W^{-}$
whose cross sections are slightly smaller. An interesting process is the
production of DM in association with two SM Higgs bosons whose cross sections
is about $1$–$3$ fb depending on the center-of-mass energy. We note that the
rates of these processes decrease as the DM mass increase, i.e. by a factor of
10 for BP2.

Finally, the production cross sections of DM in association with three SM
particles are shown in figure 6 for BP1–BP4. It is clear that these rates are
suppressed as compared to those of the DM production in association with one
SM particle and two SM particles respectively. The maximum being about $1$ fb
for $N_{R}N_{R}W^{+}W^{-}\gamma$ and $N_{R}N_{R}\gamma\gamma\gamma$ at
$\sqrt{s_{\mu\mu}}=3$ TeV. We note that the dependence on $\sqrt{s_{\mu\mu}}$
of the cross sections for the production of $N_{R}N_{R}$ in association with
three SM particles is not a strong as in case of other processes. Despite the
smallness of these cross section, these processes may have a high sensitivity
reach due to the smallness of the associated backgrounds.

$\sigma\times{\rm BR}$ [fb] (number of events) Dominant backgrounds
$\sqrt{s_{\mu\mu}}$ [TeV] $3$ $10$ $30$ $N_{R}N_{R}\gamma\gamma$ $4.97\times
10^{1}~{}(4.97\times 10^{4})$ $1.23\times 10^{1}~{}(1.23\times 10^{5})$
$2.38\times 10^{0}~{}(2.38\times 10^{5})$ $4.93\times 10^{1}~{}(4.93\times
10^{3})$ $1.28\times 10^{0}~{}(1.28\times 10^{4})$ $0.25\times
10^{0}~{}(2.53\times 10^{4})$
$\nu\bar{\nu}+\gamma\gamma,2\nu\bar{\nu}+\gamma\gamma$ $0.43\times
10^{0}~{}(0.43\times 10^{3})$ $0.17\times 10^{0}~{}(1.73\times 10^{3})$
$0.36\times 10^{-1}~{}(3.64\times 10^{3})$ $\mu\mu/VV\to H_{\rm
SM}(\to\gamma\gamma)+\nu\bar{\nu}$ $1.00\times 10^{0}~{}(1.00\times 10^{3})$
$1.87\times 10^{0}~{}(1.87\times 10^{4})$ $0.48\times 10^{0}~{}(4.85\times
10^{4})$ $N_{R}N_{R}\gamma Z(\to\ell\ell)$ $1.24\times 10^{0}~{}(1.24\times
10^{3})$ $0.49\times 10^{0}~{}(4.98\times 10^{3})$ $1.29\times
10^{-1}~{}(1.29\times 10^{4})$ $1.76\times 10^{0}~{}(1.76\times 10^{3})$
$0.76\times 10^{0}~{}(7.64\times 10^{3})$ $0.20\times 10^{0}~{}(2.02\times
10^{4})$ $\gamma Z(\to\ell\ell)+\nu\bar{\nu}$ $0.12\times 10^{0}~{}(1.23\times
10^{2})$ $9.50\times 10^{-2}~{}(9.50\times 10^{2})$ $2.80\times
10^{-2}~{}(2.80\times 10^{3})$ $\gamma
W(\to\ell\nu_{\ell})W(\to\ell\nu_{\ell})$ $0.18\times 10^{0}~{}(1.79\times
10^{2})$ $9.05\times 10^{-1}~{}(9.05\times 10^{3})$ $3.46\times
10^{-1}~{}(3.46\times 10^{4})$ $N_{R}N_{R}Z(\to\ell\ell)Z(\to\ell\ell)$
$3.53\times 10^{-2}~{}(3.53\times 10^{1})$ $2.29\times 10^{-2}~{}(2.29\times
10^{2})$ $7.21\times 10^{-3}~{}(7.21\times 10^{2})$ $0.98\times
10^{0}~{}(9.80\times 10^{2})$ $0.75\times 10^{0}~{}(7.54\times 10^{3})$
$0.24\times 10^{0}~{}(2.42\times 10^{4})$
$Z(\to\ell\ell)Z(\to\ell\ell)Z(\to\nu\bar{\nu})$ $3.23\times
10^{-2}~{}(3.23\times 10^{1})$ $7.87\times 10^{-2}~{}(7.87\times 10^{2})$
$3.08\times 10^{-2}~{}(3.08\times 10^{3})$ $H_{\rm SM}(\to ZZ^{*})Z,H_{\rm
SM}(\to WW^{*})Z$ $7.50\times 10^{-3}~{}(7.50\times 10^{0})$ $0.39\times
10^{0}~{}(3.89\times 10^{3})$ $0.30\times 10^{0}~{}(3.02\times 10^{4})$
$W(\to\ell\nu)W(\to\ell\nu)Z(\to\ell\ell)$ $N_{R}N_{R}V(\to q\bar{q})V(\to
q\bar{q})$ $1.05\times 10^{1}~{}(1.05\times 10^{4})$ $6.57\times
10^{0}~{}(6.57\times 10^{4})$ $2.02\times 10^{0}~{}(2.02\times 10^{5})$
$2.76\times 10^{0}~{}(2.76\times 10^{3})$ $2.08\times 10^{0}~{}(2.08\times
10^{4})$ $0.65\times 10^{0}~{}(6.57\times 10^{4})$ $t\bar{t},V(\to
q\bar{q})V(\to q\bar{q})\nu\bar{\nu}$ $8.90\times 10^{-2}~{}(8.90\times
10^{1})$ $2.15\times 10^{-1}~{}(2.15\times 10^{3})$ $8.30\times
10^{-2}~{}(8.30\times 10^{3})$ $H_{\rm SM}(\to gg)H_{\rm SM}(\to
gg)\nu\bar{\nu}$ $1.30\times 10^{-2}~{}(1.30\times 10^{1})$ $9.74\times
10^{-1}~{}(9.74\times 10^{3})$ $7.96\times 10^{-1}~{}(7.96\times 10^{4})$
$H_{\rm SM}(\to b\bar{b})H_{\rm SM}(\to b\bar{b})\nu\bar{\nu}$
$N_{R}N_{R}H_{\rm SM}(\to b\bar{b})H_{\rm SM}(\to b\bar{b})$ $1.21\times
10^{0}~{}(1.21\times 10^{3})$ $1.12\times 10^{0}~{}(1.12\times 10^{4})$
$3.77\times 10^{-1}~{}(3.77\times 10^{4})$ $3.95\times 10^{-1}~{}(3.95\times
10^{2})$ $5.29\times 10^{-1}~{}(5.29\times 10^{3})$ $1.88\times
10^{-1}~{}(1.88\times 10^{4})$ $t\bar{t},V(\to q\bar{q})V(\to
q\bar{q})\nu\bar{\nu}$ $1.22\times 10^{-2}~{}(1.22\times 10^{1})$ $5.32\times
10^{-2}~{}(5.32\times 10^{2})$ $2.36\times 10^{-2}~{}(2.36\times 10^{3})$
$H_{\rm SM}(\to gg)H_{\rm SM}(\to gg)\nu\bar{\nu}$ $1.40\times
10^{-3}~{}(1.40\times 10^{0})$ $2.49\times 10^{-1}~{}(2.49\times 10^{3})$
$2.27\times 10^{-1}~{}(2.27\times 10^{4})$ $H_{\rm SM}(\to b\bar{b})H_{\rm
SM}(\to b\bar{b})\nu\bar{\nu}$

Table 4:  Same as in table 3 but $N_{R}N_{R}$ in association with
$\gamma\gamma$, $\gamma Z(\to\ell\ell)$, $Z(\to\ell\ell)Z(\to\ell\ell)$,
$V(\to q\bar{q})V(\to q\bar{q}$ and $H_{\rm SM}(\to b\bar{b})H_{\rm SM}(\to
b\bar{b})$. Here $V$ refers to either $W$ or $Z$.

###  V.2 Expected event yields and dominant backgrounds

After discussing the total cross sections for all the possible production
channels of dark matter at muon colliders, it is instructive to discuss both
the total expected number of events for specific decay channels of the SM
particles and the associated backgrounds. In this subsection, we focus on two
categories of DM production channels (i) DM production in association with one
SM particle where we consider four processes: $N_{R}N_{R}\gamma$,
$N_{R}N_{R}Z(\to\ell\ell)$, $N_{R}N_{R}Z(\to q\bar{q})$ and $N_{R}N_{R}H_{\rm
SM}(\to b\bar{b})$ and (ii) DM production in association with two SM particles
where we consider five processes: $N_{R}N_{R}\gamma\gamma$, $N_{R}N_{R}\gamma
Z(\to\ell\ell)$, $N_{R}N_{R}Z(\to\ell\ell)Z(\to\ell\ell)$, $N_{R}N_{R}V(\to
q\bar{q})V(\to q\bar{q})$ and $N_{R}N_{R}H_{\rm SM}(\to b\bar{b})H_{\rm
SM}(\to b\bar{b})$. The results are shown in tables 3 and 4. The discussion
will be restricted for the following center-of-mass energies and integrated
luminosities

$\displaystyle\sqrt{s_{\mu\mu}}$ $\displaystyle=$ $\displaystyle
3,~{}10,~{}{\rm and}~{}30~{}{\rm TeV}$ $\displaystyle\int{\rm d}t{\cal L}$
$\displaystyle=$ $\displaystyle 1,~{}10,~{}{\rm and}~{}100~{}{\rm ab}^{-1},$
(34)

where we follow ref. Han:2020uak  assuming that the luminosity has a linear
scaling with the center-of-mass energy. The expected number of events for
mono–X processes is calculated using the following equation

$\displaystyle{\cal N}=\sigma_{N_{R}N_{R}X}\times{\rm BR}_{X\to
x_{1}x_{2}}\times\int{\rm d}t{\cal L}.$ (35)

For the production of DM in association with two SM particles, we have

$\displaystyle{\cal N}=\sigma_{N_{R}N_{R}XY}\times{\rm BR}_{X\to
x_{1}x_{2}}\times{\rm BR}_{Y\to y_{1}y_{2}}\times\int{\rm d}t{\cal L}.$ (36)

####  V.2.1 $\mu^{+}\mu^{-}\to N_{R}N_{R}+X$

##### $N_{R}N_{R}\gamma.$

This process leads to the final state comprising of a highly energetic photon
and a large missing transverse energy ($E_{T}^{\rm miss}$). In addition one
could have a few additional charged leptons, or photons that are emitted from
the radiation of either the initial-state muons or the final-state photon. The
dominant backgrounds for this signal process are the production of two or four
neutrinos in association with a photon. The production of two neutrinos
proceeds via muon-muon annihiliation – $\mu^{+}\mu^{-}\to
Z(\to\nu\bar{\nu})\gamma$ – and VBF – $VV\to Z(\to\nu\bar{\nu})\gamma$ – with
cross sections varying from $2.98$ pb for $\sqrt{s_{\mu\mu}}=3$ TeV to $3.27$
pb for $\sqrt{s_{\mu\mu}}=30$ TeV. The production of four neutrinos in
association with hard photon has an extremely cross section with the maximum
being $1.5$ fb for $\sqrt{s_{\mu\mu}}=30$ TeV. It is worth noting from table 3
that the signal significance can easily $5$111111The signal significance is
defined as ${\cal S}/\sqrt{{\cal B}}$ with ${\cal S}$ is the number of signal
events and ${\cal B}$ is the number of background events.. For the other
benchmark points, a more detailed selection is required to reach a signal
significance of $5$ if one can achieve an acceptance times efficiency
($A\times\epsilon$) of about $15\%$ for the signal in the signal region while
the background is having $A\times\epsilon$ of about $\mathcal{O}(10^{-3})$.

Figure 7: Production cross section of $H^{\pm}H^{\mp}+X$ as a function of the
center-of-mass energy ($\sqrt{s_{\mu\mu}}$) for the benchmark points BP1 (left
upper panel), BP2 (right upper panel), BP3 (left lower panel) and BP4 (left
lower panel). For each pane, we show the production cross section for
$H^{\pm}H^{\mp}$ plus one SM particle, plus two SM particles and in
association with three SM particles.

##### $N_{R}N_{R}Z(\to\ell\ell).$

This process leads to a very clean final state containing two same-flavour
opposite-sign (SFOS) charged leptons from the decay of the $Z$-boson in
association with large missing energy. The dominant backgrounds are found to
be the production of two $Z$-bosons with one decaying two charged leptons and
the other decaying invisibly. We note that there is another background
originated from the production of two $W$-bosons both decaying leptonically
which can significantly be reduced using the requirement of two SFOS leptons
whose invariant mass is close to the $Z$-boson mass. The cross sections for
the $ZZ$ production varies from $0.4$ fb to $26$ fb for the muon-muon
annihilation (decreases while the center-of-mass energy increases) and from
$56$ fb to about $430$ fb (increases with the center-of-mass energy). On the
other hand, the cross section for $WW$ production is larger and varies between
$8.5$ fb and $466$ fb in the muon annihilation channel and between $150$ fb
and $858$ fb in the VBF channels. The expected number of events for this
signal process is about ${\cal O}(10^{3}$–$10^{4})$. Given the differences in
the topology of the signal and backgrounds, it is easy to achieve a
significance of $5\sigma$ by suitable event selection.

#####  $N_{R}N_{R}Z(\to q\bar{q})$ and $N_{R}N_{R}H_{\rm SM}(\to b\bar{b}).$

This category of channels involves two hadronic jets in association with
missing energy. For the case of the $Z$–boson, the main decay channel is into
$q\bar{q};q=u,d,s,c,b$ with ${\rm BR}(Z\to q\bar{q})=69.911\%$
ParticleDataGroup:2022pth . For the SM Higgs boson, the main decay is into
$b\bar{b}$ with ${\rm BR}(H_{\rm SM}\to b\bar{b})=57\%$. The dominant
backgrounds to these signal processes come from $q\bar{q}$ production in
association with two neutrinos, SM Higgs boson production, $t\bar{t}$
production with one top quark decaying leptonically and the other top decaying
hadronically and $WW$ production where one $W$-boson decays leptonically and
the other decaying hadronically. In the last two backgrounds, one needs that
the charged lepton escapes the detection volume. Since the both the
hadronically decaying $Z$– and Higgs-bosons are accompanied with very large
missing energy, their decays are not always resolved as two well separated two
jets but rather a fat jet with specific characteristics. We expect for these
channels a decent statistics and the backgrounds are under control with
suitable selection.

####  V.2.2 $\mu^{+}\mu^{-}\to N_{R}N_{R}+XY$

#####  $N_{R}N_{R}\gamma\gamma$.

In this model, there is a possibility to produce DM pairs in association with
two hard photons. The expected final state would consist of two hard photons
in addition to large missing energy. Contrary to mono–$\gamma$ channel, this
process does not have large backgrounds where the main backgrounds are the
production of two photons in association with two neutrinos: non-resonant and
resonant (through the decay of the SM Higgs boson). The resonant backgrounds
can be easily suppressed via suitable requirements on the invariant mass of
the diphoton system, i.e. removing photons that are within the SM Higgs mass
window. For all the benchmark points we expect a decent statistics for the
signal events, i.e. of about ${\cal O}(10^{3}$–$10^{5})$. This would imply
that this process would one of the golden modes to probe DM at muon colliders
which will be studied in details in a future work.

#####  $N_{R}N_{R}\gamma Z(\to\ell\ell)$.

This is also one of the unique processes to probe DM at muon colliders. The
final state consists of one hard photons, two charged leptons and large
missing transverse energy. The associated background is manageable since it
consists of the production of one photon and one or two gauge bosons. The
expected number of signal events is quite large as well, i.e. of about ${\cal
O}(10^{2}$–$10^{4})$.

#####  $N_{R}N_{R}Z(\to\ell\ell)Z(\to\ell\ell)$.

This one of the most cleanest final states that can be used to probe DM at
muon colliders. The signature consists of four charged leptons in association
with missing energy. The corresponding is even smaller than for the other
signal processes. We note that enough statistics can be only achieved at
$\sqrt{s_{\mu\mu}}=30$ TeV where we expect about ${\cal O}(10^{2}$–$10^{4}$)
events. The major backgrounds arise from the production of three gauge bosons
or from the production of the SM Higgs boson decaying into $VV^{*},V=W,Z$ in
association with one or two gauge bosons.

#####  $N_{R}N_{R}V(\to q\bar{q})V(\to q\bar{q})$ and $N_{R}N_{R}H_{\rm
SM}(\to b\bar{b})H_{\rm SM}(\to b\bar{b})$.

The production of two gauge bosons or two SM Higgs bosons in association with
DM pairs lead to purely hadronic final states (either four resolved jets or
two fat jets) in association with large missing energy. The dominant
backgrounds for these signal processes consist of the production of two SM
neutrinos in association with two gauge bosons, two SM Higgs bosons, or one
Higgs boson and one gauge boson decaying hadronically. This process will be
studied in great detail in a future work.

##  VI Production of charged scalars at muon colliders

In this section, we discuss the production of charged scalar pairs at muon
colliders. Similarly to the production of DM, charged scalars can be produced
either in association with one SM particle, with two SM particles or three SM
particles. In addition, we could have the production of charged scalar pairs
with non SM particles ($H^{\pm}H^{\pm}$) or the production of four charged
scalars. An interesting feature about the production of charged scalars is
that the appearance of at least two leptons in association with missing energy
in addition to the decay products of the SM particles. For example, the
production of charged scalar pairs in association with a SM Higgs boson would
lead to two hard charged leptons, missing energy and two $b$-tagged jets (or
one fat jet). On the other hand, the charged scalar production receives
contributions from VBF thanks to their couplings to $\gamma/Z$. The results of
the production cross sections for the different processes involving charged
scalars as a function of the center-of-mass energy are shown in figure 7.
Below, we list the possible production channels for the charged scalars:

#####  $\mu\mu\to H^{\pm}H^{\mp}/H^{\pm}H^{\mp}H^{\pm}H^{\mp}$.

These processes lead to signatures of either two charged lepton and MET or
four charged leptons and MET. Charged scalar pair production proceeds through
either $s$–channel diagrams with the exchange $\gamma/Z$–bosons or $t$–channel
diagram with the exchange of the Majorana DM. The cross section for charged
scalar pair production ranges from about $10^{4}$ fb to about $10^{1}$ fb. It
is worth noting that for the benchmark point BP3 has the smallest cross
section due to the tiny mass splitting of about $2$ GeV between the charged
scalar and the DM candidate. In all the case, the number of events for this
process is quite large. The cross section for charged scalar pair production
has a $1/\sqrt{s_{\mu\mu}}$ scaling. The cross section for the production of
four charged scalars is smaller as expected due to phase suppression. It is
however quite decent as can be seen in fig. 7 and ranges from $10^{-2}$ fb to
$10^{2}$ depending on the benchmark scenarios. The most notable signatures are
$4$ muons plus MET (BP1, BP2, BP3) and $2$ muons and $2$ tau leptons (BP4).
These two channels will be studied in great detail in a future work
Jueid:2023xyz .

#####  $\mu\mu\to H^{\pm}H^{\mp}+X$.

In this case, we have three production channels: $H^{+}H^{-}\gamma$,
$H^{+}H^{-}Z$ and $H^{+}H^{-}H_{\rm SM}$. There are three contributions to
$H^{+}H^{-}\gamma$: $s$–channel contributions through $\gamma/Z$ with the
photon being emitted from the $H^{+}H^{-}$ vertex and $t$–channel contribution
through the exchange of $N_{R}$. The final-state signature for this process
consists of two charged leptons in association with one hard photon (the
kinematics is quite different from $N_{R}N_{R}\gamma Z$ production). We can
see in fig. 7 that the cross section ranges from $10^{1}$ fb to $10^{3}$ fb
depending on the center-of-mass energy and the benchmark point. Secondly, we
can have the production of charged scalar pairs in association with one
$Z$-boson which would lead to very rich signatures: $2\ell$ \+ MET, $2\ell$ \+
2 jets + MET or $4\ell$ \+ MET. The cross sections for these processes are
shown in fig. 7 where it is clear that the rates are quite important from
$10^{0}$ to $10^{2}$ fb. Finally, the charged scalar pairs can be produced in
association with a SM Higgs boson. The rates for this interesting channel are
also quite important and range between $10^{0}$ fb and $10^{2}$ fb.

#####  $\mu\mu\to H^{\pm}H^{\mp}+XY$.

For this category we have seven production channels: $H^{+}H^{-}\gamma\gamma$,
$H^{+}H^{-}\gamma Z$, $H^{+}H^{-}ZZ$, $H^{+}H^{-}W^{+}W^{-}$,
$H^{+}H^{-}H_{\rm SM}H_{\rm SM}$, $H^{+}H^{-}H_{\rm SM}Z$ and
$H^{+}H^{-}t\bar{t}$. The rates for this channels are quite smaller but still
at the noticeable level, i.e. from $10^{-2}$ fb to $10^{2}$ fb depending on
the center-of-mass energy and the benchmark point. It is worth noting that the
production of charged scalar pairs in association with two SM particles leads
to even much more richer signatures with very small backgrounds, it.e. $6$
leptons plus MET, $4$ leptons plus 4 jets plus MET and so on.

#####  $\mu\mu\to H^{\pm}H^{\mp}+XYZ$.

This is the most complicated category of processes where we can have 16
processes with many more final-state signatures. The rates for these processes
are much more smaller with the maximum being about $3$ fb for
$H^{+}H^{-}\gamma H_{\rm SM}H_{\rm SM}$ and $H^{+}H^{-}\gamma\gamma H_{\rm
SM}$ at $\sqrt{s_{\mu\mu}}=3$ TeV.

Figure 8: The production cross sections for $H^{+}H^{-}$, $H^{+}H^{-}\gamma$
and $H^{+}H^{-}H_{\rm SM}$ through VBF as a function of the center-of-mass
energy ($\sqrt{s_{\mu\mu}}$) for BP1 (solid), BP2 (dashed), BP3 (dotted) and
BP4 (dash-dotted).

We close this section by a brief discussion of the contribution of VBF to the
production of charged scalars in this model. As mentioned earlier the charged
scalar couples to the photon and the $Z$-boson and therefore may receive pure
gauge VBF contributions to the total production cross section. In this model,
we can have the production of charged scalars through $\gamma\gamma
H^{+}H^{-}$, $\gamma ZH^{+}H^{-}$, $ZZH^{+}H^{-}$, $ZZ\to H_{\rm SM}\to
H^{+}H^{-}$ and $W^{+}W^{-}\to\gamma/Z\to H^{+}H^{-}$ vertices. We take
examples of production of $H^{+}H^{-}$, $H^{+}H^{-}\gamma$ and
$H^{+}H^{-}H_{\rm SM}$ and show the corresponding results for the four
benchmark points in fig. 8. We can see that the cross sections increase with
center-of-mass energy but do not go above $2$ fb for $H^{+}H^{-}$ in BP2.
Therefore, the muon annihilation channels are the most important in our model
thanks to the $Y_{\mu N}^{4}$ dependence of the cross section.

##  VII Conclusions

In this work we have studied the production of DM and charged scalars at high
energy muon colliders within the minimal lepton portal DM model. The model
consists of extending the SM with two $SU(2)_{L}$ gauge singlets: a charged
singlet scalar and a neutral right-handed fermion (or equivalently a Majorana
fermion). We first discussed in details the phenomenology of the model at the
LHC and the corresponding constraints from direct detection, relic density
measurement, and lepton flavour violating decays of charged leptons and the SM
Higgs boson. Then we have selected a few benchmark points that define some
phenomenologically viable scenarios and which can be tested at future muon
colliders. For these benchmark points, we have calculated the cross sections
for the production of DM in association with SM particles and of charged
scalars of the models in association with SM particles as a function of the
center-of-mass energy. For DM production in association with SM particles, we
have studied the total rates of 26 possible channels for the benchmark points
considered in this study. Furthermore, we studied the total number of events
and the associated backgrounds for 9 prominent channels and found that they
are very important for the discovery DM at muon colliders for masses up to
$\sim 1$ TeV. We furthermore analysed the production of charged scalar
production in association with SM particles (about 28 channels). The potential
discovery for DM through charged scalar production at muon colliders is also
as interesting as for direct production of DM. Further investigations of this
model at muon colliders are ongoing where a full signal-to-background
optimisation will be carried out for a number of selected channels.

## Acknowledgements

The work of AJ is supported by the Institute for Basic Science (IBS) under the
project code, IBS-R018-D1.

## References

  * (1) G. Jungman, M. Kamionkowski, and K. Griest, “Supersymmetric dark matter,” Phys. Rept. 267 (1996) 195–373, arXiv:hep-ph/9506380.
  * (2) L. Bergström, “Nonbaryonic dark matter: Observational evidence and detection methods,” Rept. Prog. Phys. 63 (2000) 793, arXiv:hep-ph/0002126.
  * (3) G. Bertone, D. Hooper, and J. Silk, “Particle dark matter: Evidence, candidates and constraints,” Phys. Rept. 405 (2005) 279–390, arXiv:hep-ph/0404175 [hep-ph].
  * (4) J. L. Feng, “Dark Matter Candidates from Particle Physics and Methods of Detection,” Ann. Rev. Astron. Astrophys. 48 (2010) 495–545, arXiv:1003.0904 [astro-ph.CO].
  * (5) Planck Collaboration, P. A. R. Ade et al., “Planck 2015 results. XIII. Cosmological parameters,” Astron. Astrophys. 594 (2016) A13, arXiv:1502.01589 [astro-ph.CO].
  * (6) M. Beltran, D. Hooper, E. W. Kolb, Z. A. C. Krusberg, and T. M. P. Tait, “Maverick dark matter at colliders,” JHEP 09 (2010) 037, arXiv:1002.4137 [hep-ph].
  * (7) L. M. Carpenter, A. Nelson, C. Shimmin, T. M. P. Tait, and D. Whiteson, “Collider searches for dark matter in events with a Z boson and missing energy,” Phys. Rev. D 87 (2013) no. 7, 074005, arXiv:1212.3352 [hep-ex].
  * (8) N. F. Bell, J. B. Dent, A. J. Galea, T. D. Jacques, L. M. Krauss, and T. J. Weiler, “Searching for Dark Matter at the LHC with a Mono-Z,” Phys. Rev. D 86 (2012) 096011, arXiv:1209.0231 [hep-ph].
  * (9) A. Berlin, T. Lin, and L.-T. Wang, “Mono-Higgs Detection of Dark Matter at the LHC,” JHEP 06 (2014) 078, arXiv:1402.7074 [hep-ph].
  * (10) ATLAS Collaboration, G. Aad et al., “Search for new phenomena in events with an energetic jet and missing transverse momentum in $pp$ collisions at $\sqrt{s}$ =13 TeV with the ATLAS detector,” Phys. Rev. D 103 (2021) no. 11, 112006, arXiv:2102.10874 [hep-ex].
  * (11) ATLAS Collaboration, G. Aad et al., “Search for dark matter produced in association with a Standard Model Higgs boson decaying into b-quarks using the full Run 2 dataset from the ATLAS detector,” JHEP 11 (2021) 209, arXiv:2108.13391 [hep-ex].
  * (12) ATLAS Collaboration, G. Aad et al., “Search for new phenomena in events with two opposite-charge leptons, jets and missing transverse momentum in pp collisions at $\sqrt{\mathrm{s}}$ = 13 TeV with the ATLAS detector,” JHEP 04 (2021) 165, arXiv:2102.01444 [hep-ex].
  * (13) CMS Collaboration, A. M. Sirunyan et al., “Search for dark matter particles produced in association with a Higgs boson in proton-proton collisions at $\sqrt{\mathrm{s}}$ = 13 TeV,” JHEP 03 (2020) 025, arXiv:1908.01713 [hep-ex].
  * (14) CMS Collaboration, A. M. Sirunyan et al., “Search for dark matter produced in association with a leptonically decaying Z boson in proton-proton collisions at $\sqrt{s}=$ 13 TeV,” Eur. Phys. J. C 81 (2021) no. 1, 13, arXiv:2008.04735 [hep-ex]. [Erratum: Eur.Phys.J.C 81, 333 (2021)].
  * (15) XENON Collaboration, E. Aprile et al., “Dark Matter Search Results from a One Ton-Year Exposure of XENON1T,” Phys. Rev. Lett. 121 (2018) no. 11, 111302, arXiv:1805.12562 [astro-ph.CO].
  * (16) PandaX-II Collaboration, X. Cui et al., “Dark Matter Results From 54-Ton-Day Exposure of PandaX-II Experiment,” Phys. Rev. Lett. 119 (2017) no. 18, 181302, arXiv:1708.06917 [astro-ph.CO].
  * (17) GAMBIT Collaboration, P. Athron et al., “Status of the scalar singlet dark matter model,” Eur. Phys. J. C 77 (2017) no. 8, 568, arXiv:1705.07931 [hep-ph].
  * (18) GAMBIT Collaboration, P. Athron et al., “A global fit of the MSSM with GAMBIT,” Eur. Phys. J. C 77 (2017) no. 12, 879, arXiv:1705.07917 [hep-ph].
  * (19) M. J. Baker et al., “The Coannihilation Codex,” JHEP 12 (2015) 120, arXiv:1510.03434 [hep-ph].
  * (20) S.-M. Choi, Y.-J. Kang, H. M. Lee, and T.-G. Ro, “Lepto-Quark Portal Dark Matter,” JHEP 10 (2018) 104, arXiv:1807.06547 [hep-ph].
  * (21) K. A. Mohan, D. Sengupta, T. M. P. Tait, B. Yan, and C. P. Yuan, “Direct detection and LHC constraints on a $t$-channel simplified model of Majorana dark matter at one loop,” JHEP 05 (2019) 115, arXiv:1903.05650 [hep-ph].
  * (22) G. Belanger et al., “Leptoquark manoeuvres in the dark: a simultaneous solution of the dark matter problem and the ${R}_{D^{\left(\ast\right)}}$ anomalies,” JHEP 02 (2022) 042, arXiv:2111.08027 [hep-ph].
  * (23) B. Belfatto, D. Buttazzo, C. Gross, P. Panci, A. Strumia, N. Vignaroli, L. Vittorio, and R. Watanabe, “Dark Matter abundance via thermal decays and leptoquark mediators,” JHEP 06 (2022) 084, arXiv:2111.14808 [hep-ph].
  * (24) J. Liu, B. Shuve, N. Weiner, and I. Yavin, “Looking for new charged states at the LHC: Signatures of Magnetic and Rayleigh Dark Matter,” JHEP 07 (2013) 144, arXiv:1303.4404 [hep-ph].
  * (25) Y. Bai and J. Berger, “Lepton Portal Dark Matter,” JHEP 08 (2014) 153, arXiv:1402.6696 [hep-ph].
  * (26) S. Chang, R. Edezhath, J. Hutchinson, and M. Luty, “Leptophilic Effective WIMPs,” Phys. Rev. D 90 (2014) no. 1, 015011, arXiv:1402.7358 [hep-ph].
  * (27) P. Agrawal, Z. Chacko, and C. B. Verhaaren, “Leptophilic Dark Matter and the Anomalous Magnetic Moment of the Muon,” JHEP 08 (2014) 147, arXiv:1402.7369 [hep-ph].
  * (28) M. Garny, A. Ibarra, and S. Vogl, “Signatures of Majorana dark matter with t-channel mediators,” Int. J. Mod. Phys. D 24 (2015) no. 07, 1530019, arXiv:1503.01500 [hep-ph].
  * (29) A. Jueid, S. Nasri, and R. Soualah, “Searching for GeV-scale Majorana Dark Matter: inter spem et metum,” JHEP 04 (2021) 012, arXiv:2006.01348 [hep-ph].
  * (30) S.-I. Horigome, T. Katayose, S. Matsumoto, and I. Saha, “Leptophilic fermion WIMP: Role of future lepton colliders,” Phys. Rev. D 104 (2021) no. 5, 055001, arXiv:2102.08645 [hep-ph].
  * (31) J. Liu, X.-P. Wang, and K.-P. Xie, “Searching for lepton portal dark matter with colliders and gravitational waves,” JHEP 06 (2021) 149, arXiv:2104.06421 [hep-ph].
  * (32) J. P. Delahaye, M. Diemoz, K. Long, B. Mansoulié, N. Pastrone, L. Rivkin, D. Schulte, A. Skrinsky, and A. Wulzer, “Muon Colliders,” arXiv:1901.06150 [physics.acc-ph].
  * (33) K. Long, D. Lucchesi, M. Palmer, N. Pastrone, D. Schulte, and V. Shiltsev, “Muon colliders to expand frontiers of particle physics,” Nature Phys. 17 (2021) no. 3, 289–292, arXiv:2007.15684 [physics.acc-ph].
  * (34) H. Al Ali et al., “The muon Smasher’s guide,” Rept. Prog. Phys. 85 (2022) no. 8, 084201, arXiv:2103.14043 [hep-ph].
  * (35) A. Costantini, F. De Lillo, F. Maltoni, L. Mantani, O. Mattelaer, R. Ruiz, and X. Zhao, “Vector boson fusion at multi-TeV muon colliders,” JHEP 09 (2020) 080, arXiv:2005.10289 [hep-ph].
  * (36) R. Ruiz, A. Costantini, F. Maltoni, and O. Mattelaer, “The Effective Vector Boson Approximation in high-energy muon collisions,” JHEP 06 (2022) 114, arXiv:2111.02442 [hep-ph].
  * (37) R. Capdevilla, D. Curtin, Y. Kahn, and G. Krnjaic, “Discovering the physics of $(g-2)_{\mu}$ at future muon colliders,” Phys. Rev. D 103 (2021) no. 7, 075028, arXiv:2006.16277 [hep-ph].
  * (38) M. Chiesa, F. Maltoni, L. Mantani, B. Mele, F. Piccinini, and X. Zhao, “Measuring the quartic Higgs self-coupling at a multi-TeV muon collider,” JHEP 09 (2020) 098, arXiv:2003.13628 [hep-ph].
  * (39) T. Han, Y. Ma, and K. Xie, “High energy leptonic collisions and electroweak parton distribution functions,” Phys. Rev. D 103 (2021) no. 3, L031301, arXiv:2007.14300 [hep-ph].
  * (40) T. Han, Z. Liu, L.-T. Wang, and X. Wang, “WIMPs at High Energy Muon Colliders,” Phys. Rev. D 103 (2021) no. 7, 075004, arXiv:2009.11287 [hep-ph].
  * (41) W. Yin and M. Yamaguchi, “Muon g-2 at a multi-TeV muon collider,” Phys. Rev. D 106 (2022) no. 3, 033007, arXiv:2012.03928 [hep-ph].
  * (42) G.-y. Huang, F. S. Queiroz, and W. Rodejohann, “Gauged $L_{\mu}{-}L_{\tau}$ at a muon collider,” Phys. Rev. D 103 (2021) no. 9, 095005, arXiv:2101.04956 [hep-ph].
  * (43) R. Capdevilla, D. Curtin, Y. Kahn, and G. Krnjaic, “No-lose theorem for discovering the new physics of (g-2)$\mu$ at muon colliders,” Phys. Rev. D 105 (2022) no. 1, 015028, arXiv:2101.10334 [hep-ph].
  * (44) R. Capdevilla, F. Meloni, R. Simoniello, and J. Zurita, “Hunting wino and higgsino dark matter at the muon collider with disappearing tracks,” JHEP 06 (2021) 133, arXiv:2102.11292 [hep-ph].
  * (45) P. Asadi, R. Capdevilla, C. Cesarotti, and S. Homiller, “Searching for leptoquarks at future muon colliders,” JHEP 10 (2021) 182, arXiv:2104.05720 [hep-ph].
  * (46) M. Casarsa, M. Fabbrichesi, and E. Gabrielli, “Monochromatic single photon events at the muon collider,” Phys. Rev. D 105 (2022) no. 7, 075008, arXiv:2111.13220 [hep-ph].
  * (47) W. Liu, K.-P. Xie, and Z. Yi, “Testing leptogenesis at the LHC and future muon colliders: A Z’ scenario,” Phys. Rev. D 105 (2022) no. 9, 095034, arXiv:2109.15087 [hep-ph].
  * (48) T. Han, S. Li, S. Su, W. Su, and Y. Wu, “Heavy Higgs bosons in 2HDM at a muon collider,” Phys. Rev. D 104 (2021) no. 5, 055029, arXiv:2102.08386 [hep-ph].
  * (49) T. Han, Y. Ma, and K. Xie, “Quark and gluon contents of a lepton at high energies,” JHEP 02 (2022) 154, arXiv:2103.09844 [hep-ph].
  * (50) T. Han, W. Kilian, N. Kreher, Y. Ma, J. Reuter, T. Striegl, and K. Xie, “Precision test of the muon-Higgs coupling at a high-energy muon collider,” JHEP 12 (2021) 162, arXiv:2108.05362 [hep-ph].
  * (51) G.-S. Lv, X.-M. Cui, Y.-Q. Li, and Y.-B. Liu, “Pair production of the vectorlike top partner at future muon collider,” Nucl. Phys. B 985 (2022) 116016.
  * (52) J. Liu, Z.-L. Han, Y. Jin, and H. Li, “Unraveling the Scotogenic model at muon collider,” JHEP 12 (2022) 057, arXiv:2207.07382 [hep-ph].
  * (53) A. Azatov, F. Garosi, A. Greljo, D. Marzocca, J. Salko, and S. Trifinopoulos, “New physics in b → s$\mu$$\mu$: FCC-hh or a muon collider?,” JHEP 10 (2022) 149, arXiv:2205.13552 [hep-ph].
  * (54) J.-C. Yang, X.-Y. Han, Z.-B. Qin, T. Li, and Y.-C. Guo, “Measuring the anomalous quartic gauge couplings in the $W^{+}W^{-}\to W^{+}W^{-}$ process at muon collider using artificial neural networks,” JHEP 09 (2022) 074, arXiv:2204.10034 [hep-ph].
  * (55) Y. Bao, J. Fan, and L. Li, “Electroweak ALP searches at a muon collider,” JHEP 08 (2022) 276, arXiv:2203.04328 [hep-ph].
  * (56) S. Chen, A. Glioti, R. Rattazzi, L. Ricci, and A. Wulzer, “Learning from radiation at a very high energy lepton collider,” JHEP 05 (2022) 180, arXiv:2202.10509 [hep-ph].
  * (57) S. Baum, P. Sandick, and P. Stengel, “Hunting for scalar lepton partners at future electron colliders,” Phys. Rev. D 102 (2020) no. 1, 015026, arXiv:2004.02834 [hep-ph].
  * (58) B. Swiezewska and M. Krawczyk, “Diphoton rate in the inert doublet model with a 125 GeV Higgs boson,” Phys. Rev. D 88 (2013) no. 3, 035019, arXiv:1212.4100 [hep-ph].
  * (59) A. Arhrib, R. Benbrik, and N. Gaur, “$H\to\gamma\gamma$ in Inert Higgs Doublet Model,” Phys. Rev. D85 (2012) 095021, arXiv:1201.2644 [hep-ph].
  * (60) A. Jueid, J. Kim, S. Lee, S. Y. Shim, and J. Song, “Phenomenology of the Inert Doublet Model with a global U(1) symmetry,” Phys. Rev. D 102 (2020) no. 7, 075011, arXiv:2006.10263 [hep-ph].
  * (61) ATLAS, CMS Collaboration, G. Aad et al., “Measurements of the Higgs boson production and decay rates and constraints on its couplings from a combined ATLAS and CMS analysis of the LHC pp collision data at $\sqrt{s}=7$ and 8 TeV,” JHEP 08 (2016) 045, arXiv:1606.02266 [hep-ex].
  * (62) G. C. Branco, P. M. Ferreira, L. Lavoura, M. N. Rebelo, M. Sher, and J. P. Silva, “Theory and phenomenology of two-Higgs-doublet models,” Phys. Rept. 516 (2012) 1–102, arXiv:1106.0034 [hep-ph].
  * (63) S. Kanemura, T. Kubota, and E. Takasugi, “Lee-Quigg-Thacker bounds for Higgs boson masses in a two doublet model,” Phys. Lett. B313 (1993) 155–160, arXiv:hep-ph/9303263 [hep-ph].
  * (64) A. G. Akeroyd, A. Arhrib, and E.-M. Naimi, “Note on tree level unitarity in the general two Higgs doublet model,” Phys. Lett. B490 (2000) 119–124, arXiv:hep-ph/0006035 [hep-ph].
  * (65) I. F. Ginzburg, K. A. Kanishev, M. Krawczyk, and D. Sokolowska, “Evolution of Universe to the present inert phase,” Phys. Rev. D82 (2010) 123533, arXiv:1009.4593 [hep-ph].
  * (66) OPAL Collaboration, G. Abbiendi et al., “Search for anomalous production of dilepton events with missing transverse momentum in e+ e- collisions at s**(1/2) = 183-Gev to 209-GeV,” Eur. Phys. J. C 32 (2004) 453–473, arXiv:hep-ex/0309014.
  * (67) ATLAS Collaboration, G. Aad et al., “Search for electroweak production of charginos and sleptons decaying into final states with two leptons and missing transverse momentum in $\sqrt{s}=13$ TeV $pp$ collisions using the ATLAS detector,” Eur. Phys. J. C 80 (2020) no. 2, 123, arXiv:1908.08215 [hep-ex].
  * (68) ATLAS Collaboration, G. Aad et al., “Searches for electroweak production of supersymmetric particles with compressed mass spectra in $\sqrt{s}=$ 13 TeV $pp$ collisions with the ATLAS detector,” Phys. Rev. D 101 (2020) no. 5, 052005, arXiv:1911.12606 [hep-ex].
  * (69) L. M. Krauss, S. Nasri, and M. Trodden, “A Model for neutrino masses and dark matter,” Phys. Rev. D 67 (2003) 085002, arXiv:hep-ph/0210389.
  * (70) E. Ma, “Verifiable radiative seesaw mechanism of neutrino mass and dark matter,” Phys. Rev. D73 (2006) 077301, arXiv:hep-ph/0601225 [hep-ph].
  * (71) A. Ahriche, A. Jueid, and S. Nasri, “Radiative neutrino mass and Majorana dark matter within an inert Higgs doublet model,” Phys. Rev. D97 (2018) no. 9, 095012, arXiv:1710.03824 [hep-ph].
  * (72) S. Baumholzer, V. Brdar, and P. Schwaller, “The New $\nu$MSM ($\nu\nu$MSM): Radiative Neutrino Masses, keV-Scale Dark Matter and Viable Leptogenesis with sub-TeV New Physics,” JHEP 08 (2018) 067, arXiv:1806.06864 [hep-ph].
  * (73) D. Borah, P. S. B. Dev, and A. Kumar, “TeV scale leptogenesis, inflaton dark matter and neutrino mass in a scotogenic model,” Phys. Rev. D 99 (2019) no. 5, 055012, arXiv:1810.03645 [hep-ph].
  * (74) T. Kitabayashi, “Scotogenic dark matter and single-zero textures of the neutrino mass matrix,” Phys. Rev. D 98 (2018) no. 8, 083011, arXiv:1808.01060 [hep-ph].
  * (75) A. Ahriche, A. Arhrib, A. Jueid, S. Nasri, and A. de La Puente, “Mono-Higgs Signature in the Scotogenic Model with Majorana Dark Matter,” Phys. Rev. D 101 (2020) no. 3, 035038, arXiv:1811.00490 [hep-ph].
  * (76) S. Baumholzer, V. Brdar, P. Schwaller, and A. Segner, “Shining Light on the Scotogenic Model: Interplay of Colliders and Cosmology,” JHEP 09 (2020) 136, arXiv:1912.08215 [hep-ph].
  * (77) D. Borah, A. Dasgupta, K. Fujikura, S. K. Kang, and D. Mahanta, “Observable Gravitational Waves in Minimal Scotogenic Model,” JCAP 08 (2020) 046, arXiv:2003.02276 [hep-ph].
  * (78) M. Sarazin, J. Bernigaud, and B. Herrmann, “Dark matter and lepton flavour phenomenology in a singlet-doublet scotogenic model,” arXiv:2107.04613 [hep-ph].
  * (79) J. Hisano, T. Moroi, K. Tobe, and M. Yamaguchi, “Lepton flavor violation via right-handed neutrino Yukawa couplings in supersymmetric standard model,” Phys. Rev. D 53 (1996) 2442–2459, arXiv:hep-ph/9510309.
  * (80) E. Arganda and M. J. Herrero, “Testing supersymmetry with lepton flavor violating tau and mu decays,” Phys. Rev. D 73 (2006) 055003, arXiv:hep-ph/0510405.
  * (81) A. Ilakovac, A. Pilaftsis, and L. Popov, “Charged lepton flavor violation in supersymmetric low-scale seesaw models,” Phys. Rev. D 87 (2013) no. 5, 053014, arXiv:1212.5939 [hep-ph].
  * (82) T. Toma and A. Vicente, “Lepton Flavor Violation in the Scotogenic Model,” JHEP 01 (2014) 160, arXiv:1312.2840 [hep-ph].
  * (83) MEG Collaboration, J. Adam et al., “New constraint on the existence of the $\mu^{+}\to e^{+}\gamma$ decay,” Phys. Rev. Lett. 110 (2013) 201801, arXiv:1303.0754 [hep-ex].
  * (84) A. M. Baldini et al., “MEG Upgrade Proposal,” arXiv:1301.7225 [physics.ins-det].
  * (85) BaBar Collaboration, B. Aubert et al., “Searches for Lepton Flavor Violation in the Decays tau+- —> e+- gamma and tau+- —> mu+- gamma,” Phys. Rev. Lett. 104 (2010) 021802, arXiv:0908.2381 [hep-ex].
  * (86) Belle, Belle-II Collaboration, K. Hayasaka, “Results and prospects on lepton flavor violation at Belle/Belle II,” J. Phys. Conf. Ser. 408 (2013) 012069.
  * (87) SINDRUM Collaboration, U. Bellgardt et al., “Search for the Decay mu+ —$>$ e+ e+ e-,” Nucl. Phys. B 299 (1988) 1–6.
  * (88) A. Blondel et al., “Research Proposal for an Experiment to Search for the Decay $\mu\to eee$,” arXiv:1301.6113 [physics.ins-det].
  * (89) K. Hayasaka et al., “Search for Lepton Flavor Violating Tau Decays into Three Leptons with 719 Million Produced Tau+Tau- Pairs,” Phys. Lett. B 687 (2010) 139–143, arXiv:1001.3221 [hep-ex].
  * (90) CMS Collaboration, A. M. Sirunyan et al., “Search for lepton-flavor violating decays of the Higgs boson in the $\mu\tau$ and e$\tau$ final states in proton-proton collisions at $\sqrt{s}$ = 13 TeV,” Phys. Rev. D 104 (2021) no. 3, 032013, arXiv:2105.03007 [hep-ex].
  * (91) CMS Collaboration, V. Khachatryan et al., “Search for lepton flavour violating decays of the Higgs boson to $e\tau$ and $e\mu$ in proton–proton collisions at $\sqrt{s}=$ 8 TeV,” Phys. Lett. B 763 (2016) 472–500, arXiv:1607.03561 [hep-ex].
  * (92) Particle Data Group Collaboration, M. Tanabashi et al., “Review of Particle Physics,” Phys. Rev. D 98 (2018) no. 3, 030001.
  * (93) J. Herrero-Garcia, N. Rius, and A. Santamaria, “Higgs lepton flavour violation: UV completions and connection to neutrino masses,” JHEP 11 (2016) 084, arXiv:1605.06091 [hep-ph].
  * (94) R. J. Scherrer and M. S. Turner, “On the Relic, Cosmic Abundance of Stable Weakly Interacting Massive Particles,” Phys. Rev. D 33 (1986) 1585. [Erratum: Phys.Rev.D 34, 3263 (1986)].
  * (95) K. Griest and D. Seckel, “Three exceptions in the calculation of relic abundances,” Phys. Rev. D43 (1991) 3191–3203.
  * (96) F. Ambrogi, C. Arina, M. Backovic, J. Heisig, F. Maltoni, L. Mantani, O. Mattelaer, and G. Mohlabeng, “MadDM v.3.0: a Comprehensive Tool for Dark Matter Studies,” Phys. Dark Univ. 24 (2019) 100249, arXiv:1804.00044 [hep-ph].
  * (97) PICO Collaboration, C. Amole et al., “Dark Matter Search Results from the PICO-60 C3F8 Bubble Chamber,” Phys. Rev. Lett. 118 (2017) no. 25, 251301, arXiv:1702.07666 [astro-ph.CO].
  * (98) M. Cirelli, E. Del Nobile, and P. Panci, “Tools for model-independent bounds in direct dark matter searches,” JCAP 10 (2013) 019, arXiv:1307.5955 [hep-ph].
  * (99) R. J. Hill and M. P. Solon, “Standard Model anatomy of WIMP dark matter direct detection II: QCD analysis and hadronic matrix elements,” Phys. Rev. D 91 (2015) 043505, arXiv:1409.8290 [hep-ph].
  * (100) F. Bishara, J. Brod, B. Grinstein, and J. Zupan, “From quarks to nucleons in dark matter direct detection,” JHEP 11 (2017) 059, arXiv:1707.06998 [hep-ph].
  * (101) J. Ellis, N. Nagata, and K. A. Olive, “Uncertainties in WIMP Dark Matter Scattering Revisited,” Eur. Phys. J. C 78 (2018) no. 7, 569, arXiv:1805.09795 [hep-ph].
  * (102) XENON Collaboration, E. Aprile et al., “Projected WIMP sensitivity of the XENONnT dark matter experiment,” JCAP 11 (2020) 031, arXiv:2007.08796 [physics.ins-det].
  * (103) LZ Collaboration, D. S. Akerib et al., “Projected WIMP sensitivity of the LUX-ZEPLIN dark matter experiment,” Phys. Rev. D 101 (2020) no. 5, 052002, arXiv:1802.06039 [astro-ph.IM].
  * (104) DarkSide-20k Collaboration, C. E. Aalseth et al., “DarkSide-20k: A 20 tonne two-phase LAr TPC for direct dark matter detection at LNGS,” Eur. Phys. J. Plus 133 (2018) 131, arXiv:1707.08145 [physics.ins-det].
  * (105) J. Billard, L. Strigari, and E. Figueroa-Feliciano, “Implication of neutrino backgrounds on the reach of next generation dark matter direct detection experiments,” Phys. Rev. D 89 (2014) no. 2, 023524, arXiv:1307.5458 [hep-ph].
  * (106) M. Backović, A. Martini, O. Mattelaer, K. Kong, and G. Mohlabeng, “Direct Detection of Dark Matter with MadDM v.2.0,” Phys. Dark Univ. 9-10 (2015) 37–50, arXiv:1505.04190 [hep-ph].
  * (107) H. H. Patel, “Package-X: A Mathematica package for the analytic calculation of one-loop integrals,” Comput. Phys. Commun. 197 (2015) 276–290, arXiv:1503.01469 [hep-ph].
  * (108) G. Passarino and M. Veltman, “One Loop Corrections for e+ e- Annihilation Into mu+ mu- in the Weinberg Model,” Nucl. Phys. B 160 (1979) 151–207.
  * (109) T. Hahn and M. Perez-Victoria, “Automatized one loop calculations in four-dimensions and D-dimensions,” Comput. Phys. Commun. 118 (1999) 153–165, arXiv:hep-ph/9807565 [hep-ph].
  * (110) T. Hahn, “Generating Feynman diagrams and amplitudes with FeynArts 3,” Comput. Phys. Commun. 140 (2001) 418–431, arXiv:hep-ph/0012260 [hep-ph].
  * (111) ATLAS Collaboration, G. Aad et al., “Constraints on new phenomena via Higgs boson couplings and invisible decays with the ATLAS detector,” JHEP 11 (2015) 206, arXiv:1509.00672 [hep-ex].
  * (112) CMS Collaboration, V. Khachatryan et al., “Searches for invisible decays of the Higgs boson in pp collisions at $\sqrt{s}$ = 7, 8, and 13 TeV,” JHEP 02 (2017) 135, arXiv:1610.09218 [hep-ex].
  * (113) 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, arXiv:1809.05937 [hep-ex].
  * (114) G. Belanger, B. Dumont, U. Ellwanger, J. F. Gunion, and S. Kraml, “Global fit to Higgs signal strengths and couplings and implications for extended Higgs sectors,” Phys. Rev. D 88 (2013) 075008, arXiv:1306.2941 [hep-ph].
  * (115) J. Ellis and T. You, “Updated Global Analysis of Higgs Couplings,” JHEP 06 (2013) 103, arXiv:1303.3879 [hep-ph].
  * (116) CMS Collaboration Collaboration, “Projected performance of Higgs analyses at the HL-LHC for ECFA 2016,” Tech. Rep. CMS-PAS-FTR-16-002, CERN, Geneva, 2017.  https://cds.cern.ch/record/2266165.
  * (117) O. Cerri, M. de Gruttola, M. Pierini, A. Podo, and G. Rolandi, “Study the effect of beam energy spread and detector resolution on the search for Higgs boson decays to invisible particles at a future e+ e- circular collider,” Eur. Phys. J. C 77 (2017) no. 2, 116, arXiv:1605.00100 [hep-ex].
  * (118) D. M. Asner et al., “ILC Higgs White Paper,”  2013\.  arXiv:1310.0763 [hep-ph].  http://inspirehep.net/record/1256491/files/arXiv:1310.0763.pdf.
  * (119) CEPC Study Group Collaboration, M. Dong et al., “CEPC Conceptual Design Report: Volume 2 - Physics \& Detector,” arXiv:1811.10545 [hep-ex].
  * (120) M. Selvaggi, “Higgs measurements at the FCC-hh,” PoS ICHEP2018 (2019) 684.
  * (121) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations,” JHEP 07 (2014) 079, arXiv:1405.0301 [hep-ph].
  * (122) 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, arXiv:1108.2040 [hep-ph].
  * (123) A. Jueid and S. Nasri, “Work in progress,”.
  * (124) M. Belfkir, A. Jueid, and S. Nasri, “Work in progress,”.
  * (125) Particle Data Group Collaboration, R. L. Workman et al., “Review of Particle Physics,” PTEP 2022 (2022) 083C01.
  * (126) M. Belfkir, A. Jueid, and S. Nasri, “Work in progress,”.