# The Gas Accretion Rate of Galaxies over $z\approx 0-1.3$

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

We present here estimates of the average rates of accretion of neutral gas
onto main-sequence galaxies and the conversion of atomic gas to molecular gas
in these galaxies at two key epochs in galaxy evolution: (i) $z\approx
1.3-1.0$, towards the end of the epoch of peak star-formation activity in the
Universe, and (ii) $z\approx 1-0$, when the star-formation activity declines
by an order of magnitude. We determine the net gas accretion rate
$\rm{R_{Acc}}$ and the molecular gas formation rate $\rm{R_{Mol}}$ by
combining the relations between the stellar mass and the atomic gas mass, the
molecular gas mass, and the star-formation rate (SFR) at three epochs,
$z=1.3$, $z=1.0$, and $z=0$, with the assumption that galaxies evolve
continuously on the star-forming main-sequence. We find that, for all
galaxies, $\rm{R_{Acc}}$ is far lower than the average SFR $\rm{R_{SFR}}$ at
$z\approx 1.3-1.0$; however, $\rm{R_{Mol}}$ is similar to $\rm{R_{SFR}}$
during this interval. Conversely, both $\rm{R_{Mol}}$ and $\rm{R_{Acc}}$ are
significantly lower than $\rm{R_{SFR}}$ over the later interval, $z\approx
1-0$. We find that massive main-sequence galaxies had already acquired most of
their present-day baryonic mass by $z\approx 1.3$. At $z\approx 1.3-1.0$, the
rapid conversion of the existing atomic gas to molecular gas was sufficient to
maintain a high average SFR, despite the low net gas accretion rate. However,
at later times, the combination of the lower net gas accretion rate and the
lower molecular gas formation rate leads to a decline in the fuel available
for star-formation, and results in the observed decrease in the SFR density of
the Universe over the last 8 Gyr.

Galaxy evolution — Neutral hydrogen clouds — High-$z$ galaxies

††software:  astropy (Astropy Collaboration et al., 2013)

##  1 Introduction

The evolution of galaxies is driven by the baryon cycle in which the baryonic
constituents of galaxies and their circumgalactic mediums (CGMs) interact with
each other and convert from one form to another (e.g. Péroux & Howk, 2020).
The key processes in the baryon cycle are (i) the accretion of mostly ionized
gas from the CGM onto the “disks” of galaxies, forming neutral atomic hydrogen
(Hi) in the disks, (ii) the cooling of Hi and its conversion to molecular
hydrogen (${\rm H_{2}}$), (iii) the gravitational collapse and fragmentation
of molecular clouds to form stars, and (iv) the expulsion of gas from the
interstellar mediums (ISMs) of galaxies in outflows driven by stars or active
galactic nuclei (AGNs), with some fraction of this gas later returning to the
galaxy disk.

Measurements of the rates at which the above processes occur in galaxies, how
the rates compare to each other, and how they evolve with redshift are
critical to understanding galaxy evolution. Unfortunately, only the redshift
evolution of the stellar properties of galaxies and their star-formation rates
(SFR) are well determined today. For example, we have known for over two
decades that the SFR density of the Universe peaks in the redshift range
$z\approx 1-3$ and then declines by an order of magnitude from $z\approx 1$ to
$z\approx 0$ (e.g. Madau & Dickinson, 2014). Further, $\approx 90\%$ of the
star-formation activity of the Universe out to $z\approx 2$ occurs on the
“star-forming main sequence” (e.g. Rodighiero et al., 2011) — a tight
relationship between the SFR and the stellar mass (${\rm{M}_{*}}$) of galaxies
(e.g. Noeske et al., 2007; Whitaker et al., 2014). At a fixed stellar mass,
the SFR of galaxies on the main sequence declines by a factor of $\approx 10$
from $z\approx 1$ to the present time (e.g. Whitaker et al., 2014; Popesso et
al., 2022).

Unlike the SFR of a galaxy, it is very challenging to observationally
determine the gas accretion rate and the rate at which Hi is converted to
${\rm H_{2}}$. While there have been suggestions in the literature that a low
gas accretion rate at late times might account for various observational
results (e.g. Bouché et al., 2010; Møller et al., 2013; Walter et al., 2020;
Chowdhury et al., 2020), the lack of actual measurements of the accretion rate
has been a key limitation in our understanding of the processes that drive the
redshift evolution of the star-formation activity of the Universe. For
example, Walter et al. (2020) used the redshift evolution of the cosmological
SFR density, the cosmological Hi mass density, and the cosmological ${\rm
H_{2}}$ mass density to estimate the rate of global flow of gas, and global
conversion of Hi to ${\rm H_{2}}$. However, such global rates, averaged over
cosmological volumes, are difficult to interpret and do not provide any
information on the differences between different galaxy populations, such as
the dependence of the accretion rate on galaxy stellar mass, environment,
morphology, etc. Conversely, Bouché et al. (2010) used a toy “reservoir”
model, in which gas accretion is quenched above a halo mass of $\approx
10^{11}\,\rm{M}_{\odot}$, to argue that the decline in SFR density from
$z\approx 2$ may be driven by the decline in the gas accretion rate. This
model was found to yield the observed main-sequence and Tully-Fisher scaling
relations at different redshifts (see also Møller et al., 2013, for the mass-
metallicity relation). However, its predicted gas fractions ($\approx 30-45$%
at $z\approx 1.2$) are lower than the values recently measured ($\approx 80$%
at $z\approx 1.3$) by Chowdhury et al. (2022a).

The net gas accretion rate (i.e. the difference between the gas inflow and
outflow rates) and the ${\rm H_{2}}$ formation rate can be determined if one
knows the dependences of the Hi and ${\rm H_{2}}$ masses of galaxies on their
stellar masses (i.e. the $\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ and
${\rm{M_{Mol}}}-{\rm{M}_{*}}$ scaling relations; Scoville et al., 2017; Bera
et al., 2023a), along with the standard assumption that star-forming galaxies
evolve along the main sequence (e.g. Renzini, 2009; Peng et al., 2010;
Leitner, 2012; Speagle et al., 2014; Scoville et al., 2017). For the molecular
component, various observational tracers (primarily CO rotational lines and
dust continuum emission) have been used to determine the
${\rm{M_{Mol}}}-{\rm{M}_{*}}$ relation of galaxies out to $z\approx 5$ (e.g.
Tacconi et al., 2020). However, the weakness of the Hi 21 cm line, the only
tracer of the Hi mass of galaxies, has meant that, until very recently,
estimates of the $\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ relation were limited
to the local Universe (e.g. Catinella et al., 2018; Parkash et al., 2018).

The Hi 21 cm stacking approach, based on combining the Hi 21 cm emission
signals from a large number of galaxies with accurately known positions and
redshifts, allows one to overcome the intrinsic weakness of the Hi 21 cm line
and measure the average Hi properties of galaxy populations at cosmological
distances (Zwaan, 2000; Chengalur et al., 2001). This approach has been
recently used to measure the Hi properties of star-forming galaxies out to
$z\approx 1$ (e.g. Bera et al., 2019, 2022, 2023b; Chowdhury et al., 2020,
2021, 2022b). Recently, Chowdhury et al. (2022c) applied the Hi 21 cm stacking
approach to data from the Giant Metrewave Radio Telescope (GMRT) Cold-Hi AT
$z\approx 1$ (CAT$z1$) survey, a 510-hr Hi 21 cm emission survey of galaxies
at $z=0.74-1.45$ (Chowdhury et al., 2022b), to obtain the first measurement of
the $\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ scaling relation at $z\approx 1$.

In this Letter, we combine measurements of the
$\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ scaling relation at $z\approx 1$ from
the GMRT-CAT$z1$ survey with estimates of the star-forming main-sequence
relation and the ${\rm{M_{Mol}}}-{\rm{M}_{*}}$ scaling relation from the
literature, and the assumption of the continuity of galaxy evolution along the
main sequence, to determine the average ${\rm H_{2}}$ formation rate and the
average net gas accretion rate in star-forming galaxies over the redshift
intervals $z\approx 1.3-1.0$ and $z\approx 1.0-0$.

Throughout this work, we assume a Chabrier initial mass function (IMF) for
estimates of the stellar masses and SFRs. Further, we use a flat Lambda-cold
dark matter cosmology, with $\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$, and
$H_{0}=70$ km s-1 Mpc-1.

##  2 Determination of the Gas Accretion Rate and the ${\rm H_{2}}$ Formation
Rate

###  2.1 Formalism

The formalism used in this work to determine the net gas accretion rate and
the ${\rm H_{2}}$ formation rate of star-forming galaxies was introduced by
Scoville et al. (2017), and recently refined by Bera et al. (2023a) to include
the $\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ relation. In this approach, the
build-up of the stellar mass of a galaxy, with an initial stellar mass
(${\textrm{M}_{*,i}}$) at the initial epoch ($t_{i}$), is tracked to determine
the final stellar mass (${\textrm{M}_{*,f}}$) at the final epoch ($t_{f}$).
This is done by assuming that the galaxy remains on the star-forming main-
sequence over the entire period, such that

${\textrm{M}_{*,f}}={\textrm{M}_{*,i}}+\int_{t_{i}}^{t_{f}}(1-f_{return})\
\mathrm{SFR_{MS}}({\rm{M}_{*}},t)\ \mathrm{dt}\;$ (1)

where $\mathrm{SFR_{MS}}({\rm{M}_{*}},t)$ is the main-sequence relation at the
time $t$, and $f_{return}$ is the fraction of the stellar mass returned to the
ISM via stellar winds or supernovae ($f_{return}=0.41$ for a Chabrier IMF;
Leitner & Kravtsov, 2011; Madau & Dickinson, 2014). We define here the average
star-formation rate $\rm{R_{SF}}$ between the epochs $t_{i}$ and $t_{f}$, such
that

${\textrm{M}_{*,f}}={\textrm{M}_{*,i}}+(1-f_{return})\ {\rm{R_{SF}}}\Delta t$
(2)

where $\Delta t=t_{f}-t_{i}$ is the time interval between the two epochs of
interest. Measurements of the star-forming main-sequence (e.g. Whitaker et
al., 2014; Popesso et al., 2022) can then be used to track the evolution of
the stellar mass of a main-sequence galaxy using Equation 1, and this can be
combined with Equation 2 to determine $\rm{R_{SF}}$ using the following
relation.

${\rm{R_{SF}}}=\frac{({\textrm{M}_{*,f}}-{\textrm{M}_{*,i}})}{(1-f_{return})\Delta
t}\;$ (3)

Next, the final molecular gas reservoir111Throughout this Letter,
${\rm{M_{Mol}}}$ and ${\rm{M_{Atom}}}$ are used to refer to gas masses that
include the mass contribution of Helium, ${\rm{M_{Mol,f}}}$, of the galaxy at
time $t_{f}$ is related to the initial molecular gas reservoir,
${\rm{M_{Mol,i}}}$, at time $t_{i}$ via

${\rm{M_{Mol,f}}}={\rm{M_{Mol,i}}}-\rm{R_{SF}}\Delta t+\rm{R_{Mol}}\Delta t$
(4)

where $\rm{R_{Mol}}$ is the average molecular gas formation rate between the
two epochs of interest, i.e. the difference between the rate at which Hi is
converted to ${\rm H_{2}}$ and that at which ${\rm H_{2}}$ is dissociated to
Hi. Given that the initial and final stellar masses are known, one can use
measurements of the molecular gas mass as a function of stellar mass at both
epochs to determine ${\rm{M_{Mol,f}}}$ and ${\rm{M_{Mol,i}}}$ which can then
be combined with Equation 4 to obtain $\rm{R_{Mol}}$ using the following
relation.

