The VLA-COSMOS 3 GHz Large Project: The cosmic star formation history since z~5

TL;DR

This study uses deep 3 GHz radio observations from VLA-COSMOS to trace the evolution of star formation rates and the cosmic star formation history up to redshift 5, revealing a peak at z~2-3 and highlighting the significant contribution of dust-obscured galaxies.

Contribution

It provides the largest radio-selected sample up to z~5 and models the evolution of the radio luminosity function and star formation rate density over cosmic time.

Findings

Star formation rate density peaks between z=2 and 3.

Ultraluminous infrared galaxies contribute up to 25% of the SFRD at high redshift.

UV-based SFR estimates may underestimate true rates by 15-20%.

Abstract

We make use of the deep Karl G. Jansky Very Large Array (VLA) COSMOS radio observations at 3 GHz to infer radio luminosity functions of star-forming galaxies up to redshifts of z~5 based on approximately 6000 detections with reliable optical counterparts. This is currently the largest radio-selected sample available out to z~5 across an area of 2 square degrees with a sensitivity of rms=2.3 ujy/beam. By fixing the faint and bright end shape of the radio luminosity function to the local values, we find a strong redshift trend that can be fitted with a pure luminosity evolution L~(1+z)^{(3.16 +- 0.2)-(0.32 +- 0.07) z}. We estimate star formation rates (SFRs) from our radio luminosities using an infrared (IR)-radio correlation that is redshift dependent. By integrating the parametric fits of the evolved luminosity function we calculate the cosmic SFR density (SFRD) history since z~5. Our data suggest that the SFRD history peaks between 2<z<3 and that the ultraluminous infrared galaxies (ULIRGs; 100 Msol/yr<SFR<1000 Msol/yr) contribute up to ~25% to the total SFRD in the same redshift range. Hyperluminous infrared galaxies (HyLIRGs; SFR>1000 Msol/yr) contribute an additional <2% in the entire observed redshift range. We find evidence of a potential underestimation of SFRD based on ultraviolet (UV) rest-frame observations of Lyman break galaxies (LBGs) at high redshifts (z>4) on the order of 15-20%, owing to appreciable star formation in highly dust-obscured galaxies, which might remain undetected in such UV observations.

Tables

Table 1. Table 1: Luminosity functions of star-forming galaxies obtained with the V max subscript 𝑉 max V_{\text{max}}  method.

Redshift	$\log \left(\frac{L_{1.4\text{GHz}}}{\text{W}{\text{Hz}}^{-1}}\right)$	$\log \left(\frac{\Phi }{{\text{Mpc}}^{-3}{\text{dex}}^{-1}}\right)$	N
	21.77${}_{}{}^{+0.23}{}_{-1.1}{}^{}$	-2.85${}_{}{}^{+0.094}{}_{-0.077}{}^{}$	189
	22.15${}_{}{}^{+0.18}{}_{-0.15}{}^{}$	-2.88${}_{}{}^{+0.031}{}_{-0.029}{}^{}$	217
	22.46${}_{}{}^{+0.19}{}_{-0.14}{}^{}$	-3.12${}_{}{}^{+0.037}{}_{-0.034}{}^{}$	149
	22.77${}_{}{}^{+0.20}{}_{-0.12}{}^{}$	-3.55${}_{}{}^{+0.063}{}_{-0.055}{}^{}$	56
	23.09${}_{}{}^{+0.21}{}_{-0.12}{}^{}$	-4.05${}_{}{}^{+0.11}{}_{-0.090}{}^{}$	19
	23.34${}_{}{}^{+0.28}{}_{-0.048}{}^{}$	-4.63${}_{}{}^{+0.25}{}_{-0.22}{}^{}$	5
	22.29${}_{}{}^{+0.11}{}_{-0.31}{}^{}$	-2.97${}_{}{}^{+0.052}{}_{-0.046}{}^{}$	142
	22.54${}_{}{}^{+0.13}{}_{-0.14}{}^{}$	-3.19${}_{}{}^{+0.036}{}_{-0.033}{}^{}$	160
	22.80${}_{}{}^{+0.15}{}_{-0.12}{}^{}$	-3.33${}_{}{}^{+0.038}{}_{-0.035}{}^{}$	143
	23.04${}_{}{}^{+0.18}{}_{-0.090}{}^{}$	-3.67${}_{}{}^{+0.059}{}_{-0.052}{}^{}$	65
	23.31${}_{}{}^{+0.18}{}_{-0.096}{}^{}$	-4.32${}_{}{}^{+0.12}{}_{-0.097}{}^{}$	16
	23.68${}_{}{}^{+0.081}{}_{-0.19}{}^{}$	-5.05${}_{}{}^{+0.34}{}_{-0.30}{}^{}$	3
	22.61${}_{}{}^{+0.080}{}_{-0.23}{}^{}$	-2.89${}_{}{}^{+0.091}{}_{-0.075}{}^{}$	179
	22.84${}_{}{}^{+0.15}{}_{-0.15}{}^{}$	-3.13${}_{}{}^{+0.027}{}_{-0.025}{}^{}$	283
	23.11${}_{}{}^{+0.17}{}_{-0.12}{}^{}$	-3.47${}_{}{}^{+0.035}{}_{-0.033}{}^{}$	165
	23.40${}_{}{}^{+0.17}{}_{-0.12}{}^{}$	-3.99${}_{}{}^{+0.066}{}_{-0.057}{}^{}$	51
	23.71${}_{}{}^{+0.16}{}_{-0.13}{}^{}$	-4.68${}_{}{}^{+0.16}{}_{-0.11}{}^{}$	11
	24.06${}_{}{}^{+0.10}{}_{-0.19}{}^{}$	-5.43${}_{}{}^{+0.45}{}_{-0.37}{}^{}$	2
	22.85${}_{}{}^{+0.074}{}_{-0.17}{}^{}$	-3.01${}_{}{}^{+0.046}{}_{-0.041}{}^{}$	172
	23.05${}_{}{}^{+0.13}{}_{-0.13}{}^{}$	-3.13${}_{}{}^{+0.025}{}_{-0.024}{}^{}$	312
	23.30${}_{}{}^{+0.14}{}_{-0.12}{}^{}$	-3.45${}_{}{}^{+0.032}{}_{-0.030}{}^{}$	198
	23.54${}_{}{}^{+0.16}{}_{-0.099}{}^{}$	-3.85${}_{}{}^{+0.051}{}_{-0.046}{}^{}$	82
	23.81${}_{}{}^{+0.15}{}_{-0.11}{}^{}$	-4.31${}_{}{}^{+0.088}{}_{-0.073}{}^{}$	30
	24.11${}_{}{}^{+0.11}{}_{-0.15}{}^{}$	-4.89${}_{}{}^{+0.18}{}_{-0.17}{}^{}$	8
	23.10${}_{}{}^{+0.081}{}_{-0.21}{}^{}$	-3.19${}_{}{}^{+0.052}{}_{-0.046}{}^{}$	216
	23.31${}_{}{}^{+0.15}{}_{-0.12}{}^{}$	-3.42${}_{}{}^{+0.025}{}_{-0.024}{}^{}$	321
	23.57${}_{}{}^{+0.16}{}_{-0.12}{}^{}$	-3.86${}_{}{}^{+0.036}{}_{-0.034}{}^{}$	156
	23.84${}_{}{}^{+0.16}{}_{-0.11}{}^{}$	-4.15${}_{}{}^{+0.052}{}_{-0.046}{}^{}$	81
	24.06${}_{}{}^{+0.22}{}_{-0.051}{}^{}$	-4.74${}_{}{}^{+0.10}{}_{-0.084}{}^{}$	22
	24.38${}_{}{}^{+0.17}{}_{-0.10}{}^{}$	-5.25${}_{}{}^{+0.20}{}_{-0.19}{}^{}$	7
	23.32${}_{}{}^{+0.070}{}_{-0.16}{}^{}$	-3.36${}_{}{}^{+0.043}{}_{-0.039}{}^{}$	156
	23.53${}_{}{}^{+0.18}{}_{-0.14}{}^{}$	-3.55${}_{}{}^{+0.025}{}_{-0.024}{}^{}$	323
	23.81${}_{}{}^{+0.21}{}_{-0.10}{}^{}$	-4.10${}_{}{}^{+0.041}{}_{-0.037}{}^{}$	126
	24.15${}_{}{}^{+0.19}{}_{-0.12}{}^{}$	-4.53${}_{}{}^{+0.068}{}_{-0.059}{}^{}$	48
	24.39${}_{}{}^{+0.26}{}_{-0.053}{}^{}$	-5.30${}_{}{}^{+0.18}{}_{-0.17}{}^{}$	8
	24.82${}_{}{}^{+0.14}{}_{-0.17}{}^{}$	-5.94${}_{}{}^{+0.45}{}_{-0.37}{}^{}$	2

Equations

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

$$L_{{\nu }_1}=\frac{4\pi D_L^2(z)}{(1+z{)}^{1+\alpha }}{\left(\frac{{\nu }_1}{{\nu }_2}\right)}^{\alpha }{S}_{{\nu }_2}.$$

$$\Phi (L,z)=\frac{1}{\Delta \log L}\sum _{i=1}^N\frac{1}{V_{\text{max},i}},$$

$$V_{\text{max},i}=\sum _{z=z_{\text{min}}}^{z_{\text{max}}}[V(z+\Delta z)-V(z)]C(z),$$

$$C(z)=\frac{A_{\text{obs}}}{41253\text{leavevmode}{\text{deg}}^2}\times C_{\text{radio}}[S_{3\text{leavevmode}\text{GHz}}(z)]\times C_{\text{opt}}(z),$$

$${\sigma }_{\Phi }(L,z)=\frac{1}{\Delta \log L}\sqrt{\sum _{i=1}^N\frac{1}{V_{\text{max},i}^2}}.$$

$${\Phi }_0(L)={\Phi }_{\star }{\left(\frac{L}{L_{\star }}\right)}^{1-\alpha }\exp \left[-\frac{1}{2{\sigma }^2}{\log }^2\left(1+\frac{L}{L_{\star }}\right)\right],$$

