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

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.


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.[1][2][3][4] 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 Ω DM h 2 = 0.1198 ± 0.0015 [5].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 miss T ).Examples of these processes are mono-jet [6], mono-Z [7,8] or mono-Higgs [9] 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 [10][11][12][13][14] which were interpreted in various particle physics realizations.Furthermore, these constraints were even more stringent when the void bounds from direct detection experiments [15,16] are included [17,18].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 [19][20][21][22][23].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. [24] and was widely studied in the literature [25][26][27][28][29][30][31].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 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 [29].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 [32][33][34] and therefore competing with the future circular colliders (FCC-hh).On the other hand, these machines can achieve very high energies thanks to arXiv:2301.12524v1 [hep-ph] 29 Jan 2023 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) [35,36].Phenomenology of both the SM and beyond at muon colliders has been extensively studied in the literature (see e.g.[37][38][39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56] 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 SM → γγ and theoretical constraints.In section III we discuss the constraints from charged lepton flavour violation in α → β γ, α → 3 β and H SM → α ¯ β .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.

A. The model
We consider a minimal extension of the SM by two gauge singlet fields: a charged scalar (S) and a righthanded 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 } → {−S, −N R } and { , q, ν, Φ, V µ } → { , q, ν, Φ, V µ }.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. 1 The full Lagrangian is given by where Φ refers to the SM Higgs doublet, L S is the interaction Lagrangian for the singlet scalar (including the kinetic term), and V (Φ, S) is the scalar potential.The interaction Lagrangian for the S field is given by with D µ S = (∂ µ − ig 2 Y S B µ /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 µ and Z µ which are given, after field rotations, by where e = √ 4πα EM is the electric charge, θ W is the Weinberg mixing angle, and The most general CP-conserving, renormalizable and gauge invariant scalar potential is given by process of electroweak symmetry breaking leads to three physical scalars: H SM identified with the recently discovered 125 GeV SM Higgs boson and a pair of charged scalars denoted by H ± .Their masses are given at the lowest order in perturbation theory by with υ being the vacuum expectation value (VEV) of the SM Higgs doublet.This model involves seven additional free parameters which we parametrise as follows For convenience we define the combination of the couplings which is a very good parametrisation in case the charged leptons are assumed to be massless.

B. 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 ± where more details can be found in [29].The width of the SM Higgs boson is only affected by the rate of its decay to γγ.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 λ 3 [58][59][60] 3 .In the present work, we have used the most recent ATLAS-CMS combined measurement of |κ γ | [61] We assume the theoretical prediction to be in agreement with the experimental measurement at the 2σ-level.We found that the enhancement of |κ γ | always occur for λ 3 < 0 which excludes charged scalars with masses up to ∼ 380 GeV [60].For λ 3 > 0, we get three possible regimes: (i) large and negative contribution that implies an enhancement of κ γ , (ii) positive but small contribution which makes κ γ consistent with the experimental measurement and (iii) exact or almost exact cancellation between the H ± and the W -boson contributions which make κ γ very small.Therefore, for λ 3 > 0, charged singlet masses up to 380 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 λ 4 = λ 5 = 0.In this study, we impose constraints from vacuum stability conditions (or boundness-from-below) [62], perturbativity, perturbative unitarity [63,64] and false vacuum [65].The false vacuum condition plays a very important role in constraining the parameters λ 2 , λ 3 and M H ± .We get We found that: (i) λ 3 cannot be larger 5 for all charged scalar masses and (ii) there is a parabola in the plane defined by λ 3 and M H ± which simply tells us that the smaller is the minimum allowed value of M H ± the smaller is the maximum allowed value of λ 3 .These conclusions are mildly dependent on the choice of λ 2 and, therefore, we choose λ 2 = 2 in the remainder of this manuscript without loss of generality.
The model can be constrained by using the null results of LEP and LHC searches for supersymmetric particles [66][67][68].The OPAL collaboration of the LEP experiment has searched for charginos decaying into a charged lepton and the lightest supersymmetric neutralino using 680 pb −1 of integrated luminosity [66].Assuming that the branching ratio of H ± → µ ± N R is 100%, the production of charged singlet pairs occurs through gauge interactions (s-channel diagrams with the exchange of γ * /Z 0 ).This search constrain the mass of the charged singlet to be not heavier than 100 GeV for any value of Y µ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 ∆ = m ˜ − m χ 0 ≥ 80 GeV [67] and compressed spectra for a mass splitting as low as 0.55 GeV [68].The two searches utilized a total luminosity of 139 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 ± up to 150 GeV (see the blue and green contours of figure 1).

C. 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 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 socalled Krauss-Nasri-Trodden (KNT) three-loop radiative neutrino mass model [69].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 [70].The phenomenology of the scotogenic model has been widely studied in the literature [71][72][73][74][75][76][77][78].The relevant interaction becomes where Φ IDM = (S, (h 2 + ia 2 )/ √ 2) T , and α, β are generation indices.Identifying (7) with the first term in L S we have Y eN = h 11 , Y µN = h 21 , Y τ 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 10 F and 5F representations, the charged singlet belongs to the 10 H representation, and the right handed neutrino belongs to the singlet representation 1 α , which in this case we can write In addition to the minimal SU (5), we can obtain the first term of equation ( 2) from a flipped-SU (5) ⊗ 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 10 α representation.In this case, we have where we integrated out a heavy intermediate state with a scale Λ Λ 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 number 4 .However, the charged singlet scalar can give rise to processes violating flavor lepton numbers L α ; α = e, µ, τ at the one-loop order.These processes called charged lepton flavor violating (CLFV) processes are categorised into three categories: (i) α = β γ, (ii) α → β β ¯ β and (iii) e-µ 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 α N come from the branching ratio of µ → eγ decay.The analysis of the CLFV decays in this work are heavily based on the results of refs.[79][80][81][82].A summary of the current and future bounds on the CLFV decays is shown in Table I.
Using the most recent experimental bounds on BR( α → β γ) from the Meg [83] and BaBar [85] experiments, we can use equation (10) to derive the following bounds on the products of the couplings: Since the one-loop function varies roughly between 1/12 and 1/6, the upper bound on the coupling Y αN Y β 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 ± and become very weak for heavy H ± .This can be clearly seen in figure 2 where the maximum allowed values of , 1000, and 5000 GeV.As expected, the bounds on |Y eN Y µN | are the strongest ones while the bounds on |Y τ N Y µN | and |Y eN Y τ N are similar.We must stress out that the CLFV decays do not constraints the Y α 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 experimen-tal bounds on CLFV decays:

B.
α It is noteworthy to discuss the constraints from the CLFV decays α → β β ¯ β .These processes receive four contributions at the one-loop order: penguin diagrams with the exchange of γ, Z and H 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 [82] where B ≡ BR( α → β ν α νβ ).The contribution of the γpenguins consist of the magnetic or dipole (A M ) and the non-dipole (A ND ) contributions.The dipole contribution is the same as of BR( β → β γ) but enhanced by a factor of 16 × (log(m α /m β ) − 22)/3 which varies between 7 and 36 for τ → 3µ and τ → 3e respectively.The non-dipole contribution is given by being the one-loop function for the non-dipole γ-penguin.This function has the following limits: , and τ → 3e respectively.The Zpenguin contribution is given by where g R , g L 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 θ W is the sine of the Weinberg mixing angle, and Z ND is the momentumindependent Z-boson form factor which is given by We can see that the Z-penguin contribution involves an extra suppression by a factor of m α m β as compared to the dipole γ-contribution.Finally, the box contribution is given by where D 1,2 (x) are the one-loop box functions given by D 3 .The contribution of the box diagrams, contrarily to penguins, has an extra factor of Y 2 β N .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 α → 3 β decays (see Table II).

C. HSM → α ¯ β
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 [90,91].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 α β is given by [93] BR( with ) is the Higgs-lepton Yukawa coupling of the heavier lepton (chosen here to be α ).In this formula, the contribution of the lighter lepton is neglected.We expect the bounds from H SM → α ¯ β searches to be very weak.This can be clearly seen in from table II 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 theN R particles in section IV A and then move to a detailed analysis of the spin-independent DMnucleus scattering cross section in section IV B. 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 D.

A. 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 ∆ ≡ M H ± −M N R < 0.1×M N R while the annihilation contributes for the whole parameter space.For the annihilation, there are two major contributions: β from the exchange of the charged scalar singlet in t-and u-channels, and (ii) N R N R → X∈SM XX which arises from the exchange of the SM Higgs boson via schannel 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 [1,94,95] with H = ȧ/a, n N R is the number density of the N R particle and 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 ∝ a −3 .This equation can be solved to give approximately where σ(x f )v is the thermally-averaged annihilation cross section for the N R particle where K 1 (x) and K 2 (x) are the modified Bessel functions of the second kind and σ N R N R → α β (ŝ) is the annihilation cross section into charged lepton which is given by where κi To simplify the discussion about the relic density, we consider the annihilation cross section in the limit ŝ → 4M 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 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 ± → SM and H ± H ∓ → SM.Below, we list the individual contributions and the overall dependence of the corresponding cross section with A i , B i , C i and 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 λ 3 and Y N and may even drive the relic density to very small values (∼ orders of magnitudes smaller than the observed abundance).In general, the co-annihilation is dominated by contributions of the following two processes In the presence of co-annihilations, the Boltzmann equations become where N is the mean number of N R particles, n H ± is the number density of H ± and Γ H ± is its total width.Note that here we have replaced the thermally-averaged annihilation cross section in equation ( 16) by the effective cross section 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 [96].In figure 3, we show the values of the coupling Y 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 ± and N R is large, we need large values of the Y N .However, even for Y N near the perturbativity bound the mass splitting can not be arbitrary large: ∆ max ≈ 600 (2000) GeV for M N R = 10 (100) GeV.The relic density becomes almost independent of Y 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 N .

B. 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 by5 with S p,n being the scalar current nucleon (p/n) form factors and 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 [98][99][100][101].
The generic formula of f p,n is given by Values of the coupling Y 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 Ωh 2 ≈ 0.12 are shown for Y N = 1, 2, 3, 4, 5, 7.5 and 4π.The blue shaded area corresponds to the region where the perturbativity is broken while the shaded gray region correspond to the kinematically forbidden region MN R > M H ± in which NR is not stable and therefore not a suitable dark matter candidate.FIG. 4. The spin-independent cross section as function of the dark matter mass MN R while the colored scatter points correspond to the charged singlet scalar mass M H ± .In the same plot, we show the current bounds from Xenon 1T [15] in dashed sienna, and the future expectations from Xenon nT [102], LUX LZ [103] and DarkSide G2 [104].The shaded orange area marked by 'neutrino floor' corresponds to the backgrounds from the coherent scattering with solar neutrinos, atmospheric neutrinos and supernova neutrinos [105].The spin-independent cross section was scaled by a factor of ξ Planck = ΩN R h 2 /Ω Planck h 2 with Ω Planck h 2 ≈ 0.12.All the calculations were performed for Y N = 2 and λ3 = 4.
The parton-level scattering amplitude is where A q is connected to the non-hadronic part of the amplitude.The term ψq (p out )ψ q (p in ) should be incorpo- where N = p, n.The model-dependent non-hadronic form factor is given by here ỹ(Q 2 ≈ 0) is the effective H SM N R N R coupling computed in the low energy limit.With the help of the Package X [107], we can obtain it from equation ( 30) where is monotonous and increasing function of x in the interval [0, 1] and has the following limits lim x→0 H(x) = 0 and lim x→1 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 ± which simply means that the direct detection spin-independent cross section is always below the neutrino floor for heavy H ± .From equation ( 28) one also expect that the spin-independent cross section is always proportional to |Y N | 4 .Therefore, large Y N regions with large σ SI would also correspond to small relic density (which is proportional to 1/|Y N | 4 ) 6 and for these scenarios σ SI needs to be scaled by a factor ξ Planck ≡ Ω N R h 2 /Ω Planck h 2 .This means that the spin-independent cross section would always be consistent with the current Xenon 1T bounds [15] for most regions of the parameter space as we can see clearly in figure 4.

C. 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 6 This is consistent with our previous finding in [29] where a strong anti-correlation between σ SI and Ω N R h 2 was observed.
with ỹ is the one-loop induced effective H SM -N R -N R coupling which is given by with ), i = 0, 2 being the Passarino-Veltman three-point functions [108].The computation of the Feynman amplitudes has been performed using FeynArts, FormCalc, and LoopTools [109,110].We have used a Python interface to LoopTools to evaluate numerically the one-loop integrals 7 .We define the Higgs invisible branching ratio as where Γ SM H = 4.07 MeV.Using equations ( 29), ( 30) and (31), we can obtain bounds on the coupling Y N .The bound is analytically defined by , ) and B bound is the upper bound on BR inv .
Searches for Higgs invisible decays have been carried by the ATLAS and the CMS collaborations [111][112][113].The strongest and up-to-date stringent bound on B inv was reported by the Cms collaboration using a combination 7 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 ( ) indicate that the parameter point yields a smaller σSI than the experimental bound (present or expected) while a cross mark (X) indicates that σSI is above the experimental bound.
of previous Higgs to invisible decay searches at 7, 8 and 13 TeV, where it has been found B inv < B bound = 0.19 at 95% CL [113] 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 [114,115].Finally, several studies have been devoted to the projected sensitivities of the future collider experiments to Higgs invisible decays from HL-LHC [116], FCC-ee [117], ILC [118], CEPC [119] and FCC-hh [120].In figure 5  The Higgs diphoton rate is assumed to be equal to the SM prediction at LO.

