11institutetext: Department of Physics, School of Physics and Electronics,
Hunan Normal University, Changsha 410081, China 22institutetext: Key
Laboratory of Low Dimensional Quantum Structures and Quantum Control, Hunan
Normal University, Changsha 410081, China 33institutetext: Synergetic
Innovation Center for Quantum Effects and Applications and Institute of
Interdisciplinary Studies, Hunan Normal University, Changsha, Hunan 410081,
China 44institutetext: Yunnan Observatories, Chinese Academy of Sciences,
Kunming 650216, China 55institutetext: Center for Astronomical Mega-Science,
Chinese Academy of Sciences, Beijing 100012, China 66institutetext: Key
Laboratory for the Structure and Evolution of Celestial Objects, Chinese
Academy of Sciences, Kunming 650216, China

# Constraints on the Cosmic Star Formation History via a new modeling of the
SFG radio luminosity function

Wenjie Wang 1122 Zunli Yuan Corresponding author 1122<EMAIL_ADDRESS>Hongwei
Yu Corresponding author 112233<EMAIL_ADDRESS>Jirong Mao 445566

(Received; accepted)

###### Abstract

Context. Radio wavelengths offer a unique possibility of tracing the total
star-formation rate (SFR) in galaxies, both obscured and unobscured. To prob
the dust-unbiased star-formation history, an accurate measurement of the radio
luminosity function (LF) for star-forming galaxies (SFGs) is crucial.

Aims.  We make use of an SFG sample (5915 sources) from the Very Large Array
(VLA) COSMOS 3 GHz data to perform a new modeling of the radio LF. By
integrating the analytical LF we aim to calculate the cosmic SFR density
(SFRD) history since $z\sim 5$.

Methods.  For the first time, we use both models of the pure luminosity
evolution (PLE) and joint luminosity+density evolution (LADE) to fit the LFs
directly to the radio data using a full maximum likelihood analysis,
considering the sample completeness correction. We also incorporate the
updated observations on local radio LFs and radio source counts into the
fitting process to obtain additional constraints.

Results.  We find that the PLE model is not applicable to describe the
evolution of the radio LF at high redshift ($z>2$). By construct, our LADE
models can successfully fit a large amount of data on radio LFs and source
counts of SFGs from recent observations. Our SFRD curve shows a good fit to
the SFRD points derived by previous radio estimates. In view of that our radio
LFs are not biased in contrast to previous studies performed by fitting the
$1/V_{\rm max}$ LF points, our SFRD results should be an improvement to these
previous estimates. Below $z<2$, our SFRD matches well with the multi-
wavelength compilation from Madau & Dickinson, while our SFRD turns over at a
slightly higher redshift ($2<z<2.5$) and falls more rapidly out to high
redshift.

###### Key Words.:

galaxies: evolution – galaxies: star formation – galaxies: luminosity
function, mass function – radio continuum: galaxies

##  1 Introduction

Understanding the formation and evolution of galaxies through the cosmic time
is a major quest of modern cosmology. One of the most fundamental processes
driving the evolution of galaxies is star formation. Star formation rate
density (SFRD), defined as the amount of stars formed per year in a unit of
cosmological volume, is thus a critical parameter for galaxies (e.g. Lara-
López et al., 2010). In the past decades, a lot of research based on
multiwavelength aspects has been devoted to an accurate measurement of the
cosmic evolution of SFRD. Remarkable progress has been made on this, and the
SFRD is now well understood up to $z\sim 3$, when the universe was no more
than 3 Gyr old (see Madau & Dickinson, 2014, for an exhaustive review). A
consensus has been achieved regarding recent history, where the SFRD increases
significantly with redshift, reaching a peak at redshift $z\sim 2$, an epoch
known as “cosmic noon” (e.g., Novak et al., 2017; van der Vlugt et al., 2022;
Enia et al., 2022). The picture becomes less clear at higher redshifts. Some
studies (mainly based on UV-selected sources) show a steep decline in the SFRD
(Bouwens et al., 2015; McLeod et al., 2016; Ishigaki et al., 2018), while
studies performed at radio or submillimeter wavelengths show a flatter SFRD at
z ¿ 3 (Gruppioni et al., 2013; Rowan-Robinson et al., 2016; Novak et al.,
2017; Gruppioni et al., 2020). Such uncertainty in the evolution of SFRD at
early cosmic epochs hinders our understand of the core mechanism that governs
star formation rate (SFR) histories of individual galaxies.

The SFR of galaxies can be traced at multiple wavebands, each tracer having
its own advantages and disadvantages (Kennicutt Jr, 1998). In dust-free
environments, ultraviolet (UV) light originating primarily from young massive
stars, serves as the most direct tracer of SFR. UV light can be used to
constrain the unobscured star formation out to very high redshifts (e.g.,
McLure et al., 2013; Bowler et al., 2015; Finkelstein et al., 2015; McLeod et
al., 2015; Parsa et al., 2016; Oesch et al., 2018; Ono et al., 2018; Adams et
al., 2020; Bowler et al., 2020; Bouwens et al., 2021). However, UV
observations suffer from dust absorption, resulting in the SFR measurements
made at these wavelengths being underestimated (e.g., Smail et al., 1997;
Bouwens et al., 2009; Riechers et al., 2013; Dudzevičiūtė et al., 2020). When
the dust absorbs UV radiation, it gets heated and reradiates the energy at
far-infrared (FIR) wavelengths. Therefore, FIR emission is ideal to trace SFR
in dust-rich environments (see Kennicutt Jr, 1998). Unfortunately, FIR
observations can suffer from poor resolution and source blending.

Now it seems that a very promising tracer is using deep radio continuum
observations, which offer a unique possibility of tracing the total SFR in
galaxies, both obscured and unobscured. As such, they may provide the most
robust measurement of the star-formation history of the Universe (Jarvis et
al., 2015). Radio continuum emission, not affected by dust obscuration, is
also an end product of the formation of massive stars (e.g., van der Vlugt et
al., 2022). After these short-lived ($\tau\leqslant 3\times 10^{7}$ yr,
Matthews et al., 2021a) stars undergo supernova explosions, the expanding
remnants can accelerate the cosmic ray electrons and give rise to synchrotron
radiation at a typical frequency of $<30$ GHz (e.g. Sadler et al., 1989;
Condon, 1992; Clemens et al., 2008; Tabatabaei et al., 2017). Radio emission
triggered by the above process is empirically found to correlate well with the
FIR emission of star-forming galaxies (SFGs), known as the FIR-radio
correlation. This correlation holds over five orders of magnitude in
luminosities and extends to high redshifts (Helou et al., 1985; Yun et al.,
2001; Bell, 2003), although the redshift evolution is controversial (Jarvis et
al., 2010; Sargent et al., 2010; Magnelli et al., 2015; Calistro Rivera et
al., 2017; Delhaize et al., 2017). The FIR-radio correlation can be used to
calibrate radio luminosity as a tracer of SFR (Condon, 1992).

In the past several years deep radio surveys reaching sub-mJy detection limits
have been emerging as a powerful tool to investigate the cosmic evolution of
SFGs (e.g., van der Vlugt et al., 2022; Enia et al., 2022; Malefahlo et al.,
2022; Bonato et al., 2021a, b; Ocran et al., 2020b; Upjohn et al., 2019; Ceraj
et al., 2018; Bonato et al., 2017; Novak et al., 2017; Smolčić et al., 2009).
These studies generally measured the radio luminosity functions (LFs) of SFGs,
and the SFRD can then be estimated by taking the luminosity weighted integral
of the radio LF (e.g., van der Vlugt et al., 2022). As for the form of radio
LFs, most of them assumed pure luminosity evolution (PLE, e.g. Smolčić et al.,
2009; Novak et al., 2017; Ocran et al., 2020b; Malefahlo et al., 2022). Very
recently, van der Vlugt et al. (2022) combined the COSMOS-XS survey and Very
Large Array (VLA)-COSMOS 3 GHz data sets to constrain a radio LF with both
luminosity and density evolution. Their analytical LFs are obtained through
fitting the LF points given by the Schmidt (1968) $1/V_{\rm max}$ estimator.
This semi-parametric method has also been adopted by almost all the existing
studies (e.g., Smolčić et al., 2009; Novak et al., 2017; Cochrane et al.,
2023). However, given the ordinary precision in the $1/V_{\rm max}$ estimate,
the LF points themselves have errors, fitting to them will propagate the
uncertainties to the analytical LFs. Besides, the result would be dependent on
the choice of binning in the $1/V_{\rm max}$ method (see Fan et al., 2001). We
believe that a more reliable approach to obtain the analytical LFs is using a
full maximum likelihood analysis (e.g., Willott et al., 2001).

