# A proposed PKU-Muon experiment for muon tomography and dark matter search

Xudong Yu<EMAIL_ADDRESS>State Key Laboratory of Nuclear Physics and
Technology, School of Physics, Peking University, Beijing, 100871, China
Zijian Wang  State Key Laboratory of Nuclear Physics and Technology, School of
Physics, Peking University, Beijing, 100871, China  Cheng-en Liu  State Key
Laboratory of Nuclear Physics and Technology, School of Physics, Peking
University, Beijing, 100871, China  Yiqing Feng  State Key Laboratory of
Nuclear Physics and Technology, School of Physics, Peking University, Beijing,
100871, China  Jinning Li  State Key Laboratory of Nuclear Physics and
Technology, School of Physics, Peking University, Beijing, 100871, China
Xinyue Geng  State Key Laboratory of Nuclear Physics and Technology, School of
Physics, Peking University, Beijing, 100871, China  Yimeng Zhang  State Key
Laboratory of Nuclear Physics and Technology, School of Physics, Peking
University, Beijing, 100871, China  Leyun Gao  State Key Laboratory of Nuclear
Physics and Technology, School of Physics, Peking University, Beijing, 100871,
China  Ruobing Jiang  State Key Laboratory of Nuclear Physics and Technology,
School of Physics, Peking University, Beijing, 100871, China  Youpeng Wu
State Key Laboratory of Nuclear Physics and Technology, School of Physics,
Peking University, Beijing, 100871, China  Chen Zhou<EMAIL_ADDRESS>State Key Laboratory of Nuclear Physics and Technology, School of Physics,
Peking University, Beijing, 100871, China  Qite Li<EMAIL_ADDRESS>State Key
Laboratory of Nuclear Physics and Technology, School of Physics, Peking
University, Beijing, 100871, China  Siguang Wang  State Key Laboratory of
Nuclear Physics and Technology, School of Physics, Peking University, Beijing,
100871, China  Yong Ban  State Key Laboratory of Nuclear Physics and
Technology, School of Physics, Peking University, Beijing, 100871, China
Yajun Mao  State Key Laboratory of Nuclear Physics and Technology, School of
Physics, Peking University, Beijing, 100871, China  Qiang Li
<EMAIL_ADDRESS>State Key Laboratory of Nuclear Physics and Technology,
School of Physics, Peking University, Beijing, 100871, China

###### Abstract

We propose here a set of new methods to directly detect light mass dark matter
through its scattering with abundant atmospheric muons or accelerator beams.
Firstly, we plan to use the free cosmic-ray muons interacting with dark matter
in a volume surrounded by tracking detectors, to trace possible interaction
between dark matter and muons. Secondly, we will interface our device with
domestic or international muon beams. Due to much larger muon intensity and
focused beam, we anticipate the detector can be made further compact and the
resulting sensitivity on dark matter searches will be improved. In line with
above projects, we will develop muon tomography methods and apply them on
atmospheric and environmental sciences, archaeology and civil engineering.
Furthermore, we will measure precisely directional distributions of cosmic-ray
muons, either at mountain or sea level, and the differences may reveal
possible information of dark matter distributed near the earth. In the future,
we may also extend our study to muon on target experiments. Specifically, our
methods can have advantages over ‘exotic’ dark matters which are either muon-
philic or slowed down due to some mechanism.

##  I Introduction and Motivations

Among all the elementary particles, muon has its speciality to connect applied
studies with fundamental researches, yet both to a lesser extent compared with
its cousin, $e.g.$, electron.

Firstly, using cosmic-ray muons to image the internal structure of objects,
especially large-scale objects, has become a booming research area in particle
physics. The earliest application of muography can be traced back to 1955,
when George from Australia measured overburden of tunnel [1]. In the following
decades, muography has been widely used in the field of geology, including the
exploration of rock formations [2, 3], glaciers [4, 5], minerals [6, 7],
oceans [8] and underground carbon dioxide storage [9, 10, 11], $etc$.

In 1969, a team led by Alvarez pioneered muography by measuring the internal
structures of the Giza and Chephren pyramids in Egypt [12]. This marked the
first application of muography in archaeology, showcasing its ability to non-
destructively reveal the hidden features of ancient ruins. Since the 21st
century, muography has played a crucial role in discovering hidden chambers
within structures like the Pyramid of the Sun in Teotihuacan [16] and Khufu’s
Pyramid [13]. In 2021, researchers from Beijing Normal University employed
muography to measure the Mausoleum of Qin Shihuang, successfully
reconstructing the size and location of the underground palace tomb [14]. In
addition, muography has also aided the study of ancient earthquakes [15].

In recent years, monitoring volcanoes has emerged as a prominent research area
in muography. Numerous Japanese volcanoes, including Showa-Shinzan [17], Asama
[18], Sakurajima [19], and Stromboli in Italy [20], has been subject to
muographic studies. Similar investigations have been condected in France [21]
and the United States [22]. In addition to large natural geological
structures, there is a growing interest in combining muography with
meteorology, exemplified by its application in monitoring tropical cyclones
[23].

Due to its non-destructive nature and safety, muon imaging presents a viable
strategy for nuclear safety monitoring and detecting radioactive materials.
Applications include the visualization of reactor interiors [24], detection of
spent nuclear fuel in dry storage barrels [25] and nuclear waste [26].
However, the effectiveness of both muon radiography and muon tomography relies
heavily on precise measurements of energy spectrum and angular distribution of
cosmic-ray muons. Consequently, obtaining accurate data on muon flux is
essential. Various sea-level muon flux models have been proposed and measured
using diverse experimental methods in different geographical locations [27,
28, 29, 30]. Building on this foundation, some research teams have delved into
understanding the impact of factors such as temperature, barometric pressure
[31], and seasonal variations [32] on muon flux.

Cosmic-ray muons are characterized by low density with a mean energy of 3-4
GeV at sea level and a mean flux of approximately 1 cm-2min-1. To produce more
intense and controllable muon beams, several countries have established their
own muon sources.

For instance, the muon source at Paul Scherrer Institute (PSI) in Switzerland,
SµS, uses a continuous-wave high-current proton beam generated by two cascaded
high-current proton cyclotrons. This beam is directed successively onto
targets, providing secondary $\mu$ or $\pi$ beams, dedicated to various
experiments, such as Mu3E, MUSE, µSR [33], MuCOOL, MuSUN, HyperMu, MEG and
Mu3E [34]. TRIUMF’s muon source CMMS has three muon beams, namely M9, M15, and
M20 along with two series-connected muon targets, T1 and T2 [35]. M9A, M15,
and M20 are employed for conducting µSR experiments primarily based on surface
muon beams. M9B is dedicated to µSR experiments using decaying muon beam with
thick samples. Additionally, M11 is mainly used for detector testing. The ISIS
muon source, located at the Rutherford Appleton Laboratory in the UK, is the
world’s first pulsed muon source [36]. Equipped with a thin graphite muon
target, it has two muon beams flanking the target. These beams are dedicated
to conducting µSR experiments, leveraging both surface muon beams and decay
muon beams [37, 38]. J-PARC in Japan, houses two pulsed muon sources: MUSE,
dedicated to µSR experiments, and COMET, specializing in particle physics
experiments [39]. FNAL in the US is actively engaged in groundbreaking
experiments, notably the g-2 experiment aimed at measuring the anomalous
magnetic moment of muons [40, 41], and the Mu2e experiment, which focuses on
investigating rare decays of muons into electrons at its muon campus [42, 43].
Moreover, China and South Korea are also in the process of designing and
planning the construction of their own new muon sources.

Both cosmic-ray and man-made muons stand out as powerful tools for delving
into the mysteries of our world and the unknown. For example, the Standard
Model (SM) of particle physics, while remarkably successful, falls short in
explaining certain experimental observations, such as the neutrino mass and
Dark Matter (DM). To address these gaps, physics beyond the SM is imperative,
introducing hypothetical particles and novel interactions. Among these
hypothesized particles, understanding the nature of DM is currently a focal
point for both cosmology and particle physics. Despite our confidence in its
existence, deduced from cosmological observations, no experiments have
directly detected any form of DM particles.

There are several experiments dedicated to the direct or indirect detection of
light DM particles, whose mass ranges from sub-MeV to GeV scales. The paradigm
of Weakly Interacting Massive Particles (WIMPs) is extensively studied in the
DM sector, with numerous experiments conducted to search for them. However, no
observations of WIMPs have materialized, prompting the need to explore other
theoretically motivated scenarios. Examples include low-mass or muon-philic DM
[44, 45, 46].

