# The Tiered Radio Extragalactic Continuum (T-RECS) simulation II: HI emission
and continuum-HI cross-correlation.

Anna Bonaldi1, Philippa Hartley1, Tommaso Ronconi2,3,4, Gianfranco De
Zotti5,6, Matteo Bonato7
1 SKA Organization, Jodrell Bank, Lower Whitington, Macclesfield, SK11 9FT, UK
2 SISSA, Via Bonomea 265, 34136 Trieste, Italy
3 INAF - Osservatorio di Astrofisica e Scienza dello Spazio (OAS), Via Gobetti
93/3, I-40127 Bologna, Italy
4 CN-HPC: National Centre for HPC, Big Data and Quantum Computing - spoke 3,
Italy
5 INAF - Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5,
I-45122 Padova, Italy
6 Nicolaus Copernicus Academy, Poland
7 INAF - Istituto di Radioastronomia and Italian ALMA Regional Centre, Via
Gobetti 101, I-40129, Bologna, Italy E-mail<EMAIL_ADDRESS>
###### Abstract

In this paper we extend the Tiered Radio Extragalactic Continuum Simulation
(T-RECS) to include HI emission. The HI T-RECS model is based on the most
recent HI mass function estimates, combined with prescriptions to convert HI
mass to total integrated HI flux. It further models source size, morphology
and kinematics, including rotational velocity and HI line width. The continuum
T-RECS model is updated to improve the agreement with deeper number counts
available at 150 MHz. The model for star-forming galaxies (SFGs) is also
modified according to the most recent indications of a star formation rate
(SFR)–radio luminosity relation, which depends primarily on stellar mass
rather than redshift. We further introduce prescriptions to associate an HI
mass to the T-RECS radio continuum SFG and Active Galactic Nuclei (AGN)
populations. This gives us a way to meaningfully associate counterparts
between HI and continuum catalogues, thus building HI $\times$ continuum
simulated observations. Clustering properties of the sources in both HI and
continuum are reproduced by associating the galaxies to dark matter haloes of
a cosmological simulation. We deliver a set of mock catalogues, as well as the
code to produce them, which can be used for simulating observations and
predicting results from radio surveys with existing and forthcoming radio
facilities, such as the Square Kilometre Array (SKA).

###### keywords:

radio lines: galaxies, radio continuum: galaxies, galaxies: luminosity
function, mass function, large-scale structure of Universe, software:
simulations

††pagerange: The Tiered Radio Extragalactic Continuum (T-RECS) simulation II:
HI emission and continuum-HI cross-correlation.–B††pubyear: 2017

##  1 Introduction

The neutral gas component of galaxies can be traced by the 21-cm ?spin-flip?
transition of neutral hydrogen (HI), which can be detected over the radio
continuum emission at a rest frequency of $\nu_{21\rm{cm}}\simeq 1.420$ GHz
(Hellwig et al., 1970).

So far, due to the faintness of the signal, most untargeted HI observations
have observed the local Universe (e.g. HIPASS, Barnes et al. 2001, Wong et al.
2006, ALFALFA Giovanelli et al. 2005; Haynes et al. 2018), with a few targeted
fields to extend to slightly higher redshifts (AUDS Xi et al. 2021, BUDHIES
Gogate et al. 2023). A new generation of facilities is allowing HI surveys
that will greatly improve over their predecessors in terms of redshift range,
survey area, sensitivity and spatial resolution: CHILES on the Jansky Very
Large Array (JVLA, Fernández et al., 2016), WALLABY and DINGO on the
Australian SKA Pathfinder (ASKAP, Koribalski et al., 2020; Meyer, 2009),
LADUMA and MIGHTEE on the SKA Precursor MeerKAT (Blyth et al., 2016; Maddox et
al., 2021; Ponomareva et al., 2023) and the Apertif imaging surveys (Adams et
al., 2022). These observational efforts will culminate with the Square
Kilometre Array Observatory (SKAO) planned HI observations (Staveley-Smith &
Oosterloo, 2015).

Interstellar gas and stars, and the link between them, are at the heart of
galaxy formation and galaxy evolution studies. They together account for the
galaxy’s baryonic content, and the former fuels the latter through bursts of
star-formation activity. In this context, HI and continuum surveys combined
are a particularly powerful tool, with the former tracing neutral gas and the
latter tracing star formation (e.g. Maddox et al., 2016). In Bonaldi et al.
(2019) (hereafter, T-RECS I) we delivered a code,
T-RECS111https://github.com/abonaldi/TRECS.git, to produce mock radio
continuum catalogues consistent with the most recent data and models, with the
aim of predicting what upcoming and future radio continuum surveys could look
like, and to optimise survey specifications to science goals. This work
extends its scope to HI surveys, as well as continuum–HI combined studies.

One way to produce mock HI catalogues would be to add HI properties to the
radio continuum modelling presented in T-RECS I. However, this would produce
catalogues selected in radio continuum rather than HI, which would not be
adequate to simulate the untargeted HI surveys mentioned above. In light of
this, we have implemented two independent modules: one producing catalogues
selected in HI (described in Sec. 2) and one producing catalogues selected in
radio continuum (updated from T-RECS I as described in Sec. 3). As a further
step, those two outputs can be combined in a consistent HI $\times$ continuum
catalogue set as detailed in Sec. 4. This approach explicitly allows for
sources in one catalogue either to have or not have a counterpart in the
other, based on expected correlations between HI and radio continuum
properties and on the respective selection functions. It also gives the option
to produce either of the two simulations as stand-alone. This code
architecture also easily supports extending to other wavelengths in the future
(e.g. optical, infra-red). The issue of including realistic clustering into
the simulation is discussed in Sec 5. Finally we draw our conclusions in Sec
6.

In this work, $H_{0}=67$, $\Omega_{m}=0.32$, $\Omega_{\Lambda}=0.68$, which
are the Planck Collaboration VI (2020) $\Lambda$CDM best-fit parameters
rounded to two significant figures.

##  2 HI model

###  2.1 HI flux

Following Duffy et al. (2012), the total integrated HI flux of a source,
$f_{\rm HI}$, can be directly linked to the neutral hydrogen mass, $M_{\rm
HI}$. By adopting an optically-thin approximation for the HI, where self-
absorption is negligible, they obtain

$f_{\rm HI}=\frac{M_{\rm HI}}{49.8}d_{\rm L}^{-2}\,{\rm Jy\,Hz}$ (1)

where $d_{\rm L}$ is the luminosity distance of the source in Mpc and $M_{\rm
HI}$ is the HI mass in units of $M_{\odot}$. This relation essentially
converts between HI luminosity function and HI mass function and it is yet
another example of the intrinsic power of HI observations for probing galaxy
physics.

To generate T-RECS HI catalogues, we combined this relation with the best-fit
$z=0$ mass function from Jones et al. (2018):

$\Phi(M_{\rm HI})=\ln(10)\phi_{\star}(M_{\rm HI}/M_{\rm
HI}^{\star})^{(\alpha+1)}\exp[-M_{HI}/M_{\rm HI}^{\star}],$ (2)

with $\log(M_{\rm HI}^{\star})$=9.94 $M_{\odot}$, $\phi_{\star}=4.5\times
10^{-3}$, $\alpha=-1.25$.

The redshift evolution of the HIMF is currently poorly constrained by the
data, a situation that should drastically improve once the new-generation
surveys are completed.

Bera et al. (2022) and Paul et al. (2023) investigated the HIMF at $z\sim 0.3$
with GMRT and MeerKAT data, respectively. They both found a significant
evolution of both $M_{\rm HI}^{\star}$ and $\phi_{\star}$, which results in
far fewer massive HI galaxies at $z=0.3$ than at $z=0$.

In line with those results, we parameterize the evolution as

$\displaystyle\log(M_{\rm HI}^{\star}(z))$ $\displaystyle=$
$\displaystyle\log(M_{\rm HI}^{\star})_{z=0}+C_{\rm evol}\times z.$ (3)
$\displaystyle\log\phi_{\star}(z)$ $\displaystyle=$
$\displaystyle\log\phi_{\star z=0}+\phi_{\rm evol}\times z.$ (4)

We determine the parameters $C_{\rm evol}=-1.41$, $\phi_{\rm evol}=1.55$, by
imposing that, at $z=0.32$, $\log(M_{\rm HI}^{\star}/M_{\odot})=9.49$ and
$\phi_{\star}=1.4\times 10^{-2}$ Mpc-3 dex-1, which are the best-fit values of
Paul et al. (2023).

Given the lack of data on the HIMF at higher redshifts, we restrict the
redshift range of the simulation to $z=0$–0.5, which sets the frequency range
for HI surveys we can simulate to 1420.40–946.9 MHz. Fig.1 shows the adopted
HI mass function (black line) compared with a T-RECS realization (black
symbols) for different redshifts.

Figure 1: HI mass function from Jones et al. (2018) with redshift evolution
constrained by Paul et al. (2023, back lines), compared to the distribution of
T-RECS HI masses for an HI catalogue of 25 deg2, and flux limit $f_{\rm HI}=1$
Jy Hz (black symbols). In red we also show the distribution of the HI mass
proxy obtained for a companion continuum simulation of 25 deg2, $f_{1.4\rm
GHz}\geq 100$ nJy (see sec. 3.4).

###  2.2 Source size and morphology

There is a tight relation between HI mass and the diameter of the HI disk
$D_{\rm HI}$ in kpc, defined at an HI surface density $\Sigma_{\rm HI}$ of 1
$M_{\rm sun}\,$pc-2 (e.g., Broeils & Rhee, 1997; Verheijen & Sancisi, 2001;
Swaters et al., 2002; Noordermeer et al., 2005; Wang et al., 2014; Ponomareva
et al., 2016; Naluminsa et al., 2021), which appears to be independent of
redshift in the so far explored range (Rajohnson et al., 2022). We adopt the
recent relation by Naluminsa et al. (2021):

$\log M_{\rm HI}=(1.95\pm 0.03)\log D_{\rm HI}+(6.5\pm 0.04),$ (5)

which was derived by using a complete sample of 228 WHISP galaxies including
all morphological types. This relation is in particularly good agreement with
Broeils & Rhee (1997) and Wang et al. (2016), which have a similar
morphological selection. It is steeper than Verheijen & Sancisi (2001);
Swaters et al. (2002); Noordermeer et al. (2005), which focused on specific
morphological types.

$D_{\rm HI}$, obtained by inverting eq. (5), is finally converted from the
physical size in kpc to an angular size in arcsec depending on the angular
diameter distance of the source. For consistency with the continuum module,
where the size of SFGs is in terms of a scale radius, the T-RECS catalogue
contains the radius $R_{\rm HI}=D_{\rm HI}/2$. To further complete the
morphological description of the HI sources, the catalogue contains galaxy
inclination $i$ which, in the hypothesis of a random orientation, follows a
$\sin(i)$ distribution.

Following e.g. Holmberg (1958), galaxy inclination is linked to the axis ratio
$b/a$

${\rm cos}^{2}i=\frac{(b/a)^{2}-\kappa^{2}}{1-\kappa^{2}}$ (6)

