. S capbtabboxtable[][]

A Duet of Freeze-in and Freeze-out:
Lepton-Flavored Dark Matter and Muon Colliders

Pouya Asadi, Aria Radick, Tien-Tien Yu
Institute for Fundamental Science and Department of Physics,
University of Oregon, Eugene, OR 97403, USA

We study a Lepton-Flavored Dark Matter model and its signatures at a future
Muon Collider. We focus on the less-explored regime of feeble dark matter
interactions, which suppresses the dangerous lepton-flavor violating
processes, gives rise to dark matter freeze-in production, and leads to long-
lived particle signatures at colliders. We find that the interplay of dark
matter freeze-in and its mediator freeze-out gives rise to an upper bound of
around TeV scales on the dark matter mass. The signatures of this model depend
on the lifetime of the mediator, and can range from generic prompt decays to
more exotic long-lived particle signals. In the prompt region, we calculate
the signal yield, study useful kinematics cuts, and report tolerable
systematics that would allow for a $5\sigma$ discovery. In the long-lived
region, we calculate the number of charged tracks and displaced lepton signals
of our model in different parts of the detector, and uncover kinematic
features that can be used for background rejection. We show that, unlike in
hadron colliders, multiple production channels contribute significantly which
leads to sharply distinct kinematics for electroweakly-charged long-lived
particle signals. Ultimately, the collider signatures of this lepton-flavored
dark matter model are common amongst models of electroweak-charged new
physics, rendering this model a useful and broadly applicable benchmark model
for future Muon Collider studies that can help inform work on detector design
and studies of systematics.

###### Contents

  1. 1 Introduction
  2. 2 Lepton-Flavored Dark Matter Model
    1. 2.1 Relic Abundance
    2. 2.2 Existing Constraints
  3. 3 Production at Muon Colliders
    1. 3.1 A Future Muon Collider
    2. 3.2 Flavored Dark Matter at a Future Muon Collider
  4. 4 A Search for the Prompt Region
    1. 4.1 Kinematics and Signal Regions
    2. 4.2 Target Systematics
  5. 5 A Search for the Long-Lived Region
    1. 5.1 Long-Lived Particle Kinematics
    2. 5.2 Background and Time-of-flight
    3. 5.3 Charged Tracks and Displaced Leptons Signal
  6. 6 Conclusions
  7. A More Details on Relic Abundance Calculation
  8. B More Details on the Prompt Search
  9. C More Details on the LLP Search

##  1 Introduction

The Standard Model (SM) of particle physics is an extremely successful model,
making precise and accurate predictions for a wide range of phenomena.
However, it fails to explain several fundamental phenomena. This includes
explanations for the origin of the electroweak scale, the existence of three
flavors of matter and their mass hierarchies, and the particle nature of Dark
Matter (DM), to name a few.

Many studies are undertaken in pursuit of an answer to these questions. On the
experimental front, collider experiments are unique in their ability in
probing a multitude of different models and unveiling various aspects of
potential beyond SM (BSM) particles above the electroweak scale.

To that end, various future collider experiments are proposed to push the
energy frontier. One such option that, thanks to increasing community interest
(see e.g. [1]), has risen as a compelling contender is a Muon Collider (MuC)
[2, 3, 4, 5, 6, 7]. There are considerable technical challenges that need to
be overcome before a viable high energy MuC can be constructed, see Refs. [8,
9, 1, 10, 11, 12] for recent reviews. Nonetheless, many theoretical studies
are carried out to underscore the reach of such a machine in probing different
BSM models and a template of well-motivated targets are being developed for
searches at such a facility. The clean environment of a lepton collider,
combined with substantial parton distribution function (PDF) of electroweak
gauge bosons in a muon at TeV scales, turns such a machine into a suitable
candidate for an in-depth study of the Higgs and precision electroweak
measurements [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26],
flavorful new physics [27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42], DM candidates [43, 44, 45, 46, 47], as well as other general SM and
BSM physics studies [48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61].
There are also proposals for earlier stage facilities, such as a beam dump
with accelerated muon beams, and interesting new physics that they can probe
[62, 63].

At the moment there exists no final detector design for such a future machine.
Theoretical works tabulating well-motivated models can inform development of
such a design. In particular, there are many well-motivated DM models that can
be searched for at a MuC, including various minimal Weakly-Interacting Massive
Particle (WIMP) DM models [43]. Many of these DM models, including WIMPs, are
strongly constrained by direct detection searches.

In this work, we focus on a less-constrained – yet similarly minimal – class
of DM models, namely flavored DM [64, 65, 66, 67, 68, 69, 70, 71, 72]. In such
setups, the DM, which is a singlet fermion of the SM, interacts with the SM
only via a new scalar mediator with the same charges as one of the SM fermion
fields. We focus on the setup where the new mediator has the same charges as
right-handed (RH) leptons of the SM, and study its relic abundance calculation
and signals at a MuC.

Previous works on flavored DM models have focused on DM production through
thermal freeze-out [64, 65, 66, 67, 68, 69, 70, 71, 72]. A sub-TeV DM mass
gives rise to the correct relic abundance in such setups. DM coupling to SM
leptons should be aligned with SM lepton Yukawa matrix in order to suppress
dangerous flavor-changing neutral current (FCNC) and lepton-flavor violation
(LFV) processes. These constraints can be abated by decreasing the size of the
DM-mediator Yukawa coupling, which opens up a new mechanism for DM production:
freeze-in [73]. The smallness of the Yukawa coupling also leads to a rich set
of collider signatures, ranging from prompt to long-lived particle (LLP).

In this work, we will study in detail the production mechanism for freeze-in
flavored DM, demonstrating how the interplay between mediator and DM abundance
leads to a bounded parameter space, as well as the collider signatures at a
future MuC. In particular, given the generic nature of the production rate of
our model, our signal yield can serve as a benchmark target that will inform
future experimental studies.

In Section 2 we introduce the flavored DM model that we study in this work. We
go through its relic abundance calculation in the freeze-in regime, as well as
a review of some existing LHC bounds on its parameter space. We will then
comment on its production channels at a future MuC in Section 3. We will then
study the kinematics of this model and prospects for its discovery in parts of
its parameter space where it decays promptly in the detector in Section 4. We
will also study its kinematics and signal yield in the LLP part of the
parameter space in Section 5, before concluding in Section 6. We also provide
details about relic abundance calculation in Appendix A. Further histograms
and explanations about the prompt (LLP) signal are included in Appendix B (C)

##  2 Lepton-Flavored Dark Matter Model

We start by introducing the model under study in this paper. We augment the SM
with a scalar mediator, $\phi$, that has the same charges under SM gauge
groups as RH charged leptons, and neutral fermion $\chi$,111With the same
field content we could also have $\phi LL$ and $H\bar{\chi}L$ terms that could
potentially destabilize DM. We postulate a $\mathbb{Z}_{2}$ parity acting on
$\phi$, $\chi$, and $\bar{\chi}$ that forbids these terms.

$\mathcal{L}\supset-
m_{\chi}\bar{\chi}_{\alpha}\chi^{\alpha}-m_{\phi}^{2}|\phi|^{2}-\lambda_{i,\alpha}\phi\bar{e}_{i}\bar{\chi}_{\alpha}.$
(2.1)

Here $i$ ($\alpha$) is SM fermion (DM) flavor index, $\bar{e}$ denotes SM
charged leptons of RH chirality, $\chi$ is the fermionic DM candidate, $\phi$
is the mediator that has the same charges under SM gauge group as a RH charged
lepton, and $m_{\chi}$ ($m_{\phi}$) denotes the shared DM (mediator) mass.
While having multiple DM flavors can give rise to interesting physics [64, 65,
66, 67, 68, 69, 70, 71], we will consider only one DM flavor as multiple
flavors will not affect our collider study.

It should be noted that such setups are studied under other names such as
Effective WIMP [74, 75], Fermion Portal DM [76, 77], and colored mediator
models [78, 79] as well. Flavored DM model can be considered a more general
framework compared to these models since it also includes multiple DM flavor
indices. Thus, we use the umbrella term flavored DM for our model, even though
having multiple DM flavors, for the most part, does not change our study. We
should also note that we focus on the case where the neutral particle, $\chi$,
is a fermion. We could also consider the scenario in which this DM candidate
is scalar and its mediator to SM is a charged heavier fermion.

In the existing studies of such setups (including a recent study of its signal
at a MuC [46]), DM is assumed to be a thermal relic, freezing out of the SM
bath at some point in the early universe. Nonetheless, DM abundance could also
be determined via the freeze-in mechanism [73] instead,222The original freeze-
in proposal is now known as IR freeze-in; see Ref. [80] for what is known as
UV freeze-in. where DM never reaches thermal equilibrium with SM. While there
exists some studies of flavored DM abundance via UV freeze-in [81] or LHC
signatures in this region of the parameter space [82], a proper study of the
interplay between the DM freeze-in and its mediator’s freeze-out, and their
implications for collider searches are still undelivered.

In the upcoming sections, we will study the relic abundance of DM in such a
setup more carefully, highlighting an interesting interplay of freeze-in and
freeze-out which gives rise to an upper bound on DM mass and interesting
signals at future colliders. It should be noted that we are implicitly
assuming the reheating temperature in the early universe is above the mediator
mass $m_{\phi}$, such that there initially exists an abundance of $\phi$ in
equilibrium with SM, which subsequently freezes out before decaying to DM
$\chi$. If the reheat temperature is below $\phi$ mass, DM is only produced
via UV freeze-in [80] and it will not have an upper bound on its mass.

The original flavored DM setups required an artificial introduction of a non-
trivial flavor ansatz to suppress dangerous FCNC and LFV. In the freeze-in
regime, smallness of the $\lambda$ Yukawa coupling guarantees an automatic
suppression of these processes. We will leave explorations of a natural
realization of such small Yukawa couplings in the UV to future works and focus
on their phenomenology at colliders in this work.

###  2.1 Relic Abundance

In this section we go through the relic abundance calculation of our setup. We
will show that the interplay between the freeze-in and the freeze-out
processes bounds our viable parameter space from all directions. Further
details about our calculation of the relic abundance can be found in Appendix
A.

We start by reasserting that $\phi$ reaches equilibrium in the early universe,
i.e. the reheat temperature is higher than the $\phi$ mass. As a result, we
should solve the coupled set of Boltzmann equations governing the abundances.
Neglecting Pauli blocking and Bose enhancement, and assuming symmetric
particle–anti-particle abundances, we can rewrite the evolution equations for
$\phi$ and $\chi$ abundances, as [83, 73]

$\displaystyle\frac{dY_{\phi}}{dx}=$ $\displaystyle-\sum_{\cal
F}\frac{m_{\phi}^{3}\langle\sigma v\rangle_{\phi\phi\to{\cal
F}}}{H(m_{\phi})x^{2}}\left[Y_{\phi}^{2}-Y_{\phi,\textrm{EQ}}^{2}\right]$
$\displaystyle-\sum_{\ell}\frac{x^{3}}{H(m_{\phi})}\frac{g_{\phi}\Gamma_{\phi\to\ell\chi}}{2\pi^{2}}K_{1}(x)\left[\frac{Y_{\phi}}{Y_{\phi,\textrm{EQ}}}-\frac{Y_{\chi}}{Y_{\chi,\textrm{EQ}}}\right],$
(2.2) $\displaystyle\frac{dY_{\chi}}{dx}=$
$\displaystyle\sum_{\ell}\frac{x^{3}}{H(m_{\phi})}\frac{g_{\phi}\Gamma_{\phi\to\ell\chi}}{2\pi^{2}}K_{1}(x)\left[\frac{Y_{\phi}}{Y_{\phi,\textrm{EQ}}}-\frac{Y_{\chi}}{Y_{\chi,\textrm{EQ}}}\right],$
(2.3)

