An Empirical Estimation of the Galaxy Radio Number Counts model at 3 GHz and the Effect of Gravitational Lensing

Exploring the galaxies evolution throughout cosmic time is the central topic of astrophysics. It dates back to 70’s that astrophysicists started to pay attention to various evolutionary aspects of galaxies (e.g., Larson & Tinsley, 1974, 1978). Then, after the extensive studies by Tinsley (1980) and coevals, the evolution has been regarded as one of the most fundamental properties of galaxies. Among a huge number of relevant studies, recent deep surveys in various wavelengths  (e.g., Smolčić et al., 2009; Smolčić et al., 2017b; Ueda et al., 2014; Padovani et al., 2015) and theoretical studies (e.g., Croton et al., 2006; Hopkins et al., 2008; Somerville et al., 2008; Mancuso et al., 2015; Caplar et al., 2015, 2018) have shown a drastic evolution of physical properties of galaxies through the cosmic time.

Especially, observations indicate that galaxies and super-massive black holes (SMBHs) at the center of galaxies have evolved together (e.g., Marconi & Hunt, 2003; Silverman et al., 2009; Kiuchi et al., 2009; Kormendy & Ho, 2013; Yang et al., 2017, 2018; Krajnović et al., 2018). This is referred to as the co-evolution of galaxies and SMBHs. These studies revealed that the star formation rate density (SFRD) have a peak at $z\sim 1\mbox{--}3$, sometimes called "the Cosmic Noon" (e.g., Takeuchi et al., 2005; Bell et al., 2007; Magnelli et al., 2011; Gruppioni et al., 2013; Burgarella et al., 2013; Madau & Dickinson, 2014; Bouwens et al., 2015; Matthews et al., 2021, among others). The accretion history onto SMBHs also has been constrained from observations of AGNs (e.g., Sołtan, 1982; Marconi et al., 2004; Netzer & Trakhtenbrot, 2007; Delvecchio et al., 2018). For such approach, a luminosity function (LF) serves as a fundamental statistical tool to extract physical information from observations (e.g., Marconi et al., 2004; Kriek et al., 2007; Silverman et al., 2008; Yuan et al., 2016; Yang et al., 2017).

The radio LF has been studied for many years based on the data mainly from the VLA (e.g., Machalski, 2000; Mauch & Sadler, 2006; Smolčić et al., 2009). These studies revealed that the local radio LFs consist of two distinct components. The faint end of LF is dominated by star forming galaxies (SFG) and the bright end is dominated by active galactic nuclei (AGN). The radio continuum emission from SFG is due to synchrotron emission from relativistic electron accelerated by the shock of supernova remnants. Therefore, the radio luminosity from SFGs can be interpreted as the tracer of star forming rate (SFR) of the galaxy (Condon, 2008; Murphy et al., 2011). On the other hand, the radio emission from AGNs is dominated by the synchrotron emission from Mpc scale radio lobe and associated jet from center of a galaxy. Hence, the radio emission from AGNs can be interpreted as the mass accretion rate onto a SMBH at the center of galaxies (e.g., Sołtan, 1982; Marconi et al., 2004; Netzer & Trakhtenbrot, 2007; Körding et al., 2008; Delvecchio et al., 2018). AGNs play an important role to understand the evolution of galaxies through radio mode AGN feedback (e.g., Croton et al., 2006; Kriek et al., 2007; Merloni & Heinz, 2008; Fabian, 2012).

Since the estimation of a LF requires redshift measurements, it is not immediately available from the completion of a certain survey. Observed number counts of galaxy provides the simplest statistics to constraint their evolution. It had been used to constraint cosmological parameters in earlier studies as the value is sensitive to the geometry considered. Takeuchi et al. (2000, 2001) investigated FIR source counts and pointed out that strong luminosity evolution is necessary to explain the IRAS source counts. In radio band, the situation is slightly different as the bright end of radio number count is dominated by AGNs, expecially referred as radio loud AGNs, while most of the faint end population comes from SFGs (e.g., Condon, 1984b; Planck Collaboration et al., 2011; Padovani et al., 2015; Vernstrom et al., 2016). (see also Padovani, 2016). Additionally, probability of deflection, or $P(D)$ analysis is also applied to achieve confusion-limited number counts of galaxies below survey sensitivity limit which is pioneered by Scheuer (1957). The $P(D)$ is the probability distribution of peak flux and consists of the convolution of system noise, residual sidelobes and so on. This analysis is widely applied in various observation frequency such as in infrared band (e.g., Takeuchi & Ishii, 2004; Berta et al., 2011) and in radio and submillimeter band (e.g., Condon et al., 2012; Vernstrom et al., 2014; González-López et al., 2020; Mauch et al., 2020). Especially, the radio number counts is investigated down to sub- $\mu{\rm Jy}$ in Vernstrom et al. (2014) and Mauch et al. (2020) with $P(D)$ analysis.

