3.1. Types of galaxies in GHOST

GHOST includes four observational families that dominate the extragalactic radio sky at metre–centimetre wavelengths. Radio AGN appear in three classes – FSRQ, BLLac, and SSAGN – which differ in beaming/orientation and spectral behaviour but are often unified as the same central engine viewed at different angles (Urry & Padovani Reference Urry and Padovani1995; Heckman & Best Reference Heckman and Best2014). SFGs form a separate branch whose radio output tracks recent star formation. Below we summarise the properties and the labels used in this work.

• • FSRQ: Flat-spectrum radio quasars show core-dominated, compact radio emission with flat spectra (assumed mean $\alpha=-0.1$ Footnote ^b ), bright optical continua, and broad emission lines (Urry & Padovani Reference Urry and Padovani1995). In orientation-based unification they are relativistic jets viewed close to the line of sight, hence they have Doppler-boosted luminosities and are compact.

• • BLLac: Blazars with weak/absent emission lines and flat radio spectra, likewise interpreted as beamed (Doppler-boosted) jets at small inclination; compared to FSRQ they typically have lower line equivalent widths and similarly compact morphologies (Urry & Padovani Reference Urry and Padovani1995).

• • SSAGN: Steep-spectrum radio AGN, representing the bulk of radio-loud AGN, with mean spectral index $\zeta=-0.8$ . This class includes spectroscopic and morphological sub-taxonomies:

  1. – LERG: Jet-dominated systems with radiatively inefficient, low-Eddington accretion and weak/absent optical lines.

  2. – HERG: Radiatively efficient (‘radiative-mode’) accretors with strong optical emission lines (Heckman & Best Reference Heckman and Best2014; Bonaldi et al. Reference Bonaldi, Bonato, Galluzzi, Harrison, Massardi, Kay, De Zotti and Brown2019).

  3. – FR I: Edge-darkened morphology (brightness peaks nearer the core; jets fade outwards), more common at lower radio power (Fanaroff & Riley Reference Fanaroff and Riley1974).

  4. – FR II: Edge-brightened lobes with terminal hotspots, more common at higher radio power (Fanaroff & Riley Reference Fanaroff and Riley1974).

These schemes probe different physics and are not one-to-one (e.g. FR II LERG exist). In our catalogue the FR tag is assigned via the size-ratio parameter $R_s$ (Section 3.3), defined as the separation of peak-brightness regions divided by the total source extent; $R_s\gt0.5$ is labelled FR II and $R_s\lt0.5$ FR I.

• • SFG: Radio emission powered predominantly by recent star formation: non-thermal synchrotron from cosmic-ray electrons accelerated in supernova remnants plus thermal free–free from H ii regions (Condon Reference Condon1992; Murphy et al. Reference Murphy, Chary, Dickinson, Pope, Frayer and Lin2011). At GHz frequencies the ensemble spectral index is steep (typ. $\zeta\!\approx\!-0.7$ to $-0.8$ ) with a mild flattening from free–free towards higher frequencies. In GHOST the SFG branch is modelled with a calibrated RLF at 1.4 GHz and class-appropriate K-corrections; it is treated separately from AGN because its radio output traces star-formation physics rather than black-hole accretion.

Radio AGN typically inhabit more massive halos than SFG and therefore have a higher large-scale bias; within AGN, radiatively inefficient/jet-dominated systems tend to occupy the most massive environments, while radiatively efficient (quasar-like) systems are somewhat less clustered. Conversely, radio-selected SFG are associated with lower-mass halos and exhibit a lower bias. These trends are well-established in wide/deep radio clustering studies (e.g. Lindsay et al. Reference Lindsay2014; Magliocchetti et al. Reference Magliocchetti, Popesso, Brusa, Salvato, Laigle, McCracken and Ilbert2017; Hale et al. Reference Hale, Jarvis, Delvecchio, Hatfield and Novak2018) and are reproduced in our Section 5.3. The bias contrast between SFG and AGN is precisely the lever arm exploited by multi-tracer estimators; as such the presence of the trend in our simulation is of the upmost importance.

3.2. Star-forming galaxy radio luminosity function (SFG RLF)

We require an SFG description that is dust-insensitive and straightforward to propagate through survey selections across metre–centimetre bands. Rather than constructing an SFR function from UV/Infrared (IR) tracers and mapping $\mathrm{SFR}\!\to\!L_\nu$ , we calibrate the radio luminosity function of SFG directly at 1.4 GHz and model its smooth redshift evolution. This keeps the baseline independent of assumptions about obscured star formation and the IR-radio correlation, letting the radio data set the space density as a function of $L_\nu$ and z.

Relative to T-RECS, GHOST adopts a radio-derived baseline RLF for the SFG population, in which the 1.4 GHz RLF is fit directly and evolved smoothly with redshift. This reduces reliance on uncertain dust corrections and on the redshift evolution of the IR-radio correlation when setting the space density as a function of $L_\nu$ and z. The consequences of this modelling choice are quantified explicitly in Section 5.10, where an IR-calibrated construction and a thermal-dust SED extension are used as cross-checks. As discussed in Section 5.6, the radio-calibrated baseline tends to produce a smoother n(z) at fixed flux (less pronounced structure around $z\!\sim\!1$ -2) while preserving number counts and large-scale clustering.

For the functional form we adopt a smoothly broken (double-power-law) RLF with redshift-dependent normalisation and knee luminosity. All parameters are obtained by fitting to 1.4 GHz SFG RLF measurements in redshift slices, then validated against differential counts at 150 MHz (Franzen et al. Reference Franzen2016; Mandal et al. Reference Mandal2021), 1.4 GHz (Bondi et al. Reference Bondi, Ciliegi and Schinnerer2008; Bridle et al. Reference Bridle, Davis, Fomalont and Lequeux1972; Ciardi & Loeb Reference Ciardi and Loeb2000; Fomalont et al. Reference Fomalont, Kellermann, Cowie, Capak, Barger, Partridge, Windhorst and Richards2006; Gruppioni et al. Reference Gruppioni, Mignoli and Zamorani1999; Hopkins et al. Reference Hopkins, Afonso, Chan, Cram, Georgakakis and Mobasher2003; Ibar et al. Reference Ibar, Ivison, Biggs, Lal, Best and Green2009; Kellermann et al. Reference Kellermann, Fomalont, Mainieri, Padovani, Rosati, Shaver, Tozzi and Miller2008; Mitchell & Condon Reference Mitchell and Condon1985; Owen & Morrison Reference Owen and Morrison2008; Richards Reference Richards2000; Seymour et al. Reference Seymour2008; White et al. Reference White, Becker, Helfand and Gregg1997; de Zotti et al. Reference de Zotti, Massardi, Negrello and Wall2010) and 3 GHz (Smol $\breve{c}$ ić et al. Reference Smolčić2017) using the SFG spectral prescription applied elsewhere in this work. Fluxes at other frequencies are generated via standard K-corrections with that spectral model. The explicit evolution law and fitted coefficients are listed in Appendix A.

In the multi-tracer context the SFG branch provides the lower-bias tracer with a broadly lower-redshift kernel set by the radio RLF and survey flux limit. Together with the higher-bias AGN branch this delivers the bias contrast and distinct P(z) kernels needed for cosmic-variance cancellation on ultra-large scales. Calibrating the SFG RLF directly in radio reduces sensitivity to dust-related systematics that could otherwise imprint spurious large-scale structure in n(z) and $w(\theta)$ , improving robustness for $f_\mathrm{NL}$ forecasts.

3.2.1. Thermal dust SED term and IR-derived SFG RLF

The fiducial SFG population in GHOST is calibrated directly in radio (Section 3.2). Two additional ingredients are implemented as cross-checks and for multi-frequency completeness: an analytic thermal dust component in the SFG SED, and an IR-calibrated mapping that provides an alternative estimate of the effective $1.4\,\mathrm{GHz}$ space density. For multi-frequency catalogues, the radio SED is the sum of synchrotron and free-free components (Appendix A.6) and a thermal dust term is added,

(1) \begin{equation}{L_\nu^\mathrm{{SFG}}(\nu_\mathrm{rest})=L_\nu^\mathrm{sync}(\nu_\mathrm{rest})+L_\nu^\mathrm{ff}(\nu_\mathrm{rest})+L_\nu^\mathrm{dust}(\nu_\mathrm{rest}) .}\end{equation}

The dust term is modelled as a modified blackbody,

(2) \begin{equation}{L_\nu^\mathrm{dust}(\nu_\mathrm{rest})=A_{d}\,\nu_\mathrm{rest}^{\beta_{d}}\,B_\nu(\nu_\mathrm{rest},T_{d}),}\end{equation}

with fiducial $(T_{d},\beta_{d})=(45\,\mathrm{K},1.5)$ as adopted in the SKADS simulated skies (Wilman et al. Reference Wilman2008). The normalisation $A_{d}$ and the mapping between SFR, $L_\mathrm{{IR}}$ , and $L_{1.4}$ used in the verification are given in Appendix A.7.

3.3. AGN population

For each class we adopt the T-RECS double-power-law RLF with smooth redshift evolution in normalisation and knee luminosity; the functional form is listed in Appendix A. Coefficients are refitted where required to remain consistent with the host assignment on our simulation and with the class mix used in this work; the resulting values and uncertainties are reported in Table 7. Observed flux densities at 150 MHz, 1.4 GHz, and 3 GHz are obtained from $L_\nu$ via standard K-corrections using class-dependent spectral-index distributions (flat for FSRQ/BLLac; steep for SSAGN), as in T-RECS. We do not model variability.

Sizes and a soft edge-brightening proxy are assigned to SSAGN to provide a morphology label used only for host weighting (Section A.4); this proxy is converted to a probabilistic HERG/LERG tag to allow overlap between morphology and excitation classes. Photometry is unaffected by these labels.

Hosts are assigned using the abundance-matching scheme defined in Section 4. In brief, AGN are matched to halos within redshift shells using a monotonic, no-replacement map, with class-dependent host probabilities and stellar-mass-weighted excitation fractions inspired by Janssen et al. (Reference Janssen, Röttgering, Best and Brinchmann2012). Tying the placement to our halo supply stabilises the large-scale bias of each class and ensures that the predicted $w(\theta)$ reflects the simulation’s actual halo population.

