aainstitutetext: University of Washington, Seattle, WA, USAbbinstitutetext:
Instituto de Física y Matemáticas, Universidad Michoacana, Morelia, Michoacan,
Mexico ccinstitutetext: Istituto Nazionale di Fisica Nucleare, Sezione di
Bologna, Bologna, Italyddinstitutetext: Universidad Autónoma de Chiapas,
Tuxtla Gutiérrez, Chiapas, Mexicoeeinstitutetext: Istituto Nazionale di Fisica
Nucleare, Sezione di Roma Tor Vergata, Roma, Italyffinstitutetext: Università
degli Studi di Roma Tor Vergata, Roma, Italygginstitutetext: Benemérita
Universidad Autónoma de Puebla, Puebla, Mexicohhinstitutetext: Universidad
Mayor de San Andrés, La Paz, Boliviaiiinstitutetext: Universidad San Francisco
de Quito USFQ, Quito, Ecuadorjjinstitutetext: Escuela Nacional de Ciencias de
la Tierra, UNAM, Ciudad de Mexico, Mexico**institutetext: Corresponding
author. Email address<EMAIL_ADDRESS>
# Cosmic-ray searches with the MATHUSLA detector

mathusla-experiment.web.cern.ch
a Cristiano Alpigiani b,* , J. C. Arteaga-Velázquez g , Daniela Blanco-Lira c
, Davide Boscherini d , K. S. Caballero-Mora e,f , Paolo Camarri e , Roberto
Cardarelli i , Dennis Cazar Ramírez e , Giuseppe Di Sciascio g , Arturo
Fernández Téllez a , H. J. Lubatti d,j , O. G. Morales-Olivares h , Piter
Amador Paye Mamani b , David Rivera Rangel g , Mario Rodríguez-Cahuantzi e,f ,
Rinaldo Santonico h , Martin Alfonso Subieta Vázquez

###### Abstract

The performance of the proposed MATHUSLA detector as an instrument for
studying the physics of cosmic rays by measuring extensive air showers is
presented. The MATHUSLA detector is designed to observe and study the decay of
long-lived particles produced at the pp interaction point of the CMS detector
at CERN during the HL-LHC data-taking period. The proposed MATHUSLA detector
will be composed of many layers of long scintillating bars that cannot measure
more than one hit per bar and correctly report the hit coordinate in case of
multiple hits. This study shows that adding a layer of RPC detectors with both
analogue and digital readout significantly enhances the capabilities of
MATHUSLA to measure the local densities and arrival times of charged particles
at the front of air showers. We discuss open issues in cosmic-ray physics that
the proposed MATHUSLA detector with an additional layer of RPC detectors could
address and conclude by comparing with other air-shower facilities that
measure cosmic rays in the PeV energy range.

###### Keywords:

MATHUSLA, cosmic rays, multi-messenger astronomy, extensive air showers

##  1 Introduction

This paper describes a study of the potential contributions of the proposed
MATHUSLA detector to the understanding of the origin and nature of cosmic rays
by means of extended air showers (EAS). MATHUSLA is designed to search for new
physics beyond the standard model, in particular, decays from neutral long-
lived particles (LLPs). The detector is to be located on the surface above the
CMS (Compact Muon Spectrometer) experiment at the CERN LHC. The proposed
MATHUSLA detector has nearly 100$\%$ coverage for charged particles over 100
m2 Chou:2016lxi ; Alpigiani:2018fgd ; Alpigiani:2022Snowmass . We show in this
paper that by adding a layer of RPC detectors, the MATHUSLA detector can also
make important contributions to many unresolved issues in cosmic-ray physics.

In particular, in the next sections, we present the results of a Monte Carlo
(MC) study to analyze the sensitivity of MATHUSLA as a detector of cosmic-ray
showers. We show that the addition of a layer of RPCs combined with the many
layers of scintillating detector planes can provide unique spatial and
temporal measurements of EAS. We discuss that such a detector provides
information for studies of the energy spectrum, arrival-direction distribution
and composition of cosmic rays that can be used for the research on the origin
and nature of these particles and for tests of hadronic interaction models. We
also discuss the potential for studying inclined events.

Section 2 gives a summary of open questions in the study of cosmic rays and
Section 3 a description of the detector as designed to address the particle
physics goals. In Section 4 we discuss the addition of a layer of RPC
detectors that would allow for detailed measurements of EAS. Section 5
discusses the MC simulations of the production of EAS and the detector
response, Section 6 the event reconstruction methodology, and Section 7, the
MC event selection. Results of this study are presented in Section 8 and
discussed in Section 9. We present our conclusions in Section 10.

##  2 A glance at cosmic-ray physics

The Earth is being continuously hit by high-energy subatomic particles known
as cosmic rays. They are dominated by protons and can have energies as low as
some MeV and as high as $10^{20}\mathrm{eV}$, i.e. orders of magnitude beyond
the energies of the most powerful accelerators on Earth CosmicRays ;
CosmicRays1 ; CosmicRays2 . Their total energy spectrum can be described by a
broken power-law, with a softening at around $4\times 10^{15}$ eV (where the
spectral index $\gamma$ changes from $-2.7$ to $-3$) and a hardening close to
$5\times 10^{18}$ eV (at which the spectral index increases to $-2.6$)
CosmicRays2 . Relative abundances of elemental nuclei in the intensity of
cosmic rays are difficult to measure. In spite of this, it has been found that
the composition of cosmic rays evolves with the primary energy CosmicRays2 ;
PDG . In the region below the knee, protons and He nuclei dominates the
composition of cosmic rays, while at higher energies the heavier nuclei are
more abundant. This happens up to $\sim 10^{17}$ eV, where the composition of
cosmic rays becomes lighter. This behavior stops at energies close to the
ankle, where an additional transition to the heavy mass group is observed. The
evolution of one mass group to another one seems to be correlated with the
presence of fine structures in the energy spectrum of cosmic rays HAWCCR23 ;
Antoni2005 ; Apel2013 ; PAO_Composition_2023 ; PDG , which, in turn, could be
related with acceleration Peters61  and propagation Giacinti15  issues or the
transition between different source populations of cosmic rays Liu2022 ;
PAO_Composition_2023 .

The acceleration of cosmic rays may occur in astrophysical plasmas with strong
shock waves by means of the so called first-order Fermi acceleration mechanism
Axford77 ; Krymsky77 ; Bell78 ; Blandford78 , which predicts power-law spectra
for the intensity of cosmic rays. The propagation in space of these energetic
particles is mainly diffusive due to interactions with the random component of
magnetic fields in the interstellar and extragalactic media. As a consequence,
the angular distribution of the arrival directions of cosmic rays in the sky
seems to be mostly isotropic. There are deviations from isotropy in the sky
maps of cosmic rays at the level of $10^{-3}$ and $10^{-4}$, which creates
large- and small-scale anisotropic distributions in the sky Deligny2016 ;
Ahlers2017 ; Deligny2019 ; Ahlers2019 . The large-scale anisotropies are
thought to be due to propagation effects and the global distribution of the
sources of cosmic rays in the galaxy and in our galactic neighborhood. The
details behind the origin of the small-scale anisotropic structures are
however less clear. At TeV energies, they may be related with the heliosphere
or magnetic turbulences. At tens of EeV, they could be due to some
extragalactic sources in the direction of Centaurus and the Ursa Major, as
well as the Perseus-Pisces region PAO_TA_Aniso23 . A further implication of
the diffusive propagation of cosmic rays in the universe is that they loss
information about the position of their sources in the sky. Only at extreme
energies, the small magnetic deflections of these particles in the space may
allow to perform cosmic-ray astronomy CosmicRays2 . At energies above
$10^{18}$ eV, cosmic rays are of extragalactic origin CosmicRays3 , while
below $10^{17}$ eV, of galactic nature CRModels2019 . Somewhere, between
$10^{17}$ and $10^{18}$, a transition between galactic and extragalactic
cosmic rays is expected Alosio12 ; PAO_Composition_2023 . In our galaxy,
supernova remnants could be one of the astrophysical sources responsible for
cosmic rays of high-energy CosmicRays3 . However, there are some clues that
the galactic center CosmicRays4  and star forming regions HAWC_Cygnus_2021 ;
LHAASO_ICRC2023_010  could also accelerate cosmic rays of very high energies
in the Milky Way HAWC_Cygnus_2021 . At ultra high-energies, starburst galaxies
are promising candidates PAO_POSICRC2023_252 .

The existence of cosmic rays faces us with several compelling questions about
the universe at very high energies. The quest for satisfactory answers to
these phenomena is pushed to the limit by the combined observations of cosmic
rays, gamma rays, cosmic neutrinos and gravitational waves MultimessAJL2017 ;
MultimessAA2017 . The astroparticle messenger approach calls for scenarios
that involve both the domain of elementary particles, the so-called ”micro-
cosmos”, and the realm of the huge celestial objects, what we know as ”macro-
cosmos”.

There are several open questions as to the astroparticle messengers phenomena
such as:

  * •

What is the origin of the cosmic-ray particles of extraordinary energy?

  * •

What is the mass composition of cosmic rays?

  * •

What are the phenomena behind the small- and large-scale cosmic-ray
anisotropies?

  * •

How are cosmic-rays accelerated to high energies?

When cosmic rays of high energies arrive to the Earth, they collide with the
molecules in the atmosphere at altitudes that goes from $\sim 15$ up to $35$
km EASmodelsReview  initiating a set of chain reactions that leads to the
production of a cascade of particles that travels at relativistic velocities
to the ground. This extensive air shower, as it is known, contains a large
number of particles, which depends on the primary energy of the cosmic ray and
increases according to a power-law function. The air shower has
electromagnetic, penetrating and hadronic components. The first one is most
abundant and is composed of photons and $e^{\pm}$’s, the second one, by
$\mu^{\pm}$’s, and the last one, by nucleons, mesons, and atomic nuclei.
During their development in the air, EAS produce radio and Cherenkov radiation
CosmicRays . Also, the atmosphere emits fluorescence light from the excitation
of nitrogen molecules CosmicRays2 .

At energies above $10^{13}$ eV, cosmic rays are detected indirectly from the
measurements of EAS using arrays of particle detectors, radio antennas,
Cherenkov and fluorescence telescopes that are displayed over large surfaces
CosmicRays . At energies smaller than $10^{15}$ eV, cosmic-ray studies are
observed in a direct manner from balloon-borne and satellite instruments
CosmicRays2 ; DirectCR23 .

Data collected by many ground and underground-based observatories as well as
by sensitive radiation detectors installed on balloons and also on satellites
orbiting Earth CosmicRays13 ; CosmicRays14 ; CosmicRays15 ; CosmicRays16 ;
CosmicRays17 ; CosmicRays18 ; CosmicRays19  have provided much information on
the energy spectrum, arrival direction and mass composition of cosmic rays
that have revealed very important clues about the nature, propagation,
acceleration mechanism and the origin of these energetic particles CosmicRays3
; CosmicRays5 ; CosmicRays6 ; CosmicRays7 ; CosmicRays8 ; CosmicRays9 ;
CosmicRays10 . In addition, observations from the northern hemisphere of high-
energy gamma rays reveals many potential sources of cosmic rays within our
galaxy CosmicRays4 ; CosmicRays11 ; CosmicRays12 ; Lhaaso2021 .

Notwithstanding these excellent measurements, the abovementioned questions
still remain without conclusive answers. They will be addressed by a new
generation of ground-based observatories, some of which are located at high
altitudes ($\sim 4500\,\mbox{m}$ a.s.l) in the Earth’s southern CosmicRays20 ;
CosmicRays21  and northern hemispheres CosmicRays22 , and by new space-based
instruments Olinto21 . In addition, multiwavelength and multimessenger
observations will play an important role in these tasks MultimessAJL2017 ;
MultimessAA2017 . We demonstrate in this paper that the MATHUSLA detector
could also make important contributions to advancing the understanding of the
astrophysics and physics of high-energy cosmic rays.

Figure 1: MATHUSLA detector at CMS site Alpigiani:2022Snowmass .

##  3 The MATHUSLA Detector

The proposed MATHUSLA detector shown in Figure 1 will be located on the
surface above the CMS CMS  $pp$ interaction point (IP) located on CERN-owned
land. The proposed goal is to have MATHUSLA ready to take data when the High-
Luminosity LHC (HL-LHC) HL-LHC  turns on.

To have sufficient solid-angle coverage, the current detector concept, shown
in Fig. 2, is a 100 $\times$ 100 m2 detector consisting of a hundred
$9\,\textrm{m}\times 9\,\textrm{m}$ units. Each detector unit comprises
9-layers of scintillating-detector planes that provide position and timing
coordinates of charged particles resulting from the decay of long-lived
particles in the MATHUSLA detector decay volume. There are five detector
planes at the top, two additional planes 5 m below (to enhance the particle
position measurement precision close to the floor), and two additional planes
at the floor (to veto charged particles from the LHC and cosmic-muon
backscattering). The total height of $\sim 40$ m includes a $\sim 25$ m LLP
decay volume, 21 m of which is below ground, and 17 m above ground that houses
the tracker and the crane system used for assembly and maintenance.

Figure 2: Schematic view of the MATHUSLA detector modular concept
Alpigiani:2022Snowmass .

The layout of the building that houses the 100 m $\times$ 100 m experimental
area and the 30 m $\times$ 100 m adjacent area for the detector assembly is
shown in Figure 3. The structure, which is located on the surface near the CMS
IP fits well on CERN-owned land. Having a large part of the decay volume
underground brings it closer to the IP, which increases the solid angle in the
acceptance of LLPs generated in the collisions. To adjust to the available
land, this proposal has a 7.5 m offset to the centre of the beams. The site
allows for the detector to be as close as 68 m from the IP, which is shown in
red in Figure 3

Figure 3: MATHUSLA detector at CMS site Alpigiani:2022Snowmass .

###  3.1 Scintillating Detector Planes

Figure 4: Extruded scintillating bars produced at FNAL with one and two holes
for WLS fiber insertion.

Each of the $9\,\textrm{m}\times 9\,\textrm{m}$ detector planes consists of an
assembly of extruded scintillating bars whose length, width and thickness are
typically $4.6$ m, $4.5$ cm, $2$ cm, respectively. Each bar is extruded with a
hole at the centre into which a wave-length shifting (WLS) fibre is inserted
and connected at each end to a SiPM. The coordinate along the length of the
bar is determined by differential time measurement at the two ends of the bar
with resolution $\sigma\sim\pm$ $15$ cm. The width of the bar determines the
corresponding transverse coordinate with $\sigma=\pm$ $1.3$ cm. Figure 4 shows
some typical extruded scintillating bars produced at the Fermi National
Accelerator Laboratory (FNAL).

To facilitate installation the scintillating bars are assembled into $4$ sub-
units that comprise $\sim 100$ bars resulting in $\sim 4.5\,\textrm{m}\times
4.5\,\textrm{m}$ sub-units that allow for overlapping the $4$ sub-units
longitudinally and transversely to avoid gaps in coverage. In this arrangement
the $4$ sub-units provide hermetic surface-coverage over a $\sim
9\,\textrm{m}\times 9\textrm{m}$.