where $H(m_{\phi})=1.66\sqrt{g_{*}^{\rho}}m_{\phi}^{2}/M_{\rm Pl}$ is the
Hubble constant at $m_{\phi}$ and $x=m_{\phi}/T$. $Y_{X}=n_{X}/T^{3}$ is the
yield of particle $X$ with $n_{X}$ the number density for particle $X$, and
$Y_{X,\mathrm{EQ}}$ is the equilibrium value of $Y_{X}$. The expression for
$Y_{\phi}$ consists of two parts that contribute to the depletion of $\phi$:
the first line corresponds to the annihilation of $\phi$, where $\langle\sigma
v\rangle_{\phi\phi\to{\cal F}}$ is the thermally averaged cross-section for
$\phi^{+}\phi^{-}\to{\cal F}$ with ${\cal
F}\in[\ell\bar{\ell},q\bar{q},\gamma\gamma,ZZ]$, while the second line
corresponds to the decay of $\phi$, where $\Gamma_{\phi\to\ell\chi}$ is the
decay width, $g_{\phi}=1$ is the number of spin degrees of freedom of $\phi$,
and $K_{1}$ is the first order modified Bessel function of the second kind.
The expression for $Y_{\chi}$ has one term which corresponds to the production
of $\chi$ from $\phi\to\ell\chi$. In the calculations that follow, we assume
the couplings $\lambda_{i}=\lambda$ are all equal. In Figure 1 we see the
relevant diagrams for our calculation. In this freeze-in calculation other
diagrams involving more $\lambda$ couplings or additional external particles
are suppressed and negligible [73].

Figure 1: Some diagrams relevant for Boltzmann equations governing freeze-out
of $\phi$ (top) and the freeze-in of $\chi$ (bottom). The freeze-out diagrams
are controlled by the electroweak gauge coupling (the black dots), while the
freeze-in process is proportional to the small Yukawa coupling $\lambda$ (red
dot). We assume comparable coupling to different lepton flavors in our
calculation. Other diagrams with higher number of $\lambda$ vertices or
particle multiplicity are suppressed.

We numerically solve these coupled differential equations. In Figure 2 we see
the results for benchmark masses and a fixed mediator lifetime. We have chosen
the Yukawa coupling $\lambda$ such that we get the correct final $\chi$ relic
abundance. From this figure, we see that there are two scenarios that can give
rise to the same final $\chi$ abundance. For the lighter $\phi$, freeze-in of
$\chi$ dominates its final abundance, gaining merely a small boost when $\phi$
decays, while for a heavy $\phi$ we see that freeze-out of $\phi$ and its
subsequent decay to $\chi$ dominates the final $\chi$ abundance.

Figure 2:  Abundance of $\chi$ (orange) and $\phi$ (blue) for the two
$m_{\phi}$ values that result in today DM abundance for $\chi$ with the given
$m_{\chi}$ and $\phi$ lifetime. Solid (dashed) lines correspond to the light
(heavy) $\phi$ scenarios. The figure shows that we can have scenarios where
$\chi$ abundance today is mostly dominated by freeze-in of $\chi$ itself
(solid) or freeze-out, and subsequent decay, of the mediator $\phi$ (dashed).
It also shows that with fixed lifetime and DM mass, there are two different
mediator masses (and $\lambda$ couplings) that give rise to the right relic
abundance.

We can then use the final $\chi$ yield to calculate its expected relic
abundance. Assuming that $\chi$ constitutes all of the DM we set
$\Omega_{\chi}h^{2}=0.12$, e.g. see Ref. [84], and find the value of the
coupling $\lambda=\lambda^{\prime}$ that recovers this for each pair of $\phi$
and $\chi$ masses. In our analysis we assume the same coupling between DM and
all SM RH leptons. In Appendix A we detail our semi-analytic approximation for
calculating $\lambda^{\prime}$.

In Figure 3 we show contours of $\lambda^{\prime}$ for every point on the mass
plane. For large enough $m_{\phi}$ values the $\phi$ freeze-out cross-section
will be small enough that we are left with too many $\phi$ particles after
their freeze-out such that, after their decay to $\chi$, the universe will be
overclosed. As we decrease $m_{\phi}$, and for a fixed $m_{\chi}$, at some
point the freeze-out $\phi$ abundance will be low enough that, after decaying
to $\chi$, we find the correct relic abundance today. The contribution from DM
freeze-in to the abundance today will be sub-dominant. For lower $m_{\phi}$
values, the DM abundance from $\phi$ freeze-out and decay is augmented by a
non-negligible freeze-in contribution to get the right relic abundance today.
Stability of DM also puts a lower bound on the mediator mass
$m_{\phi}\geqslant m_{\chi}$.

The upper bound on the mediator mass (from DM relic abundance) and the lower
bound on it (from DM stability), suggests a bounded viable parameter space.
This also implies an upper bound on the DM mass. Our relic abundance
calculation shows that this upper bound is

$m_{\chi}\lesssim 3.6\leavevmode\nobreak\ \mathrm{TeV}.$ (2.4)

Figure 3: The DM Yukawa coupling, $\lambda^{\prime}$, that gives rise to the
right relic abundance today. We assume the coupling to each fermion is the
same.

###  2.2 Existing Constraints

The on-going searches at LHC already constrain parts of our parameter space.
These constraints depend strongly on the lifetime of the mediator. In Figure 4
we show the lifetime of the mediator as given by

$\tau_{\phi}=(3\Gamma_{\phi\to\ell\chi})^{-1},$ (2.5)

where the factor of three is due to the three charged lepton channels. In this
calculation we assume the same branching ratio to each final state lepton.
This lifetime affects the signal of our model at future colliders as well as
at the LHC.

To find the bounds from LHC, we note that our mediator is identical to a
right-handed slepton. Many searches for sleptons have been carried out at the
LHC, e.g. see [85, 86], however, the majority of these searches assume that
the sleptons decay promptly. To map the prompt region onto our parameter
space, we take anything that has a proper lifetime of $\tau_{\phi}\lesssim
3\times 10^{-2}$ ns as prompt (see further details in section 3.2). With this,
the prompt searches do not reach high enough $\phi$ masses to be visible on
our bounds plot, regardless of combination of slepton flavors.

On the other hand, searches for long-lived sleptons can probe some parts of
the parameter space. In Figure 4 we show the bounds from the long-lived
sleptons in Ref. [87] on our model. In applying these bounds we assumed a
branching ratio of 1 to the strongest charged lepton channel (selectrons in
this case), to be optimistic about the reach of LHC in our parameter space. We
also show the bounds from CMS search for staus as heavy stable charged
particles in the detector [88]. Figure 4 clearly shows a substantial part of
the viable parameter space remains unavailable to LHC, and this motivates
searches for this model in future colliders, which we turn to next.
Additionally, limits from LEP are always taken into account by assuming that
$m_{\phi}>100$ GeV [89, 90] and small-scale cosmological structure observables
set a limit on the mass of freeze-in DM such that $m_{\chi}\gtrsim 10$ keV
[91].

Figure 4: The charged mediator $\phi$ proper lifetime, $\tau_{\phi}$, assuming
the Yukawa couplings required to get the right relic abundance (as given in
Figure 3). Depending on this quantity, the signal, and thus the collider
search strategy, varies significantly across the parameter space. The orange
(blue) region is already ruled out by the LHC search in Ref. [87] (Ref. [88]).
This figure shows a large part of the viable parameter space is not probed by
the LHC, motivating future collider searches. The dashed purple line denotes
the boundary between our prompt (Section 4) and LLP (Section 5) search
regions, see Section 3 for details.

##  3 Production at Muon Colliders

We now turn to signals of our flavored DM model at a MuC. In our analysis we
assume there exists only one flavor of DM. We also remain inclusive over the
final state charged leptons in our analysis. Our mediator production rate is
independent of the branching ratio to different charged leptons and can be
repeated straightforwardly for models where the mediator decays with different
branching ratios.

We will calculate the signal yield of our model at a future MuC and comment on
the reach of a future MuC in its parameter space. When possible, we will also
report the tolerable systematic uncertainty that will still allow a discovery
of our model with the aim that our results can inform the detector designs and
simulations, especially with regards to LLP signals. We start with a review of
the relevant aspects of a future MuC design and then discuss the production
cross-section of our model’s mediator at such a facility.

###  3.1 A Future Muon Collider

There are many on-going experimental studies about the feasibility and
features of a future MuC, including works on a detector design; see Refs. [8,
9, 1, 10, 11, 12] for recent reviews. A final and complete detector design is
still an on-going field of research. While this prevents a detailed study of
phenomenology of BSM models at a MuC, it underscores the importance of
developing a template of canonical models that, in turn, informs the
development of the detector.

While further studies are in order, recent progress in the cooling system [92]
has revived the hopes for achieving large luminosities at a MuC. The benchmark
target luminosity for a high energy MuC with center of mass energy of
$\sqrt{s}$ is given by [6]

$\mathcal{L}\approx 10\leavevmode\nobreak\
\mathrm{ab}^{-1}\times\left(\frac{\sqrt{s}}{10\leavevmode\nobreak\
\mathrm{TeV}}\right)^{2},$ (3.1)

which we use in our study as well. We will focus on bechmark center of mass
energies $\sqrt{s}=3$ TeV and $\sqrt{s}=10$ TeV, for which Eq. (3.1) suggests
benchmark $1$ ($10$) ab-1 integrated luminosities, respectively.

In Section 5 we outline an LLP search for our model, whose primary signature
is displaced leptons at different transverse distances from the beam. This
highlights the need for a rough estimate for the size of different components
of the detector, hence the need for assuming a benchmark detector design. We
use the design sketched in Ref. [10, 11], which in turn is informed by the
existing lepton collider designs, e.g. see the CLIC [93] and ILC designs [94].
The size of different segments of the detector are reported in Table 1 (table
taken from Ref. [10]). We will use the barrel region of these segments in our
LLP analysis of Section 5.

Subsystem | Region |  $L$ dimensions [cm] |  $|Z|$ dimensions [cm] |  $\mathbf{\eta}$ bound | Material
---|---|---|---|---|---
Vertex Detector | Barrel | $3.0-10.4$ | $65.0$ | $\lesssim 2.53$ | Si
| Endcap | $2.5-11.2$ | $8.0-28.2$ | $\lesssim 1.65$ | Si
Inner Tracker | Barrel | $12.7-55.4$ | $48.2-69.2$ | $\lesssim 1.05$ | Si
| Endcap | $40.5-55.5$ | $52.4-219.0$ | $\lesssim 1.07$ | Si
Outer Tracker | Barrel | $81.9-148.6$ | $124.9$ | $\lesssim 0.76$ | Si
| Endcap | $61.8-143.0$ | $131.0-219.0$ | $\lesssim 1.21$ | Si
ECAL | Barrel | $150.0-170.2$ | $221.0$ | $\lesssim 1.08$ | W + Si
| Endcap | $31.0-170.0$ | $230.7-250.9$ | $\lesssim 1.18$ | W + Si
HCAL | Barrel | $174.0-333.0$ | $221.0$ | $\lesssim 0.62$ | Fe + PS
| Endcap | $307.0-324.6$ | $235.4-412.9$ | $\lesssim 0.71$ | Fe + PS
Solenoid | Barrel | $348.3-429.0$ | $412.9$ | $\lesssim 0.85$ | Al
Muon Detector | Barrel | $446.1-645.0$ | $417.9$ | $\lesssim 0.61$ | Fe + RPC
| Endcap | $57.5-645.0$ | $417.9-563.8$ | $\lesssim 0.79$ | Fe + RPC

Table 1: Different components of the current proposed detector design for a
high energy MuC. Table is adapted from Ref. [10]. The $L$ ($Z$) dimension
refers to the segment size in the transverse (longitudinal) direction. The
$\eta$ bound column is added to the table of Ref. [10] and shows the highest
value of pseudorapidity $\eta$ for a track that completely goes through that
region. This information will enter our LLP search in Section 5.

The closest part of the detector to the beam is a vertex detector that spans
transverse distances of up to ten centimeters. The next part is the silicon-
based tracking system, filling up the space until transverse distance of
around a meter. Outside the tracking detector reside the silicon-tungsten ECAL
with thickness of tens of centimeters and a multi-layer HCAL, followed by the
magnetic field solenoids stretching roughly up to a radius of four meters away
from the beam. Finally, a muon spectroscopy system starts at around transverse
distances of four meters and surrounds the entire system.

In addition to these, to mitigate the effect of the beam-induced background
(BIB) at a MuC, two tungsten nozzles cladded with borated polyethylene are
included in each forward direction. These nozzles limit the access to forward
directions; we will assume pseudorapidity $|\eta|\gtrsim 2.4$ is not
accessible. See Refs. [94, 10, 11] for further details about the design.