$\rm{R_{Mol}}=\rm{R_{SF}}+({\rm{M_{Mol,f}}}-{\rm{M_{Mol,i}}})/\Delta t$ (5)

Finally, the neutral atomic gas mass, ${\rm{M_{Atom,f}}}$, of the galaxy at
time $t_{f}$ is related to the initial neutral atomic gas mass,
${\rm{M_{Atom,i}}}$, at time $t_{i}$ via the relation

${\rm{M_{Atom,f}}}={\rm{M_{Atom,i}}}-\rm{R_{Mol}}\Delta t+\rm{R_{Acc}}\Delta
t$ (6)

where $\rm{R_{Acc}}$ is the _net_ accretion rate of neutral atomic gas onto
the disk of the galaxy, i.e. the difference between the neutral atomic gas
inflow and outflow rates. Again, measurements of $\rm{M_{H{\textsc{i}}}}$ as a
function of stellar mass at the two epochs can be used to determine
${\rm{M_{Atom,f}}}$ and ${\rm{M_{Atom,i}}}$, which can be combined with
Equation 6 to obtain $\rm{R_{Acc}}$ using the following relation:

$\rm{R_{Acc}}=\rm{R_{Mol}}+({\rm{M_{Atom,f}}}-{\rm{M_{Atom,i}}})/\Delta t$ (7)

We note that in the above formalism, the stellar mass returned to the ISM is
assumed to be mostly in the ionised state and is hence not explicitly included
in Equations 4-7. It is entirely possible that some fraction of this gas cools
to form Hi or ${\rm H_{2}}$ over the time interval between $t_{i}$ and
$t_{f}$. However, we emphasize that even assuming that _all_ of the stellar
mass that is returned to the ISM cools to form Hi would not significantly
alter the key conclusions of this work.

We use the above formalism to estimate $\rm{R_{SF}}$, $\rm{R_{Mol}}$, and
$\rm{R_{Acc}}$ over two redshift intervals: (a) $z\approx 1.3$ to $z\approx
1.0$, covering the end of the epoch of peak cosmic star-formation activity,
and (b) $z\approx 1.0$ to $z\approx 0$, covering the decline of star-formation
activity over the last $\approx 8$ Gyr (e.g. Madau & Dickinson, 2014). The
choice of these redshift intervals is driven by our recent measurements of the
average Hi content of stellar mass-matched samples of star-forming galaxies at
$z\approx 1.0$ and $z\approx 1.3$ (Chowdhury et al., 2022d), and the
measurement of the $\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ scaling relation over
the redshift range $0.74\leq z\leq 1.45$ (Chowdhury et al., 2022c). In the
following subsections, we briefly describe the relations that were used to
estimate SFR, ${\rm{M_{Mol}}}$, and ${\rm{M_{Atom}}}$ for star-forming
galaxies at $z\approx 1.3$, $z\approx 1.0$, and $z\approx 0$, specifically (i)
the star-forming main-sequence relation (e.g. Whitaker et al., 2014; Popesso
et al., 2022), (ii) the ${\rm{M_{Mol}}}-{\rm{M}_{*}}$ relation (Tacconi et
al., 2020), and (iii) the $\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ relation
(Catinella et al., 2018; Chowdhury et al., 2022c).

###  2.2 The Star-Forming Main-Sequence

We use the star-forming main-sequence relation of Popesso et al. (2022) to
estimate the build-up of the stellar mass of galaxies via Equation 1. These
authors converted a large number of literature measurements of the main
sequence at $0<z<6$ to a common calibration, to find that the data are
consistent with a relation that is linear at low stellar masses and flattens
above a characteristic stellar mass (see also Whitaker et al., 2014). We will
use the following main-sequence relation, obtained via a second-order
polynomial fit to the dependence of $\log[{\rm SFR}]$ on $\log[{\rm{M}_{*}}]$
(Popesso et al., 2022):

${\log}[{\rm
SFR}({\rm{M}_{*}},t)]=(a_{1}t+b_{1}){\log}[{\rm{M}_{*}}]+b_{2}{\log}^{2}[{\rm{M}_{*}}]+(b_{0}+a_{0}t)$
(8)

where $t$ is the age of the Universe at the epoch of interest and the fit
parameters are $a_{0}=0.20\pm 0.02$, $a_{1}=-0.034\pm 0.002$,
$b_{0}=-26.134\pm 0.015$, $b_{1}=4.722\pm 0.012$, and $b_{2}=-0.1925\pm
0.0011$ (Popesso et al., 2022). The main-sequence relations at $z\approx 0$,
$z\approx 1.0$, and $z\approx 1.3$, obtained using Equation 8, are shown as
the solid blue curves in Figure 1[A].

In passing, we emphasize that the main conclusions of this work are unchanged
on using the second main-sequence relation of Popesso et al. (2022), which
asymptotically approaches a maximum SFR value with increasing stellar mass, or
the main-sequence relation of Whitaker et al. (2014), which is commonly used
in the literature.

###  2.3 The ${\rm{M_{Mol}}}-{\rm{M}_{*}}$ Relation

For the ${\rm{M_{Mol}}}-{\rm{M}_{*}}$ relation, we follow Tacconi et al.
(2020), who compiled literature estimates of the molecular gas mass of 2052
galaxies with ${\rm{M}_{*}}\approx 10^{9}-10^{11.5}~{}\rm{M}_{\odot}$ at
$0\lesssim z\lesssim 5.2$. The molecular gas masses were estimated using
different tracers such as the CO rotational lines (e.g. Daddi et al., 2010;
Tacconi et al., 2013; Saintonge et al., 2017), the far-infrared dust continuum
(e.g. Santini et al., 2014), and the 1 mm dust continuum (e.g. Scoville et
al., 2016). Tacconi et al. (2020) corrected the measurements obtained from the
different tracers for any zero-point offsets, to put them on the same scale,
and obtained the following relation for the dependence of the molecular gas
mass of main-sequence galaxies on redshift and stellar mass:222We note that
the relation in Tacconi et al. (2020) also includes a term giving the
dependence on the offset from the main-sequence relation. We have dropped this
term as our study is restricted to galaxies lying on the main sequence.

$\log\left[\frac{{\rm{M_{Mol}}}}{{\rm{M}_{*}}}\right]=A+B\left[\log(1+z)-F\right]^{2}+D\left[\log{\rm{M}_{*}}-10.7\right]$
(9)

where the coefficients $A=0.06\pm 0.2$, $B=-3.33\pm 0.2$, $F=0.65\pm 0.05$,
and $D=-0.41\pm 0.03$ (note that the quoted uncertainties are $2\sigma$
errors; Tacconi et al., 2020). We will use Equation 9 to determine the
${\rm{M_{Mol}}}-{\rm{M}_{*}}$ relation at $z\approx 0$, $z\approx 1.0$, and
$z\approx 1.3$; these scaling relations are shown as the solid blue lines in
Figure 1[B].

###  2.4 The ${\rm{M_{Atom}}}-{\rm{M}_{*}}$ Relation

We determine the $\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ relation at $z\approx
1.0$ and $z\approx 1.3$ from the measurements of the average Hi mass of main-
sequence galaxies at these redshifts in the GMRT-CAT$z$1 survey (Chowdhury et
al., 2022b). This survey covered the Hi 21 cm line of 11,419 main-sequence
galaxies at $z=0.74-1.45$ with stellar masses
${\rm{M}_{*}}=10^{9}-10^{11.5}~{}\rm{M}_{\odot}$. Chowdhury et al. (2022c)
divided the sample of 11,419 galaxies into three stellar-mass subsamples, and
measured the average Hi mass of the galaxies in each subsample, and fitted a
linear relation to the $\rm{M_{H{\textsc{i}}}}$ and ${\rm{M}_{*}}$ estimates
to determine the $\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ relation at $z\approx
1$. They found that the slope of the $\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$
relation at $z\approx 1$ is consistent with that at $z\approx 0$. However, the
intercept of the relation at $z\approx 1$ is higher than that of the local
relation by a factor of $\approx 3.5$.

Chowdhury et al. (2022c) also determined the
$\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ relation for blue (NUV$-$r$<4$) galaxies
at $z\approx 0$ from $\rm{M_{H{\textsc{i}}}}$ measurements of galaxies covered
in the extended GALEX Arecibo Sloan Digital Sky Survey (xGASS; Catinella et
al., 2018): xGASS is an Arecibo telescope Hi 21 cm survey of a sample of
stellar mass-selected galaxies at $z\approx 0$, with stellar masses in the
range $10^{9.0}-10^{11.5}\,\rm M_{\odot}$. They obtained the relation
$\log[\rm{M_{H{\textsc{i}}}}]=(0.38\pm 0.05)\log[{\rm M_{*,10}}]+(9.634\pm
0.019)$, where ${\rm
M_{*,10}}=\log\left[{\rm{M}_{*}}/10^{10}\rm{M}_{\odot}\right]$. Since the
slope of the $\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ relation for the GMRT-
CAT$z$1 galaxies at $z\approx 0.74-1.45$ is the same as that of the relation
at $z\approx 0$, we will assume that the slope of the
$\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ relation is constant, and equal to
$0.38$ (i.e. the xGASS slope), over the entire redshift range $z\approx
0-1.3$.

Next, Chowdhury et al. (2022d) subdivided the DEEP2 galaxies into two redshift
bins, at $0.74<z\leq 1.25$ and $1.25<z\leq 1.45$, to measure the average Hi
mass of main-sequence galaxies at average redshifts of $\langle
z\rangle\approx 1$ and $\langle z\rangle\approx 1.3$. They obtained average Hi
masses of $\langle\rm{M_{H{\textsc{i}}}}\rangle=(3.36\pm 0.64)\times
10^{10}~{}\rm{M}_{\odot}$ at $z\approx 1.3$ and
$\langle\rm{M_{H{\textsc{i}}}}\rangle=(1.06\pm 0.19)\times
10^{10}~{}\rm{M}_{\odot}$ at $z\approx 1.0$, both for average stellar masses
of $\langle{\rm{M}_{*}}\rangle\approx 10^{10}~{}\rm{M}_{\odot}$. We determine
the intercepts of the $\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ relations at
$z\approx 1$ and $z\approx 1.3$ by requiring that the relations are consistent
with the above $\langle\rm{M_{H{\textsc{i}}}}\rangle$ measurements for
galaxies with $\langle{\rm{M}_{*}}\rangle=10^{10}~{}\rm{M}_{\odot}$ at each
redshift. This yields the scaling relations

$\log[\rm{M_{H{\textsc{i}}}}]=0.38\log[{\rm M_{*,10}}]+(10.063\pm
0.082)\;\;\;\;\;(z\approx 1)$ (10)

and

$\log[\rm{M_{H{\textsc{i}}}}]=0.38\log[{\rm M_{*,10}}]+(10.564\pm
0.087)\;\;\;\;\;(z\approx 1.3)$ (11)

