EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 683 (March 2024) A&A, 683 (2024) A174 Full HTML
Open Access
Issue
A&A
Volume 683, March 2024
Article Number A174
Number of page(s) 14
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202347746
Published online 18 March 2024

© The Authors 2024

Licence Creative CommonsOpen 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., 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, 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. 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 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 × 107 yr, Matthews et al. 2021a) stars undergo supernova explosions, the expanding remnants can accelerate the cosmic ray electrons and give rise to synchrotron radiation at a typical frequency of < 30 GHz (e.g., Sadler et al. 1989; Condon 1992; Clemens et al. 2008; Tabatabaei et al. 2017). Radio emission triggered by the above process is empirically found to correlate well with the FIR emission of star-forming galaxies (SFGs), known as the FIR–radio correlation. This correlation holds over five orders of magnitude in 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. 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 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/Vmax estimator. This semi-parametric 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/Vmax 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/Vmax 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 H0 = 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.

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 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 $ 0{{\overset{\prime\prime}{.}}}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,

r = log ( L 1.4 GHz [ W Hz − 1 ] SFR IR [ M ⊙ yr − 1 ] ) > 22 × ( 1 + z ) 0.013 . $$ \begin{aligned} r = \log \left(\frac{L_{1.4\,\mathrm{GHz}}\left[\mathrm{W\,Hz^{-1}}\right]}{\mathrm{SFR}_{\mathrm{IR} }\left[M_{\odot }\,\mathrm{yr}^{-1}\right]}\right) > 22 \times (1+z)^{0.013}. \end{aligned} $$(1)

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

f lim 1.4 GHz ( z ) = 4 π D L 2 ( 1 + z ) 1 − α ( 3 GHz 1.4 GHz ) α F lim 3 GHz , $$ \begin{aligned} f_{\mathrm{lim\,1.4\,GHz}}(z) = \frac{4 \pi D_{\rm L}^{2}}{(1+z)^{1-\alpha }}\left(\frac{3\,\mathrm{GHz}}{1.4\,\mathrm{GHz}}\right)^{\alpha }F_{\mathrm{lim\,3\,GHz}}, \end{aligned} $$(2)

thumbnail Fig. 1.

Redshift distribution (top) and the scatter plot (bottom) of our SFG sample. The red dashed curve indicates the flux limit line flim 1.4 GHz(z).

where DL represents the luminosity distance at redshift z, Flim 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:

Φ ( z , L ) = d 2 N d V d log 10 L · $$ \begin{aligned} \Phi (z,L)=\frac{d^{2}N}{\mathrm{d}V\mathrm{d}\log _{10}L}\cdot \end{aligned} $$(3)

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

S = − 2 ∑ i n ln [ Φ ( z i , L i ) p ( z i , L i ) ] + 2 ∫ ∫ W Φ ( z , L ) p ( z , L ) Ω d V d z d z d L , $$ \begin{aligned} \begin{aligned} S=&-2\!\!\sum _{i}^{n}\!\ln [\Phi (z_{i},L_{i})p(z_{i},L_{i})] \\&+2\!\!\!\int \!\!\!\!\int _{W}\!\!\!\Phi (z,L)p(z,L)\Omega \frac{\mathrm{d}V}{\mathrm{d}z}\mathrm{d}z\mathrm{d}L, \end{aligned} \end{aligned} $$(4)

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

p ( z , L ) = C radio [ F 3 GHz ( z ) ] × C opt ( z ) , $$ \begin{aligned} p(z,L)= C_{\mathrm{radio}}[F_{3\,\mathrm{GHz}}(z)] \times C_{\mathrm{opt}}(z), \end{aligned} $$(5)

where Cradio is the completeness of the VLA-COSMOS 3 GHz radio catalog as a function of the flux density F3 GHz, and Copt is the completeness owing to radio sources without assigned optical-NIR counterparts (Novak et al. 2017). We adopt the calculations of Cradio and Copt 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 two-dimensional (2D) grid of 50 × 50 in the log L − z space. For each grid point (log Li, zi), we can derive its flux density Fi from Li by assuming α = 0.7. We can then estimate the corresponding Cradio and Copt 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