We should remind the reader that Table 1 should not be construed as the final
MuC detector design. Our upcoming analysis can straightforwardly be repeated
once a concrete design is developed.

The design in Ref. [10] uses a large magnetic field as well. Nonetheless, the
effect of this magnetic field will be negligible on our study. In particular,
our proposed LLP search in Section 5 relies on transverse displacement of
displaced leptons, which is not particularly sensitive to the magnetic field
(for any feasible magnetic field size in the detector).333An adequately
precise measurement of some LLP signals, such as closest approach of the track
to the primary vertex, requires a proper inclusion of the magnetic field
effect - see Ref. [95] for a recent review. However, we will not study such
signals in this paper.

###  3.2 Flavored Dark Matter at a Future Muon Collider

In this section we go over basics of our model signals in a future MuC. Since
we are focusing on the freeze-in regime, the Yukawa coupling of the DM to SM
leptons and the mediator $\phi$ is very small, thus direct DM production is
strongly suppressed. The mediator $\phi$, on the other hand, is charged under
SM gauge groups and can be produced ubiquitously if kinematically allowed.
Subsequently, $\phi$ can only decay to the DM and a lepton, giving rise to a
charged lepton and missing energy in the detector. Since the production is
determined by SM gauge coupling, our model is a well-motivated target whose
signal yield can be used for a plethora of different BSM models, and can
inform future detector designs.

Diagrams relevant for $\phi$ production at a MuC are shown in Figure 5. The
kinematic distribution of final states will change depending on the production
channel. Similar topologies exist for other flavored DM setups (where $\phi$
has same charges as other SM fermions). They give rise to a pair of $\phi$
particles, each of which subsequently decays to a DM particle and a lepton.
The branching ratio to different charge leptons will depend on the Yukawa
couplings. Without any further structure on the flavor texture, we expect
these different branching ratios to be comparable. We remain agnostic about
the value of these branching ratios in this work.

Figure 5: Diagrams pair-producing $\phi$ at a future MuC. The initial state
can be either SM charged leptons (top), or electroweak gauge bosons (bottom).
We refer to these production channels as Drell-Yan (DY) or vector boson fusion
(VBF), respectively. The produced $\phi$ particles are on-shell and will
eventually decay to DM $\chi$ and a charged lepton. All diagrams involving an
off-shell $\phi$ or their single-production are suppressed by the small Yukawa
coupling to DM in our freeze-in setup.

The top diagram in Figure 5, which resembles a Drell-Yan (DY) production,
comes from initial state muons, whose PDF [96, 97, 98] is dominated by
longitudinal fractional momentum $x=1$, while the other three diagrams are
similar to vector boson fusion (VBF) production channels. In the VBF
production channels the original lepton typically goes undetected in the
forward direction down the beam. Nonetheless, it sometimes recoils enough to
fall within the detectable $\eta$ range and will give rise to an extra charged
lepton in the final state. We will not consider such events in our search.

We use MadGraph5 [99] for event generation and study of kinematic
distributions. Based on the existing results on various SM particles PDF in a
muon [96, 97, 98], in our simulations we assume the muons are all produced
with longitudinal fractional momentum $x=1$, while electroweak gauge bosons
are modeled with the effective vector approximation (EVA) and already included
in MadGraph5 [98]. The total cross-section for $\phi$ pair production, from
various initial states are shown in Figure 6. We can see from these figures
that for light (heavy) mediator masses the production is strongly dominated by
VBF (DY) production channels. The contribution from DY processes with other SM
fermion initial states are negligible compared to the aforementioned
production channels.

Figure 6: Cross-section for pair production of $\phi$ mediators at a future
MuC with center of mass energy of 3 TeV (blue) or 10 TeV (orange). We do this
calculation using MadGraph5 [99]. Results for different initial states are
shown with line styles, with DY as solid, $W$ VBF as dot-dashed, and
$\gamma/Z$ VBF as dashed. We find that the production is dominated by the DY
channel for the mediator masses above a few hundred GeV.

The model’s signal at a collider is substantially affected by the lifetime of
the $\phi$ mediator, see Figure 4. We divide the parameter space of the model
by $\phi$ lifetime, $\tau_{\phi}$. Specifically, we will consider two
different strategies.

  * •

For $\tau_{\phi}\lesssim 3\times 10^{-2}$ ns the mediator decays promptly to a
lepton and the missing energy in the vicinity of the interaction point.

  * •

For $\tau_{\phi}\gtrsim 3\times 10^{-2}\leavevmode\nobreak\ \mathrm{ns}$ the
produced mediator particles are LLPs, moving inside the detector for a
macroscopic distance and leaving a charged track. Depending on the lifetime,
they can either decay inside the detector (giving rise to displaced leptons),
or escape it (like heavy stable charged particles). In either case, they will
leave a detectable charged track in some parts of the detector.

It should be emphasized that our division of the parameter space into prompt
and LLP is dependent on detector geometry and exact process details; whether a
cut-and-count analysis for prompt signals is applicable to longer lifetimes,
or if an LLP search strategy to shorter lifetimes, should be studied on a
case-by-case basis and once a final detector design exists.

In the upcoming sections we detail a search strategy and calculate the signal
yield in each region. We should emphasize that our proposal is a rudimentary
analysis intended to lay out simple search strategies and to demonstrate the
power of a high enregy MuC in probing this well-motivated model. In the prompt
region we calculate the SM background and report the systematics uncertainties
that can be tolerated while still allowing for a discovery. In the LLP region
the irreducible background is generated via BIB and the detector response,
better understanding of which requires simulations beyond the current work.
Thus, in this region we simply report the signal yield and some of its unique
features, leaving further studies of the background, systematics, and the
reach of a high energy MuC to future works.

##  4 A Search for the Prompt Region

###  4.1 Kinematics and Signal Regions

As the name implies, in this part of the parameter space the signal of our
model will be a pair of charged leptons and missing energy coming from the
prompt decay of the mediator $\phi$. In this section we go over the kinematics
of events in this part of the parameter space and put forward a rudimentary
search proposal. We will show that our proposal can discover our model almost
over the entire kinematically-accessible region.

The cross-section for $\phi$ pair production, from various initial states are
shown in Figure 6. The main background for this search comes from electroweak
gauge bosons pair production and decay in Figure 7.

Figure 7: Dominant SM background for our search in the prompt region. The
$W^{+}W^{-}$ pair-production is potentially the most constraining background;
we can suppress this background by applying $M_{T2}$ cuts, see the text for
further details.

To reduce these backgrounds, we propose a cut on the event’s transverse
missing mass $M_{T2}$ [100], the lepton system invariant mass $m_{\ell\ell}$,
and the lepton system transverse momentum $p_{T,\ell\ell}$. The background
from a pair production of $Z$ bosons can be significantly reduced by
$m_{\ell\ell}$ or $p_{T,\ell\ell}$. This is due to the fact that charged
leptons originating from SM gauge bosons decays have a much lower invariant
mass $m_{\ell\ell}$ and transverse momentum $p_{T,\ell\ell}$ compared to
leptons coming from our mediator’s decay. The remaining background from pair
production of the $W$ boson can also be strongly suppressed by the use of
$M_{T2}$. To see that, we should keep in mind that in our pair production
topologies $M_{T2}$ is bounded from above by the mediator mass [100]. Thus,
$W$ pair production events have much lower $M_{T2}$ values than the signal
events from decays of $\phi$.

Figure 8: Histograms of invariant lepton mass (top), lepton transverse
momentum (middle), and $M_{T2}$ (bottom) for our model (colored) and for SM
(gray) at a MuC with center of mass energies of 3 TeV (left) and 10 TeV
(right). The bin size for each observable is 50 (200) GeV for $\sqrt{s}=3$
$(10)$ TeV. We use MadGraph5 with $10^{5}$ events to generate each histogram.
The solid (dashed) colored histograms refer to signal processes from the DY
(VBF) channel. In each figure we show the distribution for two $\phi$ masses
towards the two kinematically accessible ends of the mass ranges. We find a
sharp contrast between the SM’s and our model’s predictions, motivating use of
these cuts in our analysis.

To better justify using these cuts, in Figure 8 we show the distribution of
events in various kinematic observables, as well as the SM background
distributions. We find that the distributions of the SM background and our
model differ by orders of magnitude in certain ranges of $M_{T2}$,
$m_{\ell\ell}$, and $p_{T,\ell\ell}$, hence motivating use of these variables.
In particular, we find SM background distributions are concentrated at lower
values of these parameters, thus justifying the use of a lower bound cut on
their value in our analysis. We should also note the difference in the
distribution of events from the two signal channels (DY and VBF - see Figure
5).

We use MadGraph5 [99] to generate events for these histograms. For SM
background, we consider the processes
$\mu^{-}\mu^{+}\rightarrow\ell^{-}\ell^{\prime+}\bar{\nu}_{\ell}\nu_{\ell^{\prime}}$
with $\ell^{(^{\prime})}\in[e,\mu]$. Other background channels that come from
VBF diagrams (namely $VV\rightarrow\mu\mu\nu\nu$) are divergent at leading
order and need to be resummed. Instead of the resummation, we note that
putting cuts on $M_{T2}$ and $m_{\ell\ell}$ at above $W$ and $Z$ masses will
strongly suppress this background as well, hence we do not include it in our
study. A proper resummation and inclusion of such processes are left for
future.

We propose putting a cut on the three kinematic variables $M_{T2}$,
$m_{\ell\ell}$, and $p_{T,\ell\ell}$ and counting the number of events from
our model and from SM backgrounds. We propose different signal regions (with
different cuts on each kinematics observable) and allow for different signal
regions to have overlaps. For every point in the parameter space we choose the
signal region that has the most sensitivity to our signal.

Assuming a Gaussian distribution of SM background events in each signal
region, the statistical error bar on each signal region is simply square root
of the average number of SM events in that region, denoted by $\sqrt{B}$.
Hence, we can quantify the significance of the signal count $S$ in each signal
region by simply calculating $S/\sqrt{B}$. (For signal regions for which our
MadGraph5 simulations predict $\sqrt{B}\leqslant 2$, we use statistical error
of $\sqrt{B}=2$ instead to be conservative.) For any given mediator mass, all
regions with $S/\sqrt{B}\geqslant 5$ allow a $5\sigma$ discovery of our
model.444Since we use individual signal regions for a discovery, we do not
need to resort to more sophisticated likelihood calculations; repeating this
analysis with more complete statistical treatment is straightforward.

For $\sqrt{s}=10$ TeV, we find the cuts on $m_{\ell\ell}$, $p_{T,\ell\ell}$,
and $M_{T2}$ that optimize our sensitivity to the model with $m_{\phi}=1,4.5$
TeV:

$\displaystyle\mathrm{Optimized\leavevmode\nobreak\
Sensitivity\leavevmode\nobreak\ }m_{\phi}=1\leavevmode\nobreak\
\mathrm{TeV}\leavevmode\nobreak\ :\leavevmode\nobreak\
(m_{\ell\ell},p_{T,\ell\ell},M_{T2})=(3.50,0.83,0.20)\leavevmode\nobreak\
\mathrm{TeV}$ (4.1) $\displaystyle\mathrm{Optimized\leavevmode\nobreak\
Sensitivity\leavevmode\nobreak\ }m_{\phi}=4.5\leavevmode\nobreak\
\mathrm{TeV}\leavevmode\nobreak\ :\leavevmode\nobreak\
(m_{\ell\ell},p_{T,\ell\ell},M_{T2})=(0.99,0.20,3.50)\leavevmode\nobreak\
\mathrm{TeV}$ (4.2)

For $\sqrt{s}=3$ TeV, we find the cuts on $m_{\ell\ell}$, $p_{T,\ell\ell}$,
and $M_{T2}$ that optimize our sensitivity to the model with $m_{\phi}=1,1.3$
TeV