$$\Phi (L,z,{\alpha }_{\text{D}},{\alpha }_{\text{L}})=(1+z{)}^{{\alpha }_{\text{D}}}{\Phi }_0\left(\frac{L}{(1+z{)}^{{\alpha }_{\text{L}}}}\right),$$

$$\Phi (L,z,{\alpha }_{\text{L}},{\beta }_L)={\Phi }_0\left[\frac{L}{(1+z{)}^{{\alpha }_{\text{L}}+z{\beta }_L}}\right],$$

$$\frac{\text{SFR}}{{\text{M}}_{\odot }{\text{yr}}^{-1}}=4.5\times 10^{-37}\frac{L_{\text{TIR}}}{\text{W}},$$

$$q_{\text{TIR}}=\log \left(\frac{L_{\text{TIR}}}{3.75\times 10^{12}\text{W}}\right)-\log \left(\frac{L_{\text{1.4leavevmodenobreak GHz}}}{\text{W}{\text{Hz}}^{-1}}\right).$$

$$q_{\text{TIR}}(z)=(2.78\pm 0.02)\times (1+z{)}^{-0.14\pm 0.01}.$$

$$\frac{\text{SFR}}{{\text{M}}_{\odot }{\text{yr}}^{-1}}=f_{\text{IMF}}\times 10^{-24}10^{q_{\text{TIR}}(z)}\frac{L_{\text{1.4leavevmodenobreak GHz}}}{\text{W}{\text{Hz}}^{-1}},$$

$$\text{SFRD}={\int }_{L_{\text{min}}}^{L_{\text{max}}}\Phi (L,z,{\alpha }_{\text{D}},{\alpha }_{\text{L}})\times \text{SFR}(L)\text{d}\log L.$$

$$\log \frac{L_{\text{1.4leavevmodenobreak GHz}}}{\text{W}{\text{Hz}}^{-1}}=16.556-0.4(M_{1600,AB}-A_{UV})-q_{\text{TIR}}(z),$$

Full text

1 Introduction

One of the best methods to follow the buildup of stellar mass through cosmic times relies on inferring the cosmic star formation rate density (SFRD) history (for a review, see Madau & Dickinson, 2014). A consensus is achieved regarding recent history, where an exponential decline in SFRD by one order of magnitude from redshift $z\sim 2$ to the present day is inferred (e.g., Madau et al., 1996; Haarsma et al., 2000; Hopkins et al., 2006). On the other hand, with an increasing number of ultra-deep surveys the detection threshold is continually being pushed to higher redshifts (up to $z\sim 10$) slowly reaching the epoch of reionization (e.g., Bouwens et al., 2014a, 2015). The light of the early galaxies is a major factor in the process of reionization (e.g., Bouwens, 2016), and so accurate SFRD measurements are needed to better understand this epoch.

Although the wealth of observations has increased dramatically in the last decade, we still do not understand the core mechanism that governs star formation rate (SFR) histories of individual galaxies. This is because of our inability to actually follow these galaxies throughout their evolution. We observe galaxy populations at different cosmic epochs and try to link them in a consistent way. A picture has emerged from this method in which blue star-forming (SF) galaxies evolve into red quiescent galaxies through ways of quenching, such as rapid gas reservoir depletion after major merger interactions or active galactic nuclei (AGN) feedback (e.g., Bell et al., 2004; Schawinski et al., 2014). On the other hand, Bouché et al. (2010) presented a quenching-free model based on the cosmological decrease of accretion rates with time, which is able to reproduce the observed SFRD. Another model has also been proposed that uses simple mathematical lognormal forms for SFRD and individual SFR history to reproduce a wide range of observed relations (e.g., Gladders et al., 2013; Abramson et al., 2016). When the SFRD history is estimated with sufficient precision it can be used to further constrain semianalytical models of galaxy evolution, thereby deepening our understanding of the underlying physics.

Different SFR tracers can be used over the full electromagnetic spectrum, each with its own benefits and shortcomings (e.g., Kennicutt, 1998). The most direct tracer measures ultraviolet (UV) light from young massive stars and can be linked with the amount of star formation in the galaxy (e.g., Buat et al., 1989). The rest-frame UV emission is redshifted to optical and infrared (IR) wavelengths for the most distant galaxies; this enables the usage of very sensitive instruments, such as the Hubble Space Telescope (HST), to probe this epoch (e.g., Finkelstein et al., 2015). Currently, the SFRD in the earliest cosmic times (age of the universe less than 1 Gyr) is constrained almost exclusively with these kinds of observations (see also Behroozi et al., 2013). However, when measuring the rest-frame UV emission one must correct for dust extinction, which drastically diminishes the UV light. Well-constrained attenuation curves are needed to correct for this effect (e.g., Bouwens et al., 2009).

When dust grains absorb UV light they re-emit it at IR wavelengths. Therefore, far-IR and sub-mm traces SFR best when the dust content is high, yielding a large optical depth. These observations can suffer from poor resolution and source blending, although this was mitigated with observations with the Herschel Space Observatory. Current observations allow IR surveys to constrain the dust content and SFRs up to redshift $z<4$ (e.g., Caputi et al., 2005; Rodighiero et al., 2010; Reddy et al., 2012; Gruppioni et al., 2013). Ultraviolet and IR observations can be combined to obtain a more robust hybrid SFR estimator (e.g., Wuyts et al., 2011; Boquien et al., 2016). With the high-resolution sub-mm window opened by the Atacama Large Millimeter/submillimeter Array (ALMA), these wavelengths can be used to probe dusty submillimeter galaxies (SMGs) and their high star formation rates (e.g., Swinbank et al., 2014; Dunlop et al., 2017).

When massive stars undergo supernova explosions, the expanding remnants can accelerate the cosmic ray electrons and give rise to synchrotron radiation, which dominates the radio emission at rest-frame frequencies of $<30$ GHz. The observed nonthermal radio emission offers a dust-unbiased view at sub-arcsecond resolution of star formation processes inside the galaxy, and thus eliminates obscuration, while the high resolution assists counterpart matching (e.g., Seymour et al., 2008; Smolčić et al., 2009). However, it relies heavily on multiwavelength data to provide galaxy redshift and classification due to the featureless shape of the radio spectrum (e.g Condon, 1984). Furthermore, the SFR calibration for radio luminosities is based on the empirical IR-radio correlation to link nonthermal emission with thermal emission (e.g., Helou et al., 1985; Yun et al., 2001; Bell, 2003). This correlation continues to be valid across more than five orders of magnitudes in luminosities and holds at least up to redshift of $z<2$, albeit with some redshift evolution (e.g., Sargent et al., 2010; Magnelli et al., 2015), and it is likely to be valid even at earlier times up to $z\lesssim 5$ (Delhaize et al., 2017).

From observations and evolutionary models, we know that SF galaxies dominate the faint end of the radio counts (e.g., Condon, 1984; Gruppioni et al., 2003; Smolčić et al., 2008; de Zotti et al., 2010; Padovani, 2011; Smolčić et al., 2017a) and have strongly evolving luminosity functions (see also Rowan-Robinson et al., 1993), therefore deep surveys are needed to probe this population at early cosmic epochs. However, deep surveys have to sacrifice area in order to be feasible, which makes them more susceptible to cosmic over- and underdensities. This cosmic variance can have a strong redshift-dependent impact to any counting statistic employed (e.g., Moster et al., 2011).

The Cosmological Evolution Survey (COSMOS) 2 deg2 field (Scoville et al., 2007) is therefore well suited for our studies due to its large area, which should minimize cosmic variance, and excellent multiwavelength coverage, which allows for a precise photometric redshift determination. With the new Karl G. Jansky Very Large Array (VLA) observations obtained for the VLA-COSMOS 3 GHz Large Project (Smolčić et al., 2017b), the deepest radio survey to date given the area, we can probe the dust-unbiased SFRD up to redshift of $z\sim 5$ with $\sim$6 000 detections of SF galaxies. Our radio data best traces high-mass ($M_{\star}>10^{10}\leavevmode\nobreak\ \text{M}_{\odot}$) and highly SF galaxies (SFR $>100\leavevmode\nobreak\ \text{M}_{\odot}\,\text{yr}^{-1}$), which would also be classified as ultraluminous infrared galaxies (ULIRGs; $L_{\text{TIR, 8-1000\leavevmode\nobreak\ }\mu\text{m}}>10^{12}\leavevmode\nobreak\ \text{L}_{\odot}$, see Sanders & Mirabel 1996). At high redshift, we can also constrain even brighter hyperluminous infrared galaxy (HyLIRG; $L_{\text{TIR, 8-1000\leavevmode\nobreak\ }\mu\text{m}}>10^{13}\leavevmode\nobreak\ \text{L}_{\odot}$) populations, which have SFRs that are higher than 1 000 $\text{M}_{\odot}\,\text{yr}^{-1}$. To derive the total SFRD history of the entire radio population in the entire observed redshift range we must rely on extrapolations to lower luminosities below the sensitivity limit.

The paper is organized as follows. In Sect. 2 we briefly describe the data and selection methods used, which by itself is a topic of an accompanying paper (Smolčić et al., 2017a). We present methods of constructing luminosity functions and modeling their evolution through cosmic time and our results in Sect. 3. The calibration used to to derive SFR from radio luminosities is explained in Sect. 4 along with the cosmic SFRD history estimated from our data. We compare our results to the literature in Sect. 5. Discussion of possible systematics are given in Sect. 6. We finally summarize our findings in Sect. 7.

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

2 Data and star-forming galaxy sample

The sample of galaxies used in this work is radio selected with ancillary data from the rich multiwavelength coverage of COSMOS, which enables precise determination of redshifts and spectral energy distributions (SEDs).

2.1 Radio data

