Constraints on cosmic star formation history via a new modeling of the radio luminosity function of star-forming galaxies

Radio wavelengths offer a unique possibility to trace the total star-formation rate (SFR) in galaxies, both obscured and unobscured. To probe the dust-unbiased star-formation history, an accurate measurement of the radio luminosity function (LF) for star-forming galaxies (SFGs) is crucial. 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\sim5$ onwards. 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 updated observations of local radio LFs and radio source counts into the fitting process to obtain additional constraints. 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. 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_{\rm max}$ LF points, our SFRD results should be an improvement on these previous estimates. Below $z\sim1.5$, our SFRD matches a published multiwavelength compilation, while our SFRD turns over at a slightly higher redshift ($2


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., Lara-Ló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 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(Gruppioni et al. , 2020;;Rowan-Robinson 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 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. 2015Bowler et al. , 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 far-infrared (FIR) wavelengths.Therefore, FIR emission is ideal for tracing SFR in dust-rich 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 star-formation 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 short-lived (τ ≤ 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 star-forming galaxies (SFGs), known as the FIRradio 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 submilli-Jansky (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. 2017Bonato et al. , 2021a,b;,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 COSMOS-XS 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 maximum-likelihood analysis (e.g., Willott et al. 2001).
In the present paper, we make use of the VLA-COSMOS 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 maximum-likelihood 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 dust-unbiased 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 power-law radio spectrum for SFGs, F ν ∝ ν −α , where F ν is the flux density at frequency ν and α is the spectral index.

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 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 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 VLA-COSMOS 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 0 .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 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 rest-frame 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 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.

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 maximum-likelihood 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 VLA-COSMOS 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 optical-NIR 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 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 × 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 one-dimensional 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 = −2 ln(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 χ 2 LRLF and χ 2 SC 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 A 0 χ 2 SC is approximately equal to that of χ 2 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 Eq. ( 7), we can obtain the best-fit 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 (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.

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 large-area spectroscopic surveys.In the present work, we simultaneously use the local SFG LFs from Condon et al. (2002Condon et al. ( , 2019)), Best et al. (2005), and Mauch & Sadler (2007; see Fig. 2) to calculate χ 2 LRLF in Eq. ( 7).
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   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 F 2.5 ν (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 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 Table 1 (for the COSMOS field) and Table 2 (for the XMM-LSS 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.

Models for the luminosity function of star-forming 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 luminosity-dependent density evolution (see Singal et al. 2013Singal et al. , 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 modified-Schechter 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. 2016Yuan et al. , 2017) ) or luminosity and density evolution (LADE, e.g., Aird et al. 2010) models.

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 best-fit model parameters, and q is the number of parameters for each model.The model A174, page 4 of 14 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.

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 oneand two-dimensional 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 best-fit parameters and their 1σ uncertainties for the three models.
Figure 6 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 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-pointing 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.
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 maximum-likelihood analysis, and are therefore independent of the 1/V max estimates.

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 2017a) result is systematically lower than the others.The discrepancies could be caused by multiple factors, such as field-to-field 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.

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 COSMOS-XS survey and the shallower VLA-COSMOS 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 well-studied extragalactic fields, Elais-N1, 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 best-fit 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.
The discrepancies could be explained as follows: At any redshift bin, the best-fit 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.
A174, page 5 of 14 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.

Calculating the SFRD
Now that we have obtained the rest-frame 1.4 GHz LF, we can investigate how the SFRD evolves with redshift.We use the A174, page 6 of 14 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 rest-frame 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 FIR-to-radio luminosity ratio, which is conventionally used to parametrize the FIR-radio cor-relation in SFGs, and is defined as where L IR is the total IR luminosity (rest-frame 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 A174, page 7 of 14  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):   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 best-fit 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 A174, page 9 of 14 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.

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 maximum-likelihood 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 COSMOS-XS survey and the VLA-COSMOS 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 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.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.

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 A174, page 10 of 14 Madau & Dickinson+2014 This paper: Model B (q = 2.78 × (1 + z) −0.14 ) This paper: Model B -3σ Bell+2003 (q = 2.64) Magnelli+2015 (q = 2.35 × (1 + z) −0.12 ) Algera+2020 (q = 2.20) 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.A174, page 11 of 14 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 over-predict the SFRD at high z, while the inclusion of DE corrects for this.

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 dash-dotted 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).

Summary and conclusions
In this work, we make use of a star-forming galaxy (SFG) sample (5900 sources) from the VLA-COSMOS 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(Yuan et al. , 2022)).Based on these radio LFs, we derived the dust-unbiased star formation rate density (SFRD) up to a redshift of z ∼ 5.The main results of our work are as follows.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.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.2. 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. 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 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.4. 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 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).5.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.

Fig. 1 .
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).
as blue dash-dotted 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. (

Fig. 3 .
Fig. 3. 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 marginalized 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σ, 2σ, and 3σ contours shown by black solid lines.Our best-fitting parameters are marked by red vertical solid lines.
which is the updated version of the analysis ofDelhaize et al. (2017), using a new sample selection criteria to exclude AGN.

Fig. 6 .
Fig. 6.Radio LFs of SFGs at various redshifts compared with the previous estimates specified in the inset.The best-fit 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 dash-dotted line depicts the PLE model ofNovak et al. (2017).We also compare with the binned LFs fromGruppioni et al. (2013),Novak et al. (2017), andvan der Vlugt et al. (2022).

Fig. 7 .
Fig. 7. Comparison of our best-fit models with the Euclidean normalized 1.4 GHz source counts for SFGs observed in the literature.The blue dash-dotted line, orange dashed line, and green dotted line show our best-fit source counts of Models A-C, respectively.The source counts from Hale et al. (2023) in the COSMOS and XMM-LSS 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).

Fig. 8 .
Fig. 8. L × e 2 (z) 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).

Table 1 .
Values of AIC and BIC for Models A-C.