Traditional direct searches for DM, involving the identification of nuclear
recoils in deep underground detectors, are challenging due to insufficient
recoil energy in the sub-MeV range. Experimental searches of cold sub-GeV DM
have shifted focus towards exploring the Migdal effect [47] and interactions
with electrons [48, 49]. Proposals for low-mass DM searches using liquid
helium have also been put forward [50, 51]. Recently, experiments also exploit
the fact that a fraction of the cold DM is boosted to relativistic energies
and thus can be efficiently detected in direct detection experiments [52, 53,
54, 55, 56]. More projections on low-threshold DM direct detection in the next
decade can be found in Ref. [44, 57].

Here we extend our previous study [58] and propose a set of new methods to
directly detect light mass DM through its scattering with abundant atmospheric
muons or accelerator beams, to directly probe muon-philic DM in a model-
independent way.

Our methods can have also advantages over ‘exotic’ DMs [59] which can be well
slowed down through scattering with matter in the atmosphere or the Earth
before reaching the detector target. As shown below, methods here rely on high
speed incoming particles and thus depend less on DM velocity.

##  II Our proposal in brief

We are interested in exploiting muon detector such as RPC and GEM for both
Muon tomography and dark matter searches.

As a first step as shown in Fig. 1, we plan to use the free cosmic-ray muons
interacting with dark matter in a volume surrounded by tracking detectors, to
trace possible interaction between dark matter and muons.

Figure 1: Illustration of an experiment to detect muon-philic DM with free
leptons. The resulting kinematic shifts of leptons kicked by DM can be
measurable with tracking detectors surrounding a vacuum region. A veto region
along the chamber can be defined based on the cross-point of the in and out
tracks to suppress backgrounds.

We will then interface our device with domestic or international muon beams as
shown in Fig. 2. Due to much larger muon intensity and focused beam, we
anticipate the detector can be made further compact and the resulting
sensitivity on dark matter searches will be further improved.

Figure 2: Illustration of an experiment to detect muon-philic DM with muon
beams. Due to focused sizes of muon beams, the cylindrical detector can be
made very compact.

In line with above projects, we will develop muon tomography methods on both
detectors and algorithm parts. Depending on possible precision to be achieved,
we may also apply muon tomography on atmospheric and environmental sciences,
on archaeology and civil engineering.

Furthermore, we will measure precisely directional distributions of cosmic-ray
muons, either at mountain or sea level as shown in Fig. 3, and the differences
may reveal possible information of dark matter distributed near the earth.

Figure 3: Illustration of an experiment to detect muon-philic DM through
precisely measuring directional distributions of cosmic-ray muons, either at
mountain or sea level.

##  III Experiment Setup: GEM and RPC

Figure 4: Schematic view of a triple-GEM detector [61].

In order to meet the requirements of muon tomography and dark matter searches
as describe above, the track detectors must contain a list of conditions:
large area ($\sim$ 1 m2), high detection efficiency ($>90\%$), high spatial
resolution ($\sim$ 1 mm), high timing resolution ($\sim$ 1 ns) and low cost.

Gas Electron Multiplier (GEM) [60] detectors are widely used for high-energy
physics experiments, such as the triple-GEM detector [61] installed in the CMS
experiment, due to their excellent performance. The GEM detectors can improve
trigger capabilities and muon measurements in the forward region of CMS
because of their excellent performance: a rate capability above 10 kHZ/cm2, a
time resolution of about 8 ns and a spatial resolution of about 200 $\mu$m.

A GEM foil is a metal-clad polyimide foil which is perforated by a high
density of holes. A triple-GEM detector [61] consists of three GEM foils which
are immersed in a gas mixture and placed between a cathode drift board and an
anode readout board. The drift board routes a total of seven potentials (drift
electrode plus both faces of the three GEM foils) while the readout board is
grounded. A schematic view of a triple-GEM detector is shown in Figure 4.

A charged particle passing the triple-GEM detector produces electrons due to
ionization. The electrons experience intense electric field when moving
through the GEM holes and produce secondary ionization. This electron
avalanche process induces electric signals on the strips or pixels of the
readout board. The signals from the readout board are processed and
transmitted by front-end readout electronics and data acquisition system.

GEM detectors are distinguished by the use of GEM foils as the electron
amplification structure, where the electron avalanche process is independent
of the readout anode. This flexibility allows GEM detectors to adopt various
readout structures. However, since the signal in GEM detectors has minimal
transverse diffusion of the signal (1-3 mm), in order to obtain a good space
resolution, one has to reduce the size of the pixel and employ a large amount
of pixels to cover a required effective area, which will lead to a great
pressure on the detector construction, power consumption, spatial utilization,
etc.

The resistive anode readout method [62] can help to obtain a good spatial
resolution comparable to the pixel readout structure with an enormous
reduction of the electronic channels. By using the thick film resistor
technology, a new type of resistive structure, composed of high resistive
square pad array with low resistive narrow border strips, is developed and
applied to the readout anode of the triple GEM detector. On the basis of that,
We aim to design our exclusive readout for the specific requirements of PKU-
Muon GEM detectors.

Unlike commonly used conductive readout electrodes, where the induced charge
is promptly collected, in resistive structures, there is a continuous
diffusion process for charge induction. Charges will be collected by the
readout nodes of the adjacent cells. Using the information of the collected
charges, the hit position can be reconstructed with suitable algorithm. As
shown in Fig. 5, the black square region represents a high-resistance pad,
while the purple narrow line region represents a low-resistance strip. The
gray small square (located on the backside) serves as the signal readout
electrode. The induced charge will diffuse slowly on high-resistance pad. When
the charge diffuses into the low-resistance region, it will preferentially
spread along the low-resistance strip towards the readout electrode, as the
resistance in this area is significantly lower than that in the adjacent high-
resistance blocks. So we can calculate the hit position according to the
quantity of electric charge that 4 readout electrode collected. Using this
method, at the same time as maintaining high position resolution we can
significantly reduce the required number of electronic channels. We will carry
out further research and practical applications.

Figure 5: Structure diagram of the basic resistive anode cell.

Resistive plate chambers (RPCs) have also been used widely in high energy
physics field as particle detectors for nearly 50 years due to their
advantages such as a simple and robust structure, long-term stability, good
timing resolution, easy-maintenance and low cost. Peking University is an
early researcher in China on muon tomography by successfully building up a
prototype in 2012, i.e., a 2-dimensional glass RPC detector with a large area
based on LC delay-lines readout. The detailed design and assembly of the
prototype glass RPC were described in Ref. [63]. As shown in Fig. 6, two 2.6
mm thick float glass plates were laid in parallel with each area of $30\times
30~{}{\rm cm}^{2}$. Using several 2 mm thick gaskets, the gas gap of 2 mm
between two glass plates was set up. The graphite electrodes were coated on
the glass plates to supply high voltage. The effective area of the electrode
was $20\times 20~{}{\rm cm}^{2}$. Top and bottom were insulated with
$100~{}\mu{\rm m}$ PET layers. And then the outer layer was the printed-
circuit board (PCB), with 2.54 mm readout strips on it. To save the
electronics, we follow the delay-line, charge-division methods, which have
been proved effective in Refs. [64, 65] and help to improve the spatial
resolution as well. One end of each strip was connected to delay-lines with 4
ns delay between adjacent two strips. Signals from one end of delay-lines were
fed into 300 MHz fast preamplifiers. To get X and Y signals in the same time,
we put two prototypes in one aluminum box orthogonally, as shown in Fig. 6(b).
In addition, each layer generates a timing signal channel t, which does not
pass the delay-line. The box was then equipped with HV, signal and gas in-out
connectors. Our all RPCs used here were constructed in the same way, and were
stacked with an equal interval. It is able to decrease the positional
resolution of cosmic-ray ray muons to submillimeter, while maintaining a
detection efficiency above $90\%$.

Figure 6: (a) Layout of the prototype glass RPC. (b) One RPC layer consists of
two same structures, which get X and Y signals respectively. In addition, each
layer has a timing signal channel t not passing the delay-line.

In the near future, we plan to improve the performance of RPC detector. We
will use mixed gas
$(87\%\,\text{Freon}+8.7\%\,\mathrm{C}_{4}\mathrm{H}_{10}+4.3\%\,\mathrm{SF}_{6})$
instead of pure Freon as working gas, so that avalanche mode dominated and in
the meanwhile signal to noise ration kept a high level. Besides, to simplify
the electronic readout, we will adopt a DAQ system based on digitization ASIC,
such as AGET, PETIROC, TOFPET, $etc$.