D. 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 ≤ 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 O(1).
The benchmark points used in the discussion of the general features of DM production at Muon colliders are shown in Table II.There are four of these benchmarks and each one has its own phenomenological implications.
a. BP1.This benchmark point is characterised by a relatively light DM (M N R = 50 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 µN Y τ N > Y eN .This choice leads to a charged singlet decaying predominantly into µN R with a branching fraction approaching 100%.On the other hand, the charged lepton flavour violating decays of charged leptons are such that BR(τ → eγ) 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.
b. BP2.For this point, we choose M H ± = 500 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 µN = 1.6, Y τ N = 5 × 10 −1 and Y eN = 5 × 10 −4 .This leads to the following branching ratios BR(H ± → µ ± N R ) 91%, BR(H ± → τ ± N R ) 9% and BR(H ± → e ± N R ) 0%.The charged singlet is narrow in this case as Γ H ± /M H ± 0.04.The CLFV decays of charged leptons exhibits similar features as in BP1 with the exception that BRs of τ → µN R and µ → eγ can be probed in the future experiments as they are slightly below the current bounds.The spinindependent DM-nucleon cross section can be tested by the DarkSide G2 experiment.c.BP3.For this point, we choose the following values of the particle masses: M H ± = 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 µN , Y τ N , Y eN } = {1, 0.5, 10 −3 } which leads to the following branching fractions: 0%.On the other hand, the branching ratios of CLFV decays of µ → eγ and τ → µγ 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 λ 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.
d. BP4.Here, we choose relatively heavy DM and charged singlet scalar; M N R = 1000 GeV and M H ± = 1500 GeV.The Yukawa-type couplings are chosen such that 50% while BR(H ± → e ± N R ).Similar features to BP1 and BP2 are observed for CLFV and DM phenomenology.