χ 2 = ∑ i = 1 n ( f data i − f mod i σ data i ) 2 , $$ \begin{aligned} \chi ^{2}=\sum _{i=1}^{n}\left(\frac{f_{\mathrm{data}\,i}-f_{\mathrm{mod}\,i}}{\sigma _{\mathrm{data}\,i}}\right)^{2}, \end{aligned} $$(6)

where fdata i represents the value of the data in the ith bin, and fmod 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

S all = S + χ LRLF 2 + A 0 χ SC 2 , $$ \begin{aligned} {S\!}_{\mathrm{all}} = S + \chi ^2_{\mathrm{LRLF}} + A_0 \chi ^2_{\mathrm{SC}}, \end{aligned} $$(7)

where χ LRLF 2 $ \chi^2_{{\rm LRLF}} $ and χ SC 2 $ \chi^2_{\mathrm{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 A0 so that the value of A 0 χ SC 2 $ A_0\chi^2_{\mathrm{SC}} $ is approximately equal to that of χ LRLF 2 $ \chi^2_{{\rm LRLF}} $. This yields values of about 10−40 for our calculations. We find that varying A0 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.

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 large-area 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 χ LRLF 2 $ \chi^2_{{\rm LRLF}} $ in Eq. (7).

thumbnail 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 F ν 2.5 $ F^{2.5}_\nu $ (e.g., de Zotti et al. 2010). According to Padovani (2016) and Yuan et al. (2017), we can relate the source counts of SFGs to their LF using the following equation:

n ( F ν ) 4 π = 4 π c H 0 ∫ z min ( F ν ) z max ( F ν ) Φ ( z , L ( F ν , z ) ) D L 4 ( z ) d z ( 1 + z ) ( 3 − α ) Ω m ( 1 + z ) 3 + Ω Λ , $$ \begin{aligned} \begin{aligned} \frac{n(F_\nu )}{4\pi } = 4\pi \frac{c}{H_0} \int _{z_{\rm min}(F_\nu )}^{z_{\rm max}(F_\nu )} \frac{\Phi (z,L(F_\nu ,z))D_{\rm L}^4(z)\mathrm{d}z}{(1+z)^{(3-\alpha )}\sqrt{\Omega _\mathrm{m} (1+z)^3+\Omega _\Lambda }}, \end{aligned} \end{aligned} $$(8)

where c is the speed of light, Φ(z, L) is the LF, DL(z) is the luminosity distance, zmin and zmax 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 SCSFG,  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 star-forming galaxies

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

Φ ( z , L ) = e 1 ( z ) ϕ ( z = 0 , L / e 2 ( z ) , η j ) , $$ \begin{aligned} \Phi (z,L)=e_1(z)\phi (z=0,L/e_2(z),\eta ^j), \end{aligned} $$(9)

where e1(z) and e2(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. 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/e2(z = 0)) is described by a modified-Schechter function from Saunders et al. (1990):

ϕ ( z = 0 , L / e 2 ( z = 0 ) ) = d N d log 10 L = ϕ ⋆ ( L L ⋆ ) 1 − β exp [ − 1 2 γ 2 log 2 ( 1 + L L ⋆ ) ] , $$ \begin{aligned} \begin{aligned} \phi (z=0,L/e_2(z=0))&=\frac{\mathrm{d}N}{\mathrm{d}\log _{10}L} \\&=\phi _{\star }\left(\frac{L}{L_{\star }}\right)^{1-\beta } \exp \left[-\frac{1}{2 \gamma }^{2} \log ^{2}\left(1+\frac{L}{L_{\star }}\right)\right], \end{aligned} \end{aligned} $$(10)

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:

e 2 ( z ) = ( 1 + z ) k 1 + k 2 z . $$ \begin{aligned} e_2(z)=(1+z)^{k_{1} + k_{2} z}. \end{aligned} $$(11)

The DE function e1(z) has three different forms depending on the model: e1(z) = 1 for Model A, an exponential form

e 1 ( z ) = 10 p 1 z , $$ \begin{aligned} e_1(z)=10^{p_1 z}, \end{aligned} $$(12)

for Model B, and

e 1 ( z ) = ( 1 + z ) p 1 + p 2 z , $$ \begin{aligned} e_1(z)=(1+z)^{p_{1} + p_{2} z}, \end{aligned} $$(13)

for Model C. In the above equations, k1, k2, p1, and p2 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:

AIC = S all ( θ ̂ ) + 2 q , $$ \begin{aligned} \mathrm{AIC} = {S\!}_{\mathrm{all} }(\hat{\theta }) + 2q, \end{aligned} $$(14)

where S​all is given in Eq. (7), θ ̂ $ \hat{\theta} $ is the best-fit 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

BIC ( q ) = S all ( θ ̂ ) + q ln n , $$ \begin{aligned} \mathrm{BIC}(q) = {S\!}_{\mathrm{all} }(\hat{\theta }) + q\mathrm{ln}\,n, \end{aligned} $$(15)

where n is the sample size. When calculating the S​all values for our three models, the weight factor A0 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.

Table 1.

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 Foreman-Mackey et al. (2013). The EMCEE algorithm improves the exploration of the parameter space by using an ensemble of chains with different initial conditions. This approach helps to avoid local minima and ensures a more comprehensive exploration of the parameter space. We assume uniform priors on all the parameters. The marginalized one- and two-dimensional posterior probability distributions of the parameters for Models A–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.

thumbnail 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.

thumbnail Fig. 4.

Similar to Fig. 3, but for Model B.

thumbnail Fig. 5.

Similar to Fig. 3, but for Model C.

Table 2.

Best-fit parameters for Models A–C.

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.

thumbnail 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 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/Vmax method. Although their modeling result would be dependent on the estimation accuracy of 1/Vmax. Our analytical LFs are obtained through a full maximum-likelihood analysis, and are therefore independent of the 1/Vmax 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 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. (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.

thumbnail 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).

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 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/Vmax) 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⋆ × e2(z) and ϕ⋆ × e1(z) in our work. In Fig. 8, we show L⋆ × e2(z) and ϕ⋆ × e1(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.

thumbnail Fig. 8.

L⋆ × e2(z) (left panel) and ϕ⋆ × e1(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 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. 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/Vmax) for each redshift bin individually. At higher redshift, the LF points estimated by 1/Vmax 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/Vmax 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 rest-frame 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:

SFR M ⊙ yr − 1 = f IMF × 10 − 24 10 q IR ( z ) L 1.4 GHz WHz − 1 , $$ \begin{aligned} \frac{\mathrm{SFR} }{M_{\odot }\,\mathrm{yr}^{-1}} = f_{\mathrm{IMF} }\times 10^{-24}10^{q_{\mathrm{IR} }(z)}\frac{L_{1.4\,\mathrm{GHz}}}{\mathrm{WHz}^{-1}}, \end{aligned} $$(16)

where fIMF is a factor accounting for the IMF (fIMF = 1 for a Chabrier 2003 IMF and fIMF = 1.7 for a Salpeter 1955 IMF), and L1.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), qIR(z) is the FIR-to-radio luminosity ratio, which is conventionally used to parametrize the FIR–radio correlation in SFGs, and is defined as

q IR = log ( L IR [ W ] 3.75 × 10 12 [ Hz ] ) − log ( L 1.4 GHz [ W H z − 1 ] ) , $$ \begin{aligned} q_{\mathrm{IR} } = \log \left(\frac{L_{\mathrm{IR} }\left[\mathrm{W} \right]}{3.75\times 10^{12}\left[\mathrm{Hz} \right]}\right)-\log \left(L_{1.4\,\mathrm{GHz}}\left[\mathrm W\,Hz^{-1} \right]\right), \end{aligned} $$

where LIR is the total IR luminosity (rest-frame 8−1000 μm), and 3.75 × 1012 Hz represents the central frequency over the far-infrared domain.

Although qIR is typically taken to be a constant value derived for local galaxies, recent observations suggest that the qIR 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 qIR as a function of both stellar mass (M⋆) and redshift. These latter authors found that qIR 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):