Beamed subclasses are handled by drawing jet orientations isotropically and enforcing a maximum viewing angle per class and redshift to reproduce the observed beamed fractions. This yields flat spectra and compact sizes for the beamed subset while leaving the large-scale clustering unchanged.

In summary, we inherit the T-RECS class taxonomy, RLF functional form, and spectral prescriptions, while refitting the RLF coefficients where needed and replacing the original halo mapping with our abundance-matching scheme on the FLAMINGO light-cone.

4.1. Assumptions and limitations

We assign a single central radio source per halo (no explicit satellites). The scatter $\sigma_{\log L}$ is taken to be constant within a shell. HERG/LERG labels are proxies for host-mass weighting rather than strict spectroscopic classes. At the highest redshifts, finite halo supply can saturate the brightest SFG targets; we quantify hosted versus requested counts in Section 5.8.

4.2. Inputs and redshift shells

We partition redshift into disjoint shells: $\mathcal{Z}_i \equiv [z_i, z_{i+1}) \quad (i=0,$ $\dots, N-1),$ so that $\bigcup_{i=0}^{N-1}\mathcal{Z}_i = [z_0,z_N)$ and $\mathcal{Z}_i\cap\mathcal{Z}_j=\varnothing$ for ${i\neq j}$ . (the same grid used for n(z) and RLFs). In each shell we compute the comoving shell volume per steradian and the target numbers per luminosity bin and population from the relevant RLF under the survey selection (cf. Equation 3). Candidate hosts are all halos in the shell, each with mass $M_h$ and light-cone position. We assign a stellar-mass proxy via the redshift-dependent $M_\star(M_h, z)$ relation of Aversa et al. (Reference Aversa, Lapi, de Zotti, Shankar and Danese2015).

4.3. SFG branch: Monotonic L– $M_h$ abundance matching

For SFG we fit a smooth, monotonic mapping so that the host-assigned SFG luminosity histogram reproduces the radio SFG RLF $\Phi_\mathrm{{SFG}}(L_\nu, z)$ (Section 3.2). We adopt

(4) \begin{equation}L_\nu(M_h;\, z)\;=\;N(z)\,\Bigg[\Big(\frac{M_h}{M_\star(z)}\Big)^{\gamma}+\Big(\frac{M_h}{M_\star(z)}\Big)^{\beta}\Bigg]^{-1},\end{equation}

where $M_h$ is the halo mass, and $M_*(z)$ , $\gamma$ , $\beta$ and N(z) are free parameters. with log-normal scatter $\epsilon\sim \mathcal{N}(0,\sigma_{\log L}^2)$ . The parameters $\{N, M_\star, \gamma, \beta, \sigma_{\log L}\}$ are inferred independently in each shell by maximising a Poisson likelihood comparing the target counts (from the RLF) to the counts produced by rank-ordering halos in $M_h$ and mapping to $L_\nu$ via Equation (4). We explore the posterior with an ensemble Markov-Chain Monte-Carlo (MCMC) (emcee, Foreman-Mackey et al. Reference Foreman-Mackey, Hogg, Lang and Goodman2013) and adopt median parameters for assignment.

4.4. AGN branch: Mass-weighted host selection and monotone assignment

For radio AGN we draw hosts with a stellar-mass-dependent probability that encodes HERG/LERG tendencies and then assign luminosities monotonically within the selected subset so that the class RLFs are reproduced.

Each SSAGN is given a stochastic size from DiPompeo et al. (Reference DiPompeo, Runnoe, Myers and Boroson2013) and edge-brightening proxy $R_s$ from Lin et al. (Reference Lin, Shen, Strauss, Richards and Lunnan2010); we label $R_s\gt0.5$ as FR II and $R_s\lt0.5$ as FR I and use this as a soft proxy for excitation class. Following Janssen et al. (Reference Janssen, Röttgering, Best and Brinchmann2012), ; Bonaldi et al. (Reference Bonaldi, Bonato, Galluzzi, Harrison, Massardi, Kay, De Zotti and Brown2019), we weight hosts by stellar mass, with a high-mass cap (‘saturation’) for the LERG branch:

(5) \begin{align}P_\mathrm{LERG}(M_\star|z)\;\propto\;\begin{cases}\left(M_\star/M_0\right)^{2.5}, & M_\star \le M_\mathrm{sat},\\[3pt]\left(M_\mathrm{sat}/M_0\right)^{2.5}, & M_\star \gt M_\mathrm{sat},\end{cases} \\[-25pt] \nonumber \end{align}

(6) \begin{align}P_\mathrm{HERG}(M_\star|z)\;\propto\;\left(\dfrac{M_\star}{M_0}\right)^{1.5},\end{align}

with $M_\mathrm{sat}\simeq 10^{11.6}$ M $_\odot$ . By ‘saturating’ we mean the LERG triggering weight stops increasing above $M_\mathrm{sat}$ . We renormalise $P_t(M_\star|z)$ over candidate hosts in each redshift shell and draw hosts so that HERG/LERG totals match the RLF targets for that shell (FSRQ $\rightarrow$ HERG-like; BLLac $\rightarrow$ LERG-like; SSAGN split by $R_s$ ). Hosts are sampled without replacement. Within each class and shell we sort hosts by $M_h$ , sort luminosities by $L_\nu$ , and map by rank (with optional log-normal L scatter) to reproduce the class RLF while anchoring large-scale bias to the actual halo field. FSRQ/BL Lac orientations are drawn isotropically (uniform in $\cos\theta$ ).

Table 2.

Classification schemes of comparison surveys and the emulation used in this work.

4.5. Population-weighted host apportionment

Let $\mathcal{L}_t(s)$ be the set of luminosities requested in slice s for population $t\in\{\mathrm{FSRQ},\mathrm{BL\,Lac},\mathrm{SSAGN},\mathrm{SFG}\}$ , with size $R_t(s)=|\mathcal{L}_t(s)|$ . Let K(s) be the number of available hosts (unassigned halos) in the slice. If $K(s)\ge R_\mathrm{tot}(s)\equiv\sum_t R_t(s)$ , all requests are hosted.

When $K(s) \lt R_\mathrm{tot}(s)$ , we apportion the K(s) hosts across populations in proportion to requests using a largest-remainder rule:

(7) \begin{align}q_t \;=\; \frac{R_t(s)}{R_\mathrm{tot}(s)}\,K(s), \qquad K_t \;=\; \min\bigl(R_t(s),\,\lfloor q_t\rfloor \bigr),\end{align}

then distribute the remaining $K(s)-\sum_t K_t$ one by one to populations with the largest fractional remainders $q_t-\lfloor q_t\rfloor$ (never exceeding $R_t$ ). This guarantees $\sum_t K_t=\min (K(s),R_\mathrm{tot}(s))$ and $K_t\le R_t$ for all t.

Within each population, sort luminosities descending, $L_{t,(1)}\ge\cdots\ge L_{t,(R_t)}$ , and map rank to a host-mass index inside the quota via

(8) \begin{equation}i(r) \;=\; \Bigl\lceil r\,\frac{K_t}{R_t}\Bigr\rceil \in \{1,\dots,K_t\}.\end{equation}

Across all populations we collect these (predicted-mass, population) pairs, sort by predicted mass, and assign the top K(s) to the actual top K(s) halos (strictly preserving mass rank). Leftover requests $R_t-K_t$ are written as ‘no-halo’ rows (*NH) at the slice mid-redshift, preserving per-population totals and making host scarcity explicit.

4.6. Algorithmic summary

In each shell we compute target counts per class and luminosity bin from the relevant RLF (Equation 3) under the survey selection, gather all halos and assign $M_\star(M_h, z)$ from Aversa et al. (Reference Aversa, Lapi, de Zotti, Shankar and Danese2015), and then perform assignment by branch. For SFG we fit $L_\nu(M_h)$ via Equation (4) using an ensemble MCMC on a Poisson likelihood and assign by rank without replacement, adding log-normal scatter. For AGN we draw HERG/LERG hosts using $P_\mathrm{{HERG}/{LERG}}(M_\star)$ and, within each class, assign $L_\nu$ monotonically to match the class RLF. We then attach sky positions from hosts, apply spectral prescriptions and K-corrections, and finally apply the survey selection (flux/SNR survey selection functions and completeness surrogates).

5.1. Selections used for validation

This Section defines the cuts and clarifies how each survey separated AGN from SFG (where they did), and how we mirrored those choices in GHOST. The three emulations are: the VLA-COSMOS 1.4 GHz luminosity threshold used by Magliocchetti et al. (Reference Magliocchetti, Popesso, Brusa, Salvato, Laigle, McCracken and Ilbert2017); the second data release (DR2) of LoTSS (Hale et al. Reference Hale2024; Shimwell et al. Reference Shimwell2022) wide-area 150 MHz flux and peak–S/N survey selection functions; and the VLA–COSMOS 3 GHz peak–S/N threshold and multi-wavelength classes of Hale et al. (Reference Hale, Jarvis, Delvecchio, Hatfield and Novak2018). Table 2 summarises the survey selection functions and our emulation strategy. As emphasised in Section 5.3, these selections compress or broaden the redshift kernel n(z) and thus reshape the projected clustering amplitude/slope without altering the intrinsic bias ordering, so we present raw and selected catalogues side by side in what follows.

Recent deep 1.4 GHz surveys with MeerKAT and the VLA (e.g. MeerKAT International GHz Tiered Extragalactic Exploration Survey (MIGHTEE); Jarvis et al. Reference Jarvis2016, GOODS-N; Owen Reference Owen2018) provide valuable constraints on the faint-flux population mix and redshift distributions. We do not adopt them as primary validation targets in this paper because our emphasis is on clustering with published two-point measurements. Nevertheless, as a consistency check in the $\unicode{x03BC}$ Jy regime, we also apply a MIGHTEE-like 1.4 GHz flux cut and inspect the resulting redshift distributions and class fractions in GHOST, finding the expected shift towards an SFG-rich sample while retaining a non-negligible AGN contribution at moderate redshift (e.g. Whittam et al. Reference Whittam2024).

5.1.1. LoTSS–DR2 wide 150 MHz