###  3.2 Cosmic Rays

The dominant background to long-lived particle searches in the MATHUSLA
detector is cosmic rays. Timing information from the scintillating bars allows
for the rejection of downward-going charged cosmic rays. The two-floor
detectors provide rejection of inelastic back-scattering of cosmic-ray muons
striking the surface and muon decays that send an upward electron through the
detector planes. They also serve to reject inelastic interactions of muons
near the surface. Several simulation studies show that these backgrounds can
be rejected. A more detailed description of the possible backgrounds can be
found in Alpigiani:2018fgd .

MATHUSLA, equipped with seven scintillating detector planes at the top of the
instrument that has good space and time resolution (see section 3.1) and
nearly continuous coverage over an area of $10^{4}$ $\mbox{m}^{2}$, is also
sensitive to air showers induced by cosmic rays. However, it is limited as an
air shower detector because the long scintillating bars are not able to
provide coordinates of several charged tracks hitting the same bar at nearly
the same time. Such occurrences lead to saturation and loss of information.

##  4 The RPC detector

###  4.1 Introduction

The peculiarity of MATHUSLA, as a cosmic-ray experiment, is that its detector
is primarily conceived for searching for very long-lived particles produced by
pp interactions in the LHC. Studying its performance as a cosmic-ray detector,
and substantially improving it with a special sub-detector dedicated to the
study of extensive cosmic-ray showers, is discussed in this section.

The main constraints imposed by the primary purpose of the experiment are the
operation at CERN’s very modest altitude ($373.6$ m.a.s.l.) and the sensitive
area not exceeding $10^{4}$ m2. On the other hand, MATHUSLA, equipped with
$10$ scintillating planes that provide good spatial resolution and one ns time
resolution is an unprecedented large area cosmic-ray tracker.

Taking into account the features described above the cosmic-ray search of
MATHUSLA may be focused on the following main items:

  * •

Cosmic-ray composition, i.e. measurement of the atomic number Z of the primary
particle.

  * •

Parallel-muon bundles crossing the detector. In this case, a “pure muonic
shower” is observed. For inclined or almost horizontal showers that traverse a
larger thickness of the atmosphere where the hadronic component is absorbed,
purer muonic showers can be observed at high and very high energies.

However, a scintillating tracker such as MATHUSLA has the following limits for
detecting EAS:

  * •

Limited resolution for multiple-track detection. The limit is no more than
approximately $1$ track in a scintillating bar of area ($450\times 4.5$) cm2 =
$2025$ cm2. Multiple tracks in the same scintillating bar, in addition to
being counted as a single track, also spoil the timing information and do not
allow localization of the track position along the scintillating bar. This
severely limits the study of parallel-muon bundles, which is an excellent
field of investigation where MATHUSLA could be competitive with other cosmic-
ray experiments.

The MATHUSLA scintillating bars do not provide energy deposition measurements
that are extremely important for the shower-core detection where the hits
produced by electro-photonic component can reach a density of $10^{6}$ m-2 for
a primary cosmic-ray particle of $1$ PeV.

  * •

Absence of an obvious muon-identification strategy, which is normally based on
a heavy-material layer absorbing secondary particles except muons. This is
impossible for MATHUSLA and requires a different strategy. It should be
stressed that muon detection is important in any channel of cosmic-ray search,
but it is crucial in particular for gamma-hadron discrimination. A possible
muon-detection alternative strategy could be based on the identification of
the regular linear pattern of a muon track from the random hits produced by
the electro-photonic component of the shower. This could be possible thanks to
the tracking capability of MATHUSLA. The expected efficiency of this method in
identifying a muon would be modest. However, at about $1$ PeV a proton shower
produces an order of $10^{4}$ muons, and the identification of just one of
them would be sufficient to reject the shower as a gamma ray candidate.

###  4.2 The RPC sub-detector

The addition of a single RPC layer to the layers of MATHUSLA scintillating
detectors can significantly enhance the capability of the instrument to make
detailed measurements of EAS. In this section, we describe the RPC chambers,
the proposed layout for MATHUSLA, and the front-end electronics. The proposed
RPC configuration profits from the experience obtained during 5 years of
operation of the ARGO-YBJ experiment AIELLI200692 . The main difference
compared to ARGO-YBJ is that the RPC chambers for MATHUSLA will operate in
avalanche mode instead of the streamer mode chosen for ARGO-YBJ. Operating in
avalanche mode requires improved front-end electronics to compensate for the
smaller charge and lower discrimination threshold, which has the advantage of
substantially extending the linear response range of the analogue readout
BARTOLI201547 .

The gas-gap thickness is another difference with respect to ARGO-YBJ. At the
ARGO-YBJ altitude of $4300$ m a.s.l., the atmospheric pressure was about $600$
mbar to be compared with the CERN pressure of $980$ mbar. This makes possible
a substantial reduction of the gap with respect to the $2$-mm gap of the ARGO-
YBJ RPCs. We are testing $1$-mm gas gaps that will be used for the BI RPCs of
the ATLAS muon spectrometer ATLAS-TDR-026 . The choice of the $1$-mm standard
for the gap size of the MATHUSLA RPCs allows us to exploit the research and
development of the ATLAS Phase II upgrade. A further advantage of reducing the
gap size to $1$ mm is improved time resolution. The MATHUSLA RPC gas gap is
schematically shown in Fig. 5

Figure 5: RPC gas gap sandwiched between the strip and the “big pad” panels.
The main components are the high-pressure laminate electrode plates, the
graphite electrodes, the insulating PET foil, and the spacers.

The main RPC parameters that can be easily optimized concern the read-out
granularity. Following the ARGO-YBJ scheme there are two types of read-out
electrodes: the digitally readout strips located on top of the gas gap in Fig.
5 and and the “big pads”, located below the gas gap in Fig. 5. The pads are
large area capacitors, covering half of the gap area, and integrate the charge
of all hits located in the corresponding pad area. The analog readout of the
pad signal makes it possible to measure the shower-core multiplicity when it
is so high that the digital readout is saturated. The strips, longitudinally
oriented, are $2.2$ m long with an $11.2$-mm pitch with front-end electronics
located at both ends in order to determine the hit position along the strip
using the signal arrival time at the two ends. The area of the strips, $242$
cm2, is substantially smaller than the area of the scintillating bars, $2025$
cm2, which significantly reduces the probability of multiple hits in the same
strip. This is particularly important for the detection of parallel muon
bundles where the probability of two or more muons striking the same
scintillating bar is large as well as for selecting straight tracks from the
background of randomly distributed hits due to the electro-photonic component
of the shower.

The $0.99$ m2“big pads” ($1.1\times 0.9$) m2), are smaller than those of ARGO-
YBJ ($1.75$ m2) and provide better detail of the shower-core profile. This is
important for identifying showers generated by different kinds of primary
cosmic rays. For example, gamma showers have a very compact single core
compared to the less compact multiple-core structure of hadron showers.
Moreover, the high-Z primary showers are broader than proton showers and have
a substantially larger multiplicity of different cores.

These qualitative considerations need to be studied in detail with an adequate
simulation in order to optimize the granularity of both strips and “big pads”,
both of which can be easily tuned for specific physics channels.

The RPC characteristics are the following:

  * •

Gas gap: a single 1-mm gap of area ($2.2\times 0.9$) m2 with $2$ gas inlets
and $2$ outlets at the corners.

  * •

Strip read-out panels: $80$ longitudinal strips $2.2$ m long with $11.2$-mm
pitch. Front-end boards are located at both ends of the strip for
interpolating the position along the strip. The alternative solution with the
strips oriented along the short side of the chamber can also be considered but
it is less practical for the final assembly of the chamber.

  * •

“Big pads” for the analogue read-out of area ($1.1\times 0.9$) m2

  * •

Front-end electronics embedded in the same Faraday cage of the detector

  * •

Mechanical support structure made with honeycomb panels about $1$ cm thick.

In order to fulfill the requirements of the 9$\times$9 m2 super-modulus, 4
chambers are assembled in the same mechanical structure of 9$\times$0.9 m2
size, as shown in Fig. 6. The gas is serially distributed from one gap to the
next. The structure supporting the $4$ chambers is made of two “C”-shaped
aluminum profiles fixed at the sides of the chambers (as shown in Fig. 7).
This structure behaves as a stiff beam supported only at the ends, with a
sagitta of about 1-cm. The thickness and the height of the profiles need to be
optimized taking into account the weight of the structure. In the final
assembly of a $9\,\times\,9$ m2 super-modulus, two “C”-shaped profiles in
contact with one another are screwed together in a single “I”-shaped profile.

Figure 6: Schematic of one RPC unit with four chambers assembled in a line
within the same mechanical structure. Figure 7: Schematic side view of one RPC
unit showing details of its mechanical structure with the two “C”-shaped
profiles. All the values displayed in the figure are in mm.

###  4.3 The RPC front-end electronics

The front-end circuit was developed at the INFN laboratory of Rome Tor Vergata
for the RPC detectors of the ATLAS LHC phase-II upgrade ATLAS-TDR-026 . The
circuit consists of a charge preamplifier, a discriminator, and a Low Voltage
Differential Signal (LVDS) driver, shown schematically in Fig. 8. The
amplifier and the driver are built as discrete component circuits. The
discriminator is a full-custom circuit, conceived with the purpose of
optimizing performance and allowing for large mass production at low cost

Figure 8: Logical layout of an 8-channel RPC front-end circuit board.

####  4.3.1 Preamplifier

The preamplifier uses Silicon Bipolar Junction Transistor (BJT) technology.
The main feature of this new kind of preamplifier is fast charge integration
with the possibility of matching the input impedance to a transmission line.
Its working principle is shown in Fig. 9, which explains how the injected
signal is integrated.

Figure 9: Comparison between a bare RPC signal (upper plot) and the
corresponding amplified signal (lower plot).

The specifications of the Silicon BJT technology planned for the charge
preamplifier are listed in Table 1.

Table 1: Specifications for the preamplifier in the RPC front-end circuit. Voltage supply |  $3-5$ V
---|---
Sensitivity |  $0.2-0.4$ mV/fC
Noise (independent of the detector) |  $1000$ e- RMS
Input impedance |  $50-100$ $\Omega$
Bandwidth |  $10-100$ MHz
Power consumption |  $5$ mW/ch
Rise time for $\delta$(t) input |  $300-600$ ps

####  4.3.2 Discriminator

The requirements for the MATHUSLA cosmic-ray RPC discriminator are the
following: a) very low threshold, of the order of $6$-$7$ mV, and b) efficient
operation with very fast pulses like the ones produced by RPCs. The very low
threshold for cosmic shower detection is not required to improve the rate
capability, which is very modest in this case, but rather to improve the
linearity range in the measurement of the hit density inside the shower core.
The core of the showers can have hit densities exceeding $10^{6}$/m2, due to
the very high energy of the primaries. In this case, the linearity range
requires that the charge released in the gas by a single hit is very modest.
The above requirements are met with a new type of discriminator based on the
saturated amplifier: a signal above the discrimination threshold is amplified
by multiple stages up to saturation. The output is a square wave of fixed
amplitude independent of the input signal. Its duration is the time over the
threshold of the input signal, which allows measuring the input amplitude,
once the input shape is known. This type of measurement uses a logarithmic
scale. This is a crucial characteristic when the input charge can vary over a
wide range. In the case of MATHUSLA, the number of particles in a single
front-end channel may range over several orders of magnitude, depending on the
shower energy and the read-out strip area.

Figure 10: (a) Schematic of the discriminator circuit. (b) Signal-processing
scheme for the discriminator circuit.

The discriminator circuit is schematically shown in Fig. 10 (a). The first
stage is a simple amplifier in the BJT common-emitter configuration. The
following three stages can be represented as logic “NOT” gates and have the
task of discriminating (first logic stage) and saturating the input signal. As
shown in Fig. 10 (b), for each logic stage the output shows sharper edges and
the final output is a rectangular wave of inverse polarity with respect to the
analogue input. The discriminator features are summarised in Table 2.

Table 2: Specifications for the discriminator in the RPC front-end circuit. Technology |  $130$ nm Si-Ge BiCMOS
---|---
Voltage supply |  $1-2.5$ V
Minimum threshold |  $0.3$ $\mu$V
Minimum input-pulse width for threshold linearity |  $0.5$ ns
Bandwidth |  $10-100$ MHz
Power consumption | 10 mW/ch
Rise time for $\delta$(t) input | 300 ps
Input impedance |  $100$ $\Omega$
Double-pulse separation | 1 ns

####  4.3.3 TDC

In a system with a very large number of channels, where time measurements are
required, time digitization is one of the most relevant problems due to
cabling complexity and TDC cost. For the system proposed for MATHUSLA cosmic-
ray measurements, the TDC is embedded in the front-end electronics and is
included in the full-custom circuit under development. This is possible thanks
to new technologies, $120$-nm SiGe and CMOS, which make possible oscillators
and counting scalers that operate at $2-3$ GHz. For the detection of a large-
size cosmic shower, the intrinsic single-hit time fluctuations are in the
nanosecond range, and therefore a TDC of comparable accuracy would be suitable
for this purpose. This makes possible a drastic simplification for the TDC.
Indeed, at a frequency of $2-3$ GHz, the simple counting of the periods
covering the time interval under measurement results in an uncertainty of
$0.5-0.3$ ns that is adequate for cosmic-shower timing, and no interpolation
between two consecutive clock signals is required.

####  4.3.4 Trigger and data acquisition

This section describes the MATHUSLA cosmic-ray trigger and its implementation
in the electronic circuits that produce the trigger signal and give the Data
Acquisition the start signal. The trigger logic presented here is fully based
on the signals supplied by the RPC-layer detector, consequently, the cosmic
ray trigger is independent of the main long-lived particle trigger of the
experiment.

The MATHUSLA trigger of an EAS generated in the atmosphere by a cosmic primary
particle is based on the multiplicity of the shower hits detected by the RPC
layer. This multiplicity is approximately proportional to the EAS shower
energy on a logarithmic scale. Thus, a multiplicity threshold is equivalent to
an energy threshold, even within the large fluctuations of the energy vs.
multiplicity correlation. The main task of a multiplicity trigger is to select
the simultaneous hits due to an EAS from the accidental rate due to the
combinatorial background of the counting rate over a large-area detector.
Simultaneous hits due to a cosmic shower can be selected by a proper
coincidence among several chambers. The main challenge of this trigger is that
the time resolution of the coincidence must take into account the systematic
delay among different hits due to the inclination of the shower front with
respect to the vertical direction. For showers propagating almost horizontally
this delay is of the order of the linear detector dimension over the speed of
light, $L$/$c$. For a $(100\times 100)\,\mbox{m}^{2}$ detector like MATHUSLA,
this delay can be as large as $(100\sqrt{2}\,\mbox{m}$)/$c$, namely about
$450\,\mbox{ns}$, assuming the diagonal as the maximum distance between two
hits. This large delay does not allow for a coincidence time, which is the
main parameter for reducing the accidental rate. This problem can be solved by
dividing the detector into subdetectors of gradually increasing areas where
the coincidence time resolution can be improved according to the smaller
subdetector size. In this section, we give an example of how this trigger
concept can be implemented inside the MATHUSLA cosmic-ray detector and how the
accidental rate can be evaluated. A more precise trigger study, including the
efficiency vs. energy correlation as well as the trigger counting rate, for
each multiplicity threshold, requires a full simulation of the process
starting from the characteristics of the shower to be triggered.