A. Total cross sections
In this section, we discuss the general features of DM production at muon colliders 9 .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-γ, mono-Z and mono-Higgs (a full analysis of these channels will be done in future work [123]).
• 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 [124]).
• 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.
In figure 6 we show the total cross sections for DM production in µµ collisions as function of the centerof-mass energy ( √ s µµ ) for the four benchmark points defined in table II.Starting with mono-X processes, it is clear that the mono-γ channel has the highest rate which varies from 1 pb for √ s µµ = 3 TeV to about 80 fb for √ s µµ = 30 TeV in BP1 10 .Mono-Z production has the second highest cross sections which varies between 9 The cross sections for both DM and charged scalar production at muon colliders are computed at Leading order using Madgraph_aMC@NLO [121] with a UFO model file [122] that can be found in the FeynRules model database https://feynrules.irmp.ucl.ac.be/wiki/ModelDatabaseMainPage. 10 Note that for mono-γ, we have applied some generator-level cuts by requiring that p γ T > 25 GeV and |η γ | < 2.5.
200 fb for √ s µµ = 3 TeV and about 2 fb for √ s µµ = 30 TeV. Finally, mono-Higgs production has the lowest rates among all the mono-X processes with cross section approaching 63 fb for √ s µµ = 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 µN since the total rates are proportional to Y 4 µN .An exception to this rule is in the mono-Higgs production cross section which decrease by factors of 6-200 since it scales as λ 2  3 Y 4 µN .
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 + γγ whose cross section is between 50 fb and 2 fb.This process is followed by N R N R γ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 − γ and N R N R γγγ at √ s µµ = 3 TeV.We note that the dependence on √ s µµ 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.