The radio data were obtained with 384 hours of VLA A+C array observations in the S band (2 GHz bandwidth centered around 3 GHz) within the VLA-COSMOS 3 GHz Large Project survey. Details of the observational setup, calibration, imaging, and source extraction can be found in Smolčić et al. (2017b). Briefly, 192 pointings were used to obtain a map of the COSMOS 2 square degrees with a uniform $rms$ noise equal to 2.3 $\mu\text{Jy}\,\text{beam}^{-1}$ and an angular resolution of $0\aas@@fstack{\prime\prime}75$. Imaging was performed using the multiscale multifrequency synthesis (Rau & Cornwell, 2011) to ensure good deconvolution of both unresolved and extended sources using the entire available 2 GHz bandwidth at once. Self-calibration of pointings containing brighter sources was performed to improve the fidelity of the image. A catalog of source components with a signal to noise (S/N) greater than 5 was extracted using the Blobcat software (Hales et al., 2012), which relies on a flood fill algorithm to detect contiguous blobs of emission. After visual inspection of multicomponent sources, a final catalog of 10 830 radio sources was assembled, spanning the entire observed area of 2.6 square degrees (approximately 10 000 radio sources across the central COSMOS two square degrees). The astrometric accuracy is $0\aas@@fstack{\prime\prime}01$ at the bright end and around $0\aas@@fstack{\prime\prime}1$ for the faintest sources.

2.2 Optical and near-infrared counterparts

We use the auxiliary data from more than 30 bands in the optical, near-infrared (NIR), and near ultraviolet (NUV) available in the COSMOS field from UltraVISTA DR2, Subaru/Hyper-Suprime-Cam, and SPLASH Spitzer legacy program collected in the COSMOS2015 catalog (Laigle et al., 2016). The catalog contains $\sim$800 000 sources with reliable photometry across an area of 1.77 deg2 free of stellar contamination. Photometric redshifts were computed for all sources by SED fitting using the LePhare code (Arnouts et al., 1999; Ilbert et al., 2006) following methods described in Ilbert et al. (2013).

The counterpart matching method is fully described in Smolčić et al. (2017a), their Section 3, and briefly summarized below. Owing to high sub-arcsecond resolution of both optical and radio data, and the fact that the radio emission is usually linked with massive bright galaxies, a nearest-neighbor counterpart matching scheme was adopted in combination with a false match probability assignment using a well-constructed background model. Optical-NIR counterparts were assigned to radio sources within a 0$\aas@@fstack{\prime\prime}$8 searching radius if they were deemed reliable. Estimates of false match probabilities were drawn from simulations using a background model that takes the $m_{3.6\leavevmode\nobreak\ \mu\text{m}}$ magnitude distribution of radio counterparts into account. It was designed to consider the optical blocking effect, i.e., missing fainter optical-NIR sources in the COSMOS2015 catalog due to nearby presence of a bright radio counterpart. Given these choices, the percentage of spurious matches in the entire radio sample are negligible. Approximately 11% of radio sources out of 8 696 positioned in the unmasked optical-NIR area was not assigned a COSMOS2015 counterpart. Half of those have $S/N<6$ in the radio source catalog making them likely candidates for spurious sources. The false detection probability for radio sources can reach up to 24% for sources with $5<S/N<5.1$, and a total of $\sim 3$% of sources in the radio catalog can be considered spurious (see also Smolčić et al. 2017b). If additional optical-NIR counterpart candidates are considered from the $i$-band selected catalog (Capak et al., 2007) and the Spitzer/IRAC111Infrared Array Camera catalog (Sanders et al., 2007), then $7.6\%$ of radio sources would remain without a counterpart in the same unmasked area (see also Smolčić et al., 2017a). We limit the optical-NIR counterpart matching to the COSMOS2015 catalog for better consistency with work by Delvecchio et al. (2017) and Delhaize et al. (2017). By taking the fraction of spurious sources into account, we have an average $\sim 8\%$ incompleteness in our counterpart sample. The total number of radio sources with assigned COSMOS2015 counterparts used throughout this paper is 7 729.

We use spectroscopic redshifts from the internal COSMOS catalog (M. Salvato et al. in prep.) available for 35% of our radio sources. These redshifts were used only if the spectra were flagged as reliable and 90% of those are located at $z<1.5$. Photometric redshifts were used for the remainder of the sample. We estimate the accuracy of photometric redshifts of our radio sample by comparing them to the above mentioned spectroscopic catalog and find a median $\Delta z/(1+z_{s})=0.01$ at all redshifts, and a 4% catastrophic failure rate, defined as $\Delta z/(1+z_{s})>0.15$. At redshifts $z>1.5$ we find a slightly larger median $\Delta z/(1+z_{s})=0.04$ and a catastrophic failure rate of 12%. Laigle et al. (2016) report the photometric redshift normalized median absolute deviation of the entire COSMOS2015 catalog in $3<z<6$ to be $\sigma_{\Delta z(1+z_{s})}=0.021$ with a catastrophic failure rate of 13.2% (see their Table 5).

2.3 Removing galaxies dominated by AGN in the radio band

We are interested in measuring the amount of star formation in galaxies from radio observations, disregarding whether the galaxy is an AGN host or not. We are therefore not interested in removing all AGN host galaxies from our sample, but only those that show clear evidence of radio emission dominated by an AGN. Unlike IR observations where the photometry can be used to trace a dusty torus in AGN (e.g., Donley et al., 2012), radio emission linked to star formation and AGN cannot be disentangled without assuming some correlation with emission at other wavelengths.

In order to quantify AGN contribution in each galaxy of our sample, Delvecchio et al. (2017) in their Section 3 performed a three-component SED fit using the sed3fit code222Publicly available at http://cosmos.astro.caltech.edu/page/other-tools (Berta et al., 2013). These fits were performed on the COSMOS2015 photometry using the best redshifts available and take the energy balance between the UV light absorbed by the dust and re-emitted in the IR into account, following the approach adopted in Magphys (da Cunha et al., 2008). In addition, an AGN component including the continuum disk and the dusty torus emission was added to the fit, seeking a best-fit solution via $\chi^{2}$ minimization (Berta et al., 2013). From these fits Delvecchio et al. (2017) estimated IR luminosity arising only from the star formation processes and calculated the IR-based (8-1000 $\mu$m) star formation rate SFR${}_{\text{IR}}$ using the Kennicutt (1998) relation and the Chabrier IMF. These SFR${}_{\text{IR}}$ values should be correlated with radio luminosities $L_{\text{1.4\leavevmode\nobreak\ GHz}}$ given the existence of the IR-radio correlation. In order to quantify this correlation and find outliers, these investigators construct histograms of $r=\log(L_{\text{1.4\leavevmode\nobreak\ GHz}}/\text{SFR}_{\text{IR}})$ in different redshift bins and fit a Gaussian distribution to the histogram (see Delvecchio et al. 2017, Section 4.2). The distributions of $r$ peak at higher values with increasing redshift, and in each redshift bin they are skewed toward higher values, corresponding to higher radio luminosities. Delvecchio et al. (2017) define radio-excess sources when $r$ deviates more than $3\sigma$ from the obtained peak of the distribution as a function of redshift, i.e.,

[TABLE]

Such a cut enables us to discriminate 1 814 (23%) sources dominated by AGN emission in the radio. If we assume that the peak of the above $r$ distribution at some redshift corresponds to the ideal correlation between the $L_{\text{1.4\leavevmode\nobreak\ GHz}}$ and the SFR${}_{\text{IR}}$, and all values above the peak are due to the increasing AGN contribution, then we estimate that the above cut corresponds to at least 80% of the radio emission due to the radio AGN component. The choice of this cutoff is somewhat arbitrary and was chosen as a conservative limit by Delvecchio et al. (2017) to minimize contamination of their AGN sample by SF galaxies. Radio emission of galaxies below this $3\sigma$ threshold might still be partly contaminated by AGN emission, but not likely dominated by it. Possible biases of this selection criterion are further discussed in Sect. 6.1.

We consider 5 915 radio sources without radio excess as our main SF galaxy sample. The redshift distribution of this final SF sample as well as radio-excess sources that were removed are shown in the upper panel of Fig. 1.

3 Radio luminosity function of star-forming galaxies

Radio luminosity functions (LFs) at different cosmic epochs are used to measure the evolution of radio sources, while also providing constraints on galaxy evolution models. We first discuss methods of determining the LF from our detections, we then show how the data can be fitted with an analytical form and, finally, we present LFs for our SF galaxies up to redshift of $z\sim 5$.

3.1 Estimating the luminosity function from the data

Throughout this work we assume that radio sources exhibit a radio spectrum described as a simple power law $S_{\nu}\propto\nu^{\alpha}$, where $S_{\nu}$ is a monochromatic flux density at frequency $\nu$ and $\alpha$ is the spectral index. This leads to the standard radio K correction of $K(z)=(1+z)^{-(1+\alpha)}$. The final expression for the rest-frame radio luminosity $L_{\nu_{1}}$ at frequency $\nu_{1}$ derived from the observed flux density $S_{\nu_{2}}$ at frequency $\nu_{2}$, redshift $z$, and luminosity distance $D_{L}$ is, therefore,

[TABLE]

Luminosities calculated at the rest-frame 1.4 GHz as a function of redshift are shown in the bottom panel of Fig. 1. This frequency is chosen to simplify a comparison of our results with the literature, where 1.4 GHz observations are more common. For about $\sim$25% of the sources we were able to derive the spectral index between 1.4 GHz (Schinnerer et al., 2010) and 3 GHz, while for the remaining sources we assumed the standard $\alpha=-0.7$, which is a valid median value for SF galaxies to be expected for shock-accelerated cosmic ray electrons.

To compute the density of sources and subsequently the LF at different cosmic times (i.e., redshift bins), we employed the $V_{\text{max}}$ method (Schmidt, 1968). This method uses the maximum observable volume of each source, while satisfying all selection criteria; it is not dependent on the shape of the LF and therefore reduces the sample and selection biases. The LF $\Phi(L,z)$ gives the number of radio sources in a comoving volume per logarithm of luminosity and is obtained as

[TABLE]

where $V_{\text{max}}$,i is the maximum observable volume of the $i$-th source, $\Delta\log L$ is the width of the luminosity bin, and the sum goes over each source $i$ in a given redshift and luminosity bin. To take into account different effects and biases, such as a luminosity limited sample or nonuniform noise in the radio map, which may lead to an incompleteness of the sample, we employed a very general form for calculating the maximum observable volume $V_{\text{max}}$, i.e.,

[TABLE]