On top of above mentioned mugraphy system, we will upgrade it to the mode of
dark matter searching, as shown in Fig. 7.

Figure 7: Schematic diagram of RPC testing system

In order to minimize other factors interfering with our dark matter search, we
plan to add a chamber in the middle of RPCs to make the detection environment
at a vacuum. We will test a version with one vacuum cavity first, and later
upgrade the detection system to a version with three vacuum cavities, as shown
in Fig. 8.

Figure 8: Schematic diagram of RPC testing system with vacuum cavity

Additionally, we can combine several cosmic-ray scattering detectors as they
can function simultaneously. This reduces the limitations of individual
detectors and thus improves the accuracy of the detection so that we can
better comprehend the reaction mechanisms of matter. Meanwhile, the use of
multi-detector networking to form large-scale, wide-area cosmic-ray scattering
detection arrays for accumulating data over long periods of time can lower the
lower limit of dark matter detection.

The Active-Target of Time Projection Chamber (AT-TPC) is a type of gas
detector. It uses a working gas medium as the target material. The AT-TPC has
advantages of recording the complete kinematic information and covering a
large solid angle. Peking University has already built a small AT-TPC with a
volume of $14\times 14\times 14~{}{\rm cm}^{3}$. It consists of an electronic
field cage, double-layer GEM membrane for signal amplification and a
2-dimensional strip-readout structure. Its construction is shown as Fig. 9
[66]. In the future, we plan to increase its size to $25\times 25\times
50~{}{\rm cm}^{3}$.

Figure 9: Schematic view of the PKU AT-TPC.

A 2-dimensional strip-readout structure is shown as Fig. 10 [66]. To
characterize the performance of the TPC, it is convenient to define a
coordinate where the Z-axis is along the direction of the mean electron drift
in the field cage and the X-Y axis lies in the readout plane. As shown in Fig.
10(b), a total of 48 parallel X-strips with a width of 700 $\mu$m and 48 rows
of pad plates were placed alternately. Each pad plate was extended to the back
side of the board through a hole and collected on one of the Y-strips. In
order to improve the positional resolution, each Y-strip was divided into two
parts as shown in Fig. 10(a), using a sawtooth cut away in the middle of the
strip.

Figure 10: (a) Schematic view of readout PCB, and (b) amplified details of X-Y
readout structure.

In our vision, an AT-TPC will be placed between RPCs. With 1-dimensional
position readout plus a drift time dimension, we are able to get 2-dimensional
tracks and energy loss. As high-speed muon’s ionizing capacity varies little
before and after the collision, the variation of unit length energy loss
mostly depends on the type of the other scattering particle—a charged particle
will cause a high variation while an uncharged one (like neutrons and dark
matter, but the muon’s energy is too low to knock out a neutron) will not.
This method provides a further confirmation on whether a muon scattering event
infers the dark matter’s existence.

Therefore, we will apply AT-TPC to get the track of particles for better
precision, as shown in Fig. 11. Also, AT-TPC can work on some muon beamlines
of relatively low energy, such as MuSIC in the RCNP [67]. For the reason that
some muons in relatively low energy ranging from a few MeV to tens of MeV
cannot traverse RPC or will be scattered strongly, AT-TPC can be a suitable
detector in this situation.

Figure 11: Schematic diagram of RPC testing system with AT-TPC.

##  IV Simulations Framework

The simulation of the GEM-based detector is conducted using Geant4 [70]. The
schematic design of the triple-GEM detector can be found in Ref. [71]. The
drift cathode and readout PCB consist of 0.1 mm thick copper positioned on
both sides of the GEM. The readout strip pitch is 210 $\mu$m, and the strip
width is 150 $\mu$m. The GEM foil is composed of a 50 $\mu$m thick Kapton foil
with 5 $\mu$m thick copper cladding on both sides. In the original GEM design,
the foil is etched using photolithographic processes to form tiny bi-conical
holes arranged in a regular hexagonal pattern. In our work, we approximate the
avalanche amplification process that occurs when an electron passes through a
hole with very high electric field strength. We achieve this by simulating the
material effects of the tiny holes, reducing the density of the copper plate
equivalently. The position resolution, including electron avalanche
amplification, charge collection, readout strip pitch, $etc.$, is considered
by a smearing according to the nominal position resolution. The thicknesses of
the drift region, the three transfer regions, and the induction region are
4.8-2-2-2 mm. The upper and lower tracking systems consist of two stacked
plates, between which is the $1\times 1\times 1~{}{\rm m}^{3}$ sensitive
detection region in physics cases I and II below. The box is covered with a
100 $\mu$m thick Kapton shell. When tuned for dark matter searching mode, the
box is filled with air or vacuum. In Tomography mode, the sample to be tested
is placed inside.

###  IV.1 Dark Matter and muon scattering

The interaction of low-mass muon-philic DM interaction with free muons is
illustrated in Figs. 1, 2, and 3. In the presence of surrounding DM, these
atmospheric or accelerator muons might experience spatial shifts, providing a
distinct signal for detection. According to the standard halo model [72, 73],
the velocity distribution of DM particles in the galactic rest frame follows a
Maxwell-Boltzmann distribution with a cutoff at an escape velocity
$v_{\mathrm{esc}}$, which is

$F_{g}\left(\mathbf{v}_{\mathbf{g}}\right)=\begin{cases}k\left(\frac{1}{\sqrt{\pi}v_{0}}\right)^{3}\exp\left(-\left(v_{g}/v_{0}\right)^{2}\right)&v_{g}<v_{\mathrm{esc}}\\\
0&v_{g}\geq v_{\mathrm{esc}}\end{cases},$ (1)

where $v_{0}=\sqrt{2K_{B}T/m}$ is the most probable or modal speed of the DM
particles, and the normalization constant $k$ (1 if $v_{\mathrm{esc}}=\infty$)
is

$k=\mathrm{erf}(v_{\mathrm{esc}}/v_{0})-\frac{2}{\sqrt{\pi}}\frac{v_{\mathrm{esc}}}{v_{0}}\exp{-(v_{\mathrm{esc}}/v_{0})^{2}}.$
(2)

As the Earth moves relative to the galactic rest frame with a velocity
$v_{E}$, the velocity distribution of DM particles detected in an Earth-based
detector is not isotropic. Defining the polar angle $\phi$ with its origin in
the direction of the Earth’s motion, and the azimuth angle $\theta$ as the
angle between the observed speed $v$ and the speed of the Earth $v_{E}$, the
speed distribution is obtained by integrating over $\phi$ and $\theta$. The
distribution is given by:

$f(v)=\int_{0}^{2\pi}d\phi\int_{0}^{\pi}d\theta
v^{2}\sin(\theta)F_{g}\left(v_{g}\right).$ (3)

DM scattering experiments almost universally assume the standard halo model
with parameters $v_{0}=220~{}\mathrm{km/s}$, $v_{E}=232~{}\mathrm{km/s}$, and
$v_{\mathrm{esc}}=544~{}\mathrm{km/s}$ to ensure their results can be fairly
compared with each other. We also consider a simpler case by neglecting the
velocity distribution of DM, taking it as a constant value
$v=220~{}\mathrm{km/s}$.

In our work, we consider a model-independent elastic scattering between muons
and DMs, following Newtonian mechanics. Conservation of energy and momentum
can be used to show the muon recoil energy $E_{\mathrm{recoil}}$ obeys [72,
73]

$\displaystyle E_{\mathrm{recoil}}(v)$ $\displaystyle=\frac{{\rm
M}_{\mathrm{DM}}v^{2}}{2}\frac{1-\cos\theta}{2}\frac{4{\rm
M}_{\mathrm{DM}}{\rm M}_{\mu}}{({\rm M}_{\mathrm{DM}}+{\rm M}_{\mu})^{2}}$ (4)
$\displaystyle=\frac{{\rm M}^{\prime}_{\mathrm{DM},\mu}v^{2}}{{\rm
M}_{\mu}}(1-\cos\theta),$

where $\theta$ is the scattering angle in the center-of-mass frame and ${\rm
M}^{\prime}_{\mathrm{DM},\mu}$ is the DM-muon reduced mass. The maximum recoil
energy is achieved at a head-on collision ($\theta=\pi$):

