# Insufficient Gas Accretion Caused the Decline in Cosmic Star-Formation
Activity 8 Billion Years Ago

Aditya Chowdhury  National Centre for Radio Astrophysics, Tata Institute of
Fundamental Research, Pune, India.  Nissim Kanekar  National Centre for Radio
Astrophysics, Tata Institute of Fundamental Research, Pune, India.  Jayaram N.
Chengalur  National Centre for Radio Astrophysics, Tata Institute of
Fundamental Research, Pune, India.

###### Abstract

Measurements of the atomic hydrogen (Hi) properties of high-redshift galaxies
are critical to understanding the decline in the star-formation rate (SFR)
density of the Universe after its peak $\approx 8-11$ Gyr ago. Here, we use
$\approx 510$ hours of observations with the upgraded Giant Metrewave Radio
Telescope to measure the dependence of the average Hi mass of star-forming
galaxies at $z=0.74-1.45$ on their average stellar mass and redshift, by
stacking their Hi 21 cm emission signals. We divide our sample of 11,419 main-
sequence galaxies at $z=0.74-1.45$ into two stellar-mass ($\textrm{M}_{*}$)
subsamples, with $\textrm{M}_{*}>10^{10}~{}\textrm{M}_{\odot}$ and
$\textrm{M}_{*}<10^{10}~{}\textrm{M}_{\odot}$, and obtain clear detections, at
$>4.6\sigma$ significance, of the stacked Hi 21 cm emission in both
subsamples. We find that galaxies with
$\textrm{M}_{*}>10^{10}~{}\textrm{M}_{\odot}$, which dominate the decline in
the cosmic SFR density at $z\lesssim 1$, have Hi reservoirs that can sustain
their SFRs for only a short period, $0.86\pm 0.20$ Gyr, unless their Hi is
replenished via accretion. We also stack the Hi 21 cm emission from galaxies
in two redshift subsamples, at $z=0.74-1.25$ and $z=1.25-1.45$, again
obtaining clear detections of the stacked Hi 21 cm emission signals, at
$>5.2\sigma$ significance in both subsamples. We find that the average Hi mass
of galaxies with $\langle\textrm{M}_{*}\rangle\approx
10^{10}~{}\textrm{M}_{\odot}$ declines steeply over a period of $\approx 1$
billion years, from $(33.6\pm 6.4)\times 10^{9}~{}\textrm{M}_{\odot}$ at
$\langle z\rangle\approx 1.3$ to $(10.6\pm 1.9)\times
10^{9}~{}\textrm{M}_{\odot}$ at $\langle z\rangle\approx 1.0$, i.e. by a
factor $\gtrsim 3$. We thus find direct evidence that accretion of Hi onto
star-forming galaxies at $z\approx 1$ is insufficient to replenish their Hi
reservoirs and sustain their SFRs, thus resulting in the decline in the cosmic
SFR density 8 billion years ago.

Galaxy evolution — Radio spectroscopy — Neutral hydrogen clouds

††software: CASA (McMullin et al., 2007), calR (Chowdhury, 2021), AOFLAGGER
(Offringa et al., 2012), astropy (Astropy Collaboration et al., 2013, 2018)

##  1 Introduction

Understanding the evolution of the star-formation activity in galaxies is a
key open issue in galaxy evolution (e.g. Madau & Dickinson, 2014). The cosmic
comoving star-formation rate (SFR) density of the Universe is known to peak in
the redshift range $z\approx 1-3$, and to then decline by an order of
magnitude over the next $\approx 8$ billion years, to its value today (e.g. Le
Floc’h et al., 2005; Bouwens et al., 2010). Approximately $90\%$ of cosmic
star-formation at $z\lesssim 2.5$ takes place in galaxies lying on the galaxy
“main sequence”, a power-law relation between the SFR and the stellar mass
followed by most star-forming galaxies at any redshift (e.g. Noeske et al.,
2007; Rodighiero et al., 2011). The “amplitude” of the main sequence, i.e. the
SFR at a fixed stellar mass, has been found to decline by a factor of $\approx
10$ between $z\approx 1$ and $z\approx 0$ (e.g. Whitaker et al., 2014),
similar to the decline in the cosmic SFR density over the same period. The
cause of the decline in the star-formation activity of the Universe from
$z\approx 1$ to $z\approx 0$ remains unknown today. Addressing this issue
requires us to understand the evolution of the gas mass of galaxies, the fuel
from which the stars form.

The primary fuel for star-formation is neutral atomic hydrogen (Hi). The
atomic gas content of galaxies is best measured with the hyperfine Hi 21 cm
transition in the Hi ground state, at a wavelength of $\approx 21.11$ cm.
Unfortunately, the low Einstein A-coefficient of this transition has meant
that even deep observations with the best radio telescopes today are only able
to detect Hi 21 cm emission from individual galaxies at low redshifts,
$z\lesssim 0.4$ (e.g. Fernández et al., 2016). However, the average Hi mass of
a sample of galaxies can be measured by “stacking” the Hi 21 cm emission
signals from a large number of galaxies with known spatial positions and
redshifts (e.g. Zwaan, 2000; Chengalur et al., 2001; Lah et al., 2007; Kanekar
et al., 2016; Bera et al., 2019). Such Hi 21 cm-stacking studies have been
used to estimate the average Hi mass of star-forming galaxies out to $z\approx
1.3$ (Chowdhury et al., 2020, 2021).

Recently, Chowdhury et al. (2020) applied the Hi 21 cm-stacking technique to
galaxies in the DEEP2 fields to obtain the first detection of the average Hi
21 cm emission signal from star-forming galaxies at $z\gtrsim 1$ (see also
Chowdhury et al., 2021). Their results suggested that a decline in the Hi
reservoir of star-forming galaxies at $z\gtrsim 1$ could account for the
observed decline in the SFR density of the Universe at later times.

In this _Letter_ , we report first results from the Giant Metrewave Radio
Telescope (GMRT) Cold-Hi AT $z\approx 1$ (CAT$z$1) Survey, a 510-hr upgraded
GMRT Hi 21 cm emission survey of the DEEP2 fields, designed to directly probe
the cause of the decline in the cosmic SFR density, by measuring redshift
evolution in the average Hi mass of star-forming galaxies around $z\approx 1$
when the star-formation activity of the Universe begins to decline
rapidly.111Throughout this work, we use a flat “737” Lambda-cold dark matter
cosmology, with $\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$, and $H_{0}=70$ km
s-1 Mpc-1.

##  2 Observations and Data Analysis

The GMRT-CAT$z$1 survey is a deep Hi 21 cm emission survey of galaxies at
$z=0.74-1.45$ with the upgraded GMRT (Swarup et al., 1991; Gupta et al., 2017)
Band-4 receivers, targetting three fields of the DEEP2 Galaxy Redshift Survey
(Newman et al., 2013). The survey uses the Hi 21 cm stacking approach to carry
out the first characterisation of the Hi properties of galaxies at $z\approx
1$. The DEEP2 Survey is uniquely suited to this experiment because it provides
high-accuracy ($\approx 55$ km s-1; Newman et al., 2013) measurements of the
spectroscopic redshifts of galaxies at $z=0.74-1.45$, whose Hi 21 cm lines are
redshifted into the GMRT Band-4 frequency coverage ($\approx 550-830$ MHz),
and which are located in sky regions that are well-matched to the the GMRT
Band-4 field of view (e.g. Kanekar et al., 2016; Chowdhury et al., 2020,
2021). The survey design, the observations, and the analysis of the GMRT-
CAT$z$1 survey data are described in detail in Chowdhury et al. (2022,
submitted). We provide here a brief summary of the observations and the
analysis procedure.

The survey was carried out over three GMRT observing cycles, using 90 hours in
October 2018 – March 2019 (presented in Chowdhury et al., 2020), 170 hours in
October 2019 – March 2020, and 250 hours in May 2020 – October 2020 (proposals
35_087, 37_063, and 38_033; PI: A. Chowdhury). The total time of $\approx
510$-hours was divided approximately equally between the seven DEEP2 subfields
(21, 22, 31, 32, 33, 41 and 42); the total on-source time on each subfield was
$\approx 50-60$ hours. We used the GMRT Wideband Backend as the correlator,
with a total bandwidth of 400 MHz covering the frequency range $530-930$ MHz,
and divided into 8,192 spectral channels.

The data were analyzed in the Common Astronomy Software Applications package
(casa Version 5; McMullin et al., 2007) following standard procedures (e.g.
Chowdhury et al., 2020) to produce the final continuum images and spectral
cubes for each of the seven DEEP2 subfields. The aoflagger package (Offringa
et al., 2012) was used to excise data affected by radio frequency interference
(RFI). For each DEEP2 subfield, we chose to carry out an independent analysis
of the data from the different GMRT observing cycles. Treating the data of
different cycles independently prevents systematic errors (due to, e.g., low-
level RFI, imperfect deconvolution, etc) in the data of one cycle from
affecting the entire data set.

Our analysis yielded a total of 19 spectral cubes, two each for DEEP2
subfields 21 and 22 (which were each observed in only two GMRT cycles), and
three each for DEEP2 subfields 31, 32, 33, 41, and 42 (which were each
observed in all three GMRT cycles). The spectral cubes have a channel width of
48.8 kHz, corresponding to a velocity resolution of $18-25$ km s-1 across our
frequency coverage. The FWHMs of the synthesized beams of the cubes (i.e. the
angular resolutions of the cubes) are $4.0^{\prime\prime}-7.5^{\prime\prime}$,
corresponding to spatial resolutions of $29-63$ kpc for the redshift range
$z=0.74-1.45$.

The analysis also yielded 655 MHz continuum images of the seven DEEP2
subfields, with Root Mean Square (RMS) noise values of $\approx 5-10~{}\mu$Jy
Beam-1 in the central regions of the images, and synthesized beam FWHMs of
$\approx 3^{\prime\prime}-4^{\prime\prime}$ (Chowdhury et al., 2022,
submitted).

###  2.1 HI 21 cm Emission Stacking

The upgraded GMRT observations covered the redshifted Hi 21 cm line for 16,250
galaxies at $z=0.74-1.45$, for which the Hi 21 cm lines are redshifted to
$580-820$ MHz, and that lie within the full width at half maximum (FWHM) of
the GMRT primary beam at the redshifted Hi 21 cm line frequency of the galaxy.
The FWHM of the GMRT primary beam was assumed to be $\approx 43^{\prime}$ at
610 MHz; the primary-beam FWHM scales $\propto 1/\nu$ with the observing
frequency, $\nu$.

We restricted the Hi 21 cm stacking to blue, star-forming galaxies, and hence
excluded red galaxies and galaxies hosting active galactic nuclei (AGNs) from
the sample. We excluded 2,222 red galaxies by only retaining objects with
${\rm C}\leq 0$, where ${\rm C=U-B}+0.032\times({\rm M_{B}}+21.63)-1.014$
(Willmer et al., 2006; Chowdhury et al., 2020). For AGNs, we excluded 882
DEEP2 objects that were detected at $\geq 4\sigma$ significance in our radio-
continuum images, with 1.4 GHz luminosity L${}_{1.4\ \textrm{GHz}}\geq 2\times
10^{23}$ W Hz-1 (Condon et al., 2002). We also excluded 487 galaxies with
stellar mass $\textrm{M}_{*}<10^{9}\ \textrm{M}_{\odot}$, to ensure that our
results can be directly compared with results for the xGASS survey at
$z\approx 0$ (Catinella et al., 2018). Finally, we carried out a suite of
statistical tests on the Hi 21 cm spectra of individual galaxies, to exclude
any galaxies whose spectra are affected by non-Gaussian systematic effects
(Chowdhury et al., 2020, 2022, submitted).

After the above exclusions, the main sample of the GMRT-CAT$z$1 survey
contains 11,419 blue star-forming galaxies with accurate redshifts (redshift
quality, Q$\geq$3, in the DEEP2 DR4 catalogue; Newman et al., 2013) at
$z=0.74-1.45$, for which the Hi 21 cm line is redshifted to $\approx 580-820$
MHz (i.e. to frequencies where the Band-4 receivers have a high sensitivity),
and that lie within the FWHM of the GMRT primary beam at the redshifted Hi 21
cm frequency of the galaxy. The galaxies have stellar masses
$\textrm{M}_{*}\approx 10^{9}-10^{11}~{}\textrm{M}_{\odot}$.