where the sum starts at the beginning of a chosen redshift bin and adds together comoving volume spherical shells $\Delta V=V(z+\Delta z)-V(z)$ in small redshift steps $\Delta z=0.01$ until the end of the redshift bin is reached. The parameter $C(z)$ is the redshift-dependent geometrical and statistical correction factor that takes the observed area and sensitivity limit into account and further mitigates some of the other already mentioned completeness issues

[TABLE]

where ${A_{\text{obs}}}=1.77\leavevmode\nobreak\ \text{deg}^{2}$ corresponds to the effective unflagged area observed in the optical to NIR wavelengths, $C_{\text{radio}}$ is the completeness of the radio catalog as a function of the flux density $S_{3\leavevmode\nobreak\ \text{GHz}}$, and $C_{\text{opt}}$ is the completeness owing to radio sources without assigned optical-NIR counterpart. The area observed in the radio encompasses the entire area observed in the optical-NIR and does not have flagged (cropped) regions, therefore NIR observations set the limit for the observed area. The $C_{\text{radio}}$ factor depends on the redshift because a source with a given intrinsic luminosity changes its apparent flux density between $z_{\text{min}}$ and $z_{\text{max}}$ in $V_{\text{max}}$ calculations (see Eq. (4)).

Completeness corrections are shown and tabulated in Smolčić et al. (2017b); see their Fig. 16 and Table 2, respectively. We also show these corrections in the top panel of Fig. 2 in this work. These completeness corrections are linearly interpolated between tabulated values for any flux density below $S_{\text{3\leavevmode\nobreak\ GHz}}<100\leavevmode\nobreak\ \mu$Jy. Simulations were not optimized to probe sources with flux densities above $100\leavevmode\nobreak\ \mu$Jy and we assume a 100% completeness for such sources. Completeness corrections are based on Monte Carlo simulations of mock-source generation and extraction and take into account the nonuniform $rms$, proper derivation of flux densities for low S/N sources and the resolution bias (out-resolving and losing extended low-surface brightness radio emission). The last part, which was modeled by assuming the distribution of radio sizes, follows some functional form of flux densities, which reproduces the observed data (for details see Smolčić et al. 2017b). These corrections are a function of radio flux density only, meaning that all other physical properties are averaged out. For example, the presence of more resolved and extended (and also low-surface brightness) galaxies at lower redshift as a result of their closer proximity to us, which may introduce a redshift-dependent bias.

In Sect. 2.2 we mentioned that 11% of our radio sources were not assigned a counterpart. In the bottom panel of Fig. 2 we show the completeness of our radio catalog $C_{\text{opt}}$ due to matching with the COSMOS2015 catalog as a function of redshift. It was obtained by considering an additional optical catalog selected in the $i$ band (Capak et al., 2007). The counterpart completeness was calculated as $1-N_{i-\text{band}}/N_{i-\text{band or COSMOS2015}}$, where $N_{i-\text{band}}$ is the number of new counterparts assigned only to $i$-band selected sources (1% of the total radio sample) and $N_{i-\text{band or COSMOS2015}}$ is the number of counterparts assigned to either optical catalog. As already mentioned, we only use the COSMOS2015 catalog for consistency reasons. Our LFs results are perfectly consistent between themselves whether the actual $i$-band counterparts or the $C_{\text{opt}}$ correction curve is used. As shown in the bottom panel of Fig. 2, the counterpart sample is complete up to $z\sim 1.5$, and $\sim 90\%$ complete at $z\sim 5.5$ The addition of this completeness correction, while also considering the 3% of spurious radio sources, leaves 7% of real radio sources unaccounted for. A reliable redshift distribution is not available for these of sources, if the source follows a strong redshift trend, for example, all of these sources are located at $z>3$, it still might bias our high-redshift LFs low.

There might exist a small number of galaxies with a high radio flux density and a faint optical-NIR magnitude whose $V_{\text{max}}$ would be determined by the optical-NIR limit. Since the optical-NIR catalog was selected on the $\chi^{2}$ image, it does not have a well-defined magnitude limit and therefore we cannot apply a more precise correction on $V_{\text{max}}$. This may bias our high redshift LF low as the true $V_{\text{max}}$ would be smaller than what we have used for such sources. However, we do not expect a significant effect since there are only $\sim 10$ sources in our sample that have $K_{s}$ AB magnitudes fainter than 24.5 where the completeness of the optical-NIR catalog becomes an issue (see Laigle et al., 2016).

The $rms$ error estimate of the LF in each redshift and luminosity bin is calculated as in Marshall (1985) by weighting each galaxy by its contribution to the sum

[TABLE]

However, if there are ten sources or fewer in a luminosity bin we used the tabulated upper and lower 84% confidence intervals from Gehrels (1986). These intervals correspond to Gaussian $1\sigma$ errors so that $\sigma_{\Phi}=\Phi\times\sigma_{N}/N$, where $\sigma_{N}$ is the small-number Poissonian statistical asymmetrical error on the measured number of sources. We do not add photometric uncertainties into the error budget, but the redshift bins are chosen to be large enough to mitigate possible problems of sources falling into wrong bins. An additional contribution to the total error budget may arise from the imperfect radio SED (see also Sect. 6.4).

3.2 Local radio luminosity function and its evolution

Radio LFs of SF galaxies are usually described by four parameter analytical forms such as the power-law plus lognormal distribution from Saunders et al. (1990)

[TABLE]

where the $L_{\star}$ parameter describes the position of the turnover of the distribution, $\Phi_{\star}$ is used for the normalization, $\alpha$ and $\sigma$ are used to fit the faint and bright ends of the distribution, respectively.

Our deep COSMOS radio observations are best suited to study the high-luminosity end of the LF, especially at higher redshifts ($z>1$), where our data do not sample the faint end of the LF, but instead cover a large observed volume. If we are interested in the total amount of light emitted from SF galaxies at any redshift we must assume the shape of the LF that is not constrained by our data. These luminosities can be probed with wide and shallow low-resolution radio surveys of the local universe, such as the NVSS333National Radio Astronomy Observatory (NRAO) VLA Sky Survey. (Condon et al., 1998). There are a number of works related to the calculation of the local radio LF of SF galaxies (e.g., Condon, 1989; Condon et al., 2002; Sadler et al., 2002; Best et al., 2005; Mauch & Sadler, 2007) and they are all broadly consistent in the luminosity range of $21.5<\log L_{\text{1.4\leavevmode\nobreak\ GHz}}[\text{W}\,\text{Hz}^{-1}]<23.5$.

To obtain the local luminosity function that is used throughout this work, we performed a fit on combined volume densities from Condon et al. (2002); Best et al. (2005); Mauch & Sadler (2007) using the form given in Eq. (7). By combining the data from both wide and deep surveys we can properly constrain both the faint and bright end of the local LF. The data and fit are shown in Fig. 3. Obtained best-fit parameters are $\Phi_{\star}=3.55\times 10^{-3}\leavevmode\nobreak\ \text{Mpc}^{-3}\text{dex}^{-1}$, $L_{\star}=1.85\times 10^{21}\leavevmode\nobreak\ \text{W}\,\text{Hz}^{-1}$, and $\alpha=1.22$, $\sigma=0.63$.

We assume that the shape of the LF remains unchanged at all observed cosmic times and allows only the position of the turnover and normalization to change with redshift. This corresponds to the translation of the local LF in the $\log L-\log\Phi$ plane (Condon, 1984) and can be divided into pure luminosity evolution (horizontal shift) and pure density evolution (vertical shift). Using a simple one parameter power law for each of these evolution cases the form of the redshift evolved LF is

[TABLE]

where $\alpha_{\text{D}}$ and $\alpha_{\text{L}}$ represent pure density and pure luminosity evolution parameters, respectively, and $\Phi_{0}(L)$ is given in Eq. (7). Since our data are more sensitive to the most luminous star-forming galaxy population above the knee of the LF, these two evolution parameters may become degenerate preventing a precise estimate of the knee location, especially at higher redshifts. This choice for the LF evolution is chosen for its simplicity given that our data constrain the bright end of the LF the best (see also Sect. 6.2). In reality, all four parameters may change with redshift.

3.3 Radio luminosity functions across cosmic times

The procedure of binning sources into luminosities inherently introduces some biases due to averaging and the chosen bin sizes. To minimize possible completeness issues at the faint luminosity end within a redshift bin, all sources with luminosities below the observational luminosity limit (corresponding to $5\sigma=11.5\leavevmode\nobreak\ \mu\text{Jy}$ at 3 GHz) at $z_{\text{max}}$ of the redshift bin were put into single luminosity bin. All sources above this limit were distributed into equally wide luminosity bins spanning the observed luminosity range. The actual luminosity value of each point that we report is the median of all galaxies in a given luminosity bin, while horizontal error bars show the bin width. For easier comparison with work in the literature, we calculated each LF using the 1.4 GHz rest-frame luminosity obtained from our observations at 3 GHz. Our LFs from the $V_{\text{max}}$ method are shown in Fig. 4 as black circles and are also tabulated in Table The VLA-COSMOS 3 GHz Large Project: Cosmic star formation history since $z\sim 5$. Our data have small Poissonian error bars due to the relatively large number of sources in each bin and errors do not reflect all possible systematic effects, such as the unknown radio K correction, the error on the completeness, or the sample contamination. A comparison with LFs derived by other authors at different wavelengths is discussed in Sect. 5.

The data points were then fitted with an evolved local LF given in Eq. (8). The redshift that enters this expression is the median redshift of all galaxies in a given redshift bin. A $\chi^{2}$ minimization was performed to obtain the best fit $\alpha_{\text{L}}$ and $\alpha_{\text{D}}$ parameters. Since the LF may have asymmetric errors in sparsely populated bins due to small number Poissonian statistics, an average value of the upper and lower errors on the LF was taken for the $\chi^{2}$ computation. These parameters are degenerate when either the faint or bright ends are not sampled well, therefore a pure luminosity evolution ($\alpha_{\text{D}}=0$) was computed as well. Errors on the parameters were estimated from the $\chi^{2}$ statistics following Avni (1976). We derived the formal $1\sigma$ errors by projecting onto each parameter axis ($\alpha_{\text{L}}$ and $\alpha_{\text{D}}$) the 68% confidence contour around the minimum $\chi^{2}$.