$E_{\mathrm{recoil}}^{\mathrm{max}}=\frac{2{\rm
M}^{\prime}_{\mathrm{DM},\mu}v^{2}}{{\rm M}_{\mu}}.$ (5)

The incident cosmic-ray muons are simulated by CRY package [74]. Since the
mean energy of cosmic-ray muons at sea level is about 3-4 GeV, making them
relativistic, DM is quasi-frozen from the perspective of fast incoming muons.
This differs from other DM search experiments, such as XENON1T and PandaX,
where nuclei are quasi-static target and DMs are assumed to be incoming
particles. Notice for high speed muons, it is appropriate to treat DM as
frozen in the detector volume, and the estimated muon-DM scattering rate per
second could be, for the case of cosmic-ray muons:

$dN/dt=\rho V/{\rm M}_{\rm DM}\times\sigma_{\mu,\mathrm{DM}}\times F_{\mu}.$
(6)

Here, $F_{\mu}$ is the muon flux $\sim$ 1/60 s-1cm-2 at sea level. The local
density of DM is on the order of $\rho\sim$ 0.3 GeV/$\rm{cm}^{3}$, and the
number of DM in a box with volume $V$ is $N_{\mathrm{DM}}=\frac{\rho V}{{\rm
M}_{\mathrm{DM}}}$. For the DM mass ${\rm M}_{\mathrm{DM}}\sim 1$ (0.1, 10)
GeV, and the detector box volume $V\sim 1\,{\rm m}^{3}$, the number of DMs is
about $3\times 10^{5(6,4)}$. The sensitivity of DM-muon scattering cross
section for 1 year run will be around

$\sigma_{\mu,\mathrm{DM}}\sim 10^{-12(-13,-11)}~{}{\rm cm}^{2}.$ (7)

On the other hand, in the case of muon beams, we have:

$dN/dt=N_{\mu}\times\sigma_{\mu,{\mathrm{DM}}}\times L\times\rho/{\rm{\rm
M}_{\rm DM}}.$ (8)

Taking the example of ${\rm M}_{\mathrm{DM}}=0.03$ GeV, $L=1$ m, and
$N_{\mu}\sim 10^{6}$/s ($e.g.$, China CSNS Melody design), we get the
following estimation:

$N=10^{13}\times\sigma_{\mu,\mathrm{DM}}\times 100/{\rm cm}^{2},$ (9)

Thus the sensitivity on DM-muon scattering cross section for 1 year run will
be around

$\sigma_{\mu,\mathrm{DM}}\sim 10^{-15}~{}{\rm cm}^{2}.$ (10)

Notice, however, detector acceptance and efficiency needs to be applied, which
can change above estimated results and will be further explained in the
sections below.

##  V Physics Case I: Precise measurements of cosmic-ray muons and Muon
Tomography

The first differential and integral spectra of muons moving in the near-
vertical direction can be measured at Beijing using the devices as mentioned
above. As the first step, we put 4 RPCs in our muon tomography systems, as
illustrated in Fig. 12. When cosmic-ray muons traverse materials, they undergo
multiple deflections. The test results, depicted in Fig. 13, clearly show the
outline of the lead block in the iron shell, indicating differences in muon
scattering angles as they pass through various materials.

Figure 12: Prototype of muon imaging system based on delay line readout and
glass RPC developed by Peking University. Figure 13: Imaging results of a
$6\times 6\times 6~{}{\rm cm}^{3}$ square lead block wrapped in a $12\times
12\times 12~{}{\rm cm}^{3}$ iron shell and the letters PKU composed of
$3\times 3\times 3~{}{\rm cm}^{3}$ iron blocks by muon imaging prototype of
Peking University. The color represents the average scattering angle of the
corresponding region.

The idea of muon tomography is based on the fact that cosmic-ray muons will
undergo Coulomb scattering when they pass through materials. Therefore, by
measuring the muon scattering angular distribution, the distribution
information of substances or materials in the detection area can be
reconstructed. In 2004, Schultz et al. first conducted such experiments and
proposed Point of Closest Approach (PoCA) algorithm [75], in which multiple
consecutive scatterings of muons are approximated as one scattering, and the
scattering point (called “PoCA” point) is located at the midpoint of the
common perpendicular line between the incident and outgoing straight lines.
Subsequently, in 2006-2007, Schultz et al. proposed Maximum Likelihood
Scattering and Displacement (MLSD) algorithm [76, 77], a new iterative imaging
algorithm based on the maximum likelihood principle in statistics. This
algorithm is proposed to take into account the muons scattering angle and then
introduce the influence of the lateral displacement during the muons
penetration process into the material. However, these algorithms have certain
technical limitations. To date, new muon tomography imaging algorithms
continue to be proposed, for example, the Ratio algorithm based on the PoCA
algorithm [78].

##  VI Physics Case II: Dark matter searches in a box

For muon-philic DM, one can start with a relatively straightforward searches
using those free cosmic-ray muons, as depicted in Fig. 1. A quick estimation
of the projected sensitive on DM muon scattering cross section is shown in Eq.
7, which needs further refinement to incorporate detailed simulation effects,
as explained blow.

We simulate 5 billion incident muon events to study the backgrounds in one-
year exposure, and also simulate 10 million signal events under different DM
assumptions to estimate the signal efficiency. Due to the low cross section of
DM-muon interaction, we assume that there is only one interaction in the box.
We require a muon hit signal in each GEM detector. The MC truth hit position,
smeared by the GEM spatial resolution is taken as the reconstructed hit
position. The two hits in the upper (lower) GEM detector determine the
direction of incoming (outgoing) muons. The angle between the directions of
incoming and outgoing muon, $\theta$, is the discriminant variable. The
$\cos\theta$ distributions in signal and background samples are shown in Fig.
14. The signal distributions are scaled to the event numbers in the
background. We list some interesting results as below.

  * •

The $\cos\theta$ distribution in air has no obvious difference between that in
a vacuum. Therefore, the scattering between muons and air has no significant
effect under the current GEM spatial resolution. Considering cost and
technical difficulty, vacuuming of the boxes is not necessary in Phase I of
the project.

  * •

The $\cos\theta$ distributions in Maxwell-Bolzmann velocity distribution and a
constant velocity distribution are similar. Therefore, our signal distribution
and detection is not sensitive to the DM velocity model.

  * •

As the DM mass increases, a larger fraction occupies the region of large
scattering angles, resulting a more pronounced discrepancy between the signal
and background.

  * •

For the signal event with ${\rm M}_{\mathrm{DM}}<100$ MeV, an apparent
truncation is observed, attributed to kinematics. This truncation occurs only
when the DM mass is lower than the muon mass.

To further suppress backgrounds, we require that the reconstructed PoCA point
should be within the box, which is reasonable and supported by Punzi Figure-
of-Merit, $\mathrm{FoM}=\varepsilon/(3/2+\sqrt{B})$ [79]. Here, $\varepsilon$
is the detection efficiency, and, and $B$ is the event number of background
events. The FoM value is improved by 18%. The resulting detection efficiency
under different DM mass assumptions and the background event numbers are
listed in Table. 1.

Figure 14: The $\cos\theta$ distributions in signal and background samples. The black histogram is the background, while the other colored histograms denote signal under different DM mass assumptions, scaled to the same event number. In the large scattering region, $i.e.$, $\cos\theta<0$, signal distributions with a large DM mass and background distributions manifest a noticeable difference. Background |  Event Number ($\times 10^{9}$)
---|---
Air | 1.15
Vacuum | 1.14
DM mass (GeV) | Constant (%) | Maxwell-Bolzmann (%)
0.005 | $27.10\pm 0.01$ | $27.11\pm 0.01$
0.05 | $29.56\pm 0.01$ | $29.55\pm 0.01$
0.1 | $27.66\pm 0.01$ | $27.64\pm 0.01$
0.2 | $25.01\pm 0.01$ | $24.99\pm 0.01$
0.5 | $21.47\pm 0.01$ | $21.46\pm 0.01$
1 | $18.67\pm 0.01$ | $18.66\pm 0.01$
10 | $11.10\pm 0.01$ | $11.10\pm 0.01$
100 | $8.44\pm 0.01$ | $8.43\pm 0.01$
Table 1: Background event numbers corresponding to the integrated luminosity
of one-year exposure with the box filled with air and vacuum, along with the
signal detection efficiency under different assumptions of DM velocity
distribution and mass.