For each galaxy, we made spectral subcubes, centred on the galaxy’s position
and redshift, and covering $\pm 500$ kpc around the galaxy position, and the
velocity range $\pm 1500$ km s-1 around the galaxy redshift. Each subcube was
resampled onto a grid in the galaxy’s rest frame; the resampled subcubes have
spatial and velocity resolutions of 90 kpc and 30 km s-1, respectively. We
also convolved the subcubes to coarser spatial resolutions ($\leq 200$ kpc)
and find no evidence that the Hi 21 cm emission is spatially resolved at the
resolution of 90 kpc (Chowdhury et al., 2022, submitted).

Our observations yielded multiple Hi 21 cm subcubes, with uncorrelated
statistical noise, for nearly every galaxy of our sample. For each galaxy,
there are $2-3$ independent Hi 21 cm subcubes from the different GMRT cycles.
In addition, some galaxies lie in the overlap regions of our GMRT pointings
(Chowdhury et al., 2022, submitted), yielding additional independent Hi 21 cm
subcubes. Overall, we obtain a total of 28,993 independent Hi 21 cm subcubes
for the 11,419 blue star-forming galaxies of our sample.

For each of the 28,993 Hi 21 cm subcubes, the flux-density values were first
scaled to correct for the primary beam response at the galaxy’s location in
the GMRT primary beam. This was done using an azimuthally-symmetric polynomial
fit to measurements of the average primary beam response of the GMRT Band-4
receivers. Next, we converted each subcube from flux density (${\rm S}$) units
to luminosity density (${\rm L}$) units at the redshift of the DEEP2 galaxy by
using the relation ${\rm L}=4\pi\ {\rm S}\ {\rm D_{L}}^{2}/(1+z)$, where ${\rm
D_{L}}$ is the luminosity distance of the galaxy.

The stacking of the Hi 21 cm emission signals was done by performing a
weighted average of the spectral subcubes, using the same weight for each
spatial and velocity pixel, for different galaxy subsamples, to obtain an
average Hi 21 cm spectral cube for each set of stacked galaxies. The weights
assigned to each subcube were specific to the goal of the stacking, and were
used to control either the redshift or the stellar-mass distributions of the
galaxies in different subsamples. The weights did not take into account the
RMS noise of the individual spectral cubes. The determination of the weights
for the Hi 21 cm subcubes for the specific cases of the stellar-mass
subsamples and the redshift subsamples are discussed later.

Finally, we fitted a second-order polynomial to each spatial pixel of the
stacked Hi 21 cm cube, and then subtracted this out to obtain a residual
spectral cube. The second-order polynomial fit was performed after excluding
the central $\pm 250$ km s-1 velocity range. We note that the exclusion of the
velocity channels in the central $\pm 250$ km s-1 range from the fit results
in a slight increase in the spectral RMS noise of these channels.

The RMS noise on the stacked Hi 21 cm spectral cube was determined via Monte
Carlo simulations. In these simulations, we circularly shifted the central
velocity of each galaxy by a random offset in the range $\pm 1500$ km s-1, and
then stacked the velocity-shifted Hi 21 cm subcubes. The above procedure was
repeated to obtain 104 realizations of the stacked Hi 21 cm subcube. The RMS
noise on every spatial and velocity pixel of the final stacked Hi 21 cm cube
was then estimated from these $10^{4}$ Monte Carlo realizations.

Finally, the average Hi mass of a given subsample of DEEP2 galaxies was
obtained from the velocity-integrated Hi 21 cm line luminosity density
($\int{\rm L_{HI}\ dV}$, in units of ${\rm Jy~{}Mpc^{2}}$ km s-1), measured
from the stacked Hi 21 cm cube, via the expression
$\rm{M_{H{\textsc{i}}}}=[1.86\times 10^{4}\times\int{\rm L_{HI}\
dV}]~{}\textrm{M}_{\odot}$. The integral was carried out after smoothing the
stacked cube to a velocity resolution of 90 km s-1, and was measured from a
contiguous range of spectral channels with $\geq 1.5\sigma$ significance at
this velocity resolution.

###  2.2 Stacking the Rest-frame 1.4 GHz Continuum Emission

The average SFRs of the different galaxy subsamples were measured by stacking
their rest-frame 1.4 GHz radio luminosities (e.g. White et al., 2007; Bera et
al., 2018; Chowdhury et al., 2020), using the correlation between the far-
infrared and radio luminosities (e.g. Condon, 1992; Yun et al., 2001). The
radio continuum stacking was performed by extracting a subimage of each galaxy
from the full continuum image of its DEEP2 subfield, and then smoothing and
regridding each subimage to a spatial resolution of 40 kpc (at the galaxy
redshift) and a uniform grid with 5.2 kpc pixels, extending $\pm 260$ kpc
around each galaxy. We note that the spatial resolution of 40 kpc is much
larger than the typical size, $\approx 8$ kpc, of the star-forming regions of
galaxies at these redshifts (e.g. Trujillo et al., 2004). The flux density
values at the observed frequency in every pixel of each subimage were scaled
to the corresponding rest-frame 1.4 GHz flux densities at the galaxy’s
redshift, assuming a spectral index of $\alpha=-0.8$ (Condon, 1992), with flux
density $S_{\nu}\propto\nu^{\alpha}$. Next, we converted the rest-frame flux-
density (${\rm S_{1.4\ GHz}}$) of every pixel of each subimage to luminosity-
density (${\rm L_{1.4\ GHz}}$) units using the relation ${\rm L_{1.4\
GHz}}=4\pi\ {\rm S_{1.4\ GHz}}\ {\rm D_{L}}^{2}/(1+z)$, where ${\rm D_{L}}$ is
the luminosity distance of the galaxy. The stacked rest-frame 1.4 GHz
luminosity was computed by taking a weighted median of the rest-frame 1.4 GHz
luminosities of the corresponding spatial pixels of the different subimages,
across the subsample of galaxies; the weights were chosen such that each
galaxy has the same effective weight as in the Hi 21 cm stack of the same
subsample. The median rest-frame 1.4 GHz luminosity was then converted to a
median SFR using the relation SFR $(\textrm{M}_{\odot}/\textrm{yr})=(3.7\pm
1.1)\times 10^{-22}\times{\rm L_{1.4\ GHz}(W/Hz)}$ (Yun et al., 2001, scaled
to a Chabrier IMF). Finally, for each stacked image, we corrected the rest-
frame 1.4 GHz luminosity for a small zero-point offset ($\approx 0.4\
\mu$Jy/Beam), at the level of $2-10\%$ of the measured value; the offset was
estimated by applying the above stacking procedure to sub-images at locations
offset by $100^{\prime\prime}$ from the DEEP2 galaxies of the sample.

We note that our radio-derived SFR estimates are affected by uncertainties in
the flux density scale of our radio-continuum images. These systematic
uncertainties are typically $\lesssim 10\%$ for the GMRT, for our calibration
procedure. Thus, all quoted errors on our SFR values have been estimated by
adding, in quadrature, a $10\%$ systematic uncertainty to the $1\sigma$
statistical uncertainty.

##  3 Results and Discussion

Figure 1: [A] The average SFR and the average stellar mass of our sample of
11,419 galaxies. The green and red points show, respectively, the average SFRs
of the subsamples of galaxies at $z=0.74-1.00$ and $z=1.00-1.45$, in stellar
mass bins of width 0.5 dex. The dotted lines show the star-forming main-
sequence of Whitaker et al. (2014), interpolated to the same redshift
intervals. The average stellar mass and the average SFR of galaxies in each
subsample are seen to be consistent with the star-forming main-sequence. [B]
The contribution of star-forming galaxies to the cosmic SFR density (see also
Leslie et al., 2020). The figure shows the cosmic SFR density as a function of
the stellar mass of galaxies in two redshift bins: (i) during the epoch of
peak cosmic SFR density, $z=1.1-1.5$ (blue), and (ii) after the epoch of peak
cosmic SFR density, $z=0.8-1.1$ (orange). Star-forming galaxies with
$\textrm{M}_{*}\geq 10^{10}\ \textrm{M}_{\odot}$ (shown in the shaded region)
dominate the SFR density at both epochs and contribute $\approx 85\%$ to the
decline in the cosmic SFR density at $z\lesssim 1$.

Galaxies lying on the star-forming main-sequence dominate the SFR density of
the Universe out to $z\approx 2.5$ (Rodighiero et al., 2011). The decline in
the SFR density of the Universe at $z\lesssim 1$ is caused by a decline in the
zero-point of the star-forming main sequence (e.g. Whitaker et al., 2014;
Leslie et al., 2020). We compared the stellar masses and SFRs of our sample of
11,419 galaxies with measurements of the star-forming main-sequence at these
redshifts to determine if the sample galaxies are representative of main-
sequence galaxies. This was done by dividing the sample of 11,419 galaxies
into two redshift intervals, $z=0.74-1.00$ and $z=1.00-1.45$, and further
dividing the galaxies in each redshift interval into multiple stellar-mass
subsamples, each of width 0.5 dex. We then used the procedure of Section 2.2
to stack the rest-frame 1.4-GHz continuum emission of the galaxies in each
subsample to determine their average SFRs. For the two redshift intervals, we
used weights during the stacking process to ensure that the redshift
distributions of the stellar-mass subsamples are identical. The average
stellar mass and the average SFR of each subsample, for both redshift
intervals, are shown in Figure 1[A]. The figure also shows, for comparison,
the star-forming main-sequence in the same redshift intervals, from Whitaker
et al. (2014). Figure 1[A] shows that the average SFR and the average stellar
mass of the galaxies in each subsample are consistent with the star-forming
main sequence in both redshift intervals (see also Chowdhury et al., 2022,
submitted). Our sample of 11,419 blue star-forming galaxies at $z=0.74-1.45$
is thus representative of main-sequence galaxies at these redshifts.

Detailed characterizations of the main sequence (Whitaker et al., 2014; Leslie
et al., 2020) and the stellar mass function of high-$z$ galaxies (Davidzon et
al., 2017) have shown that the decline in the cosmic SFR density at $z\lesssim
1$ is predominantly due to the decline in the SFR density of the most massive
galaxies, with stellar masses $\textrm{M}_{*}\gtrsim 10^{10}\
\textrm{M}_{\odot}$. We combined the star-forming main sequence (Whitaker et
al., 2014) and the stellar mass function of galaxies (Davidzon et al., 2017)
to estimate the cosmic SFR density of galaxies as a function of their stellar
mass. Figure 1[B] shows the cosmic SFR density of galaxies as a function of
their stellar mass in two redshift intervals, with mean redshifts of $\langle
z\rangle=1.3$ and $\langle z\rangle=0.95$, respectively. Star-forming galaxies
with $\textrm{M}_{*}\gtrsim 10^{10}\ \textrm{M}_{\odot}$, shown as the shaded
region in Figure 1[B], contribute $\approx 85\%$ to the decline in the total
SFR density from $z\approx 1.3$ to $z\approx 0.95$. Maintaining the Hi
reservoir in such galaxies would require accretion to take place at a high
rate. If the lack of gas accretion in galaxies causes the decline of the
cosmic SFR density at $z\lesssim 1$, one would expect that the most massive
galaxies, which dominate the SFR density at these redshifts, would consume the
bulk of their Hi reservoir on a very short timescale, comparable to the
$\approx 1$-billion year interval between $z\approx 1.3$ and $z\approx 1.0$.

###  3.1 Very Short Hi Depletion Timescale in Massive Galaxies

Figure 2:  The redshift distributions of galaxies in the two stellar-mass
subsamples. The blue histograms show the number of Hi 21 cm subcubes as a
function of redshift in each of the two stellar-mass subsamples; each
histogram has been normalized by the total number of subcubes in its
subsample. The orange histogram in each panel shows the redshift distribution
of galaxies in the high stellar-mass subsample. The Hi 21 cm subcubes of the
lower stellar-mass subsample were assigned weights such that their effective
redshift distribution is the one shown in orange, i.e. is identical to that of
the higher stellar-mass subsample. The number of galaxies in each subsample is
shown above each histogram, with the number of Hi 21 cm subcubes in
parentheses.