In this paper we make use of the VLA-COSMOS 3 GHz data (Smolčić et al., 2017a)
to measure the radio LFs of SFGs. We will use both models of PLE and joint
density+luminosity evolution to fit the SFG LFs directly to the radio data
using a full maximum likelihood analysis. We aim to perform a comprehensive
parametric research on the radio LF of SFGs, by means of constraints from
multiple observational data. Finally, we can probe the dust-unbiased SFRD up
to redshift of $z\sim 5$.

The structure of the paper is outlined below. In Section 2, we briefly
describe the data used. In Section 3, we present the method of constraining
the LFs with redshift. In Section 4, we derive our radio LFs evolution through
cosmic time and compare it to those in the literature. In Section 5, we
calculate the evolution of the cosmic SFRD using the LF models we derived and
compare it to the literature. Section 6 summarizes and concludes this work.

Throughout the paper we have used the flat concordance Lambda cold dark matter
($\Lambda$CDM) cosmology with the following parameters: Hubble constant
$H_{0}=70\rm kms^{-1}Mp^{-1}$, dark energy density $\Omega_{\Lambda}=0.7$, and
matter density $\Omega_{m}=0.3$. We assume the Chabrier (2003) initial mass
function (IMF) to calculate SFRs.

##  2 Sample

Figure 1:  Redshift distribution ($top$) and the scatter plot ($bottom$) of
our SFG sample. The red dashed line indicates the flux limit of
$S_{1.4\leavevmode\nobreak\ \rm GHz}=11.5\leavevmode\nobreak\ \mu\mathrm{Jy}$.
Note that three sources are below the flux limit line, thus they are excluded
from our following analysis.

In our study, we use the same sample of SFGs as presented in Novak et al.
(2017), which was compiled from the continuum data and source catalog release
of the VLA-COSMOS 3 GHz Large Project survey (Smolčić et al., 2017a). The
sample of SFGs was selected via radio emission and complemented with ancillary
data from the comprehensive multiwavelength coverage of COSMOS. The data
analysis and multi-band association procedure are fully described in Novak et
al. (2017), and we refer readers to that publication for a complete
description. Here we summarize some key points about the sample.

The VLA-COSMOS 3 GHz Large Project survey utilized 384 hours of VLA A+C array
observations in the S band to obtain radio data. The survey covered a uniform
rms noise of 2.3 $\mu$Jy $\text{beam}^{-1}$ and had an angular resolution of
$0^{\prime\prime}_{\cdot}75$ across the 2 square degrees of COSMOS. The final
catalog contains 10 830 radio sources. Taking into account the fraction of
spurious sources, an $11\%$ incompleteness in the counterpart sample is
estimated. Therefore, a total of 7729 radio sources with assigned COSMOS 2015
counterparts were used. About 35% of these radio sources have spectroscopic
redshifts, and photometric redshifts were used for the remainder of the
sample. According to Delvecchio et al. (2017), sources were classified as
radio-excess if the value of $r$ deviates more than $3\sigma$ from the peak of
the distribution obtained as a function of redshift, i.e.,

$r=\log\left(\frac{L_{1.4\mathrm{GHz}}\left[\mathrm{W}\mathrm{Hz}^{-1}\right]}{SFR_{\mathrm{IR}}\left[M_{\odot}\mathrm{yr}^{-1}\right]}\right)>22\times(1+z)^{0.013}.$
(1)

According to this criterion, they were able to distinguish 1814 sources (23%)
that are primarily emitting due to AGN activity in the radio. By assuming a
spectral index of $\alpha=0.7$, they estimated the rest-frame 1.4 GHz
luminosity with the known 3 GHz flux density as

$\displaystyle L_{1.4\leavevmode\nobreak\ \mathrm{GHz}}=\frac{4\pi
D_{L}^{2}}{(1+z)^{1-\alpha}}\left(\frac{3\leavevmode\nobreak\
\mathrm{GHz}}{1.4\leavevmode\nobreak\
\mathrm{GHz}}\right)^{\alpha}S_{3\leavevmode\nobreak\ \mathrm{GHz}},$ (2)

where $D_{L}$ represents the luminosity distance at redshift $z$,
$S_{3\leavevmode\nobreak\ \rm GHz}$ represents the 3 GHz flux density observed
by the survey.

In summary, the sample consists of 5915 SFGs that do not exhibit radio excess.
All the sources have a (spectroscopic or photometric) redshift. The redshift
distribution of the SFG sample as well as its scatter plot are shown in Figure
1. The red dashed line indicates the flux limit of $S_{1.4\leavevmode\nobreak\
\rm GHz}=11.5\leavevmode\nobreak\ \mu\text{Jy}$. We note that three sources
are below the flux limit line, and thus they are excluded from our following
analysis.

##  3 Methods

###  3.1 LF and likelihood function

The LF $\Phi(z,L)$ is a measurement of the number of sources per unit comoving
volume per unit logarithmic luminosity interval:

$\displaystyle\Phi(z,L)=\frac{d^{2}N}{dVd\log_{10}L}.$ (3)

Given an analytical form with parameters $\mathbf{\theta}$ for the LF,
$\Phi(z,L|\mathbf{\theta})$, the maximum likelihood solution to
$\mathbf{\theta}$ is obtained by minimizing the negative logarithmic
likelihood function $S\equiv-2\ln\mathscr{L}$. Following Marshall et al.
(1983) and (Fan et al., 2001), the likelihood function can be written as

$\displaystyle\begin{aligned}
S=&-2\\!\\!\sum_{i}^{n}\\!\ln[\Phi(z_{i},L_{i})p(z_{i},,L_{i})]\\\
&+2\\!\\!\\!\int\\!\\!\\!\\!\int_{W}\\!\\!\\!\Phi(z,L)p(z,L)\Omega\frac{dV}{dz}dzdL,\end{aligned}$
(4)

where $p(z,L)$ is the selection probability of the SFG as a function of
redshift and luminosity. The inclusion of the selection probability in
equation (4) accounts for the fact that the sample is incomplete near the flux
limit. The symbol $\Omega$ represents the solid angle covered by the survey,
and $dV/dz$ denotes the differential comoving volume per unit solid angle, as
defined by Hogg (1999).

For our SFG sample, $p(z,L)$ can be estimated by

$p(z,L)=C_{\text{radio}}[S_{3\leavevmode\nobreak\ \text{GHz}}(z)]\times
C_{\text{opt}}(z),$ (5)

where $C_{\text{radio}}$ is the completeness of the VLA-COSMOS 3 GHz radio
catalog as a function of the flux density $S_{3\leavevmode\nobreak\
\text{GHz}}$, and $C_{\text{opt}}$ is the completeness owing to radio sources
without assigned optical-NIR counterparts (Novak et al., 2017). We adopt the
calculations of $C_{\text{radio}}$ and $C_{\text{opt}}$ given by Novak et al.
(2017), and refer the interested reader to their Fig. 2 for more details.