We emulate the LoTSS selection of (Shimwell et al. Reference Shimwell2022; Hale et al. Reference Hale2024) by combining an integrated-flux cut $S_{150}\ge 1.5\,\mathrm{mJy}$ with a peak–S/N threshold $\mathrm{ S/N}_\mathrm{peak}\ge 7.5$ . This corresponds to a SFR of $\approx$ $6.7\times10^2$ M $_\odot$ /yr. The image rms varies across the mosaic due to multiplicative effects, such as primary beam sensitivity, mosaicking overlap, and calibration/dynamic range, so we model the local rms $\sigma$ as log-normal in an effort to replicate this, $\ln\sigma\sim\mathcal{N}(\mu_{\ln},\sigma_{\ln}^2)$ . Anchoring the distribution to the median $\tilde{\sigma}=83\,\unicode{x03BC}\mathrm{Jy\,beam^{-1}}$ and 95th percentile $171\,\unicode{x03BC}\mathrm{Jy\,beam^{-1}}$ gives

(9) \begin{align}\mu_{\ln}=\ln\tilde{\sigma},\qquad\sigma_{\ln}=\frac{\ln\!(171 \unicode{x03BC}\mathrm{Jy})-\mu_{\ln}}{1.645}, \end{align}

which reproduces the observed quantiles by construction. For a source of integrated flux S (we take $S_\mathrm{peak}\!\approx\!S$ at 1.5 mJy; for resolved sources this is conservative), the S/N cut is equivalent to $\sigma\le\sigma_\mathrm{thr}=S/7.5$ , yielding a detection probability

(10) \begin{equation}p_\mathrm{keep}(S)=\Phi\!\left(\frac{\ln\sigma_\mathrm{thr}-\mu_{\ln}}{\sigma_{\ln}}\right),\end{equation}

with $\Phi$ the standard normal CDF. Using the class counts $dN_X/dS$ , we define the area-weighted acceptances

(11) \begin{equation}\mathrm{acc}_X=\frac{\displaystyle \int_{S_{\min}}^\infty \frac{dN_X}{dS}(S)\,p_\mathrm{keep}(S)\,dS} {\displaystyle \int_{S_{\min}}^\infty \frac{dN_X}{dS}(S)\,dS},\qquad X\in\{\mathrm{SFG},\mathrm{AGN}\},\end{equation}

and form the effective mixture entering the total n(z),

(12) \begin{equation}w_\mathrm{{SFG}}=\frac{\mathrm{acc}_\mathrm{{SFG}}\,N_\mathrm{{SFG}}}{\mathrm{acc}_\mathrm{{SFG}}\,N_\mathrm{{SFG}}+\mathrm{acc}_\mathrm{AGN}\,N_\mathrm{AGN}},\end{equation}

where $S_{\min}=1.5\,\mathrm{mJy}$ and $N_X$ are class totals from the un-gated n(z). At $S_{150}=1.5\,\mathrm{mJy}$ the flux threshold removes only a few percent of sources in both classes, so the resulting n(z) is dominated by the intrinsic SFG/AGN mix (Figure 9). Hale et al. (Reference Hale2024) does not supply a homogeneous AGN/SFG split for this clustering sample, so we emulate the selection only and compare the total n(z) and clustering; by-class predictions shown elsewhere use our internal labels to isolate population-mix effects from detection selection.

5.1.2. VLA–COSMOS 1.4 GHz

We mirror the flux density limit $S_\mathrm{1.4\,GHz}\ge 0.15$ mJy of Magliocchetti et al. (Reference Magliocchetti, Popesso, Brusa, Salvato, Laigle, McCracken and Ilbert2017) and re-classify sources using their redshift-dependent radio-luminosity boundary. This corresponds to a SFR of $\approx$ $5.3\times10^2$ M $_\odot$ /yr. Radio powers are computed as $P_{1.4}=F_{1.4}\,D_A^2\,(1+z)^{3+\alpha}$ (W Hz $^{-1}$ sr $^{-1}$ ) with $\alpha=0.7$ and $D_A$ the angular diameter distance, and the AGN/SFG division is

(13) \begin{align}\log_{10} P_\mathrm{cross}(z)=\begin{cases}21.7 + z, & z\le 1.8 \nonumber \\[3pt]23.5, & z \gt 1.8\end{cases}\end{align}

(we convert to isotropic luminosity via $L_\mathrm{cross}=4\pi P_\mathrm{cross}$ ). Objects with $L_{1.4}\ge L_\mathrm{cross}(z)$ are labelled ‘AGN-like’, the rest ‘SFG-like’. This reproduces their threshold classifier exactly and makes the comparison sensitive to population mix and redshift tails rather than to differences in detection survey selection functions; it also emphasises the overlapping mJy regime where both classes contribute (see Figure 10).

5.1.3. VLA–COSMOS 3 GHz

We emulate (Hale et al. Reference Hale, Jarvis, Delvecchio, Hatfield and Novak2018) catalogue’s peak-flux threshold $\mathrm{S/N}\ge 5.5$ with a field median rms $\tilde{\sigma}\simeq 2.3\,\unicode{x03BC}\mathrm{Jy\,beam^{-1}}$ , modelling field-to-field rms as lognormal with

(14) \begin{equation}\mu_{\ln}=\ln\tilde{\sigma},\qquad \sigma_{\ln}=0.30.\end{equation}

This corresponds to a SFR of $\approx 3.2\times10^2$ M $_\odot$ /yr. For a source of flux density S, the rms required to pass is $\sigma_\mathrm{thr}=S/5.5$ ; the area fraction that admits detection is

(15) \begin{equation}p_\mathrm{keep}(S)=\Phi\!\left(\frac{\ln\!(S/5.5)-\mu_{\ln}}{\sigma_{\ln}}\right).\end{equation}

We then compute a class-wise, area-weighted acceptance by averaging $p_\mathrm{keep}$ over the per-population 3 GHz counts histogram $H_X(S_i)$ (from the simulation),

(16) \begin{equation}\mathrm{acc}_X=\frac{\sum_i H_X(S_i)\,p_\mathrm{keep}(S_i)}{\sum_i H_X(S_i)},\qquad X\in\{\mathrm{SFG},\mathrm{AGN}\},\end{equation}

and form the effective mixture

(17) \begin{equation}w_\mathrm{{SFG}}=\frac{\mathrm{acc}_\mathrm{{SFG}}\,N_\mathrm{{SFG}}}{\mathrm{acc}_\mathrm{{SFG}}\,N_\mathrm{{SFG}}+\mathrm{acc}_\mathrm{AGN}\,N_\mathrm{AGN}}.\end{equation}

At the Hale depth the selection function is weak; in our run the acceptances are close to unity for both classes (numerically $\mathrm{acc}_\mathrm{{SFG}}\!\approx\!0.89$ , $\mathrm{ acc}_\mathrm{AGN}\!\approx\!0.94$ ), so the predicted n(z) is governed mainly by the intrinsic SFG/AGN mix. Because the Hale classification relies on deep multi-wavelength diagnostics (X-ray, infrared colour, UV–IR SED, rest-frame NUV colour, and radio excess relative to SFR), we do not attempt to reproduce their labels inside GHOST; we emulate only the S/N selection with a final flux cut of $12.5\,\unicode{x03BC}$ Jy at $3\,$ GHz and compare the total n(z), using our internal classes where a by-class illustration is helpful.

5.2. Abundance-matching outcomes

All results that follow use the host-assigned catalogue from Section 4: every mock source is either attached to a specific FLAMINGO halo or, if the slice’s host supply is exhausted, written as a ‘no-halo’ (*NH) row. When halos are scarce we apportion the available hosts across populations in proportion to their requested counts (largest-remainder rule; Section 4.5) and then assign luminosities monotonically within each population so that the class RLFs are reproduced.

At fixed flux limit the accessible radio luminosity rises with redshift,

(18) \begin{equation}L_\mathrm{lim}(z,\zeta)=4\pi\,D_L^2(z)\,S_\mathrm{lim}\,(1+z)^{\alpha-1},\end{equation}

so the observed $(L_\nu,z)$ domain is bounded jointly by Equation (18) and the evolving halo mass function. With the corrected slice projection, both the abundance-matching targets and the host-assigned catalogue populate this domain smoothly across all slices.

Figure 2 summarises the delivered mapping. SFGs occupy lower-mass halos and lower $L_\nu$ ; steep-spectrum AGN extend to higher $M_h$ with a steeper L–M trend; beamed subclasses (FSRQ, BLLac) preferentially inhabit the most massive halos and dominate the bright tail. With increasing redshift, median ridge-lines shift as expected from K-corrections and the evolving halo supply, and the per-population marginals narrow modestly.

Figure 2.

All populations – radio luminosity vs. halo mass with per-population marginals (shared density scale). Each panel shows 2D binned abundances in a fixed redshift window; coloured semi-transparent maps overlay populations (SFG solid, SSAGN dashed, FSRQ dash-dot, BLLac dotted). All density layers share the same logarithmic colourbar, so intensity is directly comparable across populations and panels. Solid lines trace the median of $\log_{10}L_\nu$ at fixed $\log_{10}M_h$ . The small top (right) axes show, for each population, the normalised marginal PDFs $p(\log_{10}M_h)$ and $p(\log_{10}L_\nu)$ – with $\int p\,{d}\log M_h=1$ and $\int p\,{d}\log L=1$ . SFG peak at lower host masses and luminosities; SSAGN extend to higher $M_h$ with a steeper $L_\nu$ – $M_h$ trend; FSRQ and BL Lac preferentially inhabit the most massive halos and dominate the bright radio tail. With increasing redshift, the medians shift to higher $M_h$ and the PDFs narrow, reflecting both selection effects and the evolving halo-occupation mix.

We track requested versus hosted counts per slice and population (Section 4.5). The hosted fraction is high across all panels; occasional *NH rows appear only at very low z and $L_\nu$ and are negligible for the selections used. In the RLF grids we show the model projected on the same (L, z) grid alongside the host-assigned estimate; the two agree within uncertainties over the full dynamic range.

The host-assigned catalogue embeds a clear bias hierarchy – SFG $\lt$ SSAGN $\lt$ FSRQ/BL Lac – set by the actual halo field under the survey selection. Because our population-weighted apportionment keeps the hosted mixture close to the requested mixture for any flux/SNR cut, changes in $w(\theta)$ across selections reflect projection kernels rather than a reversal of intrinsic bias ordering. This alignment of mixture, bias hierarchy, and selection is exactly what multi-tracer $f_\mathrm{NL}$ estimators require.

5.3. The angular correlation function