where $\kappa$ is the ratio of smallest to largest axis of an oblate spheroid
which best represents the galaxy’s 3-dimensional shape. Early measurements
indicated values of $k=0.2$ for spirals and $k=0.5$ for ellipticals. Rodríguez
& Padilla (2013) derive updated distributions for $\kappa$ by studying the
Sloan Digital Sky Survey (SDSS) Galaxy Zoo objects. A Gaussian fit to these
distributions yields $\kappa=0.267\pm 0.102$ for spirals and $\kappa=0.438\pm
0.196$ for ellipticals. We use the previous relations and the galaxy’s size
and inclination to the derive major and minor axes. Since HI-rich galaxies are
predominantly of spiral morphology, we adopt for the whole HI sample the
spiral distribution of Rodríguez & Padilla (2013).

###  2.3 Source HI line width

The width of the 21-cm emission line of neutral hydrogen is related to the
circular velocity, corrected for the galaxy’s inclination, by

$w_{50}=2\sin(i)v_{\rm max}$ (7)

(e.g. Oman, 2022), where $w_{50}$ is the full with at half maximum (FWHM) of
the 21 cm line and $v_{\rm max}$ is the maximum circular velocity.

Katz et al. (2019) derive empirical scaling relations of $v_{\rm max}$ with
the galaxy’s baryonic mass $M_{\rm b}$ (Baryonic Tully-Fisher relation, BTFR),
as well as dark matter halo mass, $M_{\rm h}$, from 120 late-type galaxies
from the SPARCS database. We use their BTFR results obtained by using a ?flat?
prior on the mass-to-light ratio, to derive $v_{\rm max}$ from the baryonic
mass. The latter has been modelled as $M_{\rm b}=M_{\rm HI}+M_{\rm
star}+M_{\rm H2}+M_{\rm HII}$, where $M_{\rm star}$ is the stellar mass,
$M_{\rm H2}$ is the mass of molecular hydrogen and $M_{\rm HII}$ is the mass
of ionised hydrogen. Although there are other baryonic components, the list
above contains the majority of baryons in a galaxy.

For $M_{\rm star}$ we use the maximum likelihood $M_{\rm star}$–$M_{\rm HI}$
relation by Naluminsa et al. (2021). Walter et al. (2020) investigate the
density of atomic and molecular hydrogen, $\rho_{\rm HI}(z)$ and $\rho_{\rm
H2}(z)$, in galaxies as a function of redshift. For $M_{\rm H2}$ we use
$M_{\rm H2}=M_{\rm HI}\,\rho_{\rm H2}(z)/\rho_{\rm HI}(z)$ where $\rho_{\rm
H2}(z)$ and $\rho_{\rm HI}(z)$ are the best-fit relations of Walter et al.
(2020) (their table 1). For $M_{\rm HII}$ we use the $\log(M_{\rm
star})$–$\log(M_{\rm HII}/M_{\rm star})$ relation by Popping et al. (2014)
(linear fit to their Figure 2). This relation is provided in the $\log M_{\rm
star}=6-13\,\log M_{\odot}$ and implies an increase of the ionised hydrogen
with decreasing stellar mass. To avoid $M_{\rm HII}$ to become unrealistically
high for low values of $M_{\rm star}$, we cap the maximum value of $M_{\rm
HII}$ to $M_{\rm HI}$. Once $v_{\rm max}$ has been obtained, by using the
$M_{\rm h}$–$v_{\rm max}$ relation by Katz et al. (2019) we add the modelling
of the dark matter mass from HI properties, $M_{\rm h,HI}$, which is used to
model source clustering in Sec 5.

In Fig. 2 we compare the distribution of the T-RECS $w_{50}$ at low redshifts
(red symbols) with the HI velocity width function (HIWF) derived from the
ALFALFA survey by Oman (2022; black lines). Despite the relative simplicity of
the model, there is a good agreement between the two. The decreasing number of
T-RECS objects below $w_{50}\sim 10^{2}$ is due to the catalogues becoming
increasingly incomplete. When including higher redshifts, the T-RECS HIWF
slowly shifts to lower $w_{50}$ values and increases in amplitude due to the
evolution of the HI mass function. We note that our modelling of the line
width relies on the presence of a plateau in the velocity curve for $v=v_{\rm
flat}$. This is not representative of irregular or dwarf galaxies (Swaters et
al., 2009; Oh et al., 2015) whose measured velocity width would be smaller.

Figure 2: Black lines: HIWF from Oman (2022) computed from the $\alpha$.100
catalogue (solid, dashed and dot-dashed lines are for the ?complete?, ?spring?
and ?fall? sky areas, respectively). Red symbols: T-RECS HIWF computed on an
HI catalogue of 25 deg2, and flux limit $f_{\rm HI}=1$ Jy Hz, over the 0–0.075
redshift range.

##  3 Radio continuum model

The radio continuum model is described in detail in T-RECS I; here we give a
summary and introduce a few modifications. T-RECS models radio continuum
sources as either star-forming galaxies (SFGs) or Radio Loud Active Galactic
Nuclei (RL AGN). There is no explicit modelling of radio-quiet (RQ) AGN (e.g.,
Kellermann et al., 2016; Padovani et al., 2015; Mancuso et al., 2017; White et
al., 2017; Hartley et al., 2019). RQ AGN can be modelled as objects where star
formation in the host galaxy is responsible for the radio emission. On the
other hand, observations have also shown that the active nucleus can produce
the radio emission. RQ AGN would contribute part of the flux of those sources
that, in T-RECS, are modelled as SFGs.

###  3.1 Star-forming galaxies (SFGs)

SFGs are modelled as in Mancuso et al. (2015) as late-type, spheroidals and
lensed spheroidals. The radio emission is based on redshift-dependent star-
formation rate (SFR) functions and a modelling of synchrotron, free-free and
thermal dust emission as a function of SFR for each of the three sub-
populations. The free-free emission and the thermal dust emission are still
modelled as in T-RECS I. The synchrotron emission, which is overall the
dominant component, is updated in this work to follow the most recent results.

Across the 150 MHz–20 GHz frequency range, we adopt

$L_{\rm synch}(\nu)=L({\rm
SFR})\,\left(\frac{\nu}{\nu_{0}}\right)^{\alpha+\delta\alpha\log(\nu/\nu_{0})},$
(8)

where $L({\rm SFR})$ follows Smith et al. (2021)

$\frac{L(\rm SFR)}{{\rm W/Hz}}=L_{0}\left(\frac{\rm SFR}{M_{\odot}/{\rm
yr}}\right)^{\beta}\,\left(\frac{M_{\rm
star}}{10^{10}\,M_{\odot}}\right)^{\gamma}.$ (9)

The frequency dependence is a power-law with spectral index that progressively
steepens towards higher frequencies; the parameters are $\alpha=-0.85$,
$\delta\alpha=-0.1$, $\nu_{0}=1.6$ GHz. The previous relation we adopted, as
in Mancuso et al. (2015), already had a steepening towards the highest
frequencies. The most recent data at 150 MHz (Siewert et al., 2019; Mandal et
al., 2021) indicate that this behaviour should be extended to the whole
frequency range. A possible mechanism for a steepening synchrotron spectral
index is synchrotron ageing.

$M_{\rm star}$ in eq. (9) is the stellar mass of the galaxy. This is modelled
starting from the SFR using the best-fit SFR–$M_{\rm star}$ relation of Aversa
et al. (2015) (their table 2, with redshift evolution). The redshift evolution
of $M_{\rm star}$ induces a redshift evolution on the radio luminosity
function. Parameters in eq. (9) are $\beta=0.850\pm 0.005$, $\gamma=0.402\pm
0.005$ (Smith et al., 2021) and $\log(L_{0}/[{\rm WHz}^{-1}])=21.28$. The
latter differs from the Smith et al. (2021) value of $\log
L_{0}/[\rm{W/Hz}]=22.10$, however our normalization is linked to that of eq.
(8) at 150 MHz $\log(L_{\rm sync}(150{\rm MHz})/\log({\rm SFR}))=0.77$, which
needs a $d\log L_{0}$=-0.77 correction. Once this is accounted for, there is
still a difference in the overall normalization of $d\log L_{0}$=-0.06,W/Hz.
Given that $L({\rm SFR})$ depends on both SFR and $M_{\rm star}$ as per eq.
(9), differences in the distribution of either quantities between the
measurements by Smith et al. (2021) and the models in Mancuso et al. (2015)
and Aversa et al. (2015) could cause some discrepancy. Furthermore, the
parameters we adopt have been chosen to achieve consistency with the available
data over the whole 150 MHz–20 GHz range rather than being optimised at 150
MHz.

In Fig. 3 we verify that the updated luminosity functions at 1.4 GHz are in
good agreement with the data. The updated comparison of differential source
counts in the 150 MH–20 GHz range is presented in Fig. 4.

As in T-RECS I, the size and shape of SFGs is modelled in terms of elliptical
Sersic profiles. As a refinement, we link galaxy inclination to galaxy
ellipticity using again eq. (6) with $\kappa=0.267\pm 0.102$ for spirals and
$\kappa=0.438\pm 0.196$ for ellipticals. We make a direct association with the
T-RECS radio continuum sub-populations as follows: late-type to spirals,
spheroids and lensed spheroids to ellipticals (Rodríguez & Padilla, 2013).

Figure 3: Comparison of continuum luminosity functions between T-RECS and the
available data from Novak et al. (2017); Butler et al. (2019); Ocran et al.
(2020); Bonato et al. (2021)

.

###  3.2 Radio Loud (RL) AGN

RL AGN are based on the Massardi et al. (2010) evolutionary model as updated
by Bonato et al. (2017). This model classifies the AGN into steep-spectrum
sources (SS-AGN), flat-spectrum radio quasars (FSRQs) and BL Lacs, with
different evolutionary properties.

The basis of our AGN model is given by luminosity functions at 1.4 GHz for
these populations and for redshift intervals ranging from $z=0$ to $z=8$.
These trivially translate to flux density at 1.4 GHz given the redshift of the
source and the adopted cosmology.

To compute flux densities at other frequencies, Bonato et al. (2017) adopt
power-law spectra $S\propto\nu^{\alpha}$ with $\alpha_{\rm FSRQ}=\alpha_{\rm
BLLac}=-0.1$, and $\alpha_{\rm steep}=-0.8$. T-RECS I modified this simple
recipe by allowing for a scatter in the spectral index between different
sources. Moreover, it introduced systematic variations with frequency of the
spectral index distributions, in order to achieve agreement with multi-
frequency observations.

An effective spectral index between the frequencies $\nu_{1}$ and $\nu_{2}$ of
sources of a given population with flux density $S_{1}$, within $dS_{1}$, at
$\nu_{1}$, was computed by finding the flux density $S_{2}$ at $\nu_{2}$ such
that $N_{1}(S_{1})dS_{1}=N_{2}(S_{2})dS_{2}$. Thus $\alpha_{\rm
eff}(\nu_{1},\nu_{2})$ is the single spectral index relating the counts at
$\nu_{1}$ to those at $\nu_{2}$. This procedure was applied to the 1.4–4.8 GHz
and 4.8–20 GHz range, using the Massardi et al. (2010) model up to 5 GHz and
the De Zotti et al. (2005) model at higher frequencies. In this work, we
extend this approach to the 150 MHz–1.4 GHz frequency range, where we had
previously applied the Bonato et al. (2017) indices plus scatter, in order to
improve the agreement with the Siewert et al. (2019) and Mandal et al. (2021)
counts at 150 MHz.