To estimate the integration term in Equation (4), one needs to find the
function values for $p(z,L)$ at given pairs of (z,L). We achieve this using an
interpolation method. Firstly, we set a two-dimensional (2D) grid of $50\times
50$ in the $\log L-z$ space. For each grid point ($\log L_{i},z_{i}$), we can
derive its flux density $S_{i}$ via Equation (2) assuming $\alpha=0.7$. Then
we can estimate the corresponding $C_{\text{radio}}$ and $C_{\text{opt}}$
through a one-dimensional linear interpolation method using the observed value
from Novak et al. (2017). Finally, we have the values for $p(z,L)$ at the
$50\times 50$ grid points, which are used to perform the 2D linear
interpolation to estimate the function value of $p(z,L)$.

Following Willott et al. (2001) and Yuan et al. (2017), we incorporate the
updated observations on local radio LFs and source counts into the fitting
process to obtain additional constraints. The local radio LF and the source
counts are one-dimensional functions, and their $\chi^{2}$ value is calculated
as

$\displaystyle\begin{aligned} \chi^{2}=\sum_{i=1}^{n}\left(\frac{f_{\text{data
}i}-f_{\bmod i}}{\sigma_{\text{data }i}}\right)^{2},\end{aligned}$ (6)

where $f_{\text{data }i}$ represents the value of the data in the $i$th bin,
$f_{\text{mod }i}$ and $\sigma_{\text{data }i}$ are the model value and data
error in the ith bin, respectively. Since $\chi^{2}$ is related to the
likelihood by $\chi^{2}=-2\ln(\text{likelihood})$ (i.e., the same form as $S$,
Willott et al., 2001), we can define a new function $S_{\text{all}}$, which
combines the constraints from all three types of data. The expression is as
follows

$\displaystyle
S_{\text{all}}=S+\chi^{2}_{\text{LRLF}}+A_{0}\chi^{2}_{\text{SC}},$ (7)

where $\chi^{2}_{\text{LRLF}}$ and $\chi^{2}_{\mathrm{SC}}$ denote the value
of $\chi^{2}$ for the local radio LFs and source counts, respectively. Because
we use three different types of data to estimate $S_{\mathrm{all}}$, we need
to balance the statistical weight for each term in Equation (7). We chose an
$A_{0}$ so that the value of $A_{0}\chi^{2}_{\mathrm{SC}}$ approximately
equals to that of $\chi^{2}_{\text{LRLF}}$. This yields values of about 10-40
for our calculations. We find that varying $A_{0}$ does not significantly bias
our final results (also see Kochanek, 1996). Using Equation (7), we can obtain
the best-fit parameters for LFs by numerically minimizing the objective
function $S_{\text{all}}$. Here we adopt a Bayesian method as in our previous
papers (e.g., Yuan et al., 2016). It enables us to determine the best
estimates for the model parameters and their probability distribution (also
see Lewis & Bridle, 2002; Yuan et al., 2017). We use uniform (so-called
“uninformative”) priors on the parameters, and employ the MCMC sampling
algorithm, available in the Python package emcee (Foreman-Mackey et al.,
2013), to estimate the best-fit parameters.

Figure 2:  Local radio LF at 1.4 GHz of SFGs from several surveys with
different observed areas and sensitivities (colored data points). The colored
lines show the fits to the combined data from our models.

###  3.2 Local LFs and Radio Source Counts

The local LFs at 1.4 GHz have been well determined for SFGs, due to the
combined use of large radio surveys, such as NVSS (NRAO VLA Sky Survey) and
FIRST (Faint Images of the Radio Sky at Twenty centimeteres), and large-area
spectroscopic surveys. In this work, we simultaneously use the local SFG LFs
from Condon et al. (2002), Best et al. (2005), Mauch & Sadler (2007), and
Condon et al. (2019) (see Figure 2) to calculate $\chi^{2}_{\text{LRLF}}$ in
Equation (7).

Besides the local LFs, the observed radio source counts can also provide an
important constraint to the modeling of SFG LFs. In the past several years
deep radio surveys reaching sub-mJy detection limits have been emerging,
enabling the faint source counts to be investigated (e.g., Smolčić et al.,
2017b; Ocran et al., 2020a; Mandal et al., 2021; van der Vlugt et al., 2021;
Matthews et al., 2021b; Hale et al., 2023). The source counts, denoted as
$n(S_{\nu})$, represent the number of sources per flux density ($S_{\nu}$) per
steradian. Its shape is closely related to the evolutionary properties of the
source as well as the geometry of the universe (Padovani, 2016). Typically,
the counts are Euclidean normalized by multiplying $S^{2.5}_{\nu}$ (e.g., de
Zotti et al., 2010). According to Padovani (2016) and Yuan et al. (2017), we
can relate the source counts of SFGs with their LF using the following
equation:

$\displaystyle\begin{aligned}
\frac{n(S_{\nu})}{4\pi}=4\pi\frac{c}{H_{0}}\int_{z_{\text{min}}(S_{\nu})}^{z_{\text{max}}(S_{\nu})}\frac{\Phi(z,L(S_{\nu},z))D_{L}^{4}(z)dz}{(1+z)^{(3-\alpha)}\sqrt{\Omega_{\mathrm{m}}(1+z)^{3}+\Omega_{\Lambda}}},\end{aligned}$
(8)

where $c$ is the speed of light, $\Phi(z,L)$ is the LF, $D_{L}(z)$ is the
luminosity distance, $z_{\text{min}}$ and $z_{\text{max}}$ represent the range
of integration in redshift and $\alpha$ is the spectral index.

In this work, we use the observed source counts from Algera et al. (2020) and
Hale et al. (2023) to provide an additional constraint in our analysis. The
Algera et al. (2020) $3$ GHz source counts are measured based on the ultra-
faint (reaching a $5\sigma$ flux limit of $\sim$2.7 $\mu$Jy beam-1 within the
center of the $3$ GHz image) radio population detected in the Karl G. Jansky
Very Large Array COSMOS-XS survey. The Hale et al. (2023) $1.4$ GHz source
counts are measured based on the continuum early science data release of the
MeerKAT International Gigahertz Tiered Extragalactic Exploration (MIGHTEE)
survey in the COSMOS and XMM-LSS fields. The MIGHTEE sources were divided into
three subsets: SFGs, AGNs and unclassified sources. Hale et al. (2023)
considered two cases: (1) the unclassified sources are assumed to be a mix of
SFGs and AGNs based on the flux density ratio of classified sources, and (2)
the unclassified sources are regarded as SFGs. The source counts for the two
cases are presented in their Table 1 (for the COSMOS field) and Table 2 (for
the XMM-LSS field). In this work, we use the data corresponding to the first
case, i.e., the $\mathrm{SC}_{\mathrm{SFG,\leavevmode\nobreak\ ratio}}$ column
in the Table 1 and Table 2 of Hale et al. (2023). For the convenience of
calculation, the above source counts are unified to $1.4$ GHz by assuming a
spectral index of 0.7.

###  3.3 Models for the SFG LFs

Without loss of generality, the SFG LF can be written as

$\displaystyle\Phi(z,L)=e_{1}(z)\phi(z=0,L/e_{2}(z),\eta^{j}),$ (9)

where $e_{1}(z)$ and $e_{2}(z)$ denote the density evolution (DE) and
luminosity evolution (LE) functions of redshift, respectively, and $\eta^{j}$
represents the parameters that determine the shape of the LF. If the values of
$\eta^{j}$ are constant, this indicates that the shape of the radio LF is
unchanged with redshift. Conversely, if $\eta^{j}$ exhibits a redshift
dependence, this implies luminosity-dependent density evolution (see Singal et
al., 2013, 2014, for a more detailed discussion). We assume the shape of the
LF to remain unchanged as many other studies (e.g., Novak et al., 2017; van
der Vlugt et al., 2022).

Following previous work (e.g. Smolčić et al., 2009; Gruppioni et al., 2013;
van der Vlugt et al., 2022), the SFG local LF $\phi(z=0,L/e_{2}(z=0))$ is
described by a modified-Schechter function from Saunders et al. (1990):