We search for excesses in the $\cos\theta$ distribution at the lower range by
performing binned maximum likelihood fits using Higgscombine [80]. As observed
yields are unavailable for this R&D study, “Asimov” data is used. We set the
upper limit (UL) using the CLs technique [81, 82]. Assuming the measurements
performed are limited by the availability of data statistics, have negligible
experimental systematic, and have excellent MC statistics, the background
uncertainty is taken as the Poisson counting uncertainty for the expected
background yield in each bin. The upper limit of the DM-muon interaction cross
section at 95% confidence level (CL) corresponding to an integrated luminosity
of one year are shown in Fig. 15. Taking into account the detection
efficiency, the upper limits on cross section for sub-GeV DM and muon
interaction are set to $10^{-7}\sim 10^{-9}\,\mathrm{cm}^{2}$.

Figure 15: Expected 95% CL upper limits on the DM-muon interaction cross
section versus DM mass. The green and yellow bands denote $1\sigma$ and
$2\sigma$ region, respectively.

##  VII Physics Case III: Dark matter searches using muon beam

As estimated with Eq. 9 and shown in Fig. 2, owing to a much larger muon
intensity and focused beam, we anticipate that the detector can be made more
compact, leading to a further improvement in sensitivity for DM searches.

To validate the correctness of estimation, we perform a MC simulation study on
DM searches using a muon beam via Geant4 [70]. As shown in Fig. 16, to suit
detection environment of the beam experiment, we adopt cylindrical GEM (CGEM)
detector structure. This innovative technique has been used in the upgrade of
the BESIII inner tracker system [83]. According to the MELODY design [84, 85],
the diameter of the beam spot, $\phi$, ranges from $10\,\mathrm{mm}$ to
$30\,\mathrm{mm}$. In our study, we choose $\phi=10\,\mathrm{mm}$. To reduce
beam background related to the size of the beam spot, the inner diameter of
our CGEM detector is designed to be $50\,\mathrm{mm}$, which is 5 times the
beam spot. The two layers of GEM detectors are stacked together, and the
structure of each GEM is the same as described above. The length of the GEM
detector is designed to be 1 m. For muon beam, we consider four different
cases with low muon momentum $P_{\mu}=10\,\mathrm{MeV}$, $100\,\mathrm{MeV}$
and high muon kinetic energy $E^{\mu}_{\mathrm{kin}}=1\,\mathrm{GeV}$ and
$10\,\mathrm{GeV}$. The momentum direction of the muon is along the $z$ axis,
which is the symmetry axis of the detector, and the profile of beam in the
$xy$ plane follows a Gaussian distribution.

Figure 16: Simulating 1 GeV muon beam hit lead plate passing through GEM
detector. Orange surfaces are drift cathodes. The blue surfaces are GEM foils.
The green surfaces are PCBs. The yellow lines are muons tracks. The red curves
are electron tracks. The green lines are photon.

For the background case, we assume that the incoming muons scatter with air in
the beam pipe, resulting in changes in their momenta. If the scattering angle
is large enough, muons may hit the detector. The two hit positions in the two
GEM detectors can reconstruct the flying direction of the outgoing muon. For
signal case, we assume that the incoming muons experience one elastic
scattering with DM particles in the detection region. The signal model is
described above, and here we only consider the case of a DM Maxwell-Bolzmann
velocity distribution. Figure 17 shows the $\cos\theta$ distributions of
signal and background samples. Notably, since $P_{\mu}=10\,\mathrm{MeV}$ is
too low, muon will lose all its energy scattering with air, resulting in no
events surviving in the background simulation. This supports the idea that the
beam pipe should be in a vacuum. However, in other cases where muon momentum
is large enough, muons will hit the GEM detector. The resulting detection
efficiency under different DM mass assumption and different muon beam energies
are listed in Table 2.

We search for excesses in the $\cos\theta$ distributions by performing binned
maximum likelihood fits using Higgscombine [80]. Since observed yields are
unavailable, “Asimov” data is used. We set the upper limit (UL) using the CLs
technique [81, 82]. Assuming the measurements performed are limited by the
availability of data statistics, have negligible experimental systematic, and
have excellent MC statistics, the background uncertainty is taken as the
Poisson counting uncertainty for expected background yield in each bin. The
upper limit of the DM-muon interaction cross section at 95% confidence level
(CL) corresponding to an integrated luminosity of one year are shown in Fig.
18. Here, we consider a muon intensity is $I=10^{6}\,\mathrm{s}^{-1}$,
consistent with the CSNS Melody design. Taking account detection efficiency,
the upper limits on the cross section for sub-GeV DM and muon interaction are
set to $10^{-9}\sim 10^{-12}\,\mathrm{cm}^{2}$ when muon momentum is 100 MeV,
$10^{-11}\sim 10^{-13}\,\mathrm{cm}^{2}$ when muon energy is 1 GeV, and
$10^{-12}\sim 10^{-13}\,\mathrm{cm}^{2}$ when muon energy is 10 GeV.

Figure 17: The $\cos\theta$ distributions with muon energies of 100 MeV, 1 GeV and 10 GeV in signal and background samples. The black histogram is the background, while the colored histograms denote the signal under different DM mass assumptions scaled to the same event number. In the large scattering region, $i.e.$, $\cos\theta<0$, the signal with a large DM mass and background distributions manifest a remarkable difference. ${\rm M}_{\rm DM}$ \$E^{\mu}_{\mathrm{kin}}$ | 100 MeV (%) | 1 GeV (%) | 10 GeV (%)
---|---|---|---
0.05 GeV | $84.29\pm 0.04$ | $74.85\pm 0.04$ | $45.93\pm 0.05$
0.1 GeV | $91.74\pm 0.03$ | $83.07\pm 0.04$ | $58.17\pm 0.05$
0.2 GeV | $94.35\pm 0.02$ | $88.16\pm 0.03$ | $68.37\pm 0.05$
0.5 GeV | $95.17\pm 0.02$ | $92.16\pm 0.03$ | $78.91\pm 0.04$
1 GeV | $95.34\pm 0.02$ | $93.88\pm 0.02$ | $84.68\pm 0.04$
10 GeV | $95.35\pm 0.02$ | $95.36\pm 0.02$ | $94.06\pm 0.02$
100 GeV | $95.43\pm 0.02$ | $95.37\pm 0.02$ | $95.37\pm 0.02$
Table 2: Signal detection efficiency under different assumptions of DM mass
and muon beam energies.

Figure 18: Expected 95% CL upper limits on the DM-muon interaction cross
section versus DM mass with muon beam energy 100 MeV, 1 GeV and 10 GeV. The
green and yellow bands denote $1\sigma$ and $2\sigma$ region, respectively.

##  VIII Physics Case IV: Dark matter searches between mountain and sea level

Directional distributions of cosmic-ray muons can be precisely measured,
either at mountain or sea level, as shown in Fig. 3, and their differences may
reveal possible information of DM distributed near the earth.

In the Geant4 framework, the dominant physics process of muon-air interaction
is ionization, in which the energy loss of muon is described by the Bethe-
Bloch formula [86]. When traveling in the atmosphere, muons scatter with
electrons of air molecule, and their momentum directions change accordingly.
Due to limitations in computing resources, we simulate $1.5\times 10^{9}$
events, which is approximately $30\%$ of the statistics for one year of data
from a 1 m2 detector. We place two layers of Gem detectors, both at the sea
level and in a mountain with a latitude of 100 m, to measure the flying
direction of cosmic-ray muons. The $\cos\theta$ distributions of muons at the
mountain and sea levels are shown in Fig. 19.

Figure 19: The $\cos\theta$ distributions at the mountain and sea levels in
the background sample.