The best-fit evolution parameters obtained are shown in Fig. 5 as a function of redshift.

3.4 A simple evolution model

In order to create a single continuous model for the evolution of star-forming LF across the entire observed cosmic time, we simultaneously fit all LF points in all redshift bins with a two parameter pure luminosity evolution described as

[TABLE]

where $\Phi_{0}$ is the local LF from Eq. (7), and we allow for an additional redshift-dependent change of the power law parametrized with $\beta_{L}$. This form follows single redshift bin fits well (see Fig. 5) and is chosen for its simplicity. Significant density evolution cannot be properly constrained by our observations, which is why we do not attempt it here. From the $\chi^{2}$ minimization fit we obtain the following values for parameters: $\alpha_{\text{L}}=3.16\pm 0.2$ and $\beta_{L}=-0.32\pm 0.07$.

4 Cosmic star formation rate density history

4.1 From radio luminosity to star formation rate

Radio emission can be used as an extinction-free tracer of star formation rate when linked to other more direct (thermal) tracers such as the IR light. The first assumption is that the UV photons of massive young stars are absorbed by the dust and re-emitted in the IR so that the total IR emission of a galaxy correlates well with its SFR, which is valid for optically thick galaxies. The conversion factor relies on estimating mass from light and was calibrated by Kennicutt (1998) assuming the Salpeter (1955) initial mass function ($\text{d}N/\text{d}M\propto M^{-2.35}$) from 0.1 to 100 $\text{M}_{\odot}$ and is given by

[TABLE]

where $L_{\text{TIR}}$ contains the total integrated IR luminosity of a galaxy between 8-1000 $\mu$m. This IMF produces more low-mass stars than are supported by observations that favor a turnover below $1\leavevmode\nobreak\ \text{M}_{\odot}$. Since low-mass stars do not contribute significantly to the total light of the galaxy, only the mass-to-light ratio is changed when the Chabrier (2003) IMF is adopted instead. This leads to a decrease in SFR by a factor of 1.7 (see also Pozzetti et al., 2007) because of there are fewer low-mass stars created. The calibration itself usually leads to a 0.3 dex scatter on a galaxy basis (see also Condon, 1992; Murphy et al., 2011; Kennicutt & Evans, 2012).

Radio observations can trace recent star formation of galaxies, and can trace these observations on timescales of up to 100 Myr (Condon, 1992). Estimation of a galaxy SFR from the radio data relies heavily on the observational IR-radio correlation that is known to span at least five orders of magnitudes (Helou et al., 1985; Yun et al., 2001). The IR-radio correlation links the radio luminosity to the TIR luminosity via the $q_{\text{TIR}}$ parameter defined as

[TABLE]

Usually, the $q_{\text{TIR}}$ parameter is taken to be a constant value derived for local galaxies. However, recent works suggest that the $q_{\text{TIR}}$ value might change with redshift (e.g., Sargent et al., 2010; Magnelli et al., 2015). In this paper, we used methods from Delhaize et al. (2017) who constrained the median of the $q_{\text{TIR}}$ as a function of redshift using a doubly censored survival analysis for a joint 3 GHz radio and IR-selected sample. They find a decrease of $q_{\text{TIR}}$ with redshift that can be parameterized with a simple power law. To be self-consistent, we ran the survival analysis on the same SF sample as utilized in this work, while also taking into account limits for IR-detected galaxies without a $5\sigma$ significant radio emission, because in their paper Delhaize et al. (2017) originally used a different sample selection criteria for excluding AGN. The obtained evolution of the IR-radio correlation for our sample can be written as

[TABLE]

The main idea behind the IR-radio correlation is that a linear relation exists between radio and IR luminosities for SF galaxies. There is a possibility that the decreasing q_{\text{TIR}}$$(z) actually mimics some complexities of the radio SED at high redshifts such as varying degrees of free-free contribution and inverse Compton losses. Inverse Compton losses off the cosmic microwave background (CMB) lead to suppression of nonthermal radio continuum emission, which would in turn increase $q_{\text{TIR}}$ with redshift, but the opposite trend was observed by Delhaize et al. (2017). In the case of a more complicated radio SED, a simple power-law K correction is not a valid assumption anymore. However, the use of a redshift-dependent $q_{\text{TIR}}$(z) parameter when calculating SFR should account for these intrinsic observational limitations under the assumption of a linear IR-radio correlation as explained in more detail in Sect. 6.3.

Finally, the expression that converts radio luminosity into SFR obtained from the steps described above can be written as

[TABLE]

where $f_{\text{IMF}}=1$ for a Chabrier IMF and $f_{\text{IMF}}=1.7$ for a Salpeter IMF.

4.2 Star formation rate density across cosmic times

Integrating below the LF first multiplied by the luminosity we can obtain the total 1.4 GHz radio luminosity density ($\text{W}\leavevmode\nobreak\ \text{Hz}^{-1}\leavevmode\nobreak\ \text{Mpc}^{-3}$) in a chosen redshift range. Similarly, if the radio luminosity is converted to SFR as given in Eq. (13) before integration, the integral yields the SFRD of a given epoch

[TABLE]

We numerically integrated the above expression by taking the analytical form of the LF in each redshift bin and using the best-fit evolution parameters shown in Fig. 5. The integral was calculated in different luminosity ranges, which are listed below (results shown in Fig. 6 and also listed in Table 1 Luminosity functions of star-forming galaxies obtained with the $V_{\text{max}}$ method.):

Entire luminosity range: This formally means setting $L_{\text{min}}=0$ and $L_{\text{max}}\rightarrow+\infty$. The integral converges and the major contribution to the SFRD arises from galaxies with luminosities around the turnover of the LF. The entire radio emission is recovered and if the LF shape and evolution is well constrained the SFRD estimate will be as well (within the SFR calibration errors). This is not the case at higher redshifts ($z>2.5$), where only the bright end of the LF is observed, therefore extrapolation to the faint end can be substantial (see Fig. 4). 2. 2.

Data constrained limits: $L_{\text{min}}$ and $L_{\text{max}}$ correspond to the lowest and highest value of the observed luminosity function. By choosing integration limits that correspond to the actual data range, any bias due to LF extrapolation toward higher or lower luminosities is removed. The shape of the local LF also does not affect this result within the fitting errors. Numbers obtained from this integration range are a very conservative lower limit on the SFRD. 3. 3.

ULIRGs: Limits that correspond to galaxies with IR luminosity of $10^{12}\leavevmode\nobreak\ \text{L}_{\odot}<L_{\text{TIR}}<10^{13}\leavevmode\nobreak\ \text{L}_{\odot}$ trace ULIRGs. The radio luminosity limits were obtained using an evolving $q_{\text{TIR}}$ parameter from Eq. (12). The integral with such a range traces SFRD of galaxies that form stars very efficiently (SFR $100-1000\leavevmode\nobreak\ \text{M}_{\odot}\,\text{yr}^{-1}$) while also being well constrained by our observational data in $0.5<z<3$ range (see also the red dotted vertical line in Fig. 4). 4. 4.

HyLIRGs: Similarly, by integrating over galaxies with radio luminosities that translate into $L_{\text{TIR}}>10^{13}\leavevmode\nobreak\ \text{L}_{\odot}$, we trace HyLIRGs that have extreme star formation, namely SFR $>1000\leavevmode\nobreak\ \text{M}_{\odot}\,\text{yr}^{-1}$.

Our errors are inferred from the LF fitting parameters uncertainties and added in quadrature with $q_{\text{TIR}}(z)$ parameter errors and do not represent the entire error budget due to LF extrapolations.

5 Comparison with the literature

To check the robustness of our LF and SFRD results presented in Figs. 4 and 6 and also to create a consistent multiwavelength picture, we compare them with work in the literature derived at radio, IR, UV, and sub-mm wavelengths. All SFR estimates were rescaled to a Chabrier IMF where necessary.

5.1 Radio and IR luminosity functions

In Fig. 4 we compare our results with the radio LFs by Smolčić et al. (2009), which are based on the VLA-COSMOS 1.4 GHz survey (Schinnerer et al., 2007). These investigators constructed LFs up to $z<1.3$ using a sample of 340 galaxies classified as star forming using optical rest-frame colors. The increase in sensitivity of the VLA-COSMOS 3 GHz survey along with a different selection method yielded $\sim 10$ times more detections of star-forming galaxies in the same redshift range. The two results generally agree with each other, although our LFs are slightly higher, most likely because of different selection criteria adopted.

We additionally plot LFs from the IR surveys to compare the validity of our results at higher redshifts as well. If the IR-radio correlation is linear, both IR and radio LFs should follow each other well. To convert the total IR (TIR) to radio luminosity function, the redshift dependent IR-radio correlation parameter $q_{\text{TIR}}$ described in Eq. (11) and (12) is used. We show the LFs by Magnelli et al. (2013) derived up to $z<2.3$ from Herschel observations of GOODS-N/S deep and GOODS-S ultra-deep fields444The Great Observatories Origins Deep Survey (GOODS).. We also show the LFs by Gruppioni et al. (2013), which were computed up to $z<4.2$ and are based on Herschel PEP/HerMES555The Herschel Guaranteed Time Observation (GTO) PACS Evolutionary Probe (PEP) Survey; The Herschel Multi tiered Extragalactic Survey (HerMES). observations. To take into account the fact that the redshift bin ranges do not necessarily coincide, we evolved the LFs of other authors using the evolution parameters they reported in their paper, if our median redshift value fell inside their redshift range. Small systematic offsets may arise when the mean redshift does not correspond to the median redshift. Our data agree well with these surveys both at low and intermediate redshifts. However, at redshift $z>2$ our LFs are systematically slightly lower than those based on IR. Some of this offset might be attributed to a higher percentage of AGN in the IR selected sample at these redshifts (Gruppioni et al., 2013). They may constitute half of the sample above $z>2.5$. However, since we start from a differently selected radio sample, exclude AGN identified with radio excess compared to IR emission, and we must rely on the redshift evolving q_{\text{TIR}}$$(z), it makes the direct comparison difficult. If a constant $q_{\text{TIR}}=2.64$ (Bell, 2003) is used for the conversion, instead of an evolving one, our observed radio LFs would actually be higher than implied by the observed IR-based LFs at high redshifts.