B. 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: and (ii) DM production in association with two SM particles where we consider five processes: where we follow ref.[40] 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 For the production of DM in association with two SM particles, we have TABLE III.The total cross sections times the branching ratio (σ × BR) and the expected number of signal events for the NRNR production in association with γ, Z(→ ), Z(→ q q), and HSM(→ b 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 refers to either an electron or a muon.
This process leads to the final state comprising of a highly energetic photon and a large missing transverse energy (E miss T ).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 -µ + µ − → Z(→ ν ν)γand VBF -V V → Z(→ ν ν)γ -with cross sections varying from 2.98 pb for √ s µµ = 3 TeV to 3.27 pb for √ s µµ = 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 √ s µµ = 30 TeV.It is worth noting from table III that the signal significance can easily 5 11 .For the other benchmark points, a more detailed selection is required to reach a signal significance of 5 if one can achieve an acceptance 11 The signal significance is defined as S/ √ B with S is the number of signal events and B is the number of background events.times efficiency (A × ) of about 15% for the signal in the signal region while the background is having A × of about O(10 −3 ).

b. N R N R Z(→
).This process leads to a very clean final state containing two same-flavour oppositesign (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 W W 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 O(10 3 -10 4 ).Given the differences in the topology of the signal and backgrounds, it is easy NRNRV (→ q q)V (→ q q) 1.05 × 10  5.29 × 10 −1 (5.29 × 10 3 ) 1.88 × 10 −1 (1.88 × 10 4 ) t t, V (→ q q)V (→ q q)ν ν 1.22 × 10   III but NRNR in association with γγ, γZ(→ ), Z(→ )Z(→ ), V (→ q q)V (→ q q and HSM(→ b b)HSM(→ b b).Here V refers to either W or Z.
to achieve a significance of 5σ by suitable event selection.
c. N R N R Z(→ q q) and N R N R H SM (→ b 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 q; q = u, d, s, c, b with BR(Z → q q) = 69.911%[125].For the SM Higgs boson, the main decay is into b b with BR(H SM → b b) = 57%.The dominant backgrounds to these signal processes come from q q production in association with two neutrinos, SM Higgs boson production, t t production with one top quark decaying leptonically and the other top decaying hadronically and W W production where one Wboson 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.
2. µ + µ − → NRNR + XY a. N R N R γγ.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-γ 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 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.
b. N R N R γZ(→ ).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  10 1 10 4 10 1 FIG. 7. Production cross section of H ± H ∓ + X as a function of the center-of-mass energy ( √ sµµ) 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 ± H ∓ plus one SM particle, plus two SM particles and in association with three SM particles.
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 O(10 2 -10 4 ).
c. N R N R Z(→ )Z(→ ).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 √ s µµ = 30 TeV where we expect about 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 V V * , V = W, Z in association with one or two gauge bosons.
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 ± H ± ) 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 γ/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: a. µµ → H ± H ∓ /H ± H ∓ H ± H ∓ .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 γ/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/ √ s µµ 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 [126].
b. µµ → H ± H ∓ + X.In this case, we have three production channels: H + H − γ, H + H − Z and H + H − H SM .There are three contributions to H + H − γ: s-channel contributions through γ/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 γ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 + MET, 2 + 2 jets + MET or 4 + 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.c. µµ → H ± H ∓ + XY .For this category we have seven production channels: H + H − γγ, H + H − γZ, H + H − ZZ, H + H − W + W − , H + H − H SM H SM , H + H − H SM Z and H + H − t 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.d. µµ → H ± H ∓ + XY Z.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 − γH SM H SM and H + H − γγH SM at √ s µµ = 3 TeV.
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 γγH + H − , γZH + H − , ZZH + H − , ZZ → H SM → H + H − and W + W − → γ/Z → H + H − vertices.We take examples of production of H + H − , H + H − γ and H + H − H 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 4 µN dependence of the cross section.FIG. 8.The production cross sections for H + H − , H + H − γ and H + H − HSM through VBF as a function of the centerof-mass energy ( √ sµµ) for BP1 (solid), BP2 (dashed), BP3 (dotted) and BP4 (dash-dotted).

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 sin-glets: 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 ∼ 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.

FIG. 2 .
FIG. 2. The maximum value of the products of the Y α N Y β N as a function of ξ for different values of the charged singlet mass M H ± .The results are shown for M H ± = 500 GeV (dashed), M H ± = 1000 GeV (dotted) and M H ± = 5000 GeV (dash-dotted).

FIG. 5 .
FIG.5.The present and future exclusions on values of Y N for M H ± = 500 GeV and λ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 MH SM = 125 GeV, and SM Higgs boson production rates.The Higgs diphoton rate is assumed to be equal to the SM prediction at LO.

FIG. 6 .
FIG.6.Production cross section of NRNR + X as a function of the center-of-mass energy ( √ sµµ) 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 NRNR plus one SM particle, plus two SM particles and in association with three SM particles.

TABLE IV .
Same as in table