We hence examined the dependence of the average Hi mass and the Hi depletion
timescale on the average stellar mass, by dividing the full sample of 11,419
galaxies into two subsamples, with [A] $\textrm{M}_{*}\leq 10^{10}\
\textrm{M}_{\odot}$ (7,834 galaxies and 19,999 independent Hi 21 cm subcubes)
and [B] $\textrm{M}_{*}>10^{10}\ \textrm{M}_{\odot}$ (3,585 galaxies and 8,994
independent subcubes). Figure 2 shows the redshift distributions of the Hi 21
cm subcubes in each of the two stellar-mass subsamples. It is clear that the
lower stellar-mass subsample contains more galaxies at lower redshifts than
the higher stellar-mass subsample; this is because the DEEP2 survey is a
magnitude-limited survey, and is hence biased toward detecting fainter (i.e.,
typically less massive) galaxies at lower redshifts. We account for this
difference in the redshift distributions of the two stellar-mass subsamples by
assigning weights to each galaxy to ensure that the effective redshift
distribution of the two stellar-mass subsamples is identical. Specifically, we
assigned weights to the Hi 21 cm subcubes in the low stellar-mass subsample to
make their effective redshift distribution identical to that of the Hi 21 cm
subcubes in the high stellar-mass subsample. These weights were then used in
stacking both the Hi 21 cm emission and the rest-frame 1.4 GHz radio continuum
of the galaxies of each subsample, following the procedures in Sections 2.1
and 2.2, to determine the average Hi mass and the average SFR of galaxies in
the two subsamples. These weights were further used in the determination of
all average quantities for the two stellar-mass subsamples. In passing, we
note that the effective redshift distribution of the two subsamples has a mean
of $\langle z\rangle\approx 1.1$.

Figure 3:  The stacked Hi 21 cm emission and the stacked rest-frame 1.4 GHz
continuum emission from star-forming galaxies with different stellar masses at
$z\approx 1$. [A] The left panels show, for the two stellar-mass subsamples,
the average Hi 21 cm emission images. The circle on the bottom left of each
panel shows the 90-kpc spatial resolution of the images. The contour levels
are at $-3.0\sigma$ (dashed), $+3.0\sigma$, $+4.0\sigma$, and $+5.0\sigma$
significance. [B] The middle panels show the average Hi 21 cm emission spectra
for the two subsamples; the channel width is 90 km s-1, and the dashed curves
show the $\pm 1\sigma$ error on each spectrum. The average Hi 21 cm emission
signals from both subsamples are clearly detected in both the images and the
spectra. [C] The right panels show the stacked rest-frame 1.4 GHz luminosity
images of the galaxies in the two subsamples. The contour levels are at
$5\sigma,\ 10\sigma,\ 20\sigma,\ 40\sigma,\ {\rm and}\ 80\sigma$ statistical
significance. The circle at the bottom left of each panel shows the 40 kpc
resolution of the images. The average rest-frame 1.4 GHz continuum emission is
detected at high ($>37\sigma$) significance in both images.  Table 1: Average
properties of galaxies in the stellar mass and the redshift subsamples.

|  | Stellar-mass Subsamples |  | Redshift Subsamples
---|---|---|---|---
|  | Low | High |  | Low | High
Stellar Mass Range ($\times 10^{9}\ \textrm{M}_{\odot}$) |  | $1.0-10$ | $10-240$ |  | $1.0-240$ | $1.0-240$
Redshift Range |  | $0.74-1.45$ | $0.74-1.45$ |  | $0.74-1.25$ | $1.25-1.45$
Number of Galaxies |  | $7,834$ | $3,585$ |  | $9,284$ | $2,135$
Number of Hi 21 cm subcubes |  | $19,999$ | $8,994$ |  | 23,630 | 5,363
Average Redshift, $\langle z\rangle$ |  | $1.1$ | $1.1$ |  | 0.98 | 1.33
Average Stellar Mass ($\times 10^{9}\ \textrm{M}_{\odot}$) |  | $4.9$ | $21.4$ |  | 10.3 | 10.3
Average Hi Mass ($\times 10^{9}\ \textrm{M}_{\odot}$) |  | $10.4\pm 2.1$ | $16.3\pm 3.5$ |  | $10.6\pm 1.9$ | $33.6\pm 6.4$
Average SFR ($\textrm{M}_{\odot}\textrm{yr}^{-1}$) |  | $5.5\pm 0.6$ | $18.9\pm 1.9$ |  | $7.5\pm 0.8$ | $15.8\pm 1.7$
Average Hi depletion timescale (Gyr) |  | $1.93\pm 0.44$ | $0.86\pm 0.20$ |  | $1.41\pm 0.31$ | $2.01\pm 0.37$

Note. — For each of the two stellar mass subsamples (columns 2 and 3) and the
two redshift subsamples (columns 4 and 5), the rows are (1) the range of
stellar masses, (2) the number of galaxies in each subsample, (3) the number
of Hi 21 cm subcubes in each subsample, (4) the average redshift of the
galaxies in each subsample, after applying the weights, (5) the average
stellar mass of the galaxies in each subsample, after applying the weights,
(6) the average Hi mass of the galaxies in each subsample, (7) the average SFR
of the galaxies in each subsample, and (8) the average Hi depletion timescale,
$\langle{\rm t}_{\rm
dep,H\textsc{i}}\rangle\equiv\langle\rm{M_{H{\textsc{i}}}}\rangle/\langle{\rm
SFR}\rangle$, of the galaxies in each subsample. See the main text for a
description of the weights that were applied to different subsamples. We note
that the distribution of the galaxies within the GMRT primary beam, as a
function of the distance from the pointing centre, is very similar for all
four subsamples.

Figure 3 shows the stacked Hi 21 cm emission image, the stacked Hi 21 cm
spectrum, and the stacked rest-frame 1.4 GHz continuum image of the galaxies
in each of the two stellar-mass subsamples. We obtain clear detections of the
average Hi 21 cm emission signal from both subsamples, each at $>4.6\sigma$
statistical significance, as well as detections of the rest-frame 1.4 GHz
luminosity, each at $>37\sigma$ significance. This allows us to determine the
average Hi mass and the average SFR of the galaxies of each stellar-mass
subsample.

The Hi depletion timescale,
$\textrm{t}_{\textrm{dep,H}\textsc{i}}=\textrm{M}_{\textrm{H}\textsc{i}}$/SFR,
quantifies the timescale on which a galaxy would convert its entire Hi
reservoir into stars, if continuing to form stars at its current SFR. Star-
formation would be quenched in a galaxy on timescales longer than
$\approx\textrm{t}_{\textrm{dep,H}\textsc{i}}$, if the Hi reservoir is not
replenished by accretion of gas from the circumgalactic medium or by minor
mergers. We combined our estimates of the average Hi mass and the average SFR
in the two stellar-mass subsamples to estimate
$\langle\textrm{t}_{\textrm{dep,H}\textsc{i}}\rangle\equiv\langle\textrm{M}_{\textrm{H}\textsc{i}}\rangle/\langle{\rm
SFR}\rangle$ for each subsample. The inferred average stellar mass, average Hi
mass, average SFR, and average Hi depletion timescale of the two stellar-mass
subsamples are listed in the last four rows of Table 1.

Figure 4 plots the average Hi depletion timescale of the two stellar-mass
subsamples against their average stellar mass. For comparison, we also plot
$\langle\textrm{t}_{\textrm{dep,H}\textsc{i}}\rangle$ against average stellar
mass for a reference sample of $z\approx 0$ galaxies, the stellar mass-
selected xGASS sample (Catinella et al., 2018). For the xGASS galaxies, we
have restricted to blue galaxies, with ${\rm NUV-r}<4$, and used weights to
obtain an effective stellar-mass distribution matched to that of the DEEP2
subsamples. We find that the average Hi depletion timescale of DEEP2 galaxies
with $\textrm{M}_{*}\leq 10^{10}\ \textrm{M}_{\odot}$ is $\approx 2$ Gyr, a
factor of $\approx 3$ smaller than that of similar galaxies in the local
Universe. However, remarkably, we find a far shorter average Hi depletion
timescale, $\langle\textrm{t}_{\textrm{dep,H}\textsc{i}}\rangle=0.86\pm 0.20$
Gyr, in the most massive DEEP2 galaxies, with $\textrm{M}_{*}>10^{10}\
\textrm{M}_{\odot}$, at $z\approx 1.1$. This is comparable to the $\approx
1.0$ Gyr period between $z=1.3$ and $z=1$ (corresponding to lookback times of
$\approx 9$ Gyr and 8 Gyr, respectively), and to the molecular gas depletion
timescale, $\textrm{t}_{\rm dep,H_{2}}\approx 0.7$ Gyr, for similar galaxies
at these redshifts (Tacconi et al., 2013). Star-forming galaxies with
$\textrm{M}_{*}>10^{10}\ \textrm{M}_{\odot}$ would thus be able to sustain
their SFR for only about a billion years unless their Hi is replenished on
this timescale via accretion from the circumgalactic medium or mergers. The
very short Hi depletion timescale in the most massive galaxies, that dominate
the decline in the cosmic SFR density at $z\lesssim 1$, is consistent with the
hypothesis that most of the Hi in such galaxies was consumed by $z\approx 1$,
and that the lack of accretion of Hi onto such galaxies led to the quenching
of their star-formation and hence, to the decline in the cosmic SFR density at
$z\lesssim 1$.

Figure 4: The average Hi depletion timescale, $\langle{\rm t}_{\rm
dep,H\textsc{i}}\rangle\equiv\langle\rm{M_{H{\textsc{i}}}}\rangle/\langle{\rm
SFR}\rangle$, of star-forming galaxies at $z\approx 1.1$, plotted against
average stellar mass. The red circles show our measurements of the average Hi
depletion timescale for the two stellar mass subsamples. The error bars along
the y-axis indicate the $1\sigma$ statistical uncertainties on the
measurements while the error bars along the x-axis show the 90% range of
stellar masses in each subsample. The blue points show the average Hi
depletion timescale, $\langle{\rm t}_{\rm
dep,H\textsc{i}}\rangle\equiv\langle\rm{M_{H{\textsc{i}}}}\rangle/\langle{\rm
SFR}\rangle$, for two subsamples of xGASS galaxies at $z\approx 0$ (Catinella
et al., 2018), each with the same stellar-mass distribution as the
corresponding DEEP2 subsample. The dashed grey line shows the ${\rm H_{2}}$
depletion timescale of star-forming galaxies at $z\approx 1.2$ (Tacconi et
al., 2013). We obtain a very short average Hi depletion timescale, $\approx
0.86$ Gyr, for star-forming galaxies with $\textrm{M}_{*}>10^{10}\
\textrm{M}_{\odot}$ at $z\approx 1.1$, comparable to the ${\rm H_{2}}$
depletion timescale (Tacconi et al., 2013) in similar galaxies at $z\approx
1.2$.

###  3.2 Rapid Redshift Evolution in the Hi mass of Star-forming Galaxies over
$z\approx 0.74-1.45$

The redshift coverage of the GMRT CAT$z$1 survey covers both the peak of
cosmic SFR density at $z\gtrsim 1$ and the decline of the SFR density at
$z\lesssim 1$. This allows us to directly test the hypothesis that
insufficient gas accretion caused the decline in the cosmic SFR density, by
measuring the Hi content of star-forming galaxies of the same stellar mass in
two redshift bins, with average redshifts $\langle z\rangle\approx 1.3$ and
$\langle z\rangle\approx 1$. If the Hi reservoir of star-forming galaxies at
$z\approx 1.3$ is indeed not replenished by accretion, we expect to observe a
significant decline in the average Hi mass of such galaxies, from $\langle
z\rangle=1.3$ to $\langle z\rangle=1$.

Figure 5: The stellar-mass distributions of galaxies in the two redshift
subsamples. The blue histogram shows the number of Hi 21 cm subcubes, in bins
of stellar mass, in the two redshift subsamples; each histogram has been
normalized by the total number of subcubes in its subsample. The orange
histograms show the stellar-mass distribution of galaxies at $z=1.25-1.45$.
The galaxies in the lower-redshift subsample were assigned weights such that
their effective stellar-mass distribution is the one shown in orange, i.e. is
identical to the stellar-mass distribution of galaxies in the higher-redshift
subsample. The number of galaxies in each subsample is shown above the
histogram, with the number of Hi 21 cm subcubes shown in parentheses.