We note that the above $\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ relations yield,
at each redshift, the mean $\rm{M_{H{\textsc{i}}}}$ value for a distribution
of galaxies at a given ${\rm{M}_{*}}$. However, the main-sequence relation of
Equation 8 and the $\rm{M_{H_{2}}}-{\rm{M}_{*}}$ relation of Equation 9 yield,
respectively, the median SFR and the median $\rm{M_{H_{2}}}$ at a given
stellar mass. For a log-normal distribution of Hi masses, which is typically
the case (e.g Saintonge & Catinella, 2022), the mean Hi mass is different from
the median Hi mass. For consistency with the main-sequence and the
${\rm{M_{Mol}}}-{\rm{M}_{*}}$ relations, we hence correct our mean
$\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ relations to obtain the median
$\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ relations (e.g. Bera et al., 2022;
Chowdhury et al., 2022c). Assuming that the Hi masses of galaxies at $z\approx
1$ and $z\approx 1.3$ follow log-normal distributions with a scatter identical
to the 0.4 dex scatter measured for the $\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$
relation at $z\approx 0$ (Catinella et al., 2018), the median
$\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ relations lie 0.184 dex lower than the
corresponding mean relations. The median $\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$
scaling relations are then:

$\log[\rm{M_{H{\textsc{i}}}}]=0.38\log[{\rm M_{*,10}}]+(9.879\pm
0.082)\;\;\;\;\;(z\approx 1),$ (12)

and

$\log[\rm{M_{H{\textsc{i}}}}]=0.38\log[{\rm M_{*,10}}]+(10.380\pm
0.087)\;\;\;\;\;(z\approx 1.3).$ (13)

Finally, we also correct the above median relations for the mass contribution
of helium, assuming ${\rm{M_{Atom}}}=1.36\ \rm{M_{H{\textsc{i}}}}$, to obtain

$\log[{\rm{M_{Atom}}}]=0.38\log[{\rm M_{*,10}}]+A$ (14)

where the values of the coefficient $A$ are $(9.583\pm 0.019)$ at $z\approx
0$, $(10.012\pm 0.082)$ at $z\approx 1.0$, and $(10.513\pm 0.087)$ at
$z\approx 1.3$. The solid blue lines in Figure 1[C] show the median
${\rm{M_{Atom}}}-{\rm{M}_{*}}$ relation at $z\approx 0$, $z\approx 1.0$, and
$z\approx 1.3$.

In passing, we note that, over the last two decades, different Hi 21 cm
surveys have measured the $\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ scaling
relation for various galaxy populations at $z\approx 0$, each with different
galaxy selection criteria (e.g. Barnes et al., 2001; Giovanelli et al., 2005;
Huang et al., 2012; Bothwell et al., 2013; Dénes et al., 2014; Catinella et
al., 2018; Parkash et al., 2018). The scaling relations obtained from “blind”
Hi 21 cm surveys (e.g. Barnes et al., 2001; Giovanelli et al., 2005) are
biased towards Hi-rich galaxies, yielding
$\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ relations that lie above those obtained
from Hi 21 cm surveys of optically-selected galaxies (e.g. Huang et al., 2012;
Catinella et al., 2018). The $\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ relation at
$z\approx 1$ from the GMRT-CAT$z1$ survey, as well as the main-sequence
relation and the $\rm{M_{H_{2}}}-{\rm{M}_{*}}$ relation used in this work, are
all obtained from optically-selected galaxy samples. We hence chose to use the
scaling relation obtained from the blue galaxies in the optically-selected
xGASS sample (Catinella et al., 2018) as the reference
$\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ relation at $z\approx 0$.

##  3 Results

Figure 1: The panels show (in blue), as a function of stellar mass, [A] the
SFR (Section 2.2), [B] the molecular gas mass (Section 2.3), and [C] the
atomic gas mass (Section 2.4), for main-sequence galaxies at $z\approx 0$,
$z\approx 1.0$ and $z\approx 1.3$. The red curve in Panel [A] shows the
evolutionary track of a galaxy with ${\rm{M}_{*}}=10^{10}~{}\rm{M}_{\odot}$ at
$z\approx 1.3$, obtained using Equation 1. The corresponding evolution of the
molecular gas mass, obtained using Equations 1 and 9, is shown in Panel [B].
Finally, the red points in Panel [C] shows the atomic gas mass of the galaxy
at the three redshifts, obtained using Equation 14. Note that the
${\rm{M_{Atom}}}-{\rm{M}_{*}}$ relation of Equation 14 is not a continuous
function of redshift; the dashed line in the panel is hence only for visual
aid, showing a linear interpolation between the points. See main text for
discussion.

We use the formalism of Section 2.1 with the main-sequence relation (Equation
8), the ${\rm{M_{Mol}}}-{\rm{M}_{*}}$ relation (Equation 9), and the
${\rm{M_{Atom}}}-{\rm{M}_{*}}$ relation (Equation 14), to estimate
$\rm{R_{SF}}$, $\rm{R_{Mol}}$, and $\rm{R_{Acc}}$ as a function of stellar
mass over the redshift intervals $z\approx 1.3-1.0$ and $z\approx 1.0-0.0$. We
emphasize that this implicitly assumes that star-forming galaxies evolve along
the main sequence (e.g. Renzini, 2009; Peng et al., 2010; Speagle et al.,
2014; Scoville et al., 2017).

###  3.1 The Evolution of a Milky-Way-like Galaxy from $z\approx 1.3$ to
$z\approx 0$

We illustrate the formalism of Section 2.1 for the case of a massive galaxy
with ${\rm{M}_{*}}=10^{10}~{}\rm{M}_{\odot}$ at $z\approx 1.3$. We use
Equation 1 and the main-sequence relation of Equation 8 to evolve the galaxy;
the evolutionary track of the galaxy in the SFR$-{\rm{M}_{*}}$ plane is shown
in the red curve in Figure 1[A]. We find that the stellar mass of the galaxy
increases to $\approx 2.1\times 10^{10}~{}\rm{M}_{\odot}$ at $z\approx 1$.
Using Equation 3, we find that the time-averaged SFR of the galaxy between
$z=1.3$ and $z=1.0$ is $\rm{R_{SF}}\approx 16~{}\rm{M}_{\odot}$ yr-1. We use
the ${\rm{M_{Mol}}}-{\rm{M}_{*}}$ relation of Equation 9 to find that the
molecular gas mass of the galaxy increases from ${\rm{M_{Mol}}}\approx
1.2\times 10^{10}~{}\rm{M}_{\odot}$ at $z=1.3$ to ${\rm{M_{Mol}}}\approx
1.3\times 10^{10}~{}\rm{M}_{\odot}$ at $z=1.0$. Using Equation 5, this yields
a time-averaged molecular gas formation rate of $\rm{R_{Mol}}\approx
17~{}\rm{M}_{\odot}{\rm yr^{-1}}$, very similar to the average SFR,
$\rm{R_{SF}}$, over the same interval. Finally, we use the
${\rm{M_{Atom}}}-{\rm{M}_{*}}$ relation of Equation 14 to find that the atomic
gas mass of the galaxy declines from ${\rm{M_{Atom}}}\approx 3.3\times
10^{10}~{}\rm{M}_{\odot}$ at $z=1.3$ to ${\rm{M_{Atom}}}\approx 1.4\times
10^{10}~{}\rm{M}_{\odot}$ at $z=1.0$. We use Equation 7 to find that the net
accretion rate of the galaxy over this interval is $\rm{R_{Acc}}\approx
1~{}\rm{M}_{\odot}$ yr-1. We thus find that the net accretion rate of a main-
sequence galaxy with a stellar mass of ${\rm{M}_{*}}=10^{10}~{}\rm{M}_{\odot}$
at $z\approx 1.3$, over the interval $z=1.3$ to $z=1.0$, is far lower than
both the average SFR and the molecular gas formation rate over the same
period.

Following the same galaxy down to $z\approx 0$, we find that the stellar mass
increases to ${\rm{M}_{*}}\approx 5\times 10^{10}~{}\rm{M}_{\odot}$ (similar
to that of the Milky Way), while the molecular gas mass and atomic gas mass
both decline, to values of ${\rm{M_{Mol}}}\approx 2\times
10^{9}~{}\rm{M}_{\odot}$ and ${\rm{M_{Atom}}}\approx 7\times
10^{9}~{}\rm{M}_{\odot}$, respectively, at $z=0$. We note that that these
stellar, atomic, and molecular gas masses are similar to those of the Milky
Way (e.g. Kalberla & Kerp, 2009). Using Equations 3–7, we obtain
$\rm{R_{SF}}\approx 7~{}\rm{M}_{\odot}$ yr-1, $\rm{R_{Mol}}\approx
5~{}\rm{M}_{\odot}$ yr-1, and $\rm{R_{Acc}}\approx 4~{}\rm{M}_{\odot}$ yr-1
over the 8-Gyr interval from $z=1.0$ to $z=0$. We thus find that both
$\rm{R_{Mol}}$ and $\rm{R_{Acc}}$ are lower than $\rm{R_{SF}}$ during the
interval $z=1.0$ to $z=0$; this is unlike the situation in the interval
$z=1.3$ to $z=1.0$, where $\rm{R_{Mol}}$ was similar to $\rm{R_{SF}}$.
Further, both $\rm{R_{SF}}$ and $\rm{R_{Mol}}$ are significantly lower in the
interval $z=1.0-0$ than at $z=1.3-1.0$; however, the time-averaged net gas
accretion rate $\rm{R_{Acc}}$ over $z=1.0-0$ is higher than that over
$z=1.3-1.0$.

###  3.2 The Net Gas Accretion and Molecular Gas Formation Rates as a Function
of Stellar Mass

Figure 2: The two panels show the redshift-averaged star-formation rate
($\rm{R_{SF}}$, in blue), molecular gas formation rate ($\rm{R_{Mol}}$, in
orange), and net gas accretion rate ($\rm{R_{Acc}}$, in green) of main-
sequence galaxies, as a function of the stellar mass, over [A] $z=1.3$ to
$z=1.0$, and [B] $z=1.0$ to $z=0$. The shaded regions around each curve show
the 68% confidence intervals. The top x-axis of panel [A] shows the initial
stellar mass of the galaxies at $z=1.3$, while the bottom x-axis shows the
final stellar mass, at $z=1.0$. This final stellar mass at $z=1.0$ is plotted
as the top x-axis of panel [B], with the bottom axis of the panel showing the
final stellar mass of the galaxies at $z=0$. The evolution of the stellar
masses was obtained using Equation 1. Panel [A] shows that the molecular gas
formation rate and star-formation rates are similar for all but the most
massive galaxies over the interval $z=1.3$ to $z=1.0$, but the net gas
accretion rate over the same period is much lower than the other two rates.
Panel [B] shows that both the average molecular gas formation rate and the
average net gas accretion rate are lower than the average SFR for the lower-
redshift interval, $z=1.0$ to $z=0$, for galaxies of all stellar masses. See
main text for discussion.

We apply a similar analysis to main-sequence galaxies with initial stellar
masses (at $z=1.3$) in the range ${\rm{M}_{*}}\approx
10^{9}-10^{11.4}~{}\rm{M}_{\odot}$, similar to the range of stellar masses
over which the GMRT-CAT$z1$ survey obtained the
$\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ relation (Chowdhury et al., 2022c).
Figure 2[A] shows, as a function of stellar mass at $z=1.0$, our estimates of
$\rm{R_{SF}}$, $\rm{R_{Mol}}$, and $\rm{R_{Acc}}$ over the interval $z=1.3$ to
$z=1.0$, while Figure 2[B] shows the same quantities for the lower redshift
interval $z=1.0$ to $z=0$, again plotted against the final stellar mass, at
$z=0$. We note that the confidence intervals shown in Figure 2 take into
account the uncertainties in the ${\rm{M_{Atom}}}-{\rm{M}_{*}}$ relation and
the ${\rm{M_{Mol}}}-{\rm{M}_{*}}$ relation, but not that in the main-sequence
relation; the uncertainties in our estimates of $\rm{R_{Mol}}$ and
$\rm{R_{Acc}}$ are dominated by those in the former two scaling relations.