For the analysis of galaxy/AGN number counts, we have to take into account yet another important effect, gravitational lensing. It is well known that the galaxy number counts in submillimeter band shows prominent magnification effect of gravitational lensing. Some previous studies show that approximately $60~{}\%$ of galaxies are strongly lensed in the bright end (e.g., Lima et al., 2009; Béthermin et al., 2012; Negrello et al., 2017; Mancuso et al., 2017) in the submillimeter and millimeter band as a result of the negative $K$-correction because of the shape of their spectral energy distribution. This effect causes a significant bias when we estimate the evolution of galaxy from blind survey data. Furthermore, the development of very long base line interferometers (VLBI) such as square kilometre array (SKA) and the Karl G. Jansky Very Large Array (VLA) allow us to perform deep survey even in the radio band. Distant galaxies are more likely to be lensed due to the long path toward us. We therefore need to estimate the statistical effect of gravitational lensing on the galaxy number counts.

In this study, we estimated the radio LFs for $0.1<z<5.5$ at 3GHz based on the data obtained from the VLA-COSMOS 3GHz large project Smolčić et al. (2017a), and construct a galaxy evolution model that is able to reproduce observed number counts with assuming the pure luminosity evolution (PLE) model that only the parameters on luminosity vary with redshift. Afterwards, we discuss about properties future SKA and other panchromatic surveys based on resultant LF evolution model.

The structure of this paper is as follows. In Section 2, we outline the VLA-COSMOS data and the classification of the sample. In Section 3, we show the method to estimate the radio LFs. In section 4, we show the results of the estimated LFs and construct our evolution model From Section  5 to 6, we discuss about the effect of gravitational lensing on galaxy number counts and cosmic star formation rate density history based on obtain LF evolution for SFG. Finally, we summarize our results in section 7.

In this study, we adopt cosmological parameter with $H_{0}=70\ [{\rm km\ s^{-1}\ Mpc^{-1}}]$, $\Omega_{\Lambda}=0.7$, $\Omega_{M}=0.3$, $\sigma_{8}=0.8$, throughout this paper.

Some literature have already investigated the evolution of radio luminosity function in the COSMOS field (e.g., Smolčić et al., 2009; Smolčić et al., 2017c; Novak et al., 2017; Novak et al., 2018; Ceraj et al., 2018; van der Vlugt et al., 2021; Malefahlo et al., 2021) on several types of SFG and AGN up to $z\lesssim 6$. In this paper, we provide evolution of luminosity function up to $z\lesssim 5.5$ estimated from two independent methods. One is a non-parametric estimator called the $C^{-}$ method (Lynden-Bell, 1971) and the other one is parametric Bayesian estimator which is formulated in (Kelly et al., 2008).

We calculated the SFRD up to $z\sim 6$ based on the radio continuum luminosity density for SFGs with our evolution model. We made use of a tight correlation between radio and IR emission (e.g., Condon, 1992; Yun et al., 2001) to convert $L_{\rm 3GHz}$ to SFRD. The relation commonly parameterized with $q\mbox{-}$parameter introduced in Helou et al. (1985).

Some studies argued that this relation shows decrease with redshift increase (e.g., Magnelli et al., 2015; Calistro Rivera et al., 2017; Ocran et al., 2020) which indicates the radio luminosity increases with redshift relative to IR luminosity. We assumed redshift dependent total infrared to radio luminosity ratio as $q_{\rm TIR}(z)=(2.78\pm 0.02)\times(1+z)^{-0.14\pm 0.01}$ obtained by Delhaize et al. (2017) for VLA-COSMOS SFG sample. For comparison, we applied constant $q_{\rm TIR}$ as $2.64\pm 0.02$ derived in Bell et al. (2003). This relation reflects the mechanisms of radio emission (Yun et al., 2001; Bell et al., 2003). As for the conversion from IR luminosity to SFR, we adopt the relation based on Kennicutt, Jr. (1998).