$\displaystyle\mathrm{Optimized\leavevmode\nobreak\
Sensitivity\leavevmode\nobreak\ }m_{\phi}=1\leavevmode\nobreak\
\mathrm{TeV}\leavevmode\nobreak\ :\leavevmode\nobreak\
(m_{\ell\ell},p_{T,\ell\ell},M_{T2})=(0.48,0.75,0.61)\leavevmode\nobreak\
\mathrm{TeV}$ (4.3) $\displaystyle\mathrm{Optimized\leavevmode\nobreak\
Sensitivity\leavevmode\nobreak\ }m_{\phi}=1.3\leavevmode\nobreak\
\mathrm{TeV}\leavevmode\nobreak\ :\leavevmode\nobreak\
(m_{\ell\ell},p_{T,\ell\ell},M_{T2})=(0.31,0.20,0.88)\leavevmode\nobreak\
\mathrm{TeV}$ (4.4)

Further details about how we arrive at these cuts are included in Appendix B.
We use the benchmark mass $m_{\phi}=1$ TeV since it is around the lowest
mediator mass inside the prompt region, see Figure 4; $m_{\phi}=4.5$ TeV
($m_{\phi}=1.3$ TeV) for $\sqrt{s}=10$ TeV ($\sqrt{s}=3$ TeV) is also around
the maximum mass ($\sqrt{s}/2$) that is kinematically accessible. We should
note that since the mediator coupling to DM is very small, the signal is
completely suppressed for the mass ranges where the mediator can not be
produced on-shell.

We use the cuts from Eqs. (4.1)–(4.4) and the average of the cuts on the two
extreme mass cases for each center of mass energy to define 27 overlapping
signal regions for each center of mass energies. For each point in the
parameter space we identify the signal region that has the best $S/\sqrt{B}$.
In Figure 9 we show signal regions, as defined in Appendix B, that are most
sensitive to each point in the parameter space and allow a $5\sigma$
discovery. The marked mass range for each signal region shows where that
signal region gives rise to a discovery, while neglecting systematics; the
small dots for each mass point indicate the best-fit signal region for that
mediator mass. It should be noted that this analysis is independent of the DM
mass $m_{\chi}$; for virtually every point in the prompt decay region the DM
can be treated as massless.

Figure 9: The range of mediator masses that each signal region could discover
with $5\sigma$ confidence at a future MuC with center of mass energies of 3
TeV (top) or 10 TeV (bottom), assuming zero systematic uncertainties. The
lines indicate the mass range in which each signal region yields a discovery,
while the dots denote the best-fit signal region for each mass point. We find
that for all masses above a TeV, and almost upto $\sqrt{s}/2$, a future MuC
can discover the signal of our model.

We find that the cuts are strong enough that essentially any part of the
parameter space where $\phi$ pairs can be produced on-shell can be probed in
future MuCs. We find that for most mediator masses there are multiple signal
regions that allow a $5\sigma$ discovery.

###  4.2 Target Systematics

We have thus far ignored various systematics uncertainties for each signal
region. There are many on-going experimental efforts for designing the MuC
detector and estimation of systematic uncertainties. Tolerable systematics on
well-motivated models, such as the ones studied in this work, can serve as a
motivated target in such studies.

In Figure 10 for each $\phi$ mass we calculate the maximum systematic
uncertainties (as a fraction of the statistical uncertainty in that region)
that still allows a discovery of our model for that point in the parameter
space. In doing so, we assume systematics and statistical uncertainties are
uncorrelated. (We also used $\sqrt{B}=2$ for signal regions that have
$\sqrt{B}\leqslant 2$ in our simulations.)

Figure 10: Tolerable systematic uncertainty for each signal region, as a
fraction of statistical uncertainty, that still allows a discovery of our
model for different mediator masses, see Figure. 19 for definitions of signal
regions. The left (right) column includes signal regions at a 3 TeV (10 TeV)
MuC. We only include signal regions with more than twenty signal events.
Signal regions not included in this figure either have fewer signal events or
do not give rise to a $5\sigma$ discovery of the signal, even in the absence
of systematics uncertainties.

This figure shows that even with systematics uncertainties that are comparable
to the statistical error in the most-sensitive signal regions, we can still
discover our model in the kinematically accessible parts of its parameter
space.

Finally, we should point out that in our analysis we used the average
background prediction in each signal region according to the results of
MadGraph5 simulations. Theoretical uncertainties on this prediction should
also be included in an actual analysis by smearing the SM prediction in each
signal region by a Gaussian distribution. We will leave a thorough treatment
of all systematics, statistical, and theoretical uncertainties for future
works.

It should also be noted that our search can be repeated for other types of
flavored DM models. In the case of quark flavored models, we expect to have a
jet, instead of the charged lepton, from the primary vertex and our search
should be updated accordingly. We expect cuts on similar set of kinematics
observables will allow us to probe the prompt region parameter space of all
such models as well.

##  5 A Search for the Long-Lived Region

For longer lifetimes, the pair of mediators can leave a charged track in the
detector and give rise to various LLP signals, see Refs. [101, 95] for recent
reviews of LLPs. In our setup, depending on the lifetime, the produced $\phi$
particles subsequently either decay to a charged lepton and a DM particle
inside the detector, or escape the detector without decaying. Further studies
are required for better understanding of the background for these searches,
see the discussion below. This prevents us from reporting the discovery reach
of a future MuC in our model’s parameter space at the moment. Instead, here we
focus on the signal yield of our model. In particular, we will calculate the
distribution of expected number of displaced leptons in different detector
components.555In our model the mediators are always produced in pairs. Thus,
each signal event includes two charged tracks and two displaced leptons that
could further enhance the search [102]. We, however, will not explicitly use
this property in our rudimentary analysis here. We also show that our model
exhibits a double-peak feature in some kinematic distributions, that
distinguishes it from SM background. Given the fact that our setup is a
theoretically-motivated target and its production rate and signals are generic
for any electroweakly-charged LLP, we hope this signal yield can serve as a
target for experimental studies and detector design efforts.

###  5.1 Long-Lived Particle Kinematics

Previous studies of LLPs at a high energy MuC focused on scenarios where the
long lifetime is an artifact of a small mass splitting between the LLP and its
neutral daughter [103, 104, 105]. In such scenarios the charged decay product,
e.g. a SM charged lepton, can carry a small momentum and escapes detection,
giving rise to a disappearing track signature. In our setup, on the other
hand, the mediator and DM have a large mass splitting and the long lifetime of
the mediator is an artifact of a small Yukawa coupling. As a result of this
large mass splitting, if the mediator decays inside the detector the charged
lepton will be easily detectable.

Figure 11: Distribution of events in $\phi$’s $\beta\gamma$ (top) and $\eta$
(bottom) for different mediator masses at a MuC with center of mass energy of
$\sqrt{s}=3$ TeV (left) or $\sqrt{s}=10$ TeV (right). The bin sizes are 0.2
and 0.5 for $\beta\gamma$ and $\eta$, respectively. The peak feature in the
$\beta\gamma$ is due to events from the DY production channel, see Figure 5.

We study the distribution of $\beta\gamma$, where $\beta=v/c$ and
$\gamma=(1-\beta^{2})^{-1/2}$, and pseudorapidity $\eta$. In Figure 11 we show
1D histograms of $\beta\gamma$ and $\eta$ for a few different points on our
parameter space and for center of mass energies of $\sqrt{s}=3$ TeV and
$\sqrt{s}=10$ TeV. In making these histograms we run MadGraph5 with $10^{5}$
events for each production channel (DY and VBF). We then weigh the number of
events from each production channel by its corresponding cross-section to get
the combined distribution.

In Figure 11 we see features that are distinct to the individual production
channels, DY and VBF, which in turn help inform the types of cuts we can use
in a search for our model. In the DY channel, the initial muons each carry
$E=\sqrt{s}/2$, which in turn gives rise to the maximum value of
$\gamma=E/m_{\phi}$ for a given $m_{\phi}$ and leads to the very sharp peak in
the $\beta\gamma$ distribution. An increase in $m_{\phi}$ results in a smaller
$\gamma$, so the peak of the $\beta\gamma$ distribution will shift to lower
values.666Currently the muon PDF in a muon beam is not implemented in
MadGraph5. Thus, in our simulations we assumed each muon enters the
interaction with $\sqrt{s}/2$ energy as an approximation. The PDF of a muon in
the beam has support for $x\neq 1$ as well [96, 97, 98]. As a result, the real
distribution of $\beta\gamma$ from this channel is slightly diffused to values
smaller than $\sqrt{s}/2$. Including the PDF will soften the peak in the
distribution of events from the DY channel. Nonetheless, given the sharp peak
in a muon PDF [96, 97, 98], we still expect a large peak in various kinematic
distributions from this channel if the real muon PDF is included in the
simulation.  We find that this channel gives rise to centralized $\phi$
particles (small $|\eta|$) as well. This motivates searches focusing on this
region. The VBF production channel gives rise to a wider spread in both
$\beta\gamma$ and $\eta$, owing to the wider spread of initial gauge bosons
PDF [96, 97, 98].

###  5.2 Background and Time-of-flight

The background for LLPs can be divided into a reducible SM background, and an
irreducible background from the detector response. While there are a handful
of SM particles that appear as LLPs at a collider (e.g. see Ref. [106]), using
various kinematic observables, such as time-of-flight or anomalous ionization
rate ($dE/dx$), allows us to distinguish them from heavy new physics LLPs.
Other techniques such as empty bunch crossing also allows us to better
understand and reject SM background; see Ref. [101] for further discussions of
SM background mitigation techniques.

A useful quantity in reducing the SM background is the time-of-
flight.777Another useful quantity in LLP searches is the anomalous ionization
rate $dE/dx$. In particular, this quantity can be used in reconstructing the
track mass and cutting away SM background. At LHC, measuring $dE/dx$ has a
smaller uncertainty when only the barrel region pixel layers are used. This
motivates LLP searches focusing on the central events at a MuC as well. A
proper calculation of this rate depends on characteristics of the detector
material and tracks interactions with the detector, which can only be
extracted from simulations. As a result, we acknowledge the possibility of
using this quantity and its importance in the LLP searches, leaving more in-
depth studies for future. In particular, the difference between the arrival
time of a $\phi$ particle and SM particles in the same direction, $\Delta
t=t_{\rm LLP}-t_{\rm SM}$, can be used to distinguish it from SM particles,
which move with $\beta=1$. In Figure 12 we show the distribution of this
variable, at the inner border of the muon system ($L=446.1$ cm transverse
distance from the beam in the design of Table 1). The time difference scales
linearly with the transverse distance $L$ at which we measure it, thus we can
use this figure to calculate the time difference at different transverse
distances in the detector.

Figure 12: Distribution of the difference between the arrival time of $\phi$
and a massless particle moving in the same direction at the onset of the muon
system ($L=446.1$ cm), as predicted in our MadGraph5 simulation for three
different mediator masses. In these histograms the log-scale x-axis is binned
in steps of 0.2 decades. The arrival time scales linearly with $L$. Events at
lower (higher) $\Delta t$ values correspond to the DY (VBF) production
channel. As the mediator mass increases, the LLPs move more slowly, shifting
the distributions to higher $\Delta t$ values.

Small $\Delta t$ events in Figure 12 are produced by the DY initial channel
(with larger $\beta$ and $\beta\gamma$ values), while events from the VBF
channel arrive later since they have lower values of $\beta$. This figure
shows that, at $L=446.1$ cm, timing resolution of around $0.1-1$ ns ($1-100$
ns) is enough to separate most events from the DY (VBF) channel from SM
background. When compared to the current timing resolution at LHC ($\sim
0.1-1$ ns for comparable $L$ values), we find that majority of events,
especially from the VBF channel or for large enough $m_{\phi}$ values, arrive
with enough delay that we can use the time-of-flight information to discern
our LLP signal from background.

In addition to the SM reducible background discussed above, interactions of
the beam, various tracks, or the BIB, with the material in the detector give
rise to an irreducible background. Currently, there are limited studies on
this topic [103]. Nonetheless, they have shown that such background events can
be mitigated by timing information, track quality, and tracks directionality
information. Further studies are in order to completely understand the
detector response and background for LLP searches at a high energy MuC and,
depending on the lifetime, different background mitigation strategies have to
be deployed.