Our goal in this section is to verify that the catalogue places each population in environments consistent with their expected large-scale bias, and to produce clustering summaries that feed the PNG analysis in Paper 2. We therefore measure the two-point angular correlation of the simulated sky, both in configuration space, $w(\theta)$ , and in harmonic space, $C_\ell$ . Given a redshift kernel $\phi(z)$ (from the n(z) of each population), these observables probe the projection of the 3D matter power spectrum modulated by the tracers’ bias, so their amplitudes and slopes are sensitive to the typical host halo masses of SFGs versus AGN subclasses and how survey selections reshape $\phi(z)$ .

Concretely, we compute $w(\theta)$ for AGN and SFG samples (and for specific survey-like cuts), fit a compact power-law model $w(\theta)=A(\theta/\theta_0)^{1-\gamma}$ over $0.02^\circ$ – $1^\circ$ to quote comparable amplitudes and slopes, and cross-check in harmonic space by fitting $C_\ell$ with fixed $\phi(z)$ templates to infer effective linear bias (both constant and gently evolving with redshift).

5.3.1. Definitions and modelling

We quantify angular clustering with the two-point correlation function $w(\theta)$ , related to the angular power spectrum $C_\ell$ by

(19) \begin{equation}w(\theta)\;=\;\sum_{\ell=0}^{\infty}\frac{2\ell+1}{4\pi}\,C_\ell\,P_\ell(\cos\theta)\;\simeq\;\int\frac{\ell\,{d}\ell}{2\pi}\,J_0(\ell\theta)\,C_\ell,\end{equation}

where the Bessel form is accurate for small angles. In the Limber regime, the auto-spectrum for a population $X\in\{\mathrm{AGN},\mathrm{SFG}\}$ is

(20) \begin{align} C_\ell^{XX}&=\int {d}z\;\underbrace{\frac{H(z)}{\chi^2(z)}\,\phi_X^2(z)}_{\text{projection}}\;\underbrace{b_X^2(k,z)\,D^2(z)}_{\text{bias/growth}}\;\nonumber\\[3pt] &\quad \times P_m\!\left(k=\frac{\ell+1/2}{\chi(z)},z\right),\end{align}

with comoving distance $\chi(z)$ , Hubble rate H(z), linear growth D(z), matter spectrum $P_m$ , and a normalised redshift kernel $\phi_X(z)$ (unit area).

Across a finite angular decade the observed slope of $w(\theta)$ reflects a weighted average of the 3D logarithmic slope $ {d}\ln[P_m\,b_X^2]/{d}\ln k $ evaluated at $k\simeq \ell/\chi(z)$ , convolved with $\phi_X^2$ . For direct comparison to the literature and to quote compact summary numbers we fit a power law

(21) \begin{equation}w(\theta)=A\left(\frac{\theta}{\theta_0}\right)^{1-\gamma},\end{equation}

using full covariances where available and a pivot $\theta_0$ taken to be the geometric mean of the fitted $\theta$ to minimise the A– $\gamma$ covariance. Unless stated otherwise we fit over $\theta\in[0.02^\circ,\,1.0^\circ]$ ( $\ell\simeq 50$ –400).

To interpret slope changes induced by selections, we also quote effective redshifts and distances computed with Limber-like weights

(22) \begin{align} W_X(z)\;&\propto\;\frac{H(z)}{\chi^2(z)}\,\phi_X^2(z)\,D^2(z),\qquad \\[-10pt] \nonumber \end{align}

(23) \begin{align} z_{\mathrm{eff}, X}&=\frac{\int \! W_X(z)\,z\,{d}z}{\int \! W_X(z)\,{d}z},\qquad \\[-10pt] \nonumber \end{align}

(24) \begin{align} \chi_{\mathrm{eff}, X}&=\frac{\int \! W_X(z)\,\chi(z)\,{d}z}{\int \! W_X(z)\,{d}z}. \\[10pt] \nonumber \end{align}

A fixed angular separation then corresponds to a characteristic projected comoving scale

(25) \begin{equation}r_\perp(\theta)=\theta\,\chi_{\mathrm{eff}, X}.\end{equation}

5.3.2. Raw GHOST sample

In Figure 3, we present the angular correlation function of the full GHOST sample showing AGN lying above SFG in amplitude across all $\theta$ , consistent with their higher bias and host-halo masses. Power-law fits give $\gamma_\mathrm{AGN}=1.886\pm0.028$ and $\gamma_\mathrm{{SFG}}=1.647\pm0.076$ . The effective depths are $(z_\mathrm{eff},\chi_\mathrm{eff})=(0.318,\,1.22 \mathrm{Gpc})$ for AGN and $(0.602,\,2.12\,\mathrm{Gpc})$ for SFG, so $r_\perp(0.05^\circ)=(1.06,\,1.85) \mathrm{Mpc}$ , respectively. The shallower SFG slope in the raw GHOST sample is naturally explained by their broader, higher-z kernel $\phi(z)$ , which mixes a wider range of $k=\ell/\chi$ in projection (Equations 20–24).

Figure 3.

Raw GHOST sample: Angular correlation function $w(\theta)$ for AGN (blue triangles) and SFG (orange circles). Solid lines show inverse-variance weighted power-law fits $w(\theta)=A(\theta/\theta_0)^{1-\gamma}$ over $\theta\in[0.02^\circ,\,1.0^\circ]$ ; shaded bands are the propagated $1\sigma$ model uncertainties from the $(A,\gamma)$ covariance. Best-fit slopes are $\gamma_\mathrm{AGN}=1.886\pm0.028$ and $\gamma_\mathrm{{SFG}}=1.647\pm0.076$ . The lower panel shows data $-$ fit residuals. The amplitude ordering AGN $\gt$ SFG reflects the intrinsic bias hierarchy, while the slightly shallower SFG slope is consistent with a broader, higher-z kernel $\phi(z)$ that mixes more $k=\ell/\chi(z)$ in projection (Equations 19 and 20).

5.3.3. VLA–COSMOS 1.4 GHz

Applying the $L_{1.4}\gtrless L_\mathrm{cross}(z)$ selection, as detailed in Section 5.1.2, compresses the SFG kernel to very low z. We find $\gamma_\mathrm{AGN}=1.858\pm0.041$ and $\gamma_\mathrm{{SFG}}=1.667\pm0.010$ , with $(z_\mathrm{eff},\chi_\mathrm{eff})=(0.339,\,1.30\,\mathrm{Gpc})$ for AGN and $(0.094,\,0.41\,\mathrm{Gpc})$ for SFG, implying $r_\perp(0.05^\circ)=(1.14,\,0.35)\,\mathrm{Mpc}$ . In Figure 4, the SFG correlation lies above the AGN over the fitted range. This does not imply a larger intrinsic SFG bias; it follows from the projection weights in the Limber integral, which scale as $\phi^2(z)/\chi^2(z)$ . Our selection compresses the SFG kernel to very low redshift ( $z_\mathrm{eff}\simeq0.094$ , $\chi_\mathrm{eff}\simeq0.41\,\mathrm{Gpc}$ ) while the AGN kernel sits higher ( $z_\mathrm{eff}\simeq0.339$ , $\chi_\mathrm{eff}\simeq1.30\,\mathrm{Gpc}$ ). The resulting geometric factor $(\chi_\mathrm{AGN}/\chi_\mathrm{{SFG}})^2\simeq(1.30/0.41)^2\!\approx\!10$ boosts the projected SFG amplitude and can outweigh the larger AGN bias, thereby producing the observed inversion in $w(\theta)$ . The residual shape mismatch – most evident in the flatter SFG slope – likely reflects blending of one- and two-halo terms when projecting such a low-z kernel over $0.02^\circ\!-\!1^\circ$ ( $r_\perp\!\sim\!0.14\!-\!7\,\mathrm{Mpc}$ ), and modelling choices that suppress small-scale power (e.g. our current radio-SFG occupation lacks an explicit satellite component). These effects account for the difference in gradient while preserving the relative bias hierarchy that underpins the multi-tracer analysis.

Figure 4.

VLA-COSMOS selection (1.4 GHz): $w(\theta)$ for AGN and SFG after the $L_{1.4}\gtrless L_\mathrm{cross}(z)$ selection function detailed in Section 5.1.2. Solid curves are our power-law fits over $\theta\in[0.02^\circ,\,1.0^\circ]$ (shaded $1\sigma$ bands); dotted lines and translucent bands show the literature reference for AGN and SFG, respectively. We obtain $\gamma_\mathrm{AGN}=1.858\pm0.041$ and $\gamma_\mathrm{{SFG}}=1.667\pm0.010$ . The SFG amplitude lies above the reference at large scales and the slope is slightly flatter. From the effective-depth summary (Table 3) the SFG kernel compresses to $z_\mathrm{eff}\simeq0.094$ and $\chi_\mathrm{eff}\simeq0.41$ Gpc, so a fixed $\theta$ probes small $r_\perp$ (e.g. $r_\perp(0.05^\circ)\approx0.35$ Mpc), flattening the angular scaling (Equations 24 and 25).

5.3.5. Picture across frequencies

The differences we see across the raw, Magliocchetti and Hale-like samples are driven primarily by how each selection reshapes the redshift kernel $\phi(z)$ , not by any change in the intrinsic bias ordering. In Equation (20) the amplitude scales as $b^2\,\phi^2/\chi^2$ and the slope reflects which 3D wavenumbers are projected ( $k\simeq \ell/\chi$ ) given the kernel. Cuts that compress the SFG kernel to lower redshift (e.g. luminosity splits or high S/N survey selection functions) increase $\phi^2/\chi^2$ and shift the map towards larger k, yielding the observed ‘lift and flattening’ of SFG $w(\theta)$ . AGN kernels move much less under the same cuts, so the intrinsic hierarchy $b_\mathrm{AGN} \gt b_\mathrm{{SFG}}$ remains intact.

This explains the apparent discrepancy between the Magliocchetti and Hale measurements: once their selection functions are emulated, both behaviours arise naturally from projection effects with a fixed bias ordering. The cross-frequency consistency (1.4 vs 3 GHz) further shows the effect is selection-driven rather than band-dependent. For multi-tracer $f_\mathrm{NL}$ applications, the lesson is operational: clustering inferences must be made with tracer-specific $\phi(z)$ that match the actual survey survey selection functions; otherwise slope and amplitude shifts can be misread as bias changes. Here we carry forward bias parameters together with the matched kernels to ensure the large-scale modes relevant for $f_\mathrm{NL}$ are interpreted consistently. The statistics of the ACF are summarised in Table 3.