To test this hypothesis, we divided our sample of 11,419 star-forming galaxies
into two subsamples based on the galaxy redshift. The high-redshift sample
includes 2,135 galaxies (5,363 independent Hi 21 cm subcubes) at $z=1.25-1.45$
(i.e. during the epoch of the peak of cosmic SFR density, with an average
redshift $\langle z\rangle\approx 1.3$), while the low-redshift sample
contains 9,284 galaxies (23,630 independent subcubes) at $z=0.74-1.25$ (i.e.
covers the decline of the cosmic SFR density, with $\langle z\rangle\approx
1.0$). The intrinsic stellar-mass distributions of the two subsamples, shown
in Figure 5, are different. We corrected for this effect by assigning weights
to each galaxy in the lower-redshift subsample to ensure that the stellar-mass
distributions of the two subsamples are identical; the weighted-mean stellar
mass of each subsample is $\approx 10^{10}\ \textrm{M}_{\odot}$.

Figure 6: The average Hi 21 cm emission and the average rest-frame 1.4 GHz
continuum emission from star-forming galaxies in the two redshift subsamples.
[A] The left panels show, for the two stellar-mass subsamples, the stacked Hi
21 cm emission images. The circle at the bottom left of each panel shows the
90-kpc spatial resolution of the images. The contour levels are at
$-3.0\sigma$ (dashed), $+3.0\sigma$, $+4.0\sigma$, and $+5.0\sigma$
significance. [B] The middle panels show the average Hi 21 cm emission
spectra; the channel width is 90 km s-1, while the dashed curves show the $\pm
1\sigma$ error on each spectrum. The stacked Hi 21 cm emission signals from
both subsamples are clearly detected in both the images and the spectra. [C]
The right panels show the stacked rest-frame 1.4 GHz luminosity images of the
galaxies in the two subsamples. The contour levels are at $5\sigma,\
10\sigma,\ 20\sigma,{\rm and}\ 40\sigma$ statistical significance. The circle
at the bottom left of each panel shows the 40 kpc resolution of the images.
The average rest-frame 1.4 GHz continuum emission is clearly detected at high
($>32\sigma$) significance in both images.

Using the above weights, we stacked the Hi 21 cm subcubes and the rest-frame
1.4 GHz continuum subimages of the galaxies in the two redshift subsamples,
following the procedures of Section 2.1 and 2.2. This yielded clear detections
of the average Hi 21 cm emission signal (with $>5.2\sigma$ statistical
significance) and the average rest-frame 1.4 GHz luminosity (with $>32\sigma$
statistical significance), from both subsamples. Figure 6 shows the stacked Hi
21 cm emission image, the stacked Hi 21 cm spectrum, and the stacked rest-
frame 1.4 GHz radio-continuum image of the galaxies in the two redshift
subsamples. We find that star-forming galaxies with
$\langle\textrm{M}_{*}\rangle\approx 10^{10}\ \textrm{M}_{\odot}$ at $\langle
z\rangle\approx 1.0$ have an average Hi mass of
$\rm{M_{H{\textsc{i}}}}=(10.6\pm 1.9)\times 10^{9}\ \textrm{M}_{\odot}$, while
those at $\langle z\rangle\approx 1.3$ have an average Hi mass of
$\rm{M_{H{\textsc{i}}}}=(33.6\pm 6.4)\times 10^{9}\ \textrm{M}_{\odot}$. We
thus find evidence, at $3.4\sigma$ statistical significance, that the average
Hi mass of star-forming galaxies with $\langle\textrm{M}_{*}\rangle\approx
10^{10}~{}\textrm{M}_{\odot}$ at $\langle z\rangle\approx 1.3$ is higher, by a
factor of $\approx 3.2$, than that of galaxies with an identical stellar-mass
distribution at $\langle z\rangle\approx 1.0$. Finally, stacking the 1.4 GHz
luminosities of the two subsamples yields average SFRs of $(15.8\pm 1.6)\
\textrm{M}_{\odot}$ yr-1 at $\langle z\rangle\approx 1.3$, but only $(7.5\pm
0.8)\ \textrm{M}_{\odot}$ yr-1 at $\langle z\rangle\approx 1.0$. We thus find
direct evidence that the Hi reservoir of star-forming galaxies is indeed
rapidly depleted between $z\approx 1.3$ and $z\approx 1.0$, i.e. on a
timescale of $\approx 1$ billion years. The decline in Hi mass is by a factor
of $(3.2\pm 0.8)$, consistent with the decline in the average SFR, by a factor
of $(2.1\pm 0.3)$, of the same galaxies over the same redshift interval.

Figure 7 plots the average Hi mass of the DEEP2 galaxies in the two redshift
subsamples against redshift, along with the average Hi mass of blue xGASS
galaxies at $z\approx 0$, again with weights chosen to yield the same stellar-
mass distribution as that of the high-$z$ galaxies. The average Hi mass222The
error on the average Hi mass of the xGASS galaxies was computed using
bootstrap resampling with replacement. of blue xGASS galaxies with
$\langle\textrm{M}_{*}\rangle\approx 10^{10}\ \textrm{M}_{\odot}$ is $(4.1\pm
0.2)\times 10^{9}\ \textrm{M}_{\odot}$. This is a factor of $\approx 2.6$
lower than the average Hi mass of our lower-redshift subsample,
$\rm{M_{H{\textsc{i}}}}=(10.6\pm 1.9)\times 10^{9}\ \textrm{M}_{\odot}$ at
$\langle z\rangle\approx 1.0$. We thus also find evidence (at $\approx
3.4\sigma$ significance) for evolution in the average Hi mass of star-forming
galaxies with the same stellar-mass distribution between $z\approx 1$ and
$z\approx 0$. Overall, Figure 7 shows that the average Hi mass in galaxies
with $\langle\textrm{M}_{*}\rangle\approx 10^{10}\ \textrm{M}_{\odot}$ rapidly
declines, by a factor of $\approx 3.2$, over a period of $\approx 1$ billion
years from $z\approx 1.3$ and $z\approx 1.0$, and then gradually declines, by
a factor of $\approx 2.6$, over the next $\approx 8$ billion years from
$z\approx 1$ to $z\approx 0$.

Finally, we investigated the possibility that our measurement of the decline
in the average Hi mass of star-forming galaxies from $z\approx 1.3$ to
$z\approx 1.0$ might be caused by incompleteness issues in the DEEP2 sample.
First, we note that the galaxies in both redshift subsamples are
representative of the main-sequence population at their redshifts (see Figure
1[A]). We tested for incompleteness issues by restricting our analysis of the
redshift evolution of the Hi mass to only the most massive galaxies of the
sample, with stellar masses greater than the median stellar mass, $6.3\times
10^{9}\ \textrm{M}_{\odot}$; such galaxies are much less likely to be affected
by incompleteness. We divided this massive-galaxy sample of 5,751 galaxies
with $\textrm{M}_{*}>6.3\times 10^{9}~{}\textrm{M}_{\odot}$ into two redshift
subsamples at $z=0.74-1.25$ (4,589 galaxies) and $z=1.25-1.45$ (1,162
galaxies) and separately stacked the Hi 21 cm subcubes of the galaxies in each
redshift subsample, following the procedures of Section 2.1, again ensuring
that the stellar-mass distribution of the two subsamples are identical. We
find that the average Hi masses of star-forming galaxies with
$\langle\textrm{M}_{*}\rangle\approx 1.6\times 10^{10}\ \textrm{M}_{\odot}$ at
$\langle z\rangle\approx 1.0$ and $\langle z\rangle\approx 1.3$ are
$\rm{M_{H{\textsc{i}}}}=(16.0\pm 3.2)\times 10^{9}\ \textrm{M}_{\odot}$ and
$\rm{M_{H{\textsc{i}}}}=(48\pm 10)\times 10^{9}\ \textrm{M}_{\odot}$,
respectively. Thus, even for the massive-galaxy subsample, we continue to find
evidence, at $\approx 3.0\sigma$ significance, for a decline in the average Hi
mass between $z\approx 1.3$ and $z\approx 1.0$. Overall, we do not find any
evidence that the measurement of the redshift evolution of the average Hi mass
of main-sequence galaxies presented in this _Letter_ might be affected by
selection biases.

Figure 7: The redshift evolution of the average Hi mass of star-forming
galaxies over $z=0.74-1.45$. The red squares show our measurements of the
average Hi mass of star-forming galaxies in the two redshift subsamples at
$\langle z\rangle\approx 1.3$ and $\langle z\rangle\approx 1.0$, with
identical stellar mass distributions, obtained by separately stacking the Hi
21 cm signals of galaxies in each subsample. The average Hi mass of blue xGASS
galaxies at $z\approx 0$ (Catinella et al., 2018), with an identical effective
stellar mass distribution, is marked with a blue circle. The average Hi mass
is seen to decline steeply from $z\approx 1.3$ to $z\approx 1.0$, and to then
decline slowly down to $z\approx 0$, for galaxies of the same average stellar
mass, $\langle\textrm{M}_{*}\rangle\approx 10^{10}\ \textrm{M}_{\odot}$.

##  4 Summary

We report first results from the GMRT-CAT$z$1 survey, a $510$-hour Hi 21 cm
emission survey of galaxies at $z=0.74-1.45$ in the DEEP2 fields. We use
stacking of the Hi 21 cm emission and the rest-frame 1.4 GHz radio-continuum
of subsamples of DEEP2 galaxies to determine the dependence of the average Hi
mass and the average Hi depletion time of main-sequence galaxies at
$z=0.74-1.45$ on their average stellar mass and redshift. We find that Hi
plays a fundamental role in the decline of the cosmic SFR density at
$z\lesssim 1$. The very short average Hi depletion timescale of $\approx 0.86$
billion years in the most massive galaxies, which dominate the cosmic SFR
density at $z>1$, implies that rapid gas accretion is needed to maintain their
Hi reservoirs and sustain their high SFRs. And the sharp decline in the
average Hi mass of galaxies with $\langle\textrm{M}_{*}\rangle\approx 10^{10}\
\textrm{M}_{\odot}$, by a factor of $\approx 3.2$ from $z\approx 1.3$ to
$z\approx 1$, indicates that such efficient gas accretion does not take place.
We conclude that the lack of accretion of Hi onto the most massive star-
forming galaxies towards the end of the epoch of peak cosmic SFR density
causes the decline in the cosmic SFR density at $z\lesssim 1$, eight billion
years ago.

We thank the staff of the GMRT who have made these observations possible. The
GMRT is run by the National Centre for Radio Astrophysics of the Tata
Institute of Fundamental Research. We thank an anonymous referee for
suggestions which improved the clarity of the paper. NK acknowledges support
from the Department of Science and Technology via a Swarnajayanti Fellowship
(DST/SJF/PSA-01/2012-13). AC, NK, $\&$ JNC also acknowledge the Department of
Atomic Energy for funding support, under project 12-R&D-TFR-5.02-0700.