Assuming DM clouds surround the earth, there is a probability that cosmic-ray
muons scatter with DM particles when traveling through the air. To simulate
signal MC samples, we extend the native Geant4 physics list by introducing the
muon-DM elastic scattering process. The cross section of muon-DM elastic
scattering process can be set manually by mean free path $L$. The relation
between muon-DM cross section and mean free path is given by
$L=1/(n_{\mathrm{DM}}\times\sigma_{\mu,\mathrm{DM}})$, where
$n_{\mathrm{DM}}=\rho_{\mathrm{DM}}/{\rm M}_{\mathrm{DM}}$ is the number
density of DM, and $\sigma_{\mu,\mathrm{DM}}$ is the muon-DM cross section.
For the case of ${\rm M}_{\mathrm{DM}}=500\,\mathrm{MeV}$, we simulate
different signal MC samples with different $L$ assumption. Each signal MC
sample include $10^{8}$ events. For $L=10^{1}\sim 10^{7}$ km, it is expected
to have $10^{6}\sim 1$ events with a large scattering angle. The $\cos\theta$
distributions at the sea level are shown in Fig. 20. We can observe a
significant discrepancy between the background distribution and signal
distributions with $L=10$ km and 100 km. For $L\geq 100$ km, no significant
differences are found. According to the $\chi^{2}$ test method, we calculate
the $\chi^{2}/\mathrm{ndf}$ between the signal histogram and background
histogram. The resulting $p$-values indicate that $L=10$ km
($\sigma_{\mu,\mathrm{DM}}=1.67\times 10^{-6}$ cm2) and 100 km ($1.67\times
10^{-7}$ cm2) signal assumptions can be rejected by more than $20\sigma$,
while the significances for $L\geq 100$ km signal assumptions are lower then
$5\sigma$. It indicates that the rejection ability for $M=500$ MeV DM is about
$\sigma_{\mu,\mathrm{DM}}=10^{-7}\sim 10^{-8}$ cm2 for 1 m2 detectors placed
at a mountain with a latitude of 100 m and the sea level after one-year data
taking.

Figure 20: The $\cos\theta$ distributions at the sea level in different signal
samples. The DM mass is 500 MeV, and the mean free path ranges from $L=10^{1}$
km to $L=10^{7}$ km. The histograms have been scaled to one year statistics.

##  IX Next: Muon on target

Lepton beams may play an important role in the detection of light DM
particles. There are some proposals using electron beams for direct detection
of sub-GeV DM particles [87, 88, 89, 90, 91]. These approaches utilize
electron beams to probe invisibly-decaying particles which couple to
electrons. Some theoretical models like Muonphilic DM and Lepton portal DM
[92, 46] show that DM states only couple to the charged leptons, and
predominantly interact with muons with a new type of interaction called muonic
force. The quest for explaining the pronounced muon magnetic moment anomaly
from the Fermilab Muon g-2 collaboration with 4.2 $\sigma$ discrepancy [93]
drives people to search such DM preferentially interacting with muons through
proton beam-dump experiments utilizing muon beams [94]. The muonphilic DM
scenario may reconcile this anomaly by introducing a force carrier particle,
muonphilic bosons [95].

Accelerator muon beams have been proposed as a potential tool for detecting DM
particles by NA64-$\mu$ experiment [96] (with first results public recently in
Ref. [97]) as well as ${\rm M}^{3}$ [98] and FNAL-$\mu$ experiment at Fermilab
[99], and also for the future muon collider [100]. The idea behind is to look
for the muon deflection caused by scattering with DM as well as the the energy
loss pattern of muons.

Accordingly, in extension of above mentioned programs in previous sections by
exploiting muon beams, we are aiming at further launching a muon missing
momentum experiment in the future, which needs further optimization R&D
studies, especially when interfacing with domestic muon beams with lower
energy. More details can be found in a forthcoming study by the authors.

##  X Summary and outlook

We propose here a set of new methods to directly detect light mass DM through
its scattering with abundant atmospheric muons or accelerator beams. Firstly,
we plan to use the free cosmic-ray muons interacting with dark matter in a
volume surrounded by tracking detectors, to trace possible interaction between
dark matter and muons. Secondly, we will interface our device with domestic or
international muon beams. Due to much larger muon intensity and focused beam,
we anticipate the detector can be made further compact and the resulting
sensitivity on dark matter searches will be improved. In line with above
projects, we will develop muon tomography methods and apply them on
atmospheric and environmental sciences, archaeology and civil engineering.
Furthermore, we will measure precisely directional distributions of cosmic-ray
muons, either at mountain or sea level, and the differences may reveal
possible information of dark matter distributed near the earth. In the future,
we may also extend our study to muon on target experiments.

Above various methods using either cosmic muons or accelerator muon beams,
will provide a first direct probe on muon and dark matter interactions, with
sensitivities spanning over different orders of magnitudes. Specifically, our
methods can have advantages over ‘exotic’ dark matters which are either muon-
philic or slowed down due to some mechanism. These searches complement and
enrich current dark matter programs, and may stimulate further efforts.

###### Acknowledgements.

This work is supported in part by the National Natural Science Foundation of
China under Grants No. 12150005, No. 12325504, No. 12075004, and No.
12061141002. Xudong Yu, Chen Zhou, Qite Li and Qiang Li contributed equally to
this work.