Table 3: Specifications for the relevant RPC trigger parameters. Effective detection area for each submodule | $(9\times 9)\,\mbox{m}^{2}$
---|---
Assumed counting rate per unit area | $1\,\mbox{kHz}/\mbox{m}^{2}$
Multiplicity threshold for EAS detection | $50\,\mbox{hits}$
Number of RPC sub-modules | $100$
Number of chambers in one RPC sub-module | $40$
Trigger majority threshold in one sub-module | $2\,\mbox{hits}$
Time resolution for a submodule majority | $50\,\mbox{ns}$
Maximum time delay for inclined showers in one sub-module | $42\,\mbox{ns}$
RPC read-out strip length | $2.2\,\mbox{m}$
Inverse signal propagation speed | $5\,\mbox{ns}/\mbox{m}$

In order to evaluate the efficiency of rejecting the accidental background
rate and to support it with some numerical calculations we assume the
following working conditions:

  * •

Multiplicity threshold of $50$ hits on the whole RPC detector. The
corresponding energy threshold should be at least an order of magnitude higher
than that of ARGO-YBJ if renormalized to the much lower altitude of MATHUSLA.

  * •

Intrinsic detector counting rate equal to $1$ kHz/m2. This is an
overestimation of the expected rate, in order to obtain a very conservative
evaluation of the accidental background rate.

  * •

The hits due to the shower are randomly distributed over the detector. This is
also a conservative hypothesis for the accidental background calculation
because the shower hits tend to be clustered thus giving a higher probability
of several hits in the same sub-detector, which is a crucial point in the
following trigger project.

The relevant parameters for the MATHUSLA RPC trigger are listed in Table 3 and
will be used in the calculations reported below. The smallest subdetector
where the shower multiplicity is considered for trigger purposes is an RPC
sub-module of about $(9\times 9)$ m2. The MATHUSLA detector, according to Fig.
2, is subdivided into $100$ such units. Each unit contains an RPC plane of the
same area, consisting of $40$ chambers arranged in groups of 4 chambers in the
same mechanical support as shown in Fig. 6. The first step of the proposed
trigger requires selecting all submodules with at least $2$ hits due to the
shower. For each sub-module, the analog sum of all output signals from the
front-end circuits is sent to a discriminator with a threshold requiring at
least two simultaneous signals. The signal shaping time is $25$ ns, resulting
in a time resolution of about $50$ ns of this majority coincidence. This
covers the uncertainties due to both the shower inclination, $42$ ns for a
$(9\times 9)$ m2 sub-detector, and the timing propagation of the signal along
the readout strips (namely $2.2$ m $\times 5$ ns$/$m $=11$ ns) combined in
quadrature. The sub-module counting rate is $81$ kHz so that the accidental
coincidence rate is $(81\times 10^{3}$ Hz$)^{2}\times(50\times 10^{-9}$
s$)=328$ Hz for each of the $100$ submodules comprising the RPC detector.

The probability of missing the trigger, because all shower hits fall in
different sub-detectors, is extremely small. For a $50$ hits shower, the
average number of hits/subdetectors is $0.5$ and the expected number of towers
showing a coincidence of at least $2$ hits is $9.02$, according to Poisson
statistics.

The next step in the trigger logic is to require the coincidence with a proper
multiplicity on the whole detector. Repeating the procedure followed for the
single units, the analog sum of the $100$ tower outputs is sent to a
discriminator with a selectable threshold requiring the coincidence of two or
more towers on the whole detector. The total rate of the analog sum is $328$
Hz $\times 100=32.8$ kHz. The coincidence time, in this case, must be $450$ ns
to take into consideration the maximum time uncertainty due to the shower
inclination over the whole detector. For a multiplicity threshold of $2$
subdetectors out of a hundred, the accidental rate is: $(32.8\times 10^{3}$
Hz$)^{2}\times(450\times 10^{-9}$ s$)=427$ Hz which is a modest trigger rate.

The probability of missing the trigger for a shower of $50$ hits is
negligible. Indeed, the probability of less than two sub-detectors firing for
an expected number of $9.02$ is $1.23\times 10^{-3}$. If needed, the trigger
rate can further be decreased by increasing the multiplicity threshold above
$2$. For a threshold of $3$, the accidental rate is $(32.8\times 10^{3}$
Hz$)^{3}\times(450\times 10^{-9}$ s$)^{2}=7.1$ Hz.

####  4.3.5 Time schedule

The front-end circuit described here is not yet fully developed. The
preamplifier and the discriminator are operational and are included in the
RPCs currently being installed for the ATLAS Phase-I upgrade. The integration
of the TDC in the circuit is in progress. It requires $1$-$2$ further foundry
runs in the next $1$-$2$ years. The current plan is to be ready for the LHC
Phase-II upgrade.

##  5 Simulations

In order to explore the performance of the RPC layer and the scintillating
planes of MATHUSLA as air-shower cosmic-ray detectors, we have prepared a set
of toy MC simulations that include the production of EAS by the primary cosmic
rays, the evolution of the particle cascade until its arrival at the Earth’s
surface and the physical response of the detectors to the passage of the
shower particles. In the following subsections, we describe each of the above
simulation steps that are used to produce our MC data.

###  5.1 Air shower simulations

For the generation and development of EAS in the atmosphere, we followed the
procedure described in Alpigiani:2018fgd . In particular, we employed the
CORSIKA v7.71 package without thinning Heck:1998vt  together with the FLUKA
2020 Ferrari:2005zk ; BOHLEN2014211  and the QGSJET-II-04 qgsjet ; qgsjet2
hadronic interaction models. FLUKA was used to simulate collisions with hadron
energies $E_{h}<200\;\mathrm{GeV}$, while QGSJET-II-04 was employed to
describe hadronic interactions with higher energies. CORSIKA simulates the
evolution of the particle cascade through the atmosphere from the production
point to the detector level. The output provides information about the
spatial, time, and momentum distributions of the shower particles at the
observation point. Without thinning, CORSIKA follows the history of each
particle in the EAS until its energy falls below a certain threshold, which is
specified at the beginning of the simulation chain. Here, we kept the CORSIKA
default values of $3\,\mathrm{MeV}$ for shower electrons and gamma-rays, and
$100\,\mathrm{MeV}$ for shower muons and hadrons.

We chose the U.S. standard atmosphere model in our simulations. In addition,
we selected the curved atmosphere option, as we will work also with very
inclined showers. Moreover, the observation level was set at $436\,\mbox{m
a.s.l.}$ ($X\sim 1000\,\mbox{g}/\mbox{cm}^{2}$) and the local geomagnetic
field was chosen as $B_{x}=22.2\,\mu T$ and $B_{z}=41.9\,\mu T$ in the
simulations. Such values111The data was taken from the NOAA (National Centers
for Environmental Information) webpage at https://www.ngdc.noaa.gov/geomag-
web/. correspond to the characteristics of the site above the ATLAS
interaction point, which was initially considered as one of the potential
places of the MATHUSLA experiment. The differences with the location of
MATHUSLA at the CMS site are small (the altitude is about $374\,\mbox{m
a.s.l.}$ and the corresponding magnitudes of the geomagnetic field are
$B_{x}=22.2\,\mu T$, $B_{z}=42.1\,\mu T$). Therefore, our conclusions
regarding the performance of the MATHUSLA detector are not expected to be
affected by such differences.

The EAS were produced with a $\sin\theta\cos\theta$ distribution for two
zenith angle intervals: $\theta=0^{\circ}-20^{\circ}$ (vertical EAS) and
$\theta=70^{\circ}-80^{\circ}$ (inclined EAS), and the shower cores were
scattered randomly on the MATHUSLA area. A flat geometry was employed. We
produced MC simulations for H and Fe primary nuclei with energies between
$1\,\mathrm{TeV}$ and $1000\,\mathrm{PeV}$ assuming an energy spectrum of the
form $E^{-2}$. In order to improve the statistics at the highest energies, we
generated MC simulations by energy decade as shown in Table 4. Each event was
weighted to achieve a unique $E^{-2}$ energy spectrum for the whole energy
range.

Table 4: Distribution of the total number of CORSIKA simulations generated by energy decade. CORSIKA version 7.71 was used together with the FLUKA and QGSJET-II-04 hadronic interaction models.  | v7.71
---|---
Energy decade |  $0^{\circ}-20^{\circ}$ |  $70^{\circ}-80^{\circ}$
1-10 TeV |  $2\times 10^{4}$ H |  $2\times 10^{4}$ H
|  $2\times 10^{4}$ Fe |  $2\times 10^{4}$ Fe
10-100 TeV |  $2\times 10^{4}$ H |  $2\times 10^{4}$ H
|  $2\times 10^{4}$ Fe |  $2\times 10^{4}$ Fe
100-1000 TeV |  $1\times 10^{4}$ H |  $1\times 10^{4}$ H
|  $1\times 10^{4}$ Fe |  $1\times 10^{4}$ Fe
1-10 PeV |  $7.5\times 10^{2}$ H |  $1\times 10^{3}$ H
|  $7.5\times 10^{2}$ Fe |  $1\times 10^{3}$ Fe
10-100 PeV |  $1\times 10^{2}$ H |  $2.5\times 10^{2}$ H
|  $1\times 10^{2}$ Fe |  $2.5\times 10^{2}$ Fe
100-1000 PeV |  $3\times 10^{1}$ H |  $2.1\times 10^{2}$ H
|  $2\times 10^{1}$ Fe |  $1.5\times 10^{2}$ Fe

###  5.2 Simulation of the MATHUSLA detector

For our study, we built a toy MC model of MATHUSLA using a program written in
ROOT CERN rootcern . In particular, we considered a simple geometry consisting
of seven extruded scintillating detector layers and the RPC layer as shown in
Fig. 1. The present work did not incorporate the last design of MATHUSLA, but
a previous one. However, we do not foresee substantial modifications to the
main results on the performance of the experiment for cosmic-ray detection. In
our simulations, the two scintillating detector layers near the bottom of the
instrument are ignored. On the other hand, the RPC chamber plane is located
$2\,\mbox{m}$ above the intermediate scintillating detector layer that is
found in the middle between the top and intermediate groups of scintillating
detector planes. Hence, we have five detector layers above the RPC and two of
them below it (see Fig. 2). The size of the detector was set to
$100\,\mbox{m}\times 100\,\mbox{m}$. The coordinate system is oriented so that
the $X$ and $Y$ axis are parallel to the horizontal dimensions of the
detector. In this coordinate frame the detector plane extends from
$x,y=0\,\mbox{m}$ to $100\,\mbox{m}$, and from $z=0\,\mbox{m}$ to
$9\,\mbox{m}$. The height is measured from the bottom of the intermediate
tracking detector layer.

The RPC geometry used assumes Big Pads with a conservative size of
$1\,\mbox{m}\times 1\,\mbox{m}$. These elements provide not only spatial and
time data of the shower particles but also local particle density
measurements. We output the coordinates of the central positions of the Big
Pads that were hit, the arrival times of the first shower particles in each
module, and the induced signals. Space-time measurements of the EAS can be
made with individual RPC strips, but to keep the simulation as simple as
possible, we have focused our studies on the data from the Big Pads. The
amplitudes of the RPC-pad signals are assumed to be proportional to the
density of charged particles as found by the ARGO-YBJ RPC measurements
Aielli2012 . We must point out that in a real simulation, it must be
considered that the Big Pads can only be used to measure hit multiplicities
with high densities due to their high capacitance. For small density
measurements and individual hit counting, the RPC strips will be used.

For the simulations, the scintillating bars have a length of $5\,\,\mbox{m}$
and a width of $4\,\,\mbox{cm}$. They are oriented along the same direction on
each plane (parallel to one of the $X$ or $Y$ axis), but they are rotated by
$90^{\circ}$ when passing from one scintillating plane to another one. In the
MC simulations, we have considered that the scintillating detector layers are
sensitive only to the passage of charged particles and that the output
saturates for more than one hit. Our assumption is conservative. In practice,
the scintillating detector bars will saturate on average at one hit per bar.
Consequently, the simulation assumes that the scintillating bars provide
information of the impact point and arrival time only for hits of single-
charged particles. In case of more than one hit per bar, we only save the
central coordinates of the corresponding bar. For more reliable estimations of
the performance of the scintillating detector layers, detailed simulations are
needed. This task will be done in a future study. Finally, in our simulations,
we have considered that the RPC and the scintillating detector planes have a
$100\,\%$ efficiency for the detection of charged particles.

To initiate the simulation, the cores of the CORSIKA air showers were randomly
distributed in the instrumented area of the detector. Then we injected the
shower particles that reach the detector level at the top layer of MATHUSLA.
We followed their trajectories through the tracking planes and recorded the
spatial coordinates of the hit modules and the corresponding hit times, the
latter only in case of one hit per bar. For the Big Pads of the RPC, we also
saved the output signals, which are proportional to the local number of
incident particles. For the simulations of the detector performance, we only
used $e^{\pm}$’s, $\mu^{\pm}$’s, $p^{\pm}$’s, charged pions and kaons. We
allowed unstable particles to decay and removed their decay products from the
simulations. The lifetimes of the decaying particles were taken from PDG . In
addition, we neglected the energy losses of the EAS particles due to
interactions with the detector and the air between the tracking layers. We
assumed that the entire surface of the detector layers of MATHUSLA are active
areas, with no gaps between contiguous scintillating bars or Big Pads.

Figure 11: Example of a simulated vertical EAS impinging on the MATHUSLA
detector. The event was induced by a primary proton with an energy of
$\log_{10}(E/\mathrm{GeV})=5.42$ and arrival direction
$(\theta,\phi)=(11.37^{\circ},34.65^{\circ})$. The azimuth is measured
counterclockwise from the $X$ direction to the horizontal component of the
momentum vector of the primary cosmic ray. The zenith is the angle between the
momentum vector of the incident nucleus and the vertical direction. The core
of the shower is located at coordinates
$(x_{c},y_{c})=(67\,\mbox{m},43\,\mbox{m})$ on the top layer. Top: Central
coordinates of the hit scintillating bars and Big Pads of MATHUSLA during the
event. Each plane corresponds to a different tracking layer. The third layer
of points, from the bottom, are associated to the RPC detector. The axis of
the EAS is shown with a continuous black line. Left bottom: Bi-dimensional
distribution of the induced signal (in arbitrary units) at the Big Pads of the
RPC. Right bottom: Arrival times of the 1st shower particles intersecting the
Big Pads of MATHUSLA RPCs. The timing information is given with respect to the
arrival time of the first charged particle of the EAS that hit the RPC plane.