## References

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bera et al. (2019) Bera, A., Kanekar, N., Chengalur, J. N., & Bagla, J. S. 2019, ApJ, 882, L7, doi: 10.3847/2041-8213/ab3656
  * Bera et al. (2018) Bera, A., Kanekar, N., Weiner, B. J., Sethi, S., & Dwarakanath, K. S. 2018, ApJ, 865, 39, doi: 10.3847/1538-4357/aad698
  * Bouwens et al. (2010) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2010, ApJ, 709, L133, doi: 10.1088/2041-8205/709/2/L133
  * Catinella et al. (2018) Catinella, B., Saintonge, A., Janowiecki, S., et al. 2018, MNRAS, 476, 875, doi: 10.1093/mnras/sty089
  * Chengalur et al. (2001) Chengalur, J. N., Braun, R., & Wieringa, M. 2001, A&A, 372, 768, doi: 10.1051/0004-6361:20010547
  * Chowdhury (2021) Chowdhury, A. 2021, calR: Robust Calibration in CASA, 3.3.1, Zenodo, doi: 10.5281/zenodo.4723688
  * Chowdhury et al. (2022, submitted) Chowdhury, A., Kanekar, N., & Chengalur, J. N. 2022, submitted, ApJ
  * Chowdhury et al. (2020) Chowdhury, A., Kanekar, N., Chengalur, J. N., Sethi, S., & Dwarakanath, K. S. 2020, Nature, 586, 369, doi: 10.1038/s41586-020-2794-7
  * Chowdhury et al. (2021) Chowdhury, A., Kanekar, N., Das, B., Dwarakanath, K. S., & Sethi, S. 2021, ApJ, 913, L24, doi: 10.3847/2041-8213/abfcc7
  * Condon (1992) Condon, J. J. 1992, ARA&A, 30, 575, doi: 10.1146/annurev.aa.30.090192.003043
  * Condon et al. (2002) Condon, J. J., Cotton, W. D., & Broderick, J. J. 2002, AJ, 124, 675, doi: 10.1086/341650
  * Davidzon et al. (2017) Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70, doi: 10.1051/0004-6361/201730419
  * Fernández et al. (2016) Fernández, X., Gim, H. B., van Gorkom, J. H., et al. 2016, ApJ, 824, L1, doi: 10.3847/2041-8205/824/1/L1
  * Gupta et al. (2017) Gupta, Y., Ajithkumar, B., Kale, H. S., et al. 2017, Current Science, 113, 707
  * Kanekar et al. (2016) Kanekar, N., Sethi, S., & Dwarakanath, K. S. 2016, ApJ, 818, L28, doi: 10.3847/2041-8205/818/2/L28
  * Lah et al. (2007) Lah, P., Chengalur, J. N., Briggs, F. H., et al. 2007, MNRAS, 376, 1357, doi: 10.1111/j.1365-2966.2007.11540.x
  * Le Floc’h et al. (2005) Le Floc’h, E., Papovich, C., Dole, H., et al. 2005, ApJ, 632, 169, doi: 10.1086/432789
  * Leslie et al. (2020) Leslie, S. K., Schinnerer, E., Liu, D., et al. 2020, ApJ, 899, 58, doi: 10.3847/1538-4357/aba044
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, Astronomical Society of the Pacific Conference Series, Vol. 376, CASA Architecture and Applications, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Newman et al. (2013) Newman, J. A., Cooper, M. C., Davis, M., et al. 2013, ApJS, 208, 5, doi: 10.1088/0067-0049/208/1/5
  * Noeske et al. (2007) Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJ, 660, L43, doi: 10.1086/517926
  * Offringa et al. (2012) Offringa, A. R., van de Gronde, J. J., & Roerdink, J. B. T. M. 2012, A&A, 539, A95, doi: 10.1051/0004-6361/201118497
  * Rodighiero et al. (2011) Rodighiero, G., Daddi, E., Baronchelli, I., et al. 2011, ApJ, 739, L40, doi: 10.1088/2041-8205/739/2/L40
  * Swarup et al. (1991) Swarup, G., Ananthakrishnan, S., Kapahi, V. K., et al. 1991, Current Science, 60, 95
  * Tacconi et al. (2013) Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74, doi: 10.1088/0004-637X/768/1/74
  * Trujillo et al. (2004) Trujillo, I., Rudnick, G., Rix, H.-W., et al. 2004, ApJ, 604, 521, doi: 10.1086/382060
  * Whitaker et al. (2014) Whitaker, K. E., Franx, M., Leja, J., et al. 2014, ApJ, 795, 104, doi: 10.1088/0004-637X/795/2/104
  * White et al. (2007) White, R. L., Helfand, D. J., Becker, R. H., Glikman, E., & de Vries, W. 2007, ApJ, 654, 99, doi: 10.1086/507700
  * Willmer et al. (2006) Willmer, C. N. A., Faber, S. M., Koo, D. C., et al. 2006, ApJ, 647, 853, doi: 10.1086/505455
  * Yun et al. (2001) Yun, M. S., Reddy, N. A., & Condon, J. J. 2001, ApJ, 554, 803, doi: 10.1086/323145
  * Zwaan (2000) Zwaan, M. A. 2000, PhD thesis, Ph.D. Thesis, Groningen: Rijksuniversiteit, 2000


#  Discovering Electroweak Interacting Dark Matter at Muon Colliders using
Soft Tracks

Rodolfo Capdevilla  Particle Theory Department, Fermi National Accelerator
Laboratory , Batavia, IL 60510, USA<EMAIL_ADDRESS>Federico Meloni
Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany
<EMAIL_ADDRESS>Jose Zurita  Instituto de Física Corpuscular, CSIC-
Universitat de València, Valencia, Spain<EMAIL_ADDRESS>
###### Abstract

Minimal Dark Matter models feature one neutral particle that serves as a
thermal relic dark matter candidate, as well as quasi-degenerate charged
states with TeV masses. When the charged states are produced at colliders,
they can decay into dark matter and a low-momentum (soft) charged particle,
which is challenging to reconstruct at hadron colliders. We demonstrate that a
3 TeV Muon Collider is capable of detecting these soft tracks, enabling the
discovery of thermal Higgsinos and similar dark matter candidates which
constitute highly motivated scenarios for future collider searches.

††preprint: FERMILAB-PUB-23-0832-T DESY-24-069 APS/123-QED

Introduction - The discovery of the Higgs boson and the multiple measurements
confirming the predictions of the Standard Model (SM) validate the successful
ongoing program at the Large Hadron Collider (LHC). A next step at the
precision frontier is to measure the Higgs properties to disentangle the
microscopic nature of the electroweak symmetry-breaking mechanism. On another
front, the current lack of new particles around the TeV scale motivates the
need for a new boost at the energy frontier. Muon colliders (MuC) are strong
candidates for explorations at the 10 TeV parton center-of-momentum (CoM)
energies. The physics program of future MuC promises both high-precision Higgs
and SM measurements Han  _et al._ (2021a); Ruiz  _et al._ (2022); Han  _et
al._ (2021b); Chen  _et al._ (2022); Forslund and Meade (2022); Chen and Liu
(2024); Ruhdorfer  _et al._ (2023); Forslund and Meade (2024); Liu  _et al._
(2024a); Stylianou and Weiglein (2024); Celada  _et al._ (2023); Cassidy  _et
al._ (2023); Han  _et al._ (2021c, 2022); Garosi  _et al._ (2023) as well as
direct production of new particles with masses well above LHC reach Al Ali
_et al._ (2022); Aime  _et al._ (2022); de Blas  _et al._ (2022); Black  _et
al._ (2024); Accettura  _et al._ (2023); Capdevilla  _et al._ (2021a); Huang
_et al._ (2021); Capdevilla  _et al._ (2022a); Liu and Xie (2021); Han  _et
al._ (2021d); Asadi  _et al._ (2021); Sen  _et al._ (2022); Bandyopadhyay  _et
al._ (2022); Liu  _et al._ (2022a); Casarsa  _et al._ (2022); Capdevilla  _et
al._ (2022b); Bao  _et al._ (2022); İnan and Kisselev (2023a); Liu  _et al._
(2022b); Lv  _et al._ (2022); Kwok  _et al._ (2023); İnan and Kisselev
(2023b); Jueid and Nasri (2023); Jueid  _et al._ (2024); Dermisek  _et al._
(2023); Belfkir  _et al._ (2024); Dermisek  _et al._ (2024); Asadi  _et al._
(2023); Liu  _et al._ (2024b).

Weakly Interacting Massive Particles (WIMPs) are well-motivated dark matter
(DM) candidates that naturally appear in many extensions of the SM. Seeking
these particles is one of the pillars for physics beyond the Standard Model
(BSM) at colliders. In typical WIMP models, the DM candidate is usually
accompanied by other states that can leave distinctive imprints at colliders.
Such is the case of Minimal Dark Matter (MDM) models, where a single
electroweak multiplet is added to the SM Cirelli  _et al._ (2006). Depending
on their quantum numbers, the thermal mass (mass value to explain the observed
relic abundance, $\Omega_{\rm DM}h^{2}=0.1198$ Ade  _et al._ (2016)) of the
MDM multiplets falls in the range of 1-200 TeV Bottaro  _et al._ (2022a, b).

The pure singlet case is phenomenologically ruled out Arkani-Hamed  _et al._
(2006); hence, the new multiplet should be at least a weak doublet. This, in
turn, means that the dark sector must include at least one electrically
charged particle $\chi^{\pm}$ quasi-degenerate with the dark matter candidate
$\chi^{0}$, with mass splittings $\Delta\lesssim$ GeV Thomas and Wells (1998).
We note that the well-studied Higgsino and Wino cases from Minimal
Supersymmetric Standard Model correspond to the Majorana fermion doublet and
triplet, respectively. Furthermore, large multiplet MDM models accounting for
only a fraction ($1-10\%$) of the relic abundance also have TeV scale masses
and are a well-motivated target. Other quasi-degenerate WIMP scenarios include
slepton-like particles, extended scalar sectors, and heavy neutral leptons
Barr and Scoville (2015); Belanger  _et al._ (2019); Jana  _et al._ (2022);
Chiang  _et al._ (2021); Calibbi  _et al._ (2021).

Figure 1: Soft tracks are originated when a heavy particle decays into an
almost mass-degenerate dark matter particle $\chi^{0}$ and a charged particle
$p^{+}$, which can be reconstructed as a $p_{T}\lesssim 1$ GeV track.

The decay process $\chi^{\pm}\to\chi^{0}p^{\pm}$, where $p^{\pm}$ denotes a
low transverse momentum or soft charged particle (Fig. 1) serves as the
primary signature of MDM models. Phenomenological investigations aimed at
discovering these models involve the use of monojet Schwaller and Zurita
(2014); Baer  _et al._ (2014); Agin  _et al._ (2024) or soft-lepton searches
suitable for $\Delta\gtrsim\text{GeV}$ Han  _et al._ (2014); Low and Wang
(2014), and disappearing tracks (DT) targeting $\Delta\lesssim\text{GeV}$
Mahbubani  _et al._ (2017); Fukuda  _et al._ (2018); Bharucha  _et al._
(2018); Calibbi  _et al._ (2018); Filimonova and Westhoff (2019); Belyaev  _et
al._ (2021) and have sparked an extensive search program at the LHC Canepa
_et al._ (2020); ATLAS Collaboration (2013); CMS Collaboration (2015); ATLAS
Collaboration (2018a); CMS Collaboration (2020); ATLAS Collaboration (2022);
Khachatryan  _et al._ (2016); ATLAS Collaboration (2018b); CMS Collaboration
(2018); ATLAS Collaboration (2020); ATLAS Collaboration (2021a); CMS
Collaboration (2022a, b); ATLAS Collaboration (2021b); CMS Collaboration
(2021). Prospects for direct searches for MDM at future colliders have been
investigated at FCC-hh Han  _et al._ (2018); Saito  _et al._ (2019), CLIC de
Blas  _et al._ (2018), and the MuC Han  _et al._ (2021e); Capdevilla  _et al._
(2021b); Bottaro  _et al._ (2021, 2022a, 2022b); Black  _et al._ (2022);
Bandyopadhyay  _et al._ (2024). Although most of these searches primarily
target DT of the charged parent $\chi^{\pm}$, proposals have also been made to
search for the soft charged particle $p^{\pm}$ at electron-proton colliders
Curtin  _et al._ (2018) (as well as the LHC Fukuda  _et al._ (2020); ATLAS
Collaboration (2019)). Furthermore, indirect searches at MuC Franceschini and
Zhao (2023); Fukuda  _et al._ (2024) provide a complementary probe.

Proposed searches for MDM models suggest that both a 100 TeV proton machine,
such as the FCC-hh, and a 10 TeV muon collider (MuC10), will have the
capability to discover the “Wino-like” state with a thermal mass of 2.7 TeV.
The prospect of discovering the “Higgsino-like” state with a thermal mass of
1.1 TeV remains uncertain (some studies suggest that MuC10 may achieve this
feat Han  _et al._ (2021e)). Larger multiplets with thermal masses beyond the
reach of any foreseeable collider, such as the five-plet with masses around
10-15 TeV Bottaro  _et al._ (2022a, b), would necessitate colliders with CoM
energies exceeding 30 TeV Han  _et al._ (2021e), or indirect probes
Franceschini and Zhao (2023); Fukuda  _et al._ (2024).

Figure 2: Signal ($f^{\pm}=\chi^{\pm}$, $f^{0}=\chi^{0}$) and main backgrounds
($f^{\pm}=\ell^{\pm},\tau^{\pm},j$, $f^{0}=\nu_{\ell},\nu_{\tau},j$). The left
diagram is the Drell-Yan-like process, whereas the right blob represents
either Vector Boson Fusion or Bhabha scattering with extra emission of a
vector boson that then splits into a pair of fermions.