$\displaystyle\begin{aligned}
\phi(z&=0,L/e_{2}(z=0))=\frac{dN}{d\log_{10}L}\\\
&=\phi_{\star}\left(\frac{L}{L_{\star}}\right)^{1-\alpha}\exp\left[-\frac{1}{2\sigma^{2}}\log^{2}\left(1+\frac{L}{L_{\star}}\right)\right],\end{aligned}$
(10)

where $L_{\star}$ determines the knee location of the LF, $\alpha$ and
$\sigma$ fit the faint and bright ends of the LF, respectively, and
$\Phi_{\star}$ is used for the normalization. In this work, we consider three
LF models, all of which adopt the same LE function:

$\displaystyle e_{2}(z)=(1+z)^{k_{1}+k_{2}z}.$ (11)

The DE function $e_{1}(z)$ has three different forms depending upon the model:
$e_{1}(z)=1$ for model A, an exponential form

$\displaystyle e_{1}(z)=10^{p_{1}z},$ (12)

for model B, and

$\displaystyle e_{1}(z)=(1+z)^{p_{1}+p_{2}z},$ (13)

for model C. In above equations, $k_{1}$, $k_{2}$, $p_{1}$ and $p_{2}$ are
free parameters. Model A is the pure luminosity evolution (PLE) model, which
is the most commonly used model for the radio LFs of SFGs in the literature
(e.g., Novak et al., 2017). Models B and C can be referred to as the mixture
evolution (e.g., Yuan et al., 2016, 2017) or luminosity and density evolution
(LADE, e.g., Aird et al., 2010) models.

###  3.4 Kernel Density Estimation

In many previous studies, the non-parametric LFs are usually derived using the
traditional $1/V_{\rm max}$ method of Schmidt (1968). The mathematics behind
$1/V_{\rm max}$ is the histogram, the simplest density estimator used mainly
for rapid visualization of results in one or two dimensions. The $1/V_{\rm
max}$ method requires the binning of data, which undoubtedly leads to
information loss, and potential biases can be caused by evolution within the
bins (Yuan et al., 2020). Besides, its result is significantly dependent on
the choice of bin center and bin width but currently there are no effective
rules to realize an “optimal” binning (Yuan & Wang, 2013). To overcome issues
surrounding the binning of LFs, Yuan et al. (2020, 2022) recently proposed a
new method for estimating LFs in the framework of kernel density estimation
(KDE), the most popular nonparametric density estimation approach developed in
modern statistics. This method does not require the binning of data and any
model assumptions, and it simultaneously has some advantages of parametric and
nonparametric methods. In this work, we apply the KDE method to derive the
non-parametric LFs of SFGs. According to Yuan et al. (2020), the KDE LF can be
estimated by

$\displaystyle\begin{aligned}
\hat{\phi}(z,L)=\frac{n\left(Z_{2}-Z_{1}\right)\hat{f}\left(x,y\mid
h_{1},h_{2}\right)}{\left(z-Z_{1}\right)\left(Z_{2}-z\right)\Omega\frac{dV}{dz}}.\end{aligned}$
(14)

where [$Z_{1},Z_{2}$] is the redshift range of sample, $n$ is the number of
objects in the redshift range, $\Omega$ is the solid angle subtended by the
sample, $dV/dz$ is the differential comoving volume per unit of solid angle,
and $\hat{f}(x,y|h_{1},h_{2})$ is the density of $(x,y)$, which corresponds to
the $(z,L)$ pair in the KDE parameter space. The calculation of $\hat{f}$ can
be seen in Equation (9) in Yuan et al. (2022). Its value depends on two
bandwidth parameters $h_{1}$ and $h_{2}$. We use the Python package “kdeLF”, a
Bayesian MCMC routine provided by Yuan et al. (2022), to determine the
posterior distributions of the bandwidth parameters (Figure 3), and then the
uncertainty estimation on the LFs. The derived LF is shown in Figure 7 as red
solid lines with orange contours.

Figure 3:  Corner plot showing the posterior probability distributions of the
bandwidths ($h_{1}$, $h_{2}$) of the KDE LF, obtained with the routine
provided by Yuan et al. (2022). Uncertainties correspond to the 16th and 84th
percentiles, while the best-fit value is the median.  Figure 4:  Corner plot
showing the one and two dimensional projections of the posterior probability
distributions of the parameters for Model A, obtained from the MCMC run. The
histograms on the diagonal show the marginalised posterior densities for each
parameter (vertical dashed lines denote the 16th and 84th percentiles). The
off-diagonal panels show the 2D joint posterior densities of all couples of
parameters, with 1$\sigma$, 2$\sigma$, and 3$\sigma$ contours shown by black
solid lines. Our best fitting parameters are marked by red vertical dashed
lines. Figure 5: Similar to Figure 4, but for Model B. Figure 6: Similar to
Figure 4, but for Model C.

##  4 Results

Figure 7:  Radio LFs of SFGs at various redshifts compared with the previous
estimates specified in the inset. The best-fit LFs for Models A, B, and C in
each redshift bin are shown by the blue, orange and green solid lines,
respectively. The light green shaded area shows the 1$\sigma$ confidence
interval for Model C. The red solid curves represent the KDE LFs. The purple
dotted line depicts the PLE model of Novak et al. (2017). We also compare with
the binned LFs from Gruppioni et al. (2013), Novak et al. (2017) and van der
Vlugt et al. (2022).  Table 1: Best-fit Parameters of The Three Models

Model | $\text{log}_{10}(\phi_{\star})$ | $\text{log}_{10}(L_{\star})$ | $\alpha$ | $\sigma$ | $k_{1}$ | $k_{2}$ | $p_{1}$ | $p_{2}$
---|---|---|---|---|---|---|---|---
A | $-2.576_{-0.023}^{+0.022}$ | $21.232_{-0.041}^{+0.040}$ | $1.306_{-0.021}^{+0.021}$ | $0.611_{-0.007}^{+0.007}$ | $4.271_{-0.027}^{+0.028}$ | $-0.548_{-0.010}^{+0.010}$ | $\ldots$ | $\ldots$
B | $-2.554_{-0.019}^{+0.019}$ | $21.284_{-0.031}^{+0.031}$ | $1.296_{-0.020}^{+0.020}$ | $0.556_{-0.005}^{+0.005}$ | $4.662_{-0.026}^{+0.027}$ | $-0.473_{-0.011}^{+0.011}$ | $-0.188_{-0.008}^{+0.008}$ | $\ldots$
C | $-2.546_{-0.019}^{+0.019}$ | $21.248_{-0.032}^{+0.032}$ | $1.296_{-0.021}^{+0.020}$ | $0.562_{-0.005}^{+0.005}$ | $4.534_{-0.028}^{+0.030}$ | $-0.366_{-0.016}^{+0.017}$ | $-0.069_{-0.033}^{+0.016}$ | $-0.397_{-0.019}^{+0.022}$

Note. Units — $\phi_{\star}$: [${\rm Mpc^{-3}dex^{-1}}$], $L_{\star}$: [${\rm
WHz^{-1}}$]. The best-fit parameters as well as their 1$\sigma$ errors for
models A, B and C.

###  4.1 Analytical LFs

The parameters in our model LFs are estimated via the MCMC algorithm, which is
performed using the Python package emcee of (Foreman-Mackey et al., 2013). The
emcee algorithm improves the exploration of the parameter space by using an
ensemble of chains with different initial conditions. This approach helps to
avoid local minima and ensures a more comprehensive exploration of the
parameter space. We assume uniform priors on all the parameters. The
marginalized one and two dimensional posterior probability distributions of
the parameters for Models A, B, and C are shown in Figures 4, 5, 6,
respectively. These corner plots illustrate that all the parameters for our
three models are well constrained. Table 1 reports the best-fit parameters and
their 1$\sigma$ uncertainties for the three models.