Figure 12: Example of a simulated inclined EAS impinging on the MATHUSLA
detector. The event was induced by a primary proton with an energy of
$\log_{10}(E/\mathrm{GeV})=7.8$ and arrival direction
$(\theta,\phi)=(70.82^{\circ},139.94^{\circ})$. The azimuth is measured
counterclockwise from the $X$ direction to the horizontal component of the
momentum vector of the primary cosmic ray. The zenith is the angle between the
momentum vector of the incident nucleus and the vertical direction. The core
of the shower is located at the coordinates
$(x_{c},y_{c})=(60\,\mbox{m},30\,\mbox{m})$ on the top layer. Top: Central
coordinates of the scintillating bars and Big Pads that are hit. The points in
each horizontal plane belongs to a different tracking layer. The third layer
of points, from the bottom, are associated to the RPC detector. The axis of
the EAS is shown with a continuous black line. Left bottom: Bi-dimensional
distribution of the induced signal (in arbitrary units) at the Big Pads of the
RPC layer. Right bottom: Arrival times of the 1st shower particles
intersecting the Big Pads of MATHUSLA RPCs. The timing information is given
with respect to the arrival time of the first charged particle of the EAS that
hit the RPC plane.

To illustrate the output of the simulation, in Figs. 11 and 12, we present two
examples: one for a vertical shower and another one for an inclined EAS,
respectively. In these figures, we display the 3D patterns of hits at the
MATHUSLA detector, the corresponding generated signals at the RPC layer, and
the arrival times of the 1st particles hitting each Big Pad of the RPC.

##  6 Event reconstruction

###  6.1 Shower core position

The reconstruction procedure of the EAS events in our analysis goes through
different phases. In the first stage, the shower core position and the arrival
direction of the EAS are estimated. To find the core position a simple
algorithm is applied. For each tracking plane, as a first approximation to the
shower core positions, we used the barycenter of the respective signal
distributions in the $XY$ space, which are calculated using the following
expressions:

$x_{c}=\frac{\sum_{i=1}^{n_{hit}}x_{i}Q_{i}}{\sum_{i=1}^{n_{hit}}Q_{i}},\hskip
28.45274pty_{c}=\frac{\sum_{i=1}^{n_{hit}}y_{i}Q_{i}}{\sum_{i=1}^{n_{hit}}Q_{i}}.$
(1)

Here, $i$ runs over all $n_{hit}$ modules with signal in the corresponding
detector plane, $(x_{i},y_{i})$ is the $XY$ coordinate of the hit $i$-th
detector element and $Q_{i}$ is the corresponding amplitude of the signal of
the Big Pad, where for the scintillating bars, $Q_{i}$ is set to $1$.

We improved the estimates of the EAS core positions by applying a second
algorithm that employs the previous core approximations as initial values. In
this procedure, for each detector layer, we projected the $XY$ signal
distributions on the $X$ and $Y$ axes. The resulting distributions immediately
reveal the position of the core, which is located close to the maximum of the
projected $X$ and $Y$ histograms. We determine the maximums by fitting the X
and Y distributions with lateral functions based on exponential distributions

$N_{\mathrm{hits},x}(x_{i})=\alpha_{1}e^{-\beta_{1}\cdot|x_{i}-x_{c}|},\hskip
28.45274ptN_{\mathrm{hits},y}(y_{i})=\alpha_{2}e^{-\beta_{2}\cdot|y_{i}-y_{c}|},$
(2)

respectively, using a $\chi^{2}$ minimisation procedure. In the above
equations, $\alpha_{1}(\alpha_{2})>0$ and $\beta_{1}(\beta_{2})>0$ are the
amplitude and the slope of the projected $X(Y)$ distribution and, as before,
$(x_{i},y_{i})$ is the coordinate of the $i$-th detector module in a given
detector plane. Finally, $x_{c}$ and $y_{c}$ are the positions of the maximums
of the projected $X$ and $Y$ histograms, which become our improved
approximations to the shower core coordinates. The advantage of this procedure
is that it is fast and can be applied for both high and low-energy EAS.

To reduce the effect of fluctuations during the fit, different bin widths are
used depending on the fraction of detector elements that are hit, $f_{hit}$.
For the RPC histograms of vertical EAS ($\theta\leq 20^{\circ}$) we used bin
widths of $10\,\mbox{m}$ and $2\,\mbox{m}$ for fractions of Big Pads
$f_{hit}\leq 0.2$ and $f_{hit}>0.2$, respectively. Meanwhile, for
$\theta>70^{\circ}$, we used only $2\,\mbox{m}$ bins for the RPC
distributions.

For the projected histograms obtained from the scintillating detector layers a
bin width of $10\,\mbox{m}$ was used for inclined EAS. This bin size was also
employed for vertical EAS that induce a fraction of hit scintillating bars
$f_{hit}\leq 0.02$. However, for larger $f_{hit}$ values associated with
events with the same zenith angle range, the bin width used was $2\,\mbox{m}$
($5\,\mbox{m}$), if the shortest (longest) dimension of the scintillating bars
coincide with the axis of projection used to obtain the histogram.

The bin width selections defined above were chosen to reduce the uncertainties
on the core position. We note that a more detailed study would be needed if
further improvement is desired.

To illustrate the core reconstruction procedure, we present on the left-hand-
side of Figs. 13 and 14 the projected histograms of the signals at the RPC
layer of the MATHUSLA detector that were used in the last step of our method
to find the EAS core of the vertical and inclined MC events presented in Figs.
11 and 12, respectively. The corresponding fits with the lateral distribution
functions of Eqs. (2), which provide the final estimations of the EAS core
positions, are also displayed in the left of Figs. 13 and 14 together with the
best fitting parameters $(x_{c},y_{c})$ for each event.

From the right-hand-side in Figs. 13 and 14, we found that if there are at
least three detector layers where the core position was successfully
reconstructed, a fit with a straight line to the corresponding
$(x_{c},y_{c},z_{c})$ coordinates of the reconstructed shower core can be
applied to get improved positions of the EAS core at each MATHUSLA layer. This
procedure could also be used to find the impact points of the shower core in
the detector planes where this observable was not successfully reconstructed,
provided that the EAS core is inside the instrumented area. This possibility
would permit exploiting the tracking capabilities of MATHUSLA, as illustrated
on the right-hand side of figs. 13 and 14, where we show the fits with a
straight line to the core positions of the events in figs. 11 and 12 estimated
in each MATHUSLA tracking layer. This technique needs to be further
investigated, therefore, it will not be applied in the present analysis.

Figure 13: Left: The projected $X$ (black) and $Y$ (gray) distributions of the
number of Big Pads at the RPC layer of MATHUSLA for the vertical MC event of
Fig. 11. The black and gray lines represent the fits with the formula (2) to
the $X$ and $Y$ projected profiles of the MC shower event, respectively. The
estimated coordinates of the EAS core are shown in the upper left part of the
figure. Right: Estimated coordinates of the impact point of the shower core of
the MC EAS of Fig. 11 at each of the MATHUSLA tracking layers. The open
squares in gray show the respective $(z_{c},y_{c})$ coordinates, and the black
circles, the $(z_{c},x_{c})$ values of the impact points of the EAS core. The
continuous line in black and the dotted line in gray are the corresponding
results of the fits with a linear relation to the $(z_{c},x_{c})$ and
$(z_{c},y_{c})$ data points to improve the estimation of the shower core
position.

Figure 14: Left: The projected $X$ (black) and $Y$ (gray) distributions of the
number of Big Pads at the RPC layer of MATHUSLA for the inclined MC event of
Fig. 12. The black and gray lines represent the fits with the formula (2) to
the $X$ and $Y$ projected profiles of the MC shower event, respectively.
Right: Estimated coordinates of the impact point of the shower core of the MC
EAS of Fig. 12 at each of the MATHUSLA tracking layers. The open squares in
gray show the respective $(z_{c},y_{c})$ coordinates, and the black circles,
the $(z_{c},x_{c})$ values of the impact points of the EAS core. The
continuous line in black and the dotted line in gray are the corresponding
results of the fits with a linear relation to the $(z_{c},x_{c})$ and
$(z_{c},y_{c})$ data points to improve the estimation of the shower core
position.

###  6.2 Arrival direction of the air shower

In a second stage, the arrival direction of the EAS is estimated using the
timing information from the RPC. As a first approximation the direction of the
EAS is found from the fit parameters $(a_{1},a_{2},a_{3})$ of a 3D plane to
the corresponding arrival time data of the shower front. These parameters are
determined by solving the following set of linear equations:

$a_{1}\cdot x_{i}+a_{2}\cdot y_{i}+a_{3}=ct_{1st,i},$ (3)

using the Moore-Penrose pseudoinverse Panos:2007 . In the above equation,
$i=1,...,n_{hit}$ runs over each detector element with signal in the layer,
where $n_{hit}$ is the total number of hit detector modules. The arrival time
of the first EAS particle to the $i$-th detector element is represented by
$t_{1st,i}$ and the respective spatial coordinates are denoted by
$(x_{i},y_{i})$, while $c$ is the velocity of the light in vacuum. The fit
parameters and the arrival direction of the cosmic ray are obtained from the
expressions

$\displaystyle\phi$ $\displaystyle=$ $\displaystyle\arctan{(a_{2}/a_{1})},$
(4) $\displaystyle\theta$ $\displaystyle=$
$\displaystyle\arcsin{\sqrt{a_{1}^{2}+a_{2}^{2}}}.$ (5)

Here, we are using the CORSIKA definitions for the zenith ($\theta$) and the
azimuth ($\phi$) angles of the shower axis Heck:1998vt : $\theta$ is the
measured angle between the momentum vector of the incident nucleus and the
vertical direction, while $\phi$ is defined as the angle between the
horizontal projection of the momentum vector of the primary cosmic ray and the
$X$ axis. The angle $\phi$ is defined as positive for counterclockwise
rotations from the $X$ direction. In this case, we need to add $\pi$ to the
result of equation (4) if $a_{1}<0$ and $a_{2}>0$, and $-\pi$ when both
$a_{1}$ and $a_{2}$ are negative.

To improve the estimation of the arrival direction of the EAS events, we fit
the measured time data with a flat EAS front using the minimum $\chi^{2}$
method and the previous estimations of $\theta$ and $\phi$ as initial values.
We minimised the following expression:

$\chi^{2}=\sum_{j=1}^{m}\left[\frac{t_{1st,j}-T_{0}+(x_{j}-x_{c})n_{x}/c+(y_{j}-y_{c})n_{y}/c}{\sigma_{t}}\right]^{2},$
(6)

where $T_{0}$, $\theta$ and $\phi$ are the fitting parameters. Here, $T_{0}$
is the arrival time of the shower core to the detector layer under
consideration, and $n_{x}$ and $n_{y}$ are the $X$ and $Y$ components of the
normalised vector $\hat{n}$ that is anti-parallel to the momentum of the
cosmic ray,

$\displaystyle\hat{n}$ $\displaystyle=$
$\displaystyle(\sin\theta\cos(\phi-\pi),\,\sin\theta\sin(\phi-\pi),\,\cos\theta),$
(7)

$(x_{c},y_{c})$ is the respective position of the EAS core (estimated with the
procedure described in section 6.1) and $\sigma_{t}$, the estimated
uncertainty in the time measurements, which we conservatively set to
$0.5\,\mbox{ns}$. The index $j$ runs only over the $m$ detector elements that
had a signal within a time window of $10^{5}\,\mbox{ns}$ measured from the
first hit recorded in the layer. For the minimization procedure, we used
MINUIT James:1994vla  from ROOT.

Final approximations of $\theta$ and $\phi$ are obtained by using a cone
instead of a plane in the $\chi^{2}$ fit of the arrival times of the EAS
front, which corrects the timing information for curvature effects. In this
case, we need to minimise the equation below:

$\displaystyle\chi^{2}$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{m^{\prime}}\left[\frac{t_{1st,k}-T_{0}+(x_{k}-x_{c})n_{x}/c+(y_{k}-y_{c})n_{y}/c-r_{k}\tan(\rho)/c}{\sigma_{t}}\right]^{2}.$
(8)

In this relation, $\rho$ is the fit parameter that controls the slope of the
cone with respect to the plane perpendicular to the shower axis, and $r_{k}$
is the radial distance of the $k$-th detector element with signal to the EAS
core in shower disk coordinates. This variable is defined as

$r_{k}=\sqrt{\left[(x_{k}-x_{c})-\Delta_{k}\times
n_{x}\right]^{2}+\left[(y_{k}-y_{c})-\Delta_{k}\times
n_{y}\right]^{2}+\left[\Delta_{k}\times n_{z}\right]^{2}},$ (9)

with

$\Delta_{k}=(x_{k}-x_{c})n_{x}+(y_{k}-y_{c})n_{y}.$ (10)

In equation (8), the sum runs only over detector modules with signals at the
RPC layer recorded during a time window of $100\,\mbox{ns}$ and
$600\,\mbox{ns}$ (taking as a reference the time of the module with the first
hit in the detector layer) for vertical and inclined EAS, respectively. The
above procedure will be also applied for inclined EAS and timing data from the
scintillating detector planes in case of events with charged particle
densities below the saturation threshold of the scintillating bars.

###  6.3 Lateral density distribution of EAS

A useful observable in EAS observatories to study the main properties of
cosmic rays is the lateral distribution of particle densities or the lateral
distribution function (LDF) at the shower front. Therefore, in the next stage
of the EAS reconstruction, the LDF of the event is built using the data from
the RPC. For this task, first, we computed the local density of deposited
charge, $\rho(r)$, in each Big Pad of the RPC detector. Then, we estimated the
radial distances, $r$, to the EAS core at shower disk coordinates (using the
estimated values for the core position and the arrival direction) of the
density measurement and, finally, we calculated the mean values of $\rho(r)$
for different bins of $r$. Here a common bin width of $\Delta r=4\,\mbox{m}$
was employed.

Figure 15: The radial distributions of the mean deposited charge per unit area
at the RPC detector of MATHUSLA (at shower disk coordinates) corresponding to
the vertical (left) and inclined (right) MC events presented in figs. 11 and
12, respectively. The LDFs are shown with data points, which were fitted with
the NKG-like function of Eq. (11) (continuous lines).

The amplitude and the steepness of the radial density distribution of an EAS
are considered to be energy and mass composition sensitive parameters of
primary cosmic rays, respectively Apel2006 ; Bartoli2017b . Therefore, it is
necessary to parameterise the LDFs obtained with the RPC measurements of the
EAS to explore the sensitivity to the primary energy and the mass of the
cosmic-ray nuclei using this information. The parameterization can be achieved
by fitting the LDF data with a Nishimura-Kamata-Greisen (NKG) formula NKG1 ;
NKG2 ; NKG3  often used in this context:

$\rho(r)=A\cdot\left(\frac{r}{r_{0}}\right)^{s-2}\left(1+\frac{r}{r_{0}}\right)^{s-3.5},$
(11)

where $A$ is the amplitude of the LDF, $s$ is a parameter (known as the
lateral shower age) that determines the flatness of the distribution and
$r_{0}$ is a radial scale, whose value we have chosen as the Molière radius
($\sim 79\,\mbox{m}$) at sea level.