The morphological model for RL AGN is the same as in T-RECS I. Following the
unified model of AGN, we assume the same parent population for all our AGN
classes. Intrinsic sizes are drawn using DiPompeo et al. (2013); BL-Lac and
FSRQ are then assigned a small viewing angle ($\theta\leq 5$ deg) and SS-AGN a
large viweing angle ($5<\theta\leq 90$ deg), thus creating the familiar
dichotomy in observed source sizes. For SS-AGN, the description is completed
by Fanaroff & Riley (1974)’s $R_{\rm s}$ parameter, modelled following Lin et
al. (2010), to provide the FRI/FRII classification.

For consistency with the SFG catalogue, we complement the morphological
description of AGN with an elliptical/spiral classification of their host
galaxy. The issue of radio vs optical morphology was investigated by Kozieł-
Wierzbowska et al. (2020) on a sample of 32,616 spectroscopically-selected
galaxies from SDSS with radio counterpart from FIRST and VLASS. For our AGN,
we adopt their result on galaxies with extended radio morphology, which have
an elliptical host in 98 % of cases.

###  3.3 Dark matter mass modelling from radio continuum properties

As in T-RECS I, we include an estimate of the dark matter (DM) mass of
galaxies based on radio continuum properties $M_{\rm h,cont}$. The main use of
this mass is to simulate the clustering properties for our sources. Special
attention has been given to ensuring that the obtained T-RECS DM mass
distribution matches well the mass function from the P-millennium (Baugh et
al., 2019) cosmological simulation, which is used to simulate clustering (see
Sec. 5 for more details).

The DM mass modelling is the same as T-RECS I. For SFGs, we used a $L_{\rm
SFR}$–$M_{\rm h}$ relation of the form:

$L(M_{\rm h})=N\times\Big{[}\Big{(}\frac{M_{\rm
h}}{M_{b}}\Big{)}^{\alpha}+\Big{(}\frac{M_{\rm
h}}{M_{b}}\Big{)}^{\omega}\Big{]}^{-1},$ (10)

where $N$, $\alpha$, $\omega$ and $M_{b}$ are free parameters which include
redshift evolution. We start from the Baugh et al. (2019) target halo mass
function as a function of redshift, and we derive the parameters so that the
resulting luminosity function matches as closely as possible the radio
luminosity function at 1.4 GHz from Bonato et al. (2017). Once the parameters
are derived, eq. (10) can be inverted to obtain $M_{\rm h}$ from the 1.4 GHz
luminosity.

For RL AGN, we start again from the halo mass function of the Baugh et al.
(2019) simulation to draw a sample of galaxies with plausible DM mass
distribution for each T-RECS redshift interval. For all galaxies in the
sample, we compute the stellar mass with the Aversa et al. (2015) $M_{\rm
star}=F(M_{\rm h})$ relation (their Table 2, including redshift evolution).

Janssen et al. (2012) derived the fraction of galaxies hosting an RL AGN as a
function of the host galaxy stellar mass, $M_{\rm star}$, separately for Low
Excitation Radio Galaxies (LERGs) and High Excitation Radio Galaxies (HERGs).
We use this result to extract subsamples of HERGs and LERGs from the initial
galaxy sample. We finally compute the $M_{\rm h}$ distribution of the two sub-
samples, which, once appropriately normalised, can be used to draw halo masses
suitable for a HERG/LERG host. We map the HERG/LERG populations into our
T-RECS observational categories as follows: FSRQs to the HERG population; BL
Lacs to the LERG population; SS-AGN morphologically classified as FR II/FR I
to the HERG/LERG population (see TRECS I for more details).

Figure 4: Comparison of differential source counts in total intensity at 150
MHz (top left) 1.4 GHz (top right), 3 GHz (bottom left) and 15 GHz (bottom
right) between T-RECS and the available data from Franzen et al. (2016); Bondi
et al. (2008); Vernstrom et al. (2016); Smolčić et al. (2017); Padovani et al.
(2015); Ibar et al. (2009); Whittam et al. (2016); Waldram et al. (2010); AMI
Consortium et al. (2011); Mandal et al. (2021); Bonato et al. (2021)

.

###  3.4 HI mass modelling from radio continuum properties

The interplay between gas and stars in the lifecycle of a galaxy results in a
correlation between HI mass, stellar mass, and SFR, which is exemplified by
the Kennicutt-Schmidt Law (Kennicutt, 1998). Scaling relations between HI mass
tracers and SFR tracers have been derived (e.g. Zhou et al., 2018; Bothwell et
al., 2013; Catinella et al., 2010; Cortese et al., 2011; Jaskot et al., 2015;
Huang et al., 2012; Sinigaglia et al., 2022, just to cite some of the most
recent); the correlation is well established, with a scatter indicative of the
many factors that regulate SFR. Among those, there is feedback from AGN
activity, that heats the surrounding gas therefore reducing star formation.
Additionally, an important role is played by galaxy type. Late-type galaxies,
associated with spiral/disk morphologies, are typically HI-rich and actively
forming stars, while early-type galaxies, with spheroidal/elliptical
morphology, are typically HI-poor and have low star formation.

The SFR–$M_{\rm{HI}}$ correlation mentioned above provides a relatively
straightforward way to add HI mass to the modelling of T-RECS star-forming
galaxies. We call this HI mass $\tilde{M_{\rm HI}}$, to distinguish it from
the HI mass modelled for HI sources in Sec. 2. It can be written as:

$\log(\tilde{M}_{\rm HI})=A\log({\rm SFR})+B$ (11)

with a scatter $\sigma$; the slope and normalization $A$ and $B$ are in
general redshift-dependent.

We derive the parameters $A$ and $B$ by imposing that the resulting
$\tilde{M}_{\rm HI}$ distribution is consistent with the HI mass function
adopted for the HI model in Sec. 2. The redshift evolution as in Paul et al.
(2023) means that the maximum HI mass changes relatively rapidly as a function
of redshift. The evolution of the SFR function in Mancuso et al. (2015) is not
as fast, and imparts a significant redshift dependence of the parameters $A$
and $B$. We find that

$\displaystyle A(z)$ $\displaystyle=$ $\displaystyle 1.1-2.46\,z+2.06\,z^{2}$
(12) $\displaystyle B(z)$ $\displaystyle=$ $\displaystyle 9.5-z$ (13)

and a scatter $\sigma\log\tilde{M}_{\rm HI}=0.2$ give a good agreement in the
considered redshift range.

In Fig.5 we show the $\log(\rm SFR)$–$\log(\tilde{M}_{\rm HI})$ relation for
the T-RECS catalogue up to $z=0.3$ compared to the best fits from Naluminsa et
al. (2021), Cluver et al. (2014) and Michałowski et al. (2015). The slope of
the correlation is consistent with the observational estimates. However, our
normalization is higher, at $\Delta\log({\rm SFR})\sim 0.3$–0.4. Using one of
the observed relations instead of the abundance-matching parameters produces
larger HI masses than what would be consistent with the HI mass function by
Jones et al. (2018) plus redshift evolution. In other words, this is
essentially a tension between the Mancuso et al. (2015) model of SF and the
Jones et al. (2018) HI mass function, which indicates that there is still some
distance to cover to get a consistent model across both observables.

Several relations have been published between $M_{\rm HI}$ and the stellar
mass, $M_{\rm star}$ (e.g. Catinella et al., 2010; Huang et al., 2012; Maddox
et al., 2015; Naluminsa et al., 2021), which could be used as an HI mass proxy
for the AGN population. The relation appears to be linear at the low-mass end,
and to flatten at the high-mass end due to the dominance of early-type, HI-
poor galaxies. We adopt the maximum likelihood relation by Naluminsa et al.
(2021) up to $M_{\rm star}=10^{10}\,M_{\odot}$ and Catinella et al. (2010) for
higher masses. The latter was derived on an unbiased sample of galaxies
selected by stellar mass and therefore should contain the gas-poor systems
that are typically missed by the current untargeted HI surveys.

In Fig. 1 we compare the HI mass function (black line and symbols) with the HI
continuum mass proxy distribution (red symbols). Overall, there is a
satisfactory agreement between the two, ensured by the abundance-matching
method used to derive the parameters in eq. (11), as well as a plausible
modelling of the HI mass proxy for the AGN population.

The discrepancy at the low-mass end is due to the different selection
functions. As mentioned, the selection in HI flux almost corresponds to a
selection in mass, and gives high completeness of the sample down to the last
mass bin. Conversely, the selection in continuum flux means that the mass
proxy sample is progressively incomplete. Departures from the two
distributions not due to completeness can be ascribed to the simplicity of eq.
(11), which cannot guarantee a perfect fit over the whole mass range. At the
high-mass end, the mass-proxy distribution tends to over-predict the number of
objects. This region of the distribution contains very few objects, and
therefore its weight in the abundance-matching scheme is relatively low.

Figure 5: Relation between SFR and $\tilde{M}_{\rm HI}$ for the $z\leq 0.3$
continuum catalogue, compared to Naluminsa et al. (2021) (dashed line), Cluver
et al. (2014) (dot-dashed line) and Michałowski et al. (2015) (solid line).
Contours include are 99%, 90% and 50% of all galaxies.

##  4 Continuum $\times$ HI model

Figure 6: _Left:_ Distribution of $M_{\rm HI}$ for an HI catalogue of 25 deg2
$f_{\rm HI}\geq 1$ Jy/Hz and $z=$0.475–0.525 (solid line) and of
$\tilde{M}_{\rm HI}$ for a continuum catalogue with the same sky area and
redshift range and $S_{\rm 1.4\,GHz}\geq 100\,$nJy (dot-dashed line). _Right:_
distribution of $M_{\rm HI}$ for the cross-matched catalogue (red solid line);
HI unmatched sources (black solid line) and continuum un-matched sources
(black dot-dashed line).  Figure 7: Selection function as a function of HI
mass for the HI catalogue (black solid line) and the matched continuum
counterparts (black diamonds and dot-dashed line) for the example in Fig. 6.
Superimposed in red is the theoretical selection function corresponding to a
Gaussian error $\sigma\log M=0.65$ and an offset $\Delta\log M=0.3$, which is
a good fit to the observed selection.

So far, the T-RECS HI and continuum produced with the modules described in
Secs. 2 and 3 are independent. Each produces statistically plausible samples
of sources, but the two samples are uncorrelated and have different sky
coordinates.

If the two catalogues were used to represent the same area of the sky, one
would expect a certain number of positional cross-matches between the two, as
the same galaxy is visible in both bands. In this section, we discuss how a
consistent HI $\times$ continuum catalogue can be constructed by associating
counterparts between the two catalogues and rewriting one of the two
coordinate sets to match.