In this letter, we demonstrate how a dedicated search for the soft charged
particle $p^{\pm}$ can serve as the discovery channel for MDM models, even at
the early stage of the Muon Collider program with a CoM energy of 3 TeV
(MuC3). Our Soft Track (ST) analysis, indicates that MuC could discover the
Higgsino-like thermal relic. Our analysis utilize MuC detector simulations and
a rigorous (conservative) treatment of the beam-induced background (BIB) from
the decay of the muons as they travel through the machine lattice. As a bonus,
our ST analysis can provide unique insights into the dark sector by probing
the mass splittings between the DM and its companions, potentially
distinguishing MDM from other scenarios. This result highlights the importance
of MuC3 not only as a stage towards MuC10 but also as a powerful discovery
machine.

Setup - We consider a simplified model composed of an electrically charged
particle, $\chi^{\pm}$, and a neutral one, $\chi^{0}$, with a small mass gap,
$\Delta=m_{\chi^{\pm}}-m_{\chi^{0}}\ll m_{\chi^{\pm},\chi^{0}}\equiv
m_{\chi}$. As a proof of concept, in this letter we focus on a search for a
Higgsino-like state at the 3 TeV MuC (hence we will refer to $\chi^{\pm}$ and
$\chi^{0}$ as charginos and neutralinos respectively).

Our focus on this benchmark model is twofold: on one hand, this is one of the
most challenging cases for MDM searches at future colliders. On the other
hand, its thermal mass (1.1 TeV) is low enough for it to be reachable at a
MuC3. To go beyond our benchmark, in our last section we will indicate how to
recast our results to other models.

The $\chi$ particles are odd under a $Z_{2}$ symmetry that renders $\chi^{0}$
stable (hence a good dark matter candidate) and forces these particles to be
pair-produced at colliders. Fig. 2 illustrates the main production modes for
both signal and background processes in our analysis. Signal production, for
which $f^{\pm}=\chi^{\pm}$ and $f^{0}=\chi^{0}$, is dominated by Drell-Yan-
like process (DY). Background production, for which
$f^{\pm}=\ell^{\pm},\tau^{\pm},j$, where $\ell^{\pm}$ stands for an electron
or muon, and $f^{0}=\nu_{\ell},\nu_{\tau},j^{\prime}$, is dominated by two
processes: Vector Boson Fusion (VBF), where the initial state muons radiate
vector bosons that interact, producing the desired final state, and
$t$-channel Bhabha-like Scattering (BS) with the extra emission of a vector
boson splitting into a pair of fermions.

$\mu^{+}\mu^{-}\to\gamma+{\rm X}\,\,(+\,Z\to\nu\nu)$
---
X |  $\,\,\,\,\sigma(\gamma{\rm X})\,$ [fb] |  $\,\,\,\,\sigma(\gamma{\rm X}Z)\,$ [fb]
$\ell^{+}\ell^{-}\nu_{\ell}\bar{\nu}_{\ell}$ | 242.0 | 2.828
$\ell^{+}\ell^{-}\mu^{+}\mu^{-}$ | 60.45 | 0.012
$e^{+}\nu_{e}\mu^{-}\bar{\nu}_{\mu}$ \+ CP | 226.6 | 2.710
$\tau^{+}\tau^{-}\nu_{\ell}\bar{\nu}_{\ell}$ | 6.493 | 0.058
$\tau^{+}\tau^{-}\mu^{+}\mu^{-}$ | 30.86 | 0.006
$\tau^{+}\nu_{\tau}\mu^{-}\bar{\nu}_{\mu}$ \+ CP | 226.2 | 2.722
$jj\nu_{\ell}\bar{\nu}_{\ell}$ | 104.5 | 0.904
$jj\mu^{+}\mu^{-}$ | 30.63 | 0.019
$jj\mu^{-}\bar{\nu}_{\mu}$ \+ CP | 1215. | 11.57
Table 1: Background processes for the 3 TeV MuC. The third column corresponds
to extra emission of a $Z$ boson that decays invisibly. The symbol “+ CP”
indicates that the charge-conjugated process is included. These results were
obtained using the cuts $p_{T}^{\gamma}\geq 20$ GeV, $|\eta_{\gamma}|<2.44$,
$p_{T}^{\ell}\geq 0.1$ GeV, $p_{T}^{j}\geq 0.1$ GeV, $p_{T}^{\ell}\geq 0.1$
GeV and no cuts on $\eta_{\ell}$ and $\eta_{j}$.

Signal - The expected total thermal Higgsino production cross section at
$\sqrt{s}=3$ TeV is 12.53 fb.

The decay $\chi^{\pm}\to\chi^{0}p^{\pm}$ gives rise to a ST and large missing
energy. The distinctive signature we pursue involves a soft charged particle,
produced near the interaction point (IP), which traverses most layers of the
tracker. Objects that can meet these requirements include $\ell,\tau$ or
hadrons. In the case of thermal Higgsinos, the charged state has a 90%
branching ratio (BR) to pions, with only a 10% BR to light leptons. The
characteristic transverse momentum of our ST is approximately
$p_{T}\sim\gamma\Delta$ in the laboratory frame, where $\Delta\sim 0.354$ GeV
and $\gamma$ is the boost of the parent chargino from which the ST originates.
At MuC3, a pair of 1.1 TeV charginos is predominantly produced with a low
boost. This implies a $p_{T}$ spectrum that peaks around few hundreds of MeV
with no events beyond 1 GeV. This challenging signature is the focus of our
study.

Our signal region of interest includes 2 STs with a polar angle
$\theta\in(\pi/3,2\pi/3)$ and at least one hard photon with
$p_{T}^{\gamma}>40$ GeV. The two STs must have opposite electric charge and a
$p_{T}\in(0.2,0.75)$ GeV. The photon cut is pivotal for effective signal-to-
background discrimination. In order to balance the large $p_{T}^{\gamma}$, the
charged particles in the background processes acquire a substantial $p_{T}$.
However, the signal processes are largely unaffected: the impact on the ST
spectrum is minimal because the considerable transverse energy in the event is
balanced by the neutralinos. Finally, we require the ST to have an angular
distance of $\Delta R>1.5$ to help eliminate background from soft photons
splitting into light fermions.

Background - The cross sections of the main background processes are listed in
Table 1. These correspond to $\mu^{+}\mu^{-}\to\gamma$ X, where X is one of
the final states in the left column. The “+ CP” mark indicates that the quoted
cross section includes the charge-conjugated process. Very soft photons
splitting into a pair of light leptons is a particularly large source of
background. We can see this by comparing how the process
$\ell^{+}\ell^{-}\nu_{\ell}\bar{\nu}_{\ell}$ has a much larger cross section
than the analogous processes with taus or quarks. (Our analysis includes
processes with full collinear initial state photons, not shown in Table 1,
simulated with the improved Weizsacker-Williams approximation implemented in
MadGraph Frixione  _et al._ (1993). These backgrounds are easily removed by
our $p_{T}^{\gamma}$ cut.)

For all final states in the left column, we also considered an additional
$Z\to\nu\nu$ emission, which can alter the transverse kinematics and thus
pollute our SR. Emitting this extra particle significantly reduces the cross
section. However, some processes have a cross section comparable to the signal
even after the emission of the extra $Z$ boson. The reported values include
the invisible branching ratio penalty for the $Z$ boson. We note that adding
background processes with an extra invisibly-decaying $Z$ boson adds the
complication of increasing the particle multiplicity in the final state. This
can slow down the calculation beyond a reasonable point. In order to speed up
the calculation of said processes, we use the “mixed beam” (MB) approximation,
where one of the incoming muon beams is replaced by an electron Ruiz  _et al._
(2022). This removes in a gauge invariant manner diagrams with neutral current
that do not correspond to VBF topologies and contribute less than $1\%$ of the
total cross section.

The presence of BIB affects the analysis in two ways: 1) by decreasing the
tracking reconstruction efficiency $\epsilon_{\text{trk}}$ of the ST and 2) by
creating random hit patterns that can be reconstructed by the tracking
algorithm in our detector simulation as extra ST. Our results indicate that
for the most favourable $\theta$ region, the ST reconstruction efficiency
reaches about 45% for $p_{T}\geq 500$ MeV, as shown in Fig. 3. The presence of
BIB causes a decrease of a factor 2–4 in the ST reconstruction efficiency when
compared to a background-free scenario. Fake tracks from BIB hits mainly
populate the forward regions of the tracker. We select tracks with at least 8
hits in the tracking detectors Accettura  _et al._ (2023) and no more than 3
expected-but-missing hits on the fitted trajectory. STs must have a small
transverse ($<3$ mm) and longitudinal ($<15$ mm) impact parameters. These
requirements reduce the rate of fake tracks by more than five orders of
magnitude to a negligible level.

Additionally, BIB energy deposits in the calorimeter could be reconstructed as
fake photons that can be paired with uncorrelated ST from the IP producing
fake signal events. We verified using the full detector simulation described
below that by appropriately raising the calorimeter cell energy thresholds the
photon fake rate for photons $p_{T}\geq 20$ GeV becomes negligible. We expect
future reconstruction techniques to achieve similar rejection power by
exploiting timing information and shower shapes, and have neglected fake
photons in this study.

Figure 3: Track reconstruction efficiency at MuC3 in presence of beam-induced
backgrounds, determined in thermal Higgsino events with a realistic beamspot
distribution Accettura  _et al._ (2023). The efficiency is around 45% for
central tracks with $p_{T}$ between 500 MeV and 1 GeV. The null efficiency in
the highest-$p_{T}$ bins is due to lack of signal events, not reconstruction
effects. A conservative way to extend this map for recasting purposes would be
to extrapolate the last non-null value to higher $p_{T}$.

Finally, additional sources of backgrounds could emerge from beam dynamics.
For example, lepton pairs from incoherent production via two photons can be
generated at the IP. The spectrum of these pairs is very soft ($p_{T}\leq 200$
MeV) so they do not contaminate the SR Calzolari  _et al._ (shed).

Analysis and Detector Simulation - Monte Carlo samples for signal and
backgrounds are simulated, targeting an effective integrated luminosity of 10
ab-1, using MadGraph5_aMC@NLO 3.5.1 Alwall  _et al._ (2014) interfaced to
Pythia 8.307 Sjöstrand  _et al._ (2015) for the parton showering and
hadronisation. The response of the detector, described in Accettura  _et al._
(2023), is modelled using GEANT 4 Agostinelli  _et al._ (2003). Events were
reconstructed and analyzed using the MuonColliderSoft software within key4hep
Sailer  _et al._ (2020); Ganis  _et al._ (2022).

Tracks are reconstructed with the ACTS toolkit Ai  _et al._ (2022) using a
combinatorial Kalman filter algorithm seeded by hit triplets. The tracking
algorithm is configured with the default parameters Accettura  _et al._
(2023), with the exception of the minimum estimated track $p_{T}$ being
lowered to 200 MeV and the width of the $\chi^{2}$ hit search window at each
layer being increased to 30. Furthermore, in order to save computational time,
no hit triplets were build outside of the barrel region of the Vertex
detector, causing the sharp drop in efficiency for $\theta<0.9$ in Fig. 3.

Dedicated particle gun samples including simulated BIB data from the MARS15
package Mokhov and James (2017) for $\sqrt{s}=1.5$ TeV collisions were used to
derive track reconstruction efficiency corrections that are applied to all
reconstructed events when computing the final predictions. The $\sqrt{s}=1.5$
TeV BIB simulation is expected to be a pessimistic approximation based on the
studies available in FLUKA Ferrari  _et al._ (2005); Böhlen  _et al._ (2014)
at $\sqrt{s}=10$ TeV.

Results and Discussion - Fig. 4 shows the $p_{T}$ distribution of the leading
ST after our events selection at MuC3, and compare it with the main
backgrounds. By selecting $p_{T}<750$ MeV we can obtain a sensitivity larger
than $5$ sigma for the thermal dark matter benchmark. An improved pattern
recognition algorithm could substantially improve the sensitivity by
increasing the ST reconstruction efficiency for a similar fake rate. This
result demonstrates that in spite of the BIB, charged soft particles can be
looked for with sufficient signal efficiencies in the central region of the
detector. To our fortune, chargino production at the MuC is dominated by the
DY process, whose angular distribution peaks at the central region, while
background processes (and BIB-related tracks) peak at the forward regions of
the detector.