Figure 2[A] shows that, except for the most massive galaxies, the time-
averaged molecular gas formation rate is similar to the net average SFR over
the higher-redshift interval $z=1.3$ to $z=1.0$. Both $\rm{R_{SF}}$ and
$\rm{R_{Mol}}$ are $\gtrsim 10~{}\rm{M}_{\odot}~{}\textrm{yr}^{-1}$ for
galaxies with ${\rm{M}_{*}}\gtrsim 10^{10}~{}\rm{M}_{\odot}$ at $z\approx
1.3$. Remarkably, Figure 2[A] also shows that $\rm{R_{Acc}}$ for main-sequence
galaxies over the entire stellar mass range is far lower than both
$\rm{R_{SF}}$ and $\rm{R_{Mol}}$. Indeed, our estimate of $\rm{R_{Acc}}$ is
consistent with no net accretion of neutral gas onto the disks of star-forming
galaxies over the redshift range $z\approx 1.3$ to $z\approx 1$. This
essentially implies that the Hi that is consumed in molecular gas formation
(and later, star-formation) is not replenished in main-sequence galaxies
towards the end of the epoch of galaxy assembly.

For the most massive galaxies, with ${\rm{M}_{*}}\gtrsim
10^{10.5}~{}\rm{M}_{\odot}$ at $z\approx 1.3$, Fig. 2[A] indicates that the
net gas accretion rate is negative, i.e. that outflows dominate infall in such
galaxies over the redshift range $z\approx 1.3$ to $z\approx 1.0$. Evidence
has indeed been found for stronger outflows in higher stellar-mass galaxies at
these redshifts: for example, Weiner et al. (2009) find that the outflow
velocity in DEEP2 galaxies at $z\approx 1.4$ scales $\propto{\rm SFR}^{0.3}$,
implying that the outflowing material may escape from high-mass galaxies.
However, we emphasize that the result that outflows dominate infall in high-
mass galaxies remains tentative, due to the uncertainties in the different
scaling relations.

In the lower-redshift interval ($z=1$ to $z=0$), Figure 2[B] shows that
$\rm{R_{SF}}\approx 4-12~{}\rm{M}_{\odot}$ yr-1, typically far lower (except
for the lowest stellar-mass galaxies) than $\rm{R_{SF}}$ for the same galaxy
in the higher redshift interval, $z=1.3$ to $z=0$. This is a direct
consequence of the redshift evolution of the main-sequence relation (e.g.
Whitaker et al., 2014; Popesso et al., 2022). Further, it can be seen from the
figure that $\rm{R_{Mol}}$ is lower than $\rm{R_{SF}}$ for galaxies over the
entire stellar-mass range. This is different from the behavior in the higher-
redshift interval, where Figure 2[A] shows that
$\rm{R_{Mol}}\approx\rm{R_{SF}}$ for all, except the most massive, galaxies.
Finally, we find that the net average gas accretion rate of main-sequence
galaxies over the lower redshift interval is comparable to $\rm{R_{Mol}}$ for
galaxies with ${\rm{M}_{*}}\gtrsim 10^{10.4}\ \rm{M}_{\odot}$ today; however,
$\rm{R_{Acc}}$ continues to be lower than the average SFR, especially for
galaxies today with stellar masses similar to or larger than that of the Milky
Way.

###  3.3 The Build-up of Baryons in Main-Sequence Galaxies over the past
$\approx 9$ Gyr

Figure 3: The ratio of the baryonic mass of main-sequence galaxies at
$z\approx 1.3$ to that at $z=0$, as a function of their stellar mass at
$z\approx 0$. The black dot indicates the position on the relation of a galaxy
with the stellar mass of the Milky Way today. The dashed horizontal line is at
a ratio of 1, indicating no change in the baryonic mass between $z=1.3$ and
$z=0$. It is clear that the highest-mass (${\rm{M}_{*}}\gtrsim
10^{10.5}~{}\rm{M}_{\odot}$ at $z\approx 1.3$) galaxies had the bulk of their
baryonic mass in place by $z\approx 1.3$; the baryonic mass of these galaxies
has declined over the past $\approx 9$ Gyr. Conversely, the baryonic mass of
galaxies with stellar mass ${\rm{M}_{*}}\lesssim 10^{10}~{}\rm{M}_{\odot}$ at
$z\approx 1.3$ has increased since $z\approx 1.3$, by approximately $30\%$ for
galaxies with ${\rm{M}_{*}}\approx 10^{9}~{}\rm{M}_{\odot}$ at $z\approx 1.3$.

Our formalism can be used to track the build-up of baryons in galaxies333We
emphasize that this does not include the baryons in the CGM. from $z\approx
1.3$ to the present epoch. For this, we combine the estimates of the stellar
mass, molecular gas mass, and atomic gas mass to estimate the baryonic mass of
main-sequence galaxies at $z\approx 1.3$, $z\approx 1$, and $z\approx 0$.
Again, we emphasize that we assume that the ionized gas mass within galaxies
can be neglected relative to the stellar and neutral gas masses, and that we
are not considering the CGM here. Figure 3 shows, as a function of the stellar
mass of galaxies today, the fraction of the baryonic mass in individual
galaxies today that was already within the galaxy by $z\approx 1.3$. It can be
seen that most of the baryonic mass in star-forming galaxies with
${\rm{M}_{*}}\gtrsim 10^{9}~{}\rm{M}_{\odot}$ at $z\approx 1.3$ was already in
place by $z\approx 1.3$. Indeed, Figure 3 indicates that the baryonic mass in
the disk of the most massive galaxies, with ${\rm{M}_{*}}\gtrsim
10^{10.4}~{}\rm{M}_{\odot}$ at $z\approx 1.3$, has actually _declined over the
past $\approx 9$ Gyr_, by $\approx 20$% in galaxies with ${\rm{M}_{*}}\gtrsim
10^{11}~{}\rm{M}_{\odot}$ at $z\approx 1.3$. The bulk of the evolution in
massive galaxies over the last $\approx 9$ Gyr thus appears to involve the
conversion of the baryonic material from atomic gas to molecular gas, and
thence to stars, rather than acquisition of fresh gas from the CGM.
Conversely, gas accretion from the CGM continues to play an important role in
the evolution of low stellar-mass galaxies, with ${\rm{M}_{*}}\lesssim
10^{10}~{}\rm{M}_{\odot}$ at $z\approx 1.3$. Such galaxies have acquired
significant amounts of baryonic material from the CGM over the past $\approx
9$ Gyr, with the largest fractional increase in the baryonic mass ($\approx
30$%) taking place in the lowest stellar-mass galaxies, with
${\rm{M}_{*}}\approx 10^{9}~{}\rm{M}_{\odot}$ at $z\approx 1.3$

##  4 Discussion

Our results show that the net average gas accretion rate in main-sequence
galaxies is lower than the average SFR in both redshift intervals, $z=1.3$ to
$z=1.0$, and $z=1.0$ to $z=0$. This implies that the gas accretion rate does
not keep pace with the SFR in main-sequence galaxies at $z\lesssim 1.3$: the
star-formation activity in such galaxies is driven by the Hi that is present
in their disks from earlier times. Our results thus suggest that the net gas
accretion was significantly more efficient at high redshifts, $z\gtrsim 1.3$,
causing the build-up of Hi in the disks of main-sequence galaxies. It is
plausible that this efficient accretion is via the cold mode, with gas flowing
along filaments onto the disks of galaxies (e.g. Kereš et al., 2005; Dekel et
al., 2009). In this scenario, the large Hi reservoirs resulted in a very
efficient conversion of Hi to ${\rm H_{2}}$ at $z\approx 1-3$, consistent with
our estimate of a high $\rm{R_{Mol}}$, comparable to $\rm{R_{SF}}$, over
$z\approx 1.3-1.0$. This efficient Hi-to-${\rm H_{2}}$ conversion, in turn,
triggers the high star-formation activity in galaxies at these redshifts,
resulting in the observed era of high star-formation activity over $z\approx
1-3$, the epoch of galaxy assembly (Madau & Dickinson, 2014)

At later times, $z\lesssim 1.3$, a combination of physical processes are
likely to have led to the observed low net gas accretion rate of Figure 2[B].
These include a transition to inefficient hot-mode accretion in massive
galaxies, due to the increase in the masses of their dark matter halos (e.g.
Dekel & Birnboim, 2006), an increased efficiency of stellar- and AGN-driven
outflows in massive galaxies that remove gas from galaxies (e.g. Weiner et
al., 2009; Veilleux et al., 2020), interactions between neighbouring galaxies
that disrupt the smooth cold-mode filamentary flows, etc. The low net gas
accretion rate at $z\lesssim 1.3$, combined with the very high SFR and high
molecular gas formation rate, implies that the Hi reservoir is consumed on a
very short timescale of $\lesssim 1$ Gyr (Chowdhury et al., 2022b). This
results in a decline in the average Hi content of galaxies by $z\approx 1$
(Chowdhury et al., 2022d), and in a lower average molecular gas formation rate
at later times, as seen in Figure 2[B] at $z<1$. The resulting decline in the
amount of molecular gas available for star formation causes a decline in the
star-formation activity in individual galaxies, especially in the more massive
systems that dominate the SFR density of the Universe at $z\gtrsim 1$. The low
net gas accretion rate at $z\lesssim 1.3$ is thus the primary cause of the
observed decline in the SFR density of the Universe at $z\lesssim 1$. We
emphasize that the above scenario is consistent with the measurements of the
large Hi and ${\rm H_{2}}$ reservoirs of main-sequence galaxies at $z\gtrsim
1$, and the decline in the neutral-gas content of these galaxies at later
times (Chowdhury et al., 2022d, b, a; Tacconi et al., 2020).

Our estimates of $\rm{R_{SF}}$, $\rm{R_{Mol}}$, and $\rm{R_{Acc}}$ over
$z\approx 0-1$ are the time-averaged rates over a relatively long, $\approx 8$
Gyr, interval. This implies that there could be significant evolution in each
of these quantities between $z\approx 1$ and $z\approx 0$ that is not captured
by the averages in Figure 2[B]. Indeed, Bera et al. (2023a) apply a similar
formalism to find that galaxies over the last $\approx 4$ Gyr, from $z\approx
0.35$ to $z\approx 0$, have $\rm{R_{Acc}}\approx\rm{R_{SF}}$, while
$\rm{R_{Mol}}$ is lower than the other two444We note, for ease of comparison,
that Figure 4 of Bera et al. (2023a) shows the rate of change of the molecular
gas content of galaxies, i.e. $({\rm{M_{Mol,f}}}-{\rm{M_{Mol,i}}})/\Delta
t\equiv(\rm{R_{Mol}}-\rm{R_{SF}})$, rather than $\rm{R_{Mol}}$ as defined in
Equation 4; similar conclusions are reached using either metric.. They hence
conclude that the decline in the star-formation activity over the past
$\approx 4$ Gyr is primarily due to inefficient conversion of Hi to ${\rm
H_{2}}$ (Bera et al., 2023a). Combining the two results, it then appears that
for the first few Gyrs after $z\approx 1$, the star-formation activity
continues to be limited by the accretion of Hi onto galaxies, but at later
times, even though the relatively low SFRs of galaxies are now comparable to
their gas accretion rates, the star-formation activity continues to decline,
but now primarily due to the low molecular gas formation rate. However, we
caution that. the ranges of stellar mass explored in the study of Bera et al.
(2023a) and this work are very different, with the typical stellar masses here
being far higher than those in Bera et al. (2023a). In addition, as pointed
out by Bera et al. (2023a), the small cosmic volume covered in their study
implies that their $\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ scaling relation may
be affected by cosmic variance.