First of all, counterparts are only explored between galaxies that belong to
the same redshift slice between the two simulations. We consider a sample
$\\{i\\}$ of continuuum sources and $\\{j\\}$ of HI sources, which have been
selected with $S^{i}>S1$ and $f^{j}_{\rm HI}>F1$, where $S^{i}$ are the fluxes
in continuum, $f^{j}_{\rm HI}$ the fluxes in HI, and $S1$, $F1$ are the
respective detection limits. If we ignore redshift variation within a slice,
from eq. (1) we can rewrite the HI selection in terms of HI mass, as
$M^{j}_{\rm HI}>M1$ where $M^{j}_{\rm HI}$ and $M1$ are the HI masses
corresponding to $f^{j}_{\rm HI}$ and $F1$ at the given redshift.

The counts of both the $\\{i\\}$ and $\\{j\\}$ samples come from the
appropriately normalised distributions (luminosity function, SFR function and
HI mass function) and need to be conserved. Therefore, all objects are
retained, and in general only some of them have a counterpart. Based on the HI
selection function, continuum sources belonging to $\\{i\\}$ would be visible
also in HI if their HI mass is above M1. Assessing this requires a modelling
of the HI mass for the continuum sources, $\tilde{M}^{i}_{\rm HI}$, which we
have introduced and described in Sec. 3.4.

We note that $\tilde{M}^{i}_{\rm HI}$ does not come from sampling an expected
HI mass distribution and so selecting a subsample of potential counterparts as
$\tilde{M}^{i}_{\rm HI}>M1$ would not automatically reproduce HI source counts
and could lead to objects that are mismatched with $\\{j\\}$. However, if we
have compatible $\tilde{M}_{\rm HI}$ and $M_{\rm HI}$ distributions, the
differential counts are similar $N(M_{\rm HI})\sim N(\tilde{M_{\rm HI}})$ and
so are the number of sources above the same mass limit. In this case,
assigning counterparts based on matching $M_{\rm HI}$ from the HI catalogue
with $\tilde{M}_{\rm HI}$ from the continuum catalogue produces a similar
selection function.

As already discussed in Sec. 3.4, our modelling of $\tilde{M}_{\rm HI}$ for
the SFGs, which are the vast majority of the continuum sources, has been
explicitly designed to reproduce the HI mass function as best as possible, by
using an abundance-matching technique to constrain its parameters (see also
Fig 1). The mass match ensures that broadly compatible objects between the two
catalogues are assigned, and it propagates, although with additional scatter,
the modelled correlations between continuum and HI properties. For example,
the SFR–$M_{\rm HI}$ correlation that has been modelled for the continuum
catalogue shown in Fig. 5 is conserved when we consider the HI $\times$
continuum catalogue.

The HI selection $M^{j}_{\rm HI}>M1$ maps to a selection in
$\tilde{M}^{j}_{\rm HI}$ which is smoother than the original step function,
due to the various uncertainties introduced by the additional modelling.

The HI mass-matching itself is performed with a nearest neighbour
implementation. Initially, for each galaxy in the HI catalogue, the 20 nearest
neighbours from the continuum catalogue are retained, as counterpart
candidates. The HI galaxies are then reprocessed in order of decreasing HI
mass, and a unique counterpart is assigned to each of them, if available, as
the best of the 20 candidates that has not been already matched. The 20
nearest neighbours have been chosen as a trade-off between computational speed
and having enough samples for a good mass-matching. The reason for proceeding
in inverse mass order is that high-mass objects are more rare and therefore
more difficult to match accurately, but at the same time brighter and
therefore more likely to have a counterpart.

Figure 8: Median T-RECS $M_{\rm h,cont}$–$M_{\rm HI}$ (T-RECS 1) and $M_{\rm
h,HI}$–$M_{\rm HI}$ (T-RECS 2) compared to HI–halo mass relation from
observations (Guo et al., 2020; Rhee et al., 2023, central galaxies only) and
theoretical modelling (Baugh et al., 2019; Chauhan et al., 2021). The contours
are the full distribution for the T-RECS 1 model, including 99%, 90% and 10%
of the galaxies.

In the left panel of Fig. 6 we show the histograms of the HI masses $M_{\rm
HI}$ for an HI catalogue with $f_{\rm HI}>1$ Jy Hz and $\tilde{M}_{\rm HI}$
for a continuum catalogue with $S_{1.4\,\rm{GHz}}>$100 nJy, over 25 deg2 and
$z=0.5$. The HI catalogue has a mass limit $\log M1=6.45$. Similar numbers of
galaxies are available in both catalogues as a function of $M_{\rm HI}$ for
$\log M_{\rm HI}\geq 7$. In the $\log M_{\rm HI}=6$ mass bin, the counts
$N(M_{\rm HI})$ and $N(\tilde{M_{\rm HI}})$ start to diverge, and there is an
excess of sources in the continuum catalogue.

In the right panel of Fig. 6 we show the mass distribution of the cross-
matched catalogue with the red line, and the un-matched sources with the black
solid line and the black dash-dotted line for HI and continuum, respectively.
In this case, the continuum is deeper than the HI and most HI sources have a
continuum counterpart. When considering another couple of catalogues with
different detection limits, the proportion of sources in the matched, HI-only
and continuum-only categories would change. The $M_{\rm HI}$–$\tilde{M}_{\rm
HI}$ relation for the matched objects is very tight, with $\log{M_{\rm
HI}}-\log\tilde{M}_{\rm HI}=2.4\pm 9.4\times 10^{-3}$.

In Fig. 7 we compare the selection in mass of the HI catalogue ($\log
M1>6.45$) with that of the matched continuum counterparts (diamonds and dot-
dashed line). The error bars are consistent with Poisson statistcs and reflect
the very different number of sources that are present in the different mass
bins.

In the Gaussian approximation, for a scatter $\sigma_{M}$ on the mass $M$ and
a detection limit $M1$, the selection function $F(M)$ is an error function

$F(M)=\frac{1}{2}\left[1-{\rm
erf}\left(\frac{M_{1}-M}{\sigma_{M}\sqrt{2}}\right)\right].$ (14)

The red line in Fig. 7 is an error function similar to the measured selection,
which corresponds to a mass limit $\log\tilde{M}_{1}=6.75\pm 0.65$. The
scatter is due to that of the SFR–$\tilde{M}_{\rm HI}$ relation and that of
the mass matching. The offset $\Delta\log M=0.3$ is due to differences in the
distributions of $\tilde{M}_{\rm HI}$ and $M_{\rm HI}$, in particular the
excess of continuum sources in the $\log M_{\rm HI}=6$ mass bin. At the
highest masses, the observed selection departs from the error function, with a
few radio continuum galaxies not having an HI counterpart. This is due to
discrepancies between the high-mass end of the HIMF and the HI mass proxy
distribution (see also Fig 1).

The difficulty in accurately reproducing the expected HI mass function
starting from a radio continuum catalogue is one of the reasons behind our
choice to model the HI independently and then cross-match with the continuum
catalogue. In this way those discrepancies, rather than affecting the number
of HI sources for an HI selection, only affect the number of continuum
counterparts.

We note that the method we described does not necessarily identify pairs of
galaxies that are the most consistent in terms of properties other than the HI
mass. To achieve this, a better approach would be to simultaneously match
multiple catalogue attributes at once, for example $M_{\rm h}$, $M_{\rm star}$
and $M_{\rm HI}$. However, this approach would increase the scatter between
the matched $\tilde{M_{\rm HI}}$ and $M_{\rm HI}$ which would result in a
selection of the continuum sources largely independent from the HI signal
strength.

We point out that the selection functions adopted in T-RECS are an
approximation of what actually achieved by real observations. In radio
continuum, the ability to detect a source of a given integrated flux density
depends, amongst other things, on how much of that flux is contained within
the instrumental point spread function. Similarly, completeness for a given HI
flux depends on the 3D size and shape of the source, particularly the line
width, with sources having narrower line widths being easier to detect (see,
e.g., Haynes et al., 2011). Instrumental specifications, such as spatial and
spectral resolution, need to feature in an accurate description of selection
functions. Additionally, an important role is played by the adopted source
finding algorithms, which perform spatial and spectral filtering of the
sources to maximise detection (see, e.g. Bonaldi et al., 2021; Hartley et al.,
2023, for a comparison of different source finding methods on the same data).
If the goal is to accurately reproduce the selection functions of specific
surveys, T-RECS catalogues should be generated with flux limits consistent
with the deepest achievable fluxes, and further filtering of the source lists
should be performed on additional source properties, as appropriate.

##  5 Clustering

As described in T-RECS I, the clustering properties of the sources are
simulated by associating each galaxy with a dark matter sub-halo of the
P-millennium (Baugh et al., 2019) cosmological simulation. The cosmological
parameters adopted in this simulation are: $H_{0}=67.77$ km s-1 Mpc-1,
$\Omega_{\Lambda}=0.693$, $\Omega_{\rm M}$=0.307, $\sigma_{8}=0.8288$ (Planck
Collaboration XVI, 2014).

The reason for choosing this simulation is its very high mass resolution,
which means that sub-haloes are tracked individually down to $M_{h}=1.6\times
10^{9}\,M_{\odot}$, but also a large box size (800 comoving Mpc/h). The latter
allows producing a $5\times 5$ deg2 lightcone up to $z=8$, which is therefore
the maximum T-RECS survey size for which we can simulate clustering. This
approach essentially matches galaxies to haloes by means of the dark mass
$M_{\rm h}$, separately for each redshift slice, and overwrites the T-RECS
latitude, longitude and redshift coordinates with those of the matched haloes.
The clustering is therefore simulated in 3 dimensions.

Two predictions for the dark matter mass of T-RECS galaxies are available:
$M_{\rm h,HI}$ as modelled from the HI properties in Sec 2, and $M_{\rm
h,cont}$ as modelled from the radio continuum properties for SFGs and AGN in
Sec 3. Given the uncertainties related with both modelling and measurement of
dark matter mass, both quantities are retained as alternative models in an HI
$\times$ continuum catalogue.

Several works have explored the total HI mass in dark matter haloes (HI-halo
mass relation, HIHM), which is important to understand the role of HI in
galaxy formation and its connection to structure formation. This relation has
been investigated extensively in different theoretical models (Barnes &
Haehnelt, 2014; Paul et al., 2018; Obuljen et al., 2019; Kim et al., 2017;
Villaescusa-Navarro et al., 2018) and observationally, via spectral stacking
(Guo et al., 2020; Rhee et al., 2023). Theoretical studies predict a positive
correlation between $M_{\rm h}$ and $M_{\rm HI}$ in both the low-$M_{\rm h}$
regime, which is dominated by central galaxies, and the high-$M_{\rm h}$
regime, dominated by groups and satellite galaxies. The overall relation has a
flattening at masses in between these two regimes. Observational estimates
have so far measured a much flatter relation, with $M_{\rm HI}$ only mildly
depending on $M_{\rm h}$. Possible observational biases to account for these
discrepancies are discussed in details in Chauhan et al. (2021).