Two examples of LDFs obtained from the RPC data for the simulated air showers
presented in figs. 11 and 12 are shown in the left and right panels of Fig.
15, respectively. The radial distributions have been fitted with the NKG
function (11) for $r\leq 50\,\mbox{m}$. The fitted functions are also
displayed in the above figures. As observed from the plots in Fig. 15, the NKG
expression seems to describe reasonably well the LDF of the MC data.

##  7 Selection criteria

For this work, we have selected a set of quality cuts to reduce the magnitude
of the systematic uncertainties in the reconstructed EAS that are used to
study the sensitivity of MATHUSLA to high-energy cosmic rays.

We require that the number of elements with signals in each detector layer
satisfies the condition $n_{hit}\geq 100$ for vertical events, and
$n_{hit}\geq 50$, for inclined EAS. These cuts eliminate low-energy events
with large systematic uncertainties. Note that we have reduced the minimum
number of hits for inclined events in order to lower the energy threshold in
our analysis. In addition, for the analysis of the RPC data, we accepted only
EAS that passed the core and arrival direction reconstruction procedure and
have reconstructed EAS cores within a central area of $90\;\mbox{m}\times 90\
\mbox{m}$ at the RPC layer to avoid border effects. Similar cuts were also
applied for the study of inclined events recorded with the tracking layers of
MATHUSLA. For vertical data from the scintillating layers, we removed the cut
that requires the successful reconstruction of the shower direction at the
respective layer, since we did not estimate the shower direction of vertical
EAS with these detector elements.

##  8 Results

###  8.1 Performance of the RPC

####  8.1.1 Reconstruction efficiency

The selection-cut efficiency of the proposed RPC layer of MATHUSLA to
vertical, $\theta\leq 20^{\circ}$, and inclined,
$\theta=[70^{\circ},80^{\circ}]$, EAS induced by protons and iron nuclei are
presented in Fig. 16. The efficiencies were estimated using our MC simulations
generated with QGSJET-II-04. From these graphs, we observe that MATHUSLA
reaches maximum efficiency at $10^{14}\,\mathrm{eV}$ for vertical air showers,
and increases to $10^{16}\,\mathrm{eV}$ for larger zenith angles, which is
expected due to atmospheric attenuation of the EAS. We did not simulate a
particular trigger condition in our MC simulations. For EAS detection, the RPC
is expected to provide the trigger in MATHUSLA, which will be the same for
vertical and inclined events. A more detailed EAS analysis with a specific
trigger condition for MATHUSLA will be included in a further study.

Figure 16: The selection-cut efficiency of the MATHUSLA RPC layer, estimated
using the QGSJET-II-04 model, for shower events induced by protons (circles)
and iron nuclei (open squares) Left: for vertical EAS with $\theta\leq
20^{\circ}$ and right: for inclined showers with
$\theta=[70^{\circ},80^{\circ}]$.

####  8.1.2 Systematic uncertainties in EAS reconstruction

The performance of the proposed RPC layer as a cosmic-ray detector was
additionally investigated by calculating the systematic uncertainties expected
from the EAS reconstruction described in section 6. First, we have estimated
the bias and resolution of the arrival angle of cosmic-ray EAS as a function
of the primary energy. For this aim, the data is divided in different energy
intervals of width $\Delta\log_{10}(E/\mbox{GeV})=0.2$ and for each bin, the
distributions of the angular systematic uncertainty are calculated. This
systematic is defined as the angle between the reconstructed and true shower
axes of the EAS events. Then, the mean and the $68\%$ containment region of
these distributions are obtained, which are defined as the angular bias
$\Delta\alpha$ and resolution $\Delta\alpha_{68}$ of the RPC layer to EAS
events, respectively, and are plotted versus the central energy of the
corresponding energy bin. The results are shown in Fig. 17 for vertical (left
plot) and inclined events (right graph). The uncertainties were estimated for
H and Fe-induced EAS. The left plot of Fig. 17 shows that both $\Delta\alpha$
and $\Delta\alpha_{68}$ diminish with the primary energy from $E=3.2\times
10^{13}\,\mathrm{eV}$ to $2.5\times 10^{17}\,\mathrm{eV}$. At low energies
($\lesssim 3\times 10^{16}\,\mathrm{eV}$), the angular bias and the resolution
vary between $1^{\circ}$ and $2^{\circ}$, while at high energies, both
uncertainties decrease rapidly to values close to $\sim 0.3^{\circ}$. On the
other hand, the right plot of Fig. 17 shows that for inclined events in the
interval $E=[3.2\times 10^{16}\,\mathrm{eV},2.5\times 10^{17}\,\mathrm{eV}]$,
the angular resolution is not as good. At the RPC layer, the bias
$\Delta\alpha$ is larger for inclined events than for vertical EAS. From MC
simulations, at energies between $3.2\times 10^{16}\,\mathrm{eV}$ and
$2.5\times 10^{17}\,\mathrm{eV}$, we estimated that
$\Delta\alpha=[1.3^{\circ},2.7^{\circ}]$ and
$\Delta\alpha_{68}=[1.3^{\circ},2.9^{\circ}]$ for inclined showers. In this
case, we still observe that both uncertainties decrease at high energies as
for vertical EAS, albeit the reduction is not as pronounced as observed for
the case of vertical showers.

Figure 17: The bias and resolution on the arrival direction of EAS expected
for the MATHUSLA RPC layer as estimated from MC events induced by H and Fe
nuclei. The simulations were produced with QGSJET-II-04. Quality cuts were
applied. Symbols represent the angular bias for protons (circles) and iron
nuclei (open squares). The vertical errors are standard uncertainties of the
mean. Continuous and dashed lines are used to represent the angular
resolutions for H and Fe induced events, respectively. Left: Estimations for
vertical EAS with $\theta\leq 20^{\circ}$. Right: Results for inclined showers
within the interval $\theta=[70^{\circ},80^{\circ}]$.

Figure 18: The bias and resolution on the shower core expected for the
MATHUSLA RPC layer as estimated from MC events induced by H and Fe nuclei. The
simulations were produced with QGSJET-II-04. Quality cuts were applied.
Symbols represent the bias on the shower core location for events produced by
protons (circles) and iron nuclei (open squares). The vertical uncertainty
errors are standard uncertainties of the mean. Continuous and dashed lines are
used to represent the core resolutions for H and Fe-induced events,
respectively. Left: Estimations for vertical EAS with $\theta\leq 20^{\circ}$.
Right: Results for inclined showers within the interval
$\theta=[70^{\circ},80^{\circ}]$.

We also estimated the bias $\Delta R$ and resolution $\Delta R_{68}$ on the
EAS core position expected for the RPC layer at different primary energies.
These quantities are defined as the mean and the $68\%$ containment region of
the distribution for the distance between the reconstructed and the true
shower core locations, respectively, for each energy bin. $\Delta R$ and
$\Delta R_{68}$ are presented in Fig. 18 for vertical and inclined EAS. We
employed a bin width of $\Delta\log_{10}(E/\mbox{GeV})=0.2$ . The plots were
estimated for H and Fe events using QGSJET-II-04. The left plot in Fig. 18
shows that the reconstruction uncertainties for the shower core position of
vertical events decrease from a few tens of meters at $3.2\times
10^{13}\,\mathrm{eV}$ up to a fraction of a meter between $\sim
10^{15}\,\mathrm{eV}$ and $10^{16}\,\mathrm{eV}$, above which the
reconstruction uncertainties remain almost constant with the primary energy.
In particular, the bias (resolution) of the EAS core position at $3.2\times
10^{13}\,\mathrm{eV}$ is $\Delta R=8.6\,\mbox{m}$ ($\Delta
R_{68}=8.6\,\mbox{m}$) and $35.4\,\mbox{m}$ ($42.1\,\mbox{m}$) for H and Fe
events, respectively. At energies larger than $10^{15}\,\mathrm{eV}$, $\Delta
R$ and $\Delta R_{68}$ are about $0.4\,\mbox{m}$ for vertical showers created
by protons and $0.5\,\mbox{m}$ above $10^{16}\,\mathrm{eV}$ for iron nuclei
with $\theta<20^{\circ}$.

On the other hand, from Fig. 18, right, we see that for energies greater than
$3.2\times 10^{15}\,\mathrm{eV}$, the shower-core bias and resolution for
inclined data are of the order of some tens of meters based on our shower
reconstruction method of MC simulations for the MATHUSLA RPC layer. We also
observe that both $\Delta R$ and $\Delta R_{68}$ are reduced at higher
energies, however, their corresponding magnitudes in the high-energy domain do
not reach the same values that we found for vertical events. The reason is the
reduction in the number of hit Big Pads close to the EAS core for inclined
events due to attenuation effects in the atmosphere. For inclined EAS and
energies between $3.2\times 10^{15}\,\mathrm{eV}$ and $2.5\times
10^{17}\,\mathrm{eV}$, $\Delta R$ decreases from $36\,\mbox{m}$ to $\sim
14\,\mbox{m}$, while the core resolution is reduced from $44\,\mbox{m}$ to
approximately $16\,\mbox{m}$.

####  8.1.3 Fraction of hit Big Pads at the RPC layer

We have also investigated the mean fraction of hit Big Pads $f_{hit}$ at the
RPC layer of MATHUSLA that is expected from EAS events of different energies,
see Fig. 19. As before, calculations were performed with our MC simulations
for vertical and inclined showers, and for protons and iron nuclei as
primaries. We used energy bins of size $\Delta\log_{10}(E/\mathrm{GeV})=0.2$.
For the estimations, we computed $f_{hit}$ event-by-event, then, for a given
primary energy bin and cosmic-ray nuclei, we determined the corresponding
average values of $f_{hit}$. From these plots, we can see that $f_{hit}$
increases for vertical showers up to an energy threshold of $\sim 8\times
10^{15}\,\mathrm{eV}$ at which point the mean hit fraction saturates. For
inclined events with $\theta>70^{\circ}$, there is no saturation of $f_{hit}$
at least up to energies of $2.5\times 10^{17}\,\mathrm{eV}$.

Figure 19: The mean fraction of hits in the Big Pads at the MATHUSLA RPC layer
obtained from MC simulations with QGSJET-II-04 for proton (circles) and iron
(open squares) cosmic-ray primaries. The results are presented for two zenith
angle intervals: left column: $\theta=[0^{\circ},20^{\circ}]$ and right
column: $\theta=[70^{\circ},80^{\circ}]$. Standard uncertainties on the mean
were estimated, but they are smaller than the size of the symbols.

According to Fig. 19, below the energy threshold for the saturation of
$f_{hit}$, the expected fraction of hits is different for proton and iron-
induced events. For $\theta<20^{\circ}$, $f_{hit}$ is greater for light
primaries than for heavy nuclei because proton primaries have a larger
production of charged shower particles. For events with $\theta>70^{\circ}$,
where the atmospheric absorption of the electromagnetic component is more
important and the fractional abundance of the muon component increases, the
fraction of hits associated with iron events tends to be slightly larger than
for proton-induced showers.

The fraction of hits increases with the primary energy of the EAS below the
$f_{hit}$ saturation region for both vertical and inclined events. For
inclined EAS, this dependence is linear in a logarithmic scale. These results
suggest that $f_{hit}$ could be used to estimate the cosmic-ray energy below
the threshold for $f_{hit}$ saturation.

####  8.1.4 Maximum density of hit particles at the Big Pads due to EAS events

Individual saturation at the Big Pads of the RPC of MATHUSLA is expected to
appear for more than $\sim 10^{6}$ hits of charged particles per module. These
densities can be found for air showers of very high energies. Therefore, there
could be some EAS energy threshold for which, we would begin to see some
saturation of signals at the Big Pads due to shower events. In order to
investigate this shower energy threshold, we have used our MC simulations to
calculate the mean of the maximum number of hit-charged particles per Big Pad
for EAS as a function of the primary energy. The plots for H and Fe-induced
showers are presented in Fig. 20 for vertical and inclined events. From the
results of Fig. 20, we do not find any potential saturation at the Big Pads
due to hadronic-induced EAS in the energy range of interest.

Figure 20: The mean maximum number of charged particles per Big Pad at the RPC
layer of MATHUSLA induced by EAS events as a function of the primary energy.
The plots were estimated from MC data after applying our selection cuts and
employing the QGSJET-II-04 model. Two primary nuclei were simulated: H
(circles) and Fe (open squares). MC simulations for vertical EAS ($\theta\leq
20^{\circ}$) are presented on the left panel and for inclined showers
($\theta=[70^{\circ},80^{\circ}]$), on the right one. Vertical error bars
correspond to statistical uncertainties of the mean.

####  8.1.5 Spatial structure of EAS at the RPC layer

Figure 21: The 2D charge distributions of the RPC MATHUSLA layer for
individual MC events with $\theta=7^{\circ}$ produced by H and Fe primaries at
three different energies. The simulations were produced using QGSJET-II-04.
The shower cores hit the center of the instrument on the top scintillation
layer. The column on the left was calculated for H primaries, while the column
on the right, for Fe nuclei. The top row contains events with energy
$\log_{10}(E/\mathrm{GeV})=5.0$, the middle row, with
$\log_{10}(E/\mathrm{GeV})=5.5$ and the bottom row, with
$\log_{10}(E/\mathrm{GeV})=6.0$. The bins have dimensions $1\,\mbox{m}\times
1\,\mbox{m}$.

One of the advantages of instrumenting MATHUSLA with an RPC layer is the
possibility of studying the individual morphology of EAS in detail. To
illustrate this point, we present in Fig. 21 some MC examples that show the
bi-dimensional signal distributions at the Big Pads produced by H and Fe
vertical events with shower cores at the centre of the detector for three
different primary energies: $\log_{10}(E/\mathrm{GeV})=5.0$, $5.5$ and $6.0$.
The observed patterns in each of these plots show noticeable differences
depending on the energy and mass of the primary nuclei. For example, the
shower core has a sharper edge, larger amplitudes, and fewer fluctuations at
high energies than at low energies. From Fig. 21, we observe that for energies
in the interval from $10^{5}\,\mathrm{GeV}$ to $10^{6}\,\mathrm{GeV}$, the
proton-induced EAS are more compact and less clumpy than the showers produced
by iron nuclei. The distinctive morphological features of the EAS measured
with the RPC layer are potential observables for distinguishing the nature of
the incident cosmic-ray particle and its primary energy.

Another advantage of the RPC layer is that it clearly shows the position of
the EAS core at very high energies. For showers produced by iron nuclei of low
energies (smaller than a few $10^{5}\,\mathrm{GeV}$), the EAS core is more
difficult to locate from the signal patterns at the RPC because these types of
EAS fronts are clumpy, flatter (as we will see later) and are subject to
strong fluctuations. This is the reason for the rapid increase of the
systematic uncertainty on the shower core position, which we have observed in
Fig. 18 for Fe-induced events in the low-energy regime.

Further details about the physical differences of the EAS fronts can be
obtained with the RPC layer from the corresponding measurements of the lateral
density profiles, which are presented in the next subsection.

####  8.1.6 Lateral distribution of EAS from the RPC data