where $\rho_{\rm IR}$ is IR luminosity density per unit comoving volume $[\rm Mpc^{3}]$ at redshift $z$ determined by integrating $\int^{L_{\rm max}}_{L_{\rm min}}L\phi(L,z)d\log L$. We set $(L_{\rm min},L_{\rm max})=(0,\infty)$ in the integration. While the integration range covers unrealistic values, the bulk of contribution to luminosity density comes from characteristic luminosity. Note that this assumption requires considerable extrapolation in LF at faint end especially in the high redshift as the data is limited to the bright end.

where $\beta_{\rm SFG}$ is the spectral index for SFG radio spectrum. The value of the coefficient $f_{\rm IMF}$ depends on IMF. In this study, we assumed Chabrier IMF Chabrier (2003), thus $f_{\rm IMF}=1$.

The resultant cosmic SFRD is displayed in Fig. 4 with blue solid line along with the observational results from other wavelength. The shaded area corresponds to $95\%$ confidence interval of the LF resultant parameters. The SFRDs obtained by (Smolčić et al., 2009; Karim et al., 2011; Novak et al., 2017) are derived through radio surveys. Smolčić et al. (2009) SFRD, red diamonds, is based on LF in $1.4~{}{\rm GHz}$ estimated from 340 radio selected SFG up to $z=1.3$ detected in the VLA-COSMOS $1.4~{}{\rm GHz}$ survey (Schinnerer et al., 2006a). Karim et al. (2011) constructed SFRD evolution model derived from stacked $1.4~{}{\rm GHz}$ continuum for $3.6~{}{\rm\mu m}$ selected SFG in the COSMOS field at $0.2<z<3$ which is plotted as green squares. While our result is baed on the same sample as Novak et al. (2018) (blue circles), their LF is based on $1/V_{\rm max}$ method. In order to compare our result with SFRD derived through LF in various wave band, we plot result obtained by Burgarella et al. (2013); Bouwens et al. (2015); Koprowski et al. (2017), yellow triangle, green stars and black dot-dashed line, respectively, and analytical form of the cosmic SFRD history obtained by Madau & Dickinson (2014) with red dashed line. We also compare the result from UV band. Bouwens et al. (2015) derived rest-frame dust corrected UV LF obtained by over $10^{4}$ Lyman break galaxies at $4<z<10$ in the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS) field identified by HST observation and Burgarella et al. (2013) combination of FIR (Herschel/PEP and Herschel/HerMES) and FUV LFs (VIMOS-VLT Deep Survey). Koprowski et al. (2017) estimated SFRD with SCUBA-2 and ALMA imaging data in sub-mm/mm up to $z\sim 5$ and adding the SFRD contribution derived by UV luminosity. Madau & Dickinson (2014) compiled IR and UV data from multiple literature and derived parametric SFRD evolution function.

The resultant monotonically growing SFRD up to $z\sim 2.5$ agrees well with the literature. Especially, the SFRD by Novak et al. (2018) has a good agreement with our result for all redshift range considered and confirmed that the LF from independent method gives consistent evolution. The difference between our SFRD and Karim et al. (2011) becomes larger at $z>2$. The difference comes from the fact that they applied constant $q_{\rm TIR}$ parameter derived by Bell et al. (2003) which makes the IR luminosity larger at high redshift. The brown dotted line indicates SFRD obtained assuming constant $q_{\rm TIR}$. While the result is has good agreement with Karim et al. (2011) and Burgarella et al. (2013) up to $z\sim 3$, the values at higher redshift is too high comparing to the literature. Adding to this, the model takes peak value at $z\sim 4$ which is too large. Based on this result, redshift dependent $q_{\rm TIR}$ is gives SFRD which is consistent with the result from observations in other wavelength. We also plot SFRD whose integral base is assumed to be $L_{\rm lim}(z,S_{\rm lim})$ instead of $L_{\rm lim}=0$ with light-blue dotted line as a lower limit of prediction. Note that there is still uncertainty in the faint end slope to be constrained through future observation.