Table 3.

Projection summary for AGN and SFG by selection. Errors are $1\sigma$ from power-law fits over $\theta\in[0.02^\circ,0.5^\circ]$ . Effective depths use the Limber-like weights of Equation (24). We quote the characteristic projected comoving separation $r_\perp(\theta)=\theta\,\chi_\mathrm{eff}$ at three angles.

Table 4.

Bias fits from $C_\ell$ over $\ell\in[40\,,\,250]$ . We fit a constant-bias model ( $b=$ const) and an evolving form $b(z)=b_0(1+z)^\zeta$ . Quoted $\chi^2$ values are for the best fit of each model. The last column gives the constant-bias ratio $b/b_\mathrm{{SFG}}$ .

5.4. Angular power spectra and bias fits

We model the auto-spectra of each radio population $X\in\{\mathrm{SFG},\mathrm{SSAGN},\mathrm{BLLac},\mathrm{FSRQ}\}$ in the Limber regime using Equation (20). We construct, for each tracer, a template $T_\ell^X$ by evaluating Equation (20) with a unit, redshift-independent bias (i.e. $b(z) \equiv 1$ ). The measured spectrum in a given $\ell$ -band then determines an overall amplitude related to the bias.

We fit two bias models over $\ell\in[\ell_{\min},\ell_{\max}]$ (default 40–250), after sanitising n(z) (Section 5.6) and optionally subtracting a shot component from the auto-correlations when available.

5.4.1. Constant bias

We assume $b_X(z) = b_X$ and fit an amplitude $a_X = b_X^2$ against the fixed template $T_\ell^X$ . Using (optionally) Knox weights $w_\ell=[2/((2\ell+1)f_\mathrm{sky})]/\sigma^2_\ell$ or per-band errors we solve the weighted least-squares problem

(26) \begin{equation}\widehat{a}_X \;=\; \frac{\sum_\ell w_\ell\,T_\ell^X\,C_\ell^{\mathrm{obs}}}{\sum_\ell w_\ell\,(T_\ell^X)^2},\qquad\widehat{b}_X\;=\;\sqrt{\widehat{a}_X}\,,\end{equation}

with error propagation $\sigma^2_{a}\simeq 1/\!\sum_\ell w_\ell(T_\ell^X)^2$ and $\sigma_b\simeq \sigma_a/(2\,\widehat{b})$ . Goodness of fit is summarised by $\chi^2=\sum_\ell w_\ell(C_\ell^\mathrm{obs}-\widehat{a}\,T_\ell)^2$ .

Figure 5.