Figure 4: Signal and background $p_{T}$ distribution of the leading ST at MuC3
for events passing all signal region selections, except that on the leading ST
$p_{T}$. The signal, shown as a red histogram, corresponds to the thermal
Higgsino-like state.

Obtaining the expected sensitivity of our ST analysis for arbitrary models is
straightforward. Since the background distribution is shown in Fig. 4, only
the signal ST $p_{T}$ distribution is required, which is obtainable by
reweighing the $p_{T}$ parton level distribution with the efficiency map shown
in Fig. 3, and applying the selection criteria outlined in our analysis above.
Following this recipe, Fig. 5 shows the expected sensitivity of a Higgsino-
like state as a function of the mass and mean proper lifetime of $\chi^{\pm}$.
For a given mass, increasing the lifetime means decreasing $\Delta$, which
consequently modifies the chargino decay channels. We accounted for all these
variations by utilizing the results in Ibe  _et al._ (2024) as we scanned the
parameter space. Our scan includes lifetimes as low as $1.1\times 10^{-4}$ ns
($\Delta\sim 1.5$ GeV). Below this point, charmed-flavored decay channels
start to appear and the theoretical calculation of the BR suffers from large
hadronic uncertainties Ibe  _et al._ (2024). Our thermal Higgsino-optimized
selection criteria are applied to all masses and lifetimes, causing a loss of
sensitivity for $c\tau\lesssim\sim 8\times 10^{-4}$ ns. A mass-lifetime
dependent optimization, enlarging the covered parameter space, is foreseeable.

For illustration purposes, we present results from DT Capdevilla  _et al._
(2021b) and two signal benchmarks corresponding to triplet (fiveplet) MDM
states with a mass splitting of 166 MeV and a mass of 774 GeV (1.02 TeV),
respectively, whose neutral component account for 10% (1%) of the DM in the
Universe Mitridate  _et al._ (2017). This demonstrates how a combination of DT
and ST analyses can reveal these MDM scenarios.

Figure 5: Reach on the parameter space of MDM models. The blue region will be
covered by the high luminosity stage of the LHC. Solid (dashed) red and green
lines refer to 5$\sigma$ (2$\sigma$) contours from DT and ST analysis at the
MuC3, respectively. The star is the Higgsino thermal target of 1.1 TeV. The
square (triangle) represents a 774 GeV (1.02 TeV) Wino-like (5-plet) state
with a mass splitting of 166 MeV whose neutral component corresponds to 10%
(1%) of the DM in the Universe.

It can be confidently concluded that the soft-track analysis presented in this
study is apt at successfully detecting the thermal Higgsino-like Minimal Dark
Matter particle. Moreover, our technique proves instrumental in characterizing
the dark sector by facilitating the measurement of mass splittings,
potentially enabling the distinction of MDM from other models. It is
noteworthy that our generic approach is versatile and applicable to a range of
other beyond the Standard Model benchmarks, extending beyond the scope of MDM.
The realization that a Muon Collider can uncover highly plausible dark matter
scenarios, even in its early operational phase at an energy of 3 TeV, stands
as a significant milestone for the MuC program.

Acknowledgments - The authors thank Daniel Schulte and Daniele Calzolari for
useful discussions, Priscilla Pani, Andrea Wulzer, Zhen Liu, Simone Pagan
Griso, Donatella Lucchesi, Karol Krizka, and Karri Di Petrillo for their
thorough discussion of the draft, and the International Muon Collider
Collaboration (IMCC) for fostering this effort. The authors are thankful to
Theorist of the Month program at DESY, which enabled the conditions for the
development of this project at its early stage. This manuscript has been
authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359
with the U.S. Department of Energy, Office of High Energy Physics. JZ is
supported by the Generalitat Valenciana (Spain) through the plan GenT program
(CIDEGENT/2019/068), by the Spanish Government (Agencia Estatal de
Investigación) and ERDF funds from European Commission
(MCIN/AEI/10.13039/501100011033, Grant No. PID2020-114473GB-I00). This work
has benefited from computing services provided by the German National Analysis
Facility (NAF).