Having constructed LF evolution model, we estimated the effect of gravitational lensing on galaxy number counts. The magnification $\mu$ due to gravitational lensing is given as the Jacobian between image coordinate and lens coordinate. Magnification profile $\mu(\theta)$ is calculated as below.

where $\kappa$ is convergence, $\gamma$ is shear and $\theta$ denotes the separate angle from the center of the dark matter halo. These values depend on dark matter halo properties such as the distance to lens object $z_{l}$, mass $M$ and ellipticity $\epsilon$, and the distance to source object. We therefore adopt the parametric form of function from Sheth & Tormen (1999) for galaxy halo mass function and the NFW profile (Navarro et al., 1997) as the radial density profile of each dark matter halo. We adopted the concentration parameter $c_{\rm vir}\equiv R_{\rm vir}/r_{s}$ in the NFW profile is given as a function of virial mass $M_{\rm vir}$ and redshift that is given by Bullock et al. (2001). We also applied the analytical form for $\kappa$ and $\gamma$ obtained by Takada & Jain (2003a) and Takada & Jain (2003b). Now we considered the ellipticity of DM halos $\epsilon$. It is well known that the ellipticity of the lens object modifies magnification profile. The ellipticity is defined as the ratio of projected semi-major $a$ and semi-minor axis $b$ as $\epsilon=1-b/a$. In this elliptical coordinate, convergence and shear is modified as below.

The elliptical coordinates are given as $x=\sqrt{x^{2}_{1\epsilon}+x^{2}_{2\epsilon}}$, $x_{1\epsilon}=\sqrt{a_{1\epsilon}}x_{1}$, $x_{2\epsilon}=\sqrt{a_{2\epsilon}}x_{2}$ and $\phi_{\epsilon}=\arctan{(x_{2\epsilon}/x_{1\epsilon})}$ which are defined in Golse & Kneib (2002), where $a_{1\epsilon}$ and $a_{2\epsilon}$ are the parameters defining the ellipticity in the lens coordinate. In this study, we applied $a_{1\epsilon}=1-\epsilon$ and $a_{2\epsilon}+\epsilon$. From this, we modify convergence and shear in eq. (18) as $\kappa\to\kappa_{\epsilon}$ and $\gamma\to\gamma_{\epsilon}$. We applied empirical value of ellipticity as $\langle\epsilon_{\rm DM}\rangle=0.482\pm 0.028$ estimated from HST cluster sample which is obtained by Okabe et al. (2020). In this study, we assumed triaxiality of dark matter halos does not depend on their mass for simplicity although they have mentioned that $\epsilon_{\rm DM}$ has lower value for lower $M_{\rm DM}$ based on simulation in Okabe et al. (2018). In this study, we estimated effect of gravitational lensing on the statistical values of galaxy. We adopted the method from Lapi et al. (2012) and Lima et al. (2009). The lensing probability $f(z_{s})$ is computed as,

where $d^{2}V/dzd\Omega$ is the comoving volume element per unit solid angle and $dN/dM_{\rm H}$ is the galaxy halo mass function, $\epsilon$ is the ellipticity of dark matter halos and $P(\epsilon)$ is the probability distribution of the ellipticity. We assumed that $P(\epsilon)$ to be $\delta\mbox{-}$function (i.e., $P(\epsilon)=\delta(\epsilon-\langle\epsilon_{\rm DM}\rangle)$) that only takes value at mean value obtained from observation. $\Delta\Omega$ in the third integration denotes the lensing cross section of a lens halo with $M_{\rm H}$ at $z_{l}$ for a source galaxy at $z_{s}$ whose flux is magnified by larger than $\mu$ times the original. The lensing cross section is derived by integrating the inverse of magnification with image plane coordinate considering the fact that the magnification is the Jacobian of source and image plane,

where $\beta$ and $\theta$ denotes the coordinates of the source plane and image plane respectively. Our logarithmic $f(z_{s})$ is displayed in Fig. 5 as a 2D map on source redshift $z_{s}$ and magnification $\mu$. Lensing probability decrease with $z_{s}$ as the possibility of having dark matter halo between the observer and source always decrease as the light path goes shorter. Especially, $f_{\mu}$ is significantly small The $f_{\mu}$ at $\mu=2$ is approximately $1\times 10^{-2}$ and monotonically decrease with $\mu$ for above $z_{s}=0.2$.

Finally, the lensed differential number counts would be,