VLA-COSMOS selection (3 GHz): Same as Figure 4 but for the 3 GHz cut (Section 5.1.3. Our fits give $\gamma_\mathrm{AGN}=1.903\pm0.037$ and $\gamma_\mathrm{{SFG}}=1.662\pm0.010$ . AGN closely follow the band from Hale et al. (Reference Hale, Jarvis, Delvecchio, Hatfield and Novak2018), while SFG sit high and slightly flatter. The projection summary (Table 3) shows $z_\mathrm{eff,\,SFG}\simeq0.161$ and $\chi_\mathrm{eff,\,SFG}\simeq0.67$ Gpc (so $r_\perp(0.05^\circ)\approx0.59$ Mpc), again explaining the slope change as a selection-driven projection effect rather than frequency-dependent systematics. Residuals (lower panel) are shown relative to the fitted power laws.

Figure 6.

$C_\ell$ bias fits by population. Points show measured bandpowers; dashed lines are best-fit constant-bias models ( $b=$ const) and solid lines are evolving models $b(z)=b_0(1+z)^\zeta$ , both fit over $\ell\in[40,250]$ . Numerical values are listed in Table 4. The amplitude ordering SFG $\lt$ SSAGN $\lesssim$ FSRQ $\simeq$ BLLac and the positive $\zeta$ for AGN are consistent with the halo–mass distributions and the configuration-space results (Section 5.3).

Figure 6. Long description

Panel A: A line graph titled SFG. The x-axis is labeled with the variable l and ranges from 40 to 200. The y-axis is labeled G1 and ranges from 0 to 10^-6. The graph includes two lines: a dashed orange line representing a constant-bias model with const b=1.55 and a solid green line representing an evolving model with evol b0=1.87, alpha=1.09. Data points are marked with blue symbols. Panel B: A line graph titled SSAGN. The x-axis is labeled with the variable l and ranges from 40 to 200. The y-axis is labeled G1 and ranges from 0 to 10^-6. The graph includes two lines: a dashed orange line representing a constant-bias model with const b=1.95 and a solid green line representing an evolving model with evol b0=3.05, alpha=3.41. Data points are marked with blue symbols. Panel C: A line graph titled BLLac. The x-axis is labeled with the variable l and ranges from 40 to 200. The y-axis is labeled G1 and ranges from 0 to 10^-6. The graph includes two lines: a dashed orange line representing a constant-bias model with const b=2.32 and a solid green line representing an evolving model with evol b0=4.35, alpha=4.19. Data points are marked with blue symbols. Panel D: A line graph titled FSRQ. The x-axis is labeled with the variable l and ranges from 40 to 200. The y-axis is labeled G1 and ranges from 0 to 10^-6. The graph includes two lines: a dashed orange line representing a constant-bias model with const b=2.27 and a solid green line representing an evolving model with evol b0=4.64, alpha=4.71. Data points are marked with blue symbols.

5.4.2. Evolving bias

We test a minimal two-parameter form

(27) \begin{equation}b_X(z)=b_{0,X}\,(1+z)^{\zeta_X}.\end{equation}

Within our pipeline this redshift dependence is folded into the shape of the projection kernel by replacing $\phi_X(z)\!\to\!\phi_X(z)\,{(1+z)^{\zeta_X}}$ (and re-normalising to unit area), while keeping the overall amplitude as $b_0^2$ . For each tracer we form the profile likelihood in $\zeta$ by optimising $b_0$ at fixed $\zeta$ using Equation (26), that is, $\widehat b_0(\zeta)=\arg\min_{b_0}\chi^2(b_0,\zeta)$ . We then select $(\widehat b_0(\widehat\zeta),\widehat\zeta)$ with $\widehat\zeta=\arg\min_\zeta \chi^2(\widehat b_0(\zeta),\zeta)$ .

5.4.4. Implications

The bias ordering mirrors the halo–mass distributions inferred from the catalogue: SFG peak at lower masses ( $\sim$ $10^{12}$ M $_\odot$ ) while AGN subclasses peak at higher masses (few $\times 10^{12.5-13}$ M $_\odot$ ), naturally producing larger large-scale clustering amplitudes. The mild preference for positive $\zeta$ in AGN ( $\zeta\sim3$ –5) indicates bias growth with redshift, aligning with the picture of AGN inhabiting rarer halos at earlier times. Overall, the $C_\ell$ fits and the configuration-space results (Section 5.3) are mutually consistent: amplitudes trace $b^2$ while the $w(\theta)$ slope variations are explained by projection effects rather than changes in the underlying bias hierarchy.

Figure 7.

Host halo mass PDFs by population. normalised distributions of $\log_{10}(M_h/M_\odot)$ (per dex) for SFG, SSAGN, BLLac, and FSRQ. Distributions for the AGN subclasses peak near $10^{12.8-13}$ M $_\odot$ , while SFG peak near $10^{12}$ M $_\odot$ . The mass ordering maps onto the bias hierarchy measured in Section 5.4.

5.5. Host halo mass distributions

To connect clustering amplitudes with halo occupation, we examine the distribution of host halo masses $M_h$ for each population. Figure 7 shows per-population PDFs of $\log_{10}(M_h/M_\odot)$ (normalised per dex). Table 5 lists the median, the central 68% interval, and the sample size.

The ordering is clear: SFG inhabit lower-mass halos (median $\log_{10}M_h \simeq 11.97$ ), while all AGN subclasses peak near $\sim$ $10^{12.6-13.0}$ M $_\odot$ (SSAGN $\simeq12.65$ , BLLac $\simeq12.95$ , FSRQ $\simeq12.92$ ). The interquartile widths (IQR) are similar ( $\sim$ 0.63–0.68 dex), indicating that the dominant difference is a shift in typical mass rather than a change in spread.

Because the large-scale bias b(M, z) is a monotonic function of halo mass (e.g. in peak-background split models),

(28) \begin{align}\begin{split}b(M,z)\; &=\;1+\frac{a\nu^2-1}{\delta_c}+\frac{2p}{\delta_c\,[1+(a\nu^2)^p]}\,,\qquad\\[2pt]\nu &\equiv \frac{\delta_c}{\sigma(M,z)} ,\end{split}\end{align}

(with $\delta_c\simeq1.686$ and $(a,p)\approx(0.707,0.3)$ as a common choice), the $\sim$ 0.7–1.0 dex upward shift from SFG to radio-loud AGN naturally implies larger b at the redshifts that dominate our projection. This is precisely what we measure from $C_\ell$ in Section 5.4: the constant-bias fits obey $b_\mathrm{{SFG}} \lt b_\mathrm{{SSAGN}} \lesssim b_\mathrm{{FSRQ}} \simeq b_\mathrm{BLLac}$ , and the evolving model prefers positive $\zeta$ for AGN, consistent with rarer halos at earlier times.

Table 5.

Summary of host halo mass distributions per population. We report the sample size N, the median $\mathrm{med}[\log_{10}M_h/M_\odot]$ , the central 68% interval (16th–84th percentiles), the IQR width, and the median offset relative to SFG.

Overall, the mass PDFs corroborate the bias/clustering hierarchy found in both configuration and harmonic space: AGN populate more massive halos than SFG, hence exhibit higher large-scale clustering amplitudes, while the small changes in $w(\theta)$ slope across selections are driven by projection (Section 5.3) rather than a reversal of the intrinsic bias-mass ordering.

5.6. Redshift distributions n(z)

Figure 8 shows normalised redshift distributions for the full GHOST catalogue and for AGN and SFG sub-samples (area under each curve equals unity). Before any survey selection, the total population peaks near $z\!\sim\!0.8$ with a long tail to $z\gt3$ (median $z=1.21$ , IQR $0.74$ - $1.81$ ). The AGN distribution is narrower (median $z=0.85$ , IQR $0.50$ – $1.31$ ), while SFG peak at slightly higher redshift (median $z=1.23$ , IQR $0.75$ – $1.83$ ) and maintain the broader high-z tail. This is consistent with radio-loud AGN dominating the bright counts at low z and SFG contributing increasingly towards the sub-mJy regime.

At our pre-selection stage with the raw GHOST catalogue, the n(z) for SFG and AGN look consistent with the behaviour from deep $\unicode{x03BC}$ Jy radio surveys: SFG dominate at faint fluxes and exhibit a broad distribution that extends to $z\approx4$ , while AGN show a more concentrated, lower z peak. Similar behaviour is seen in the acsLoTSS Deep fields (Best et al. Reference Best2023; Bonzini et al. Reference Bonzini, Padovani, Mainieri, Kellermann, Miller, Rosati, Tozzi and Vattakunnel2013; Smol $\breve{c}$ ić et al. Reference Smolčić2017). Disagreements that emerge once survey cuts are applied should therefore be read chiefly as mixture/evolution issues rather than a failure of the overall redshift density.

Figure 8.

Normalised PDFs P(z) (unit area) for the full GHOST catalogue (top) and by class (AGN, SFG; middle/bottom). Vertical lines mark medians (solid) and IQRs (dashed). The total peaks near $z\!\sim\!0.8$ with a long high-z tail, establishing the baseline cosmic distribution before survey selection.

Our SFG p(z) shape – a broad peak around $z\!\sim\!1$ with a long, shallow high-z tail – is qualitatively consistent with the component in empirical SKAO telescope forecasting work (e.g. Mancuso et al. Reference Mancuso2015, their Figure 10). Two differences help explain why our curves are smoother and show no pronounced ‘bump’ near $z\!\sim\!1$ at $\sim\!\unicode{x03BC}$ Jy depths:

1. 1. Population partitioning. In those forecasts, SFGs are split into at least two SFG families with distinct evolution: ‘late-type’ discs (normal+starburst) that dominate at $z\!\lesssim\!1$ –1.5, and high-z dusty star-forming systems (often tagged as ‘proto-spheroids’). The sum of two or more evolving families can imprint a mild feature near the transition redshift (Mancuso et al. Reference Mancuso2015). In GHOST we treat SFG as a single, smoothly evolving family, which naturally washes out such structure.

2. 2. How dust enters. Although radio SFR tracers are intrinsically insensitive to dust attenuation, many forecasting frameworks first construct an SFR function from UV, H $\alpha$ , and IR data (the IR term capturing the obscured star formation) and then map SFR $\!\to\!L_{\nu}$ via radio-SFR calibrations (Madau & Dickinson Reference Madau and Dickinson2014; Bonato et al. Reference Bonato2017; Mancuso et al. Reference Mancuso2015). In such models, a distinctly evolving, dust-rich SFG population that grows rapidly near the cosmic SFR peak ( $z\!\sim\!1$ –2) can imprint mild features in n(z) at fixed flux (e.g. late-type versus proto-spheroidal contributions) (Mancuso et al. Reference Mancuso2015). By contrast, our implementation fits the SFG RLF directly at 1.4 GHz with smooth redshift evolution, so a given flux limit projects to a correspondingly smooth redshift distribution.

5.6.2. LoTSS-like selection ( $S_{150}\ge1.5$ mJy)

We emulate the LoTSS DR2 (Hale et al. Reference Hale2024) wide-area selection with an integrated-flux threshold $S_\mathrm{150\,MHz}\ge1.5$ mJy and a peak-S/N threshold, as described in Section 5.1.1. At this depth the sample is still expected to contain a mixed population of radio-loud AGNs and SFGs, and the relative contribution is expected to evolve with redshift because a fixed flux limit implies a rapidly increasing luminosity threshold $L_\mathrm{lim}(z,\alpha)$ ; at $z\gtrsim1$ the selection progressively favours intrinsically luminous systems, which are typically AGNs. We therefore quantify the implied population mixture directly from the selected n(z) and its flux-threshold dependence.

Table 6.

Effect of increasing the 150 MHz flux threshold on the mixture weight ( $w_\mathrm{{SFG}}$ ) and the total-sample median redshift.

Figure 9.

Redshift distributions after applying the flux threshold $S_{150}\!\ge\!1.5\,\mathrm{mJy}$ (Section 5.1). Top: total P(z) from GHOST (black). The green curve reconstructs the total as $w_\mathrm{{SFG}}\,P_\mathrm{{SFG}}(z)+w_\mathrm{AGN}\,P_\mathrm{AGN}(z)$ using the intrinsic flux-limited class fractions $(w_\mathrm{{SFG}},w_\mathrm{AGN})$ ; the blue curve uses the same $P_\mathrm{{SFG}}(z),P_\mathrm{AGN}(z)$ but reweights the class fractions by the acsLoTSS peak–S/N $\ge 7.5$ area acceptance from a log-normal rms model ( $\tilde{\sigma}=83\,\unicode{x03BC}\mathrm{Jy\,beam^{-1}}$ , 95th percentile $171\,\unicode{x03BC}\mathrm{Jy\,beam^{-1}}$ ). The dashed line shows the acsLoTSS-DR2 parametric model. The green and blue reconstructions are essentially identical, indicating that the peak–S/N cut only weakly perturbs the class mix at 1.5 mJy. Middle/bottom: AGN and SFG PDFs (unit area) with medians and IQRs indicated. Relative to acsLoTSS, the simulated total places excess weight at $z\!\lesssim\!0.4$ and underpowers the $z\!\gtrsim\!0.6$ shoulder, suggesting a slightly SFG-heavy mix.

After replicating the $1.5$ mJy sensitivity the selection function removes only a few percent of sources, so the resulting n(z) is driven by the intrinsic class mixture, shown in Figure 9. Where shown, the ‘LoTSS model’ is the parametric p(z) derived from Deep Fields. The total p(z) is low-z dominated (median $z=0.35$ , IQR $0.15$ – $0.83$ ), with a conspicuous SFG contribution at $z\lesssim0.3$ (SFG median $z=0.14$ , IQR $0.09$ – $0.20$ ). The AGN component peaks broadly around $z\sim0.5$ –1 (median $z=0.75$ , IQR $0.47$ – $1.08$ ). Using the measured 150 MHz counts-by-population to set the mixture yields $w_\mathrm{{SFG}}\approx0.42$ . Figure 9 compares the resulting PDFs to the LoTSS model. Relative to the LoTSS model, the total median is lower and the high-z shoulder is underpowered, indicating a SFG-heavy mix at the mJy level and an under-representation of moderate-z AGNs in the total.

Increasing the flux cut shifts the prediction as expected (more AGN-dominated, higher median z), yet even at 10 mJy the median remains below the LoTSS curve (Table 6). The direction and size of the offset are consistent with the RLF diagnostics (Section 5.8); possible population-level origins are deferred to Section 5.6.

5.6.3. VLA-COSMOS split at 1.4 GHz

To mirror VLA-COSMOS (Magliocchetti et al. Reference Magliocchetti, Popesso, Brusa, Salvato, Laigle, McCracken and Ilbert2017), we apply a flux cut $F_{1.4}\ge0.15$ mJy and partition sources using a redshift-dependent luminosity boundary, as detailed in Section 5.1.

At $F_{1.4}\ge150\,\unicode{x03BC}\mathrm{Jy}$ the sample is still expected to contain a mixed population of radio-loud AGN and SFGs, that is, the selection lies in a transition regime rather than the more purely SFG-dominated depths $\lesssim$ $100\,\unicode{x03BC}\mathrm{Jy}$ . In our framework this follows directly from the redshift-dependent luminosity threshold $L_\mathrm{ lim}(z,\alpha)$ implied by a fixed flux limit: $L_\mathrm{lim}$ rises rapidly with z and progressively favours intrinsically luminous systems at $z\gtrsim1$ , which are typically AGN. For this reason, the population fractions inferred from the selected n(z) are expected to evolve with redshift, and we therefore report the predicted AGN/SFG fractions for this flux cut (both integrated and in broad redshift bins) and present the analogous assessment for the 150 MHz and 3 GHz selections.

Figure 10.

VLA-COSMOS-like n(z) with 1.4 G Hz selection. We apply $F_{1.4}\!\ge\!0.15$ mJy and classify sources via $L_{1.4}\gtrless L_\mathrm{cross}(z)$ with $L_\mathrm{cross}=4\pi P_\mathrm{cross}$ and $\log_{10}P_\mathrm{cross}(z)=21.7+z$ ( $z\!\le\!1.8$ ) or $23.5$ otherwise (Section 5.1.2). Top: total P(z) with median and IQR; middle/bottom: AGN-like and SFG-like PDFs. The expected qualitative in redshift – with SFG tightly clustered and AGN more broad – is recovered, but the integrated SFG-like fraction is high, again pointing to a slightly SFG-heavy mix at the mJy level.

The total resulting distribution, as shown in Figure 10 remains low-z dominated (median $z=0.37$ , IQR $0.17$ – $0.87$ ). The AGN-like PDF peaks at $z\sim0.6$ –1 with a tail to $z\sim2$ , while the SFG-like PDF is tightly concentrated at $z\lesssim0.3$ . Integrating above/below $L_\mathrm{cross}(z)$ gives an SFG fraction for this selection of

(29) \begin{equation}f_\mathrm{{SFG}} \equiv \frac{N_\mathrm{{SFG}}}{N_\mathrm{{SFG}}+N_\mathrm{AGN}} \;=\; 0.49 .\end{equation}

Figure 11.

VLA–COSMOS-like n(z) at 3 GHz. GHOST total (top) and class PDFs (middle/bottom) for a $\mathrm{S/N}\!\ge\!5.5$ selection function (Section 5.1.3); the top panel also shows the effective counts-weighted mixture (green). Medians and IQRs are indicated. The total median is low and the SFG fraction is high compared to Hale et al. (Reference Hale, Jarvis, Delvecchio, Hatfield and Novak2018), consistent with the trend that GHOST overweights nearby SFG relative to moderate-z radio-loud AGN.

The selection function recovers the expected qualitative separation in n(z) prevalent in Magliocchetti et al. (Reference Magliocchetti, Popesso, Brusa, Salvato, Laigle, McCracken and Ilbert2017), with SFG tightly clustered at low redshift and a broader AGN distribution, but the integrated SFG share is high for this selection, again pointing to a SFG-heavy mix around the mJy level. Because this diagnostic does not depend on peak-vs-integrated flux (unlike LoTSS), it isolates the population modelling: too many low-z SFG and/or too few moderate-z radio-loud AGNs.

Relative to the COSMOS-based fractions of (Magliocchetti et al. Reference Magliocchetti, Popesso, Brusa, Salvato, Laigle, McCracken and Ilbert2017), our $f_\mathrm{{SFG}}$ is larger in the same sense indicated by the LoTSS test. A hard L-threshold is an imperfect classifier (e.g. it can move luminous SFG at $z\sim1$ into the AGN-like bin), so exact numerical agreement is not expected; the direction of the difference is robust and consistent across selections.

5.6.4. VLA-COSMOS-like selection at 3 GHz

We emulate the $\mathrm{S/N}\ge5.5$ 3 GHz selection (median rms $\simeq2.3\,\unicode{x03BC}$ Jy beam $^{-1}$ ) and compare to the multi-wavelength SFG/AGN classes of Hale et al. (Reference Hale, Jarvis, Delvecchio, Hatfield and Novak2018).

At the Hale et al. (Reference Hale, Jarvis, Delvecchio, Hatfield and Novak2018) depth the applied flux threshold and selection function detailed in Section 5.1 is weak: numerically the acceptances are close to unity for both classes (in our run $\mathrm{acc}_\mathrm{{SFG}}\!\approx\!0.89$ , $\mathrm{acc}_\mathrm{AGN}\!\approx\!0.94$ ), so the predicted n(z) is governed mainly by the intrinsic SFG/AGN mix. The total has median $z=0.38$ . Class medians are $z_\mathrm{med}^\mathrm{{SFG}}=0.30$ (IQR $0.19$ – $0.45$ ) and $z_\mathrm{med}^\mathrm{AGN}=0.80$ (IQR $0.48$ – $1.22$ ). The effective SFG fraction is $w_\mathrm{{SFG}}=0.69$ .

Figure 12.

Normalised redshift PDFs under a MIGHTEE-like 1.4 GHz flux cut: GHOST (solid) compared to the MIGHTEE catalogue (dashed), shown for the total sample and AGN and SFG subsets. Vertical lines indicate medians (solid) and interquartile ranges (dashed).

At this depth and frequency the selection is intrinsically SFG-rich, but the simulated mixture is still somewhat SFG-heavy relative to Hale et al. (Reference Hale, Jarvis, Delvecchio, Hatfield and Novak2018), again implicating the population balance rather than the survey selection function. Hale et al. (Reference Hale, Jarvis, Delvecchio, Hatfield and Novak2018) classify $\sim$ 79% of sources via multi-wavelength diagnostics and report a SFG share of $\simeq$ $3\,704/3\,704+2\,937)\approx0.56$ in the masked field, together with a higher total median z (see their Figure 3). Our Hale-like selection therefore echoes the acsLoTSS and Magliocchetti results: GHOST places relatively more weight in nearby SFG and underweights moderate-z radio-loud AGN. The sign and scale of the difference are consistent with the RLF and clustering comparisons.