5.2 UV luminosity functions

Our radio data are good tracers of highly star-forming and dusty galaxies (ULIRGs and HyLIRGs), but lack the sensitivity to probe fainter sources at high redshifts. We make use of the work performed by Bouwens et al. (2015) in an attempt to constrain the faint end of the luminosity functions of SF galaxies with actual detections and to simultaneously test their dust corrections. Bouwens et al. (2015) utilize HST observations of more than ten deep and wide surveys covering $\sim 1000$ arcmin2 to derive the rest-frame UV LFs between $4<z<10$ using a sample of more than 10 000 Lyman break galaxies (LBGs). The rest-frame UV light correlates strongly with the SFR, unless the galaxy is very dusty. Therefore we can make a broad comparison with our SF galaxy sample.

The SFR calibrations from Kennicutt (1998) are self-consistent, meaning that all tracers (radio, IR, and UV) should provide the same SFR estimate, thus enabling the link between radio and UV luminosities via the SFR. Although this correlation likely has a large scatter when applied to a specific galaxy, if used on larger samples as a statistical conversion factor, it should allow the conversion of UV magnitudes into radio luminosities. The conversion is needed to compare LFs at these two different wavelengths. The expression for this conversion using the Kennicutt (1998) calibration, Chabrier IMF, and the redshift-dependent IR-radio correlation from Eq. (12) is

[TABLE]

where $M_{1600,AB}$ is the rest-frame UV absolute magnitude reported in Bouwens et al. (2015) and $A_{UV}$ is the extinction needed to calculate the dust-corrected magnitudes. The dust extinction, obtained from the IRX-$\beta$ relationship (correlation between the ratio of FIR to UV fluxes with the UV spectral slope $\beta$; see Meurer et al. 1999), is given in the form of $A_{UV}=4.43-1.99\beta$ and tabulated as a function of UV magnitudes in Bouwens et al. (2014b).

The bottom panels in Fig. 4 show dust-uncorrected (green right triangles) and dust-corrected (dark green circles) LFs from Bouwens et al. (2015) for their $z\sim 4$ and $z\sim 5$ dropouts. Several results can be noted in the plots:

Significant dust corrections are needed at high luminosities and the rest-frame UV cannot be used to detect dustier ULIRGs and HyLIRGs, which are easily observed in the radio at high redshifts. Our radio detections, also available across more than $6\times$ larger area, can therefore provide an independent test for these dust corrections. 2. 2.

The bright end of the UV LF is lower than our radio LF, with the discrepancy being larger at $z\sim 4$ than at $z\sim 5$. Our result is also broadly consistent with the result of Heinis et al. (2013) in which they perform stacking of Herschel images of UV selected galaxies at $z\sim 1.5$. These authors found that the IR luminosity of bright galaxies ($L_{\text{IR}}>10^{11}\leavevmode\nobreak\ \text{L}_{\odot}$) obtained via stacking can recover as low as 52-65% of the total SFRD derived from the IR-selected samples. 3. 3.

Even if we disregard the dust correction, the density of faint galaxies two decades in luminosity below our detection limit is very high. The dust-corrected UV LFs, are in broad agreement with our pure-luminosity fit extrapolation more than two decades below the lowest observed radio luminosity at $z\sim 5$. These arguments can be used to rule out most of the gray-shaded area in the bottom two panels of Fig. 4 arising from significant negative density evolution (see also lower panel of Fig. 5). Rest-frame UV observations are in favor of higher densities of galaxies than what would be obtained if a turnover in the radio LF is introduced immediately below the faintest observed radio luminosity. 4. 4.

There is a discrepancy between our pure luminosity evolution model and the UV LF at the faintest observed end. Since our radio data cannot constrain such low luminosities, pure luminosity evolution is probably not the best possible model for extrapolating our observed LFs below our detection limits. Indeed, a continuous steepening with redshift of the faint end slope of the LF has recently been proposed (e.g., Parsa et al., 2016).

We further discuss the UV data in the context of our radio estimates of SFRD in Sect. 5.6.

5.3 Total SFRD estimates

Throughout all panels in Fig. 6 we show our total SFRD derived by integrating the pure luminosity evolved LF in individual redshift bins as black filled circles. In panel A we compare our SFRD with the curve from the review by Madau & Dickinson (2014) where the fit was performed on a collection of previously published UV and IR data (red line). Below $z<2$ our data agree well with this compilation of data, but show a turnover at higher redshift ($z\sim 2.5$) with a shallower decline yielding up to 2-3 times larger SFRD at $z\gtrsim 4$. We also plot a slightly different Behroozi et al. (2013) fit to the data compilation in the same panel.

If we allow for both luminosity and density evolution there is a degeneracy of parameters leading to large uncertainties in the total SFRD estimate; the gray shaded area in panel A of Fig. 6 is obtained with fit parameters taken from the $1\sigma$ significant region in $\alpha_{\text{D}}$ and $\alpha_{\text{L}}$ parameter space. We do not fit a pure density evolution because it would increase the normalization of the LF to very high densities. The SFRD estimates would consequently be significantly higher, making our data completely inconsistent with other works in the literature at intermediate and high redshifts. In the same panel we also show very strict lower limits constrained by the data with blue triangles that demonstrate the amount of extrapolation needed to obtain the total SFRD. Although the extrapolation is be significant, especially at higher redshifts, we note that the UV LFs support the need for such large extrapolations.

5.4 Comparison with previous radio SFRD

In panel B of Fig. 6 we show two radio estimates based on the VLA-COSMOS 20 cm survey (Schinnerer et al., 2007, 2010). Smolčić et al. (2009) calculated the SFRD by integrating the pure luminosity evolution fit of a local LF taken from Sadler et al. (2002) and their results are shown with blue squares. Also, these estimates were obtained in the COSMOS field and therefore they represent a good consistency check at low redshift. A different approach was taken by Karim et al. (2011) who performed stacking on mass selected galaxies (shown as orange diamonds). They obtained a monotonous rise in the SFRD up to $z\sim 3$. Although the field is the same, their method of estimating SFRD is significantly different from ours since it depends on stacking individually undetected sources. Our estimates are slightly lower than theirs, with the difference increasing with redshift. This offset is primarily due to a different IR-radio correlation used. They adopt a calibration from Bell (2003), which yields higher SFRD at higher redshifts.

5.5 Comparison with IR SFRD

In panel C of Fig. 6 the pink shaded area shows the $1\sigma$ uncertainty for the SFRD derived from the integrated total IR LF by Gruppioni et al. (2013). This LF has a rising trend up to $z\sim 1.1$ and then flattens out. The highest redshift estimate should be considered as a lower limit because the PEP selection might miss high-z sources. Our results are in broad agreement with theirs. Discrepancies at some redshifts might be attributed to different sample selection, since we are excluding AGN host galaxies classified as such only in the radio (see Sect. 6.1). Additionally, the agreement between the IR and the radio SFRD is better at $z>2$ than at $z\sim 1.5$, while the opposite is true when comparing IR to radio LF (see Fig. 4). The reason for this effect is that the normalization of the Gruppioni et al. (2013) IR LF is slightly higher than ours. However, because of the significant negative density evolution and the unchanging faint end slope, this higher normalization is being progressively compensated in their SFRD integral by decreased densities of the faint end. Differing contributions of the faint and the bright end to the total SFRD as a function of redshift lead to apparent agreement between IR and radio SFRD estimates, even though the actual observed LFs do not match perfectly.

We also show the results from the recent work by Rowan-Robinson et al. (2016) as purple plus signs in the plot. Using SED fitting on $\sim 3000$ Herschel sources from 20.3 square degrees of the sky they derive an IR-based SFRD since $z\sim 6$. Even though their result has large uncertainties, the finding supports a much flatter SFRD trend at high redshifts. It is still consistent, however, with our findings within the error bars. In the same panel we additionally show SFRD results of an extended halo model estimated by Planck Collaboration et al. (2014) from the measured power spectra of the cosmic IR background anisotropies as orange-red shaded area. They report several possible reasons for such high values of SFRD at $z>2$, and it is important to note that all measurements rely on some form of extrapolation toward fainter galaxies. These results might therefore be considered as upper limits.

5.6 UV addition to SFRD estimates

In panel D of Fig. 6 we show the recent rest-frame UV estimates from Bouwens et al. (2015) as dark-green filled squares (dust-corrected) and green open circles (dust-uncorrected). The SFRD is scaled to Chabrier IMF and Kennicutt (1998) calibration. Ultraviolet observations like these are well suited to study the early universe owing to the ability to probe exceptionally high redshifts $z\sim 10$, as also reviewed in Madau & Dickinson (2014).

To simultaneously model both the faint and bright ends of the SF LFs at high redshifts in an attempt to better constrain the SFRD of that epoch we use the dust-corrected UV LFs from Bouwens et al. (2015) along with our own radio LFs and perform a fit on the combined data as explained below. The UV dust corrections are more severe at high luminosities and the LBG selection criteria can easily miss the most massive and dusty galaxies with significant SFRs. On the other hand, the radio emission is an excellent tracer of such SF galaxies. Therefore, we disregard the three most luminous UV LF points at redshifts $z\sim 4$ and $z\sim 5$ and fit an analytical form given in Eq. (7) to the remaining UV points combined with all of our radio LFs at the same redshift. The combined data span more than four decades in luminosities. Our results are shown in Fig. 7, where we show the SFR on the x-axis instead of the usual luminosity. Ultraviolet luminosities were scaled to SFR according to Kennicutt (1998) relation, while our radio luminosities were scaled using the redshift-dependent $q_{\text{TIR}}$ given in Eq. (12). It is not our intention to obtain the best SF LF at these redshifts, but rather to calculate an estimate of the missed SFRD in the LBG sample from the radio perspective. Still, for completeness we report here the best-fit parameters obtained. They are $\Phi_{\star}=9.35\times 10^{-3}\leavevmode\nobreak\ \text{Mpc}^{-3}\text{dex}^{-1}$, $L_{\star}=1.81\times 10^{22}\leavevmode\nobreak\ \text{W}\,\text{Hz}^{-1}$, $\alpha=1.62$, and $\sigma=0.83$ at $z\sim 4$ and $\Phi_{\star}=1.23\times 10^{-3}\leavevmode\nobreak\ \text{Mpc}^{-3}\text{dex}^{-1}$, $L_{\star}=1.26\times 10^{23}\leavevmode\nobreak\ \text{W}\,\text{Hz}^{-1}$, $\alpha=1.76$, and $\sigma=0.67$ at $z\sim 5$.