where $dn/dS$ is differential source counts at $z=z_{s}$ whose flux is magnified by $\mu$ and $\langle\mu\rangle$ indicates average $\mu$ weighted by lensing probability. Fig. 6 shows the result of analysis of the effect of gravitational lensing on the galaxy differential number counts. This result shows number counts is biased by gravitational lensing approximately $1\%$ at most in radio band with a peak around $S_{\rm 3GHz}\sim 100~{}{\rm\mu Jy}$. For SFGs, the effect would be diminished and becomes negligible towards bright end. This reflects the fact that the number counts of SFGs in bright end is dominated by contribution from nearby universe and lensing probability goes small at low $z_{s}$. We compare our result with lensed SFG differential number counts model obtained by Mancuso et al. (2017) which is constructed based on redshift dependent SFR function. While the peak flux agrees well with their result, our result shows subtler effect on lensed SFG counts especially towards faint population. On the other hand, the behaviour of lensed AGN differential number counts keeps growing trend because the population from high redshift is still considerable as discussed in Sec. 4.2 for AGN. Nevertheless, the effect is smaller than that of SFG.

In this paper, we have estimated the 3GHz LFs up to $z\sim 5.5$ for sources brighter than $S_{\rm 3GHz}>11.0\ {\rm\mu Jy}$ with making use of the data from the VLA-COSMOS 3GHz Large project, one of the deepest radio survey. The LFs are estimated through two independent method, a non-parametric $C^{-}$ method and parametric MCMC method. The sample contains 5410 SFGs and 1908 AGNs classified based on panchromatic counterparts and SFR based radio excess. The evolution of the LF is described by PLE model. The resultant redshift dependency of the characteristic luminosity is, $L_{*,{\rm SFG}}(z)\propto(1+z)^{(3.36\pm 0.05)-(0.31\pm 0.02)z}$ for SFG and $L_{*,{\rm AGN}}(z)\propto(1+z)^{(2.91^{+0.18}_{-0.19})-(0.99\pm 0.08)z}$ for AGN. The 0th and 1st order term of PLE function showed a strong negative correlation. As expected from previous works, SFG shows stronger positive evolution comparing to that of AGN. Our resultant LFs agreed well with that derived by Novak et al. (2018) with MCMC to determine the parameters of total LF obtained via non-parametric $1/V_{\rm max}$ method for overall redshift (Fig. 1). The Euclidean normalized number counts estimated from our PLE model had a good agreement with various observations including ultra-deep surveys below $S_{\rm lim}$ in the VLA-COSMOS 3 GHz large survey while we assumed fixed faint end slope for LF (Fig. 3).

We discussed about cosmic star formation rated density history estimated from radio luminosity. Our SFRD model showed monotonic increase towards $z\sim 3$ and then compared our model with other SFRD models constructed from various wavelength in Fig. 4. We examined SFRD evaluated from constant radio to IR luminosity fraction for all redshift and found that constant $q_{\rm TIR}$ gives higher SFRD at $z>3.5$. Matthews et al. (2021) introduced luminosity dependent non-linear $q$ parameter in order to exclude the bias caused by flux limited sample. The FIR/radio correlation is a critical quantity to draw SFRD from radio luminosity. It remained as an issue to be discussed further.

In the end, we assessed the statistical effect of gravitational lensing due to the dark matter halos with the estimated redshift distribution of galaxies. In the Euclidean normalized number counts, the lensed SFG population shows a peak at $S_{\rm 3GHz}\sim 100\ {\rm\mu Jy}$ with a contribution of about $1\%$ while lensed AGN counts found to be less effective (Fig. 6). Additionally, we performed Fisher analysis to predict the capability of PLE parameters for SFG and AGN parameter constraint with SKA I MID surveys. The direction of main axes of the error elliptical agrees with that is obtained for VLA-COSMOS sample (Fig. 7) and we find that the parameter will be constrain about 10 times tighter in the SKA era (Tab.3).

This work has been supported by the Japan Society for the Promotion of Science (JSPS) Grants-in-Aid for Scientific Research (19H05076 and 21H01128). This work has also been supported in part by the Collaboration Funding of the Institute of Statistical Mathematics “New Development of the Studies on Galaxy Evolution with a Method of Data Science”.

The data appeared in this article will be shared on reasonable request to the corresponding author.