5.7. Differential and Euclidean-normalised number counts

This section tests the total source density in GHOST against the literature via Euclidean-normalised differential counts, $S^{2.5}\,dN/dS$ , at 150 MHz, 1.4 GHz, and 3 GHz. Because number counts integrate over luminosity and redshift, they primarily probe the combination of luminosity-function normalisations and radio SEDs, largely independent of the host-halo assignment used for clustering.

Figure 13 shows that over nearly seven decades in flux density the simulation follows the observed trends at all three frequencies. At the bright end ( $S\gtrsim10$ mJy), the counts are close to Euclidean ( $S^{2.5}\,dN/dS\!\approx\!\mathrm{const}$ ), consistent with a population dominated by radio-loud AGN.

Figure 13.

Euclidean-normalised differential counts, $S^{2.5}\,{d}N/{d}S$ , for GHOST (black markers) compared to measurements at 150 MHz (Franzen et al. Reference Franzen2016; Mandal et al. Reference Mandal2021), 1.4 GHz (Bondi et al. Reference Bondi, Ciliegi and Schinnerer2008; Bridle et al. Reference Bridle, Davis, Fomalont and Lequeux1972; Ciardi & Loeb Reference Ciardi and Loeb2000; Fomalont et al. Reference Fomalont, Kellermann, Cowie, Capak, Barger, Partridge, Windhorst and Richards2006; Gruppioni et al. Reference Gruppioni, Mignoli and Zamorani1999; Hopkins et al. Reference Hopkins, Afonso, Chan, Cram, Georgakakis and Mobasher2003; Ibar et al. Reference Ibar, Ivison, Biggs, Lal, Best and Green2009; Kellermann et al. Reference Kellermann, Fomalont, Mainieri, Padovani, Rosati, Shaver, Tozzi and Miller2008; Mitchell & Condon Reference Mitchell and Condon1985; Owen & Morrison Reference Owen and Morrison2008; Richards Reference Richards2000; Seymour et al. Reference Seymour2008; White et al. Reference White, Becker, Helfand and Gregg1997; de Zotti et al. Reference de Zotti, Massardi, Negrello and Wall2010) and 3 GHz (Smol $\breve{c}$ ić et al. Reference Smolčić2017). Error bars are $1\sigma$ measurement uncertainties, incorporating the estimations of Gehrels (Reference Gehrels1986) if $N\lt20$ for any bin; shaded bands (if shown) indicate the field-to-field variance from GHOST mock realisations with the same area and flux cuts as the data. Residual panels show $(\text{model}-\text{data})/\text{data}$ .

A tension appears around the classical sub-mJy upturn, where the data rise due to the emergence of SFG. In our model this rise is weaker at $S\sim0.1$ – $0.3$ mJy. This behaviour is fully consistent with the adopted SFG RLF: the radio-first, smoothly evolving SFG luminosity function underpredicts the luminous tail of SFG at $z\gtrsim1$ , which suppresses the upturn when projected to flux. Since total counts depend on the RLFs and SEDs rather than on host assignment, abundance matching does not drive this feature.

5.8. RLF

We describe the class-dependent RLF by $\Phi_c(L_\nu,z)$ in units of $\mathrm{Mpc}^{-3}\,\mathrm{dex}^{-1}$ , where $c\in\{\mathrm{SFG},\mathrm{SSAGN},\mathrm{FSRQ},\mathrm{BLLac}\}$ . Observed flux density at $\nu_\mathrm{obs}$ is related to rest-frame luminosity through the standard K-correction,

(30) \begin{equation} L_\nu(z,\alpha_c) = 4\pi D_L^2(z)\,S_{\nu_\mathrm{obs}}\, \Big(\frac{\nu_\mathrm{ref}}{\nu_\mathrm{obs}}\Big)^{-\alpha_c}(1+z)^{1+\alpha_c}, \end{equation}

so that, for a survey flux limit $S_\mathrm{lim}$ , the luminosity boundary within a redshift slice is

(31) \begin{equation} L_{\mathrm{lim}}(z,\alpha_c) = 4\pi D_L^2(z)\,S_{\mathrm{lim}}\, \Big(\frac{\nu_\mathrm{ref}}{\nu_\mathrm{obs}}\Big)^{-\alpha_c}(1+z)^{1+\alpha_c},\end{equation}

with $L_\nu\propto\nu^{-\alpha}$ .

Forward (projected) model. When comparing to data taken over a finite redshift bin $z\in[z_1,z_2]$ we project the intrinsic RLF under the same selection,

(32) \begin{equation}\begin{aligned}v c&\Phi^\mathrm{proj}_c(L\,|\,z_1{:}z_2) = \\&\frac{\int_{z_1}^{z_2} \Phi_c(L,z)\, \Theta \big[L-L_\mathrm{lim}(z,\alpha_c)\big]\, ({{d}V}/{d}z)\,{d}z} {\int_{z_1}^{z_2} ({{d}V}/{d}z)\,{d}z} .\end{aligned}\end{equation}

These orange dashed curves in Figures 16 and 17 should be regarded as the baseline ‘model-to-data’ comparison.

Binned estimator from the catalogue ( $V_\mathrm{eff}$ ). The blue points are a direct estimate from the catalogue in the $[z_1,z_2]\times[L_1,L_2]$ cell using an effective volume that accounts for partial visibility of the L-bin across the slice. Define the bin width in dex by $\Delta\log_{10}L=\log_{10}\!(L_2/L_1)$ . At any redshift z, the fraction of the luminosity bin that is observable is

(33) \begin{align}&f_\mathrm{vis}(z;\,L_1,L_2)= \nonumber \\[2pt] &\begin{cases}1 &L_\mathrm{lim}(z,\alpha_c)\le L_1,\\[2pt]\displaystyle\frac{\log_{10}\!\big(L_2/L_\mathrm{lim}(z,\alpha_c)\big)} {\log_{10}\!\big(L_2/L_1\big)} &L_1\lt L_\mathrm{lim}(z,\alpha_c) \lt L_2,\\[2pt]0 &L_\mathrm{lim}(z,\alpha_c)\ge L_2\, .\end{cases}\end{align}

Figure 14.

Fraction of (L, z) cells with zero effective volume, $V_\mathrm{eff}=0$ , in each redshift bin for the raw selection. Bars show the average masked fraction across AGN and SFG classes; these are indistinguishable within $\lesssim$ $0.5\%$ in every bin because the mask is set by the survey flux limit (and K-correction) rather than population physics. The fraction rises monotonically from $\lesssim$ 15% at $z\le0.3$ to $\simeq$ 40% by $z\sim1.8$ , reflecting the shrinking accessible luminosity range in a flux-limited survey. Cells with $V_\mathrm{eff}=0$ are excluded from the panel-averaged RLFs.

The corresponding effective comoving volume for the cell is

(34) \begin{equation} V_\mathrm{eff}(L_1,L_2;\,z_1{:}z_2) \;=\; \Omega \int_{z_1}^{z_2} \frac{{d}V}{{d}z}(z)\, f_\mathrm{vis}(z;\,L_1,L_2)\,{d}z, \end{equation}

with survey solid angle $\Omega$ (sr). If $N_c$ catalogue objects of class c fall in the cell, the minimum-variance, Poisson estimator is