To be able to discover our model, we need to estimate this irreducible
background as well. This can only be done with extensive simulations, beyond
the scope of this work. Nonetheless, the signal yield of our model sets a
reasonable target precision. Given the theoretical motivation of our model,
and how generic its signal yield is for all models of electroweakly-generated
LLPs, we hope our results in the upcoming section inform the on-going studies
in developing the detector design.

###  5.3 Charged Tracks and Displaced Leptons Signal

In this section we calculate the signal yield of our model on the LLP part of
the parameter space. We will show that our model gives rise to many events
across the entire kinetically-accessible part of the parameter space. We will,
in particular, highlight a double-peak feature in the distribution of events
in the detector, and argue that it is a shared signature of all electroweakly-
charged LLPs at a high energy MuC.

The kinematic quantities $\beta\gamma$ and $\eta$ directly enter the
calculation of the trajectory of LLPs inside the detector. In LLP searches we
are in particular interested in the transverse direction away from the beam,
$L$. It can be shown that for realistic magnetic field values, the track
curvature for $m_{\phi}\gtrsim 1$ TeV is very small and will not affect the
LLP observables under study here.888For precise measurements of some LLP
observables, such as closest approach of the track to the primary vertex, one
needs to include the effect of the magnetic field. However, this is not the
case for the displaced leptons and charged tracks signals studied in this
work. Thus, we neglect the effect of the magnetic field and assume the tracks
move on straight lines after their production.

We can show that as a function of time, the transverse distance from the beam
at time $t$ (in the lab frame) is given by

$L(t,\beta,\eta)=\frac{\beta t}{\cosh\eta}.$ (5.1)

Using this we can calculate the transverse displacement for every event at
$t=\gamma\tau_{\phi}$. The time $t=\gamma\tau_{\phi}$, is the average time at
which the $\phi$ particle decays to a charged lepton and missing energy,
giving rise to a displaced lepton.

In Figure 13 we show the distribution of transverse displacement $L$ at
$t=\gamma\tau_{\phi}$ for a few different points in our parameter space. The
sharp peak in large $L/\tau_{\phi}$ values in the distributions corresponds to
events from the DY channel that all have a very large $\beta\gamma$ value. We
also find that for smaller mediator masses the distribution stretches to
larger values of $L/\tau_{\phi}$, giving rise to more displaced leptons
further away from the beam.

Figure 13: Distribution of events in $\phi$’s transverse displacement at
$t=\gamma\tau_{\phi}$ normalized by the lifetime, $L/\tau_{\phi}$ for
different mediator masses at a MuC with center of mass energy of $\sqrt{s}=3$
TeV (left) or $\sqrt{s}=10$ TeV (right). The bin sizes are 6 cm/ns. The shape
of the $L/\tau_{\phi}$ can be inferred from distributions in Figure 11 and Eq.
(5.1). For different points in the parameter space the mediator lifetime
$\tau_{\phi}$ can be read off of Figure 4; combined with these plots, we can
then find the transverse displacement for events at any point in the parameter
space for every DM mass. The double-peak feature in the right figure can give
rise to interesting LLP signals, as elaborated in the text.

The $L/\tau_{\phi}$ distribution exhibits and interesting double-peak feature;
the peak at large values arises from the DY production while the peak at lower
values is an artifact of the $\eta$ dependence in Eq. (5.1) and arise from VBF
channel. Since these production channels universally apply to all
electroweakly-charged LLPs, we expect a similar double-peak feature in the
distribution of $L/\tau_{\phi}$ in all such models.999Once the muon PDF is
properly included in the simulations, we expect this peak to spread out to
some extent. Nonetheless, the sharp peak at $x\rightarrow 1$ in this PDF [96,
97, 98] still will give rise to a peak in the $\beta\gamma$ and in
$L/\tau_{\phi}$ distributions, maintaining this double-peak feature.  The main
irreducible background for LLP searches come from secondary interactions of
particles with the detector material and this double-peak feature may not be
manifested in this background. As a result, this feature could be a smoking
gun signature of electroweakly-charged LLPs at a MuC and could be used for an
efficient background rejection.

Figure 14: Distribution of events in the $\eta-\beta\gamma$ plane for
different $\phi$ masses at a 3 (10) TeV MuC on left (right). The bin sizes are
0.5 for $\eta$. For $\beta\gamma$, the bin size is 0.05 except for
$\sqrt{s}=10$ TeV and $m_{\phi}\leq 1$ TeV, where we use 0.2. Contours of
constant $L/\tau_{\phi}$ [cm/ns] are shown as well (see Eq. (5.1)), see the
text for further explanation on the shape of the contours. The cluster of
events at maximum $\beta\gamma$ values are from the DY production channel. We
find that majority of events from the DY and the VBF channel appear with
different kinematics.

In Figure 14 we show the joint distribution of events on the plane of
$\eta-\beta\gamma$ for a few different mediator masses. In each plot we can
see a cluster of events at highest physically possible value of $\beta\gamma$
which correspond to the DY-generated events. The VBF-generated events have a
wider distribution in both $\eta$ and $\beta\gamma$, as expected.

We also show contours of constant $L/\tau_{\phi}$ on the $\beta\gamma-\eta$
plane in Figure 14. The shape of these contours can be inferred from Eq.
(5.1). For a fixed value of $\tau_{\phi}$, these contours tell us at what
transverse distances the $\phi$ particle will (on average) decay and inform
the signal yield of our model. For $\tau_{\phi}=1$ ns these contours show $L$
in cm, while for larger lifetimes, the contours of constant $L$ move down on
the plots. We can repeat this calculation for all mediator and DM masses to
calculate the rate for displaced leptons appearing in any transverse distance
segment of the detector, inclusive over other kinematic variables such as
$\eta$ and $\beta\gamma$. The rate for displaced leptons appearing at a given
transverse distance is our model’s main signature in the LLP region.

The transverse distance can be mapped to detector components as defined in
Table 1. We show the average number of displaced leptons in the barrel region
of each detector component on our model parameter space in Figure 15 (Figure
16) for $\sqrt{s}=3$ TeV ($\sqrt{s}=10$ TeV).101010In the compressed region,
$m_{\phi}\approx m_{\chi}$, the daughter lepton from $\phi$ decay will be soft
and will typically escape detection, giving rise to a disappearing track if
they decay inside the detector. However, the lifetime $\tau_{\phi}$ is very
long in this region of the parameter space (see Figure 4), so our signal will
still be a heavy stable charged track. The disappearing track signals for LLPs
that decay inside a MuC are studied extensively in Ref. [103]. Endcap regions
are neglected since, similar to LHC, it is reasonable to expect more SM
background in endcaps with large track masses that are difficult to
distinguish from LLP signals, and that $dE/dx$ measurements in the endcaps are
slightly less accurate compared to their barrel region counterparts. Our
analysis can be straightforwardly repeated for the endcap regions as well.

Figure 15: Average rate of displaced leptons in different barrel regions of
the detector for a 3 TeV MuC, assuming the tentative design of Table 1, as
well as average number of detector-stable charged tracks. The DM Yukawa
coupling to the mediator is chosen to get the right relic abundance today. The
left (right) cluster corresponds to events from the initial DY (VBF) channel
of Figure 5, see the text for further explanation about the shape of contours.
The gray region corresponds to $m_{\chi}>m_{\phi}$ and is not
phenomenologically viable. Our model’s signal yield can serve as a benchmark
in studies of tolerable systematics in LLP searches at a MuC. The region below
the orange (blue) line is already ruled out by the LHC search in Ref. [87]
(Ref. [88]).

Figure 16: Similar to Figure 15, but for a MuC with center of mass energy of
$\sqrt{s}=10$ TeV. We again find that virtually the entire kinematically-
accessible part of the parameter space gives rise to large numbers of charged
tracks and displaced leptons. The region below the orange (blue) line is
already ruled out by the LHC search in Ref. [87] (Ref. [88]).

In each detector component, the distribution of events clearly divides up into
two separate regions. The upper cluster of events - at lower $m_{\chi}$ and
higher $m_{\phi}$ values - corresponds to events originating from initial DY
channel, while the other cluster of events come from the initial VBF channel.
These two separate clusters of events are a manifestation of the previously
explained “double-peak” feature seen in Figure 13.

The features of the signal regions in Figures 15 and 16 can be understood
through the shape of the $L$ contours in Figure 14. Recall that for a fixed
value of $\tau_{\phi}$, the contours in Figure 14 tell us at what transverse
distances the $\phi$ particle will (on average) decay, $i.e.$
$L(t=\gamma\tau_{\phi})$. At very short lifetimes, the contours associated
with each barrel region in Figure 14 are above the entire distribution of
events. As the lifetime increases, contours of constant $L$ move down on
Figure 14 (see Figure 4 for contours of constant lifetime). For fixed
$m_{\phi}$, $\tau_{\phi}$ increases with higher $m_{\chi}$ values, and so we
find that the events at low $m_{\chi}$ are predominantly from the DY channel
(as they are on top of 2D histograms in Figure 14). As we increase $m_{\chi}$
(and $\tau_{\phi}$), the VBF channel increasingly contributes. Since
distributions of DY- and VBF-generated events peak at different $\beta\gamma$
values in Figure 14, we find two separate cluster of events in each detector
component in Figures 15 and 16. The slope of each cluster of events is
inherited from contours of constant lifetime on the mass plane, see Figure 4.

The last panel in both Figures 15 and 16 corresponds to tracks that are
expected to go outside the detectors without decaying. We only include the
tracks with small enough $\eta$s that go through all barrel regions, i.e.
$\eta\lesssim 0.61$, see Table 1. These panels clearly show a large number of
such stable tracks are expected on a large part of the viable parameter space,
allowing us to look for the signals of this model on virtually the entire
kinematically accessible parameter space.

Our results should be contrasted to the current bounds from LHC on Figure 4.
Putting the numbers from all panels of Figures 15 and 16 together, we find
that we have at least $10^{3}$ events in at least one part of the detector for
almost the entire kinematically-accessible parts of our parameter space, see
Appendix C. Once a reasonable degree of systematics is achieved, future high
energy MuCs can sift through a much larger part of the parameter space than
LHC. This clearly underscores the supremacy of a high energy MuC in searching
for this DM model. Crucial in this supremacy are the stable tracks that go
through the entire detector without decaying, see the last panel of Figures
15–16. For most of the parameter space this is the dominant LLP signal and can
be used for probing the parameter space inaccessible by displaced lepton
signals in the detector. The double-peak feature, which in turn is an artifact
of two distinct production channels, suggests for many points in the parameter
space we can have many displaced leptons in more than one component of the
detector.

We should emphasize that Figures 15 and 16 show the average number of decays
and displaced leptons (or stable charged tracks). The true number of events in
each barrel region will be drawn from a Poisson distribution around this
average. Further details about kinematics of our model in the LLP region are
included in Appendix C.

##  6 Conclusions

We initiated a detailed study of the freeze-in flavored DM setup in which
flavored DM coupling to the slepton-like mediator is very small. We showed
that the interplay of direct freeze-in and the mediator freeze-out gives rise
to an interesting abundance calculation, such that the dark matter mass should
be below a few TeV in order to avoid overclosing the universe. Combined with
the existing collider and astrophysical bounds, this gives rise to a bounded
viable parameter space for this model.

We also studied signals of our model at a future high energy Muon Collider
(MuC). As a result of the feeble Yukawa coupling of DM, all signals of the
model in a MuC originate from the on-shell pair production of the mediator.
The range of mediator masses and lifetimes in the viable parameter space is
such that we can have prompt, as well as various long-lived particle signals
at a MuC. We divided our parameter space into two broad ranges according to
the mediator’s lifetime, namely prompt and LLP regions.

In the prompt decay region, we proposed a rudimentary analysis and identified
the target systematics that ought to be achieved in order to discover our
model at a MuC. Our analysis uses cuts on invariant leptonic system mass
$m_{\ell\ell}$, missing transverse mass parameter $M_{T2}$, and the leptonic
system transverse momentum $p_{T,\ell\ell}$. We found that in the part of the
parameter space with prompt mediator decay, if systematic uncertainties are
comparable or smaller than SM background, our proposed analysis can discover
this model at effectively the entire kinematically-accessible parameter space.