Figure 7 shows our best-fit LFs for Models A (solid blue lines) , B (solid
orange lines), and C (solid green lines). All the LFs are measured at the
rest-frame 1.4 $\rm GHz$. We also compare our result with the binned LFs from
Gruppioni et al. (2013), Novak et al. (2017) and van der Vlugt et al. (2022),
which are represented by orange left triangles, dark blue, and sky blue
squares with error bars, respectively. At lower redshifts of $z<1.0$, our
three models can be hardly distinguishable. Above $z\sim 1.0$, model A
gradually separates from the other two models, and also significantly deviates
from the binned and KDE LFs. This indicates that model A is not applicable to
describe the evolution of LFs at higher redshift. By contrast, Models B and C
are in good agreement with the binned and KDE LFs for all redshift intervals.

In Figure 7, the purple dotted lines depict the PLE model of Novak et al.
(2017), which is generally in agreement with our Model A. Nevertheless, the
difference between the two results increases with redshift. This is not
surprising, since the PLE model of Novak et al. (2017) is constrained through
simultaneously fitting the LF points in all redshift bins, while these LF
points are estimated by the traditional $1/V_{\rm max}$ method. Obviously,
their modeling result would be dependent on the estimation accuracy of
$1/V_{\rm max}$. Our analytical LFs are obtained through a full maximum
likelihood analysis, thus being independent of the $1/V_{\rm max}$ estimates.

Figure 8: Comparison of our best-fit models with the Euclidean normalized
$1.4$ GHz source counts for SFGs observed in the literature. The orange,
green, and blue dotted dashed lines show our best-fit source counts of Models
A, B, and C, respectively. The source counts from Hale et al. (2023) in the
COSMOS and XMM-LSS fields are shown as purple and brown right triangles,
respectively. Also shown are the observed source counts from Smolčić et al.
(2017a) (green squares), Algera et al. (2020) (orange pentagons).

Figure 9: $L_{\star}\times e_{2}(z)$ (left panel) and $\phi_{\star}\times
e_{1}(z)$ (right panel) for our three LF models shown in different colours,
compared with the inference (in all fields) from Cochrane et al. (2023). The
light shaded areas take into account the 3$\sigma$ error bands. Our
$L_{\star}$ has been converted to $150$ MHz by assuming a spectral index of
0.7. The purple circles represent the variation of $L_{\star}$ (left panel)
and $\phi_{\star}$ (right panel) as functions of redshift inferred by Cochrane
et al. (2023).

###  4.2 Fitting the Observed source Counts

We calculate the source counts for each of our three model LFs using
Equation(8). The result is shown in Figure 8 as orange, green, and blue dotted
dashed lines for models A, B, and C, respectively. In the Figure, we compare
our models with the Euclidean normalized $1.4$ GHz source counts for SFGs
measured from Hale et al. (2023), Smolčić et al. (2017a) and Algera et al.
(2020). All our three models can reproduce the Algera et al. (2020) and Hale
et al. (2023) source counts fairly well. But the Smolčić et al. (2017a) result
is systematically lower than the others. The discrepancies could arise from
multiple reasons, such as field to field variation, differences in the
assumptions used to calculate completeness, spectral index assumptions
(converting the source counts from $3$ GHz to $1.4$ GHz) and resolution bias
(see Hale et al., 2023). Note that at the bright region ($S>1$ mJy), the
measurements of Hale et al. (2023) display huge uncertainties. This is due to
the contamination of AGNs at these fluxes, preventing an easy classification
between SFGs and AGNs.

###  4.3 LADE versus PLE

Due to the limitation of survey depth, most of the existing radio studies can
hardly reach the knee of the SFG LF at z ¿ 1. If one fits the LF using a LADE
model, the density and luminosity evolution parameters may become degenerate
(also see van der Vlugt et al., 2022). Therefore, previous studies usually
assumed a PLE model for their SFG LFs. Thanks to combining data set from the
ultra-deep COSMOS-XS survey and the shallower VLA-COSMOS 3 GHz large project,
van der Vlugt et al. (2022) was able to jointly constrain the luminosity and
density evolution, finding evidence for significant density evolution. From
Figure 7, we also find that the LADE model is superior to the PLE model.

Very recently, Cochrane et al. (2023) measured the 150 MHz LFs of SFGs
utilizing data from the Low Frequency Array (LOFAR) Two Metre Sky Survey in
three well-studied extragalactic fields, Elais-N1, Bo$\ddot{o}$es, and the
Lockman Hole. By fixing the faint and bright end shape of the radio LF to the
local values (equivalent to fixing $\sigma=0.49$ and $\alpha=1.12$ in Equation
10), they fitted the LF points (via $1/V_{\rm max}$) to find the best-fit
$L_{\star}$ and $\phi_{\star}$ at each individual redshift bin. In their
Figure 7, the variation of $L_{\star}$ and $\phi_{\star}$ as functions of
redshift was displayed. Cochrane et al. (2023)’s method is equivalent to
assuming a LADE model for their LFs. The variation of $L_{\star}$ and
$\phi_{\star}$ as functions of redshift in their analysis correspond to
$L_{\star}\times e_{2}(z)$ and $\phi_{\star}\times e_{1}(z)$ in our work. In
Figure 9, we show $L_{\star}\times e_{2}(z)$ and $\phi_{\star}\times e_{1}(z)$
for our three LF models, compared with the inference from Cochrane et al.
(2023). Below $z<2$, our two LADE models, especially Model C, agree well with
the $L_{\star}$ and $\phi_{\star}$ evolutions given by Cochrane et al. (2023).
Above $z>2$, their $L_{\star}$ evolution curve is significantly higher than
those of our models, while their $\phi_{\star}$ evolution falls more rapidly
than those of our models out to high redshift.

The discrepancies could be explained as follows: At any redshift bin, the
best-fit $L_{\star}$ and $\phi_{\star}$ are negatively correlated with each
other (e.g., see the right panel of Figure 4 in Cochrane et al. (2023)),
implying that the luminosity and density evolution parameters are degenerate.
The degeneracy should be stronger at higher redshift where the knee location
of LFs is increasingly difficult to identified. Recall that Cochrane et al.
(2023) fitted the LF points (via $1/V_{\rm max}$) for each redshift bin
individually. At higher redshift, the LF points estimated by $1/V_{\rm max}$
usually have larger uncertainty. Fitting these discrete LF points to find a
precise knee location would be very difficult. This will inevitably bias the
inferred $L_{\star}$ and $\phi_{\star}$ values. Therefore, their $L_{\star}$
and $\phi_{\star}$ evolutions for high redshift are subject to uncertainties
due to dual factors. Our LADE modeling is also subject to the uncertainty due
to degeneracy, but it is free from the $1/V_{\rm max}$ estimates.

Cochrane et al. (2023) found that their $\phi_{\star}$ remains roughly
constant back to $z\sim 0.8$ but then falls steeply at higher redshifts. Our
Model C displays a similar trend where the density evolution falls more
rapidly at higher redshifts. A comparison between Cochrane et al. (2023)’s
inference with our models are meaningful, implying that Model C could be taken
as our prefered model.

Figure 10: Cosmic star formation rate density (SFRD) history. Our SFRD history
results are shown as blue, red, and green solid lines for Models A, B, and C,
respectively. The light shaded areas take into account the $3\sigma$ error
bands. The compilation of Madau & Dickinson (2014) is shown as a black dashed
line in all panels. All data shown for comparison are indicated in the legend
of each panel; see text for details.

##  5 Cosmic star formation rate density history

###  5.1 Calculating the SFRD

Now that we have obtained the rest-frame $\mathrm{1.4\leavevmode\nobreak\
GHz}$ LF, we can investigate how the SFRD evolves with redshift. We use the
functional form provided in Delvecchio et al. (2021) to convert the radio
luminosity into an SFR:

$\displaystyle\frac{\mathrm{SFR}}{M_{\odot}\mathrm{yr}^{-1}}=f_{\mathrm{IMF}}\times
10^{-24}10^{q_{\mathrm{IR}}(z)}\frac{L_{1.4\mathrm{GHz}}}{\mathrm{WHz}^{-1}}$
(15)

