Issue 
A&A
Volume 683, March 2024



Article Number  A174  
Number of page(s)  14  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/202347746  
Published online  18 March 2024 
Constraints on cosmic star formation history via a new modeling of the radio luminosity function of starforming galaxies
^{1}
Department of Physics, School of Physics and Electronics, Hunan Normal University, Changsha 410081, PR China
email: yzl@hunnu.edu.cn, hwyu@hunnu.edu.cn
^{2}
Key Laboratory of Low Dimensional Quantum Structures and Quantum Control, Hunan Normal University, Changsha 410081, PR China
^{3}
Synergetic Innovation Center for Quantum Effects and Applications, and Institute of Interdisciplinary Studies, Hunan Normal University, Changsha, Hunan 410081, PR China
^{4}
Yunnan Observatories, Chinese Academy of Sciences, Kunming 650216, PR China
^{5}
Center for Astronomical MegaScience, Chinese Academy of Sciences, Beijing 100012, PR China
^{6}
Key Laboratory for the Structure and Evolution of Celestial Objects, Chinese Academy of Sciences, Kunming 650216, PR China
Received:
17
August
2023
Accepted:
12
January
2024
Context. Radio wavelengths offer a unique possibility to trace the total starformation rate (SFR) in galaxies, both obscured and unobscured. To probe the dustunbiased starformation history, an accurate measurement of the radio luminosity function (LF) for starforming galaxies (SFGs) is crucial.
Aims. We make use of an SFG sample (5900 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 history of the cosmic SFR density (SFRD) from z ∼ 5 onwards.
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 maximumlikelihood analysis, considering the sample completeness correction. We also incorporate updated observations of local radio LFs and radio source counts into the fitting process to obtain additional constraints.
Results. We find that the PLE model cannot be used 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. The Akaike information criterion (AIC) also demonstrates that the LADE model is superior to the PLE model. We therefore conclude that density evolution is genuinely indispensable in modeling the evolution of SFG radio LFs. Our SFRD curve shows a good fit to the SFRD points derived by previous radio estimates. In view of the fact that our radio LFs are not biased, as opposed those of previous studies performed by fitting the 1/V_{max} LF points, our SFRD results should be an improvement on these previous estimates. Below z ∼ 1.5, our SFRD matches a published multiwavelength compilation, while our SFRD turns over at a slightly higher redshift (2 < z < 2.5) and falls more rapidly out to high redshift.
Key words: Galaxy: evolution / galaxies: luminosity function / mass function / galaxies: star formation / radio continuum: galaxies
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
Understanding the formation and evolution of galaxies through 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 therefore a critical parameter for galaxies (e.g., LaraLópez et al. 2010). In recent decades, a lot of research based on multiwavelength aspects has been devoted to the accurate measurement of the cosmic evolution of the SFRD. Remarkable progress has been made, and the SFRD is now well understood up to z ∼ 3, when the Universe was no more than 3 Gyr old (see Madau & Dickinson 2014 for an exhaustive review). A consensus has been reached regarding recent history, where the SFRD increases significantly with redshift, reaching a peak at redshift z ∼ 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 UVselected 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, 2020; RowanRobinson et al. 2016; Novak et al. 2017). Such uncertainty in the evolution of the SFRD at early cosmic epochs hinders our understanding of the core mechanism that governs the 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 1998). In dustfree 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, 2020; Finkelstein et al. 2015; McLeod et al. 2015; Parsa et al. 2016; Oesch et al. 2018; Ono et al. 2018; Adams et al. 2020; Bouwens et al. 2021). However, UV observations suffer from dust absorption, which means the SFR measurements made at these wavelengths are 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 farinfrared (FIR) wavelengths. Therefore, FIR emission is ideal for tracing SFR in dustrich environments (see Kennicutt 1998). Unfortunately, FIR observations can suffer from poor resolution and source blending.
Deep radio continuum observations are now believed to be very promising tracers; they offer a unique possibility to trace the total SFR in galaxies, both obscured and unobscured. As such, they may provide the most robust measurement of the starformation history of the Universe (Jarvis et al. 2015). Radio continuum emission, which is 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 shortlived (τ ≤ 3 × 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 starforming galaxies (SFGs), known as the FIR–radio correlation. This correlation holds over five orders of magnitude in luminosity 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 submilliJansky (mJy) detection limits have emerged 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. 2017, 2021a,b; Ocran et al. 2020a; Upjohn et al. 2019; Ceraj et al. 2018; Novak et al. 2017; Smolčić et al. 2009b). These studies generally measured the radio luminosity functions (LFs) of SFGs; 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. 2009a; Novak et al. 2017; Ocran et al. 2020a; Malefahlo et al. 2022). Very recently, van der Vlugt et al. (2022) combined the COSMOSXS survey and Very Large Array (VLA)COSMOS 3 GHz data sets to constrain a radio LF with both luminosity and density evolution. The analytical LFs from these studies are obtained through fitting the LF points given by the Schmidt (1968) 1/V_{max} estimator. This semiparametric method has also been adopted by almost all the existing studies (e.g., Smolčić et al. 2009b; Novak et al. 2017; Cochrane et al. 2023). However, given the ordinary precision in the 1/V_{max} estimate, the LF points themselves have errors, and fitting to them will propagate the uncertainties to the analytical LFs. In addition, the result would be dependent on the choice of binning in the 1/V_{max} method (see Fan et al. 2001). We believe that a more reliable approach to obtain the analytical LFs is to use a full maximumlikelihood analysis (e.g., Willott et al. 2001).
In the present paper, we make use of the VLACOSMOS 3 GHz data (Smolčić et al. 2017a) to measure the radio LFs of SFGs. We use both models of PLE and joint density+luminosity evolution to fit the SFG LFs directly to the radio data using a full maximumlikelihood analysis. We aim to perform a comprehensive parametric study of the radio LF of SFGs by means of constraints from multiple observational data. Finally, we can probe the dustunbiased SFRD up to a redshift of z ∼ 5.
The structure of the present paper is outlined below. In Sect. 2, we briefly describe the data used. In Sect. 3, we present the method used to constrain the LFs with redshift. In Sect. 4, we derive our radio LF evolution through cosmic time and compare it to those in the literature. In Sect. 5, we calculate the evolution of the cosmic SFRD using the LF models we derived and compare it to the literature. In Sect. 6, we summarize our findings and conclusions.
Throughout the paper, we use the flat concordance Λ cold dark matter (ΛCDM) cosmology with the following parameters: Hubble constant H_{0} = 70 km s^{−1} Mpc^{−1}, dark energy density Ω_{Λ} = 0.7, and matter density Ω_{m} = 0.3. We assume the Chabrier (2003) initial mass function (IMF) to calculate SFRs. We assume a simple powerlaw radio spectrum for SFGs, F_{ν} ∝ ν^{−α}, where F_{ν} is the flux density at frequency ν and α is the spectral index.
2. Sample
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 VLACOSMOS 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 multiband 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 VLACOSMOS 3 GHz Large Project survey utilized 384 h of VLA A+C array observations in the S band to obtain radio data. The survey covered a uniform rms noise of 2.3 μJy beam^{−1} and had an angular resolution of 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 radioexcess if the value of r deviates by more than 3σ from the peak of the distribution obtained as a function of redshift, that is,
According to this criterion, they were able to distinguish 1814 sources (23%) that are primarily emitting due to AGN activity in the radio. The sample consists of 5915 SFGs that do not exhibit radio excess. All the sources have a (spectroscopic or photometric) redshift and the restframe 1.4 GHz luminosity. The redshift distribution of the SFG sample as well as its scatter plot are shown in Fig. 1. The red dashed curve indicates the 1.4 GHz flux limit line defined as
Fig. 1. Redshift distribution (top) and the scatter plot (bottom) of our SFG sample. The red dashed curve indicates the flux limit line f_{lim 1.4 GHz}(z). 
where D_{L} represents the luminosity distance at redshift z, F_{lim 3 GHz} = 11.5 μJy is the 5σ detection limit of the survey at 3 GHz, and the spectral index α is set to 0.7. We have excluded all sources below the flux limit line, and the total number of sources used in this work is 5900.
3. Methods
3.1. Luminosity function and likelihood function
The LF Φ(z, L) is a measurement of the number of sources per unit comoving volume per unit logarithmic luminosity interval:
Given an analytical form with parameters θ for the LF, Φ(z, Lθ), the maximumlikelihood solution to θ is obtained by minimizing the negative logarithmic likelihood function S. Following Marshall et al. (1983) and (Fan et al. 2001), S can be written as
where p(z, L) is the selection probability of the SFG as a function of redshift and luminosity, and W is the survey region. The inclusion of the selection probability in Eq. (4) accounts for the fact that the sample is incomplete near the flux limit. The symbol Ω 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
where C_{radio} is the completeness of the VLACOSMOS 3 GHz radio catalog as a function of the flux density F_{3 GHz}, and C_{opt} is the completeness owing to radio sources without assigned opticalNIR counterparts (Novak et al. 2017). We adopt the calculations of C_{radio} and C_{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 Eq. (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 twodimensional (2D) grid of 50 × 50 in the log L − z space. For each grid point (log L_{i}, z_{i}), we can derive its flux density F_{i} from L_{i} by assuming α = 0.7. We can then estimate the corresponding C_{radio} and C_{opt} through a onedimensional linear interpolation method using the observed value from Novak et al. (2017). Finally, we have the values for p(z, L) at the 50 × 50 grid points, which are used to perform the 2D linear interpolation to estimate the function value of p(z, L).
Following the method of Willott et al. (2001) and Yuan et al. (2017), we incorporate the most recent observations of the local radio LFs and source counts (see Sect. 3.2) into the fitting process to obtain additional constraints. The local radio LF (LRLF) and the source counts (SCs) are onedimensional functions, and their χ^{2} value is calculated as
where f_{data i} represents the value of the data in the ith bin, and f_{mod i} and σ_{data i} are the model value and data error in the ith bin, respectively. As χ^{2} is related to a likelihood by χ^{2} = −2ln(likelihood) (i.e., the same form as S; Willott et al. 2001), we can define a new function S_{all}, which combines the constraints from all three types of data (i.e., the SFG sample, LRLF, and SC data). The expression is as follows
where and denote the value of χ^{2} for the local radio LFs and source counts, respectively. Because we use three different types of data to estimate S_{all}, we need to balance the statistical weight for each term in Eq. (7). We chose an A_{0} so that the value of is approximately equal to that of . 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 Eq. (7), we can obtain the bestfit parameters for LFs by numerically minimizing the objective function S_{all}. Here we adopt a Bayesian method as in our previous papers (e.g., Yuan et al. 2016). This latter 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 (socalled “uninformative”) priors on the parameters, and employ the MCMC sampling algorithm available in the Python package EMCEE (ForemanMackey et al. 2013) to estimate the bestfit parameters.
3.2. Local luminosity functions and radio source counts
The local LFs at 1.4 GHz have been well determined for SFGs thanks 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 centimeters), and largearea spectroscopic surveys. In the present work, we simultaneously use the local SFG LFs from Condon et al. (2002, 2019), Best et al. (2005), and Mauch & Sadler (2007; see Fig. 2) to calculate in Eq. (7).
Fig. 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. 
In addition to local LFs, the observed radio source counts can provide an important constraint to the modeling of SFG LFs. In the past several years, deep radio surveys have emerged that reach submJy detection limits, enabling investigation of the faint source counts (e.g., Smolčić et al. 2017b; Ocran et al. 2020b; Mandal et al. 2021; van der Vlugt et al. 2021; Matthews et al. 2021b; Hale et al. 2023). The source counts, denoted n(F_{ν}), represent the number of sources per flux density (F_{ν}) per steradian. The shape of n(F_{ν}) 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 by (e.g., de Zotti et al. 2010). According to Padovani (2016) and Yuan et al. (2017), we can relate the source counts of SFGs to their LF using the following equation:
where c is the speed of light, Φ(z, L) is the LF, D_{L}(z) is the luminosity distance, z_{min} and z_{max} represent the range of integration in redshift, and α is the spectral index.
In this work, we use the observed source counts from Algera et al. (2020a) and Hale et al. (2023) to provide an additional constraint in our analysis. The Algera et al. (2020a) 3 GHz source counts are measured based on the ultrafaint (reaching a 5σ flux limit of ∼2.7 μJy beam^{−1} within the center of the 3 GHz image) radio population detected in the Karl G. Jansky Very Large Array COSMOSXS 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 XMMLSS 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 Table 1 (for the COSMOS field) and Table 2 (for the XMMLSS field) presented by these latter authors, respectively. In this work, we use the first case, where the unclassified or unmatched sources are assumed to have the same split between SFGs and AGN as the classified sources at the given flux density. The source counts are shown in the SC_{SFG, ratio} column in Tables 1 and 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 luminosity function of starforming galaxies
Without loss of generality, the SFG LF can be written as
where e_{1}(z) and e_{2}(z) denote the density evolution (DE) and luminosity evolution (LE) functions of redshift, respectively, and η^{j} represents the parameters that determine the shape of the LF. If the values of η^{j} are constant, this indicates that the shape of the radio LF is unchanged with redshift. Conversely, if η^{j} exhibits a redshift dependence, this implies luminositydependent density evolution (see Singal et al. 2013, 2014 for a more detailed discussion). We assume the shape of the LF to remain unchanged (i.e., η^{j} is constant) as in many other studies (e.g., Novak et al. 2017; van der Vlugt et al. 2022).
Following previous work (e.g., Smolčić et al. 2009a; Gruppioni et al. 2013; van der Vlugt et al. 2022), the SFG local LF ϕ(z = 0, L/e_{2}(z = 0)) is described by a modifiedSchechter function from Saunders et al. (1990):
where L_{⋆} determines the location of the knee in the LF, β and γ fit the faint and bright ends of the LF, respectively, and Φ_{⋆} is used for the normalization. In this work, we consider three LF models, all of which adopt the same LE function:
The DE function e_{1}(z) has three different forms depending on the model: e_{1}(z) = 1 for Model A, an exponential form
for Model B, and
for Model C. In the 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. Model selection
In order to evaluate which model is a better fit to the data, a helpful tool is the information criterion (Takeuchi 2000). The Akaike Information Criterion (AIC; Akaike 1974) is one of the most widely used information criterion. For our problem, the AIC can be written as:
where S_{all} is given in Eq. (7), is the bestfit model parameters, and q is the number of parameters for each model. The model with the smallest value of AIC is considered to be the most accurate. Another commonly used criterion is the Bayesian information criterion (BIC; Schwarz 1978), which can be written as
where n is the sample size. When calculating the S_{all} values for our three models, the weight factor A_{0} in Eq. (7) is set to 1. The AIC and BIC values are listed in Table 1. We find that the AIC and BIC are consistent with each other, both indicating that the LADE model is superior to the PLE model. The AIC value of Model B is slightly smaller than that of Model C, implying that Model B could be taken as our preferred model.
Values of AIC and BIC for Models A–C.
4. Results
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 ForemanMackey 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 twodimensional posterior probability distributions of the parameters for Models A–C are shown in Figs. 3–5, respectively. These corner plots illustrate that all the parameters for our three models are well constrained. Table 2 reports the bestfit parameters and their 1σ uncertainties for the three models.
Fig. 3. Corner plot showing the one and twodimensional projections of the posterior probability distributions of the parameters for Model A obtained from the MCMC run. The histograms on the diagonal show the marginalized posterior densities for each parameter (vertical dashed lines denote the 16th and 84th percentiles). The offdiagonal panels show the 2D joint posterior densities of all couples of parameters, with 1σ, 2σ, and 3σ contours shown by black solid lines. Our bestfitting parameters are marked by red vertical solid lines. 
Bestfit parameters for Models A–C.
Figure 6 shows our bestfit LFs for Models A (solid blue lines), B (solid orange lines), and C (solid green lines). All the LFs are measured at the restframe 1.4 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 leftpointing triangles, dark blue circles, and sky blue squares with error bars, respectively. At lower redshifts of z < 1.0, our three models are barely distinguishable. As redshift increases, Model A begins to diverge, predicting larger number densities at L < L_{⋆} than Models B and C. This deviation increases towards higher redshifts, and the disagreement is in excess of the 3σ confidence intervals at z > ∼1. We also note that, at z > ∼2, the binned LFs of van der Vlugt et al. (2022) present a decline in number density at the faint end, while Model A cannot reproduce the behavior. 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 LFs for all redshift intervals.
Fig. 6. Radio LFs of SFGs at various redshifts compared with the previous estimates specified in the inset. The bestfit LFs for Models A–C in each redshift bin are shown by the blue, orange, and green solid lines, respectively. The light green shaded area shows the 3σ confidence interval for Model C. The red dashed curves represent the KDE LFs (see Appendix A), with the 3σ confidence interval shown by the pink shaded area. The purple dashdotted 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). 
In Fig. 6, 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 given that 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 using the traditional 1/V_{max} method. Although their modeling result would be dependent on the estimation accuracy of 1/V_{max}. Our analytical LFs are obtained through a full maximumlikelihood analysis, and are therefore independent of the 1/V_{max} estimates.
4.2. Fitting the observed source counts
We calculate the source counts for each of our three model LFs using Eq. (8). The result is shown in Fig. 7 as blue dashdotted line, orange dashed line, and green dotted line for Models A–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. (2020a). All our three models can reproduce the Algera et al. (2020a) 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 be caused by multiple factors, such as fieldtofield variation, differences in the assumptions used to calculate completeness, and resolution bias (see Hale et al. 2023). We note that, at the bright region (F_{ν} > 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.
Fig. 7. Comparison of our bestfit models with the Euclidean normalized 1.4 GHz source counts for SFGs observed in the literature. The blue dashdotted line, orange dashed line, and green dotted line show our bestfit source counts of Models A–C, respectively. The source counts from Hale et al. (2023) in the COSMOS and XMMLSS fields are shown as purple right triangles and light blue pentagons, respectively. Also shown are the observed source counts from Smolčić et al. (2017a, green squares), Algera et al. (2020a, orange circles). 
4.3. LADE versus PLE
Due to the limitation of survey depth, most of the existing radio studies barely reach the knee of the SFG LF at z > 1. If fitting the LF using a LADE model, the DE and LE parameters may become degenerate (also see van der Vlugt et al. 2022). On several occasions, this has led previous authors to assume a PLE model for their SFG LFs. By combining data from the ultradeep COSMOSXS survey and the shallower VLACOSMOS 3 GHz large project, van der Vlugt et al. (2022) was able to jointly constrain the LE and DE, finding evidence for significant DE. From Fig. 6, 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 using data from the Low Frequency Array (LOFAR) Two Metre Sky Survey in three wellstudied extragalactic fields, ElaisN1, Boötes, and the Lockman Hole. By fixing the faint and bright end shape of the radio LF to the local values (equivalent to fixing γ = 0.49 and β = 1.12 in Eq. (10)), these latter authors fitted the LF points (via 1/V_{max}) to find the bestfit L_{⋆} and ϕ_{⋆} at each individual redshift bin. Their Fig. 7 shows the variation of L_{⋆} and ϕ_{⋆} as functions of redshift. The method of Cochrane et al. (2023) is equivalent to assuming a LADE model for their LFs. The variation of L_{⋆} and ϕ_{⋆} as functions of redshift in their analysis correspond to L_{⋆} × e_{2}(z) and ϕ_{⋆} × e_{1}(z) in our work. In Fig. 8, we show L_{⋆} × e_{2}(z) and ϕ_{⋆} × e_{1}(z) for our three LF models compared with the inference from Cochrane et al. (2023). We converted our L_{⋆} to 150 MHz by assuming a spectral index of 0.7. Below z < 1.5, our two LADE models – especially Model C – agree well with the L_{⋆} and ϕ_{⋆} evolutions given by Cochrane et al. (2023). Above z > 1.5, the L_{⋆} evolution curve obtain by these authors is significantly higher than those of our models (outside the 3σ uncertainties), while their ϕ_{⋆} evolution falls more rapidly than those of our models out to high redshift.
Fig. 8. L_{⋆} × e_{2}(z) (left panel) and ϕ_{⋆} × 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σ error bands. Our L_{⋆} and ϕ_{⋆} have been converted to 150 MHz by assuming a spectral index of 0.7. The purple circles represent the variation of L_{⋆} (left panel) and ϕ_{⋆} (right panel) as functions of redshift inferred by Cochrane et al. (2023). 
The discrepancies could be explained as follows: At any redshift bin, the bestfit L_{⋆} and ϕ_{⋆} are negatively correlated with each other (e.g., see the right panel of Fig. 4 in Cochrane et al. 2023), implying that the LE and DE parameters are degenerate. The degeneracy should be stronger at higher redshift where the knee location of LFs is increasingly difficult to identify. We highlight the fact that Cochrane et al. (2023) fitted the LF points (via 1/V_{max}) for each redshift bin individually. At higher redshift, the LF points estimated by 1/V_{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_{⋆} and ϕ_{⋆} values. Therefore, their L_{⋆} and ϕ_{⋆} 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 is free from the 1/V_{max} estimates.
Cochrane et al. (2023) found that their ϕ_{⋆} remains roughly constant back to z ∼ 0.8 but then falls steeply at higher redshifts. Our Model C displays a similar trend. A comparison between the inference of Cochrane et al. (2023) and our models lends strong support to the efficacy of our LADE models.
5. History of the cosmic star formation rate density
5.1. Calculating the SFRD
Now that we have obtained the restframe 1.4 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:
where f_{IMF} is a factor accounting for the IMF (f_{IMF} = 1 for a Chabrier 2003 IMF and f_{IMF} = 1.7 for a Salpeter 1955 IMF), and L_{1.4 GHz} is restframe 1.4 GHz radio luminosity. Following Novak et al. (2017), we use the Chabrier IMF in the following analysis. In Eq. (16), q_{IR}(z) is the FIRtoradio luminosity ratio, which is conventionally used to parametrize the FIR–radio correlation in SFGs, and is defined as
where L_{IR} is the total IR luminosity (restframe 8−1000 μm), and 3.75 × 10^{12} Hz represents the central frequency over the farinfrared domain.
Although q_{IR} is typically taken to be a constant value derived for local galaxies, recent observations suggest that the q_{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_{IR} as a function of both stellar mass (M_{⋆}) and redshift. These latter authors found that q_{IR} primarily evolves with M_{⋆}, and only weakly with redshift. This finding implies that using radio emission as an SFR tracer requires M_{⋆}dependent conversion factors, but its robustness still needs to be verified with further study. In this work, we use the expression given by Novak et al. (2017):
which is the updated version of the analysis of Delhaize et al. (2017), using a new sample selection criteria to exclude AGN.
The SFRD of a given epoch can then be estimated by the following integral:
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 bestfit evolution parameters presented in Table 2 and Fig. 6. 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 →∞. In Fig. 9a, our SFRD history results are shown as blue, orange, and green solid lines for Models A–C, respectively. The light shaded areas take into account the 3σ error bands. Our three models coincide at lower redshift of z < 1. Model A gradually separates from the other two models towards higher redshift. Models B and C start to diverge at z > ∼2, but the difference is within the 3σ confidence intervals.
Fig. 9. History of the cosmic SFRD. Our SFRD history results are shown as blue, orange, and green solid lines for Models A–C, respectively. The light shaded areas take into account the 3σ 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.2. Comparison with the literature
In Fig. 9, 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), who performed a fit 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 Fig. 9a, we find that below z ∼ 1.5, the rate at which the SFRD increases with redshift in our three models shows striking similarity to the trend observed by Madau & Dickinson (2014), although their SFRD is slightly higher. The discrepancy in magnitude is probably due to the assumed FIR–radio relation in our result, which calibrates radio emission as a tracer of SFR. Above z ∼ 1.5, 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 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. A similar behavior was observed in recent radio estimates by van der Vlugt et al. (2022).
In Fig. 9a, we also show our Model A SFRD compared to the radio estimates from Smolčić et al. (2009b) and Novak et al. (2017). Smolčić et al. (2009b) 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 an agreement with the Smolčić et al. (2009b) estimates within the permissible error ranges, and therefore their result provides a good consistency check for 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, because our analysis is based on the sample studied by these latter authors. Overall, the curve of our Model A SFRD seems to be a good fit to their SFRD points. Nevertheless, there are two points at z < 2 and one point at z > 2 that seem to disagree with our Model A at the 3σ level. Because we use the same q_{IR}(z) evolution to calculate the SFRD as Novak et al. (2017) did, any discrepancy between the two results can only arise from the difference in LF (see Fig. 6). As discussed in Sect. 4.1, the analytical LFs of Novak et al. (2017) are obtained by fitting the 1/V_{max} LF points in all redshift bins, while our LFs are obtained through a full maximumlikelihood analysis, incorporating 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 an improvement on their estimates.
In Fig. 9b, we show the SFRD derived from our Models B and C compared to the radio estimates from van der Vlugt et al. (2022). These latter authors also assumed a LADE LF – similar to our Models B and C – to calculate their SFRD based on the combined data set from the ultradeep COSMOSXS survey and the VLACOSMOS 3 GHz large project. Although adopting a different q_{IR}(z) evolution from that used here, their result is in good agreement with those of our Models B and C – especially our Model C – within the error bars. We note that three points of the estimates of van der Vlugt et al. (2022) show a slightly elevated SFRD compared to our model prediction. The discrepancies could be attributed to the different q_{IR}(z) evolution used by these authors, or the uncertainties in their LF measurement propagated from the 1/V_{max} estimator. Similar to Novak et al. (2017), van der Vlugt et al. (2022) also obtained their analytical LFs by fitting the 1/V_{max} LF points in all redshift bins.
Figure 9c shows 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 approaches taken; these latter authors performed stacking on massselected galaxies and used a nonevolving FIR–radio correlation established by Bell (2003). The measurements of Cochrane et al. (2023) are systematically higher than ours. We find that a vertical shift of our Model B SFRD curve will match the Cochrane et al. (2023) data points over the whole redshift range. This is equivalent to multiplying our SFRD by ∼1.5, shown as the orange dashed curve. Different from the q_{IR}(z)based L_{radio} − SFR calibration used in the present 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 MHz} − SFR relation. This may explain the discrepancy between their measurements and ours. As noted by Leslie et al. (2020), the impact of different L_{radio} − SFR calibrations is significant.
5.3. Density evolution is indispensable
According to Yuan et al. (2016), the evolution of a LF may be regarded as a vector E, and can be written as
where E_{d} and E_{l} are the base vectors of DE and LE, respectively; and e_{1} and e_{2} are DE and LE functions as mentioned in Eq. (9). DE carries a physical meaning, and can tell us whether the sources are more or less numerous than those of today, while the LE can tell us whether the sources are systematically more or less luminous than those of today. In all three of our models, the LE function has a peak, indicating that SFGs are, on average, most luminous in radio at z ∼ 3 − 4. The DE function, according to our LADE models, monotonically decreases with redshift, implying that, in the radio view, SFGs are less numerous in earlier epochs. From Eq. (18), we speculate that the shape of the SFRD curve is jointly determined by the form of DE and LE. Comparing Figs. 8 and 9, we note that below z ∼ 2, for all three of our models, the effect of positive LE is dominant over the DE. Therefore, the SFRD curves display a monotonically increasing trend. Above z ∼ 2, the LE gradually begins to turn over, and the effect of DE begins to show up. As Model C has the strongest negative DE, its SFRD falls more rapidly out to high redshift than Models B and A. One of the main findings of this work is that a DE is genuinely indispensable in modeling the evolution of SFG radio LFs. This finding is further confirmed by the picture that, the assumption of pure LE seems to overpredict the SFRD at high z, while the inclusion of DE corrects for this.
5.4. The effect of IR–radio correlation
As shown in Eq. (18), the calculation of SFRD relies on two components: the LF and the derived SFR. Although our LF has been well constrained (see Sect. 4), the calibration of SFR using different scaling factors may also affect the final result. In this section, we show that the SFRD derived from our Model B LF can change depending on the choice of q_{IR}(z). The first q_{IR}(z) we tested is that from Bell (2003), where a constant q_{IR} value of 2.64 was assumed. In Fig. 9d, the derived SFRD is shown as the gray solid curve. The model seems to significantly overpredict the SFRD above z ∼ 1.5. The second q_{IR}(z) tested is that from Magnelli et al. (2015), where q_{FIR}(z) = 2.35 × (1 + z)^{−0.12}; this relation can be scaled as log(L_{FIR}) = log(L_{IR})−log(2) to obtain the q_{IR}(z). The SFRD based on this q_{IR}(z) is presented in the same panel as the gray dashed line, which is generally consistent with our result. Finally, we considered q_{IR}(z) = 2.20 from Algera et al. (2020b). The resultant SFRD is shown as the gray dashdotted line, and is consistent with our result at z > 3.5, but is significantly lower at z < 3.5. In conclusion, we show that the assumed IR–radio correlation has a significant impact on the derived SFRD and is crucial to accurately constrain the q_{IR} value at all observed redshifts (also see Novak et al. 2017).
6. Summary and conclusions
In this work, we make use of a starforming galaxy (SFG) sample (5900 sources) from the VLACOSMOS 3 GHz data (Smolčić et al. 2017a) to measure the radio luminosity functions (LFs). For the first time, we use both models of pure luminosity evolution (PLE) and joint luminosity+density evolution (LADE) to fit the LFs directly to the radio data using a full maximumlikelihood analysis. We 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 applied the kernel density estimation (KDE) method described in our previous work (Yuan et al. 2020, 2022). Based on these radio LFs, we derived the dustunbiased star formation rate density (SFRD) up to a redshift of z ∼ 5. The main results of our work are as follows.

Our radio LFs are fitted using three models, assuming a modifiedSchechter 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. 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, and well reproduce the latest radio source counts from Algera et al. (2020a) and Hale et al. (2023). The Akaike information criterion (AIC) also demonstrates that the LADE model is superior to the PLE model.

We find that the luminosity and density evolutions of our LADE models are broadly consistent with the L_{⋆} and ϕ_{⋆} evolutions inferred by Cochrane et al. (2023) within the uncertainty due to degeneracy of LE and DE parameters. As the inference of Cochrane et al. (2023) is model independent, the fact that their results are consistent with ours lends strong support to the efficacy of our LADE models, and we conclude that density evolution is genuinely indispensable in modeling the evolution of SFG radio LFs.

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_{max} LF points in all redshift bins, and possibly suffer from bias in the 1/V_{max} estimate itself. Our analytical LFs are independent of the 1/V_{max} estimates and should be more accurate, which would make our SFRD results an improvement on these previous estimates.

Below z ∼ 1.5, the SFRD from both our PLE and LADE models is a good match to the multiwavelength compilation from Madau & Dickinson (2014), if considering the calibration uncertainty of the FIR–radio relation. Their SFRD curve shows a very similar gradient to ours. This provides a good consistency check for radio emission as a tracer of SFR at low redshift. Comparing with Madau & Dickinson (2014), the SFRD predicted by our three models turns over at a slightly higher redshift (2 < z < 2.5) and falls more rapidly out to high redshift. This was also observed for recent radio estimates by van der Vlugt et al. (2022).

The very recent measurements from Cochrane et al. (2023) are systematically higher than ours (by up to 0.4 dex). This discrepancy could be due to the different L_{radio} − SFR calibration and its scatter correction in the work of these latter authors. From the comparison with Cochrane et al. (2023), we highlight that the L_{radio} − SFR calibration and its scatter, or choice of q_{IR} and its evolution, has a significant impact and remains one of the biggest uncertainties in the radio estimates of SFRD.
Acknowledgments
We thank the anonymous reviewer for the many constructive comments and suggestions, leading to a clearer description of these results. 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 Nos. 12073069, 12075084, and 12393813). Z.Y. is supported by the Xiaoxiang Scholars Programme of Hunan Normal University. J.M. is supported by the China Manned Space Project (CMSCSST2021A06), and the Yunnan Revitalization Talent Support Program (YunLing Scholar project). We would like to thank Yaqian Yu and Yang Liu for helpful discussions.
References
 Adams, N. J., Bowler, R. A. A., Jarvis, M. J., et al. 2020, MNRAS, 494, 1771 [NASA ADS] [CrossRef] [Google Scholar]
 Aird, J., Nandra, K., Laird, E. S., et al. 2010, MNRAS, 401, 2531 [Google Scholar]
 Akaike, H. 1974, IEEE Trans. Autom. Control, 19, 716 [Google Scholar]
 Algera, H. S. B., van der Vlugt, D., Hodge, J. A., et al. 2020a, ApJ, 903, 139 [Google Scholar]
 Algera, H. S. B., Smail, I., Dudzevičiūtė, U., et al. 2020b, ApJ, 903, 138 [Google Scholar]
 Bell, E. F. 2003, ApJ, 586, 794 [Google Scholar]
 Best, P., Kauffmann, G., Heckman, T., & Ivezić, Ž. 2005, MNRAS, 362, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Bonato, M., Negrello, M., Mancuso, C., et al. 2017, MNRAS, 469, 1912 [Google Scholar]
 Bonato, M., Prandoni, I., De Zotti, G., et al. 2021a, A&A, 656, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bonato, M., Prandoni, I., De Zotti, G., et al. 2021b, MNRAS, 500, 22 [Google Scholar]
 Bouwens, R. J., Illingworth, G. D., Franx, M., et al. 2009, ApJ, 705, 936 [Google Scholar]
 Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34 [Google Scholar]
 Bouwens, R. J., Oesch, P. A., Stefanon, M., et al. 2021, AJ, 162, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Bowler, R. A. A., Dunlop, J. S., McLure, R. J., et al. 2015, MNRAS, 452, 1817 [Google Scholar]
 Bowler, R. A. A., Jarvis, M. J., Dunlop, J. S., et al. 2020, MNRAS, 493, 2059 [Google Scholar]
 Calistro Rivera, G., Williams, W. L., Hardcastle, M. J., et al. 2017, MNRAS, 469, 3468 [Google Scholar]
 Ceraj, L., Smolčić, V., Delvecchio, I., et al. 2018, A&A, 620, A192 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
 Clemens, M. S., Vega, O., Bressan, A., et al. 2008, A&A, 477, 95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cochrane, R. K., Kondapally, R., Best, P. N., et al. 2023, MNRAS, 523, 6082 [NASA ADS] [CrossRef] [Google Scholar]
 Condon, J. J. 1992, ARA&A, 30, 575 [Google Scholar]
 Condon, J. J., Cotton, W. D., & Broderick, J. J. 2002, AJ, 124, 675 [Google Scholar]
 Condon, J. J., Matthews, A. M., & Broderick, J. J. 2019, ApJ, 872, 148 [NASA ADS] [CrossRef] [Google Scholar]
 de Zotti, G., Massardi, M., Negrello, M., & Wall, J. 2010, A&ARv, 18, 1 [CrossRef] [Google Scholar]
 Delhaize, J., Smolčić, V., Delvecchio, I., et al. 2017, A&A, 602, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delvecchio, I., Smolčić, V., Zamorani, G., et al. 2017, A&A, 602, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delvecchio, I., Daddi, E., Sargent, M. T., et al. 2021, A&A, 647, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2020, MNRAS, 494, 3828 [Google Scholar]
 Enia, A., Talia, M., Pozzi, F., et al. 2022, ApJ, 927, 204 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, X., Strauss, M. A., Schneider, D. P., et al. 2001, AJ, 121, 54 [NASA ADS] [CrossRef] [Google Scholar]
 Finkelstein, S. L., Ryan, R. E., Jr., Papovich, C., et al. 2015, ApJ, 810, 71 [NASA ADS] [CrossRef] [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
 Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23 [Google Scholar]
 Gruppioni, C., Béthermin, M., Loiacono, F., et al. 2020, A&A, 643, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hale, C. L., Whittam, I. H., Jarvis, M. J., et al. 2023, MNRAS, 520, 2668 [NASA ADS] [CrossRef] [Google Scholar]
 Helou, G., Soifer, B. T., & RowanRobinson, M. 1985, ApJ, 298, L7 [Google Scholar]
 Hogg, D. W. 1999, ArXiv eprints [arXiv:astroph/9905116] [Google Scholar]
 Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2018, ApJ, 854, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Jarvis, M. J., Smith, D. J. B., Bonfield, D. G., et al. 2010, MNRAS, 409, 92 [Google Scholar]
 Jarvis, M., Seymour, N., Afonso, J., et al. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 68 [Google Scholar]
 Karim, A., Schinnerer, E., MartínezSansigre, A., et al. 2011, ApJ, 730, 61 [Google Scholar]
 Kennicutt, R. C., Jr. 1998, ARA&A, 36, 189 [Google Scholar]
 Kochanek, C. S. 1996, ApJ, 473, 595 [NASA ADS] [CrossRef] [Google Scholar]
 LaraLópez, M., Cepa, J., Bongiovanni, A., et al. 2010, A&A, 521, L53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leslie, S. K., Schinnerer, E., Liu, D., et al. 2020, ApJ, 899, 58 [Google Scholar]
 Lewis, A., & Bridle, S. 2002, Phys. Rev. D, 66, 103511 [Google Scholar]
 Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [Google Scholar]
 Magnelli, B., Ivison, R., Lutz, D., et al. 2015, A&A, 573, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Malefahlo, E. D., Jarvis, M. J., Santos, M. G., et al. 2022, MNRAS, 509, 4291 [Google Scholar]
 Mandal, S., Prandoni, I., Hardcastle, M. J., et al. 2021, A&A, 648, A5 [EDP Sciences] [Google Scholar]
 Marshall, H., Tananbaum, H., Avni, Y., & Zamorani, G. 1983, ApJ, 269, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Matthews, A. M., Condon, J. J., Cotton, W. D., & Mauch, T. 2021a, ApJ, 914, 126 [CrossRef] [Google Scholar]
 Matthews, A. M., Condon, J. J., Cotton, W. D., & Mauch, T. 2021b, ApJ, 909, 193 [NASA ADS] [CrossRef] [Google Scholar]
 Mauch, T., & Sadler, E. M. 2007, MNRAS, 375, 931 [Google Scholar]
 McLeod, D. J., McLure, R. J., Dunlop, J. S., et al. 2015, MNRAS, 450, 3032 [NASA ADS] [CrossRef] [Google Scholar]
 McLeod, D. J., McLure, R. J., & Dunlop, J. S. 2016, MNRAS, 459, 3812 [NASA ADS] [CrossRef] [Google Scholar]
 McLure, R. J., Dunlop, J. S., Bowler, R. A. A., et al. 2013, MNRAS, 432, 2696 [Google Scholar]
 Novak, M., Smolčić, V., Delhaize, J., et al. 2017, A&A, 602, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ocran, E. F., Taylor, A. R., Vaccari, M., et al. 2020a, MNRAS, 491, 5911 [NASA ADS] [CrossRef] [Google Scholar]
 Ocran, E. F., Taylor, A. R., Vaccari, M., IshwaraChandra, C. H., & Prandoni, I. 2020b, MNRAS, 491, 1127 [Google Scholar]
 Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labbé, I., & Stefanon, M. 2018, ApJ, 855, 105 [Google Scholar]
 Ono, Y., Ouchi, M., Harikane, Y., et al. 2018, PASJ, 70, S10 [Google Scholar]
 Padovani, P. 2016, A&ARv, 24, 13 [Google Scholar]
 Parsa, S., Dunlop, J. S., McLure, R. J., & Mortlock, A. 2016, MNRAS, 456, 3194 [Google Scholar]
 Riechers, D. A., Bradford, C. M., Clements, D. L., et al. 2013, Nature, 496, 329 [Google Scholar]
 RowanRobinson, M., Oliver, S., Wang, L., et al. 2016, MNRAS, 461, 1100 [Google Scholar]
 Sadler, E. M., Jenkins, C. R., & Kotanyi, C. G. 1989, MNRAS, 240, 591 [NASA ADS] [Google Scholar]
 Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
 Sargent, M. T., Schinnerer, E., Murphy, E., et al. 2010, ApJS, 186, 341 [Google Scholar]
 Saunders, W., RowanRobinson, M., Lawrence, A., et al. 1990, MNRAS, 242, 318 [Google Scholar]
 Schmidt, M. 1968, ApJ, 151, 393 [Google Scholar]
 Schwarz, G. 1978, Ann. Stat., 6, 461 [Google Scholar]
 Singal, J., Petrosian, V., Lawrence, A., et al. 2013, ApJ, 764, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Singal, J., Ko, A., & Petrosian, V. 2014, ApJ, 786, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Smail, I., Ivison, R. J., & Blain, A. W. 1997, ApJ, 490, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, D. J. B., Haskell, P., Gürkan, G., et al. 2021, A&A, 648, A6 [EDP Sciences] [Google Scholar]
 Smolčić, V., Zamorani, G., Schinnerer, E., et al. 2009a, ApJ, 696, 24 [CrossRef] [Google Scholar]
 Smolčić, V., Schinnerer, E., Zamorani, G., et al. 2009b, ApJ, 690, 610 [CrossRef] [Google Scholar]
 Smolčić, V., Novak, M., Bondi, M., et al. 2017a, A&A, 602, A1 [Google Scholar]
 Smolčić, V., Delvecchio, I., Zamorani, G., et al. 2017b, A&A, 602, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tabatabaei, F. S., Schinnerer, E., Krause, M., et al. 2017, ApJ, 836, 185 [Google Scholar]
 Takeuchi, T. T. 2000, Astrophys. Space Sci., 271, 213 [Google Scholar]
 Upjohn, J. E., Brown, M. J. I., Hopkins, A. M., & Bonne, N. J. 2019, PASA, 36, e012 [Google Scholar]
 van der Vlugt, D., Algera, H. S. B., Hodge, J. A., et al. 2021, ApJ, 907, 5 [NASA ADS] [CrossRef] [Google Scholar]
 van der Vlugt, D., Hodge, J. A., Algera, H. S. B., et al. 2022, ApJ, 941, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Willott, C. J., Rawlings, S., Blundell, K. M., Lacy, M., & Eales, S. A. 2001, MNRAS, 322, 536 [NASA ADS] [CrossRef] [Google Scholar]
 Yuan, Z., & Wang, J. 2013, Ap&SS, 345, 305 [NASA ADS] [CrossRef] [Google Scholar]
 Yuan, Z., Wang, J., Zhou, M., & Mao, J. 2016, ApJ, 820, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Yuan, Z., Wang, J., Zhou, M., Qin, L., & Mao, J. 2017, ApJ, 846, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Yuan, Z., Jarvis, M. J., & Wang, J. 2020, ApJS, 248, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Yuan, Z., Zhang, X., Wang, J., Cheng, X., & Wang, W. 2022, ApJS, 260, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Yun, M. S., Reddy, N. A., & Condon, J. J. 2001, ApJ, 554, 803 [Google Scholar]
Appendix A: Kernel density estimation
In many previous studies, the nonparametric LFs are usually derived using the traditional 1/V_{max} method of Schmidt (1968). The mathematics behind 1/V_{max} is the histogram, the simplest density estimator used mainly for rapid visualization of results in one or two dimensions. The 1/V_{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). In addition, 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 or any model assumptions, and simultaneously has some advantages and parametric and nonparametric methods. In this work, we apply the KDE method to derive the nonparametric LFs of SFGs. According to Yuan et al. (2020), the KDE LF can be estimated using
where [Z_{1}, Z_{2}] is the redshift range of the sample, n is the number of objects in the redshift range, Ω is the solid angle subtended by the sample, dV/dz is the differential comoving volume per unit of solid angle, and is the density of (x, y), which corresponds to the (z, L) pair in the KDE parameter space. The calculation of 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 A.1), and then the uncertainty estimation on the LFs. The derived LFs are shown in the panels of Figure 6 as red dashed lines with pink contours. The KDE LFs seem to be in good agreement with the binnd LFs, while with smaller uncertainties.
Fig. A.1. Corner plot showing the posterior probability distributions of the bandwidths (h_{1}, h_{2}) of the KDE LF, which were obtained with the routine provided by Yuan et al. (2022). Uncertainties correspond to the 16th and 84th percentiles, while the bestfit value is the median. 
All Tables
All Figures
Fig. 1. Redshift distribution (top) and the scatter plot (bottom) of our SFG sample. The red dashed curve indicates the flux limit line f_{lim 1.4 GHz}(z). 