For longer lifetimes, we will have a charged track in the detector which,
depending on the lifetime, can either decay to a lepton and missing energy or
decay outside the detector. While the SM background can be reduced away using
quantities such as time-of-flight or anomalous ionization rate, detector
response can give rise to an irreducible background. A proper study of these
effects requires extensive simulations and is left for future works. In the
absence of concrete simulations, here we reported the signal yield and various
kinematics distributions. In particular, we showed that the displaced lepton
distribution has a double-peak signature that could potentially be used for
efficient background rejection. This feature appears as two large sets of
signal events, produced through the two distinct channels (DY and VBF), in
different detector components.

Our work can be extended in many directions. Foremost among them is a detailed
study of the detector design that enables detailed studies of the background.
Given the great reach of future MuCs in our model’s parameter space, the
universality of its production rate, and its theoretical motivation, we
believe our model can serve as a motivated benchmark target for future studies
on development of the detector by setting a target for systematics of such
works.

Our analysis can also be repeated for other flavored DM models, i.e. mediators
with same charges as other SM fermions. All these setups have the added bonus
that the artificial alignment of the DM Yukawa and SM Yukawas in flavored DM
models, required for evading bounds from flavor-changing neutral currents and
lepton flavor violation, are completely avoided thanks to the feeble DM
coupling. While the relic abundance calculation resembles the calculation
here, we expect a different collider phenomenology for some of them thanks to
their different production and decay channels.

Furthermore, as alluded to earlier, the muon PDF is not yet included in
MadGraph5, thus we neglected its effect in the DY production channel. With the
ever-growing interest in BSM signals at a MuC, it is interesting to work on
proper embedding of this PDF in the event-generation pipeline.

It will also be interesting to consider signals of such models in other
experiments, such as various astrophysical or direct detection searches. In
particular, there is a large part of the parameter space inaccessible to
colliders in foreseeable future, motivating searches in the aforementioned
complementary fronts.

High energy MuCs are strongly motivated for their ability in probing different
models of the Higgs boson UV-completion, flavorful new physics, and DM models.
Our study outlines a simple and interesting DM model with unique signatures
that can serve as a target for searches at a high energy MuC and inform its
detector design.

## Acknowledgement

We thank Rodolfo Capdevilla, Spencer Chang, Jochen Heinrich, Samuel Homiller,
Graham Kribs, and Ben Lillard for helpful discussions. We are especially
grateful for numerous illuminating discussions with Laura Jeanty about LLPs.
We also thank Samuel Homiller for constructive comments on the draft. The work
of P.A. is supported by the U.S. Department of Energy under grant number DE-
SC0011640. A.R. and T.-T.Y. are supported in part by NSF CAREER grant
PHY-1944826.

The code to produce the plots and results in this paper is available here. We
also use the following software:

Julia [107], DataFrames.jl [108], CSV.jl [109], DifferentialEquations.jl
[110], HCubature.jl [111], Optim.jl [112], PythonCall.jl [113], QuadGK.jl
[114], Roots.jl [115], matplotlib [116].

##  Appendix A More Details on Relic Abundance Calculation

The relic abundance is directly related to the final $\chi$ yield,
$Y_{\chi}^{\infty}$,

$\displaystyle\Omega_{\chi}$
$\displaystyle=2m_{\chi}Y_{\chi}^{\infty}\frac{T_{0}^{3}}{30}\frac{8\pi
G}{3H_{0}^{2}}$ (A.1)

where $T_{0}=2.7255$ K and $H_{0}=100h$ km s-1 Mpc-1 are the temperature and
Hubble constant today, respectively, and $G$ is the gravitational constant.
The factor of $2$ accounts for $\chi$ and $\bar{\chi}$ both contributing to
the final relic abundance. To be fully accurate, one should run the Boltzmann
equation solver out to a sufficiently long time to obtain $Y_{\chi}^{\infty}$.
However, we can also estimate this by noticing that $Y_{\phi}$ follows its
freeze-out value until decaying, and $Y_{\chi}$ follows its freeze-in value
(assuming $Y_{\phi}$ is in equilibrium) until $\phi$ decays. This information
leads to the approximation that

$\displaystyle Y_{\chi}^{\infty}\approx
Y_{\chi}^{\textrm{FI}}+Y_{\phi}^{\textrm{FO}}$ (A.2)

This approximation is useful because $Y_{\phi}^{\textrm{FO}}$ only depends on
$\phi$’s mass, its couplings to the SM are fixed by its charge, and we can
analytically calculate $Y_{\chi}^{\textrm{FI}}$. We have checked that for the
majority of the parameter space, this approximation has less than $1\%$ error
in the abundance calculation.

Starting with $\chi$’s Boltzmann equation we set
$Y_{\phi}=Y_{\phi,\textrm{EQ}}$ and start from $Y_{\chi}=0$. The Boltzmann
equation for $\chi$ will be [73]

$\displaystyle\frac{dY_{\chi}}{dx}$
$\displaystyle=\frac{x^{3}}{H(m_{\phi})}\frac{g_{\phi}3\Gamma_{\phi\to\ell\chi}}{2\pi^{2}}K_{1}(x),$
(A.3)

and we can simply integrate this from $x=0$ to $x=\infty$ to find

$\displaystyle
Y_{\chi}^{\textrm{FI}}=\frac{9g_{\phi}\Gamma_{\phi\to\ell\chi}}{4\pi
H(m_{\phi})},$ (A.4)

with the decay width for $\phi\to\ell\chi$ given by

$\displaystyle\Gamma_{\phi\to\ell\chi}$
$\displaystyle=\frac{\lambda^{2}m_{\phi}}{16\pi}\left(1-\frac{m_{\chi}^{2}}{m_{\phi}^{2}}\right)^{2}.$
(A.5)

We can plug in the approximation to the relic abundance to get

$\displaystyle\Omega_{\chi}$
$\displaystyle=2m_{\chi}\left(Y_{\chi}^{\textrm{FI}}+Y_{\phi}^{\textrm{FO}}\right)\frac{T_{0}^{3}}{30}\frac{8\pi
G}{3H_{0}^{2}}$ (A.6)

To find the particular $\lambda$ that gives the correct relic abundance
$\Omega_{\chi}h^{2}\approx 0.12$, we insert our expression for
$Y_{\chi}^{\textrm{FI}}$ and rearrange to get the Yukawa coupling that gives
rise to the right relic abundance today

$\displaystyle\lambda^{\prime}$
$\displaystyle=\frac{8\pi}{3(1-m_{\chi}^{2}/m_{\phi}^{2})}\sqrt{\frac{H(m_{\phi})}{g_{\phi}m_{\chi}m_{\phi}}}\sqrt{\frac{45\Omega_{\chi}H_{0}^{2}}{8\pi
GT_{0}^{3}}-m_{\chi}Y_{\phi}^{\textrm{FO}}}.$ (A.7)

We plot this quantity in Figure 3.

##  Appendix B More Details on the Prompt Search

In this appendix we provide further results of our MadGraph5 simulations that
inform the cuts we used in our prompt region search.

To better determine what cuts can optimize the reach of a MuC collider, in
Figure 17 (Figure 18) we show 2D histograms of the ratio of the signal ($S$)
to square root of the background ($B$) for a few different $\phi$ masses and
for center of mass energy of 3 TeV (10 TeV). These histograms inform us about
what lower bounds on each kinematic parameter yields the best reach for any
given mediator mass and collider’s center of mass energy.

Figure 17: Contours of significance of the signal (signal count $S$ over
square root of background $\sqrt{B}$) for different minimum values on each
observable. We show the effect of cuts on events from two different $\phi$
masses, 1 TeV (left) and 1.3 TeV (right), at a MuC with center of mass energy
of $\sqrt{s}=3$ TeV. In each panel we show the effect of varying two cuts,
while fixing the third one to the value that maximizes the significance. We
find large significances are attainable. These figures inform the cuts used in
definition of the signal regions.

Figure 18: Contours of significance of the signal (signal count $S$ over
square root of background $\sqrt{B}$) for different minimum values on each
observable. This is the same as Figure 17, but for $\phi$ masses 1 TeV (left)
and 4.5 TeV (right), and at a MuC with center of mass energy of $\sqrt{s}=10$
TeV. We find large significances are attainable. These figures inform the cuts
used in definition of the signal regions.

Based on these results, we can define a set of signal regions that together
are sensitive to all values of $m_{\phi}$. The definition of signal regions we
use and their background $B$ is included in Figure 19.

Figure 19: Signal regions (SRs) used in our proposed search for $\sqrt{s}=3$
TeV (left) and for $\sqrt{s}=10$ TeV (right). The cuts used on each kinematic
variable and the SM background in each region is reported. Determination of
the cuts are informed by the most sensitive cuts on various extreme masses in
the parameter space from Eqs. (4.1)–(4.4).

To further illustrate the effect of each cut as a function of $\phi$ mass, in
Figure 20 we show the signal yield and the background for four sample signal
regions with very similar cuts. These figures indicate that an $M_{T2}$ cut is
the lower bound on the mediator masses that populate a signal region, as
expected [100]. The $p_{T,\ell\ell}$ and $m_{\ell\ell}$ cuts also reduce the
SM and total signal yield in a signal region. Increasing the $p_{T,\ell\ell}$
cut shifts the distribution to higher $m_{\phi}$ masses as well.

Figure 20: Signal yield and the SM background of four sample signal regions
with similar cuts, for different mediator masses. For each mass if a signal
region has $S\geqslant 5\sqrt{B}$ signal events it can give rise to a
discovery (neglecting systematics). We find that the $M_{T2}$ cut controls the
lowest mediator mass each signal region is sensitive too. Different
$m_{\ell\ell}$ cuts affect the total number of events, while leaving the
overall dependence on mediator mass unchanged. Comparing the signal regions
with similar cuts we find that increasing the $p_{T,\ell\ell}$ cut shifts the
signal distribution to higher $\phi$ masses.

##  Appendix C More Details on the LLP Search

For completeness, in this appendix we include the 2D histograms of events in
different kinematic observables in the LLP search. In particular, in Figure 21
(Figure 22) we show the 2D histograms of the distribution of events in the
transverse distance $L(t=\gamma\tau_{\phi})/\tau_{\phi}$ and their
pseudorapidity $\eta$ (their $\beta\gamma$), for a few mass points in the
parameter space. The double-peak feature in Figures 11-14 are still visible
here as well. The figure also shows that the smaller the mediator mass, the
more separated in $L/\tau_{\phi}$ the two peaks in the distribution.

Figure 21: Joint distribution of events in $\eta-L/\tau_{\phi}$ for three
different mediator masses (different rows) at a 3 TeV MuC (left) or a 10 TeV
MuC (right). The bin sizes are all 0.5 for $\eta$ and 1.5 cm/ns for
$L/\tau_{\phi}$ (except for $\sqrt{s}=10$ TeV, $m_{\phi}=0.5$ and $1$ TeV
where they are 6 cm/ns). Events at the largest $L/\tau_{\phi}$ values are due
to the DY initial channel. The double-peak feature in the distribution (at
large and small $L$ values) is manifested as well.

Figure 22: Joint distribution of events in $\beta\gamma-L/\tau_{\phi}$ for
three different mediator masses (different rows) at a 3 TeV MuC (left) or a 10
TeV MuC (right). Bin sizes are ($L/\tau_{\phi}$, $\beta\gamma$) = (1.5 cm/ns,
0.05) (except $\sqrt{s}=10$ TeV, $m_{\phi}=0.5$ and $1$ TeV where they are
($L/\tau_{\phi}$, $\beta\gamma$) = (6 cm/ns, 0.2)). Events at the largest
$L/\tau_{\phi}$ values are due to the DY initial channel. The double-peak
feature in the distribution (at large and small $L$ values) is manifested as
well.

Figure 21 also reiterates the concentration of events at small $|\eta|$
values. This allows us to introduce stringent cuts on this variable to cut
down the SM background, while preserving majority of signal events. The fact
that the DY-generated events are more concentrated at smaller $|\eta|$ values
is also an underlying cause of the separation of the cluster of events from
different channels in Figures 15–16.

In Figure 23 we show the maximum number of decays across different detector
segments. To make these plots, for each point in the parameter space we simply
report the maximum event count from among all panels in Figures 15–16. We find
that for most points in the parameter space we have at least $10^{3}-10^{4}$
events in at least one part of the detector (including events that are
detector stable and give rise to stable charged tracks).