where $f_{\text{IMF}}$ is a factor accounting for the IMF ( $f_{\text{IMF}}=1$
for a Chabrier (2003) IMF, and $f_{\text{IMF}}=1.7$ for a Salpeter (1955)
IMF), and $L_{1.4\mathrm{GHz}}$ is rest-frame $1.4\,\mathrm{GHz}$ radio
luminosity. Following Novak et al. (2017), we use the Chabrier IMF in the
following analysis. In Equation (15), $q_{\mathrm{IR}}(z)$ is the FIR-to-radio
luminosity ratio which is conventionally used to parametrize the FIR-radio
correlation in SFGs, defined as

$q_{\mathrm{IR}}=\log\left(\frac{L_{\mathrm{IR}}\left[\mathrm{W}\right]}{3.75\times
10^{12}\left[\mathrm{Hz}\right]}\right)-\log(L_{1.4\mathrm{GHz}}[\mathrm{W\,Hz^{-1}}]),$

where $L_{\mathrm{IR}}$ is the total IR luminosity (rest-frame 8-1000
$\mu\text{m}$), and $3.75\times 10^{12}\mathrm{Hz}$ represents the central
frequency over the far-infrared domain.

Although $q_{\mathrm{IR}}$ is typically taken to be a constant value derived
for local galaxies, recent observations suggest that the $q_{\mathrm{IR}}$
value probably changes with redshift (e.g. Sargent et al., 2010; Magnelli et
al., 2015; Delhaize et al., 2017; Calistro Rivera et al., 2017). Recently, for
the first time, Delvecchio et al. (2021) calibrated $q_{\mathrm{IR}}$ as a
function of both stellar mass ($M_{\star}$) and redshift. They found that
$q_{\mathrm{IR}}$ primarily evolves with $M_{\star}$, and only weakly with
redshift. This finding implies that using radio emission as an SFR tracer
requires $M_{\star}$-dependent conversion factors, but its robustness still
needs to be verified by more following studies. In this work, we use the
expression given by Novak et al. (2017):

$\displaystyle q_{\mathrm{IR}}(z)=(2.78\pm 0.02)\times(1+z)^{-0.14\pm 0.01},$
(16)

which is the updated version of Delhaize et al. (2017)’s analysis by using a
new sample selection criteria for excluding AGN.

The SFRD of a given epoch can then be estimated by the following integral:

$\displaystyle\text{SFRD}=\int_{L_{\text{min}}}^{L_{\text{max}}}\Phi(L,z)\times\text{SFR}(L_{1.4\,\text{GHz}})\,\mathrm{d}\log_{10}L\,.$
(17)

To obtain the SFRD for a given epoch, we performed a numerical integration of
the analytical form of the LF in each redshift bin, employing the best-fit
evolution parameters presented in Table 1 and Figure 7. The resulting integral
provides an estimate of the SFRD. Unless otherwise specified, our reported
SFRD values correspond to the integral of the fitted LF from 0.0 to
$\to\infty$. In Figure 10(a), our SFRD history results are shown as blue, red,
and green solid lines for Models A, B, and C, respectively. The light shaded
areas take into account the $3\sigma$ error bands. Our three models coincide
at lower redshift of $z<1$. After then, model A gradually separates from the
other two models. The difference of Models B and C occurs at $\sim z>2$, but
the significance is less than 3 $\sigma$ confidence level.

###  5.2 Comparison with the Literature

In Figure 10, we show the SFRD history derived from our three models compared
with those in the literature derived at different wavelengths. The SFRD from
the review by Madau & Dickinson (2014), where the fit was performed on a
collection of previously published UV and IR data, is shown as the black
dashed curve in all panels for ease of comparison. From Figure 10(a), we find
that below $z\sim 2$, the SFRD from our three models match well with that of
Madau & Dickinson (2014), with the rising shape of their SFRD curve being
almost the same as ours. Above $z\sim 2$, our PLE model predicts a
significantly higher SFRD, while the SFRD given by our LADE models is lower
than that of Madau & Dickinson (2014). The discrepancy in magnitude is
probably due to the assumed FIR–radio relation, which has a significant impact
on our derived SFRD (see van der Vlugt et al., 2022). The SFRD derived from
all our three models turns over at a slightly higher redshift ($2<z<2.5$) and
falls more rapidly than that of Madau & Dickinson (2014) out to high redshift.
This has also been observed by recent radio estimates from van der Vlugt et
al. (2022).

In Figure 10(b), we show our Model A SFRD compared to the radio estimates from
Smolčić et al. (2009) and Novak et al. (2017). Smolčić et al. (2009) derived
the SFRD out to $z=1.3$ by assuming a PLE LF and a nonevolving FIR–radio
correlation established by Bell (2003). We find a good agreement with the
Smolčić et al. (2009) estimates, and therefore their result provides a good
consistency check of our models at low redshift. The result of Novak et al.
(2017), who also assumed a PLE LF, is the key comparative object for our model
A, since our analysis is based on the same sample with them. Unsurprisingly,
the curve of our Model A SFRD seems to be a good fit to their SFRD points. The
only discrepancy occurs at the data point with the highest redshift, their
SFRD being significantly higher than the prediction of our Model A. Because
Novak et al. (2017) and we use the same $q_{\mathrm{IR}}(z)$ evolution to
calculate the SFRD, any discrepancy between the two results can only arise
from the difference in LFs (see Figure 7). As discussed in section 4.1, the
analytical LFs of Novak et al. (2017) are obtained by fitting the $1/V_{\rm
max}$ LF points in all redshift bins, while our LFs are obtained through a
full maximum likelihood analysis and also utilizing the additional constraints
from source counts and local LFs. Therefore our LFs should be more accurate
than that of Novak et al. (2017), making our SFRD to be an improvement to
their estimates.

In Figure 10(c), we show the SFRD derived from our Models B and C compared to
the radio estimates from van der Vlugt et al. (2022) and Smolčić et al.
(2009). van der Vlugt et al. (2022) also assumed a LADE LF, similar to our
Models B and C, to calculate their SFRD based on the combined data set from
the ultra-deep COSMOS-XS survey and the VLA-COSMOS 3 GHz large project.
Although adopting a different $q_{\mathrm{IR}}(z)$ evolution with us, their
result is in good agreement with those of our Models B and C, especially our
Model C, within the error bars. Note that three points of their estimates show
a slightly elevated SFRD compared to our model prediction. The discrepancies
could be attributed to the different $q_{\mathrm{IR}}(z)$ evolution they used,
or the uncertainties in their LF measurement. Similar to Novak et al. (2017),
van der Vlugt et al. (2022) also obtained their analytical LFs by fitting the
$1/V_{\rm max}$ LF points in all redshift bins.

In Figure 10(d), we show the SFRD derived from our Model B compared to the
radio estimates from Karim et al. (2011) and Cochrane et al. (2023). Our
estimates are slightly lower than that of Karim et al. (2011), with the
difference increasing with redshift. The discrepancies could be attributed to
the different approach, where they performed stacking on mass selected
galaxies and used a nonevolving FIR–radio correlation established by Bell
(2003). The measurements of Cochrane et al. (2023) are systematically higher
than ours. If we multiply our SFRD by 1.5, the new curve (red dashed) will
match the Cochrane et al. (2023) data over the whole redshift range. Different
from the $q_{\mathrm{IR}}(z)$-based $L_{\mathrm{radio}}-\mathrm{SFR}$
calibration used in this work, Cochrane et al. (2023) used the calibrated
relation between 150 MHz radio luminosity and SFR from Smith et al. (2021),
and also constrained and corrected the scatter in the
$L_{150\leavevmode\nobreak\ \mathrm{MHz}}-\mathrm{SFR}$ relation. This can
explain the discrepancy between their measurements and ours. As noted by
Leslie et al. (2020), the impact of different
$L_{\mathrm{radio}}-\mathrm{SFR}$ calibrations is significant.

##  6 Summary and Conclusions