## References

  * [1] EP George, Commonwealth Engineer 455,1955
  * [2] L. Malmqvist, G. Jönsson, K. Kristiansson, L. Jacobsson, Geophysics 44, no. 9, 1549-1569(1979).
  * [3] F. Hivert, I. Lázaro Roche, J. B. Decitre,J. Brunner, J. Busto, S. Gaffet, Near Surface Geophysics,15, 2017.
  * [4] R. Nishiyama, A. Ariga, T. Ariga, S. Käser, F. Schlunegger, Geophysical Research Letters 44, no. 12, 2017.
  * [5] R. Nishiyama, A. Ariga, T. Ariga, A. Lechmann, D. Mair, C. Pistillo, P. Scampoli, PG. Valla, M. Vladymyrov, A. Ereditato, et al. Scientific Reports 9, no. 1, 6970,(2019).
  * [6] A. Bonneville, R. Kouzes, J. Yamaoka, A. Lintereur, J. Flygare, G. S. Varner, I. Mostafanezhad, E. Guardincerri, C. Rowe, R. Mellors, Philosophical Transactions of the Royal Society A 377, no. 2137, 20180060(2019).
  * [7] Z. X. Zhang, T. Enqvist, M. Holma, P. Kuusiniemi, Rock Mechanics and Rock Engineering 1-15(2020).
  * [8] H. K. M. Tanaka, M. Aichi, C. Bozza, R. Coniglione, J. Gluyas, N. Hayashi, M. Holma, O. Kamoshida, Y. Kato, T. Kin et al. Scientific Reports 11, no. 1, 19485, (2021).
  * [9] V. A. Kudryavtsev, N. J. C. Spooner, J. Gluyas, C. Fung, M. Coleman, International Journal of Greenhouse Gas Control 11, 21-24(2012).
  * [10] J. Klinger, S. J. Clark, M. Coleman, J. G. Gluyas, V. A. Kudryavtsev, D. L. Lincoln, S. Pal, S. M. Paling, N. J. C. Spooner, S. Telfer et al. International Journal of Greenhouse Gas Control 42, 644-654(2015).
  * [11] J. Gluyas, L. Thompson, D. Allen, C. Benton, P. Chadwick, S. Clark, J. Klinger, V. Kudryavtsev, D. Lincoln, B. Maunder et al. Philosophical Transactions of the Royal Society A 377, no. 2137, 20180059(2019).
  * [12] L. W. Alvarez, J. A. Anderson, F. E. Bedwei, J. Burkhard, A. Fakhry, A. Girgis, A. Goneid, F. Hassan, D. Iverson, G. Lynch et al. Science 167, no. 3919, 832-839(1970).
  * [13] K. Morishima, M. Kuno, A. Nishio, N. Kitagawa, Y. Manabe, M. Moto, F. Takasaki, H. Fujii, K. Satoh, H. Kodama et al. Nature 552, no. 7685, 386-390(2017).
  * [14] S. Ning, L. Yuanyuan, W. Li, C. Jianping, Acta Physica Sinica,2020.
  * [15] Tanaka, H. K. M., Sumiya, K., and Oláh, L.: Muography as a new tool to study the historic earthquakes recorded in ancient burial mounds, Geosci. Instrum. Method. Data Syst., 9, 357–364, https://doi.org/10.5194/gi-9-357-2020, 2020.
  * [16] S. Aguilar, R. Alfaro, E. Belmont, 2013\.
  * [17] H. K. M. Tanaka, T. Nakano, S. Takahashi, J. Yoshida, H. Ohshima, T. Maekawa, H. Watanabe, K. Niwa, Geophysical Research Letters 34, no. 22(2007).
  * [18] H. K. M. Tanaka, T. Uchida, M. Tanaka, M. Takeo, J. Oikawa, T. Ohminato, Y. Aoki, E. Koyama, H. Tsuji, Geophysical Research Letters 36, no. 17(2009).
  * [19] L. Oláh, H. K. M. Tanaka, T. Ohminato, D. Varga, Scientific Reports 8, no. 1, 3207(2018).
  * [20] V. Tioukov, A. Alexandrov, C. Bozza, L. Consiglio, N. D’Ambrosio, G. De Lellis, C. De Sio, F. Giudicepietro, G. Macedonio, S. Miyamoto et al. Scientific Reports 9, no. 1, 6695(2019).
  * [21] F. Ambrosino, A. Anastasio, A. Bross, S. Béné, P. Boivin, L. Bonechi, C. Cârloganu, R. Ciaranfi, L. Cimmino, Ch. Combaret et al. Journal of Geophysical Research: Solid Earth 120, no. 11, 7290-7307(2015).
  * [22] I. D. Guerrero, D. F. Cabrera, J. C. Paz, J. D. Estrada, C. A. Villota, E. A. Velasco, F. E. Fajardo, O. Rodriguez, J. Rodriguez, D. Arturo et al. Journal of Physics: Conference Series 1247, no. 1, 012020(2019).
  * [23] H. K. M. Tanaka, J. Gluyas, M. Holma, J. Joutsenvaara, P. Kuusiniemi, G. Leone, D. Lo Presti, J. Matsushima, L. Oláh, S. Steigerwald et al. Scientific Reports 12, no. 1, 16710(2022).
  * [24] K. Takamatsu, H. Takegami, C. Ito, K. Suzuki, H. Ohnuma, R. Hino, T. Okumura, Annals of Nuclear Energy 78, 166-175(2015).
  * [25] D. Poulson, J. M. Durham, E. Guardincerri, C. L. Morris, J. D. Bacon, K. Plaud-Ramos, D. Morley, A. A. Hecht, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 842, 48-53(2017).
  * [26] D. Mahon, A. Clarkson, S. Gardner, D. Ireland, R. Jebali, R. Kaiser, M. Ryan, C. Shearer, G. Yang, Philosophical Transactions of the Royal Society A 377, no. 2137, 20180048(2019).
  * [27] J.-L. Autran, D. Munteanu, T. S. Saad, S. Moindjie, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 903, 77-84(2018).
  * [28] N. Su, Y. Liu, L. Wang, B. Wu, J. Cheng, Frontiers in Energy Research 9, 750159(2021).
  * [29] M. Abbrescia, C. Avanzini, L. Baldini, R. Baldini Ferroli, G. Batignani, M. Battaglieri, S. Boi, E. Bossini, F. Carnesecchi, D. Cavazza, and others, The European Physical Journal C 83, no. 4, 293(2023).
  * [30] Elio Conte, Nicoletta Sala, Marco Arcani, Symmetry 15 no. 4, 659(2023).
  * [31] I. Alekseev, V. Belov, V. Brudanin, A. Bystryakov, M. Danilov, V. Egorov, D. Filosofov, M. Fomina, S. Kazartsev, A. Kobyakin, et al. The European Physical Journal C 82, no. 6, 515(2022).
  * [32] F. P. An, A. B. Balantekin, H. R. Band, M. Bishai, S. Blyth, D. Cao, G. F. Cao, J. Cao, Y. L. Chan, J. F. Chang, et al. Journal of Cosmology and Astroparticle Physic 2018, no. 01, 001(2018).
  * [33] E. Morenzoni, H. Glückler, T. Prokscha, HP. Weber, EM. Forgan, TJ. Jackson, H. Luetkens, Ch. Niedermayer, M. Pleines, M. Birke, et al. Physica B: Condensed Matter,289, 653-657 (2000).
  * [34] A. Bravar, The Mu3e Experiment at PSI,2016.
  * [35] TRIUMF/CMMS website, http://cmms.triumf.ca/
  * [36] ISIS muon source website, https：//www.isis.stfc.ac.uk/Pages/Muons.aspx
  * [37] A. D. Hillier, D. J. Adams, P. J. Baker, A. Bekasovs, F. C. Coomer, S. P. Cottrell, S. D. Higgins, S. J. S. Jago, K. G. Jones, J. S. Lord, et al. Journal of Physics: Conference Series,551,no.1, 012067 (2014).
  * [38] T. Matsuzaki, K. Nagamine, M. Tanase, M. Kato, K. Kurosawa, H. Sugai, K. Ishida, S. N. Nakamura, I. Watanabe, G. H. Eaton, Nuclear Instruments and Methods in Physics Research A, 480,no.2-3, 814-827 (2002).
  * [39] J-PARC/MUSE website, http：//www.j-parc.jp/MatLife/en/index.html
  * [40] Fermilab website, https://www.fnal.gov/
  * [41] A. Chapelain, Y. Foka, N. Brambilla, V. Kovalenko, Epj Web of Conferences,137,08001(2017).
  * [42] L. Bartoszek, E. Barnes, J. P. Miller, J. Mott, J. L. Popp, 2015\.
  * [43] R. P. Litchfield, 2014\.
  * [44] R. Essig, G. K. Giovanetti, N. Kurinsky, D. McKinsey, K. Ramanathan, K. Stifter and T. T. Yu, [arXiv:2203.08297 [hep-ph]].
  * [45] P. Harris, P. Schuster and J. Zupan, [arXiv:2207.08990 [hep-ph]].
  * [46] Y. Bai and J. Berger, JHEP 08, 153 (2014) doi:10.1007/JHEP08(2014)153 [arXiv:1402.6696 [hep-ph]].
  * [47] M. Ibe, W. Nakano, Y. Shoji, and K. Suzuki, JHEP 03, 194 (2018), arXiv:1707.07258 [hep-ph]
  * [48] Q. Arnaud et al. (EDELWEISS), Phys. Rev. Lett. 125, 141301 (2020), arXiv:2003.01046 [astro-ph.GA]
  * [49] L. Barak et al. (SENSEI), Phys. Rev. Lett. 125, 171802 (2020), arXiv:2004.11378 [astro-ph.CO].
  * [50] J. Liao, Y. Gao, Z. Liang, Z. Peng, Z. Ouyang, L. Zhang, L. Zhang and J. Zhou, [arXiv:2103.02161 [astro-ph.IM]].
  * [51] J. Liao, Y. Gao, Z. Jiang, Z. Liang, Z. OuYang, Z. Peng, F. Zhang, L. Zhang and J. Zhou, doi:10.1140/epjp/s13360-023-03747-2 [arXiv:2209.02320 [astro-ph.IM]].
  * [52] K. Abe et al. [Super-Kamiokande], Phys. Rev. Lett. 130, no.3, 031802 (2023) doi:10.1103/PhysRevLett.130.031802 [arXiv:2209.14968 [hep-ex]].
  * [53] X. Cui et al. [PandaX-II], Phys. Rev. Lett. 128, no.17, 171801 (2022) doi:10.1103/PhysRevLett.128.171801 [arXiv:2112.08957 [hep-ex]].
  * [54] T. Bringmann and M. Pospelov, Phys. Rev. Lett. 122, no.17, 171801 (2019) doi:10.1103/PhysRevLett.122.171801 [arXiv:1810.10543 [hep-ph]].
  * [55] R. Plestid, V. Takhistov, Y. D. Tsai, T. Bringmann, A. Kusenko and M. Pospelov, Phys. Rev. D 102, 115032 (2020) doi:10.1103/PhysRevD.102.115032 [arXiv:2002.11732 [hep-ph]].
  * [56] P. K. Hu, A. Kusenko and V. Takhistov, Phys. Lett. B 768, 18-22 (2017) doi:10.1016/j.physletb.2017.02.035 [arXiv:1611.04599 [hep-ph]].
  * [57] G. Elor, R. McGehee and A. Pierce, Phys. Rev. Lett. 130, no.3, 031803 (2023) doi:10.1103/PhysRevLett.130.031803 [arXiv:2112.03920 [hep-ph]].
  * [58] A. Ruzi, C. Zhou, X. Sun, D. Wang, S. Wang, Y. Ban, Y. Mao, Q. Li and Q. Li, Int. J. Mod. Phys. A 38 (2023) no.29n30, 2350154 doi:10.1142/S0217751X23501543 [arXiv:2303.18117 [hep-ph]].
  * [59] D. McKeen, D. E. Morrissey, M. Pospelov, H. Ramani and A. Ray, Phys. Rev. Lett. 131, no.1, 011005 (2023) doi:10.1103/PhysRevLett.131.011005 [arXiv:2303.03416 [hep-ph]].
  * [60] F. Sauli, Nucl. Instrum. Meth. A 386, 531-534 (1997) doi:10.1016/S0168-9002(96)01172-2
  * [61] CMS Collaboration, A. M. Sirunyan, et al., Tech. Design Rep. CERN-LHCC-2017-012, CMS-TDR-016 (2017) [https://cds.cern.ch/record/2283189].
  * [62] Ju Xu-Dong, Dong Ming-Yi, Zhou Chuan-Xing, Dong Jing, Zhao Yu-Bin, Zhang Hong-Yu, Qi Hui-Rong, Ouyang Qun, Acta Phys. Sin., 2017, 66(7): 072902 doi: 10.7498/aps.66.072902
  * [63] Qite Li and Yanlin Ye, et al. Nucl. Instrum. Methods. Phys. Res. A 663 (2012) 22-25
  * [64] Li Q.T, et al. Chinese Phys. C37(2013)016002
  * [65] Chen S.D., et al. JINST9(2014)C10022
  * [66] Yang,LS., Xu,JY., Li,QT.et al. NUCL.SCI.TECH 32,85(2021)
  * [67] Cook, S. and D’Arcy, R. and Edmonds, A.et al. Phys.Rev.Accel.Beams20,3,030101,10(2017) doi:10.1103/physrevaccelbeams.20.030101 [arXiv:1610.07850v1 [physics.acc-ph]]
  * [68] M. Abbas, M. Abbrescia, H. Abdalla, A. Abdelalim, S. AbuZeid, A. Agapitos, A. Ahmad, A. Ahmed, W. Ahmed and C. Aimè, et al. Nucl. Instrum. Meth. A 1034, 166716 (2022) doi:10.1016/j.nima.2022.166716 [arXiv:2203.12037 [physics.ins-det]].
  * [69] A. Pellecchia et al. [CMS Muon], Nucl. Instrum. Meth. A 1046, 167618 (2023) doi:10.1016/j.nima.2022.167618 [arXiv:2207.09906 [physics.ins-det]].
  * [70] S. Agostinelli et al. (GEANT4 Collaboration), Nucl. Instrum. Meth. A 506, 250-303 (2003)
  * [71] https://ep-news.web.cern.ch/content/gems-cms
  * [72] J. D. Lewin and P. F. Smith, Astropart. Phys. 6 (1996), 87-112 doi:10.1016/S0927-6505(96)00047-3
  * [73] J. Aalbers, en. 9789462339842, 2018. ISBN: 978-94-6233-984-2
  * [74] C. Hagmann, D. Lange, D. Wright IEEE Nuclear Science Symposium Conference Record.IEEE, doi:10.1109/NSSMIC.2007.4437209.
  * [75] L. J. Schultz, K. N. Borozdin, J. J. Gomez, G. E. Hogan, J. A. McGill, C. L. Morris, W. C. Priedhorsky, A. Saunders and M. E. Teasdale, Nucl. Instrum. Meth. A 519, 687-694 (2004) doi:10.1016/j.nima.2003.11.035
  * [76] L. J. Schultz et al., 2006 IEEE Nuclear Science Symposium Conference Record, San Diego, CA, USA, 2006, pp. 2574-2577, doi: 10.1109/NSSMIC.2006.354434.
  * [77] L. J. Schultz et al., IEEE Transactions on Image Processing, vol. 16, no. 8, pp. 1985-1993, Aug. 2007, doi: 10.1109/TIP.2007.901239.
  * [78] L. Cheng-Ming, W. Qun-Gang, Z. Zhi-Yong and H. Guang-Shun, Radiat. Detect. Technol. Methods 4, 263-268 (2020) doi:10.1007/s41605-020-00179-9 [arXiv:1907.08928 [physics.ins-det]].
  * [79] G. Punzi, eConf C030908, MODT002 (2003) [arXiv:physics/0308063 [physics]].
  * [80] https://cms-analysis.github.io/HiggsAnalysis-CombinedLimit/
  * [81] T. Junk, Nucl. Instrum. Meth. A 434 (1999), 435-443 doi:10.1016/S0168-9002(99)00498-2
  * [82] A. L. Read, J. Phys. G 28 (2002), 2693-2704 doi:10.1088/0954-3899/28/10/313
  * [83] A. Bortone [CGEM-IT Working Group], Nucl. Instrum. Meth. A 1048, 167957 (2023) doi:10.1016/j.nima.2022.167957 [arXiv:2212.10896 [physics.ins-det]].
  * [84] J. Y. Tang, S. N. Fu, H. T. Jing, H. Q. Tang, J. Wei and H. H. Xia, Chin. Phys. C 34, 121-125 (2010) doi:10.1088/1674-1137/34/1/022
  * [85] Y. Bao, J. Chen, C. Chen, H. Cheng, C. Deng, R. Fan, Y. Guo, N. He, H. Hu and Q. Li, et al. J. Phys. Conf. Ser. 2462, no.1, 012034 (2023) doi:10.1088/1742-6596/2462/1/012034
  * [86] R.L. Workman et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2022, 083C01 (2022) and 2023 update.
  * [87] R. Essig, A. Manalaysay, J. Mardon, P. Sorensen and T. Volansky, Phys. Rev. Lett. 109, 021301 (2012) doi:10.1103/PhysRevLett.109.021301 [arXiv:1206.2644 [astro-ph.CO]].
  * [88] M. Battaglieri, P. Bisio, M. Bondí, A. Celentano, P. L. Cole, M. De Napoli, R. De Vita, L. Marsicano, G. Ottonello and F. Parodi, et al. Eur. Phys. J. C 81, no.2, 164 (2021) doi:10.1140/epjc/s10052-021-08957-5 [arXiv:2011.10532 [physics.ins-det]].
  * [89] A. Berlin, P. deNiverville, A. Ritz, P. Schuster and N. Toro, Phys. Rev. D 102, no.9, 095011 (2020) doi:10.1103/PhysRevD.102.095011 [arXiv:2003.03379 [hep-ph]].
  * [90] T. Åkesson et al. [LDMX], [arXiv:1808.05219 [hep-ex]].
  * [91] P. Agnes et al. [DarkSide], Phys. Rev. Lett. 130, no.10, 101002 (2023) doi:10.1103/PhysRevLett.130.101002 [arXiv:2207.11968 [hep-ex]].
  * [92] H. Al Ali, N. Arkani-Hamed, I. Banta, S. Benevedes, D. Buttazzo, T. Cai, J. Cheng, T. Cohen, N. Craig and M. Ekhterachian, et al. Rept. Prog. Phys. 85, no.8, 084201 (2022) doi:10.1088/1361-6633/ac6678 [arXiv:2103.14043 [hep-ph]].
  * [93] B. Abi et al. [Muon g-2], Phys. Rev. Lett. 126, no.14, 141801 (2021) doi:10.1103/PhysRevLett.126.141801 [arXiv:2104.03281 [hep-ex]].
  * [94] D. Forbes, C. Herwig, Y. Kahn, G. Krnjaic, C. Mantilla Suarez, N. Tran and A. Whitbeck, [arXiv:2212.00033 [hep-ph]].
  * [95] S. N. Gninenko, N. V. Krasnikov and V. A. Matveev, Phys. Rev. D 91, 095015 (2015) doi:10.1103/PhysRevD.91.095015 [arXiv:1412.1400 [hep-ph]].
  * [96] C. Y. Chen, J. Kozaczuk and Y. M. Zhong, JHEP 10, 154 (2018) doi:10.1007/JHEP10(2018)154 [arXiv:1807.03790 [hep-ph]].
  * [97] Y. M. Andreev, D. Banerjee, B. B. Oberhauser, J. Bernhard, P. Bisio, N. Charitonidis, P. Crivelli, E. Depero, A. V. Dermenev and S. V. Donskov, et al. [arXiv:2401.01708 [hep-ex]].
  * [98] Y. Kahn, G. Krnjaic, N. Tran and A. Whitbeck, JHEP 09, 153 (2018) doi:10.1007/JHEP09(2018)153 [arXiv:1804.03144 [hep-ph]].
  * [99] C. Y. Chen, M. Pospelov and Y. M. Zhong, Phys. Rev. D 95, no.11, 115005 (2017) doi:10.1103/PhysRevD.95.115005 [arXiv:1701.07437 [hep-ph]].
  * [100] C. Cesarotti, S. Homiller, R. K. Mishra and M. Reece, Phys. Rev. Lett. 130, no.7, 071803 (2023) doi:10.1103/PhysRevLett.130.071803 [arXiv:2202.12302 [hep-ph]].