The SFRD integral of the best LF fit of the combined dust-corrected UV and radio data is 0.08 dex higher at $z\sim 4$ and 0.06 dex higher at $z\sim 5$ than the values obtained from the UV data alone in the same luminosity range. These integrated values are also plotted as blue diamonds in panel D of Fig. 6. Assuming the dust corrections calculated by Bouwens et al. (2014b) are correct and start to become significant only at higher luminosities and SFRs (for details, see also Wang & Heckman, 1996), this suggests a 15-20% underestimation of highly obscured SFR estimated from the rest-frame UV observations. Since our radio LFs are slightly lower than the IR LFs at $z\sim 4$ (see Fig. 4), this underestimation could be considered a lower limit. Also, our pure radio SFRD estimate is likely underestimated at $z\sim 4$ due to a rather flat faint end slope, while at $z\sim 5$ it is actually higher than the combined UV plus radio estimate owing to a higher normalization of the pure evolution fit.

Our radio LFs are in very good agreement with the work done by Mancuso et al. (2016). These authors used the continuity equation approach with the main sequence star formation timescales to conclude that the number density of SF galaxies at high redshifts ($z\lesssim 7$) cannot be reliably estimated from the UV-data alone, even when corrected for dust extinction. Their results also imply the existence of a high-redshift heavily dust-obscured galaxy population with SFRs larger than $100\leavevmode\nobreak\ \text{M}_{\odot}\,\text{yr}^{-1}$.

In their work, Burgarella et al. (2013) attempted to constrain the SFRD by taking into account dust obscuration using combined IR and UV LFs reported in Gruppioni et al. (2013) and Cucciati et al. (2012), respectively. We show their results as magenta crosses in panel D of Fig. 6. It is interesting to note a good agreement in SFRD at $z\simeq 4$ between substantially different approaches such as the pure UV-based data, IR plus UV data, and the radio plus UV data. They are all consistent within $\sim$20%, but at the same time higher than previously reported SFRD fits (Madau & Dickinson, 2014; Behroozi et al., 2013). Work carried out by Dunlop et al. (2017) is another example that aims at a complete dust-obscured and unobscured (UV + FIR) SFRD census at high redshifts utilizing ALMA observations of the Hubble Ultra Deep Field (HUDF) at 1.3 mm. These investigators estimate UV contribution to the total SFR from evolving luminosity functions given in Parsa et al. (2016) and Bouwens et al. (2015). Dunlop et al. (2017) find SFRD (shown as red squares in panel D of Fig. 6) consistent with Behroozi et al. (2013) in the redshift range $2.5<z<4.5$. They also find a transition in the dominant SF population from dust obscured to dust unobscured at $z\gtrsim 4$. Given the wide distribution and uncertainties of calculated SFRD arising in insufficient knowledge of dust corrections, we believe that the inclusion of radio observations as a dust-unbiased tracer can help achieve a better consensus.

5.7 ULIRGs and HyLIRGs

In panel E of Fig. 6 we decompose our total SFRD to focus on galaxies that form stars efficiently (SFR $>100\leavevmode\nobreak\ \text{M}_{\odot}\,\text{yr}^{-1}$). Our SFRD estimates for ULIRGs and HyLIRGs are shown as purple asterisks and red crosses, respectively. As previously, the first consistency check is to compare our SFRD for ULIRGs with those estimated by Smolčić et al. (2009), shown as blue downward triangles. Our values are slightly higher than theirs, as was the case with the LFs, in the redshift range sampled by Smolčić et al. (2009). Ultraluminous IR galaxies are well constrained by our data in the redshift range $0.6<z<3.3$, with little to no extrapolation needed. The data show that ULIRGs contribute above 16% to the total SFRD at $z>1$ with a peak of $\sim$25% at a redshift of $z\sim 2.5$, while HyLIRGs contribute an additional $\lesssim 2\%$ in the entire observed range.

Caputi et al. (2007) inferred the bolometric IR ($5-1000\leavevmode\nobreak\ \mu$m) luminosity density up to $z\sim 2$ using Spitzer 24 $\mu$m selected galaxies in the GOODS fields. We show their SFRD results for ULIRGs only as green dot-dashed line. The agreement is good up to $z\sim 1$, but it gets worse at higher redshifts where their estimates are significantly higher than ours. Discrepancies may be caused by a different star-forming galaxy sample selection as mentioned in the previous section where we compared IR- and radio-based SFRD. Additionally, the luminosity integration limits needed to calculate contributions from ULIRGs are directly scaled (see Sect. 4.2) by the redshift-dependent $q_{\text{TIR}}$ parameter. The total scaling effect of the q_{\text{TIR}}$$(z) on the SFRD integral is further discussed in 6.3.

Additionally, in the same panel, we show SFRD based on $\sim 100$ ALMA LESS (ALESS666LABOCA Extended Chandra Deep Field South (ECDFS) Submm Survey (LESS).) submillimeter galaxies ($S_{870\mu\text{m}}>1$ mJy) by Swinbank et al. (2014) as the magenta dashed line and $1\sigma$ errors as dotted lines. These SMGs are highly dust obscured and have large SFR. Since our lower integration limit for ULIRGs ($100\leavevmode\nobreak\ \text{M}_{\odot}\,\text{yr}^{-1}$) is slightly higher than theirs ($80\leavevmode\nobreak\ \text{M}_{\odot}\,\text{yr}^{-1}$), 0.11 dex should be added to our ULIRG SFRD values prior to comparison with Swinbank et al. (2014) results. There is a broad agreement within the error bars between these two results. However, there are some additional complications in comparing these results because their observations are less sensitive to hotter than average dust temperatures, and they report up to a factor of 2 uncertainty due to missing these ULIRGs. Therefore, their results represent lower limits. Also, both the Swinbank et al. (2014) results and our results rely on non-negligible extrapolations to fainter flux densities.

6 Potential biases

Here we summarize some critical assumptions and associated possible systematic effects on our results. While the biggest uncertainties arise from extrapolations, there are a number of additional redshift-dependent and independent effects that may scale our LFs and SFRD history.

6.1 AGN contamination

In this paper we adopted an IR-radio-based discrimination of galaxy populations since our goal was to estimate LFs of star-forming galaxies and the total SFRD history from the radio perspective. We assume that the IR is a good tracer of SFR in our radio detected galaxies and that SF galaxies follow a IR-radio correlation with some intrinsic scatter. Because the observed scatter is nonsymmetrical, i.e., there is a tail of sources with large radio fluxes compared to IR measurements, we conclude that AGN contribution to the radio emission is large in such galaxies. The radio-excess method described in Sect. 2.3 is therefore good at selecting galaxies dominated by AGN emission in the radio band. The $3\sigma$ cut given in Eq. (1) ensures that only $\sim$0.15% of removed galaxies are SF, giving us a high level of completeness of our SF sample. On the other hand, by counting sources below the $3\sigma$ radio-excess limit and above the best-fitting symmetric profile in all redshift bins (see Sect. 2.3) we estimate that the integrated radio emission can be contaminated by some AGN contribution for around 1 000 sources (17% of the sample), and this contribution is limited to a maximum of 80%. This potential AGN contribution is mitigated when calculating the SFRD integral by using a properly calibrated $q_{\text{TIR}}$ relation (see Sect. 6.3). When AGN enter the sample, they increase the density in the LF, but at the same time lower the $q_{\text{TIR}}$ parameter (see Delhaize et al., 2017). If a smaller, less than $3\sigma$, cut were used, then more and more SF galaxies would be removed from the sample trading completeness for purity.

In an attempt to obtain a clean SF sample, free of AGN hosts, we employed a different selection method explained in Smolčić et al. (2017a). We start from the full radio sample with optical-NIR counterparts. The first step in removing AGN includes the use of a cutoff in the X-ray full band (rest-frame $0.5-8$ keV) luminosity (see Szokoly et al., 2004). In the second step, a warm dusty torus signature around the supermassive black hole is found in the MIR using a cut in the four IRAC bands as prescribed in Donley et al. (2012). The third step uses SED fits with AGN templates (da Cunha et al., 2008; Berta et al., 2013) to exclude galaxies with significant AGN contribution (see also Delvecchio et al., 2017). These three criteria remove moderate-to-high radiative luminosity AGN from the sample. The next step uses rest-frame optical colors $M_{NUV}-M_{r}$ corrected for internal dust extinction to select red quiescent galaxies ($M_{NUV}-M_{r}>3.5$, Ilbert et al. 2010) that may host an AGN detected in the radio. If galaxies with such colors do not have a $5\sigma$ detection in the Herschel image, they are then classified as low-to-moderate radiative luminosity AGN and excluded from our sample. The remaining 4555 green and blue galaxies ($M_{NUV}-M_{r}\leq 3.5$) without a $3\sigma$ radio excess (see Sect. 2.3) are considered a clean sample of SF galaxies based on available AGN diagnostics. Since this sample does not account for the star formation in AGN hosts, it represents a conservative lower limit to SF LFs and SFRD. When the analysis is repeated with this clean SF sample, we find a median decrease of 0.12 dex in number densities across all observed epochs without significant redshift trends. For consistency, the $q_{\text{TIR}}$ was recalculated for the clean SF galaxy sample. It gives slightly higher values at all redshifts, but agrees with that given in Eq. (12) within $1.5\sigma$. The total SFRD integral median decrease is 0.035 dex, which is within the uncertainties of our nominal sample. The SFRD median decrease is not as significant as the LF median density decrease because we are still fitting the pure luminosity evolution to newly derived LFs, meaning that the faint end remains mostly unchanged; we also recalibrated the $q_{\text{TIR}}$ parameter to match the clean SF sample.