## References

  * Han  _et al._ (2021a) T. Han, D. Liu, I. Low, and X. Wang, Phys. Rev. D 103, 013002 (2021a), arXiv:2008.12204 [hep-ph] .
  * Ruiz  _et al._ (2022) R. Ruiz, A. Costantini, F. Maltoni, and O. Mattelaer, JHEP 06, 114 (2022), arXiv:2111.02442 [hep-ph] .
  * Han  _et al._ (2021b) T. Han, W. Kilian, N. Kreher, Y. Ma, J. Reuter, T. Striegl, and K. Xie, JHEP 12, 162 (2021b), arXiv:2108.05362 [hep-ph] .
  * Chen  _et al._ (2022) S. Chen, A. Glioti, R. Rattazzi, L. Ricci, and A. Wulzer, JHEP 05, 180 (2022), arXiv:2202.10509 [hep-ph] .
  * Forslund and Meade (2022) M. Forslund and P. Meade, JHEP 08, 185 (2022), arXiv:2203.09425 [hep-ph] .
  * Chen and Liu (2024) M. Chen and D. Liu, Phys. Rev. D 109, 075020 (2024), arXiv:2212.11067 [hep-ph] .
  * Ruhdorfer  _et al._ (2023) M. Ruhdorfer, E. Salvioni, and A. Wulzer, Phys. Rev. D 107, 095038 (2023), arXiv:2303.14202 [hep-ph] .
  * Forslund and Meade (2024) M. Forslund and P. Meade, JHEP 01, 182 (2024), arXiv:2308.02633 [hep-ph] .
  * Liu  _et al._ (2024a) Z. Liu, K.-F. Lyu, I. Mahbub, and L.-T. Wang, Phys. Rev. D 109, 035021 (2024a), arXiv:2308.06323 [hep-ph] .
  * Stylianou and Weiglein (2024) P. Stylianou and G. Weiglein, Eur. Phys. J. C 84, 366 (2024), arXiv:2312.04646 [hep-ph] .
  * Celada  _et al._ (2023) E. Celada, T. Han, W. Kilian, N. Kreher, Y. Ma, F. Maltoni, D. Pagani, J. Reuter, T. Striegl, and K. Xie, (2023), arXiv:2312.13082 [hep-ph] .
  * Cassidy  _et al._ (2023) M. E. Cassidy, Z. Dong, K. Kong, I. M. Lewis, Y. Zhang, and Y.-J. Zheng, (2023), arXiv:2311.07645 [hep-ph] .
  * Han  _et al._ (2021c) T. Han, Y. Ma, and K. Xie, Phys. Rev. D 103, L031301 (2021c), arXiv:2007.14300 [hep-ph] .
  * Han  _et al._ (2022) T. Han, Y. Ma, and K. Xie, JHEP 02, 154 (2022), arXiv:2103.09844 [hep-ph] .
  * Garosi  _et al._ (2023) F. Garosi, D. Marzocca, and S. Trifinopoulos, JHEP 09, 107 (2023), arXiv:2303.16964 [hep-ph] .
  * Al Ali  _et al._ (2022) H. Al Ali _et al._ , Rept. Prog. Phys. 85, 084201 (2022), arXiv:2103.14043 [hep-ph] .
  * Aime  _et al._ (2022) C. Aime _et al._ , (2022), arXiv:2203.07256 [hep-ph] .
  * de Blas  _et al._ (2022) J. de Blas _et al._ , (2022), arXiv:2203.07261 [hep-ph] .
  * Black  _et al._ (2024) K. M. Black _et al._ , JINST 19, T02015 (2024), arXiv:2209.01318 [hep-ex] .
  * Accettura  _et al._ (2023) C. Accettura _et al._ , Eur. Phys. J. C 83, 864 (2023), arXiv:2303.08533 [physics.acc-ph] .
  * Capdevilla  _et al._ (2021a) R. Capdevilla, D. Curtin, Y. Kahn, and G. Krnjaic, Phys. Rev. D 103, 075028 (2021a), arXiv:2006.16277 [hep-ph] .
  * Huang  _et al._ (2021) G.-y. Huang, F. S. Queiroz, and W. Rodejohann, Phys. Rev. D 103, 095005 (2021), arXiv:2101.04956 [hep-ph] .
  * Capdevilla  _et al._ (2022a) R. Capdevilla, D. Curtin, Y. Kahn, and G. Krnjaic, Phys. Rev. D 105, 015028 (2022a), arXiv:2101.10334 [hep-ph] .
  * Liu and Xie (2021) W. Liu and K.-P. Xie, JHEP 04, 015 (2021), arXiv:2101.10469 [hep-ph] .
  * Han  _et al._ (2021d) T. Han, S. Li, S. Su, W. Su, and Y. Wu, Phys. Rev. D 104, 055029 (2021d), arXiv:2102.08386 [hep-ph] .
  * Asadi  _et al._ (2021) P. Asadi, R. Capdevilla, C. Cesarotti, and S. Homiller, JHEP 10, 182 (2021), arXiv:2104.05720 [hep-ph] .
  * Sen  _et al._ (2022) C. Sen, P. Bandyopadhyay, S. Dutta, and A. KT, Eur. Phys. J. C 82, 230 (2022), arXiv:2107.12442 [hep-ph] .
  * Bandyopadhyay  _et al._ (2022) P. Bandyopadhyay, A. Karan, R. Mandal, and S. Parashar, Eur. Phys. J. C 82, 916 (2022), arXiv:2108.06506 [hep-ph] .
  * Liu  _et al._ (2022a) W. Liu, K.-P. Xie, and Z. Yi, Phys. Rev. D 105, 095034 (2022a), arXiv:2109.15087 [hep-ph] .
  * Casarsa  _et al._ (2022) M. Casarsa, M. Fabbrichesi, and E. Gabrielli, Phys. Rev. D 105, 075008 (2022), arXiv:2111.13220 [hep-ph] .
  * Capdevilla  _et al._ (2022b) R. Capdevilla, D. Curtin, Y. Kahn, and G. Krnjaic, JHEP 04, 129 (2022b), arXiv:2112.08377 [hep-ph] .
  * Bao  _et al._ (2022) Y. Bao, J. Fan, and L. Li, JHEP 08, 276 (2022), arXiv:2203.04328 [hep-ph] .
  * İnan and Kisselev (2023a) S. C. İnan and A. V. Kisselev, J. Phys. G 50, 105002 (2023a), arXiv:2207.03325 [hep-ph] .
  * Liu  _et al._ (2022b) J. Liu, Z.-L. Han, Y. Jin, and H. Li, JHEP 12, 057 (2022b), arXiv:2207.07382 [hep-ph] .
  * Lv  _et al._ (2022) G.-S. Lv, X.-M. Cui, Y.-Q. Li, and Y.-B. Liu, Nucl. Phys. B 985, 116016 (2022).
  * Kwok  _et al._ (2023) T. H. Kwok, L. Li, T. Liu, and A. Rock, (2023), arXiv:2301.05177 [hep-ph] .
  * İnan and Kisselev (2023b) S. C. İnan and A. V. Kisselev, (2023b), arXiv:2301.08585 [hep-ph] .
  * Jueid and Nasri (2023) A. Jueid and S. Nasri, Phys. Rev. D 107, 115027 (2023), arXiv:2301.12524 [hep-ph] .
  * Jueid  _et al._ (2024) A. Jueid, T. A. Chowdhury, S. Nasri, and S. Saad, Phys. Rev. D 109, 075011 (2024), arXiv:2306.01255 [hep-ph] .
  * Dermisek  _et al._ (2023) R. Dermisek, K. Hermanek, N. McGinnis, and S. Yoon, Phys. Rev. D 108, 055019 (2023), arXiv:2306.13212 [hep-ph] .
  * Belfkir  _et al._ (2024) M. Belfkir, T. A. Chowdhury, and S. Nasri, Phys. Lett. B 852, 138605 (2024), arXiv:2307.16111 [hep-ph] .
  * Dermisek  _et al._ (2024) R. Dermisek, K. Hermanek, T. Lee, N. McGinnis, and S. Yoon, Phys. Rev. D 109, 095003 (2024), arXiv:2311.05078 [hep-ph] .
  * Asadi  _et al._ (2023) P. Asadi, A. Radick, and T.-T. Yu, (2023), arXiv:2312.03826 [hep-ph] .
  * Liu  _et al._ (2024b) D. Liu, L.-T. Wang, and K.-P. Xie, JHEP 04, 084 (2024b), arXiv:2312.09117 [hep-ph] .
  * Cirelli  _et al._ (2006) M. Cirelli, N. Fornengo, and A. Strumia, Nucl. Phys. B 753, 178 (2006), arXiv:hep-ph/0512090 .
  * Ade  _et al._ (2016) P. A. R. Ade _et al._ (Planck), Astron. Astrophys. 594, A13 (2016), arXiv:1502.01589 [astro-ph.CO] .
  * Bottaro  _et al._ (2022a) S. Bottaro, D. Buttazzo, M. Costa, R. Franceschini, P. Panci, D. Redigolo, and L. Vittorio, Eur. Phys. J. C 82, 31 (2022a), arXiv:2107.09688 [hep-ph] .
  * Bottaro  _et al._ (2022b) S. Bottaro, D. Buttazzo, M. Costa, R. Franceschini, P. Panci, D. Redigolo, and L. Vittorio, Eur. Phys. J. C 82, 992 (2022b), arXiv:2205.04486 [hep-ph] .
  * Arkani-Hamed  _et al._ (2006) N. Arkani-Hamed, A. Delgado, and G. F. Giudice, Nucl. Phys. B 741, 108 (2006), arXiv:hep-ph/0601041 .
  * Thomas and Wells (1998) S. D. Thomas and J. D. Wells, Phys. Rev. Lett. 81, 34 (1998), arXiv:hep-ph/9804359 .
  * Barr and Scoville (2015) A. Barr and J. Scoville, JHEP 04, 147 (2015), arXiv:1501.02511 [hep-ph] .
  * Belanger  _et al._ (2019) G. Belanger _et al._ , JHEP 02, 186 (2019), arXiv:1811.05478 [hep-ph] .
  * Jana  _et al._ (2022) S. Jana, N. Okada, and D. Raut, Eur. Phys. J. C 82, 927 (2022), arXiv:1911.09037 [hep-ph] .
  * Chiang  _et al._ (2021) C.-W. Chiang, G. Cottin, Y. Du, K. Fuyuto, and M. J. Ramsey-Musolf, JHEP 01, 198 (2021), arXiv:2003.07867 [hep-ph] .
  * Calibbi  _et al._ (2021) L. Calibbi, F. D’Eramo, S. Junius, L. Lopez-Honorez, and A. Mariotti, JHEP 05, 234 (2021), arXiv:2102.06221 [hep-ph] .
  * Schwaller and Zurita (2014) P. Schwaller and J. Zurita, JHEP 03, 060 (2014), arXiv:1312.7350 [hep-ph] .
  * Baer  _et al._ (2014) H. Baer, A. Mustafayev, and X. Tata, Phys. Rev. D 89, 055007 (2014), arXiv:1401.1162 [hep-ph] .
  * Agin  _et al._ (2024) D. Agin, B. Fuks, M. D. Goodsell, and T. Murphy, Phys. Lett. B 853, 138597 (2024), arXiv:2311.17149 [hep-ph] .
  * Han  _et al._ (2014) Z. Han, G. D. Kribs, A. Martin, and A. Menon, Phys. Rev. D 89, 075007 (2014), arXiv:1401.1235 [hep-ph] .
  * Low and Wang (2014) M. Low and L.-T. Wang, JHEP 08, 161 (2014), arXiv:1404.0682 [hep-ph] .
  * Mahbubani  _et al._ (2017) R. Mahbubani, P. Schwaller, and J. Zurita, JHEP 06, 119 (2017), [Erratum: JHEP 10, 061 (2017)], arXiv:1703.05327 [hep-ph] .
  * Fukuda  _et al._ (2018) H. Fukuda, N. Nagata, H. Otono, and S. Shirai, Phys. Lett. B 781, 306 (2018), arXiv:1703.09675 [hep-ph] .
  * Bharucha  _et al._ (2018) A. Bharucha, F. Brümmer, and N. Desai, JHEP 11, 195 (2018), arXiv:1804.02357 [hep-ph] .
  * Calibbi  _et al._ (2018) L. Calibbi, L. Lopez-Honorez, S. Lowette, and A. Mariotti, JHEP 09, 037 (2018), arXiv:1805.04423 [hep-ph] .
  * Filimonova and Westhoff (2019) A. Filimonova and S. Westhoff, JHEP 02, 140 (2019), arXiv:1812.04628 [hep-ph] .
  * Belyaev  _et al._ (2021) A. Belyaev, S. Prestel, F. Rojas-Abbate, and J. Zurita, Phys. Rev. D 103, 095006 (2021), arXiv:2008.08581 [hep-ph] .
  * Canepa  _et al._ (2020) A. Canepa, T. Han, and X. Wang, Ann. Rev. Nucl. Part. Sci. 70, 425 (2020), arXiv:2003.05450 [hep-ph] .
  * ATLAS Collaboration (2013) ATLAS Collaboration, Phys. Rev. D 88, 112006 (2013), arXiv:1310.3675 [hep-ex] .
  * CMS Collaboration (2015) CMS Collaboration, JHEP 01, 096 (2015), arXiv:1411.6006 [hep-ex] .
  * ATLAS Collaboration (2018a) ATLAS Collaboration, JHEP 06, 022 (2018a), arXiv:1712.02118 [hep-ex] .
  * CMS Collaboration (2020) CMS Collaboration, Phys. Lett. B 806, 135502 (2020), arXiv:2004.05153 [hep-ex] .
  * ATLAS Collaboration (2022) ATLAS Collaboration, Eur. Phys. J. C 82, 606 (2022), arXiv:2201.02472 [hep-ex] .
  * Khachatryan  _et al._ (2016) V. Khachatryan _et al._ (CMS), Phys. Lett. B 759, 9 (2016), arXiv:1512.08002 [hep-ex] .
  * ATLAS Collaboration (2018b) ATLAS Collaboration, Phys. Rev. D 97, 052010 (2018b), arXiv:1712.08119 [hep-ex] .
  * CMS Collaboration (2018) CMS Collaboration, Phys. Lett. B 782, 440 (2018), arXiv:1801.01846 [hep-ex] .
  * ATLAS Collaboration (2020) ATLAS Collaboration, Phys. Rev. D 101, 052005 (2020), arXiv:1911.12606 [hep-ex] .
  * ATLAS Collaboration (2021a) ATLAS Collaboration, Eur. Phys. J. C 81, 1118 (2021a), arXiv:2106.01676 [hep-ex] .
  * CMS Collaboration (2022a) CMS Collaboration, JHEP 04, 147 (2022a), arXiv:2106.14246 [hep-ex] .
  * CMS Collaboration (2022b) CMS Collaboration, JHEP 04, 091 (2022b), arXiv:2111.06296 [hep-ex] .
  * ATLAS Collaboration (2021b) ATLAS Collaboration, Phys. Rev. D 103, 112006 (2021b), arXiv:2102.10874 [hep-ex] .
  * CMS Collaboration (2021) CMS Collaboration, JHEP 11, 153 (2021), arXiv:2107.13021 [hep-ex] .
  * Han  _et al._ (2018) T. Han, S. Mukhopadhyay, and X. Wang, Phys. Rev. D 98, 035026 (2018), arXiv:1805.00015 [hep-ph] .
  * Saito  _et al._ (2019) M. Saito, R. Sawada, K. Terashi, and S. Asai, Eur. Phys. J. C 79, 469 (2019), arXiv:1901.02987 [hep-ph] .
  * de Blas  _et al._ (2018) J. de Blas _et al._ (CLIC),  3/2018 (2018), 10.23731/CYRM-2018-003, arXiv:1812.02093 [hep-ph] .
  * Han  _et al._ (2021e) T. Han, Z. Liu, L.-T. Wang, and X. Wang, Phys. Rev. D 103, 075004 (2021e), arXiv:2009.11287 [hep-ph] .
  * Capdevilla  _et al._ (2021b) R. Capdevilla, F. Meloni, R. Simoniello, and J. Zurita, JHEP 06, 133 (2021b), arXiv:2102.11292 [hep-ph] .
  * Bottaro  _et al._ (2021) S. Bottaro, A. Strumia, and N. Vignaroli, JHEP 06, 143 (2021), arXiv:2103.12766 [hep-ph] .
  * Black  _et al._ (2022) K. Black, T. Bose, Y. Chen, S. Dasu, H. Jia, D. Pinna, V. Sharma, N. Venkatasubramanian, and C. Vuosalo, in  _Snowmass 2021_ (2022) arXiv:2205.10404 [hep-ex] .
  * Bandyopadhyay  _et al._ (2024) P. Bandyopadhyay, S. Parashar, C. Sen, and J. Song, (2024), arXiv:2401.02697 [hep-ph] .
  * Curtin  _et al._ (2018) D. Curtin, K. Deshpande, O. Fischer, and J. Zurita, JHEP 07, 024 (2018), arXiv:1712.07135 [hep-ph] .
  * Fukuda  _et al._ (2020) H. Fukuda, N. Nagata, H. Oide, H. Otono, and S. Shirai, Phys. Rev. Lett. 124, 101801 (2020), arXiv:1910.08065 [hep-ph] .
  * ATLAS Collaboration (2019) ATLAS Collaboration, “Performance of tracking and vertexing techniques for a disappearing track plus soft track signature with the ATLAS detector,” (2019).
  * Franceschini and Zhao (2023) R. Franceschini and X. Zhao, Eur. Phys. J. C 83, 552 (2023), arXiv:2212.11900 [hep-ph] .
  * Fukuda  _et al._ (2024) H. Fukuda, T. Moroi, A. Niki, and S.-F. Wei, JHEP 02, 214 (2024), arXiv:2310.07162 [hep-ph] .
  * Frixione  _et al._ (1993) S. Frixione, M. L. Mangano, P. Nason, and G. Ridolfi, Phys. Lett. B 319, 339 (1993), arXiv:hep-ph/9310350 .
  * Calzolari  _et al._ (shed) D. Calzolari, D. Schulte, _et al._ , (To be published.).
  * Alwall  _et al._ (2014) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, JHEP 07, 079 (2014), arXiv:1405.0301 [hep-ph] .
  * Sjöstrand  _et al._ (2015) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands, Comput. Phys. Commun. 191, 159 (2015), arXiv:1410.3012 [hep-ph] .
  * Agostinelli  _et al._ (2003) S. Agostinelli _et al._ (GEANT4), Nucl. Instrum. Meth. A 506, 250 (2003).
  * Sailer  _et al._ (2020) A. Sailer, G. Ganis, P. Mato, M. Petrič, and G. A. Stewart, EPJ Web Conf. 245, 10002 (2020).
  * Ganis  _et al._ (2022) G. Ganis, C. Helsens, and V. Völkl, Eur. Phys. J. Plus 137, 149 (2022), arXiv:2111.09874 [hep-ex] .
  * Ai  _et al._ (2022) X. Ai _et al._ , Comput. Softw. Big Sci. 6, 8 (2022), arXiv:2106.13593 [physics.ins-det] .
  * Mokhov and James (2017) N. V. Mokhov and C. C. James,  (2017), 10.2172/1462233.
  * Ferrari  _et al._ (2005) A. Ferrari, P. R. Sala, A. Fasso, and J. Ranft,  (2005), 10.2172/877507.
  * Böhlen  _et al._ (2014) T. T. Böhlen, F. Cerutti, M. P. W. Chin, A. Fassò, A. Ferrari, P. G. Ortega, A. Mairani, P. R. Sala, G. Smirnov, and V. Vlachoudis, Nucl. Data Sheets 120, 211 (2014).
  * Ibe  _et al._ (2024) M. Ibe, Y. Nakayama, and S. Shirai, JHEP 03, 012 (2024), arXiv:2312.08087 [hep-ph] .
  * Mitridate  _et al._ (2017) A. Mitridate, M. Redi, J. Smirnov, and A. Strumia, JCAP 05, 006 (2017), arXiv:1702.01141 [hep-ph] .