##  5 Summary

In this Letter, we have combined our measurement of the
$\rm{M_{H{\textsc{i}}}}-{\rm{M}_{*}}$ scaling relation in main-sequence
galaxies at $z\approx 0.74-1.45$ from the GMRT-CAT$z$1 survey, with estimates
of the main-sequence relation and the ${\rm{M_{Mol}}}-{\rm{M}_{*}}$ relation,
and the assumption of the continuity of main-sequence evolution, to determine
the net average gas accretion rate and the molecular gas formation rate in
main-sequence galaxies over two interesting intervals, $z=1.3$ to $z=1.0$, at
the end of the epoch of galaxy assembly, and $z=1.0$ to $z=0$, the period when
the SFR density of the Universe declines by an order of magnitude. We find
that the average molecular gas formation rate is comparable to the average SFR
in the former interval, but that the net gas accretion rate is significantly
lower than both $\rm{R_{Mol}}$ and $\rm{R_{SF}}$. The low net gas accretion
rate results in a depletion of the Hi reservoir in massive galaxies by
$z\approx 1$, and causes the decline in the SFR density of the Universe at
lower redshifts. In the lower-redshift interval, both $\rm{R_{Acc}}$ and
$\rm{R_{Mol}}$ are lower than $\rm{R_{SF}}$ for all stellar masses; this
results in the consumption of both the Hi and ${\rm H_{2}}$ reservoirs, and a
continuing decline in the SFR density of the Universe down to the present
epoch. We also find that most of the baryonic contents of massive main-
sequence galaxies, with ${\rm{M}_{*}}\gtrsim 10^{10.5}~{}\rm{M}_{\odot}$ at
$z\approx 1.3$, was already in place $\approx 9$ Gyr ago, while low-mass
galaxies, with ${\rm{M}_{*}}\lesssim 10^{10}~{}\rm{M}_{\odot}$ at $z\approx
1.3$, have built up approximately $30\%$ of their baryonic content over this
period.

We thank an anonymous referee for a report that improved the clarity of this
paper. We also 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. AC and NK thank Apurba Bera for many
discussions about Hi scaling relations and gas accretion that have contributed
to this work. 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
  * Barnes et al. (2001) Barnes, D. G., Staveley-Smith, L., de Blok, W. J. G., et al. 2001, MNRAS, 322, 486, doi: 10.1046/j.1365-8711.2001.04102.x
  * 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. (2022) —. 2022, ApJ, 940, L10, doi: 10.3847/2041-8213/ac9d32
  * Bera et al. (2023a) —. 2023a, ApJ, 956, L15, doi: 10.3847/2041-8213/acf71a
  * Bera et al. (2023b) —. 2023b, ApJ, 950, L18, doi: 10.3847/2041-8213/acd0b3
  * Bothwell et al. (2013) Bothwell, M. S., Maiolino, R., Kennicutt, R., et al. 2013, MNRAS, 433, 1425, doi: 10.1093/mnras/stt817
  * Bouché et al. (2010) Bouché, N., Dekel, A., Genzel, R., et al. 2010, ApJ, 718, 1001, doi: 10.1088/0004-637X/718/2/1001
  * 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 et al. (2022a) Chowdhury, A., Kanekar, N., & Chengalur, J. N. 2022a, ApJ, 935, L5, doi: 10.3847/2041-8213/ac8150
  * Chowdhury et al. (2022b) —. 2022b, ApJ, 937, 103, doi: 10.3847/1538-4357/ac7d52
  * Chowdhury et al. (2022c) —. 2022c, ApJ, 941, L6, doi: 10.3847/2041-8213/ac9d8a
  * Chowdhury et al. (2022d) —. 2022d, ApJ, 931, L34, doi: 10.3847/2041-8213/ac6de7
  * 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
  * Daddi et al. (2010) Daddi, E., Bournaud, F., Walter, F., et al. 2010, ApJ, 713, 686, doi: 10.1088/0004-637X/713/1/686
  * Dekel & Birnboim (2006) Dekel, A., & Birnboim, Y. 2006, MNRAS, 368, 2, doi: 10.1111/j.1365-2966.2006.10145.x
  * Dekel et al. (2009) Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Nature, 457, 451, doi: 10.1038/nature07648
  * Dénes et al. (2014) Dénes, H., Kilborn, V. A., & Koribalski, B. S. 2014, MNRAS, 444, 667
  * Giovanelli et al. (2005) Giovanelli, R., Haynes, M. P., Kent, B. R., et al. 2005, AJ, 130, 2598, doi: 10.1086/497431
  * Huang et al. (2012) Huang, S., Haynes, M. P., Giovanelli, R., & Brinchmann, J. 2012, ApJ, 756, 113, doi: 10.1088/0004-637X/756/2/113
  * Kalberla & Kerp (2009) Kalberla, P. M. W., & Kerp, J. 2009, ARA&A, 47, 27, doi: 10.1146/annurev-astro-082708-101823
  * Kereš et al. (2005) Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2, doi: 10.1111/j.1365-2966.2005.09451.x
  * Leitner (2012) Leitner, S. N. 2012, ApJ, 745, 149, doi: 10.1088/0004-637X/745/2/149
  * Leitner & Kravtsov (2011) Leitner, S. N., & Kravtsov, A. V. 2011, ApJ, 734, 48, doi: 10.1088/0004-637X/734/1/48
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * Møller et al. (2013) Møller, P., Fynbo, J. P. U., Ledoux, C., & Nilsson, K. K. 2013, MNRAS, 430, 2680, doi: 10.1093/mnras/stt067
  * Noeske et al. (2007) Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJ, 660, L43, doi: 10.1086/517926
  * Parkash et al. (2018) Parkash, V., Brown, M. J. I., Jarrett, T. H., & Bonne, N. J. 2018, ApJ, 864, 40, doi: 10.3847/1538-4357/aad3b9
  * Peng et al. (2010) Peng, Y.-j., Lilly, S. J., Kovač, K., et al. 2010, ApJ, 721, 193, doi: 10.1088/0004-637X/721/1/193
  * Péroux & Howk (2020) Péroux, C., & Howk, J. C. 2020, ARA&A, 58, 363, doi: 10.1146/annurev-astro-021820-120014
  * Popesso et al. (2022) Popesso, P., Concas, A., Cresci, G., et al. 2022, MNRAS, doi: 10.1093/mnras/stac3214
  * Renzini (2009) Renzini, A. 2009, MNRAS, 398, L58, doi: 10.1111/j.1745-3933.2009.00710.x
  * Rodighiero et al. (2011) Rodighiero, G., Daddi, E., Baronchelli, I., et al. 2011, ApJ, 739, L40, doi: 10.1088/2041-8205/739/2/L40
  * Saintonge & Catinella (2022) Saintonge, A., & Catinella, B. 2022, arXiv e-prints, arXiv:2202.00690.  https://arxiv.org/abs/2202.00690
  * Saintonge et al. (2017) Saintonge, A., Catinella, B., Tacconi, L. J., et al. 2017, ApJS, 233, 22, doi: 10.3847/1538-4365/aa97e0
  * Santini et al. (2014) Santini, P., Maiolino, R., Magnelli, B., et al. 2014, A&A, 562, A30, doi: 10.1051/0004-6361/201322835
  * Scoville et al. (2016) Scoville, N., Sheth, K., Aussel, H., et al. 2016, ApJ, 820, 83, doi: 10.3847/0004-637X/820/2/83
  * Scoville et al. (2017) Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150, doi: 10.3847/1538-4357/aa61a0
  * Speagle et al. (2014) Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15, doi: 10.1088/0067-0049/214/2/15
  * Tacconi et al. (2020) Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157, doi: 10.1146/annurev-astro-082812-141034
  * 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
  * Veilleux et al. (2020) Veilleux, S., Maiolino, R., Bolatto, A. D., & Aalto, S. 2020, A&A Rev., 28, 2, doi: 10.1007/s00159-019-0121-9
  * Walter et al. (2020) Walter, F., Carilli, C., Neeleman, M., et al. 2020, ApJ, 902, 111, doi: 10.3847/1538-4357/abb82e
  * Weiner et al. (2009) Weiner, B. J., Coil, A. L., Prochaska, J. X., et al. 2009, ApJ, 692, 187, doi: 10.1088/0004-637X/692/1/187
  * 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
  * Zwaan (2000) Zwaan, M. A. 2000, PhD thesis, Ph.D. Thesis, Groningen: Rijksuniversiteit, 2000


# Learned Spike Encoding of the Channel Response for Low-Power Environment
Sensing

Eleonora Cicciarella∗†, Riccardo Mazzieri†, Jacopo Pegoraro†, Michele Rossi§†
∗Corresponding author<EMAIL_ADDRESS>†These authors are
with the University of Padova, Department of Information Engineering. §This
author is with the University of Padova, Department of Mathematics “Tullio
Levi-Civita”. This work was partially supported by the European Union under
the Italian National Recovery and Resilience Plan (NRRP) of NextGenerationEU,
partnership on “Telecommunications of the Future” (PE0000001 - program
“RESTART”).

###### Abstract

Radio Frequency (RF) sensing holds the potential for enabling pervasive
monitoring applications. However, modern sensing algorithms imply complex
operations, which clash with the energy-constrained nature of edge sensing
devices. This calls for the development of new processing and learning
techniques that can strike a suitable balance between performance and energy
efficiency. Spiking Neural Networks (SNNs) have recently emerged as an energy-
efficient alternative to conventional neural networks for edge computing
applications. They process information in the form of sparse binary spike
trains, thus potentially reducing energy consumption by several orders of
magnitude. Their fruitful use for RF signal processing critically depends on
the representation of RF signals in the form of spike signals. We underline
that existing spike encoding algorithms to do so generally produce inaccurate
signal representations and dense (i.e., inefficient) spike trains.

In this work, we propose a lightweight neural architecture that learns a
tailored spike encoding representations of RF channel responses by jointly
reconstructing the input and its spectral content. By leveraging a tunable
regularization term, our approach enables fine-grained control over the
performance-energy trade-off of the system. Our numerical results show that
the proposed method outperforms existing encoding algorithms in terms of
reconstruction error and sparsity of the obtained spike encodings.

###### Index Terms:

spiking neural networks, edge computing, radio frequency sensing, radar
sensing, energy efficiency.

##  I Introduction

Radio Frequency (RF) sensing has gathered a huge interest as a privacy-
preserving and ubiquitous alternative to cameras for environment monitoring,
pervasive healthcare, and touchless human-machine interaction, among others
[1, 2]. Depending on the adopted technology, existing RF sensing systems can
use dedicated commercial radar sensors adopting Frequency-Modulated
Continuous-Wave (FMCW) or Ultra-Wideband (UWB) signals [2], reuse the signals
transmitted by existing Wi-Fi routers or by cellular base stations [3] (termed
Integrated Sensing and Communication (ISAC)). RF sensing entails the
estimation of the wireless channel, as the latter embeds information about the
location and movement of objects and people in the physical surroundings. Due
to the complexity of such channel, Deep Learning (DL) solutions are commonly
employed to enable advanced sensing applications such as target recognition,
movement analysis, and human activity recognition [4, 3]. However, DL is well-
known to be energy-hungry and, in turn, unsuited for implementation on
resource-constrained edge devices and routers, where the RF data are collected
[5]. This is further exacerbated by the extremely high sampling rates of
modern analog-to-digital converters, which cause RF sensors to generate a huge
amount of data compared to other sensing modalities. Hence, a widespread
adoption of RF sensing, and ISAC in particular, raises concerns regarding the
large carbon footprint of the resulting sensor networks.