In this work, we make use of an SFG sample (5915 sources) from the VLA-COSMOS
3 GHz data (Smolčić et al., 2017a) to measure the radio LFs. For the first
time, we use both models of PLE and joint luminosity+density evolution (LADE)
to fit the LFs directly to the radio data using a full maximum likelihood
analysis. We have fully considered the effect of completeness correction for
the sample. We also incorporate the updated observations on local radio LFs
and source counts into the fitting process to obtain additional constraints.
The parameters of fitting are determined through the Bayesian Markov Chain
Monte Carlo (MCMC) approach. In addition, to provide an alternative
nonparametric LF estimate for SFGs, we have applied the kernel density
estimation (KDE) method described in our previous work (Yuan et al., 2020,
2022). Based on these radio LFs, we have derived the dust-unbiased SFRD up to
redshift of $z\sim 5$. The main results of our work are as follows.

  1. 1.

Our radio LFs are fitted using three models, assuming a modified-Schechter
function with PLE (Model A) and LADE (Models B and C), respectively. Below
$z<2$, the PLE model can fit the radio LFs well, while at higher redshift it
gradually deviates from the observations, indicating that PLE is not
applicable to describe the evolution at higher redshift. By contrast, our LADE
models can successfully fit these observed radio LFs from previous studies
(Gruppioni et al., 2013; Novak et al., 2017; van der Vlugt et al., 2022) over
the whole redshift range. All our three models can reproduce the latest radio
source counts from Algera et al. (2020) and Hale et al. (2023) fairly well.

  2. 2.

We find that the luminosity and density evolutions of our LADE models are
broadly consistent with the $L_{\star}$ and $\phi_{\star}$ evolutions inferred
by Cochrane et al. (2023), within the uncertainty due to degeneracy of
luminosity and density evolution parameters. In view of that Cochrane et al.
(2023)’s inference is model independent, thus the consistency with their
result is a strong support to our LADE models, and we conclude that a density
evolution is genuinely indispensable in modeling the evolution of SFG radio
LFs.

  3. 3.

The SFRD curves derived from our PLE and LADE models show a good fit to the
SFRD points derived by Novak et al. (2017) and van der Vlugt et al. (2022),
respectively. Both their analytical LFs were obtained by fitting the $1/V_{\rm
max}$ LF points in all redshift bins, possibly suffering from extra errors
propagated from the $1/V_{\rm max}$ estimate itself. Our analytical LFs are
independent of the $1/V_{\rm max}$ estimates and should be more accurate,
making our SFRD results to be an improvement to these previous estimates.

  4. 4.

Below $z\sim 2$, the SFRD from both our PLE and LADE models matchs well with
the multi-wavelength compilation from Madau & Dickinson (2014), with the
rising shape of their SFRD curve being almost the same as ours. Comparing with
Madau & Dickinson (2014), our SFRD predicted by all three models turns over at
a slightly higher redshift ($2<z<2.5$) and falls more rapidly out to high
redshift. This has also been observed by recent radio estimates from van der
Vlugt et al. (2022).

  5. 5.

The very recent measurements from Cochrane et al. (2023) are systematically
higher than ours (by up to 0.4 dex). The discrepancy coud be due to the
different $L_{\text{radio}}-\text{SFR}$ calibration and its scatter correction
in their work. From the comparison with Cochrane et al. (2023), we highlight
that the $L_{\text{radio}}-\text{SFR}$ calibration and its scatter has a
significant impact and remains one of the biggest uncertainties in the radio
estimates of SFRD.

###### Acknowledgements.

This work is supported by the National Key R&D Program of China
(2023YFE0101200). We acknowledge the financial support from the National
Natural Science Foundation of China (grants No. 12073069 and No. 12075084).
Z.Y. is supported by the Xiaoxiang Scholars Programme of Hunan Normal
University. J.M is supported by the National Science Foundation of China (NSFC
11673062), the China Manned Space Project (CMS-CSST2021-A06), and the Yunnan
Revitalization Talent Support Program (YunLing Scholar Award). We would like
to thank Yaqian Yu and Yang Liu for helpful discussions.