6.2 The choice of the local LF

The choice of the analytical shape of the LF can have a significant effect on the total SFRD due to extrapolation toward unobserved luminosities. A compilation of local LF is shown in the top panel Fig. 8. Where necessary, according to Ascasibar et al. (2002), LF was corrected for the change of cosmology by scaling $L(z)\propto d_{m}^{2}(z)\propto h^{-2}$ and $\Phi(z)\propto d_{m}^{-3}(z)\propto h^{3}$. The bottom panel of Fig. 8 shows the contribution to the SFRDs as a function of luminosity for the various LFs; although all the LFs show a peak at a similar radio luminosity, the positions of these peaks can differ by up to $\sim$0.3 dex. There is also no physical argument for the shape of the LF being fix and evolving in redshift by simple translation. However, our data cannot constrain the full luminosity range required to obtain the most significant bulk of the SFRD integral at all redshifts, so this way of extrapolation was chosen for its simplicity. For a more complex handling of the LF evolution, see for example Fotopoulou et al. (2016), where they used a Bayesian approach to model and constrain the shape of the AGN LF as a function of redshift.

6.3 IR-radio correlation

The most significant factor in our SFRD estimates is the $q_{\text{TIR}}$ parameter since it directly scales our integrated radio luminosities as a function of redshift. Throughout this work we used q_{\text{TIR}}$$(z) estimated on the same SF galaxy sample with the methods from Delhaize et al. (2017). The following are a few underlying assumptions when using such an evolving $q_{\text{TIR}}$(z) value:

The IR emission is an accurate tracer of SFR at all redshifts and radio emission originates mostly in SF processes. Extragalactic radio observations can properly trace emission from SF processes in a galaxy when cosmic ray electrons are not allowed to escape it. The escape scenario is possible for small sized galaxies with $L_{\text{1.4\leavevmode\nobreak\ GHz}}\lesssim 2\times 10^{21}\leavevmode\nobreak\ \text{W}\,\text{Hz}^{-1}$ (e.g., Bell, 2003), which is far below our observational limit. However, the nonthermal radio emission needs a proxy to derive the actual SFR value and the assumption is that the IR emission is a good proxy. 2. 2.

Infrared-radio correlation is linear, meaning that it can be represented as a single line with a slope of one in the log-log plot of radio and IR luminosities. 3. 3.

Radio spectrum is a simple power law in frequency. This is a widely used approximation and is often taken for granted because of insufficient radio data, however, it plays an important role, especially at high redshifts.

Within the framework of these assumptions it is correct to use an evolving $q_{\text{TIR}}$(z) when calculating the SFR of a galaxy from radio emission. Even if the second or the third assumption was not correct, for example, because of various free-free contributions in the radio spectrum or the luminosity dependence of the IR-radio correlation, which might cause a difference between the IR and radio LF evolution, the q_{\text{TIR}}$$(z) evolution takes these wrong assumptions into account and produces a correct SFR value on average because it was calibrated using both the radio and IR data.

To demonstrate the scaling effect of the $q_{\text{TIR}}$ parameter on our SFRDs we integrate our continuous simple evolution model from Sect. 3.4 and show the results with a blue line in panel F of Fig. 6, while the shaded area corresponds to the $1\sigma$ uncertainty owing to the errors on the fit parameters added in quadrature with the q_{\text{TIR}}$$(z) uncertainty. If we instead take the standard constant local value of $q_{\text{TIR}}=2.64$ from Bell (2003) and apply it to our simple LF evolution model, we would obtain three times larger SFRD estimates at $z\sim 4$ (see gray dotted line in the same panel). Observations however do not favor this choice. Another analysis of the IR-radio correlation was performed through stacking by Magnelli et al. (2015). They obtained $q_{\text{FIR}}(z)=2.35\times(1+z)^{-0.12}$. This relation can be scaled as $\log(L_{\text{FIR}})=\log(L_{\text{TIR}})-\log(2)$ to obtain the $q_{\text{TIR}}$(z) needed for our conversion, which is valid in terms of median statistics; see also Delhaize et al. (2017). The SFRD obtained from this expression is shown as a red dot-dashed line in the same panel and is similar to ours. To summarize, the trend in the cosmic SFRD history that we obtain from our simple LF evolution model is linked with the trend in the $q_{\text{TIR}}$ and it is important for this value to be well constrained at all observed redshifts.

6.4 Radio spectral indices

Regarding the accuracy of the computed rest-frame 1.4 GHz luminosity, the highest uncertainty, especially at high redshifts, lies in the insufficient knowledge of the radio K correction. For example, a rather large photometric error of $\Delta z=0.3$ would result in a 0.05 dex error in luminosity at $z\sim 5$. However, an uncertainty in spectral index of just $\Delta\alpha=0.1$ would produce an error of 0.1 dex in luminosity. It is known that a symmetric spread of the spectral indices around the mean value for SF galaxies is approximately $\sim 0.4$ (e.g., Kimball & Ivezić, 2008). Even though this spread would produce a significant uncertainty on the estimated $L_{\text{1.4\leavevmode\nobreak\ GHz}}$ of each single galaxy, these variations approximately cancel each other and give a valid average luminosity because of the symmetry of the distribution of spectral indices. It is widespread in the literature to assume a single spectral index for radio SED, where the usual values are $\alpha=-0.8$ or $-0.7$.

Approximately 75% of our radio sources were only detected at 3 GHz. The use of the measured spectral indices for the remaining 25% can introduce a small bias toward steeper spectra (therefore higher luminosities), since our survey is currently the deepest radio survey of the COSMOS field. For example, a point source at the limit of our sensitivity ($rms=2.3\leavevmode\nobreak\ \mu\text{Jy}\,\text{beam}^{-1}$) would have to have a spectral index steeper than -1.9 to be observed in the previous deep 1.4 GHz survey ($rms=10\leavevmode\nobreak\ \mu\text{Jy}\,\text{beam}^{-1}$; Schinnerer et al. 2010). The median spectral index of sources detected in both surveys is $\alpha=-0.85$.

To assess the impact of the used spectral indices on our results, we repeated the analysis two times: the first time with the standard $\alpha=-0.7$ and the second time with $\alpha=-0.8$ for all sources regardless of the observed radio spectrum. When a single and identical spectral index is used for each source, the pure luminosity evolution of the chosen analytical local function from Eq. 7 better fits (smaller $\chi^{2}$) the derived LFs at all luminosities and redshifts. Specifically, when $\alpha=-0.7$ is used, the best pure luminosity fit evolution remains essentially unchanged from that presented in Sect. 3.3, which is unsurprising given that 75% of the spectral indices remained unchanged as well. When $\alpha=-0.8$ is used, a stronger luminosity evolution is obtained, which is described by an increase of 0.16 in the evolution parameter $\alpha_{L}$ as given in Eq. (8) and previously shown in Fig. 5. This increase is still within the model uncertainties as given in Sect. 3.4. Before deriving the SFRD values, for consistency, the $q_{\text{TIR}}$ parameter was again recalculated using different spectral indices and obtained expressions are within the $1.5\sigma$ of the nominal values given in Eq. (12). Derived total SFRDs are within the uncertainties of the nominal sample in both cases, which strengthens the robustness of our results.

Finally, assuming a simple power-law radio spectra might be an overly simplistic approach given the unknown contribution of free-free (Bremsstrahlung) emission to the total SED. Additional deep radio observations at higher frequencies are needed to properly model the radio SED and mitigate this limitation.

7 Summary

We studied a radio-selected sample of star-forming galaxies from deep VLA-COSMOS 3 GHz observations (Smolčić et al., 2017b). Galaxy classifications were performed by relying heavily on the optical-NIR COSMOS2015 catalog (Laigle et al., 2016) and SED fits by Delvecchio et al. (2017). The final sample contains 5 915 galaxies, where the radio emission is not dominated by an AGN. Using this sample we derived radio LFs up to $z\sim 5$. By comparing them with LFs derived using IR and UV selected samples we checked their robustness and found that our radio LF can be very well described by a local LF with a power-law plus lognormal form evolved only in luminosity as ${L_{\text{1.4\leavevmode\nobreak\ GHz}}\propto(1+z)^{(3.16\pm 0.2)-(0.32\pm 0.07)z}}$. However, we do not observe the faint end of the LF at all redshifts to properly constrain a more complex evolution. The difference between radio and UV LFs suggests an underestimation of dust corrections obtained from UV slopes by Bouwens et al. (2014b). We converted our radio luminosities to SFR using a redshift-dependent IR-radio correlation where $q_{\text{TIR}}$ parameter decreases with increasing redshift (Delhaize et al., 2017). An accurate constraint on this parameter is the most important factor for estimating SFR from radio observations in the early universe. Our data suggest that the peak of the total SFRD history occurs at $2<z<3$. We find that the total SFRD estimates using only LBG galaxies (e.g., Bouwens et al., 2015), even if corrected for dust extinction, are still likely to miss up to 15-20% of SFR in highly obscured galaxies at $z\gtrsim 4$. By integrating LF fits in various luminosity limits we estimated SFRDs of the total SF sample and the subpopulations of the sample, such as ULIRGs and HyLIRGs. We find that ULIRGs contribute at maximum up to $\sim$25% of the total SFRD at $z\sim 2.5$, where this population of galaxies is well constrained by our data. Even though HyLIRGs can have very large SFRs (several 1000 $\text{M}_{\odot}\,\text{yr}^{-1}$), we find that they contribute less than 2% to the total SFRD at all redshifts owing to their low volume density.

Acknowledgements.

This research was funded by the European Unions Seventh Frame-work program under grant agreement 337595 (ERC Starting Grant, ’CoSMass’). ES warmly acknowledges support from the National Radio Astronomy Observatory for visits to the Array Operations Center in Socorro, NM, in conjunction with this project. AK acknowledges support by the Collaborative Research Council 956, sub-project A1, funded by the Deutsche Forschungsgemeinschaft (DFG). MB and PC acknowledge support from the PRIN-INAF 2014.