A possible solution to tackle the energy efficiency problem of RF sensing at
the network edge has been identified in replacing standard DL architectures
with Spiking Neural Networks [6]. SNNs closely mimic the natural functioning
of the brain, which is a strikingly energy-efficient processor operating
within $\sim 12-20$ W of power. They encode information into event-based spike
trains, which can be represented as sparse time-domain binary signals [6].
Their energy consumption is directly related to the number of processed spikes
and can be orders of magnitude smaller than than of standard GPU accelerators
[7]. As such, the way in which the input (channel) information is encoded into
spikes is key in determining the trade-off between performance on the learning
task and sparsity of the resulting spike train, which is tied to the energy
consumption. However, most available works on SNNs for RF sensing overlook the
spike encoding step and rely on simple, general purpose techniques [8, 9].
Specifically, these approaches perform standard radar signal processing
operations to estimate wireless channel features and then feed them to
threshold- or moving-window-based encoding methods [10]. These are not
tailored for RF sensing data, and offer limited control on the sparsity of the
spike trains. Other preliminary works use the raw RF samples as the SNN input
[11, 12], thus avoiding the use of encoding algorithms. Notably, [12] proposes
an SNN-based ISAC system using short UWB pulses, which are naturally
represented as spikes after receiver sampling. These solutions achieve a
better integration of the SNN into the RF signal processing chain, but do not
address the minimization of the energy consumption.

In this work, we propose a new approach to gain control over the performance-
energy trade-off in RF sensing with SNNs. Our key insight is that channel
representations in RF sensing can be cast under a sum of sinusoids model, in
which the main features of interest are the frequency components and their
amplitudes. We propose a lightweight neural architecture that learns a
tailored spike encoding for RF channel responses (treated as input signals),
preserving their spectral content. This architecture is trained to jointly
optimize the fidelity to the original spectrum (hence the sensing performance)
and the energy consumption, by promoting the sparsity of the encoded spike
train.

The main contributions of this work are summarized next.

1) We identify a sum of sinusoids signal model to represent the channel
response in most RF sensing tasks of interest (like, e.g., range, Doppler, and
Angle of Arrival (AoA) estimation). Hence, we leverage this model to formulate
a spike encoding learning mechanism that achieves a desirable trade-off
between fidelity to the original channel properties and energy consumption.
This is achieved by minimizing the number of spikes in the encoding.

2) We design and validate a Convolutional Autoencoder (CAE) that leverages a
Straight-Through Estimator (STE) to obtain a spike train feature encoding. Two
output branches are used. One reconstructs the original channel
representation, pushing the spike train to accurately represent the input. The
other is a SNN performing a regression on the frequency components of the
channel representation. To the best of our knowledge, this is the first
encoding method to learn spike patterns from RF sensing data.

3) We train the proposed network on a synthetic dataset, showing its excellent
capabilities to effectively encode RF channel responses. Specifically, our
method encodes channel features (and its spectral content in particular) by
achieving spike trains that are more than $70\%$ sparser than those obtained
from existing encoding algorithms.

##  II Background and signal model

In both radar and ISAC, sensing is performed by estimating the wireless
propagation channel, which embeds information about the location and movement
of objects and people in the surroundings. The key motivation behind our
learned encoding design and signal model (see Section II-A) lies in the common
representation of the channel response for both radar and ISAC processing,
which takes a sum of sinusoids form [13, 4, 3]. For example, FMCW radars
estimate distance, Doppler speed, and AoA of the sensing targets by processing
the echo of a probe waveform, obtaining the so-called intermediate frequency
signal [13]. This has a complex sinusoidal shape in the delay, Doppler, and
spatial domains, which allows sensing parameters estimation via 3D- Discrete
Fourier Transform (DFT). Although different waveforms and processing are
employed, also ISAC systems based on Orthogonal Frequency Division
Multiplexing (OFDM) or single-carrier transmission obtain similar channel
representations, as detailed, e.g., in [4, 3]. Hence, to design low-power RF
sensing technology based on SNNs it is fundamentally important to efficiently
encode sum of sinusoids channels into spike trains that: (i) preserve the
frequency content of the channel response, as this contains the sensing
parameters of interest and (ii) contain as few spikes as possible, to minimize
the energy consumption.

In the following, we introduce the signal model representing the channel
response as a sum of sinusoids and provide useful background on SNNs.

###  II-A Signal model

As anticipated above, a sum of sinusoids model is considered to represent the
channel response in the delay, Doppler, or spatial domains. The number of
sinusoidal components is denoted by $M$, which correspond to the channel
reflectors. Moreover, $a_{m},\phi_{m},f_{m}$ respectively represent the
amplitude (accounting for the propagation and reflection losses), the phase,
and the frequency of the $m$-th component. The sampling time is denoted by $T$
and a processing window of $K$ subsequent samples is considered, indexed by
$k=0,\dots,K-1$. The signal model is written as

$x[k]\triangleq x(kT)=\sum_{m=1}^{M}a_{m}e^{j(2\pi f_{m}kT+\phi_{m})}+w(kT),$
(1)

where $w$ is a complex Gaussian noise with variance $\sigma_{w}^{2}$. We
stress that the above model has huge applicability in RF sensing, as it can
represent the channel response in terms of reflectors distance, Doppler speed,
or AoA of the reflector, depending on the domain in which the signal is
sampled (fast-time, slow-time, or spatial) [13, 4, 14]. In each case, the
frequency $f_{m}$ assumes a different meaning and it is mapped onto a
propagation delay, a Doppler speed, or an angle.

###  II-B Spiking Neural Networks and spike encoding methods

####  II-B1 Spiking neuron model

We employ the well-known Leaky Integrate and Fire (LIF) neuron [15], which is
characterized by an internal state called membrane potential. Each neuron acts
as a leaky integrator of the input signal (called current) and fires an output
spike whenever its membrane potential exceeds a threshold. In the absence of
input, the membrane potential decays exponentially to zero with rate $\beta$.
After generating a spike, the value of the threshold is subtracted from the
membrane potential. SNNs are constituted by a network of LIF neurons,
typically organized into layers. The output spikes of layer $\ell$ are
multiplied by learnable parameters called synaptic weights before being fed to
layer $\ell+1$. This enables optimizing the SNN using a surrogate gradient
algorithm, which approximates the standard backpropagation used in standard
neural networks [6]. Specifically, the non-differentiable threshold function
that determines the emission of the spikes is approximated by a differentiable
surrogate function during the backward pass, while in the forward pass the
original thresholding is applied.

####  II-B2 Spike encoding

The vast majority of available data for practical applications is not
collected in the form of spike trains, as in the case of RF sensing. To
efficiently process this kind of data with SNNs, the need to develop spike
encoding techniques arises [16]. The most common spike encoding approach for
sequential data is temporal contrast encoding [10]. This class of methods
keeps track of temporal changes in the signal and produces positive or
negative spikes according to a thresholding criterion. In this work, we
consider the following spike encoding algorithms:

  * •

Threshold-Based Representation (TBR) generates a positive or negative spike
when the absolute difference between two consecutive signal values (first-
order difference) is higher than the threshold, and the sign of the generated
spike is the sign of the difference. Given an input sequence, $x[k]$, the
value of the threshold is set to $\mu(x^{\prime})+\delta\sigma(x^{\prime})$,
where $x^{\prime}$ is the sequence of first-order differences of $x$, $\mu$
and $\sigma$ denote the mean and standard deviation, respectively, and
$\delta$ is a parameter of the algorithm.

  * •

Step-Forward (SF) compares the difference between consecutive values with a
dynamic baseline. When the absolute difference deviates from the baseline by
more than the threshold, a positive or negative spike is generated, and the
baseline is updated to the current signal value.

  * •

Moving-Window (MW) is equivalent to SF, but the value of the baseline is set
equal to the moving average of the signal at each step. Besides the threshold,
this algorithm requires choosing the width of the sliding window used to
compute the moving average.

##  III Learned spike encoding

Next, we describe the proposed learning-based encoding method for RF sensing
signals, shown in Fig. 1. The neural network is composed of three main blocks:
an encoder, a decoder, and an SNN, as introduced in the following.

  1. 1.

After a preprocessing step, the encoder block extracts salient features of the
input channel response (see Section III-B), which are then fed to a threshold
function to produce the spike encoding (Section III-C).

  2. 2.

The decoder block reconstructs the input channel response by approximately
inverting the encoding (Section III-D). The network comprising the encoder and
decoder blocks is referred to as CAE in the following.

  3. 3.

The SNN block (Section III-E) takes the spike encoding produced by the encoder
as input and performs regression to estimate: (i) the frequency components of
the input, $f_{m}$, and (ii) their amplitudes, $a_{m}$. By adding this block,
we push the encoder to output spike encodings that preserve the spectral
content of the channel.

The whole network is jointly trained as described in Section III-F. A detailed
description of each block follows.

###  III-A Channel response preprocessing

The input channel response is modeled as in Eq. (1) and processed in windows
of $K$ complex-valued samples, i.e., $[x[0],\dots,x[K-1]]^{T}$. Before feeding
it to the encoder, each processing window is transformed into a real-valued
tensor, $\boldsymbol{X}$, with shape $(2,K)$, where the first dimension is the
number of convolution channels of the CAE. We call $X_{c}[k]$ the $c$-th
convolution channel of $\boldsymbol{X}$ at time $k$, with $c\in\\{1,2\\}$.

Figure 1: Block diagram of the proposed architecture.

###  III-B Encoder

The encoder takes as input the channel response processing window
$\boldsymbol{X}\in\mathbb{R}^{2\times K}$ and outputs a feature representation
$\boldsymbol{Y}\in\mathbb{R}^{2\times K}$. Let
$E_{\xi}(\cdot):\mathbb{R}^{2\times K}\rightarrow\mathbb{R}^{2\times K}$ be
the parameterized function representing the encoder, where $\xi$ is the set of
weights of the encoder. Then, $\boldsymbol{Y}=E_{\xi}(\boldsymbol{X})$.

In the proposed architecture, the encoder employs $3$ one-dimensional (1D)
convolutional layers having $128$, $128$, and $2$ feature maps, respectively,
and filters of length $7$ each. A stride of $1$ and zero-padding are used to
preserve the input dimension. Finally, batch normalization and the hyperbolic
tangent function are applied at each hidden layer, except for the last layer
which has no activation function.

###  III-C Spike encoding

From the feature representation $\boldsymbol{Y}$, we obtain the spike encoding
$\boldsymbol{Z}$, which lies in the space $\mathcal{Z}=\\{-1,0,1\\}^{2\times
K}$. This is done by applying an element-wise hard-threshold function,
$\mathcal{H}_{\tau}$, to $\boldsymbol{Y}$. Using $Y_{c}[k],Z_{c}[k]$ to
respectively denote the $c$-th convolution channel of $\boldsymbol{Y}$ and
$\boldsymbol{Z}$ at time $k$, we have