In Fig. 8 we compare the $M_{\rm h,HI}$–$M_{\rm HI}$ (T-RECS B) and $M_{\rm
h,cont}$–$M_{\rm HI}$ (T-RECS A) relations with some of the HIMH in the
literature. They can only be qualitatively compared, since T-RECS’s are
galaxy-galaxy relations and the MIMH is the total HI content of haloes
including multiple galaxies. The comparison is more meaningful at the lowest
halo masses, where a single galaxy occupancy of haloes is prevalent, or when
considering the central galaxies’ contribution to the total HIMH. The filled
circles show the median T-RECS A and T-RECS B relations; black contours are
the full T-RECS A $M_{\rm h,cont}$ model including 99%, 90% and 50% of the
galaxies. Diamonds and squares are the central galaxies’ observational
relations by Rhee et al. (2023) and Guo et al. (2020). Dashed and dash-dotted
lines are the theoretical relations by Baugh et al. (2019) and Chauhan et al.
(2021). Error bars on the observational estimates have been omitted for figure
clarity.

The T-RECS A model is in qualitatively good agreement with both the
observational estimates and with the Chauhan et al. (2021) model over the mass
range $10^{11}<M_{\rm h}<10^{13}$ $M_{\odot}$. It also manages to reproduce
the flattening of the relation towards high masses, which in T-RECS is the
result of associating some HI galaxies to AGN sources, which typically have a
larger halo mass. The T-RECS B model does not present any flattening and it
predicts a higher HI mass for the same dark mass. This is because the relation
used to model $M_{\rm h,HI}$ has been derived for late-type galaxies, which
have higher HI fraction. As such, this relation fails to capture the most
massive, early-type galaxies which are HI-poor. A lower DM mass for these
sources would result in lower clustering. For this reason, in the clustering
analysis that follows, the dark matter mass of the continuum counterpart
$M_{\rm h,cont}$ has been used, whenever present.

The implementation of the DM mass matching is based on the same nearest-
neighbour approach used for the $M_{\rm HI}$ mass matching in Sec. 4. Given
the very large number of DM sub-haloes contained in the lightcone, however,
comparing every galaxy to every halo is a significant bottleneck in terms of
execution time. At low $M_{\rm h}$ especially, there is an over-abundance of
DM haloes with respect to galaxies, due to the shape of the DM mass function
and selection effects for the T-RECS catalogues. On those mass ranges, and
depending on the sizes of the samples to be compared, a subsample of the
available sub-haloes can be randomly selected for the subsequent analysis.
This has very significant advantages in terms of computational cost, at the
price of a negligible increase of the scatter between the matched masses.

###  5.1 Validation for continuum sources

Several observational estimates of clustering for AGN and SFGs are present in
the literature (Magliocchetti et al., 2017; Hale et al., 2018; Chakraborty et
al., 2020; Bonato et al., 2017). Details of the source selection for those
studies are collected in Table 1.

Table 1: Comparison on the selection between different radio continuum clustering estimates in the literature Observation |  $\nu$[GHz] |  $S_{{\rm limit},\nu}$[mJy] | Ref.
---|---|---|---
COSMOS | 3 | 0.013 | Hale et al. (2018)
COSMOS | 1.4 | 0.15 | Magliocchetti et al. (2017)
LH | 1.4 | 0.15 | Bonato et al. (2021)
EN1 | 0.612 | 0.05 | Chakraborty et al. (2020)

The observational estimates adopted the standard power-law shape for the two-
point angular correlation function: $w(\theta)=A\theta^{1-\gamma}$. Since the
data did not allow an accurate determination of both $A$ and $\gamma$ for each
source population, they fixed $\gamma$ and fit for the normalization.

Figure 9: Two-point angular correlation function $w(\theta)$ yielded by our
simulation (symbols with error bars) for radio AGN (left) and SFGs (right)
compared to the observational results in Table 1 (coloured lines).

To produce the T-RECS correlation functions we use the same flux limits of the
observational estimates and a sky area of 4 $\times$4 square degrees. We use
the Hamilton (1993) estimator:

$w(\theta)=\frac{DD\cdot RR}{DR\cdot DR}-1,$ (15)

where $DD$, $RR$ and $DR$ are the number of data-data, random-random and data-
random pairs separated by $\theta$.

The random catalogue is constructed by allocating a random position of all
sources within the simulated sky area. The $w(\theta)$ is computed for 500
realizations. Each realization uses a different generation of the continuum
catalogue and a different 4$\times$4 deg2 portion of the DM lightcone, which
is extracted from the 5$\times$5 deg2 with a random uniform shift from the
field centre in both coordinates.

In Fig. 9 we show the mean observational estimates as lines of different
colours. Uncertainties have not been included, in order to improve figure
clarity. The different selection criteria are responsible for some of the
differences in the measured clustering. At the same selection frequency, a
deeper flux limit would correspond to less luminous sources, which are
expected to be found in less massive, weakly clustered haloes. Therefore,
deeper catalogues find a lower normalization of the correlation function. A
different selection frequency would favour different sub-populations, which
could have different clustering properties.

T-RECS mean values and their dispersion for catalogues with flux limits as in
Magliocchetti et al. (2017) and Hale et al. (2018) are shown as symbols and
error bars with the same colour as the observational estimates. For the same
source selection, our simulations give amplitudes of the angular correlation
function consistently lower than the observational ones.

The higher amplitudes of the observationally estimated $w(\theta)$ imply
higher average bias factors, i.e. higher halo masses. As already discussed in
T-RECS I, various factors can influence the comparison between the
observational determinations and the T-RECS simulations. One of those is the
observational classification between AGN and SFG and how this compares with
the definition of the two populations in T-RECS. Another factor is the
availability of high-mass haloes in the Baugh et al. (2019) cosmological
simulation adopted in T-RECS. The mass function of the cosmological simulation
and that inferred from the luminosity function are somewhat different, which
means that there is a deficit of suitable haloes in some mass ranges and a
surplus in others. Allowing for some scatter between the predicted and the
associated mass alleviates the problem; however, it typically favours
association to smaller halo masses, given the shape of the mass function.

T-RECS results are still consistent with the observations using the same flux
limit, with the exception of Hale et al. (2018) SFGs. In the latter case the
discrepancy is of 2.5 $\sigma$, where $\sigma$ is the quadratic sum of errors
of the observational estimates and of the simulations.

Despite the differences noted above, and taking into account the uncertainties
in both determinations, the agreement between the T-RECS clustering and the
empirically-determined one is reasonably good.

###  5.2 Validation for HI sources

Observational estimates of the real-space correlation functions for galaxies
selected in HI are presented by Passmoor et al. (2011), Martin et al. (2012)
and Papastergis et al. (2013). Details of the survey used are reported in
Table 2. The observational estimates adopt a power-law model for the
correlation function of the form $\xi(r)=(r/r_{0})^{-\gamma}$, and fit for
both the $r_{0}$ and the $\gamma$ parameters.

We select T-RECS sources with a mass limit of $\log M_{\rm HI}=7.5\,M_{\odot}$
and a redshift limit of $z=0.05$, in an attempt to be comparable with the
observational estimates. However, the observed samples have a complex
selection function and they are neither flux- or mass-limited, which makes the
comparison of the results more complicated.

In order to get good statistics for the computation of $\xi(r)$, we use the
maximum possible area for a T-RECS clustered simulation, of 25 deg2.
Observational estimates were obtained within larger surveys, of $\sim$400 deg2
for ALFALFA and 29,000 deg2 for HIPASS.

To produce T-RECS correlation functions we use again the estimator in eq.
(15); we checked that using the Landy & Szalay (1993) estimator as in Passmoor
et al. (2011), Martin et al. (2012) and Papastergis et al. (2013) gives
comparable results. For the random catalogue we use the unclustered T-RECS
catalogue, which therefore has a plausible redshift distribution and random
angular coordinates within the simulated sky area. The clustered catalogue
contains the same sources, but the angular and redshift coordinates are
replaced with those of dark haloes in the same redshift range that have
similar mass.

The $\xi(r)$ is computed for 500 realizations. Each realization uses a
different generation of both the HI and continuum catalogue (which is also
used in the dark mass matching process, as previously explained).

In Fig. 10 we show the comparison between T-RECS and the observational
estimates (uncertainties in the latter were not included to improve figure
clarity). T-RECS results reproduce the slope of the observations and are lower
in amplitude; the discrepancy is of 1.5–2.5 $\sigma$ depending on the dataset.

As previously mentioned, this comparison is made more difficult by the
different source selection. T-RECS sources are complete down to $M_{\rm
HI}=10^{7.5}$. This is the minimum mass achieved by the observational samples
(specifically, Papastergis et al., 2013), which could in part explain the
lower normalization. The HI catalogue would also most likely inherit the same
issue of the continuum catalogue, that has to do with the availability of
high-mass haloes in the cosmological simulation.

Table 2: Details of the samples used for HI clustering estimates in the literature. Values of mass denoted with * are average masses rather than lowest masses. Observation | $v_{\rm max}{\rm[km/s]}$ | $z_{\rm max}$ |  $M_{\rm HI,min}$ [$M_{\odot}$] | Ref.
---|---|---|---|---
ALFALFA | 15,000 | 0.05 | $10^{7.5}$ | Papastergis et al. (2013)
ALFALFA | 15,000 | 0.06 |  | Martin et al. (2012)
HIPASS | 12,700 | 0.06 | 3.24 $10^{9}$* | Passmoor et al. (2011)
ALFALFA | 18,000 | 0.06 | 2.48 $10^{9}$* | Passmoor et al. (2011)
Figure 10: Real-space correlation function $\xi(r)$ yielded by our simulation
(symbols with error bars) compared to the observational results in Table 2
(coloured lines).

##  6 Conclusions

This work updated the T-RECS radio continuum model presented in Bonaldi et al.
(2019) to improve the agreement with the most recent observations and extended
it to include HI emission from extra-galactic sources as well as the cross-
correlation between the two.

The HI model reproduces the Jones et al. (2018) HI mass function and the Oman
(2022) HI velocity width function at $z=0$. Evolution of the HI mass function
has been modelled as a linear dependence of the $\log(M^{*}_{\rm HI})$ and
$\log\phi\star$ parameters with redshift, with amplitudes calibrated to give
results consistent with Paul et al. (2023) at $z=0.3$. Uncertainties in this
modelling are due to a lack of high-redshift data and currently limit our HI
simulation to a maximum redshift of 0.5. The HI model further includes source
size (modelled with the $\log M_{\rm HI}$–$\log D_{\rm HI}$ relation by
Naluminsa et al., 2021) and morphology (galaxy type, inclination, ellipticity
using Holmberg 1958 and Rodríguez & Padilla 2013).

The updated radio continuum model makes use of a steepening synchrotron
spectral index and the Smith et al. (2021) $L({\rm SFR})$–$M_{\rm star}$
relation. The goodness of the model has been evaluated by comparing T-RECS
realizations with observed luminosity functions at 1.4 GHz as well as the most
recent compilations of number counts in the 150 MHz–15 GHz. We also modelled
an HI mass proxy for the continuum sources, $\tilde{M}_{\rm HI}$, by using the
SFR–$M_{\rm HI}$ correlation for star-forming galaxies and published relations
between the HI mass and the stellar mass for the AGN.

We showed that, by matching the HI mass proxy in the continuum catalogue with
the HI mass in the HI catalogue, we are able to assign plausible counterparts
between the otherwise independent T-RECS continuum and HI simulations. Our
method ensures a correct number of counterparts by applying an HI selection
function to the radio continuum sources consistent with that of the HI
catalogue, within an error $\sigma_{\log M}\sim 0.65$. It also propagates
modelled correlations between the radio continuum and HI properties via their
link to $M_{\rm HI}$.