The radial density profile of EAS as measured with the RPC Big Pads of
MATHUSLA could be a useful observable for studying primary particles on an
event-by-event basis. In Fig. 22, we show the mean radial density profiles,
$\bar{\rho}(r)$, calculated from our RPC simulations for proton and iron
primaries at different energy ranges and for both vertical and inclined
events. The lateral distributions were estimated at shower disk coordinates.
For the calculation, we divided the shower front plane in concentric disks of
width $\Delta r=4\,\mbox{m}$ from $r=0\,\mbox{m}$ to $200\,\mbox{m}$. For each
disk with radius $r_{i}$, we estimated the average, $\bar{\rho}(r_{i})$, of
the density measurements of the Big Pads, whose centers were located inside
the bin $r_{i}$. The local density at each Big Pad was estimated by dividing
the respective signal by the projected area of the module
$A_{BP}(\theta)=1\,\mbox{m}^{2}\times\cos(\theta)$, where $\theta$ is the
zenith angle of the shower axis.

Figure 22: The mean lateral distribution for charged particles of EAS at high
energies expected from our simulations of the MATHUSLA RPC detector. On the
left, for events with zenith angles in the range
$\theta=[0^{\circ},20^{\circ}]$ for three primary energy intervals:
$E=[100\,\mathrm{TeV},1\,\mathrm{PeV}]$, $[1\,\mathrm{PeV},10\,\mathrm{PeV}]$,
$[10\,\mathrm{PeV},100\,\mathrm{PeV}]$, the latter are represented by squares,
triangles, and circles, respectively. On the right, for
$\theta=[70^{\circ},80^{\circ}]$ and $E=[1\,\mathrm{PeV},10\,\mathrm{PeV}]$,
$[10\,\mathrm{PeV},100\,\mathrm{PeV}]$ and
$[100\,\mathrm{PeV},1\,\mathrm{EeV}]$ . The energy intervals are represented
by squares, triangles, and circles, respectively. Solid symbols show the
results for EAS induced by protons, and hollowed points, for iron nuclei.

It is clear from the plots in Fig. 22 that there are differences among the
expected radial density profiles of EAS at the RPC layer depending on the
energy, nature and arrival direction of the incident cosmic ray. For example,
we observe that the amplitude of the distributions increases with the cosmic-
ray energy but decreases with the zenith angle of observation. In addition,
the slopes of the density profiles are flatter for low-energy showers than for
high-energy showers, as well as for vertical events produced by heavy cosmic-
ray nuclei in comparison with vertical EAS due to light mass primaries. From
these results, we expect that the RPC data of the lateral density
distributions of EAS will provide information on the primary energy and the
composition of high-energy cosmic rays. Hence, we should have some shower
observables associated to the radial density profiles that could be employed
as cosmic-ray energy and composition estimators. We have looked for some of
these observables and will present the results in the next subsections.

####  8.1.7 Potential estimators of the cosmic-ray energy from the RPC data

Using our MC simulations, we investigated the amplitude $A$ of the lateral
shower density distributions and the value of the charge density at a fixed
radial distance from the centre of the shower core as potential estimators of
the primary energy of the EAS events. First, we calculated the dependence of
the mean logarithm of the primary energy $\log_{10}(E/\mathrm{GeV})$ on
$\log_{10}(A)$ for H and Fe events and vertical and inclined EAS. The results
are presented in Fig. 23. We used MC events with $f_{hit}>0.03$ for vertical
data and $f_{hit}>0.003$ in the case of inclined showers. This selection was
chosen to reduce the influence of low-energy events in the results, which have
large reconstruction uncertainties.

Figure 23: Mean primary energy as a function of the amplitude of the lateral
distribution $A$ as defined by equation (11) for H (circles) and Fe nuclei
(open squares) simulated with QGSJET-II-04. The plot on the left was built for
EAS with $\theta\leq 20^{\circ}$, while the plot on the right was estimated
for showers within the interval $\theta=[70^{\circ},80^{\circ}]$. Vertical
error bars represent the uncertainty of the mean. In general, the
uncertainties are smaller than the size of the data points.

The plots from Fig. 23 show an approximately linear relation between
$\log_{10}(E/\mathrm{GeV})$ and $\log_{10}(A)$ for both proton and iron
primaries in the energy region of maximum reconstruction efficiency, which
indicates that for vertical and inclined air shower events, the amplitude $A$
can be used to estimate their primary energy. It is important to point out
that, according to Fig. 23, left, a mass dependence of the energy calibration
function based on the $A$ observable for vertical events is expected to
introduce systematic effects when reconstructing the energy spectrum of cosmic
rays with this energy scale. Such dependence, however, decreases at large
zenith angles, see Fig. 23, right.

As an alternative energy estimator, we also explored the value of the lateral
density distribution of the event at a fixed radial distance from the core, in
shower disk coordinates. In this case, as an example, we selected the mean
density at the radial bin $r=18\,\mbox{m}$, which we denoted as
$\rho_{18}=\rho(r=18\,\mbox{m})$. The mean logarithm of the primary energy
versus $\log_{10}(\rho_{18})$, as computed with MC data, is shown in the plots
of Fig. 24 for protons and iron nuclei and within the intervals,
$\theta=[0^{\circ},20^{\circ}]$ and $\theta=[70^{\circ},80^{\circ}]$. We have
applied again the cuts $f_{hit}>0.03$ and $f_{hit}>0.003$ for vertical and
inclined EAS, respectively. In both cases, and independently of the primary
mass composition, we observe from the plots in Fig. 24 that the curves for the
average $\log_{10}(E/\mathrm{GeV})$ vs $\log_{10}(\rho_{18})$ dependence is
approximately linear in the energy region of maximum reconstruction
efficiency, which implies that $\rho_{18}$ could be also employed as a primary
energy estimator. Nevertheless, for vertical EAS, we observe a dependence on
the mass of the primary nuclei. This dependence is also observed for values of
$\rho(r)$ up to distances of at least $r=48\,\mbox{m}$ for
$\theta<20^{\circ}$.

Figure 24: Primary energy as a function of the measured lateral density at the
radial bin centred at $18\,\mbox{m}$ with the MATHUSLA RPC layer for EAS due
to H (circles) and Fe (open squares) nuclei simulated with QGSJET-II-04 and
after applying the quality cuts. The plot on the left was obtained for showers
with $\theta\leq 20^{\circ}$, while the plot on the right was estimated for
events within the interval $\theta=[70^{\circ},80^{\circ}]$. Vertical error
bars represent the standard uncertainty of the mean.

####  8.1.8 Sensitivity of the lateral shower age to the mass of the cosmic
ray

We have previously noted, in section 8.1.5, that the shape of the lateral
distribution in MATHUSLA could be sensitive to the type of the primary
particle. One of the parameters of the lateral distribution of EAS that is
sensitive to the primary composition is the lateral shower age, $s$, as
defined, for example, in equation (11). This parameter measures the flatness
of the lateral distribution of the EAS and it is related to the depth of the
maximum development of the air shower in the atmosphere, $X_{max}$, Apel2006 ;
Bartoli2017a . It is generally known that, at a given energy, light nuclei
create air showers with steeper lateral distributions and larger values of
$X_{max}$ than heavy primaries. Based on such differences, we say that the EAS
associated with the first group are young and the second set, old. On average,
the young EAS have relative $s$ values smaller than the old ones. Thereby the
age parameter has the potential to be used in mass composition studies (see,
for example, HAWCCR23 ).

Figure 25: The mean lateral shower age, $s$, defined in eq. (11), against the
logarithm of the amplitude $A$ of the lateral density distribution of EAS
expected from measurements with the RPC of MATHUSLA. The graphs were estimated
from MC simulations using QGSJET-II-04. The circles and open squares represent
the results for protons and iron primaries, respectively. The plots on the
left are valid for showers with $\theta\leq 20^{\circ}$, while those on the
right, for $\theta=[70^{\circ},80^{\circ}]$. Vertical uncertainties show the
standard deviations of the measurements.

To investigate the sensitivity of the $s$ parameter to the primary cosmic-ray
nature in MATHUSLA, we plot in Fig. 25 the mean age parameter versus the
amplitude $A$ as calculated from our MC simulations for the RPC of the
apparatus. The plots were computed with our MC simulations for vertical and
inclined EAS using both H and Fe primaries. By comparing the graphs for both
nuclei for vertical EAS (c.f. Fig. 25, left), we found that there is a mass
classification capability in the age parameter above some $100\,\mathrm{TeV}$
in MATHUSLA, as the separation of both curves is larger than one standard
deviation in this regime. It is noteworthy that, for vertical events, the mass
sensitivity of the age parameter is, in general, enhanced at high energies, at
least up to $10^{17}\,\mathrm{eV}$, because the uncertainty bars shrink.

For inclined showers, the mass sensitivity of the age parameter in MATHUSLA
decrease as seen from the right plot in Fig. 25, since EAS fluctuations are
larger and the mean shower age values for H and Fe are very similar. This
means that for composition analysis of the primary nuclei with
$\theta>70^{\circ}$, alternative mass-sensitive parameters must be explored.

###  8.2 Performance of the scintillating layers

The analyses shown in the previous sections show that the RPC detector layer
proposed for MATHUSLA could be well suited for studying high-energy cosmic
rays in the experiment. However, the capabilities of the instrument as a
cosmic-ray detector could be enhanced by adding data from the scintillating
layers of MATHUSLA. Before exploring some of these possibilities, we need to
investigate the performance of the scintillating detectors for EAS detection.
For this aim, we carried out additional analyses, which will be presented in
the following subsections.

Figure 26: The selection cut efficiency of the top MATHUSLA scintillating
layer, estimated using the QGSJET-II-04 model, for shower events induced by
protons (black circles) and iron nuclei (open squares) Left: for vertical EAS
with $\theta\leq 20^{\circ}$ and right: for inclined showers with
$\theta=[70^{\circ},80^{\circ}]$. The uncertainty bars show the statistical
uncertainties. The upper and lower limits of the statistical uncertainties
show the values of the $68.2\,\%$ and $31.8\,\%$ containment regions,
respectively, of the efficiency distributions at a given energy, which were
estimated by dividing the MC data into sub-samples according to the energy
interval and calculating the selection cut efficiency for each of them.

####  8.2.1 Reconstruction efficiency

We start the analyses by investigating the energy range of sensitivity of the
scintillating layers of MATHUSLA for cosmic-ray observations. For this task,
using MC simulations, we calculated the selection cut efficiency of the
scintillating planes as a function of the primary energy. Plots for vertical
and inclined events at the top scintillating layer are shown in Fig. 26 to
illustrate the results. By comparing these graphs with those in Fig. 16, we
observe that the energy thresholds for maximum reconstruction efficiency of
EAS at the top scintillating layer are similar to those for the RPC detector
plane. In particular, for vertical EAS, we found that the energy threshold of
the scintillating detector layer is close to $10^{14}\,\mathrm{eV}$, while for
inclined events, it is approximately $10^{16}\,\mathrm{eV}$. Similar
conclusions apply to the other MATHUSLA scintillating detector layers since
their efficiency plots are similar to the ones presented in Fig. 26.

####  8.2.2 Systematic uncertainties in EAS reconstruction

We continue the study of the performance of the MATHUSLA scintillating
detector layers with an estimation of the systematic uncertainties that could
be expected for the reconstruction of EAS events at the scintillating planes
using the procedures described in section 6. To evaluate the systematic
uncertainties, we have calculated the bias $\Delta R$ of the shower core
position for vertical and inclined events and the bias $\Delta\alpha$ in the
arrival direction of inclined showers. We did not calculate $\Delta\alpha$ for
vertical EAS due to saturation problems at some scintillating bars as
explained in section 6. In Fig. 27, we present the shower core bias expected
for the MATHUSLA scintillating layers vs the primary energy. On the other
hand, in Fig. 28, we show the corresponding bias on the arrival direction
$\Delta\alpha$ and the core position $\Delta R$ for inclined EAS. For the
calculation of the bias, we divided the primary energy range in bins of
$\Delta\log_{10}(E/\mathrm{GeV})=0.2$ and followed the procedure described in
section 8.1.2.

Figure 27: The bias on the shower core position of vertical EAS with
$\theta\leq 20^{\circ}$ expected for the MATHUSLA scintillating layers, as
estimated from a selected sample of MC events induced by H (circles) and Fe
(open squares) nuclei. The simulations were produced with the QGSJET-II-04
model. Quality cuts were applied. Vertical uncertainty bars represent the
statistical uncertainties on the mean. They are smaller than the size of the
data points.

Figure 28: The bias on the EAS core position (left) and arrival direction
(right) of inclined showers within the interval
$\theta=[70^{\circ},80^{\circ}]$ expected for the MATHUSLA scintillating
layers. The plots were estimated from a selected sample of MC events using H
(circles) and Fe (open squares) nuclei as primaries. The simulations were
produced with QGSJET-II-04 model. Quality cuts were applied. Vertical
uncertainty bars represent the statistical uncertainties on the mean. They are
smaller than the size of the data points.

In Fig. 27, we note that for vertical events, the bias on the shower core
position $\Delta R$ at the scintillating layers is less than $41\,\mbox{m}$
for H and Fe primaries, respectively, in the region from $3.2\times
10^{13}\,\mathrm{eV}$ to $2.5\times 10^{17}\,\mathrm{eV}$. From the same
figure, we also find that the average $\Delta R$ decreases in the energy
interval from $E=3.2\times 10^{13}\,\mathrm{eV}$ to $3.2\times
10^{14}\,\mathrm{eV}$ ($1.3\times 10^{15}\,\mathrm{eV}$) for proton and iron
nuclei, respectively, and then increases at large energies to $41\,\mbox{m}$.
The reduction in $\Delta R$ at $\mathrm{TeV}$ energies is due to an increment
in the number of hit scintillating bars with the primary energy, while the
increase appears because at very high energies the EAS core covers most of the
surface of the scintillating layers with a large number of tracks (see Fig.
29, left). When all the scintillating layers are covered by tracks it is
difficult to locate with sufficient precision the shower core with our
reconstruction method (see subsection 6.1). It is worth noting that, from Fig.
27, left, the bias $\Delta R$ for EAS produced by iron nuclei is larger than
for proton events at $\mathrm{TeV}$ energies. This difference is related to
the fact that for iron-induced EAS, the shower core is not as well defined as
for proton-induced events with the same primary energy.

For showers with $\theta>70^{\circ}$ (see Fig. 28, left), saturation effects
in the fraction of hit scintillating bars are no longer important for
$E=3.2\times 10^{15}\,\mathrm{eV}-2.5\times 10^{17}\,\mathrm{eV}$. However,
the quality of the reconstruction of the EAS core diminishes mainly because
the shower core is elongated due to the large inclination of the shower axis,
and since the lateral particle density is flatter than for vertical EAS (see
Fig. 29, right). In addition, the shower fluctuations tend to bias the
reconstruction of the EAS core. From Fig. 28, left, we observe that the bias
of the EAS core of inclined events decreases with the primary energy. We also
notice that at $3.2\times 10^{15}\,\mathrm{eV}$, the systematic uncertainty is
equal to $28\,\mbox{m}$ and $35\,\mbox{m}$ for proton and iron events, and
decreases to approximately $17\,\mbox{m}$ at $2.5\times 10^{17}\,\mathrm{eV}$.