(35) \begin{align} \widehat{\Phi}_c(L\,|z_1{:}z_2) \; &=\; \frac{N_c}{V_\mathrm{eff}(L_1,L_2;\,z_1{:}z_2)\,\Delta\log_{10}L}\,, \\[-10pt] \nonumber \end{align}

(36) \begin{align} \sigma^2\!\left[\widehat{\Phi}_c\right] \;&=\; \frac{N_c}{\big[V_\mathrm{eff}\,\Delta\log_{10}L\big]^2}. \\[12pt] \nonumber \end{align}

If $N\lt20$ for any bin, we use the estimators of Gehrels (Reference Gehrels1986) in constructing the confidence interval. We plot the estimate at the bin centre with vertical error bars (from Equation 36) and omit bins where $V_\mathrm{eff}=0$ . Horizontal bars indicate $\Delta\log_{10}L$ . For empty bins with $V_\mathrm{eff}\gt0$ we show one-sided upper limits (Gehrels-style) when useful.

Masked fraction of ( L , z ) cells. Figure 14 shows, for each redshift bin, the fraction of (L, z) cells that are masked because their effective survey volume is zero, $V_\mathrm{eff}(z,L)=0$ . A cell is masked precisely when the visibility fraction in Equation (33) vanishes, $f_\mathrm{vis}(z;\,L_1,L_2)=0$ ; in those cells we set $\Phi$ to NaN and exclude them from the RLF averages. As expected for a flux-limited sample, the masked fraction is $\simeq 0$ in the nearest slice $(0 \lt z\le0.1)$ and then rises monotonically with redshift as the flux boundary $L_\mathrm{lim}(z,\zeta)\propto D_L^2(z)\,(1+z)^{\zeta-1}$ moves to higher luminosities. Numerically, the fractions are $\approx$ 15% at $0.10\text{-}0.30$ , $\approx$ $23\text{-}28\%$ at $0.30\text{-}0.70$ , $\approx$ $32\text{-}38\%$ at $0.70\text{-}1.40$ , and $\approx$ $40\text{-}42\%$ by $z\simeq1.8$ . The AGN and SFG panels are nearly identical because the $V_\mathrm{eff}$ mask is set by survey geometry and the adopted spectral index used in the visibility correction, not by population physics. Thus the increasing masked fraction simply quantifies the shrinking accessible L-range with redshift and explains the growing white space in the high-z RLF panels; it is a geometric selection effect, not a loss of simulated sources nor a model shortcoming.

Figure 15.

Comparison of 1.4 GHz radio luminosity-SFR calibrations. Dashed curves show literature relations (Molnár et al. Reference Molnár2021; Heesen et al. Reference Heesen2022; Cook et al. Reference Cook2024; Matthews et al. Reference Matthews2024; Davies et al. Reference Davies2017; Murphy et al. Reference Murphy, Chary, Dickinson, Pope, Frayer and Lin2011; Bell et al. Reference Bell, McIntosh, Katz and Weinberg2003). The solid red line (‘this work’) is our adopted mapping. At high SFR the slope of our relation is shallower than several published fits, compressing the bright tail when mapping a given $\phi_\mathrm{ SFR}(z)$ into $L_{1.4}$ and thereby yielding fewer SFG at $L_{1.4}\!\gtrsim\!10^{23\text{-}24}\,\mathrm{W\,Hz^{-1}}$ . This behaviour contributes to the lower bright-end SFG space densities seen at $z\gtrsim1$ in Figure 17. The halo mass range and abundance matching limits the upper SFR range in practice.

Figure 16.

AGN radio luminosity functions at 1.4 GHz in redshift slices (raw cut). Each panel shows $\Phi(L_\nu)\,(\mathrm{Mpc}^{-3}\,\mathrm{dex}^{-1})$ for the stated redshift bin. Blue circles: RLF measured from the simulation using the effective-volume method, $\Phi = N/(V_\mathrm{eff}\,\Delta\log_{10}L)$ , where $V_\mathrm{eff}$ accounts for the flux limit within the bin (Section 5.8; Equation 36). Error bars are $1\sigma$ Poisson uncertainties propagated through $V_\mathrm{eff}$ ; cells with $V_\mathrm{eff}=0$ are omitted. Purple dashed: forward (projected) model $\Phi^\mathrm{proj}(L|z_1:z_2)$ obtained by averaging the intrinsic AGN RLF over the finite $\Delta z$ of the panel under the same flux selection (Equation 32). Points in the background: literature measurements at 1.4 GHz (legend in the top-left panel). Across the ten slices the simulation (blue) and the forward model (purple) show excellent agreement in both normalisation and shape; a mild feature near $L_\nu\!\sim\!10^{26}\,\mathrm{W\,Hz^{-1}}$ reflects the combined contribution of flat-spectrum (FSRQ/BLLac) and steep-spectrum (SSAGN) components with luminosity-dependent evolution (see Section 5.8). Axis limits in each panel are set by the extent of the observational points to aid readability: Squares indicate radio loud AGN; diamonds indicate radio quiet AGN. Filled markers are fully within the z-bin; open markers partially overlap the bin.

Figure 17.

SFG radio luminosity functions at 1.4 GHz in redshift slices (raw cut). As in the AGN figure, blue circles show the simulation RLF recovered with the effective-volume estimator $\Phi=N/(V_\mathrm{eff}\,\Delta\log_{10}L)$ (with $1\sigma$ Poisson errors), and orange dashed curves are the forward (projected) model evaluated over the same $\Delta z$ and flux selection. Background points: literature measurements at 1.4 GHz (legend in the top-left panel). At $z\!\lesssim\!1$ the model and simulation track the observed faint end well. At higher redshift the bright tail in the simulation/model sits below some compilations: this follows naturally from our adopted L-SFR mapping (shallower high-SFR slope), the tightening luminosity threshold imposed by the flux limit within each $\Delta z$ slice, and the resulting reduction of accessible high-SFR systems (see Section 5.8). As before, panel limits are set by the extent of the observational points; filled markers are fully within the z-bin; open markers partially overlap the bin.

Results at $\boldsymbol{1.4}$ GHz. Figures 16 and 17 show excellent agreement between the projected model (orange dashed) and the catalogue-based, $V_\mathrm{eff}$ -corrected estimates (blue points):

• • AGN: The match is tight across all redshift slices and over the full dynamic range where the survey is complete. This validates both the evolutionary form of the AGN RLF and the class mix encoded in the model.

• • SFG: At $z\gtrsim 1$ the blue points fall below some literature bright-end measurements. This is a model outcome: our adopted $L_{1.4}$ –SFR calibration is shallower at high SFR (see Figure 15), compressing the bright tail and yielding fewer $L_{1.4}\gtrsim 10^{23{-}24}\,\mathrm{W\,Hz^{-1}}$ systems when an evolving SFR function is mapped into radio luminosity. The same behaviour is consistent with the divergence of the total counts at the faint flux end where SFG dominate.

A shallow ‘bump/inflection’ occasionally visible in the orange curves near $10^{26}\,\mathrm{W\,Hz^{-1}}$ is naturally produced by the sum of AGN subclasses and their evolution; the double-power-law knee of the steep-spectrum tail, combined with the beamed (FSRQ/BLLac) component beginning to dominate the bright end with a slightly flatter local slope, and the finite-width slice average in Equation (32). A per-class decomposition of $\Phi^\mathrm{proj}$ confirms that the total slope locally reflects the crossover of SSAGN and beamed contributions. The feature is therefore not required by the data and does not affect the excellent overall agreement.

5.9. The effect of observational parameters

The results discussed in this work use the GHOST intrinsic catalogue prior to applying survey-specific observational effects. In real radio-continuum analyses, going from an intrinsic sky to a detected catalogue depends on the synthesised beam, imaging/source-finding choices, and the confusion environment, which vary between instruments and releases; we therefore treat a fully general observational model as future work, and here summarise the main consequences of its omission.

5.10. Impact of SFG modelling choices

The impact of alternative SFG modelling choices on number counts is quantified through two cross-checks: an IR-calibrated construction used to derive an effective radio RLF, and the inclusion of thermal dust emission as an SED component when generating multi-frequency fluxes (Section 3.2.1).

Figures 22–23 compare the redshift and flux distributions predicted by the IR-calibrated construction (Appendix A.7) against the fiducial radio-calibrated prescription. Over the redshift range where most SFGs contributing to typical radio selections lie ( $z\lesssim 2$ ), the two constructions yield broadly similar N(z) and N(S) shapes, with differences at the factor-of-few level that decrease with redshift. The dominant qualitative difference is the high-redshift behaviour: in the tested IR-calibrated prescription the SFG density turns over more sharply at $z\gtrsim 3$ –4, whereas the radio-calibrated model retains a higher high-z tail. This behaviour reflects the fact that the IR-calibrated construction inherits assumptions about IR luminosity evolution and the IR-radio mapping (Yun et al. Reference Yun, Reddy and Condon2001; Cai et al. Reference Cai2013), while the fiducial GHOST model is explicitly calibrated to radio observables.

Figure 25 quantifies the effect of adding thermal dust emission as an SED component (Appendix A.7) on SFG number counts at high observing frequency. Across most redshifts the impact on N(z) is negligible, consistent with synchrotron/free-free dominating the SFG spectrum below $\sim$ 10–20 GHz (Appendix A.6). In Bonaldi et al. (Reference Bonaldi, Bonato, Galluzzi, Harrison, Massardi, Kay, De Zotti and Brown2019), the rest-frame SFG SED is dust-dominated only at $\nu \gtrsim 100$ GHz, implying that dust becomes important in the observer frame only at frequencies of a few tens of GHz for high-redshift galaxies; they implement this by adding dust emission using the model of Cai et al. (Reference Cai2013) and the same SFR- $L_\mathrm{{IR}}$ calibration. Our results follow the same trend: the dust term affects only the extreme high-z tail at the highest frequencies, where the observer-frame band probes substantially higher rest-frame frequencies and the modified blackbody emission rises rapidly.

This conclusion is also consistent with (Wilman et al. Reference Wilman2008), which adopts a modified blackbody dust component with $(T,\beta)=(45\,\mathrm{K},1.5)$ normalised via the IR-radio correlation and finds a discernible impact on the observed SED only at very high redshift and at the highest output frequency. For the metre-centimetre regimes emphasised in this work, omitting thermal dust emission therefore has negligible impact on the reported results, while including it improves SED completeness for multi-frequency catalogues.