$Z_{c}[k]=\mathcal{H}_{\tau}\left(Y_{c}[k]\right)=\begin{cases}\mathrm{sign}(Y_{c}[k])\,\,&\mathrm{if}\,|Y_{c}[k]|\geq\tau,\\\
\,\,0&\mathrm{otherwise,}\end{cases}$ (2)

where $\tau>0$ is a threshold parameter. The two channels of $\boldsymbol{Z}$
represent the encoding of the real and imaginary part of the channel response,
respectively. It is worth noting that, although the encoder does not reduce
the dimensionality of the input, it compresses it by discretizing the signal
values in $\mathcal{Z}$. This leads to a loss of information that prevents the
CAE from simply learning the identity function.

Remark on the use of $\mathcal{H}_{\tau}$: The use of the threshold function
$\mathcal{H}_{\tau}$ prevents a direct application of the backpropagation
algorithm to train the network. Indeed, the derivative of
$\mathcal{H}_{\tau}(Y_{c}[k])$ evaluates to zero everywhere, except in
$Y_{c}[k]=\pm\tau$, where it is not defined. This causes the gradients to
vanish during backpropagation. To address this problem, we employ an STE [17],
which consists in replacing $\mathcal{H}_{\tau}$ with the hard hyperbolic
tangent, which clamps the input value to the range $[-1,1]$, during the
backward propagation of the gradients.

###  III-D Decoder

We represent the decoder with a parameterized function
$D_{\psi}(\cdot):\mathcal{Z}\rightarrow\mathbb{R}^{2\times K}$, where $\psi$
is the set of weights of the decoder. The reconstructed channel response is
$\hat{\boldsymbol{X}}=D_{\psi}(\boldsymbol{Z})$. The decoder is composed of
$3$ _transposed_ 1D convolutional layers with $128,128$, and $2$ filters,
respectively. Filter size, stride, padding, batch normalization, and
activation function are the same as those used in the encoder, except for the
last layer which implements the sigmoid function.

###  III-E Spiking Neural Network

The SNN processes the input spike encoding, $\boldsymbol{Z}$, to regress the
frequencies and amplitudes of the sinusoidal components in the input signal,
denoted by $\hat{\bm{f}}=[\hat{f}_{1},\dots,\hat{f}_{M}]^{T}\in\mathbb{R}^{M}$
and
$\hat{\boldsymbol{a}}=[\hat{a}_{1},\dots,\hat{a}_{M}]^{T}\in\mathbb{R}^{M}$,
respectively. Denoting by $\gamma$ the set of parameters of the SNN, we
represent the network with the function
$S_{\gamma}:\mathcal{Z}\rightarrow\mathbb{R}^{M}\times\mathbb{R}^{M}$, such
that $(\hat{\bm{f}},\hat{\boldsymbol{a}})=S_{\gamma}(\boldsymbol{Z})$. To do
so, the SNN interprets $\boldsymbol{Z}$ as a $K$ timesteps-long sequence of
$2$-dimensional vectors, each containing the real and imaginary spike
components. Specifically, we define the input vector to the SNN at time $k$ as
$\boldsymbol{z}[k]\triangleq[Z_{1}[k],Z_{2}[k]]^{T}$.

The SNN is composed of LIF layers. Each layer takes as input the spike train
produced by the previous layer and behaves as a leaky integrator until a
firing threshold is reached. At layer $\ell$, denote by
$\boldsymbol{u}_{\ell}[k]$ the vector of membrane potentials at time $k$, by
$\boldsymbol{\theta}_{\ell}$ the vector containing the firing thresholds of
each neuron, by $\boldsymbol{W}_{\ell}$ the learnable weights matrix, and by
$\boldsymbol{\beta}_{\ell}$ the vector of membrane potential decay constants
[6]. The input-output equation of the $\ell$-th LIF layer is

$\boldsymbol{u}_{\ell}[k]=\boldsymbol{\beta}_{\ell}\boldsymbol{u}_{\ell}[k-1]+\boldsymbol{W}_{\ell}\boldsymbol{s}_{\ell-1}[k]-\boldsymbol{\theta}_{\ell}\boldsymbol{s}_{\ell}[k-1],$
(3)

where vector $\boldsymbol{s}_{\ell}[k]$ is the output spike vector of layer
$\ell$ at time $k$. Note that for the first layer we set
$\boldsymbol{s}_{0}[k]=\boldsymbol{z}[k]$.

Our SNN architecture comprises an input layer with $128$ neurons, a hidden
layer with $64$ neurons, and two parallel output layers, which we denote by
$\mathcal{O}_{f}$ and $\mathcal{O}_{a}$, with $5$ neurons each. The output
layers are specifically designed to reconstruct the channel frequencies and
the amplitudes, respectively. The output spikes from the hidden layer are
processed in parallel by $\mathcal{O}_{f}$ and $\mathcal{O}_{a}$ and then the
membrane potential of each output layer is passed through the sigmoid function
to obtain $\hat{\bm{f}}$ and $\hat{\bm{a}}$, respectively.

The membrane potential decay rates, $\boldsymbol{\beta}_{\ell}$, and the
firing thresholds, $\boldsymbol{\theta}_{\ell}$, are learned during training,
hence the set of learnable parameters of the SNN, $\gamma$, contains
$\boldsymbol{W}_{\ell},\boldsymbol{\beta}_{\ell},\boldsymbol{\theta}_{\ell}$
for the four layers.

###  III-F Network training and hyperparameter optimization

In this section, we describe the training process of the learned spike
encoding network.

####  III-F1 Loss function and sparsity regularization

The network is trained end-to-end through the minimization of a loss function
$\mathcal{L}$ which depends on the outputs of the decoder and the SNN, as well
as on the spike encoding itself. Specifically, let $\boldsymbol{f}$ and
$\boldsymbol{a}$ be the vectors of the true frequency and amplitude components
of the channel response, respectively, with components $f_{m},a_{m}$. We
define $3$ reconstruction loss terms using the Mean Squared Error (MSE) as

$\mathcal{L}_{1}\triangleq\mathrm{MSE}(\boldsymbol{X},\hat{\boldsymbol{X}})=\frac{1}{2K}\sum_{c=1}^{2}\sum_{k=0}^{K-1}(X_{c}[k]-\hat{X}_{c}[k])^{2},$
(4)

$\mathcal{L}_{2}\triangleq\mathrm{MSE}(\boldsymbol{f},\hat{\bm{f}})=\sum_{m=1}^{M}(f_{m}-\hat{f}_{m})^{2}/M$,
and
$\mathcal{L}_{3}\triangleq\mathrm{MSE}(\boldsymbol{a},\hat{\boldsymbol{a}})=\sum_{m=1}^{M}(a_{m}-\hat{a}_{m})^{2}/M$.
The frequencies, amplitudes, and values of $\boldsymbol{X}$ are normalized in
the range $[0,1]$, as usual practice with neural networks. Note that this
implies that all the error terms also lie in $[0,1]$.

Additionally, a regularization term $\Omega(\boldsymbol{Z})$, which depends on
the number of spikes, is added to the loss to penalize non-sparse spike
trains. The use of this regularization is motivated by the objective of
obtaining _sparse_ spike trains since the number of processed spikes is tied
to the energy consumption of the system. Recalling that
$\boldsymbol{Z}\in\\{-1,0,1\\}^{2\times K}$, we define

$\Omega(\boldsymbol{Z})\triangleq\frac{1}{2K}\sum_{c=1}^{2}\sum_{k=0}^{K-1}|Z_{c}[k]|,$
(5)

where the absolute value is due to the presence of negative spikes. We use a
regularization coefficient $\lambda\in[0,1]$, that allows us to tune the
importance of the regularization term in the loss and gain control of the
energy content of the spike encoding. Larger values of $\lambda$ lead to
higher sparsity (lower energy consumption), while $\lambda=0$ means no
regularization is applied (higher energy consumption). The total loss function
$\mathcal{L}$ is obtained as
$\mathcal{L}=\mathcal{L}_{1}+\mathcal{L}_{2}+\mathcal{L}_{3}+\lambda\,\Omega(\boldsymbol{Z})$.

####  III-F2 Training

We train the network in batches of $256$ samples, using backpropagation with
the Adam optimizer with learning rate $\eta=10^{-3}$. In computing the
gradients of the SNN parameters, we use the surrogate gradient approach to
avoid the non-differentiability of the LIF layer [6]. Specifically, we adopt
the arctangent approximation of the firing threshold function during the
backward pass. The training process is stopped upon convergence of the loss on
a validation set.

####  III-F3 Hyperparameter optimization

We optimize the network hyperparameters by performing a greedy search, using
the loss on a validation set as the optimization objective. Notably, we obtain
that the best threshold function for the spike encoding layer is $\tau=0.1$.

##  IV Simulation and numerical results

In this section, we compare the results obtained by the proposed architecture,
referred to as Learned Spike Encoding (LSE) in the following, and the encoding
methods from the literature, presented in Section II-B. All the code is
developed in Python, using the PyTorch and the snnTorch frameworks [6].

###  IV-A Synthetic dataset generation

We create a synthetic dataset using the model presented in Section II-A. To
simulate realistic RF channel responses, we used typical parameters of RF
sensing systems for Doppler spectrum estimation in the $60$ GHz unlicensed
band, see e.g. [4]. In this setting, frequencies $f_{m}$ represent Doppler
shifts induced by the motion of objects in the environment. This is just one
of the possible applications of our proposed encoding algorithm, which is
general enough to be also applied to range or AoA estimation.

Following [4], we set $T=0.27$ ms, and $K=64$, while the carrier frequency is
$f_{c}=60$ GHz. This choice of $T$ leads to a maximum observable Doppler
velocity of $v_{\rm max}=c/(4f_{c}T)\approx 4.48$ m/s, where $c$ is the speed
of light.

To generate channel responses with a variable number of sinusoidal components,
we set $M=5$, to be intended as the maximum number of components in the
synthetic data. Then, channel responses with fewer than $M$ components are
obtained by setting to $0$ the corresponding amplitudes $a_{m}$. We generate
$3,000$ windows for each possible number of sinusoidal components in
$\\{1,\dots,M\\}$, for a total of $15,000$. The dataset is randomly split into
training, validation, and test sets with the proportion $75$:$15$:$10$. We
generate each channel response window by applying the following pipeline:

  1. 1.

the frequency $f_{m}$ of each sinusoidal component is randomly generated in
the interval $[0,1/(2T)]$;

  2. 2.

the phase $\phi_{m}$ of each component is randomly generated in the interval
$[0,2\pi)$;

  3. 3.

the amplitudes are randomly generated and normalized with respect to the
maximum one so that it equals $1$;

  4. 4.

the complex Gaussian noise is generated such that the signal-to-noise ratio
(SNR) belongs to the set $\\{5,10,15,20\\}$ dB.

###  IV-B Comparison with state-of-the-art algorithms

In this section, we compare our learned encoding algorithm to the state-of-
the-art methods introduced in Section II-B. Recall that the quality of the
final encoding provided by the TBR, SF, and MW algorithms heavily depends on
the choice of the input parameters. Hence, for each algorithm, we performed a
grid search to find the values of the parameters that minimize
$\mathrm{MSE}(\boldsymbol{X},\hat{\boldsymbol{X}})$. For each algorithm, the
specific search space and final optimal values are reported in Tab. I.