Finally, in Fig. 28, right, we notice that the bias $\Delta\alpha$ in the
arrival direction of inclined EAS decreases from $2.9^{\circ}$ at about
$E=3.2\times 10^{15}\,\mathrm{eV}$ up to $1.2^{\circ}$ at $2.5\times
10^{17}\,\mathrm{eV}$ at the scintillating layers. In this case, the magnitude
and energy dependence of $\Delta\alpha$ is similar for inclined EAS produced
by H and Fe primaries.

Figure 29: Mean projections on the $Y$ direction of the number of hit bars at
the second scintillating detector layer of MATHUSLA (counted from the top of
MATHUSLA) expected for proton induced EAS with different energies according to
QGSJET-II-04. The shower cores are at the center of the top MATHUSLA layer.
Left: Vertical events with energies of $100\,\mathrm{TeV}$ (upward triangles),
$320\,\mathrm{TeV}$ (circles), $1\,\mathrm{PeV}$ (squares) and
$10\,\mathrm{PeV}$ (downward triangles). Right: Inclined EAS with energies
$E=10\,\mathrm{PeV}$ (upward triangles), $100\,\mathrm{PeV}$ (circles) and
$320\,\mathrm{PeV}$ (downward triangles). Projections were build for bins in
$y$ of $5\,\mbox{m}$. Vertical uncertainty bars represent uncertainties on the
mean.

###  8.3 Fraction of detector hits per layer

To continue with the study of the performance of the scintillating layers of
MATHUSLA for cosmic ray detection, we also estimated the mean fraction
$f_{hit}$ of scintillating bars activated in each MATHUSLA detector layer by
the EAS at different primary energies using MC data and following the
procedure described in section 8.1.3. We used again energy bins of size
$\Delta\log_{10}(E/\mathrm{GeV})=0.2$. The results for $f_{hit}$ are presented
in Fig. 30. From here, we observe an energy dependence of $f_{hit}$ with the
primary energy at the scintillating layers analogous to that observed at the
RPC system (c.f. Fig. 19). However, for vertical EAS all the detector elements
of the scintillating layers are activated at higher EAS energies ($4\times
10^{16}\,\mathrm{eV}$ ) than for the RPC detector. For inclined EAS, as in the
case of the RPC system, we found that $f_{hit}$ for the scintillating layers
is not saturated below $2.5\times 10^{17}\,\mathrm{eV}$ and that $f_{hit}$
increases with the primary energy following a power-law with a small
sensitivity to the primary mass.

Figure 30: The mean fraction of hit scintillating bars at the MATHUSLA
detector layers obtained from the MC simulations with QGSJET-II-04 for proton
(circles) and iron (open squares) cosmic-ray primaries. The results are
presented for two zenith angle intervals: left:
$\theta=[0^{\circ},20^{\circ}]$ and right: $\theta=[70^{\circ},80^{\circ}]$.
Standard uncertainties on the mean were estimated, but they are smaller than
the size of the symbols.

###  8.4 Saturation of the scintillation detector layers of MATHUSLA

Another key point in the performance of the scintillating detectors of
MATHUSLA studied here is the saturation of the individual scintillating bars.
Saturation at each bar occurs, on average, for particle densities greater than
$\rho_{th}\sim 5\,\mbox{charged part.}/\mbox{m}^{2}$. Therefore, if the
particle density of the EAS exceeds this limit, we will lose information on
the distribution of charged particles at the shower front at the scintillation
detector layers of MATHUSLA.

To investigate the EAS energy range that will be affected in MATHUSLA by the
saturation of the scintillating bars, we have compared the average density
threshold of this effect (after applying a geometrical correction factor
$1/cos(\theta)$ due to the inclination of the shower axis) with MC estimations
for the lateral distributions of charged particles in air showers (at shower
disk coordinates) from primary protons of distinct energies. For vertical
events, we used $E=100\,\mathrm{TeV}$, $320\,\mathrm{TeV}$, $1\,\mathrm{PeV}$
and $10\,\mathrm{PeV}$, while for inclined EAS, we employed
$E=10\,\mathrm{PeV}$, $100\,\mathrm{PeV}$ and $320\,\mathrm{PeV}$. The
results, presented in Fig. 31, left, indicate that the saturation of the
scintillation bars begins to affect the measurements of vertical EAS induced
by protons at energies of $\sim 100\,\mathrm{TeV}$ near the shower core and
that, at higher energies, the area affected by the saturation grows to cover
the whole surface of the detector layers at about $10\,\mathrm{PeV}$. In Fig.
31, right, the saturation of the scintillation bars is not dominant in
MATHUSLA for inclined EAS with $\theta=[70^{\circ},80^{\circ}]$ up to primary
energies of about $3.2\times 10^{17}\,\mathrm{eV}$.

Figure 31: Mean radial densities of charged particles in the EAS front of
events produced by protons of different energies. Left: Vertical events with
energies of $100\,\mathrm{TeV}$ (upward triangles), $320\,\mathrm{TeV}$
(circles), $1\,\mathrm{PeV}$ (squares) and $10\,\mathrm{PeV}$ (downward
triangles). Right: Inclined EAS with energies $E=10\,\mathrm{PeV}$ (upward
triangles), $100\,\mathrm{PeV}$ (circles) and $320\,\mathrm{PeV}$ (downward
triangles). Distributions were calculated in shower disk coordinates. The
densities were estimated for the MATHUSLA site using QGSJET-II-04. The dashed
line shows the value of the average particle density for which the
scintillating bars of MATHUSLA are saturated ($\sim 5\,\mbox{charged
part.}/\mbox{m}^{2}$ for vertical events and $\sim 5\,\mbox{charged
part.}/[\cos(70^{\circ})\times\mbox{m}^{2}]$, for inclined showers). Vertical
uncertainty bars represent statistical uncertainties on the mean. Bins of
$40\,\mbox{m}$ in $r$ were employed to build the plots.

###  8.5 Muon bundles in MATHUSLA

####  8.5.1 Events from highly inclined EAS

At large zenith angles, the relative content of muons in EAS becomes
increasingly larger due to the stronger absorption of electromagnetic and
hadronic particles in the atmosphere. When the arrival direction of the events
has zenith angles greater than $65^{\circ}$, muons become the most abundant
component of EAS EASmodelsReview . This phenomenon produces muon bundles
Bogdanov20 , which are shower events dominated by multiple muon hits at the
detectors. These types of showers are quite interesting to study because muons
are sensitive to the physics of hadronic interactions in EAS, the primary
energy, and the composition of primary cosmic rays EASmodelsReview . At the
location of MATHUSLA, inclined EAS from primary nuclei with $\theta$ between
$70^{\circ}$ and $80^{\circ}$ have an average content of muons of about
$(80\pm 1)\%$ in the region of maximum reconstruction efficiency. This can be
seen in Fig. 32, where we present QGSJET-II-04 predictions of the relative
ratio of muons to charged particles in hadronic EAS for highly inclined events
as a function of the primary energy. Therefore, local measurements with
MATHUSLA of the total shower particle content of such events would be quite
sensitive to the penetrating component of air showers.

Figure 32: Predictions for the average ratio of muons to the total number of
charged particles in hadronic EAS from inclined directions
($\theta=[70^{\circ},80^{\circ}]$) against the primary energy. The MC
simulations were produced for H (circles) and Fe (open squares) primaries
using QGSJET-II-04. Standard uncertainties on the mean were estimated and were
displayed with vertical error bars, however, they can not be seen because they
are smaller than the size of the symbols.

To study the performance of MATHUSLA to these inclined events, we have
estimated the expected spectrum and average value of the local muon density
$D$ for inclined EAS at the RPC layer using the post-LHC models QGSJET-II-04,
EPOS-LHC and SIBYLL 2.3c. The estimations were calculated as a function of the
primary energy for H and Fe nuclei. We define $D$ following Bogdanov20 , as
the ratio of the total number of particle hits per event to the area of the
detector layer. For this study, we carried out a small production of MC events
with EPOS-LHC and SIBYLL 2.3c consisting of $3.3\times 10^{4}$ and $4.1\times
10^{4}$ inclined events, respectively, for each primary nucleus following the
procedure described in section 5. We also weighted the primary energy spectra
of these models to have the same number of events and overall spectral shape
($E^{-2}$) as QGSJET-II-04. The model predictions for the average values and
the spectra of $D$ are shown on the left and right plots of Fig. 33,
respectively. We used bin widths of $\Delta\log_{10}(E/\mathrm{GeV})=0.2$ and
$\Delta\log_{10}(D/\mbox{m}^{2})=0.32$ for the corresponding plots of the
average value of $D$ and its spectrum. From Fig. 33, left, we observe that, at
energies above a few $10^{16}\,\mathrm{eV}$, the local magnitude of $D$
increases linearly with the primary energy in logarithmic scale, and it is
larger for heavy primaries than for light cosmic-ray nuclei. The differences
of $D$ for distinct primaries decrease at low energies due to the loss of
reconstruction efficiency. We note that the sensitivity of $D$ to the hadronic
interaction model is not clear due to the small statistics of the MC data
sets. However, it suggests that QGSJET-II-04 could produce inclined events
with larger values of $D$ at the MATHUSLA RPCs. Whereas the plot on the right-
hand-side of Fig. 33 seems to indicate that the spectra of $D$ for EAS with a
high content of muons are higher and harder than for events with a low
multiplicity of muons according to the results with QGSJET-II-04, EPOS-LHC
Pierog_2015  and SIBYLL 2.3c SIBYLL2.3c_2017 . These types of curves could be
used to test the predictions of high-energy hadronic interaction models with
MATHUSLA RPCs Bogdanov20 . In particular, if the hadronic models predict fewer
muons in EAS above $10^{16}\,\mathrm{eV}$ than the actual measurements as
indicated by some experiments Whisp , then we may expect that the predicted
$D$ spectra, as estimated with a more realistic cosmic-ray composition model,
e.g. Dembinski , are smaller than the measured data.

Figure 33: Left: The mean of the local muon density $D$ expected at the RPC of
MATHUSLA for inclined events against the primary energy of the EAS according
to QGSJET-II-04 (circles), EPOS-LHC (squares) and SIBYLL 2.3c (triangles)
predictions for protons (filled symbols) and iron nuclei (hollowed symbols).
Vertical bars represent standard uncertainty of the mean. Right: Spectra of
the local muon density $D$ of EAS predicted at the RPC of MATHUSLA and
multiplied by a factor $D$. The curves were estimated with the same post-LHC
high-energy hadronic interaction models and primary nuclei as in the plot on
the left. The primary energy spectra of the models were weighted to have the
same number of events and energy spectrum ($\propto E^{-2}$). Selection cuts
were applied to the left and right plots.

Muon bundles induced by cosmic rays at large zenith angles are expected to
have low-particle densities at the MATHUSLA site, with average values smaller
than the threshold for density saturation at the scintillating bars for
primary energies below a few $10^{17}\,\mathrm{eV}$ (see Fig. 31). Under these
circumstances, MATHUSLA could be used not only to perform independent
measurements of local muon densities $D$ but also as a tracking detector to
improve the pointing direction of the instrument. According to Fig. 20,
fluctuations at the shower front are expected to saturate occasionally some
scintillating bars. Despite that, such modules could be removed in a joint EAS
analysis involving the scintillating tracking detector and the RPC detector
layer.

####  8.5.2 Events in tandem with the CMS detector

As a stand-alone cosmic-ray detector, MATHUSLA will sample the front of
hadronic EAS with its tracking system and the RPC detector layer. However,
during the High-Luminosity LHC runs, measurements in combination with the
underground detector CMS may also be carried out. This possibility is still
under exploration with the CMS group. However, if approved, the CMS detector
would make possible EAS measurements of high-energy muons ($E_{\mu}\gtrsim
61\,\mathrm{GeV}$ for vertical incidence222Using the analytical formula for
the range of muons in matter from LipariPRD44 , cosmic-ray muons of
$61\,\mathrm{GeV}$ that propagate vertically in $90\,\mbox{m}$ of standard
rock will arrive to the top of the CMS cavern with $\sim 10\,\mathrm{GeV}$ of
energy, which is enough to traverse the CMS detector CMS_cosmicmuons_2010 .
The value for the muon energy threshold at ground level could be smaller since
the material above the CMS cavern is not completely composed of rock
CMS_cosmicmuons_2010b .) in MATHUSLA at different lateral distances on the
shower plane, ranging from $r=0\,\mbox{m}$ to $168\,\mbox{m}$, depending of
the zenith and azimuth of the shower axis and the impact point of the EAS core
on the surface of MATHUSLA. The data may be taken with independent triggers,
then time correlations could be searched to find EAS events simultaneously
recorded in both instruments. With the CMS data on shower muons, MATHUSLA
would have local, high-precision information about the multiplicity $N_{\mu}$,
radial density $\rho_{\mu}$, charged ratio $\mu^{+}/\mu^{-}$, momentum
$p_{\mu}$ and energy spectrum of muons for large threshold energies in the EAS
front.

For events observed in MATHUSLA and CMS, the most interesting ones would be
those in which the shower cores traverse through the instrumented areas of
both experiments. The muon data from these golden events would be collected
directly from the forward region of the EAS, which is difficult to measure in
dedicated particle physics experiments and is not well constrained by hadronic
interaction models at high energies. In this case. The statistics are expected
to be reduced due to the limited physical area of the CMS detector ($\sim
315\,\mbox{m}^{2}$) and the reduced solid angle of observation for golden
events. If we only consider the field of view subtended by the RPC detector
layer of MATHUSLA from the CMS location, the zenith angle for such events is
restricted from approximately $40^{\circ}$ to $64^{\circ}$ and the width of
the azimuth interval to $\approx 71^{\circ}$ in the CMS reference system
(using the geometry presented in figs. 1 and 3). Hence, we could expect an
acceptance of around $77\,\mbox{m}^{2}\,\mbox{sr}$ in CMS for these events,
which implies a reduction in statistics by a factor of approximately
$2.5\times 10^{-3}$ with respect to MATHUSLA.

Simultaneous EAS measurements with MATHUSLA and CMS also could help to
investigate the connection between hadronic cosmic rays and the properties of
muon bundles in underground experiments. These are air-shower events, which
produce groups of penetrating muons with variable multiplicity in detectors
buried underground under several meters of rock and soil. First measurements
of muon bundles were reported by early deep underground experiments like LVD
LVD_cosmicmuon_90 , Frejus Frejus_cosmicmuon_89  and MACRO MacroBundle .
Detailed studies of muon bundles with underground detectors began with the LEP
experiments ALEPH aleph , L3+C L3C_cosmicmuons_02  and DELPHI delphi . The
results from the LEP experiments showed that the measured rate of muon bundles
with high multiplicity ($N_{\mu}\gtrsim 100$) from cosmic-ray showers was
larger than the predictions of pre-LHC hadronic interaction models cosmoLep
like QGSJET01 qgsjet01 . In recent measurements at the LHC aliceref2 , the
ALICE collaboration has shown that the frequency of muon bundle events with
$\theta<50^{\circ}$ and multiplicities above $100$ are in better agreement
with QGSJET-II-04 expectations if it is assumed that at primary energies
around $10^{17}\,\mathrm{eV}$ the composition of cosmic rays is dominated by
heavy nuclei aliceref1 . Currently, there is no other experiment that confirms
the ALICE result. Data from MATHUSLA and CMS in tandem, if they become
available, would provide an independent check of the ALICE observations.