In the text 
Fig. 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. 

In the text 
Fig. 3. Corner plot showing the one and twodimensional projections of the posterior probability distributions of the parameters for Model A obtained from the MCMC run. The histograms on the diagonal show the marginalized posterior densities for each parameter (vertical dashed lines denote the 16th and 84th percentiles). The offdiagonal panels show the 2D joint posterior densities of all couples of parameters, with 1σ, 2σ, and 3σ contours shown by black solid lines. Our bestfitting parameters are marked by red vertical solid lines. 

In the text 
Fig. 4. Similar to Fig. 3, but for Model B. 

In the text 
Fig. 5. Similar to Fig. 3, but for Model C. 

In the text 
Fig. 6. Radio LFs of SFGs at various redshifts compared with the previous estimates specified in the inset. The bestfit LFs for Models A–C in each redshift bin are shown by the blue, orange, and green solid lines, respectively. The light green shaded area shows the 3σ confidence interval for Model C. The red dashed curves represent the KDE LFs (see Appendix A), with the 3σ confidence interval shown by the pink shaded area. The purple dashdotted 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). 

In the text 
Fig. 7. Comparison of our bestfit models with the Euclidean normalized 1.4 GHz source counts for SFGs observed in the literature. The blue dashdotted line, orange dashed line, and green dotted line show our bestfit source counts of Models A–C, respectively. The source counts from Hale et al. (2023) in the COSMOS and XMMLSS fields are shown as purple right triangles and light blue pentagons, respectively. Also shown are the observed source counts from Smolčić et al. (2017a, green squares), Algera et al. (2020a, orange circles). 

In the text 
Fig. 8. L_{⋆} × e_{2}(z) (left panel) and ϕ_{⋆} × 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σ error bands. Our L_{⋆} and ϕ_{⋆} have been converted to 150 MHz by assuming a spectral index of 0.7. The purple circles represent the variation of L_{⋆} (left panel) and ϕ_{⋆} (right panel) as functions of redshift inferred by Cochrane et al. (2023). 

In the text 
Fig. 9. History of the cosmic SFRD. Our SFRD history results are shown as blue, orange, and green solid lines for Models A–C, respectively. The light shaded areas take into account the 3σ 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. 

In the text 
Fig. A.1. Corner plot showing the posterior probability distributions of the bandwidths (h_{1}, h_{2}) of the KDE LF, which were obtained with the routine provided by Yuan et al. (2022). Uncertainties correspond to the 16th and 84th percentiles, while the bestfit value is the median. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.