Figure 23: Maximum rate of displaced leptons across different detector regions
or stable charged tracks for every point in the parameter space in a 3 TeV
(left) and a 10 TeV (right) MuC. The region below the orange (yellow) line is
already ruled out by the LHC search in Ref. [87] (Ref. [88]).

## References

  * [1] Muon Collider Collaboration, D. Stratakis et al., “A Muon Collider Facility for Physics Discovery,” arXiv:2203.08033 [physics.acc-ph].
  * [2] C. M. Ankenbrandt et al., “Status of muon collider research and development and future plans,” Phys. Rev. ST Accel. Beams 2 (1999) 081001, arXiv:physics/9901022.
  * [3] M.-H. Wang, Y. Nosochkov, Y. Cai, and M. Palmer, “Design of a 6 tev muon collider,” Journal of Instrumentation 11 no. 09, (Sep, 2016) P09003. https://dx.doi.org/10.1088/1748-0221/11/09/P09003.
  * [4] M. Boscolo, J.-P. Delahaye, and M. Palmer, “The future prospects of muon colliders and neutrino factories,” Rev. Accel. Sci. Tech. 10 no. 01, (2019) 189–214, arXiv:1808.01858 [physics.acc-ph].
  * [5] D. Neuffer and V. Shiltsev, “On the feasibility of a pulsed 14 TeV c.m.e. muon collider in the LHC tunnel,” JINST 13 no. 10, (2018) T10003, arXiv:1811.10694 [physics.acc-ph].
  * [6] 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].
  * [7] V. Shiltsev and F. Zimmermann, “Modern and Future Colliders,” Rev. Mod. Phys. 93 (2021) 015006, arXiv:2003.09084 [physics.acc-ph].
  * [8] Muon Collider Collaboration, J. de Blas et al., “The physics case of a 3 TeV muon collider stage,” arXiv:2203.07261 [hep-ph].
  * [9] Muon Collider Collaboration, S. Jindariani et al., “Promising Technologies and R&D Directions for the Future Muon Collider Detectors,” arXiv:2203.07224 [physics.ins-det].
  * [10] Muon Collider Collaboration, N. Bartosik et al., “Simulated Detector Performance at the Muon Collider,” arXiv:2203.07964 [hep-ex].
  * [11] K. M. Black et al., “Muon Collider Forum Report,” arXiv:2209.01318 [hep-ex].
  * [12] C. Accettura et al., “Towards a muon collider,” Eur. Phys. J. C 83 no. 9, (2023) 864, arXiv:2303.08533 [physics.acc-ph].
  * [13] E. Eichten and A. Martin, “The Muon Collider as a $H/A$ Factory,” Phys. Lett. B 728 (2014) 125–130, arXiv:1306.2609 [hep-ph].
  * [14] N. Chakrabarty, T. Han, Z. Liu, and B. Mukhopadhyaya, “Radiative Return for Heavy Higgs Boson at a Muon Collider,” Phys. Rev. D 91 no. 1, (2015) 015008, arXiv:1408.5912 [hep-ph].
  * [15] D. Buttazzo, D. Redigolo, F. Sala, and A. Tesi, “Fusing Vectors into Scalars at High Energy Lepton Colliders,” JHEP 11 (2018) 144, arXiv:1807.04743 [hep-ph].
  * [16] 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].
  * [17] T. Han, D. Liu, I. Low, and X. Wang, “Electroweak couplings of the Higgs boson at a multi-TeV muon collider,” Phys. Rev. D 103 no. 1, (2021) 013002, arXiv:2008.12204 [hep-ph].
  * [18] P. Bandyopadhyay and A. Costantini, “Obscure Higgs boson at Colliders,” Phys. Rev. D 103 no. 1, (2021) 015025, arXiv:2010.02597 [hep-ph].
  * [19] D. Buttazzo, R. Franceschini, and A. Wulzer, “Two Paths Towards Precision at a Very High Energy Lepton Collider,” JHEP 05 (2021) 219, arXiv:2012.11555 [hep-ph].
  * [20] W. Liu and K.-P. Xie, “Probing electroweak phase transition with multi-TeV muon colliders and gravitational waves,” JHEP 04 (2021) 015, arXiv:2101.10469 [hep-ph].
  * [21] T. Han, S. Li, S. Su, W. Su, and Y. Wu, “Heavy Higgs bosons in 2HDM at a muon collider,” Phys. Rev. D 104 no. 5, (2021) 055029, arXiv:2102.08386 [hep-ph].
  * [22] R. Franceschini and M. Greco, “Higgs and BSM Physics at the Future Muon Collider,” Symmetry 13 no. 5, (2021) 851, arXiv:2104.05770 [hep-ph].
  * [23] M. Chiesa, B. Mele, and F. Piccinini, “Multi Higgs production via photon fusion at future multi-TeV muon colliders,” arXiv:2109.10109 [hep-ph].
  * [24] 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].
  * [25] S. Spor, “Probe of the anomalous neutral triple gauge couplings in photon-induced collision at future muon colliders,” Nucl. Phys. B 991 (2023) 116198, arXiv:2207.11585 [hep-ph].
  * [26] H. Amarkhail, S. C. Inan, and A. V. Kisselev, “Probing anomalous $\gamma\gamma\gamma\gamma$ couplings at a future muon collider,” arXiv:2306.03653 [hep-ph].
  * [27] R. Capdevilla, D. Curtin, Y. Kahn, and G. Krnjaic, “Discovering the physics of $(g-2)_{\mu}$ at future muon colliders,” Phys. Rev. D 103 no. 7, (2021) 075028, arXiv:2006.16277 [hep-ph].
  * [28] D. Buttazzo and P. Paradisi, “Probing the muon $g-2$ anomaly with the Higgs boson at a muon collider,” Phys. Rev. D 104 no. 7, (2021) 075021, arXiv:2012.02769 [hep-ph].
  * [29] W. Yin and M. Yamaguchi, “Muon g-2 at a multi-TeV muon collider,” Phys. Rev. D 106 no. 3, (2022) 033007, arXiv:2012.03928 [hep-ph].
  * [30] G.-y. Huang, F. S. Queiroz, and W. Rodejohann, “Gauged $L_{\mu}{-}L_{\tau}$ at a muon collider,” Phys. Rev. D 103 no. 9, (2021) 095005, arXiv:2101.04956 [hep-ph].
  * [31] 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 no. 1, (2022) 015028, arXiv:2101.10334 [hep-ph].
  * [32] N. Chen, B. Wang, and C.-Y. Yao, “The collider tests of a leptophilic scalar for the anomalous magnetic moments,” arXiv:2102.05619 [hep-ph].
  * [33] G.-y. Huang, S. Jana, F. S. Queiroz, and W. Rodejohann, “Probing the RK(*) anomaly at a muon collider,” Phys. Rev. D 105 no. 1, (2022) 015013, arXiv:2103.01617 [hep-ph].
  * [34] 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].
  * [35] P. Bandyopadhyay, A. Karan, R. Mandal, and S. Parashar, “Distinguishing signatures of scalar leptoquarks at hadron and muon colliders,” Eur. Phys. J. C 82 no. 10, (2022) 916, arXiv:2108.06506 [hep-ph].
  * [36] S. Qian, C. Li, Q. Li, F. Meng, J. Xiao, T. Yang, M. Lu, and Z. You, “Searching for heavy leptoquarks at a muon collider,” JHEP 12 (2021) 047, arXiv:2109.01265 [hep-ph].
  * [37] S. Homiller, Q. Lu, and M. Reece, “Complementary signals of lepton flavor violation at a high-energy muon collider,” JHEP 07 (2022) 036, arXiv:2203.08825 [hep-ph].
  * [38] A. Azatov, F. Garosi, A. Greljo, D. Marzocca, J. Salko, and S. Trifinopoulos, “New physics in $b\rightarrow s\mu\mu$: FCC-hh or a muon collider?,” JHEP 10 (2022) 149, arXiv:2205.13552 [hep-ph].
  * [39] J.-L. Yang, C.-H. Chang, and T.-F. Feng, “The leptonic di-flavor and di-number violation processes at high energy $\mu^{\pm}\mu^{\pm}$ colliders,” arXiv:2302.13247 [hep-ph].
  * [40] W. Altmannshofer, S. A. Gadam, and S. Profumo, “Probing New Physics with $\mu^{+}\mu^{-}\to bs$ at a Muon Collider,” arXiv:2306.15017 [hep-ph].
  * [41] S. Jana and S. Klett, “Muonic Force and Neutrino Non-Standard Interactions at Muon Colliders,” arXiv:2308.07375 [hep-ph].
  * [42] N. Ghosh, S. K. Rai, and T. Samui, “Search For a Leptoquark and Vector-like Lepton in a Muon Collider,” arXiv:2309.07583 [hep-ph].
  * [43] T. Han, Z. Liu, L.-T. Wang, and X. Wang, “WIMPs at High Energy Muon Colliders,” Phys. Rev. D 103 no. 7, (2021) 075004, arXiv:2009.11287 [hep-ph].
  * [44] S. Bottaro, A. Strumia, and N. Vignaroli, “Minimal Dark Matter bound states at future colliders,” JHEP 06 (2021) 143, arXiv:2103.12766 [hep-ph].
  * [45] 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].
  * [46] A. Jueid and S. Nasri, “Lepton portal dark matter at muon colliders: Total rates and generic features for phenomenologically viable scenarios,” Phys. Rev. D 107 no. 11, (2023) 115027, arXiv:2301.12524 [hep-ph].
  * [47] M. Belfkir, A. Jueid, and S. Nasri, “Boosting dark matter searches at muon colliders with Machine Learning: the mono-Higgs channel as a case study,” arXiv:2309.11241 [hep-ph].
  * [48] L. Di Luzio, R. Gröber, and G. Panico, “Probing new electroweak states via precision measurements at the LHC and future colliders,” JHEP 01 (2019) 011, arXiv:1810.10993 [hep-ph].
  * [49] 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].
  * [50] W. Liu, K.-P. Xie, and Z. Yi, “Testing leptogenesis at the LHC and future muon colliders: A Z’ scenario,” Phys. Rev. D 105 no. 9, (2022) 095034, arXiv:2109.15087 [hep-ph].
  * [51] H. Al Ali et al., “The muon Smasher’s guide,” Rept. Prog. Phys. 85 no. 8, (2022) 084201, arXiv:2103.14043 [hep-ph].
  * [52] M. Casarsa, M. Fabbrichesi, and E. Gabrielli, “Monochromatic single photon events at the muon collider,” Phys. Rev. D 105 no. 7, (2022) 075008, arXiv:2111.13220 [hep-ph].
  * [53] Y. Bao, J. Fan, and L. Li, “Electroweak ALP searches at a muon collider,” JHEP 08 (2022) 276, arXiv:2203.04328 [hep-ph].
  * [54] S. C. İnan and A. V. Kisselev, “Probe of axion-like particles in vector boson scattering at a muon collider,” J. Phys. G 50 no. 10, (2023) 105002, arXiv:2207.03325 [hep-ph].
  * [55] 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.
  * [56] M. Chen and D. Liu, “Top Yukawa Coupling at the Muon Collider,” arXiv:2212.11067 [hep-ph].
  * [57] T. H. Kwok, L. Li, T. Liu, and A. Rock, “Searching for Heavy Neutral Leptons at A Future Muon Collider,” arXiv:2301.05177 [hep-ph].
  * [58] S. C. İnan and A. V. Kisselev, “Probe of a Randall-Sundrum-like model from muon pair production at high energy muon collider,” arXiv:2301.08585 [hep-ph].
  * [59] T. A. Chowdhury, A. Jueid, S. Nasri, and S. Saad, “Probing Zee-Babu states at Muon Colliders,” arXiv:2306.01255 [hep-ph].
  * [60] M. Belfkir, T. A. Chowdhury, and S. Nasri, “Doubly-charged scalars of the Minimal Left-Right Symmetric Model at Muon Colliders,” arXiv:2307.16111 [hep-ph].
  * [61] S. Chigusa, S. Girmohanta, Y. Nakai, and Y. Zhang, “Aiming for Tops of ALPs with a Muon Collider,” arXiv:2310.11018 [hep-ph].
  * [62] C. Cesarotti, S. Homiller, R. K. Mishra, and M. Reece, “Probing New Gauge Forces with a High-Energy Muon Beam Dump,” Phys. Rev. Lett. 130 no. 7, (2023) 071803, arXiv:2202.12302 [hep-ph].
  * [63] C. Cesarotti and R. Gambhir, “The New Physics Case for Beam-Dump Experiments with Accelerated Muon Beams,” arXiv:2310.16110 [hep-ph].
  * [64] P. Agrawal, S. Blanchet, Z. Chacko, and C. Kilic, “Flavored Dark Matter, and Its Implications for Direct Detection and Colliders,” Phys. Rev. D 86 (2012) 055002, arXiv:1109.3516 [hep-ph].
  * [65] B. Batell, T. Lin, and L.-T. Wang, “Flavored Dark Matter and R-Parity Violation,” JHEP 01 (2014) 075, arXiv:1309.4462 [hep-ph].
  * [66] P. Agrawal, B. Batell, D. Hooper, and T. Lin, “Flavored Dark Matter and the Galactic Center Gamma-Ray Excess,” Phys. Rev. D 90 no. 6, (2014) 063512, arXiv:1404.1373 [hep-ph].
  * [67] P. Agrawal, M. Blanke, and K. Gemmler, “Flavored dark matter beyond Minimal Flavor Violation,” JHEP 10 (2014) 072, arXiv:1405.6709 [hep-ph].
  * [68] P. Agrawal, Z. Chacko, C. Kilic, and C. B. Verhaaren, “A Couplet from Flavored Dark Matter,” JHEP 08 (2015) 072, arXiv:1503.03057 [hep-ph].
  * [69] P. Agrawal, Z. Chacko, E. C. F. S. Fortes, and C. Kilic, “Skew-Flavored Dark Matter,” Phys. Rev. D 93 no. 10, (2016) 103510, arXiv:1511.06293 [hep-ph].
  * [70] P. Agrawal, C. Kilic, S. Swaminathan, and C. Trendafilova, “Secretly Asymmetric Dark Matter,” Phys. Rev. D 95 no. 1, (2017) 015031, arXiv:1608.04745 [hep-ph].
  * [71] N. Desai, C. Kilic, Y.-P. Yang, and T. Youn, “Suppressed flavor violation in Lepton Flavored Dark Matter from an extra dimension,” Phys. Rev. D 101 (2020) 075043, arXiv:2001.00720 [hep-ph].
  * [72] H. Acaroğlu, P. Agrawal, and M. Blanke, “Lepton-flavoured scalar dark matter in Dark Minimal Flavour Violation,” JHEP 05 (2023) 106, arXiv:2211.03809 [hep-ph].
  * [73] L. J. Hall, K. Jedamzik, J. March-Russell, and S. M. West, “Freeze-In Production of FIMP Dark Matter,” JHEP 03 (2010) 080, arXiv:0911.1120 [hep-ph].
  * [74] S. Chang, R. Edezhath, J. Hutchinson, and M. Luty, “Effective WIMPs,” Phys. Rev. D 89 no. 1, (2014) 015011, arXiv:1307.8120 [hep-ph].
  * [75] S. Chang, R. Edezhath, J. Hutchinson, and M. Luty, “Leptophilic Effective WIMPs,” Phys. Rev. D 90 no. 1, (2014) 015011, arXiv:1402.7358 [hep-ph].
  * [76] Y. Bai and J. Berger, “Fermion Portal Dark Matter,” JHEP 11 (2013) 171, arXiv:1308.0612 [hep-ph].
  * [77] Y. Bai and J. Berger, “Lepton Portal Dark Matter,” JHEP 08 (2014) 153, arXiv:1402.6696 [hep-ph].
  * [78] A. DiFranzo, K. I. Nagao, A. Rajaraman, and T. M. P. Tait, “Simplified Models for Dark Matter Interacting with Quarks,” JHEP 11 (2013) 014, arXiv:1308.2679 [hep-ph]. [Erratum: JHEP 01, 162 (2014)].
  * [79] M. Garny, A. Ibarra, S. Rydbeck, and S. Vogl, “Majorana Dark Matter with a Coloured Mediator: Collider vs Direct and Indirect Searches,” JHEP 06 (2014) 169, arXiv:1403.4634 [hep-ph].
  * [80] F. Elahi, C. Kolda, and J. Unwin, “UltraViolet Freeze-in,” JHEP 03 (2015) 048, arXiv:1410.6157 [hep-ph].
  * [81] G. D’Ambrosio, S. Chatterjee, R. Laha, and S. K. Vempati, “Freezing In with Lepton Flavored Fermions,” SciPost Phys. 11 (2021) 006, arXiv:2103.03886 [hep-ph].
  * [82] J. A. Evans and J. Shelton, “Long-Lived Staus and Displaced Leptons at the LHC,” JHEP 04 (2016) 056, arXiv:1601.01326 [hep-ph].
  * [83] E. W. Kolb and M. S. Turner, The Early Universe, vol. 69.  1990\.
  * [84] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [85] ATLAS Collaboration, G. Aad et al., “Search for direct pair production of sleptons and charginos decaying to two leptons and neutralinos with mass splittings near the W-boson mass in $\sqrt{s}$ = 13 TeV pp collisions with the ATLAS detector,” JHEP 06 (2023) 031, arXiv:2209.13935 [hep-ex].
  * [86] 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 no. 2, (2020) 123, arXiv:1908.08215 [hep-ex].
  * [87] ATLAS Collaboration, G. Aad et al., “Search for Displaced Leptons in $\sqrt{s}=13$ TeV $pp$ Collisions with the ATLAS Detector,” Phys. Rev. Lett. 127 no. 5, (2021) 051802, arXiv:2011.07812 [hep-ex].
  * [88] CMS Collaboration, V. Khachatryan et al., “Search for long-lived charged particles in proton-proton collisions at $\sqrt{s}=$ 13 TeV,” Phys. Rev. D 94 no. 11, (2016) 112004, arXiv:1609.08382 [hep-ex].
  * [89] DELPHI Collaboration, J. Abdallah et al., “Searches for supersymmetric particles in e+ e- collisions up to 208-GeV and interpretation of the results within the MSSM,” Eur. Phys. J. C 31 (2003) 421–479, arXiv:hep-ex/0311019.
  * [90] ALEPH Collaboration, A. Heister et al., “Search for charginos nearly mass degenerate with the lightest neutralino in e+ e- collisions at center-of-mass energies up to 209-GeV,” Phys. Lett. B 533 (2002) 223–236, arXiv:hep-ex/0203020.
  * [91] F. D’Eramo and A. Lenoci, “Lower mass bounds on FIMP dark matter produced via freeze-in,” JCAP 10 (2021) 045, arXiv:2012.01446 [hep-ph].
  * [92] MICE Collaboration, M. Bogomilov et al., “Demonstration of cooling by the Muon Ionization Cooling Experiment,” Nature 578 no. 7793, (2020) 53–59, arXiv:1907.08562 [physics.acc-ph].
  * [93] CLICdp Collaboration, N. Alipour Tehrani et al., “CLICdet: The post-CDR CLIC detector model,”.
  * [94] ILC Collaboration, G. Aarons et al., “ILC Reference Design Report Volume 4 - Detectors,” arXiv:0712.2356 [physics.ins-det].
  * [95] S. Knapen and S. Lowette, “A guide to hunting long-lived particles at the LHC,” arXiv:2212.03883 [hep-ph].
  * [96] T. Han, Y. Ma, and K. Xie, “High energy leptonic collisions and electroweak parton distribution functions,” Phys. Rev. D 103 no. 3, (2021) L031301, arXiv:2007.14300 [hep-ph].
  * [97] 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].
  * [98] 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].
  * [99] J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer, and T. Stelzer, “MadGraph 5 : Going Beyond,” JHEP 06 (2011) 128, arXiv:1106.0522 [hep-ph].
  * [100] C. G. Lester and D. J. Summers, “Measuring masses of semiinvisibly decaying particles pair produced at hadron colliders,” Phys. Lett. B 463 (1999) 99–103, arXiv:hep-ph/9906349.
  * [101] L. Lee, C. Ohm, A. Soffer, and T.-T. Yu, “Collider Searches for Long-Lived Particles Beyond the Standard Model,” Prog. Part. Nucl. Phys. 106 (2019) 210–255, arXiv:1810.12602 [hep-ph]. [Erratum: Prog.Part.Nucl.Phys. 122, 103912 (2022)].
  * [102] O. Buchmueller, A. De Roeck, K. Hahn, M. McCullough, P. Schwaller, K. Sung, and T.-T. Yu, “Simplified Models for Displaced Dark Matter Signatures,” JHEP 09 (2017) 076, arXiv:1704.06515 [hep-ph].
  * [103] 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].
  * [104] S. Bottaro, D. Buttazzo, M. Costa, R. Franceschini, P. Panci, D. Redigolo, and L. Vittorio, “Closing the window on WIMP Dark Matter,” Eur. Phys. J. C 82 no. 1, (2022) 31, arXiv:2107.09688 [hep-ph].
  * [105] S. Bottaro, D. Buttazzo, M. Costa, R. Franceschini, P. Panci, D. Redigolo, and L. Vittorio, “The last complex WIMPs standing,” Eur. Phys. J. C 82 no. 11, (2022) 992, arXiv:2205.04486 [hep-ph].
  * [106] J. Alimena et al., “Searching for long-lived particles beyond the Standard Model at the Large Hadron Collider,” J. Phys. G 47 no. 9, (2020) 090501, arXiv:1903.04497 [hep-ex].
  * [107] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, “Julia: A Fresh Approach to Numerical Computing,” SIAM Review 59 no. 1, (Jan., 2017) 65–98. https://epubs.siam.org/doi/10.1137/141000671. Publisher: Society for Industrial and Applied Mathematics.
  * [108] M. Bouchet-Valat and B. Kamiński, “Dataframes.jl: Flexible and fast tabular data in julia,” Journal of Statistical Software 107 no. 4, (2023) 1–32. https://www.jstatsoft.org/index.php/jss/article/view/v107i04.
  * [109] J. Quinn, M. Bouchet-Valat, N. Robinson, B. Kamiński, G. Newman, A. Stukalov, C. Vogt, cjprybol, T. Holy, A. Noack, T. Kelman, E. Davies, ExpandingMan, Ian, L. O. Hafner, M. Piibeleht, R. Finnegan, evalparse, A. M. Silberstein, A. Heimerson, P. Anthony Blaom, B. Lungwitz, B. König, C. de Graaf, C. Woodfield, D. Barton, D. Aluthge, E. Saba, F. Noronha, and kragol, “Juliadata/csv.jl: v0.10.11,” June, 2023.  https://doi.org/10.5281/zenodo.8004128.
  * [110] C. Rackauckas and Q. Nie, “Differentialequations.jl – a performant and feature-rich ecosystem for solving differential equations in julia,” The Journal of Open Research Software 5 no. 1, (2017) . https://app.dimensions.ai/details/publication/pub.1085583166andhttp://openresearchsoftware.metajnl.com/articles/10.5334/jors.151/galley/245/download/. Exported from https://app.dimensions.ai on 2019/05/05.
  * [111] S. G. Johnson, “The HCubature.jl package for multi-dimensional adaptive integration in Julia.” https://github.com/JuliaMath/HCubature.jl, 2017\.
  * [112] P. K. Mogensen and A. N. Riseth, “Optim: A mathematical optimization package for Julia,” Journal of Open Source Software 3 no. 24, (2018) 615.
  * [113] C. Rowley, “Pythoncall.jl: Python and julia in harmony,” 2022.  https://github.com/JuliaPy/PythonCall.jl.
  * [114] S. G. Johnson, “QuadGK.jl: Gauss–Kronrod integration in Julia.” https://github.com/JuliaMath/QuadGK.jl, 2013.
  * [115] J. Verzani, “Roots.jl: Root finding functions for julia.” https://github.com/JuliaMath/Roots.jl, 2020.
  * [116] J. D. Hunter, “Matplotlib: A 2d graphics environment,” Computing in Science & Engineering 9 no. 3, (2007) 90–95.