## References

  * Adams et al. (2020) Adams, N. J., Bowler, R. A. A., Jarvis, M. J., et al. 2020, MNRAS, 494, 1771
  * Aird et al. (2010) Aird, J., Nandra, K., Laird, E. S., et al. 2010, MNRAS, 401, 2531
  * Algera et al. (2020) Algera, H. S. B., van der Vlugt, D., Hodge, J. A., et al. 2020, ApJ, 903, 139
  * Bell (2003) Bell, E. F. 2003, ApJ, 586, 794
  * Best et al. (2005) Best, P., Kauffmann, G., Heckman, T., & Ivezić, Ž. 2005, Monthly Notices of the Royal Astronomical Society, 362, 9
  * Bonato et al. (2017) Bonato, M., Negrello, M., Mancuso, C., et al. 2017, MNRAS, 469, 1912
  * Bonato et al. (2021a) Bonato, M., Prandoni, I., De Zotti, G., et al. 2021a, A&A, 656, A48
  * Bonato et al. (2021b) Bonato, M., Prandoni, I., De Zotti, G., et al. 2021b, MNRAS, 500, 22
  * Bouwens et al. (2009) Bouwens, R. J., Illingworth, G. D., Franx, M., et al. 2009, ApJ, 705, 936
  * Bouwens et al. (2015) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34
  * Bouwens et al. (2021) Bouwens, R. J., Oesch, P. A., Stefanon, M., et al. 2021, AJ, 162, 47
  * Bowler et al. (2015) Bowler, R. A. A., Dunlop, J. S., McLure, R. J., et al. 2015, MNRAS, 452, 1817
  * Bowler et al. (2020) Bowler, R. A. A., Jarvis, M. J., Dunlop, J. S., et al. 2020, MNRAS, 493, 2059
  * Calistro Rivera et al. (2017) Calistro Rivera, G., Williams, W. L., Hardcastle, M. J., et al. 2017, MNRAS, 469, 3468
  * Ceraj et al. (2018) Ceraj, L., Smolčić, V., Delvecchio, I., et al. 2018, A&A, 620, A192
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Clemens et al. (2008) Clemens, M. S., Vega, O., Bressan, A., et al. 2008, A&A, 477, 95
  * Cochrane et al. (2023) Cochrane, R. K., Kondapally, R., Best, P. N., et al. 2023, MNRAS, 523, 6082
  * Condon (1992) Condon, J. J. 1992, ARA&A, 30, 575
  * Condon et al. (2002) Condon, J. J., Cotton, W. D., & Broderick, J. J. 2002, AJ, 124, 675
  * Condon et al. (2019) Condon, J. J., Matthews, A. M., & Broderick, J. J. 2019, ApJ, 872, 148
  * de Zotti et al. (2010) de Zotti, G., Massardi, M., Negrello, M., & Wall, J. 2010, A&A Rev., 18, 1
  * Delhaize et al. (2017) Delhaize, J., Smolčić, V., Delvecchio, I., et al. 2017, A&A, 602, A4
  * Delvecchio et al. (2021) Delvecchio, I., Daddi, E., Sargent, M. T., et al. 2021, A&A, 647, A123
  * Delvecchio et al. (2017) Delvecchio, I., Smolčić, V., Zamorani, G., et al. 2017, A&A, 602, A3
  * Dudzevičiūtė et al. (2020) Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2020, MNRAS, 494, 3828
  * Enia et al. (2022) Enia, A., Talia, M., Pozzi, F., et al. 2022, ApJ, 927, 204
  * Fan et al. (2001) Fan, X., Strauss, M. A., Schneider, D. P., et al. 2001, The Astronomical Journal, 121, 54
  * Finkelstein et al. (2015) Finkelstein, S. L., Ryan, Russell E., J., Papovich, C., et al. 2015, ApJ, 810, 71
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306
  * Gruppioni et al. (2020) Gruppioni, C., Béthermin, M., Loiacono, F., et al. 2020, A&A, 643, A8
  * Gruppioni et al. (2013) Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23
  * Hale et al. (2023) Hale, C. L., Whittam, I. H., Jarvis, M. J., et al. 2023, MNRAS, 520, 2668
  * Helou et al. (1985) Helou, G., Soifer, B. T., & Rowan-Robinson, M. 1985, ApJ, 298, L7
  * Hogg (1999) Hogg, D. W. 1999, arXiv preprint astro-ph/9905116
  * Ishigaki et al. (2018) Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2018, ApJ, 854, 73
  * Jarvis et al. (2015) Jarvis, M., Seymour, N., Afonso, J., et al. 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14), 68
  * Jarvis et al. (2010) Jarvis, M. J., Smith, D. J. B., Bonfield, D. G., et al. 2010, MNRAS, 409, 92
  * Karim et al. (2011) Karim, A., Schinnerer, E., Martínez-Sansigre, A., et al. 2011, The Astrophysical Journal, 730, 61
  * Kennicutt Jr (1998) Kennicutt Jr, R. C. 1998, Annual Review of Astronomy and Astrophysics, 36, 189
  * Kochanek (1996) Kochanek, C. S. 1996, ApJ, 473, 595
  * Lara-López et al. (2010) Lara-López, M., Cepa, J., Bongiovanni, A., et al. 2010, Astronomy & Astrophysics, 521, L53
  * Leslie et al. (2020) Leslie, S. K., Schinnerer, E., Liu, D., et al. 2020, ApJ, 899, 58
  * Lewis & Bridle (2002) Lewis, A. & Bridle, S. 2002, Physical Review D, 66, 103511
  * Madau & Dickinson (2014) Madau, P. & Dickinson, M. 2014, ARA&A, 52, 415
  * Magnelli et al. (2015) Magnelli, B., Ivison, R., Lutz, D., et al. 2015, Astronomy & Astrophysics, 573, A45
  * Malefahlo et al. (2022) Malefahlo, E. D., Jarvis, M. J., Santos, M. G., et al. 2022, MNRAS, 509, 4291
  * Mandal et al. (2021) Mandal, S., Prandoni, I., Hardcastle, M. J., et al. 2021, A&A, 648, A5
  * Marshall et al. (1983) Marshall, H., Tananbaum, H., Avni, Y., & Zamorani, G. 1983, The Astrophysical Journal, 269, 35
  * Matthews et al. (2021a) Matthews, A. M., Condon, J. J., Cotton, W. D., & Mauch, T. 2021a, ApJ, 914, 126
  * Matthews et al. (2021b) Matthews, A. M., Condon, J. J., Cotton, W. D., & Mauch, T. 2021b, ApJ, 909, 193
  * Mauch & Sadler (2007) Mauch, T. & Sadler, E. M. 2007, MNRAS, 375, 931
  * McLeod et al. (2016) McLeod, D. J., McLure, R. J., & Dunlop, J. S. 2016, MNRAS, 459, 3812
  * McLeod et al. (2015) McLeod, D. J., McLure, R. J., Dunlop, J. S., et al. 2015, MNRAS, 450, 3032
  * McLure et al. (2013) McLure, R. J., Dunlop, J. S., Bowler, R. A. A., et al. 2013, MNRAS, 432, 2696
  * Novak et al. (2017) Novak, M., Smolčić, V., Delhaize, J., et al. 2017, A&A, 602, A5
  * Ocran et al. (2020a) Ocran, E. F., Taylor, A. R., Vaccari, M., Ishwara-Chandra, C. H., & Prandoni, I. 2020a, MNRAS, 491, 1127
  * Ocran et al. (2020b) Ocran, E. F., Taylor, A. R., Vaccari, M., et al. 2020b, MNRAS, 491, 5911
  * Oesch et al. (2018) Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labbé, I., & Stefanon, M. 2018, ApJ, 855, 105
  * Ono et al. (2018) Ono, Y., Ouchi, M., Harikane, Y., et al. 2018, PASJ, 70, S10
  * Padovani (2016) Padovani, P. 2016, Astronomy and Astrophysics Review, 24, 13
  * Parsa et al. (2016) Parsa, S., Dunlop, J. S., McLure, R. J., & Mortlock, A. 2016, MNRAS, 456, 3194
  * Riechers et al. (2013) Riechers, D. A., Bradford, C. M., Clements, D. L., et al. 2013, Nature, 496, 329
  * Rowan-Robinson et al. (2016) Rowan-Robinson, M., Oliver, S., Wang, L., et al. 2016, MNRAS, 461, 1100
  * Sadler et al. (1989) Sadler, E. M., Jenkins, C. R., & Kotanyi, C. G. 1989, MNRAS, 240, 591
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161
  * Sargent et al. (2010) Sargent, M. T., Schinnerer, E., Murphy, E., et al. 2010, ApJS, 186, 341
  * Saunders et al. (1990) Saunders, W., Rowan-Robinson, M., Lawrence, A., et al. 1990, Monthly Notices of the Royal Astronomical Society, 242, 318
  * Schmidt (1968) Schmidt, M. 1968, ApJ, 151, 393
  * Singal et al. (2014) Singal, J., Ko, A., & Petrosian, V. 2014, The Astrophysical Journal, 786, 109
  * Singal et al. (2013) Singal, J., Petrosian, V., Lawrence, A., et al. 2013, The Astrophysical Journal, 764, 43
  * Smail et al. (1997) Smail, I., Ivison, R. J., & Blain, A. W. 1997, ApJ, 490, L5
  * Smith et al. (2021) Smith, D. J. B., Haskell, P., Gürkan, G., et al. 2021, A&A, 648, A6
  * Smolčić et al. (2009) Smolčić, V., Zamorani, G., Schinnerer, E., et al. 2009, The Astrophysical Journal, 696, 24
  * Smolčić et al. (2017a) Smolčić, V., Novak, M., Bondi, M., et al. 2017a, A&A, 602, A1
  * Smolčić et al. (2017b) Smolčić, V., Delvecchio, I., Zamorani, G., et al. 2017b, A&A, 602, A2
  * Smolčić et al. (2009) Smolčić, V., Schinnerer, E., Zamorani, G., et al. 2009, ApJ, 690, 610
  * Tabatabaei et al. (2017) Tabatabaei, F. S., Schinnerer, E., Krause, M., et al. 2017, ApJ, 836, 185
  * Upjohn et al. (2019) Upjohn, J. E., Brown, M. J. I., Hopkins, A. M., & Bonne, N. J. 2019, PASA, 36, e012
  * van der Vlugt et al. (2021) van der Vlugt, D., Algera, H. S. B., Hodge, J. A., et al. 2021, ApJ, 907, 5
  * van der Vlugt et al. (2022) van der Vlugt, D., Hodge, J. A., Algera, H. S. B., et al. 2022, ApJ, 941, 10
  * Willott et al. (2001) Willott, C. J., Rawlings, S., Blundell, K. M., Lacy, M., & Eales, S. A. 2001, Monthly Notices of the Royal Astronomical Society, 322, 536
  * Yuan et al. (2020) Yuan, Z., Jarvis, M. J., & Wang, J. 2020, The Astrophysical Journal Supplement Series, 248, 1
  * Yuan & Wang (2013) Yuan, Z. & Wang, J. 2013, Ap&SS, 345, 305
  * Yuan et al. (2016) Yuan, Z., Wang, J., Zhou, M., & Mao, J. 2016, The Astrophysical Journal, 820, 65
  * Yuan et al. (2017) Yuan, Z., Wang, J., Zhou, M., Qin, L., & Mao, J. 2017, The Astrophysical Journal, 846, 78
  * Yuan et al. (2022) Yuan, Z., Zhang, X., Wang, J., Cheng, X., & Wang, W. 2022, The Astrophysical Journal Supplement Series, 260, 10
  * Yun et al. (2001) Yun, M. S., Reddy, N. A., & Condon, J. J. 2001, ApJ, 554, 803