Table I: Optimization of TBR, SF, and MW parameters. Method | Search interval | Step value | Optimal value
---|---|---|---
TBR (factor $\delta$) | $[0.001,0.01]$ | 0.001 | 0.005
SF (threshold) | $[0.05,0.3]$ | 0.05 | 0.2
MW (threshold) | $[0.01,0.6]$ | 0.05 | 0.06
MW (window) | $[2,4]$ | 1 | 3
Table II: Performance comparison (normalized reconstruction errors and average sparsity) between our method (LSE) and threshold-based methods. We report the mean and standard deviation on the test set. Method | Rec. RMSE | DFT mag. RMSE | Sparsity
---|---|---|---
TBR [10] | 0.374 $\pm$ 0.075 | 0.043 $\pm$ 0.010 | 0.015
SF [10] | 0.222 $\pm$ 0.044 | 0.029 $\pm$ 0.009 | 0.433
MW [10] | 0.262 $\pm$ 0.059 | 0.039 $\pm$ 0.011 | 0.202
LSE | 0.133 $\pm$ 0.028 | 0.017 $\pm$ 0.004 | 0.736

Figure 2: Qualitative comparison between the reconstructed channel and the
original one (top), and the learned spike encoding (bottom), where the real
part of the channel (Re) is shown on the left and the imaginary one (Im) on
the right.

####  IV-B1 Reconstruction error

To evaluate the quality of the reconstructed channel, we adopt the per-window
Root Mean Squared Error (RMSE) computed over the test set, obtained as
$\mathrm{RMSE}(\boldsymbol{X},\hat{\boldsymbol{X}})=\sqrt{\mathrm{MSE}(\boldsymbol{X},\hat{\boldsymbol{X}})}$.
All errors are then averaged over all the windows of the test set, to obtain
the mean reconstruction error. As shown in Tab. II, LSE achieves a remarkably
low reconstruction error compared to other methods. Indeed, temporal contrast
techniques are hindered by the use of a fixed threshold to encode the signal,
which prevents the model from capturing fine-grained, dynamic channel
features. An example of a test window reconstructed by LSE, along with the
obtained spike train, is shown in Fig. 2.

####  IV-B2 Spectral components preservation

We evaluate and compare the encoding algorithms on their capability of
preserving the channel spectral content by computing the per-window RMSE of
the reconstructed channel response DFT magnitude with respect to the original
one. As shown in Tab. II, our learned encoding demonstrates superior
performance due to its higher representation power. As an additional result,
we report that our method also achieves good results in the estimation of the
channel frequency components through the SNN, achieving an average per-window
RMSE of $0.1$ for the frequencies and $0.16$ for the amplitudes.

####  IV-B3 Sparsity of the encoding

We define the sparsity of a spike encoding as the fraction of zeros in
$\boldsymbol{Z}$, i.e. $\rho(\boldsymbol{Z})=1-\Omega(\boldsymbol{Z})$.
$\rho(\boldsymbol{Z})$ directly correlates with the efficiency of the SNN used
for processing the spike trains. The sparsity of each method is then obtained
by averaging the sparsity of all the encoded test windows, as reported in Tab.
II. LSE exhibits $70\%$ higher sparsity than the best algorithm from the
literature, which is obtained by SF ($0.73$ vs. $0.43$). This demonstrates its
efficiency in representing key features with minimal redundancy (this result
is obtained using $\lambda=0.2$). Other standard encoding methods, especially
TBR, produce denser spike trains. Fig. 3 highlights the trade-off between
encoding sparsity and reconstruction accuracy. Our LSE offers a substantial
improvement over existing approaches in encoding RF sensing data. Moreover,
the proposed approach allows directly controlling the sparsity-accuracy trade-
off by tuning the regularization parameter $\lambda$. This is shown in Fig. 4,
where we plot the sparsity and the corresponding validation reconstruction MSE
for different values of $\lambda$ (averaged over $14$ training with different
initialization). Higher $\lambda$ leads to higher sparsity (lower energy
consumption) but implies performance degradation (higher MSE).

Figure 3: Reconstruction RMSE for different spike encoding sparsity levels of
the channel response (left) and its DFT magnitude (right).

Figure 4: Sparsity level for different values of $\lambda$ using the learned
spike encoding.

###  IV-C Robustness to noise

Lastly, we analyze the robustness of each method to noise. To this end, we
generated multiple test sets, with $10,000$ samples each, with different SNR
values ($5,10,15$, and $20\leavevmode\nobreak\ \mathrm{dB}$), following the
pipeline outlined in Section II-A. Hence, for each test set and encoding
method, we computed the reconstruction and DFT magnitude errors. The results
are shown in Fig. 5. LSE obtains the lowest errors and improves the most at
high SNR. Notably, LSE also has the lowest standard deviation across different
SNR values, which makes it more reliable.

##  V Concluding remarks

In this work, we have designed and implemented the first neural network that
learns a _tailored_ spike encoding for RF channel responses in RF sensing
applications. To this end, we proposed a CAE to obtain the spike encoding of
the input data, enforcing sparsity by penalizing dense spike trains. In
addition, we used an SNN block to perform a regression on both the frequency
and the amplitude components of the input data, to push the spike train to
preserve the original spectral content. The proposed network is trained end-
to-end on a synthetic dataset. Our results show that it learns to produce a
better spike encoding than existing thresholding-based methods, both in terms
of fidelity to the spectral content of the original signal and sparsity of the
spike trains. While our primary objective is to shift from traditional DL
architectures to SNN for energy efficiency reasons, it is important to
highlight the compactness of the proposed DL-based encoder architecture (fewer
than $120$ K parameters and less than $2$ MB of size in memory), which ensures
a lower inference time, compared to that of the SNN, and allows its deployment
on low-memory edge devices. Moreover, the proposed technique allows
controlling the performance-energy trade-off through a regularization
parameter that trades information content of the spike encoding with its
sparsity, which is directly connected to the network energy consumption. This
provides high flexibility in adapting it to different applications and
hardware implementations.

Future developments of this work include: (i) an evaluation of the trained
encoder on real RF sensing data, (ii) an in-depth analysis of the energy
consumption of the encoding during inference, (iii) the possibility of
integrating the spike encoding step directly into the SNN, and (iv) the
training of the whole network to autonomously recognize the different
sinusoidal components without having to set a maximum number a priori.

Figure 5: Robustness to noise when reconstructing the channel response (left)
and its DFT magnitude (right). Shaded areas represent the standard deviation.

## References

  * [1] J. Liu, H. Liu, Y. Chen, Y. Wang, and C. Wang, “Wireless sensing for human activity: A survey,” IEEE Comm. Surveys & Tutorials, vol. 22, no. 3, pp. 1629–1645, 2019.
  * [2] A. Shastri, N. Valecha, E. Bashirov, H. Tataria, M. Lentmaier, F. Tufvesson, M. Rossi, and P. Casari, “A review of millimeter wave device-based localization and device-free sensing technologies and applications,” IEEE Comm. Surveys & Tutorials, vol. 24, no. 3, pp. 1708–1749, 2022.
  * [3] Y. Ma, G. Zhou, and S. Wang, “WiFi sensing with channel state information: A survey,” ACM Computing Surveys (CSUR), vol. 52, no. 3, pp. 1–36, 2019.
  * [4] J. Pegoraro, J. O. Lacruz, F. Meneghello, E. Bashirov, M. Rossi, and J. Widmer, “RAPID: Retrofitting IEEE $802.11$ay Access Points for Indoor Human Detection and Sensing,” IEEE Transactions on Mobile Computing, pp. 1–18, 2023.
  * [5] S. Lahmer, A. Khoshsirat, M. Rossi, and A. Zanella, “Energy consumption of neural networks on Nvidia edge boards: an empirical model,” in 20th International Symposium on Modeling and Optimization in Mobile, Ad hoc, and Wireless Networks (WiOpt), pp. 365–371, IEEE, 2022.
  * [6] J. K. Eshraghian, M. Ward, E. O. Neftci, X. Wang, G. Lenz, G. Dwivedi, M. Bennamoun, D. S. Jeong, and W. D. Lu, “Training spiking neural networks using lessons from deep learning,” Proc. of the IEEE, 2023.
  * [7] C. Frenkel, M. Lefebvre, J.-D. Legat, and D. Bol, “A 0.086-mm212.7-pJ/SOP 64k-synapse 256-neuron online-learning digital spiking neuromorphic processor in 28-nm CMOS,” IEEE Transactions on biomedical circuits and systems, vol. 13, no. 1, pp. 145–158, 2018.
  * [8] A. Safa, F. Corradi, L. Keuninckx, I. Ocket, A. Bourdoux, F. Catthoor, and G. G. Gielen, “Improving the accuracy of spiking neural networks for radar gesture recognition through preprocessing,” IEEE Transactions on Neural Networks and Learning Systems, 2021.
  * [9] A. Safa, A. Bourdoux, I. Ocket, F. Catthoor, and G. G. Gielen, “On the use of spiking neural networks for ultralow-power radar gesture recognition,” IEEE Microwave and Wireless Components Letters, vol. 32, no. 3, pp. 222–225, 2021.
  * [10] B. Petro, N. Kasabov, and R. M. Kiss, “Selection and optimization of temporal spike encoding methods for spiking neural networks,” IEEE transactions on neural networks and learning systems, vol. 31, no. 2, pp. 358–370, 2019.
  * [11] M. Arsalan, A. Santra, and V. Issakov, “Spiking neural network-based radar gesture recognition system using raw ADC data,” IEEE Sensors Letters, vol. 6, no. 6, pp. 1–4, 2022.
  * [12] J. Chen, N. Skatchkovsky, and O. Simeone, “Neuromorphic integrated sensing and communications,” IEEE Wireless Communications Letters, vol. 12, pp. 476–480, March 2022.
  * [13] S. M. Patole, M. Torlak, D. Wang, and M. Ali, “Automotive radars: A review of signal processing techniques,” IEEE Signal Processing Magazine, vol. 34, no. 2, pp. 22–35, 2017.
  * [14] M. Liu, F. Gao, Z. Cui, S. Pollin, and Q. Liu, “Sensing with OFDM Waveform at mmWave Band based on Micro-Doppler Analysis,” arXiv preprint arXiv:2303.16559, 2023.
  * [15] W. Gerstner, W. M. Kistler, R. Naud, and L. Paninski, Neuronal dynamics: From single neurons to networks and models of cognition.  Cambridge University Press, 2014.
  * [16] D. Auge, J. Hille, and E. Mueller, “A survey of encoding techniques for signal processing in spiking neural networks,” Neural Processing Letters, vol. 53, no. 3, pp. 4693–4710, 2021.
  * [17] Y. Bengio, N. Léonard, and A. Courville, “Estimating or propagating gradients through stochastic neurons for conditional computation,” arXiv preprint arXiv:1308.3432, 2013.

  *[RF]: Radio Frequency
  *[FMCW]: Frequency-Modulated Continuous-Wave
  *[UWB]: Ultra-Wideband
  *[ISAC]: Integrated Sensing and Communication
  *[DL]: Deep Learning
  *[SNNs]: Spiking Neural Network
  *[SNN]: Spiking Neural Network
  *[AoA]: Angle of Arrival
  *[CAE]: Convolutional Autoencoder
  *[STE]: Straight-Through Estimator
  *[DFT]: Discrete Fourier Transform
  *[OFDM]: Orthogonal Frequency Division Multiplexing
  *[LIF]: Leaky Integrate and Fire
  *[TBR]: Threshold-Based Representation
  *[SF]: Step-Forward
  *[MW]: Moving-Window
  *[MSE]: Mean Squared Error
  *[LSE]: Learned Spike Encoding
  *[SNR]: signal-to-noise ratio
  *[RMSE]: Root Mean Squared Error