q IR ( z ) = ( 2.78 ± 0.02 ) × ( 1 + z ) − 0.14 ± 0.01 , $$ \begin{aligned} q_{\mathrm{IR} }(z)=(2.78 \pm 0.02) \times (1+z)^{-0.14 \pm 0.01}, \end{aligned} $$(17)

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:

SFRD = ∫ L min L max Φ ( L , z ) × SFR ( L 1.4 GHz ) d log 10 L . $$ \begin{aligned} \mathrm{SFRD} = \int _{L_{\mathrm{min}}}^{L_{\mathrm{max}}} \Phi (L,z) \times \mathrm{SFR}(L_{1.4\,\mathrm{GHz}}) \, \mathrm{d} \log _{10} L. \end{aligned} $$(18)

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 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.

thumbnail 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 qIR(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/Vmax 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 qIR(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 qIR(z) evolution used by these authors, or the uncertainties in their LF measurement propagated from the 1/Vmax estimator. Similar to Novak et al. (2017), van der Vlugt et al. (2022) also obtained their analytical LFs by fitting the 1/Vmax 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 qIR(z)-based Lradio − 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 L150 MHz − SFR relation. This may explain the discrepancy between their measurements and ours. As noted by Leslie et al. (2020), the impact of different Lradio − 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

E = e 1 E d + e 2 E l , $$ \begin{aligned} \boldsymbol{E} = e_1{\boldsymbol{E}}_{\rm d}+e_2{\boldsymbol{E}}_{\rm l}, \end{aligned} $$(19)

where Ed and El are the base vectors of DE and LE, respectively; and e1 and e2 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 over-predict 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 qIR(z). The first qIR(z) we tested is that from Bell (2003), where a constant qIR value of 2.64 was assumed. In Fig. 9d, the derived SFRD is shown as the gray solid curve. The model seems to significantly over-predict the SFRD above z ∼ 1.5. The second qIR(z) tested is that from Magnelli et al. (2015), where qFIR(z) = 2.35 × (1 + z)−0.12; this relation can be scaled as log(LFIR) = log(LIR)−log(2) to obtain the qIR(z). The SFRD based on this qIR(z) is presented in the same panel as the gray dashed line, which is generally consistent with our result. Finally, we considered qIR(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 qIR value at all observed redshifts (also see Novak et al. 2017).

6. 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 maximum-likelihood 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 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/Vmax LF points in all redshift bins, and possibly suffer from bias in the 1/Vmax estimate itself. Our analytical LFs are independent of the 1/Vmax 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 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).

  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 Lradio − 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 Lradio − SFR calibration and its scatter, or choice of qIR 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 (CMS-CSST2021-A06), and the Yunnan Revitalization Talent Support Program (YunLing Scholar project). We would like to thank Yaqian Yu and Yang Liu for helpful discussions.

References

  1. Adams, N. J., Bowler, R. A. A., Jarvis, M. J., et al. 2020, MNRAS, 494, 1771 [NASA ADS] [CrossRef] [Google Scholar]
  2. Aird, J., Nandra, K., Laird, E. S., et al. 2010, MNRAS, 401, 2531 [Google Scholar]
  3. Akaike, H. 1974, IEEE Trans. Autom. Control, 19, 716 [Google Scholar]
  4. Algera, H. S. B., van der Vlugt, D., Hodge, J. A., et al. 2020a, ApJ, 903, 139 [Google Scholar]
  5. Algera, H. S. B., Smail, I., Dudzevičiūtė, U., et al. 2020b, ApJ, 903, 138 [Google Scholar]
  6. Bell, E. F. 2003, ApJ, 586, 794 [Google Scholar]
  7. Best, P., Kauffmann, G., Heckman, T., & Ivezić, Ž. 2005, MNRAS, 362, 9 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bonato, M., Negrello, M., Mancuso, C., et al. 2017, MNRAS, 469, 1912 [Google Scholar]
  9. Bonato, M., Prandoni, I., De Zotti, G., et al. 2021a, A&A, 656, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Bonato, M., Prandoni, I., De Zotti, G., et al. 2021b, MNRAS, 500, 22 [Google Scholar]
  11. Bouwens, R. J., Illingworth, G. D., Franx, M., et al. 2009, ApJ, 705, 936 [Google Scholar]
  12. Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34 [Google Scholar]
  13. Bouwens, R. J., Oesch, P. A., Stefanon, M., et al. 2021, AJ, 162, 47 [NASA ADS] [CrossRef] [Google Scholar]
  14. Bowler, R. A. A., Dunlop, J. S., McLure, R. J., et al. 2015, MNRAS, 452, 1817 [Google Scholar]
  15. Bowler, R. A. A., Jarvis, M. J., Dunlop, J. S., et al. 2020, MNRAS, 493, 2059 [Google Scholar]
  16. Calistro Rivera, G., Williams, W. L., Hardcastle, M. J., et al. 2017, MNRAS, 469, 3468 [Google Scholar]
  17. Ceraj, L., Smolčić, V., Delvecchio, I., et al. 2018, A&A, 620, A192 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
  19. Clemens, M. S., Vega, O., Bressan, A., et al. 2008, A&A, 477, 95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Cochrane, R. K., Kondapally, R., Best, P. N., et al. 2023, MNRAS, 523, 6082 [NASA ADS] [CrossRef] [Google Scholar]
  21. Condon, J. J. 1992, ARA&A, 30, 575 [Google Scholar]
  22. Condon, J. J., Cotton, W. D., & Broderick, J. J. 2002, AJ, 124, 675 [Google Scholar]
  23. Condon, J. J., Matthews, A. M., & Broderick, J. J. 2019, ApJ, 872, 148 [NASA ADS] [CrossRef] [Google Scholar]
  24. de Zotti, G., Massardi, M., Negrello, M., & Wall, J. 2010, A&ARv, 18, 1 [CrossRef] [Google Scholar]
  25. Delhaize, J., Smolčić, V., Delvecchio, I., et al. 2017, A&A, 602, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Delvecchio, I., Smolčić, V., Zamorani, G., et al. 2017, A&A, 602, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Delvecchio, I., Daddi, E., Sargent, M. T., et al. 2021, A&A, 647, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2020, MNRAS, 494, 3828 [Google Scholar]
  29. Enia, A., Talia, M., Pozzi, F., et al. 2022, ApJ, 927, 204 [NASA ADS] [CrossRef] [Google Scholar]
  30. Fan, X., Strauss, M. A., Schneider, D. P., et al. 2001, AJ, 121, 54 [NASA ADS] [CrossRef] [Google Scholar]
  31. Finkelstein, S. L., Ryan, R. E., Jr., Papovich, C., et al. 2015, ApJ, 810, 71 [NASA ADS] [CrossRef] [Google Scholar]
  32. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
  33. Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23 [Google Scholar]
  34. Gruppioni, C., Béthermin, M., Loiacono, F., et al. 2020, A&A, 643, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Hale, C. L., Whittam, I. H., Jarvis, M. J., et al. 2023, MNRAS, 520, 2668 [NASA ADS] [CrossRef] [Google Scholar]
  36. Helou, G., Soifer, B. T., & Rowan-Robinson, M. 1985, ApJ, 298, L7 [Google Scholar]
  37. Hogg, D. W. 1999, ArXiv e-prints [arXiv:astro-ph/9905116] [Google Scholar]
  38. Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2018, ApJ, 854, 73 [NASA ADS] [CrossRef] [Google Scholar]
  39. Jarvis, M. J., Smith, D. J. B., Bonfield, D. G., et al. 2010, MNRAS, 409, 92 [Google Scholar]
  40. Jarvis, M., Seymour, N., Afonso, J., et al. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 68 [Google Scholar]
  41. Karim, A., Schinnerer, E., Martínez-Sansigre, A., et al. 2011, ApJ, 730, 61 [Google Scholar]
  42. Kennicutt, R. C., Jr. 1998, ARA&A, 36, 189 [Google Scholar]
  43. Kochanek, C. S. 1996, ApJ, 473, 595 [NASA ADS] [CrossRef] [Google Scholar]
  44. Lara-López, M., Cepa, J., Bongiovanni, A., et al. 2010, A&A, 521, L53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Leslie, S. K., Schinnerer, E., Liu, D., et al. 2020, ApJ, 899, 58 [Google Scholar]
  46. Lewis, A., & Bridle, S. 2002, Phys. Rev. D, 66, 103511 [Google Scholar]
  47. Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [Google Scholar]
  48. Magnelli, B., Ivison, R., Lutz, D., et al. 2015, A&A, 573, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Malefahlo, E. D., Jarvis, M. J., Santos, M. G., et al. 2022, MNRAS, 509, 4291 [Google Scholar]
  50. Mandal, S., Prandoni, I., Hardcastle, M. J., et al. 2021, A&A, 648, A5 [EDP Sciences] [Google Scholar]
  51. Marshall, H., Tananbaum, H., Avni, Y., & Zamorani, G. 1983, ApJ, 269, 35 [NASA ADS] [CrossRef] [Google Scholar]
  52. Matthews, A. M., Condon, J. J., Cotton, W. D., & Mauch, T. 2021a, ApJ, 914, 126 [CrossRef] [Google Scholar]
  53. Matthews, A. M., Condon, J. J., Cotton, W. D., & Mauch, T. 2021b, ApJ, 909, 193 [NASA ADS] [CrossRef] [Google Scholar]
  54. Mauch, T., & Sadler, E. M. 2007, MNRAS, 375, 931 [Google Scholar]
  55. McLeod, D. J., McLure, R. J., Dunlop, J. S., et al. 2015, MNRAS, 450, 3032 [NASA ADS] [CrossRef] [Google Scholar]
  56. McLeod, D. J., McLure, R. J., & Dunlop, J. S. 2016, MNRAS, 459, 3812 [NASA ADS] [CrossRef] [Google Scholar]
  57. McLure, R. J., Dunlop, J. S., Bowler, R. A. A., et al. 2013, MNRAS, 432, 2696 [Google Scholar]
  58. Novak, M., Smolčić, V., Delhaize, J., et al. 2017, A&A, 602, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Ocran, E. F., Taylor, A. R., Vaccari, M., et al. 2020a, MNRAS, 491, 5911 [NASA ADS] [CrossRef] [Google Scholar]
  60. Ocran, E. F., Taylor, A. R., Vaccari, M., Ishwara-Chandra, C. H., & Prandoni, I. 2020b, MNRAS, 491, 1127 [Google Scholar]
  61. Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labbé, I., & Stefanon, M. 2018, ApJ, 855, 105 [Google Scholar]
  62. Ono, Y., Ouchi, M., Harikane, Y., et al. 2018, PASJ, 70, S10 [Google Scholar]
  63. Padovani, P. 2016, A&ARv, 24, 13 [Google Scholar]
  64. Parsa, S., Dunlop, J. S., McLure, R. J., & Mortlock, A. 2016, MNRAS, 456, 3194 [Google Scholar]
  65. Riechers, D. A., Bradford, C. M., Clements, D. L., et al. 2013, Nature, 496, 329 [Google Scholar]
  66. Rowan-Robinson, M., Oliver, S., Wang, L., et al. 2016, MNRAS, 461, 1100 [Google Scholar]
  67. Sadler, E. M., Jenkins, C. R., & Kotanyi, C. G. 1989, MNRAS, 240, 591 [NASA ADS] [Google Scholar]
  68. Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
  69. Sargent, M. T., Schinnerer, E., Murphy, E., et al. 2010, ApJS, 186, 341 [Google Scholar]
  70. Saunders, W., Rowan-Robinson, M., Lawrence, A., et al. 1990, MNRAS, 242, 318 [Google Scholar]
  71. Schmidt, M. 1968, ApJ, 151, 393 [Google Scholar]
  72. Schwarz, G. 1978, Ann. Stat., 6, 461 [Google Scholar]
  73. Singal, J., Petrosian, V., Lawrence, A., et al. 2013, ApJ, 764, 43 [NASA ADS] [CrossRef] [Google Scholar]
  74. Singal, J., Ko, A., & Petrosian, V. 2014, ApJ, 786, 109 [NASA ADS] [CrossRef] [Google Scholar]
  75. Smail, I., Ivison, R. J., & Blain, A. W. 1997, ApJ, 490, L5 [NASA ADS] [CrossRef] [Google Scholar]
  76. Smith, D. J. B., Haskell, P., Gürkan, G., et al. 2021, A&A, 648, A6 [EDP Sciences] [Google Scholar]
  77. Smolčić, V., Zamorani, G., Schinnerer, E., et al. 2009a, ApJ, 696, 24 [CrossRef] [Google Scholar]
  78. Smolčić, V., Schinnerer, E., Zamorani, G., et al. 2009b, ApJ, 690, 610 [CrossRef] [Google Scholar]
  79. Smolčić, V., Novak, M., Bondi, M., et al. 2017a, A&A, 602, A1 [Google Scholar]
  80. Smolčić, V., Delvecchio, I., Zamorani, G., et al. 2017b, A&A, 602, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Tabatabaei, F. S., Schinnerer, E., Krause, M., et al. 2017, ApJ, 836, 185 [Google Scholar]
  82. Takeuchi, T. T. 2000, Astrophys. Space Sci., 271, 213 [Google Scholar]
  83. Upjohn, J. E., Brown, M. J. I., Hopkins, A. M., & Bonne, N. J. 2019, PASA, 36, e012 [Google Scholar]
  84. van der Vlugt, D., Algera, H. S. B., Hodge, J. A., et al. 2021, ApJ, 907, 5 [NASA ADS] [CrossRef] [Google Scholar]
  85. van der Vlugt, D., Hodge, J. A., Algera, H. S. B., et al. 2022, ApJ, 941, 10 [NASA ADS] [CrossRef] [Google Scholar]
  86. Willott, C. J., Rawlings, S., Blundell, K. M., Lacy, M., & Eales, S. A. 2001, MNRAS, 322, 536 [NASA ADS] [CrossRef] [Google Scholar]
  87. Yuan, Z., & Wang, J. 2013, Ap&SS, 345, 305 [NASA ADS] [CrossRef] [Google Scholar]
  88. Yuan, Z., Wang, J., Zhou, M., & Mao, J. 2016, ApJ, 820, 65 [NASA ADS] [CrossRef] [Google Scholar]
  89. Yuan, Z., Wang, J., Zhou, M., Qin, L., & Mao, J. 2017, ApJ, 846, 78 [NASA ADS] [CrossRef] [Google Scholar]
  90. Yuan, Z., Jarvis, M. J., & Wang, J. 2020, ApJS, 248, 1 [NASA ADS] [CrossRef] [Google Scholar]
  91. Yuan, Z., Zhang, X., Wang, J., Cheng, X., & Wang, W. 2022, ApJS, 260, 10 [NASA ADS] [CrossRef] [Google Scholar]
  92. 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/Vmax method of Schmidt (1968). The mathematics behind 1/Vmax is the histogram, the simplest density estimator used mainly for rapid visualization of results in one or two dimensions. The 1/Vmax 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

ϕ ̂ ( z , L ) = n ( Z 2 − Z 1 ) f ̂ ( x , y ∣ h 1 , h 2 ) ( z − Z 1 ) ( Z 2 − z ) Ω dV dz , $$ \begin{aligned} \hat{\phi }(z, L)=\frac{n\left(Z_{2}-Z_{1}\right) \hat{f}\left(x, { y} \mid h_{1}, h_{2}\right)}{\left(z-Z_{1}\right)\left(Z_{2}-z\right) \Omega \frac{d V}{d z}}, \end{aligned} $$(A.1)

where [Z1, Z2] 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 f ̂ ( x , y | h 1 , h 2 ) $ \hat{f}(x,\mathit{y}|h_{1},h_{2}) $ is the density of (x, y), which corresponds to the (z, L) pair in the KDE parameter space. The calculation of f ̂ $ \hat{f} $ can be seen in Equation (9) in Yuan et al. (2022). Its value depends on two bandwidth parameters, h1 and h2. 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.

thumbnail Fig. A.1.

Corner plot showing the posterior probability distributions of the bandwidths (h1, h2) 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 best-fit value is the median.

All Tables

Table 1.

Values of AIC and BIC for Models A–C.

In the text
Table 2.

Best-fit parameters for Models A–C.

In the text

All Figures

thumbnail Fig. 1.

Redshift distribution (top) and the scatter plot (bottom) of our SFG sample. The red dashed curve indicates the flux limit line flim 1.4 GHz(z).
In the text
thumbnail 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
thumbnail 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.
In the text
thumbnail Fig. 4.

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

Similar to Fig. 3, but for Model C.
In the text
thumbnail 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 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
thumbnail 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).
In the text
thumbnail Fig. 8.

L⋆ × e2(z) (left panel) and ϕ⋆ × e1(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
thumbnail 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
thumbnail Fig. A.1.

Corner plot showing the posterior probability distributions of the bandwidths (h1, h2) 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 best-fit value is the median.
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.