####  8.5.3 Searching for exotic shower events with MATHUSLA

ALICE results on muon bundles allow for the possibility that these underground
events are due to exotic mechanisms in EAS crCERN . Some theoretical models
propose that the high muon multiplicity events reported by ALICE are produced
by strangelets, hypothetical nuclear matter composed of doublets or triplets
of strange quarks Farhi84 , with possible astrophysical origin Kan2017 . In
this context, the muon bundles detected by ALICE are thought to be created by
a mixture of strange quark matter that loses mass in many successive
interactions with air nuclei when penetrating the atmosphere. According to
Maciej19 , strangelets could induce muon bundles in the atmosphere that could
behave as EAS induced by atomic nuclei with an atomic mass larger than that of
iron at a given energy. Hence, air shower events produced by strangelets would
be characterised by smaller $X_{max}$ values and larger muon contents than
hadronic air showers of the same primary energy. By exploiting the differences
between exotic and hadronic EAS, it could be possible to perform an event-by-
event search for strangelet-induced events in the MATHUSLA cosmic-ray data.

Figure 34: Predictions for the average of the lateral density distributions
measured in the RPCs for muon bundles induced by strangelet events in the
scenario of Maciej19  for primary energies of $5\,\mathrm{PeV}$ (black
diamonds), $45\,\mathrm{PeV}$ (black crosses) and $150\,\mathrm{PeV}$ (black
stars) and zenith angles $\theta<20^{\circ}$. For comparison, the averages of
the lateral density profiles for EAS expected for primary protons and iron
nuclei at different energy intervals according to QGSJET-II-04 are presented.
Results for hadronic EAS for the energy intervals
$100\,\mathrm{TeV}-1\,\mathrm{PeV}$, $1\,\mathrm{PeV}-10\,\mathrm{PeV}$ and
$10\,\mathrm{PeV}-100\,\mathrm{PeV}$ are displayed with gray squares,
triangles and circles. For proton and iron nuclei closed and hollow symbols
are used, respectively. For hadronic events, we also used EAS with zenith
angles below $20^{\circ}$.

To investigate the potential of MATHUSLA to detect possible muon bundles from
an exotic origin in the atmosphere, we simulated the signals of muon bundles
at the MATHUSLA RPCs that are created by strangelets and we estimated the
lateral density distributions of these events. These profiles were compared
with the corresponding predictions for protons and iron-induced EAS in order
to look for possible differences in MATHUSLA between exotic and hadronic
events. To simulate the production of muon bundles due to strangelets, we took
the shower muon content of some QGSJET-II-04 EAS produced by iron nuclei with
vertical incidence ($\theta<20^{\circ}$) and with primary energies
$1\,\mathrm{PeV}$, $10\,\mathrm{PeV}$ and $30\,\mathrm{PeV}$. These shower
muon particles were passed to our MATHUSLA simulation and reconstruction
programs, which are described in sections 5 and 6. Then we applied the
selection cuts described in section 7 and, next, with the chosen data we built
the average of the muon lateral density distributions of the events at
different energies. The magnitudes of the lateral distributions were
multiplied by an appropriate factor to reproduce the muon numbers that we
should expect in strangelet events that have a similar $X_{max}$. To find this
factor we used the predictions for the mean $X_{max}$ and muon content of EAS
induced by protons, iron nuclei and strangelets given in Maciej19 . These
plots indicate that EAS produced by iron nuclei with energies
$1\,\mathrm{PeV}$, $10\,\mathrm{PeV}$ and $30\,\mathrm{PeV}$ have on average
the same $X_{max}$ as muon bundles with energies $5\,\mathrm{PeV}$,
$45\,\mathrm{PeV}$ and $150\,\mathrm{PeV}$, respectively, and that the muon
content of these exotic events has approximately $1.3$ more particles than Fe
induced events with the same primary energy. For QGSJET-II-04, we found
$N_{\mu}\approx 7.3\times 10^{-2}(E/\mathrm{GeV})^{0.92}$ for Fe nuclei with
$\theta<20^{\circ}$ and $E=[1\,\mathrm{PeV},1\,\mbox{EeV}]$ at the MATHUSLA
detector. Hence, we need to apply a scale factor333The procedure has several
drawbacks but it gives us a naive estimation of the effect that we are looking
for in MATHUSLA. In Maciej19 , for example, the EAS were simulated for the
location of the Pierre Auger Observatory, which is situated at an altitude of
$1400\,\mbox{m}$ a.s.l. PAOInclined2014 . In addition, for EAS simulations in
Maciej19 , a modified version of SHOWERSIM showersim84  was employed. Also
$X_{max}$ for hadronic EAS from vertical directions is dominated by photons
and electrons and not by muons EASmodelsReview . of
$1.3\times(E_{str}/E_{Fe})^{0.92}$ to the lateral density of muons simulated
for Fe nuclei of energy $E_{Fe}$ in order to obtain the respective profile for
strangelets of primary energy $E_{str}$.

The average of the lateral density profiles for the simulated muon bundles are
presented in Fig. 34 in comparison with QGSJET-II-04 predictions for vertical
EAS resulting from protons and iron primaries. We note that the lateral
distributions of muon bundles from strangelets decrease less rapidly as a
function of $r$ than the respective distributions for hadronic EAS in the
scenario described in Maciej19 . A more dedicated analysis is underway.
However, this preliminary analysis suggests that MATHUSLA measurements could
test the hypothesis on the strangelet origin of the high multiplicity muon
bundles observed by ALICE by looking for anomalous events with flatter lateral
density distributions than the cosmic-ray-induced EAS. The tests can be more
restrictive if these analyses also include the golden events measured in
correlation with the CMS detector in the hypothetical case that simultaneous
measurements with the underground detector become available. Under these
circumstances, by analysing MATHUSLA data from EAS correlated with muon
bundles detected at the CMS detector, we could directly verify or discard the
exotic nature of these underground events in the context of models like
Maciej19 .

##  9 Discussion

The purpose behind the proposal to install an RPC layer in MATHUSLA is to
exploit the advantages for EAS studies of combining a large area and full
coverage detector that offers a wide dynamical range for charged particle
densities, with the benefits of the cosmic-ray tracking capability of
MATHUSLA. The RPC layer would provide detailed information of cosmic-ray
induced EAS about the spatial and time structure of the shower front in EAS
events with different inclinations, while the tracking detector system of
MATHUSLA would mainly help to measure the arrival times and radial density
profiles of shower particles from inclined EAS. The combination of MATHUSLA
with the RPC detector layer would result in an improved tracking detector for
studying muon bundles from very inclined directions and for measuring vertical
EAS events of TeV energies.

According to the analyses of selection cut efficiency of the RPC layer
presented in subsection 8.1.1 the instrument would be sensitive to EAS from
primary cosmic rays with energies as low as $10^{14}\,\mathrm{eV}$. The
particular energy threshold will depend on the selection cuts and arrival
direction of the shower. For the selection cuts here employed, the proposed
design would have an energy threshold of $10^{14}\,\mathrm{eV}$ to vertically
incident cosmic rays, while for inclined events with zenith angles between
$70^{\circ}$ and $80^{\circ}$ the threshold would be $10^{16}\,\mathrm{eV}$.
For higher energies, EAS detection at the RPC would be restricted by the
reduction in the intensity of cosmic rays and the limited sampling size of the
detector. This point is illustrated in Fig. 35 that shows the expected event
rate of cosmic-ray showers at MATHUSLA as a function of the primary energy
considering only events with shower cores impinging on the instrumented area
and assuming full detector efficiency. From the rate shown on this plot and
considering a collection time of around $3\,\mbox{yr}$, which is similar to
the projected time of the HL-LHC run, we estimate that MATHUSLA would collect
$\sim 530$ air showers with $\theta\leq 70^{\circ}$ at energies
$E>10^{17}\,\mathrm{eV}$ and almost $52$ events with
$\theta=[70^{\circ},80^{\circ}]$ in the same energy range. Therefore, from the
previous results, we could expect that the RPC detector layer would be mainly
constrained to the primary energy range
$E=[10^{14}\,\mathrm{eV},10^{17}\,\mathrm{eV}]$ of the cosmic-ray spectrum,
i.e. the region around the knee. Based on similar analyses for the MATHUSLA
scintillating layers, see section 8.2, we found analogous energy thresholds to
those for the RPC detector.

Figure 35: Number of cosmic-ray events per year above an energy $E$ expected
at MATHUSLA from the zenith angle interval $\theta=[0^{\circ},90^{\circ}]$ and
for an instantaneous field of view of $\pi$ sr considering an horizontal
detector geometry. These results assume that cores of the cosmic-ray induced
EAS are located inside the $100\,\mbox{m}^{2}\times 100\,\mbox{m}^{2}$
instrumented area and that the detector has full efficiency. In the
calculations, the intensities of cosmic rays for different nuclei were taken
from the Global Spline Fit model that was tuned to the latest measurements of
cosmic rays in the energy range from 3 GeV to $10^{11}$ GeV Dembinski . The
black line represents the predictions for the sum of all-particle cosmic-ray
nuclei and the colour lines, the expectations for four representative mass
groups.

From Fig. 18, 17 and 21, we observe that the proposed RPC layer of MATHUSLA
will make possible precise measurements of the core location and arrival
direction of vertical events. The performance of the instrument, however,
would be reduced when increasing the inclination of the shower axes of the
events. In this case, our estimations for the angular resolution of the RPC
layer might be improved by using the tracking capabilities of the MATHUSLA
scintillating detectors. This possibility needs further study, but it should
apply also to vertical events in the $\mathrm{TeV}$ energy regime that have
particle densities below the saturation level of the scintillating bars (see
section 8.4). A good angular resolution is important when building precise
maps of the arrival direction of cosmic rays. These measurements are the basis
of studies of the anisotropy CosmicRays2 ; Deligny2016 ; Ahlers2017  and
clustering of these particles Antoni2004 ; Kang2015  that are done to get a
better understanding of the propagation of this radiation in the interstellar
medium and to constrain models on the spatial distribution of galactic cosmic-
ray sources Pohl_2013 . Shower core measurements could also be possible with
the MATHUSLA scintillating detectors for TeV showers and inclined events (see
figs. 27 and 28). However, we did not observe an improvement here in
comparison with the results for the RPC detector.

From figures 11, 12 and 21, we observe that the proposed RPC detector layer of
MATHUSLA would provide detailed information on the spatial and temporal
structure of the EAS, which could be used to estimate the primary energy of
the event, to measure the distribution of charged particles at the shower
front, and to investigate the primary mass of cosmic rays. This unique
characteristic of the planned RPC detector could help to uncover hints of the
internal hadronic processes that occur inside the EAS. Searches for clusters
or multiple cores in the shower front of an event multicoreArgo2009 ;
multicoreArgo2011 , as well as measurements of the amount and distribution of
charged particles at the forward direction in EAS could be easily done with
the RPC layer and would contribute to testing the validity of modern hadronic
interaction models at high energies qgsjet ; qgsjet2 ; Pierog_2015 ;
SIBYLL2.3c_2017 .

We found that the amplitude of the lateral distribution of the deposited
charge (see Fig. 23) and the value of the charge density at a fixed distance
from the core (for example, Fig. 24) are the measured RPC variables that could
be used to estimate the primary energy of the EAS. The disadvantage here is
that the energy calibration will be model dependent and some influence from
the shower fluctuations may be expected in the result because the measurements
will be made close to ground level. Notwithstanding, cross-calibration with
other cosmic-ray observatories by using the knee feature of the energy
spectrum or the intensity of cosmic rays could be exploited to accomplish this
task, as it is done, for example, in Dembinski . Alternative energy
calibration methods, like neural networks ICECUBE19 , need to be also
investigated. In this case, information on the time structure of the EAS from
the RPC and on the spatial structure of the shower front from the
scintillating planes may increase the precision in the energy assignment of
the shower events.

Another observable that also depends on the primary energy of the EAS is the
fraction of Big Pads with signal, $f_{hit}$, as seen in Fig. 19, left and
right. The sensitivity of this observable, however, seems to be constrained to
energies smaller than $8\times 10^{15}\,\mathrm{eV}$ for vertical showers,
since at higher energies, $f_{hit}$ at the RPC is saturated. In contrast, for
inclined EAS, we do not see saturation effects on $f_{hit}$ at least up to
$2.5\times 10^{17}\,\mathrm{eV}$. In addition, this observable shows a smaller
dependence on the nature of the primary particle for inclined events. Thus,
$f_{hit}$ could also provide an energy estimate with reduced sensitivity to
the mass of the primaries for data with large zenith angles. However, for
vertical events, the mass sensitivity of the fraction hit increases.
Consequently, the fraction hit from the RPC layer may be also used to study
the mass composition of primary nuclei at least in the low-energy regime
($E<10^{15}\,\mathrm{eV}$). From Fig. 30, we notice that the scintillating
layers of MATHUSLA may also provide an estimation of the primary energy of the
events by means of the fraction of scintillating bars with signals, which
exhibits a larger energy threshold for vertical EAS (at around $\sim 4\times
10^{16}\,\mbox{eV}$) than the corresponding one for the RPC detector.

For the research of the composition of cosmic rays with the proposed RPC layer
of MATHUSLA, a mass sensitive observable could be derived from the analysis of
the slope of the lateral distribution of EAS. In this case, the slope of the
lateral density of particles can be parameterised in terms of the lateral
shower age using an NKG-like function (c.f. section 6.3). Based on our studies
in section 8.1.8, the lateral shower age seems to be sensitive to the type of
primary nuclei in MATHUSLA (see Fig. 25) and therefore could be employed for
analyses of the relative abundances of cosmic rays at least for vertical EAS.

Further sensitivity to the mass composition could be gained by looking at the
2D morphology of the air shower, as measured with the RPC detector, by
exploiting the fine-grained spatial resolution of the instrument along with
its large coverage. According to our MC simulations, there appear to be mass-
dependent differences in the 2D shape of the shower front of EAS as measured
with the RPC (see section 8.1.5). Fig. 21 illustrates this point. For example,
the clumpiness in iron-induced EAS seems to be more pronounced than for
protons of the same energy, and proton-induced events seem to produce more
particles in the forward direction than iron nuclei of identical energy. A
neural network or deep learning algorithm could be trained with this detailed
information on the 2D spatial structure of the EAS to improve the mass
classification of cosmic-ray events. In this regard, additional information
from the scintillating layers on the morphology of the EAS front could be
useful. At TeV energies, the MATHUSLA tracking detector could give additional
data on the particle densities at the EAS front in regions outside the shower
core, where the scintillating bars are not saturated (see Fig. 31). At higher
energies ($>1\,\mathrm{PeV}$), however, the measurements of the lateral
particle densities with the tracking detector will be compromised because of
the saturation of both the scintillating bars and the fraction of detector