The $M_{\rm h}$–$M_{\rm HI}$ relation for T-RECS galaxies is qualitatively in
good agreement with the HI–halo mass relation from data (Guo et al., 2020;
Rhee et al., 2023) and theoretical modelling (Chauhan et al., 2021).

Clustering properties were simulated by associating the galaxies with haloes
of the P-millennium (Baugh et al., 2019) simulation of similar dark halo mass.
We assessed the clustering properties of both the HI and the continuum
catalogue and validated them against the empirical estimates of $w(\theta)$
(Magliocchetti et al., 2017; Bonato et al., 2021; Chakraborty et al., 2020)
and $\xi(r)$ (Passmoor et al., 2011; Martin et al., 2012; Papastergis et al.,
2013). The T-RECS simulation correctly reproduces the slope of both
correlation functions. The normalization is in both cases lower (of 1.5–2.5
$\sigma$). This could be explained at least in part by the difficulty in
reproducing the selection functions of the observational estimates; issues
with the availability of high-mass haloes in the Baugh et al. (2019) dark
matter simulation we derived the clustering from could also be responsible.

The code used to produce these results is available on github
(https://github.com/abonaldi/TRECS.git) and can be used to produce mock
observations of radio continuum and HI observations, as well as their cross-
correlation. We believe these can be very useful to plan for future surveys.
This code has been used as a basis of the second science data challenge (SDC2
Hartley et al., 2023) organised by the SKA Observatory to prepare the science
community to analyse SKA HI data.

## Acknowledgements

We thank M. Massardi, V. Galluzzi, A. Lapi, I. Prandoni for useful suggestions
on how to update the T-RECS continuum model. We thank S. Blyth for her useful
comments on the paper. We also thank the anonymous referee for their review,
that has led to a substantial improvement of this work.

We acknowledge the use of computational resources at the SKA Observatory.

TR is supported by the PRIN MIUR 2017 prot. 20173ML3WW, ’Opening the ALMA
window on the cosmic evolution of gas, stars and supermassive black holes’,
and by the Fondazione ICSC - Spoke 3 Astrophysics and Cosmos Observations -
National Recovery and Resilience Plan Project ID CN-00000013 ’Italian Research
Center on High-Performance Computing, Big Data and Quantum Computing’ - Next
Generation EU.

We acknowledge usage of the Fortran (Backus & Heising, 1964) and Python (Van
Rossum et al., 2007) programming languages, Astropy (The Astropy Collaboration
et al., 2013), NumPy (Harris et al., 2020), Scikit-learn (Pedregosa et al.,
2011), LAPACK (Anderson et al., 1999), GSL (Galassi et al., 1996), HEALPix
(Gorski et al., 2005), CFitsIO (Pence, 2010).

## Data availability

The data underlying this article are available at the link
https://tinyurl.com/TRECS2222The full link is
https://www.dropbox.com/sh/52pbxscr8pkrdn7/AAA-RPe6QcTeQrcqEvNpSWkoa?dl=0. The
code used to generate the data is available on Gihub
(https://github.com/abonaldi/TRECS.git)

## References

  * Adams et al. (2022) Adams, E. A. K., Adebahr, B., de Blok, W. J. G., et al., 2022, A&A, 667, A38, arXiv:2208.05348
  * AMI Consortium et al. (2011) AMI Consortium, Davies, M. L., Franzen, T. M. O., et al., 2011, MNRAS, 415, 2708, arXiv:1012.3659
  * Anderson et al. (1999) Anderson, E., Bai, Z., Bischof, C., et al., 1999, LAPACK users’ guide, SIAM
  * Aversa et al. (2015) Aversa, R., Lapi, A., de Zotti, G., Shankar, F., Danese, L., 2015, ApJ, 810, 74, arXiv:1507.07318
  * Backus & Heising (1964) Backus, J. W., Heising, W. P., 1964, IEEE Transactions on Electronic Computers, , 4, 382
  * Barnes et al. (2001) Barnes, D. G., Staveley-Smith, L., de Blok, W. J. G., et al., 2001, MNRAS, 322, 3, 486
  * Barnes & Haehnelt (2014) Barnes, L. A., Haehnelt, M. G., 2014, MNRAS, 440, 3, 2313, arXiv:1403.1873
  * Baugh et al. (2019) Baugh, C. M., Gonzalez-Perez, V., Lagos, C. d. P., et al., 2019, MNRAS, 483, 4, 4922, arXiv:1808.08276
  * Bera et al. (2022) Bera, A., Kanekar, N., Chengalur, J. N., Bagla, J. S., 2022, ApJ, 940, 1, L10, arXiv:2210.17018
  * Blyth et al. (2016) Blyth, S., Baker, A. J., Holwerda, B., et al., 2016, in MeerKAT Science: On the Pathway to the SKA, 4
  * Bonaldi et al. (2021) Bonaldi, A., An, T., Brüggen, M., et al., 2021, MNRAS, 500, 3, 3821, arXiv:2009.13346
  * Bonaldi et al. (2019) Bonaldi, A., Bonato, M., Galluzzi, V., et al., 2019, MNRAS, 482, 1, 2, arXiv:1805.05222
  * Bonato et al. (2017) Bonato, M., Negrello, M., Mancuso, C., et al., 2017, MNRAS, 469, 2, 1912, arXiv:1704.05459
  * Bonato et al. (2021) Bonato, M., Prandoni, I., De Zotti, G., Brienza, M., Morganti, R., Vaccari, M., 2021, MNRAS, 500, 1, 22, arXiv:2010.08748
  * Bondi et al. (2008) Bondi, M., Ciliegi, P., Schinnerer, E., et al., 2008, ApJ, 681, 1129, arXiv:0804.1706
  * Bothwell et al. (2013) Bothwell, M. S., Maiolino, R., Kennicutt, R., et al., 2013, MNRAS, 433, 2, 1425, arXiv:1304.4940
  * Broeils & Rhee (1997) Broeils, A. H., Rhee, M. H., 1997, A&A, 324, 877
  * Butler et al. (2019) Butler, A., Huynh, M., Kapińska, A., et al., 2019, A&A, 625, A111, arXiv:1904.12507
  * Catinella et al. (2010) Catinella, B., Schiminovich, D., Kauffmann, G., et al., 2010, MNRAS, 403, 2, 683, arXiv:0912.1610
  * Chakraborty et al. (2020) Chakraborty, A., Dutta, P., Datta, A., Roy, N., 2020, MNRAS, 494, 3, 3392, arXiv:2002.12383
  * Chauhan et al. (2021) Chauhan, G., Lagos, C. d. P., Stevens, A. R. H., et al., 2021, MNRAS, 506, 4, 4893, arXiv:2102.12203
  * Cluver et al. (2014) Cluver, M. E., Jarrett, T. H., Hopkins, A. M., et al., 2014, ApJ, 782, 2, 90, arXiv:1401.0837
  * Cortese et al. (2011) Cortese, L., Catinella, B., Boissier, S., Boselli, A., Heinis, S., 2011, MNRAS, 415, 2, 1797, arXiv:1103.5889
  * De Zotti et al. (2005) De Zotti, G., Ricci, R., Mesa, D., et al., 2005, A&A, 431, 893, astro-ph/0410709
  * DiPompeo et al. (2013) DiPompeo, M. A., Runnoe, J. C., Myers, A. D., Boroson, T. A., 2013, ApJ, 774, 24, arXiv:1307.1731
  * Duffy et al. (2012) Duffy, A. R., Meyer, M. J., Staveley-Smith, L., et al., 2012, MNRAS, 426, 4, 3385, ISSN 0035-8711, https://academic.oup.com/mnras/article-pdf/426/4/3385/3338965/426-4-3385.pdf
  * Fanaroff & Riley (1974) Fanaroff, B. L., Riley, J. M., 1974, MNRAS, 167, 31P
  * Fernández et al. (2016) Fernández, X., Gim, H. B., van Gorkom, J. H., et al., 2016, ApJ, 824, 1, L1, arXiv:1606.00013
  * Franzen et al. (2016) Franzen, T. M. O., Jackson, C. A., Offringa, A. R., et al., 2016, MNRAS, 459, 3314, arXiv:1604.03751
  * Galassi et al. (1996) Galassi, M., Davies, J., Theiler, J., et al., 1996, No. Release, 2
  * Giovanelli et al. (2005) Giovanelli, R., Haynes, M. P., Kent, B. R., et al., 2005, AJ, 130, 6, 2598, arXiv:astro-ph/0508301
  * Gogate et al. (2023) Gogate, A. R., Verheijen, M. A. W., van der Hulst, J. M., Jaffé, Y. L., 2023, MNRAS, 519, 3, 4279, arXiv:2303.06618
  * Gorski et al. (2005) Gorski, K. M., Hivon, E., Banday, A. J., et al., 2005, The Astrophysical Journal, 622, 2, 759
  * Guo et al. (2020) Guo, H., Jones, M. G., Haynes, M. P., Fu, J., 2020, ApJ, 894, 2, 92, arXiv:2004.04762
  * Hale et al. (2018) Hale, C. L., Jarvis, M. J., Delvecchio, I., et al., 2018, MNRAS, 474, 4133, arXiv:1711.05201
  * Hamilton (1993) Hamilton, A. J. S., 1993, ApJ, 417, 19
  * Harris et al. (2020) Harris, C. R., Millman, K. J., Van Der Walt, S. J., et al., 2020, Nature, 585, 7825, 357
  * Hartley et al. (2023) Hartley, P., Bonaldi, A., Braun, R., et al., 2023, MNRAS, arXiv:2303.07943
  * Hartley et al. (2019) Hartley, P., Jackson, N., Sluse, D., Stacey, H. R., Vives-Arias, H., 2019, MNRAS, 485, 3, 3009, arXiv:1901.05791
  * Haynes et al. (2018) Haynes, M. P., Giovanelli, R., Kent, B. R., et al., 2018, ApJ, 861, 1, 49, arXiv:1805.11499
  * Haynes et al. (2011) Haynes, M. P., Giovanelli, R., Martin, A. M., et al., 2011, AJ, 142, 5, 170, arXiv:1109.0027
  * Hellwig et al. (1970) Hellwig, H., Vessot, R. F. C., Levine, M. W., Zitzewitz, P. W., Allan, D. W., Glaze, D. J., 1970, IEEE Transactions on Instrumentation Measurement, 19, 4, 200
  * Holmberg (1958) Holmberg, E., 1958, Meddelanden fran Lunds Astronomiska Observatorium Serie II, 136, 1
  * Huang et al. (2012) Huang, S., Haynes, M. P., Giovanelli, R., Brinchmann, J., 2012, ApJ, 756, 2, 113, arXiv:1207.0523
  * Ibar et al. (2009) Ibar, E., Ivison, R. J., Biggs, A. D., Lal, D. V., Best, P. N., Green, D. A., 2009, MNRAS, 397, 281, arXiv:0903.3600
  * Janssen et al. (2012) Janssen, R. M. J., Röttgering, H. J. A., Best, P. N., Brinchmann, J., 2012, A&A, 541, A62, arXiv:1206.0578
  * Jaskot et al. (2015) Jaskot, A. E., Oey, M. S., Salzer, J. J., Van Sistine, A., Bell, E. F., Haynes, M. P., 2015, ApJ, 808, 1, 66, arXiv:1506.05081
  * Jones et al. (2018) Jones, M. G., Haynes, M. P., Giovanelli, R., Moorman, C., 2018, MNRAS, 477, 1, 2, arXiv:1802.00053
  * Katz et al. (2019) Katz, H., Desmond, H., McGaugh, S., Lelli, F., 2019, MNRAS, 483, 1, L98, arXiv:1810.12347
  * Kellermann et al. (2016) Kellermann, K. I., Condon, J. J., Kimball, A. E., Perley, R. A., Ivezić, Ž., 2016, ApJ, 831, 168, arXiv:1608.04586
  * Kennicutt (1998) Kennicutt, J., Robert C., 1998, ApJ, 498, 2, 541, arXiv:astro-ph/9712213
  * Kim et al. (2017) Kim, H.-S., Wyithe, J. S. B., Baugh, C. M., Lagos, C. d. P., Power, C., Park, J., 2017, MNRAS, 465, 1, 111, arXiv:1603.02383
  * Koribalski et al. (2020) Koribalski, B. S., Staveley-Smith, L., Westmeier, T., et al., 2020, Ap&SS, 365, 7, 118, arXiv:2002.07311
  * Kozieł-Wierzbowska et al. (2020) Kozieł-Wierzbowska, D., Goyal, A., Żywucka, N., 2020, ApJS, 247, 2, 53, arXiv:1912.09959
  * Landy & Szalay (1993) Landy, S. D., Szalay, A. S., 1993, ApJ, 412, 64
  * Lin et al. (2010) Lin, Y.-T., Shen, Y., Strauss, M. A., Richards, G. T., Lunnan, R., 2010, ApJ, 723, 1119, arXiv:1006.5452
  * Maddox et al. (2021) Maddox, N., Frank, B. S., Ponomareva, A. A., et al., 2021, A&A, 646, A35, arXiv:2011.09470
  * Maddox et al. (2015) Maddox, N., Hess, K. M., Obreschkow, D., Jarvis, M. J., Blyth, S. L., 2015, MNRAS, 447, 2, 1610, arXiv:1412.0852
  * Maddox et al. (2016) Maddox, N., Jarvis, M. J., Oosterloo, T. A., 2016, MNRAS, 460, 4, 3419, arXiv:1605.04962
  * Magliocchetti et al. (2017) Magliocchetti, M., Popesso, P., Brusa, M., et al., 2017, MNRAS, 464, 3271, arXiv:1606.08286
  * Mancuso et al. (2015) Mancuso, C., Lapi, A., Cai, Z.-Y., et al., 2015, ApJ, 810, 72, arXiv:1507.08411
  * Mancuso et al. (2017) Mancuso, C., Lapi, A., Prandoni, I., et al., 2017, ApJ, 842, 95, arXiv:1705.06539
  * Mandal et al. (2021) Mandal, S., Prandoni, I., Hardcastle, M. J., et al., 2021, A&A, 648, A5, arXiv:2011.08829
  * Martin et al. (2012) Martin, A. M., Giovanelli, R., Haynes, M. P., Guzzo, L., 2012, ApJ, 750, 1, 38, arXiv:1202.6005
  * Massardi et al. (2010) Massardi, M., Bonaldi, A., Negrello, M., Ricciardi, S., Raccanelli, A., de Zotti, G., 2010, MNRAS, 404, 532, arXiv:1001.1069
  * Meyer (2009) Meyer, M., 2009, in Panoramic Radio Astronomy: Wide-field 1-2 GHz Research on Galaxy Evolution, 15, arXiv:0912.2167
  * Michałowski et al. (2015) Michałowski, M. J., Gentile, G., Hjorth, J., et al., 2015, A&A, 582, A78, arXiv:1508.03094
  * Naluminsa et al. (2021) Naluminsa, E., Elson, E. C., Jarrett, T. H., 2021, MNRAS, 502, 4, 5711, arXiv:2101.04514
  * Noordermeer et al. (2005) Noordermeer, E., van der Hulst, J. M., Sancisi, R., Swaters, R. A., van Albada, T. S., 2005, A&A, 442, 1, 137, arXiv:astro-ph/0508319
  * Novak et al. (2017) Novak, M., Smolčić, V., Delhaize, J., et al., 2017, A&A, 602, A5, arXiv:1703.09724
  * Obuljen et al. (2019) Obuljen, A., Alonso, D., Villaescusa-Navarro, F., Yoon, I., Jones, M., 2019, MNRAS, 486, 4, 5124, arXiv:1805.00934
  * Ocran et al. (2020) Ocran, E. F., Taylor, A. R., Vaccari, M., et al., 2020, MNRAS, 491, 4, 5911, arXiv:1912.00934
  * Oh et al. (2015) Oh, S.-H., Hunter, D. A., Brinks, E., et al., 2015, AJ, 149, 6, 180, arXiv:1502.01281
  * Oman (2022) Oman, K. A., 2022, MNRAS, 509, 3, 3268, arXiv:2108.08856
  * Padovani et al. (2015) Padovani, P., Bonzini, M., Kellermann, K. I., Miller, N., Mainieri, V., Tozzi, P., 2015, MNRAS, 452, 1263, arXiv:1506.06554
  * Papastergis et al. (2013) Papastergis, E., Giovanelli, R., Haynes, M. P., Rodríguez-Puebla, A., Jones, M. G., 2013, ApJ, 776, 1, 43, arXiv:1308.2661
  * Passmoor et al. (2011) Passmoor, S. S., Cress, C. M., Faltenbacher, A., 2011, MNRAS, 412, 1, L50
  * Paul et al. (2018) Paul, N., Choudhury, T. R., Paranjape, A., 2018, MNRAS, 479, 2, 1627, arXiv:1712.04469
  * Paul et al. (2023) Paul, S., Santos, M. G., Chen, Z., Wolz, L., 2023, arXiv e-prints, arXiv:2301.11943, arXiv:2301.11943
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al., 2011, the Journal of machine Learning research, 12, 2825
  * Pence (2010) Pence, W. D., 2010, Astrophysics Source Code Library, ascl–1010
  * Planck Collaboration VI (2020) Planck Collaboration VI, 2020, A&A, 641, A6, arXiv:1807.06209
  * Planck Collaboration XVI (2014) Planck Collaboration XVI, 2014, A&A, 571, A16
  * Ponomareva et al. (2023) Ponomareva, A. A., Jarvis, M. J., Pan, H., et al., 2023, arXiv e-prints, arXiv:2304.13051, arXiv:2304.13051
  * Ponomareva et al. (2016) Ponomareva, A. A., Verheijen, M. A. W., Bosma, A., 2016, MNRAS, 463, 4, 4052, arXiv:1609.00378
  * Popping et al. (2014) Popping, G., Somerville, R. S., Trager, S. C., 2014, MNRAS, 442, 3, 2398, arXiv:1308.6764
  * Rajohnson et al. (2022) Rajohnson, S. H. A., Frank, B. S., Ponomareva, A. A., et al., 2022, MNRAS, 512, 2, 2697, arXiv:2203.06149
  * Rhee et al. (2023) Rhee, J., Meyer, M., Popping, A., et al., 2023, MNRAS, 518, 3, 4646, arXiv:2210.09697
  * Rodríguez & Padilla (2013) Rodríguez, S., Padilla, N. D., 2013, MNRAS, 434, 3, 2153, arXiv:1306.3264
  * Siewert et al. (2019) Siewert, T. M., Hale, C., Bhardwaj, N., et al., 2019, arXiv e-prints, arXiv:1908.10309, arXiv:1908.10309
  * Sinigaglia et al. (2022) Sinigaglia, F., Rodighiero, G., Elson, E., et al., 2022, ApJ, 935, 1, L13, arXiv:2208.01121
  * Smith et al. (2021) Smith, D. J. B., Haskell, P., Gürkan, G., et al., 2021, A&A, 648, A6, arXiv:2011.08196
  * Smolčić et al. (2017) Smolčić, V., Novak, M., Bondi, M., et al., 2017, A&A, 602, A1, arXiv:1703.09713
  * Staveley-Smith & Oosterloo (2015) Staveley-Smith, L., Oosterloo, T., 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14), 167, arXiv:1506.04473
  * Swaters et al. (2009) Swaters, R. A., Sancisi, R., van Albada, T. S., van der Hulst, J. M., 2009, A&A, 493, 3, 871, arXiv:0901.4222
  * Swaters et al. (2002) Swaters, R. A., van Albada, T. S., van der Hulst, J. M., Sancisi, R., 2002, A&A, 390, 829, arXiv:astro-ph/0204525
  * The Astropy Collaboration et al. (2013) The Astropy Collaboration, Robitaille, Thomas P., Tollerud, Erik J., et al., 2013, A&A, 558, A33
  * Van Rossum et al. (2007) Van Rossum, G., et al., 2007, in USENIX annual technical conference, vol. 41, 1–36, Santa Clara, CA
  * Verheijen & Sancisi (2001) Verheijen, M. A. W., Sancisi, R., 2001, A&A, 370, 765, arXiv:astro-ph/0101404
  * Vernstrom et al. (2016) Vernstrom, T., Scott, D., Wall, J. V., et al., 2016, MNRAS, 462, 2934, arXiv:1603.03085
  * Villaescusa-Navarro et al. (2018) Villaescusa-Navarro, F., Genel, S., Castorina, E., et al., 2018, ApJ, 866, 2, 135, arXiv:1804.09180
  * Waldram et al. (2010) Waldram, E. M., Pooley, G. G., Davies, M. L., Grainge, K. J. B., Scott, P. F., 2010, MNRAS, 404, 1005, arXiv:0908.0066
  * Walter et al. (2020) Walter, F., Carilli, C., Neeleman, M., et al., 2020, ApJ, 902, 2, 111, arXiv:2009.11126
  * Wang et al. (2014) Wang, J., Fu, J., Aumer, M., et al., 2014, MNRAS, 441, 3, 2159, arXiv:1401.8164
  * Wang et al. (2016) Wang, J., Koribalski, B. S., Serra, P., et al., 2016, MNRAS, 460, 2, 2143, arXiv:1605.01489
  * White et al. (2017) White, S. V., Jarvis, M. J., Kalfountzou, E., et al., 2017, MNRAS, 468, 217, arXiv:1702.00904
  * Whittam et al. (2016) Whittam, I. H., Riley, J. M., Green, D. A., et al., 2016, MNRAS, 457, 1496, arXiv:1601.00282
  * Wong et al. (2006) Wong, O. I., Ryan-Weber, E. V., Garcia-Appadoo, D. A., et al., 2006, MNRAS, 371, 4, 1855, arXiv:astro-ph/0607491
  * Xi et al. (2021) Xi, H., Staveley-Smith, L., For, B.-Q., et al., 2021, MNRAS, 501, 3, 4550, arXiv:2012.09516
  * Zhou et al. (2018) Zhou, Z., Wu, H., Zhou, X., Ma, J., 2018, PASP, 130, 991, 094101, arXiv:1807.09955

##  Appendix A Content of the catalogues

The structure of each T-RECS HI and continuum catalogues is presented in Table
A1 and A2, respectively. The continuum catalogues include the AGN and SFG
populations. A few of the modelled quantities apply to only one of the
populations and are left blank for the other. A flag value of -100 is used in
this case.

HI $\times$ continuum catalogues contain all columns in Table A1 and Table A2
with the exception of the second set of coordinates (x_coord, y_coord,
latitude, longitude, redshift) which are rewritten after plausible
counterparts between the two catalogues have been identified. Other quantities
that have the same name in both catalogues are both kept, and the second set
is renamed (e.g. bmaj_1 instead of bmaj). This is because they all contain
unique information, (e.g., the size as traced by the stars and by the neutral
gas).

Table 3: Structure of the HI catalogues released with this paper. Catalogues produced with the T-RECS code will have the same format. Tag Name |  Units |  Description
---|---|---
ID_HI |  |  Numerical identifier for the source
MHI |  log($M_{\odot}$) |  HI mass
HI flux |  mJY Hz |  HI flux
Mh |  log($M_{\odot}$) |  Dark halo mass
Mstar |  log($M_{\odot}$) |  Stellar mass
x_coord |  degs |  First angular coordinate in the flat-sky approximation (see T-RECS I for more details)
y_coord |  degs |  Second angular coordinate in the flat-sky approximation (see T-RECS I for more details)
latitude |  degs |  Latitude spherical coordinate for a chosen centre of the field
longitude |  degs |  Longitude spherical coordinate for a chosen centre of the field
redshift |  |  redshift
HI size |  arcsec |  Apparent size of the HI source
inclination |  degs |  galaxy inclination
axis ratio |  |  galaxy axis ratio, defined as the ratio between major and minor axis
bmaj |  arcsec |  major axis
bmin |  arcsec |  minor axis
PA |  degs |  position angle
OptClass |  |  Number identifying the optical type: 1 for elliptical, 2 for spiral
Table 4: Structure of the radio continuum catalogues released with this paper. Catalogues produced with the T-RECS code will have the same format except for the number and list of frequencies and the possibility to optionally output the luminosities for each frequency as additional columns. Tag Name |  Units |  Description
---|---|---
ID_cont |  |  Numerical identifier for the source
Lum1400 |  log(erg/s/Hz) |  Luminosity at 1.4 GHz (AGN only)
logSFR |  log($M_{\rm sun}$)/yr |  SFR (SFG only)
Ifreq |  mJy |  Total intensity flux density of the source at frequency _freq_ for each of a list of $N_{\rm freqs}$ frequencies as specified by the user.
Pfreq |  mJy |  Polarized flux density of the source at frequency _freq_ for each of a list of $N_{\rm freqs}$ frequencies as specified by the user. See T-RECS I for details of the polarization model.
Mh |  log($M_{\rm sun}$) |  Dark halo mass
Mstar |  log($M_{\rm sun}$) |  Stellar mass
MHI_pred |  log($M_{\rm sun}$) |  HI mass proxy
x_coord |  degs |  First angular coordinate in the flat-sky approximation (see T-RECS I for more details)
y_coord |  degs |  Second angular coordinate in the flat-sky approximation (see T-RECS I for more details)
latitude |  degs |  Latitude spherical coordinate for a chosen centre of the field
longitude |  degs |  Longitude spherical coordinate for a chosen centre of the field
redshift |  |  redshift
size |  arcsec |  Apparent size of the source. This is the maximum size of the core+jet emission for AGN, and the scale radius of a Sersic profile for SFG.
inclination |  degs |  galaxy inclination
axis ratio |  |  galaxy axis ratio, defined as the ratio between major and minor axis
bmaj |  arcsec |  major axis
bmin |  arcsec |  minor axis
PA |  degs |  position angle
Rs |  |  Ratio between the distance between the spots and the total size of the jets, for the FR I /FR II classification (steep-spectrum AGN only)
RadioClass |  |  Number identifying the sub-population: 1 for late-type; 2 for spheroids; 3 for lensed spheroids; 4 for FSRQ, 5 for BL Lac and 6 for SS-AGN.
OptClass |  |  Number identifying the optical type: 1 for elliptical, 2 for spiral.
Lfreq |  erg/s/Hz |  OPTIONAL OUTPUT: Luminosity of the source at frequency _freq_ for each of a list of $N_{\rm freqs}$ frequencies as specified by the user.

##  Appendix B Running the complete T-RECS pipeline

Extending T-RECS to include HI and HI-continuum cross-correlation has meant an
increase in the code complexity and an evolution in its architecture. Along
with these modifications, we have implemented a new build system and a user
interface we detail in this appendix.

When building the library, separate FORTRAN modules are compiled, each of
which deals with a different modular component of the work-flow. The work-flow
is controlled by a master script, accessed as terminal command, which executes
calls to the different components. The modules called by the master script are
listed below:

  * •

continuum sampler (trecs_sampler_continuum) produces a sample of star forming
galaxies and AGN to whom a dark mass and a HI mass are associated in order to
allow for cross-matching with clustering and HI properties, as detailed in
Sec. 3;

  * •

HI sampler (trecs_sampler_hi) produces a sample of HI galaxies, with the total
HI mass, line width, and morphological properties described in Sec. 2; this
catalogue is still un-correlated to the one produced by the continuum module;

  * •

cross-matching of HI properties with continuum properties (trecs_xmatch_hi)
implements the algorithm matching continuum and HI sources, detailed in Sec.
4. It is used when the outputs of the two modules have to be used as a
consistent set on the same sky area. Since the input of this step is given by
the HI and continuum catalogues, these are to be computed in advance.

  * •

cross-match of the DM-halo masses with those of the P-Millenium light-cone
(Sec. 5) delivered with the input data-set coupled with the library
(trecs_xmatch_clustering). The coordinates of haloes in the P-Millenium (Baugh
et al., 2019) light-cone embed the spatial properties of dark matter haloes
and sub-haloes in a Planck cosmology. If for a source in the T-RECS catalogue
a DM-halo counter-part is identified in the light-cone, the sky coordinates of
the source are overwritten with those of the associated halo. This cross-match
can be applied to each output catalogue of the T-RECS samplers.

  * •

wrapping of the produced catalogues in just one file (trecs_wrapper). For
parallelisation and RAM-usage purposes, the T-RECS functions are launched and
operate on redshift slices that are stored in separate files. The list of
redshift slices is reported in table B1. This final module collates the raw
outputs into a single catalogue file, rotating the field of view towards user-
defined coordinates.

All of the components listed above can be run separately and are independent
one from the other with some obvious exceptions (e.g. cross-matching, whether
it is between continuum and HI properties or between the sampled catalogues
and the DM haloes in the light-cone, requires to have first produced the
relevant catalogues to match).

Running the full code requires to provide a parameter file with user-defined
properties that will apply to the whole run (e.g. field of view, redshift
limits, what output catalogues to cross-match). By changing the initial random
seed one can either produce different realisations or reproduce the exact
results of some previous run.

As already mentioned, all of the functionalities of the library are wrapped by
a command line interface that can be launched by running the

    $ trecs [OPTIONAL FLAGS]
        -p|--params /path/to/parameter_file.ini

command. The optional flags control the overall behaviour and which module(s)
to run. An example working call to the command would be:

    $ trecs --continuum --HI
        --xmatch --clustering
        --wrap HI_continuum
        --params /path/to/parameter_file.ini

the command executes continuum sampling and HI sampling then it cross-matches
the two catalogues and assigns coordinates in the light-cone to as many
sources as possible; finally, the raw catalogues are wrapped into one single
file with sky coordinates rotated to match the user-defined field of view.

Table 5: Redshift slices adopted in T-RECS Central $z$ | Minimum $z$ | Maximum $z$
---|---|---
0.01 | 0.00 | 0.01
0.02 | 0.01 | 0.03
0.05 | 0.03 | 0.07
0.10 | 0.07 | 0.12
0.15 | 0.12 | 0.17
0.20 | 0.17 | 0.22
0.25 | 0.22 | 0.27
0.30 | 0.27 | 0.32
0.35 | 0.32 | 0.37
0.40 | 0.37 | 0.42
0.45 | 0.42 | 0.47
0.50 | 0.47 | 0.52
0.55 | 0.52 | 0.57
0.60 | 0.57 | 0.62
0.65 | 0.62 | 0.67
0.70 | 0.67 | 0.72
0.75 | 0.72 | 0.77
0.80 | 0.77 | 0.82
0.85 | 0.82 | 0.87
0.90 | 0.87 | 0.92
0.95 | 0.92 | 0.97
1.00 | 0.97 | 1.10
1.20 | 1.10 | 1.30
1.40 | 1.30 | 1.50
1.60 | 1.50 | 1.70
1.80 | 1.70 | 1.90
2.00 | 1.90 | 2.10
2.20 | 2.10 | 2.30
2.40 | 2.30 | 2.50
2.60 | 2.50 | 2.70
2.80 | 2.70 | 2.90
3.00 | 2.90 | 3.10
3.20 | 3.10 | 3.30
3.40 | 3.30 | 3.50
3.60 | 3.50 | 3.70
3.80 | 3.70 | 3.90
4.00 | 3.90 | 4.10
4.20 | 4.10 | 4.30
4.40 | 4.30 | 4.50
4.60 | 4.50 | 4.70
4.80 | 4.70 | 4.90
5.00 | 4.90 | 5.10
5.20 | 5.10 | 5.30
5.40 | 5.30 | 5.50
5.60 | 5.50 | 5.70
5.80 | 5.70 | 5.90
6.00 | 5.90 | 6.10
6.20 | 6.10 | 6.30
6.40 | 6.30 | 6.50
6.60 | 6.50 | 6.70
6.80 | 6.70 | 6.90
7.00 | 6.90 | 7.10
7.20 | 7.10 | 7.30
7.40 | 7.30 | 7.50
7.60 | 7.50 | 7.70
7.80 | 7.70 | 7.90
8.00 | 7.90 | 8.10
Table 6: Optional flags of the T-RECS terminal command (trecs) and the different type of simulation they produce if added. The [TAG] argument in the last column determines which raw catalogue to wrap (e.g. HI for the output of the HI sampler, HI_continuum for the cross-matched catalogue). The pipe symbol (|) in the column labels specifies that the flag can be called in two equivalent ways. Simulation type | \--HI|-i | \--continuum|-c | \--xmatch|-x | \--clustering|-C | \--wrap|-w [TAG]
---|---|---|---|---|---
HI, no clustering | X |  |  |  | HI
HI, with clustering | X |  |  | X | HI_clustering
Continuum, no clustering |  | X |  |  | continuum
Continuum, with clustering |  | X |  | X | continuum_clustering
Continuum $\times$ HI, no clustering | X | X | X |  | HI_continuum
Continuum $\times$ HI, with clustering | X | X | X | X | HI_continuum_clustering

Tab. 6 summarizes the available command line options and the different types
of simulations that are generated by calling the trecs command with different
combinations of arguments. A manual page with an explanation of all available
optional flags is shown on screen by calling the trecs command with the
-h|--help flag.