The optically-selected 1.4-GHz quasar luminosity function below 1 mJy

TL;DR

This paper introduces a Bayesian stacking method to measure the radio luminosity function of optically-selected quasars below 1 mJy, revealing the emergence of radio-quiet quasars and potential star formation contributions.

Contribution

The study develops and applies a novel Bayesian stacking technique to accurately reconstruct the quasar radio luminosity function below detection thresholds, extending previous analyses to fainter fluxes.

Findings

RLF agrees with deeper data at low redshift

Radio-loud RLF flattens and steepens at specific luminosities

Radio-quiet quasars emerge at lower luminosities, possibly linked to star formation

Abstract

We present the radio luminosity function (RLF) of optically-selected quasars below 1~mJy, constructed by applying a Bayesian-fitting stacking technique to objects well below the nominal radio flux-density limit. We test the technique using simulated data, confirming that we can reconstruct the RLF over three orders of magnitude below the typical $5\sigma$ detection threshold. We apply our method to 1.4-GHz flux-densities from the Faint Images of the Radio Sky at Twenty-cm survey (FIRST), extracted at the positions of optical quasars from the Sloan Digital Sky Survey (SDSS) over seven redshift bins up to $z=2.15$ {and measure the RLF down to two orders of magnitude below the FIRST detection threshold}. In the lowest redshift bin ($0.2<z<0.45$), we find that our measured RLF agrees well with deeper data from the literature. The RLF for the radio-loud quasars flattens below $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}] \approx 25.5$ and becomes steeper again below $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}] \approx 24.8$, where radio-quiet quasars start to emerge. The radio luminosity where radio-quiet quasars emerge coincides with the luminosity where star-forming galaxies are expected to start to dominate the radio source counts. This implies that there could be a significant contribution from star formation in the host galaxies, but additional data is required to investigate this further. The higher-redshift bins show a similar behaviour as for the lowest-$z$ bin, implying that the same physical process may be responsible.

Tables

Table 1. Table 1: The redshift bins used to separate the sources, along with the absolute-magnitude cut, the number of sources ( N 𝑁 N ), and the number of sources with extracted FIRST flux-densities, N FIRST subscript 𝑁 FIRST N_{\mathrm{FIRST}} (most of these sources are undetected in FIRST). The RQF (‘radio-quiet’ fraction) of the quasars is calculated by integrating over the low-luminosity part of the radio luminosity function described in Sec  3.3 .

Redshift bin	max($M_i$)	$N$	$N_{\mathrm{FIRST}}$	RQF (%)
	-23.0	1234	1222	96.2
	-24.1	1437	1424	91.8
	-24.9	2401	2359	86.3
	-25.4	4534	4472	93.8
	-25.8	5967	5879	93.7
	-26.2	4988	4923	91.4
	-26.6	3250	3211	91.4

Table 2. Table 2: Assumed priors. L 5 ​ σ subscript 𝐿 5 𝜎 L_{5\sigma} is the luminosity corresponding to the 5 ​ σ n 5 subscript 𝜎 𝑛 5\sigma_{n} flux-density cut for a given redshift.

Parameter	Prior
${\alpha }_1,{\beta }_1,{\alpha }_2,{\beta }_2,\delta$	uniform $\in [-5,5]$
${\sigma }_{LF}$	uniform $\in [0.1,2]$
${\log }_{10}\left[L_{{\min }_{\{1,2\}}}/\left({\mathrm{WHz}}^{-1}\right)\right]$	uniform $\in [20,30]$
${\log }_{10}\left[L_{{\max }_{\{1,2\}}}/\left({\mathrm{WHz}}^{-1}\right)\right]$	uniform $\in [20,30]$
${\log }_{10}\left[{\varphi }_{\{1,2\}}^{\ast }/\left({\mathrm{Mpc}}^{-3}{\mathrm{mag}}^{-1}\right)\right]$	uniform $\in [-12,-2]$
${\log }_{10}\left[L_1^{\ast }/\left({\mathrm{WHz}}^{-1}\right)\right]$	uniform $\in [{\log }_{10}\left(L_{5\sigma }\right),30]$
${\log }_{10}\left[L_2^{\ast }/\left({\mathrm{WHz}}^{-1}\right)\right]$	uniform $\in [20,{\log }_{10}\left(L_{5\sigma }\right)]$

Table 3. Table 3: The relative log 10 − limit-from subscript 10 \log_{10}- evidence, Δ ​ log 10 ⁡ Z Δ subscript 10 𝑍 \Delta\log_{10}Z , for the models using single radio luminosity functions (what we call Model A ′ superscript 𝐴 ′ A^{\prime} , B ′ superscript 𝐵 ′ B^{\prime} and C ′ superscript 𝐶 ′ C^{\prime} ) relative to the model with the lowest evidence in each case, applied to the full SKADS sample and the high-luminosity sample with different noise levels.

	15 $\mu$Jy	150 $\mu$Jy
Model	$\Delta {\log }_{10}Z$	$\Delta {\log }_{10}Z$
$A^{\prime }$	$0.0\pm 0.00$	$0.45\pm 0.17$
$B^{\prime }$	$139.1\pm 0.22$	$0.00\pm 0.00$
$C^{\prime }$	$141.8\pm 0.21$	$0.35\pm 0.17$
$A^{\prime }$	$0.00\pm 0.00$	$0.87\pm 0.14$
$B^{\prime }$	$141.5\pm 0.20$	$0.00\pm 0.00$
$C^{\prime }$	$147.2\pm 0.22$	$0.71\pm 0.14$

Table 4. Table 4: Relative evidence of the different models for each redshift bin in the FIRST data. The reference evidence is the model with the lowest value for each redshift bin, and the winning model has the highest relative evidence (in bold).

	$△{\log }_{10}𝒵$	$△{\log }_{10}𝒵$	$△{\log }_{10}𝒵$	$△{\log }_{10}𝒵$	$△{\log }_{10}𝒵$	$△{\log }_{10}𝒵$	$△{\log }_{10}𝒵$
Model
A	$4.4\pm 0.20$	$0.0\pm 0.00$	$2.0\pm 0.21$	$10.1\pm 0.22$	$0.1\pm 0.23$	$0.0\pm 0.00$	$0.0\pm 0.00$
B	$5.6\pm 0.20$	$0.7\pm 0.19$	$2.7\pm 0.21$	$11.9\pm 0.22$	$0.8\pm 0.23$	$3.0\pm 0.22$	$2.0\pm 0.22$
C	$0.0\pm 0.00$	$0.6\pm 0.19$	$0.0\pm 0.00$	$0.0\pm 0.00$	$0.0\pm 0.00$	$1.9\pm 0.22$	$1.9\pm 0.25$

Table 5. Table 5: The MAP posterior parameters for the double power-law – the winning model – for the quasar RLF, in each of the redshift bins and their 2 ​ σ 2 𝜎 2\sigma . The units of the parameters are as shown in Table 2 .

Parameter
${\log }_{10}\left[L_{{\min }_1}\right]$	${20.83}_{-0.79}^{+1.27}$	${20.69}_{-0.65}^{+1.69}$	${21.94}_{-1.87}^{+0.85}$	${20.76}_{-0.71}^{+1.97}$	${20.45}_{-0.40}^{+2.68}$	${22.00}_{-1.90}^{+1.38}$	${22.38}_{-2.35}^{+0.40}$
${\log }_{10}\left[L_{{\max }_1}\right]$	${28.75}_{-3.82}^{+1.14}$	${24.18}_{-0.91}^{+5.65}$	${24.69}_{-0.07}^{+3.25}$	${29.44}_{-3.98}^{+0.43}$	${29.35}_{-3.92}^{+0.52}$	${26.30}_{-0.85}^{+3.57}$	${25.92}_{-0.29}^{+3.94}$
${\log }_{10}\left[L_{{\min }_2}\right]$	${20.60}_{-0.44}^{+2.85}$	${23.43}_{-3.29}^{+0.94}$	${22.69}_{-2.54}^{+1.67}$	${21.26}_{-3.05}^{+3.89}$	${18.63}_{-0.39}^{+6.43}$	${22.31}_{-2.11}^{+2.54}$	${22.71}_{-2.51}^{+2.55}$
${\log }_{10}\left[L_{{\max }_2}\right]$	${29.65}_{-2.18}^{+0.32}$	${28.91}_{-1.36}^{+1.05}$	${29.03}_{-0.91}^{+0.94}$	${29.76}_{-1.49}^{+0.22}$	${29.55}_{-1.49}^{+0.42}$	${28.12}_{-0.09}^{+1.83}$	${29.92}_{-1.26}^{+0.07}$
${\log }_{10}\left[{\Phi }_1^{\ast }\right]$	$-{8.01}_{-3.50}^{+1.57}$	$-{9.02}_{-0.16}^{+0.68}$	$-{8.79}_{-0.35}^{+0.72}$	$-{8.10}_{-0.12}^{+0.26}$	$-{8.03}_{-0.26}^{+0.10}$	$-{8.58}_{-0.18}^{+0.42}$	$-{8.39}_{-0.35}^{+0.16}$
${\log }_{10}\left[L_1^{\ast }\right]$	${25.25}_{-0.55}^{+3.52}$	${27.74}_{-2.72}^{+0.08}$	${27.36}_{-2.39}^{+0.37}$	${25.38}_{-0.13}^{+0.70}$	${25.78}_{-0.08}^{+1.41}$	${27.61}_{-0.82}^{+0.22}$	${25.04}_{-0.01}^{+2.29}$
${\alpha }_1$	${0.82}_{-2.45}^{+3.80}$	${1.16}_{-1.10}^{+2.74}$	${2.27}_{-2.09}^{+1.21}$	${0.46}_{-0.11}^{+0.39}$	${0.50}_{-0.08}^{+0.70}$	${2.46}_{-1.56}^{+1.00}$	${0.27}_{-0.04}^{+0.70}$
${\beta }_1$	$-{3.26}_{-1.60}^{+3.79}$	${0.18}_{-4.70}^{+0.01}$	${0.24}_{-4.90}^{+0.10}$	$-{3.30}_{-1.45}^{+2.84}$	$-{1.31}_{-1.19}^{+1.30}$	${0.11}_{-0.56}^{+0.11}$	$-{3.44}_{-1.29}^{+3.53}$
${\log }_{10}\left[{\Phi }_2^{\ast }\right]$	$-{6.73}_{-4.91}^{+0.04}$	$-{6.94}_{-0.16}^{+0.00}$	$-{6.97}_{-0.12}^{+0.15}$	$-{6.68}_{-0.02}^{+0.11}$	$-{6.53}_{-0.09}^{+0.07}$	$-{6.55}_{-0.12}^{+0.10}$	$-{6.73}_{-0.19}^{+0.09}$
${\log }_{10}\left[L_2^{\ast }\right]$	${22.55}_{-1.36}^{+0.78}$	${22.74}_{-0.14}^{+0.20}$	${23.25}_{-0.42}^{+0.22}$	${23.54}_{-0.27}^{+0.02}$	${23.53}_{-0.11}^{+0.26}$	${23.76}_{-0.20}^{+0.23}$	${24.00}_{-0.32}^{+0.07}$
${\alpha }_2$	${0.68}_{-0.24}^{+3.37}$	${1.05}_{-0.40}^{+0.06}$	${1.21}_{-0.65}^{+0.54}$	${1.07}_{-0.22}^{+0.14}$	${1.05}_{-0.17}^{+0.47}$	${1.31}_{-0.40}^{+0.47}$	${1.67}_{-0.80}^{+0.26}$
${\beta }_2$	$-{1.92}_{-2.97}^{+1.92}$	$-{4.21}_{-0.69}^{+2.67}$	$-{1.82}_{-3.09}^{+0.42}$	$-{1.48}_{-3.31}^{+0.06}$	$-{3.56}_{-1.33}^{+2.22}$	$-{4.05}_{-0.85}^{+2.67}$	$-{2.75}_{-2.17}^{+0.96}$

Equations

$$M_i=m_i-5{\log }_{10}[{\text{d}}_L(z_{\mathrm{up}})/10]-K(z),$$

$$S=\min (1.40S_{\mathrm{F}},S_{\mathrm{F}}+0.25\mathrm{mJy}),$$

$$S_m=S_{\mathrm{F}}+n,$$

$$\text{centering}\mathcal{P}(\Theta |D,H)=\frac{\mathcal{L}(D|\Theta ,H)\Pi (\Theta |H)}{\mathcal{Z}},\text{@add@centering}$$

$$\mathcal{Z}=\int \mathcal{L}\Pi {\text{d}}^n\Theta .$$

$${\mathcal{L}}_i\left(k_i|\mathbf{\Theta }\right)=\frac{I_i^{k_i}{\text{e}}^{-I_i}}{k_i!},$$

$$I_i={\int }_{S_{min}}^{S_{max}}dS\frac{dN(S)}{dS}{\int }_{S_{m_i}}^{S_{m_i}+\Delta S_{m_i}}d{S}_m\frac{1}{{\sigma }_n\sqrt{2\pi }}{\mathrm{e}}^{-\frac{{\left({\mathrm{S}}_{\mathrm{F}}-{\mathrm{S}}_{\mathrm{m}}\right)}^2}{2{\sigma }_{\mathrm{n}}^2}}.$$

$$I_i=$$

$$\frac{1}{2}\left\{\mathrm{erf}\left(\frac{S_{\mathrm{F}}-S_{m_i}}{{\sigma }_n\sqrt{2}}\right)-\mathrm{erf}\left(\frac{S_{\mathrm{F}}-(S_{m_i}+\Delta S_{m_i})}{{\sigma }_n\sqrt{2}}\right)\right\}.$$

$$\mathcal{L}\left(\mathbf{k}|\mathbf{\Theta }\right)=\prod _{i=1}^N{\mathcal{L}}_i\left(k_i|\mathbf{\theta }\right).$$

$$L_{\nu }=4\pi D_L^2(1+z{)}^{\alpha -1}{S}_{\nu },$$

$$\Phi (L_{\nu })=\frac{dN}{dVdm}=\frac{dN}{dVd{L}_{\nu }}\frac{d{L}_{\nu }}{dm}=\ln (10^{0.4})L_{\nu }\rho (L_{\nu }).$$

$$\Phi (L{)}_A=\frac{{\Phi }_1^{\ast }}{(L/L_1^{\ast }{)}^{{\alpha }_1}+(L/L_1^{\ast }{)}^{{\beta }_1}}+\frac{{\Phi }_2^{\ast }}{(L/L_2^{\ast }{)}^{{\alpha }_2}}.$$

$$\Phi (L{)}_B=\frac{{\Phi }_1^{\ast }}{(L/L_1^{\ast }{)}^{{\alpha }_1}+(L/L_1^{\ast }{)}^{{\beta }_1}}+\frac{{\Phi }_2^{\ast }}{(L/L_2^{\ast }{)}^{{\alpha }_2}+(L/L_2^{\ast }{)}^{{\beta }_2}}.$$

$$\Phi (L{)}_C=$$

$$+{\Phi }_2^{\ast }{\left(\frac{L}{L_2^{\ast }}\right)}^{1-\delta }\exp \left[-\frac{1}{2{\sigma }_{LF}^2}{\log }_{10}^2\left(1+\frac{L}{L_2^{\ast }}\right)\right].$$

$$\frac{dN}{dS}$$

$$=\rho (L)V_{i}4\pi D_L^2(1+z_i{)}^{\alpha -1}$$

$$=\frac{\Phi (L)V_i}{L\ln (10^{0.4})}4\pi D_L^2(1+z_i{)}^{\alpha -1},$$

$$\Phi (L_{\nu })=\frac{1}{\Delta m}\sum _{i=1}^N{\left(\frac{1}{V_{\max }}\right)}_i,$$

$$\sigma (\Phi )=\frac{1}{\Delta m}{\left[\sum _{i=1}^N{\left(\frac{1}{V_{\max }}\right)}_i^2\right]}^{1/2},$$

Full text

††pagerange: The optically-selected 1.4-GHz quasar luminosity function below 1 mJy–The optically-selected 1.4-GHz quasar luminosity function below 1 mJy††pubyear: 2019

1 Introduction

The evolution of quasars has been a subject of interest right since their discovery (Schmidt, 1963). Quasars have been of particular interest over the past decade due to the role that they — and active galactic nuclei (AGN) in general — play in galaxy evolution. For example, feedback from AGN may expel or heat gas in a galaxy, thereby quenching star formation (SF) in the host galaxies (e.g. Granato et al., 2004; Scannapieco & Oh, 2004; Croton et al., 2006; Hopkins et al., 2008; Antonuccio-Delogu & Silk, 2008), or feasibly in the wider environment (e.g. Rawlings & Jarvis, 2004; Hatch et al., 2014). This may be a major contributor to establishing the observed relationship between supermassive black holes (SMBHs) and the central bulge properties in a galaxy (e.g. Ferrarese & Merritt 2000; Hopkins et al. 2006).

They were originally discovered as strong radio sources and later also found to be bright in the optical (e.g. Schmidt, 1963). However only $\sim 10$ per cent of optically-selected quasars were detected in large-area radio surveys (e.g. Strittmatter et al. 1980). The sources that were detected in these surveys were termed ‘radio-loud quasars’, while the remaining 90 per cent of the quasar population, which are fainter in the radio, were referred to as ‘radio-quiet quasars’. The radio emission from radio-loud quasars is known to be mainly dominated by synchrotron radiation from electrons accelerated by powerful jets, while the source of radio-quiet quasars is still debated. One suggestion is that the radio emission from radio-quiet quasars is a result of synchrotron radiation from supernova explosions associated with star formation in the host galaxy, rather than being the result of AGN processes (e.g. Terlevich et al. 1987, 1992; Padovani et al. 2011; Kimball et al. 2011; Bonzini et al. 2013; Condon et al. 2013; Kellermann et al. 2016; Gürkan et al. 2018; Stacey et al. 2018). However, some authors suggest the radio emission in radio-quiet quasars is still dominated by AGN-related processes such as low-power jets (e.g. Falcke & Biermann, 1995; Wilson & Colbert, 1995; Hartley et al., 2019), accretion disk winds (e.g. Jiang et al., 2010; Zakamska & Greene, 2014), coronal disk emissions (Laor & Behar, 2008; Laor et al., 2019) or a combination of these process (see Panessa et al. 2019 for a review), with factors such as different accretion rates (Fernandes et al., 2011), SMBH spin (Blandford & Znajek, 1977; Schulze et al., 2017), SMBH mass (Dunlop et al., 2003; McLure & Jarvis, 2004), host-galaxy morphology (Bessiere et al., 2012), galactic environments (Fan et al., 2001), or a combination of these, being responsible for lack of powerful jets .

One of the ways to study quasars and their source of radio emission is through luminosity functions (LFs, i.e. the number of sources with a certain luminosity in a given volume and luminosity bin). It is now accepted that SMBHs accrete most of their mass during the active-galaxy phase, when they are radiating at quasar luminosities (Salpeter, 1964; Zel’dovich & Novikov, 1965; Lynden-Bell, 1969; Soltan, 1982). Therefore, with accurate measurements of the quasar LF and its evolution, one can map out the SMBH accretion history (e.g.  Shankar et al. 2009; Shankar 2010; Shen 2009; Shen & Kelly 2012), constrain the formation history of SMBHs (e.g. Rees 1984, Haiman et al. 2012), and potentially determine the contribution of quasars to feedback.

The radio luminosity function (RLF) of radio-loud quasars is well-studied (e.g. Schmidt, 1970; Willott et al., 1998; Jiang et al., 2007), but the faint (radio-quiet) end is not well-explored, as these fainter sources lie below the detection threshold of most wide-area radio surveys. There are various methods used in the literature to study radio-quiet populations. One such method is through deep-narrow radio surveys (e.g. Condon et al., 2003; Kellermann et al., 2008; Padovani et al., 2009; Padovani et al., 2011; Miller et al., 2013). Such surveys have contributed to our understanding of the radio emission from the radio-quiet population. For instance, Padovani et al. (2015) found that emission from radio-quiet AGNs has a contribution from black-hole activity as well as emission related to star formation. However, very few genuinely luminous quasars are detected in these deep-narrow survey ($\sim$ 15 quasars per deg2).

The most popular means of studying $\mu$Jy sources in the past two decades have involved some form of ‘stacking’ (Ivezić et al., 2002; White et al., 2007; Hodge et al., 2008; Mitchell-Wynne et al., 2014; Roseboom & Best, 2014; Zwart et al., 2015a). There are a number of different versions and definitions of stacking seen in the literature (see Zwart et al. 2015a for an overview). Usually stacking involves using positional information of a source population that is selected (and classified) from an auxiliary survey, and then extracting the flux density at those positions in the survey of interest (where they are above or below the detection threshold). In most cases stacking is used to explore the average (mean, median or weighted versions thereof) properties of sources below the detection threshold (i.e. $<<5\sigma$). For example, stacking can be employed to infer average SF rates (e.g. Dunne et al. 2009; Karim et al. 2011, Zwart et al. 2014), where 1.4-GHz radio flux-densities are extracted at positions of sources selected by stellar mass.

Traditional stacking techniques have added a great deal to our understanding of $\mu$Jy source populations. However, they only return a single statistic, and new techniques have been developed that extract more information from the stacked data. Mitchell-Wynne et al. (2014) went beyond stacking by combining stacking with maximum-likelihood methods to fit a source-count model to the stacked sources. Roseboom & Best (2014) adopted a similar approach to Mitchell-Wynne et al. (2014) by fitting a luminosity-function model to stacked star-forming galaxies. Zwart et al. (2015b) then extended the technique of Mitchell-Wynne et al. (2014) to a fully-Bayesian framework (bayestack), which allows for model selection. Chen et al. (2017) extended the technique of Zwart et al. (2015b) by including the effects of the point spread function and source confusion, an approach that incorporates some of the reasoning from Vernstrom et al. (2014), which combined a traditional $P(D)$ analysis with a Bayesian likelihood model fitting.

In this work we measure the RLF of optically-selected quasars below 1 mJy by building on the work of Roseboom & Best (2014) and Zwart et al. (2015b). We use a set of models for the RLF and fit directly to the radio data using a full Bayesian approach. We apply the technique to a large sample of quasars from the Sloan Digital Sky Survey (SDSS; York et al. 2000) Data release 7 (DR7) quasar catalogue (Shen et al., 2011), using flux densities taken from the Faint Images of The Radio Sky at Twenty-centimeters (FIRST; Becker et al. 1995).

In Section 2 we describe the optical and radio data used in this study. We then outline our technique for making measurements below the noise level using bayestack (Section 3). In Section 4 we test the technique, and our results are given in Section 5. We discuss the results and compare them to the literature in Section 6, finally concluding in Section 7. Throughout the work, unless stated otherwise, we use AB magnitudes and the positions are in the J2000 equinox. We set the spectral index, defined as $\alpha\equiv\log(S/S_{0})/\log(\nu_{0}/\nu)$, to $\alpha=0.7$ (e.g. Kukula et al., 1998) when converting flux density to luminosity and one reference frequency to another. We assume a $\Lambda$CDM cosmology, with $H_{0}=70$ km*-1* Mpc*-1*, $\Omega_{\Lambda}=0.7$ and $\Omega_{m}=0.3$.

2 Data

In a stacking experiment, where we try to extract information from undetected sources in a given survey, one needs data from another survey in which the sources have already been identified. In this paper we use optically-selected quasars from the SDSS and radio data from the FIRST survey.

2.1 The optical quasar sample

The optical data are drawn from the quasar catalogue (Schneider et al., 2010) of the SDSS seventh data release (DR7, Abazajian et al. (2009)). In SDSS, quasars are mainly identified using colour selection for objects in the magnitude range $15<i<19.1$ (Richards et al., 2002; Richards, 2006). Quasars are then differentiated from galaxies and stars by their unique colours in multi-dimensional colour-colour space (Fan, 1999): SDSS’s candidate quasars are primarily outliers from stellar regions in colour-colour space (Richards et al., 2001), and the regions having large stellar contamination were avoided. The quasar sample includes additional sources that are selected because they have a FIRST counterpart. A source is targeted for spectroscopic follow-up if it is within 2 arcsec from a source in the FIRST catalogue. The final catalogue contains 105,783 spectroscopically-confirmed quasars, all brighter than $M_{i}=-22$ with at least one emission line with full width at half-maximum greater than 1,000 km/s or a relevant absorption feature.

We use a subsample consisting of 59,932 quasars selected across the survey area (purely colour-selected sources with the flag UNIFORM=1 from Shen et al. 2011) for the purpose of having a homogeneous sample of quasars. This sample covers an effective area of 6,248 deg2 (Shen et al., 2011). Fig. 1 shows the distribution of sources in absolute magnitude and redshift.

We divide the sources into seven redshift bins (see Table 1) reducing the total to 48,046 sources. Since each redshift bin has a non-negligible width, we apply an absolute-magnitude cut to each redshift bin (corresponding to a minimum luminosity cut per bin) to ensure all quasars in the bin are observed within the sensitivity limit of the survey. This reduces the total number of sources to 24,003. The maximum absolute magnitude in each redshift bin corresponds to the optical flux limit at the highest redshift in that bin given by

[TABLE]

where $m_{i}=18.7$ is just above the magnitude completeness limit ($m_{i}=19.1$) for DR7, $\rm{d}_{L}(z_{\rm{up}})$ is the luminosity distance (in pc) at the upper redshift of the bin and $K(z)$ is the $K$-correction from Richards (2006).

2.2 Radio data

The FIRST survey (Becker et al. 1995) was carried out with the Very Large Array (VLA; Thompson et al. 1980) in its ‘B’ configuration at 20 cm (1.4 GHz), yielding a synthesized beam size of $5.4^{\prime\prime}$ (FWHM). It covered 8,444 deg2 in the North Galactic cap and 2,131 deg2 in the South Galactic cap giving a total coverage of 10,575 deg2. The survey footprint overlaps with the area that SDSS covered in the North Galactic Cap, as well as with a smaller $\simeq$ 2.5 deg2 wide strip along the Celestial equator. The maps have a rms of $\approx 150\,\mu$Jy/beam. The survey catalogue contains more than 800,000 sources above the detection limit of 1 mJy, and includes peak and integrated flux-densities calculated by fitting a two-dimensional Gaussian to each source. The survey is 95 per cent complete at 2 mJy and 80 per cent complete at 1 mJy. The maps are stored as FITS images and have 1.8*′′* pixels.

2.3 Cross-matching catalogues

We first matched the SDSS quasars with detected sources from the FIRST catalogue. The allowed separation between the coordinates of the two catalogues should be as small as possible to avoid random matching with other sources, but also large enough to ensure real matches are not omitted because of slight random offsets in position between the optical and radio data.

Fig. 2 shows the results of matching our sample to the FIRST catalogue. We choose a limiting separation of 1.8 arcsec based on Fig. 2, which is the pixel size of the FIRST images. From the original 105,783 quasars we made 3,815 matches ($\sim$3 per cent), which is consistent with the low number of optical-to-radio matches found by Paris et al. (2012) and Pâris et al. (2017). We find 2,381 ($\sim$ 10 per cent) matches from our sample of 24,003 SDSS quasars.

The FIRST catalogue only contains sources with flux densities above the detection threshold of $\sim 1$ mJy. In order to obtain sources with flux densities below the FIRST detection threshold we extracted $11\times 11$ pixel stamps ($19.8\times 19.8$ arcsec) from the FIRST maps, centered on the SDSS quasar positions, and used the central pixel value as the radio flux-density of the quasar. 23,490 of our quasars have fluxes densities, the rest fall outside FIRST coverage.

In Fig 3 we compare the catalogued peak flux-densities and the extracted flux-densities for 2,381 detected sources. Most of the extracted flux-densities are in good agreement with peak flux-densities, with the exception of fluxes densities below 10 mJy, which underestimate the peak flux densities. There are also about 10 sources with high scatter from the peak flux-densities. The difference between the extracted flux-densities and peak flux-densities at low flux densities will affect the results and therefore needs to be accounted for and understood. Note, however, that there could be a difference in the effect on extraction of high signal-to-noise (detected) sources compared to the undetected ones. For instance, detected sources could be more extended and therefore slightly resolved by the FIRST restoring beam. Other possible contributions to the difference in the flux densities are clean and snapshot biases.

clean bias is a systematic effect that decreases the peak flux-density of a source above the detection limit and redistributes it around the map. This phenomenon is associated with the non-linear clean process (Condon et al. 1994) and affects large-area radio surveys such as FIRST and the NRAO VLA Sky Survey (NVSS; Condon et al. 1998). The bias is additive and has an approximately constant magnitude, with a value of 0.25 mJy beam*-1* for FIRST (Becker et al., 1995). White et al. (2007) discovered another bias that affects sub-threshold sources (which are not cleaned) and suggested that it is associated with the sidelobes of the beam pattern. This snapshot bias behaves differently from the one associated with clean as it is multiplicative (i.e. the higher the flux density the higher the bias). The proposed total bias correction summarized by White et al. (2007) is

[TABLE]

where $S$ should be the intrinsic flux-density of the source and $S_{F}$ is the noiseless flux-density from FIRST. This is an idealised case, where $S_{\rm F}$ is meant to incorporate the calibration effects in the FIRST data. It is important to note that this correction can only be applied when the noise can be neglected, that is, for the case of detected sources or the stacked median flux as done in White et al. (2007). With the correction, the low (detected) flux-densities are in good agreement with the catalogue flux-densities (Fig 3). Since, in our case, we are also dealing with noise dominated sources, these corrections need to be incorporated directly into the likelihood function in a forward model way, as described in Section 3.2. Throughout the paper, we use $S_{m}$ to indicate the measured fluxes from FIRST, which will include noise, as opposed to the ideal noiseless case above, that is:

[TABLE]

where $n$ represents the noise distribution.

To proceed with the analysis, the sources in each redshift bin in Table 1, were further binned in terms of the measured radio flux-density ($S_{m}$, Fig. 4). This includes both detected and undetected sources. One can see from the negative side of the flux-density distributions that the noise is Gaussian to a good approximation, whilst there is a tail on the positive side of the distributions, which shows the contribution from faint real sources. There is an offset in the noisy flux-densities because the average ‘true’ flux-density of these faint sources is comparable to the noise. The flux-density distribution is more Gaussian-like if the noise is much larger than the true flux-density of the faint sources. This effect is more-pronounced for the higher-redshift bins, where a greater fraction of sources are undetected.

3 Bayestack framework

Our stacking analysis is based on a Bayesian formalism that can probe the quasar RLF below the FIRST detection threshold, down to sub-mJy levels. We made use of a modified version of the software bayestack (Zwart et al., 2015b). The idea is to start with a model for the RLF for a given redshift bin. We then translate that into a source-count model, take into account the FIRST bias correction, and fit to the number of sources per flux-density bin, as extracted from the data. Below we review the basics of the method.

3.1 Bayesian analysis

The fitting approach centres on Bayes’ theorem,

[TABLE]

where $\mathcal{P}$ is the posterior distribution of the parameters $\Theta$, given the data $D$ and model $H$. $\mathcal{L}$ is the likelihood, the probability distribution of the data given the model and parameters, and $\Pi$ is the prior, the known constraints on the parameters. $\mathcal{Z}$ is the Bayesian evidence, which normalizes $\mathcal{P}$ and can be written as an integral of $\mathcal{L}$ and $\Pi$ over the $n$-dimensional parameter space $\Theta$,

[TABLE]

A model has high evidence when a large portion of its prior parameter space is likely (i.e. large likelihood), and small evidence when a large portion of its parameter space has a small likelihood, irrespective of how peaked the likelihood function is. This therefore automatically encapsulates Occam’s razor (e.g. Feroz et al. 2009b).

In order to compute this posterior distribution, one needs to sample from it. Sampling has always been one of the most computationally expensive parts of model selection because it involves solving the multidimensional integral in Eq. 5. Nested sampling (Skilling 2004) was created for its efficiency in calculating the evidence, with an added bonus of producing posterior inferences as a by-product. MultiNest  (Feroz et al., 2009b; Feroz et al., 2009a; Buchner et al., 2014) is a robust implementation of nested sampling, returning the full posterior distribution from which the uncertainty analysis can be correctly undertaken.

In Bayesian model selection, one compares the evidences of two models, A and B. This is quantified by considering the ratio of their evidences $\mathcal{Z}_{A}/\mathcal{Z}_{B}$ (equivalent to the difference of their log-evidence, $\ln[\mathcal{Z}_{B}-\mathcal{Z}_{A}$]), known as the Bayes factor. Jeffreys (1961) introduced a way to conclude how much better Model A is compared to B using the Bayes factor: $\Delta\ln\mathcal{Z}<1$ is ‘not significant’, $1<\Delta\ln\mathcal{Z}<2.5$ is ‘significant’, $2.5<\Delta\ln\mathcal{Z}<5$ is ‘strong’, and $\Delta\ln\mathcal{Z}>5$ is ‘decisive’. We adopt this scale in our analysis and use it to compare different luminosity models.

3.2 Assumed Likelihood

To proceed with our Bayesian analysis we need a likelihood for the data, which in our case comprises the extracted flux-densities, $S_{m}$. As explained in Section 2.3, this flux density is a combination of the biased FIRST flux density ($S_{\rm F}$) of the source and the noise distribution (Eq. 3). The noise is assumed to follow a Gaussian distribution, centered at zero with a variance $\sigma_{n}^{2}$. This is a good assumption considering that the flux-density distribution of source-less (‘random sky’) FIRST extracted flux-densities is well approximated by a Gaussian.

Since we are working with binned flux-densities, the likelihood of finding $k_{i}$ objects in the $i^{th}$ flux-density bin [$S_{mi},S_{mi}+\triangle S_{m}$] follows a Poisson distribution,

[TABLE]

where $I_{i}$ is the theoretically-expected number of sources in the $i^{th}$ measured bin, $[S_{m_{i}},S_{m_{i}}+\Delta S_{m}]$, given by the modified equation taken from Mitchell-Wynne et al. (2014),

[TABLE]

Here $dN/dS$ is the source-count model (number of sources per flux-density bin), $\sigma_{n}$ is the mean noise of the data and $S$ is again the intrinsic flux density of the source. The noiseless FIRST flux density $S_{\rm F}$ is related to $S$ by Eq. 2. Therefore, in order to compare to the measured flux, we need to first apply the expected bias effects of the FIRST observation, e.g. $S_{\rm F}=\max\{S/1.4,(S-0.25\rm{\,mJy})\}$. This approach naturally takes into account sample variance (at the Poisson level) since it does not fix the total number of predicted sources to the observed number (e.g. other regions of the sky could have a different total number). This will have implications for the allowed minimum and maximum flux-density values of our fits, as we will see later. The fitting will have large variance at the low flux-density level (because of the noise) and at the high flux-density level (because of Poisson fluctuations due to the low number of sources). Solving the second integral, Eq. 7 becomes

[TABLE]

The total likelihood for the $N$ bins is given by the product of the likelihood in each bin, assuming that the bins are independent,

[TABLE]

As we aim to fit models that describe the radio luminosity function, we need to convert those luminosity-function models to source counts, $dN/dS$, and compare to the binned flux-densities in data space where the noise is Gaussian. As a final detail, we would like to point out that we include bins with zero sources at the low flux-density (negative) end. This means we do not actually see any galaxies below a certain flux-density level (including noise) and models that predict galaxies in those flux-density bins should be penalized. At the high flux-density level this is not done as the maximum flux-density cutoff is our choice, and models that predict some sources above that should not be penalized. However, such models will likely over-predict sources in our highest flux-density bin and will therefore have a lower probability. In any case, such a choice has very little impact on the low-flux-density stacked sample that we are targeting in our analysis.

3.3 Models for the radio luminosity functions

The luminosity per unit frequency (luminosity density) of a radio source, $L_{\nu}$, can be related to the observed flux-density at the same frequency, $S_{\nu}$, through

[TABLE]

where $D_{L}$ is the luminosity distance, $\alpha$ is the spectral index of the source, and $z$ is the redshift of the source.

The luminosity function (LF), $\rho(L_{\nu})$, is the number density of sources per luminosity density bin, e.g. $\rho(L_{\nu})=dN/(dLdV)$ (where $dV$ is comoving volume). Another common definition of the LF ($\Phi$), which we use here, involves binning the source counts in magnitude, that is, $m-m_{0}=-2.5\log_{10}(L/L_{0})$. The relationship between these two definitions is then

[TABLE]

We define parametric models for the quasar RLF consisting of two functions, one for the luminous sources and the other for faint sources (using subscripts 1 and 2 respectively). The radio-loud quasar RLF has been shown to follow a double power-law (see e.g. Boyle et al. 1988), so we parameterize the luminous part of the RLF as a double power-law for all the models considered here. The shape of the quasar RLF at low luminosities is still uncertain, so for that we consider 3 models: a power-law, a double-power-law and a log-normal power-law.

Model A is the simplest overall form for the quasar RLF – a double power-law for the high luminosities (the detected sources) and a single power-law to describe the RLF at low luminosities:

[TABLE]

Note that $L_{2}^{*}$ and $\Phi_{2}^{*}$ will be degenerate here, but we keep this form for convenience.

Model B has a double power-law for both the high- and low-luminosity sources:

[TABLE]

Model C has a double power-law for the luminous sources and a log-normal power-law, which has earlier been used for star-forming galaxies (Tammann et al. 1979), for low-luminosity sources:

[TABLE]

Finally, we note that each of the model functions will be bounded: $L_{{\rm min}_{1}}\leq L\leq L_{{\rm max}_{1}}$ for the high-luminosity end and $L_{{\rm min}_{2}}\leq L\leq L_{{\rm max}_{2}}$ for the low-luminosity end. The boundaries are allowed to overlap since there might be a contribution from both populations. We also consider a different set of models when fitting to simulations which include only the part used for the low luminosity region (parameters with subscript ”2”). We call it Model $A^{\prime}$ (single power law), Model $B^{\prime}$ (one double power-law) and Model $C^{\prime}$ (a log-normal power-law).

The likelihood (Eq. 8) is computed in flux-density space, which means that our LF models, $\Phi(L)$, have to be converted into source-count models, $dN/dS$:

[TABLE]

where $V_{i}$ is the volume of the survey for the redshift bin $i$ and $z_{i}$ is the mean redshift for that bin.

3.4 Priors

Priors play an important role in Bayesian inference as they define the sampled parameter space. A uniform prior is the simplest form, providing an equal weighting of the parameter space. We assign a uniform prior to the slopes $\alpha_{1,2}$, $\beta_{1,2}$, and $\delta$. $\sigma_{LF}$ also has a uniform prior. To avoid degeneracy in the slopes for the double power law, we also impose $\alpha_{1,2}\geq\beta_{1,2}$. $L_{1,2}^{*}$, $L_{{\rm min}_{1,2}}$, $L_{{\rm max}_{1,2}}$ and $\phi_{1,2}^{*}$ all have uniform priors in log-space. The priors are summarised in Table  2.

Combining Eq. 9 with the priors shown in Table 2, and substituting into Eq. 4, one can determine the posterior probability distribution as well as the evidence. We use a Python implementation (Buchner et al. 2014) of MultiNest  (PyMultiNest) to fit the models with evidence_tolerence=0.5 and sampling_efficiency=0.1.

4 Tests on simulated data

We first test our technique by applying it to the Square Kilometre Array Design Studies SKA Simulated Skies (SKADS-S3) simulations (see Wilman et al. 2008; Wilman et al. 2010). SKADS is a semi-empirical simulation of the extragalactic radio continuum sky, covering a sky area of $20\times 20$ deg2 with $\approx 320$ million sources out to a redshift of $z=20$ and flux density of 10 nJy.

We took $\sim 555,000$ sources contained within a 8 deg2 patch of the simulation in the redshift range $1.0<z<1.3$. 223,457 of those sources have radio luminosities between $20.5<\log_{10}[L/{\rm W\,Hz}^{-1}]<24.5$ and we call this the full sample. In order to test how a higher luminosity cut may alter our fits, we also consider a brighter sample of 91,458 which lie between $21.5<\log_{10}[L/{\rm W\,Hz}^{-1}]<24.5$. Such luminosity cut could arise due to the input optical sample being being flux limited, and if there is a correlation between the optical emission and radio emission, this would in turn lead to a downturn in the measured RLF that may not happen if one could measure it directly from a purely radio-selected sample.

We added random noise – generated from a Gaussian distribution, with standard deviation $\sigma_{n}$ that corresponds to the FIRST rms of 150 $\mu$Jy ($S_{n}=S+N[150$ $\mu$Jy, 0]) – to both the high luminosity sample ($\log_{10}[L/(\rm{WHz}^{-1})]>21.4$) and the full sample ($\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]>20.4$), to emulate the observed data (i.e. the ‘noisy’ sources in FIRST). We bin the noisy SKADS sources in flux density and apply our technique for fitting the three models (setting $S=S_{\rm F}$ as there are no calibration biases in the simulated data). In this case we only fit a single function from each model (either a power-law, double power-law or log-normal) to the faint SKADS sources, to test the technique on sources around and below the detection threshold (what we call model $A^{\prime}$, $B^{\prime}$ and $C^{\prime}$). We repeat this using the extreme case of a noise level of 15 $\mu$Jy which will allow for detected sources.

We note that with real data, if the parent catalogue (in the case of the SDSS quasar sample we use in this paper) is flux limited, this may naturally lead to a lower limit in radio luminosity that we can probe, if there is a correlation between optical and radio luminosity (e.g. Serjeant et al., 1998; White et al., 2017).

MultiNest returns the Bayesian evidence of the model and the posterior distribution for all the fitted parameters. The ‘relative evidence’ for a model is the difference between the model evidence and the reference-model evidence (where the reference model is the model with the lowest evidence). We show the relative evidence for the SKADS samples in Table 3, where the winning model, the one with the highest relative evidence, is in bold. From the relative evidences it is clear that the data prefers the log-normal function (Model $C^{\prime}$) for the 15 $\mu$Jy noise levels in both samples and the power-law (Model $A^{\prime}$) for the 150$~{}\mu$Jy noise levels, although the evidence is marginal between models for the 150$~{}\mu$Jy noise level. The evidence also suggests that the power-law function (Model $A^{\prime}$) is a significantly poor fit compared to the other models for the 15 $\mu$Jy noise-levels.

In Fig. 5 we show the one-dimensional (1-D) and two-dimensional (2-D) posterior distributions for the fits of the various models to the ‘noisy’ low-luminosity SKADS sample. The 1-D posterior distribution is the marginalization for each parameter, located at the end of each row in Fig. 5. The peaks in 1D do not always do justice to the 2-D posteriors as they are not just simple Gaussians. They show distorted ‘banana-like’ shapes, with some having long tails. The limits on each plot are the maximum and minimum values from the posterior distribution. Some of these parameters are unconstrained and therefore limited by the assumed priors. Note, however, that this has no effect on the reconstructed RLFs.

Along with the posterior distribution, MultiNest returns three values to summarize each parameter: the mean, maximum-likelihood and maximum-a-posteriori (MAP, maximizing the product of the likelihood and prior) values. Obtaining a single value for a parameter is straightforward if the 1-D posterior is Gaussian, as the mean, MAP and maximum likelihood are the same or very close to each other. It is clear that some of the posteriors in Fig. 5 are not Gaussian, which would mean that the three summaries are likely to be different from each other. We use all three parameters to reconstruct the LFs in turn. (The maximum likelihood gives the same value as the MAP for models with uniform priors, so we just quote the MAP.) Although they are good estimates, they still do not fully describe the complex nature of the posterior, as clearly shown in Fig. 5.

In Fig. 6 we show the reconstruction of radio luminosity functions (RLFs) of the noisy SKADS sources, using both the $21.5<\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]<24.5$ and $20.5<\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]<24.5$ samples. We also show the average total RLF, MAP functions and the 95 per cent confidence interval for each model fit and noise level. Such a choice is not unique as the models that span the 95 per cent confidence interval do not necessarily give a continuous region in terms of the RLF curves. For plotting such a region, we chose a set of luminosity bins and calculated all the values of the RLFs in each bin corresponding to all the models in the posterior to determine the 95 per cent confidence limits.

Since this is a simulation, we can calculate the true underlying RLF by converting flux density to luminosity and bin in luminosity and volume (given the sky area and redshift bin). We therefore show the comparison to these RLFs in Fig. 6. We see that the reconstructed RLFs with 150-$\mu$Jy noise levels have a large scatter but are in good agreement with the true SKADS RLF. As expected, using lower noise levels produces RLF reconstructions with better fits to the SKADS RLF and smaller 95-per-cent confidence regions (the 15 $\mu$Jy noise level panels in Fig. 6). Thus the fitting method works for our current noise levels and for those that will be obtained by future radio surveys. We see that the fit is unbiased, though (of course) if the model is quite poor, the fitting will also be poor (as in the case of the power law). Moreover, the fitting is not affected when we use the sample that includes sources down to $\log_{10}[L/(W\,Hz^{-1})]=20.5$, although (as one would expect) the uncertainty increases as we move to lower luminosities.

5 Results

Having now illustrated the effectiveness of the bayestack algorithm, we apply it to the observational data from SDSS and FIRST. We apply the technique for all three models to each of the redshift bins shown in Table 1.

For the high signal-to-noise (detected) flux-densities we can calculate the luminosity function directly by converting flux density to luminosity (neglecting noise) and binning the number of sources in luminosity. Although the source populations are volume-limited in the optical (i.e no brightness cutoff), in the radio some of the sources might not be above the radio flux-density threshold if they are placed at the highest redshift for a given bin. Therefore, we need to apply the $1/V_{\rm{max}}$ correction: the spectral RLF of $N$ sources in a logarithmic bin of width $\Delta m$, using the $1/V_{\rm{max}}$ method (Schmidt, 1968) is given by,

[TABLE]

with an uncertainty

[TABLE]

where $V_{\rm{max}}$ is the maximum comoving volume at which the source is detected.

5.1 Quasars at $0.2<z<0.45$

We start with the lowest redshift sample because it allows a direct comparison to the work of Kellermann et al. (2016).

We use the Bayesian technique to fit the three RLF models to all the sources (radio detected and undetected) from a volume-limited sample defined by $M_{i}<-23$ at $0.2<z<0.45$. The first column of Table 4 shows the relative evidence of the fit to the models. From the relative evidence we conclude that the data significantly prefer Model B, which consists of a double power-law for the luminous sources and a second double power-law for the low-luminosity and undetected sources.

Fig. 7 shows the posterior distributions for the winning Model B. Like with the SKADS sources, the boundary parameters $L_{\rm{max}_{1,2}}$ along with $L_{\rm{min}_{1,2}}$ exceed the prior limit and are unconstrained. Note, however, that this has very little impact on the actual observed numbers (and that the uncertainty increases due to noise and/or Poisson fluctuations). The other parameters have well-defined peaks, except for the faint-end slopes $\beta_{1,2}$ which span a large range below 0.

In Fig. 8 we show the optically-selected quasar RLF across the full luminosity and redshift range from our sample. The black circles denote the RLF determined using the $1/V_{\max}$ method, which is only possible for those detected above a certain flux-density threshold (we use $5\sigma$), whereas the lines and shaded regions show the full RLF distribution from the Bayesian modelling. Concentrating on the lowest redshift bin (top-left panel of Fig. 8), we find that the number density of radio-bright quasars increases with decreasing radio luminosity in all redshift bins, as expected.

We compare our inferred RLFs for optically-selected quasars with similar RLFs from the literature: Condon et al. (2013) and Kellermann et al. (2016). The optical data are the same volume-limited sample from SDSS’s DR7. The radio data are all from the VLA but each sample was observed with different configurations and depths. Our data are from FIRST, which was observed in the ‘B’ configuration, with a resolution of $5^{\prime\prime}$ and rms of 0.15 mJy, corresponding to a detection threshold of 1 mJy. The Condon et al. sample is from the NRAO VLA Sky Survey (NVSS; Condon et al. 1998) observed using the compact ‘D’ and ‘DnC’ configurations, with a resolution of $45^{\prime\prime}$ and rms of 0.45 mJy (a detection threshold of 2.4 mJy). Kellermann et al. observed a complete sub-sample of these quasars, over a reduced redshift range ($0.2<z<0.3$) at 6 GHz using the Karl G. Jansky Very Large Array (JVLA) in the ‘C’ configuration, with a resolution of $3.5^{\prime\prime}$ and rms as low as $6\mu$Jy for the fainter sources. In order to enable a direct comparison to the results in our lowest redshift bin, the 6 GHz luminosities of the Kellermann et al. sources are converted to 1.4 GHz luminosities using a spectral index of $\alpha=0.7$ and their number density is increased by $\log_{10}[\Phi/(\rm{Mpc}^{-3}\rm{mag}^{-1})]=0.2$ to correct for evolution (Condon et al. 2013) when comparing the RLF over the redshift range $0.2<z<0.3$ with that over $0.2<z<0.45$.

The RLF above the nominal 5-$\sigma$ threshold for our sample is in good agreement with both Condon et al. (2013) and Kellermann et al. (2016) RLFs between radio luminosities $25<\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]<26$ but is less consistent with the Condon et al. (2013) RLF towards the low-luminosity end of where we have direct detections ($23.6<\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]<25$). Furthermore, our RLF has large uncertainties above $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]\sim 26$. These are both likely due to the fact that only 7 of the 26 sources observed in NVSS are compact (Condon et al. 2013) and the rest are extended sources that have emission resolved out by FIRST (hence the sources occupy lower luminosity bins below $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]\sim 26$), and as such lead to the discrepancy with the Condon et al. (2013) study and reduce the numbers in the highest luminosity bins.

Each of the RLFs in this redshift bin also show a flattening in the number density between $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]\approx 25.5$ and $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]\approx 24$. Below $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]=24.4$ our RLF is higher than that of Condon et al. but still in good agreement with Kellermann et al.. The difference between our RLF and that from Condon et al. is most probably due to the difference in resolution of the radio data, which results in sources moving into lower-luminosity bins due to some emission being resolved out.

Given the likely underestimation of extended emission using the FIRST survey, we use the Condon et al. flux-densities for sources found in both NVSS and FIRST in the RLF fit (shown in Fig. 7 and Fig. 8). However, we note that extended emission may still be resolved out for sources below the flux-density limit but we have no way of estimating this. Although we could potentially use the NVSS data here, we would then have to deal with confusion issues due to the larger synthesised beam. We therefore continue to use the FIRST data, but the issue of extended emission should be borne in mind.

The reconstruction of the RLF below the detection threshold continues to follow the slope (of the double power-law) established from $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]\geq 24.4$, dropping at $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]\approx 22.4$. This therefore measures the RLF $\sim 2$ orders of magnitude below $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]=23.4$ ($=5\sigma_{n}$). The steep drop-off in the RLF at $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]<22$ is due to the optical limit of the quasar sample, meaning that there are no optically-selected quasars with $i<19.1$ contributing to this part of the RLF.

Comparing the reconstructed RLF to the Kellermann et al. (2016) individually-observed sources, we find that the two measurements are in complete agreement.

5.2 Higher-redshift bins

In Sec. 5.1 we demonstrated that the technique is able to reconstruct the RLF below the detection threshold in the lowest-$z$ sample, for which there are deeper radio data. In this section we present the result using our algorithm and the three models describing the RLF to the higher redshift bins. The relative evidence of the models for each redshift bin are shown in Table 4. The data prefer Model B (a double power-law for the low-luminosity sources) for all the redshift bins. The posterior distributions for the winning models are shown in the online version of the paper.

The optically-selected quasar RLF mirrors the general shape seen in the lowest redshift bin over all redshifts. In all cases we see that the bright-end of the RLF increases steeply as the radio luminosity decreases towards $L_{1}^{*}\sim 10^{25}$ W Hz*-1* and then turns over. Just below this luminosity we see the second (faint-end) double power law starting to dominate the RLF, where we find a steep increase as the radio luminosity decreases towards $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]\sim{23}$. Our reconstructed RLF also follows the 1/$V_{\rm{max}}$ points very well where we are able to measure them.

This flattening of the bright-end of the RLF and subsequent increase below $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]<26$ is also observed in optically-selected quasar RLFs studies (e.g. Condon et al., 2013; Kellermann et al., 2016; Hwang et al., 2018). A similar flattening is also observed in the RLF of other optically-selected-AGN samples (e.g. Rush et al., 1996; Padovani et al., 2015). A clear change of the slope in the number density is also observed in radio-selected AGN RLFs (e.g Willott et al. 2001, Smolčić et al. 2009, McAlpine et al. 2013). Indeed, our fitted values for $L_{1}^{*}$ are in good agreement with the RLF derived using the deep VLA-3 GHz survey from Smolčić et al. (2017).

At radio luminosities below where the flattening takes place, the reconstructed RLF steeply increases towards lower luminosities, with a slope established above $5\sigma_{n}$ for all redshift bins. The preferred model for all redshift bins is Model B (Table 4, the double power-law). With this model, the RLF in all redshift bins has a peak at $L_{2}^{*}$, dropping below $L_{2}^{*}$. This fall-off is due to the hard absolute magnitude cut-off in the parent sample, and essentially means that there is no significant evidence for any radio continuum emission from our quasar sample below $L_{2}^{*}$.

6 Discussion

The definition of radio-loudness varies in the literature, as some objects can be classified as ‘radio-quiet’ in one definition and ‘radio-loud’ in another, e.g. either by considering the ratio of optical to radio emission (e.g. Kellermann et al., 1989) or by just using a radio lumninosity threshold (e.g. Miller et al., 1990). In this paper, we do not explicitly classify our quasars as radio-loud or radio-quiet, instead using the shape of the RLF to infer where these populations dominate. In all of our RLFs (Fig. 8) there is a clear change in behaviour at or around $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]\sim 25$. We define the ‘radio-loud’ population as the quasars that are described by a bright-end double power-law (parameters with subscript ‘1’ in the modelling). The faint end (radio-quiet quasars) is parameterised by the power-law, double power-law or log-normal function. For this study all redshift bins had the double power-law as the winning model (Table 4).

6.1 Radio-loud quasars

The radio emission from radio-loud quasars are powered by processes associated with the accretion on to the central supermassive black hole. Falling within the AGN orientation-based unification model (e.g. Barthel, 1989; Antonucci, 1993; Urry & Padovani, 1995), these radio-loud objects have been shown to require a supermassive black-hole of mass $M_{\rm BH}>10^{8}$ M*⊙* (McLure & Jarvis, 2004), whereas their radio-quiet counterparts can have lower-mass black holes. By integrating under the two double power-law models, representing the bright- and faint-end of the RLF, we find that the radio-loud fraction of quasars make up $\approx$ 10 per cent of the total quasar population in our sample at $z>0.7$ (Table 1). However, we find that the radio-loud fraction drops to $\approx$ 7 and 4 per cent of the total quasars in the two lowest redshift bins (Table 1). This lower fraction of radio-loud quasars towards lower redshifts reflects the fact that we have a much fainter optical magnitude limit at low redshift, and if radio-loudness is linked to the combination of accretion rate and black-hole mass, then lower-optical luminosity quasars are more likely to be radio quiet. These fractions are in line with previous studies of radio-loud and radio-quiet quasars with a variety of classification schemes (e.g. White et al., 2007; Cirasuolo et al., 2005; Baloković et al., 2012).

However, one of the differences is that we actually find a much more pronounced flattening than the studies based purely on radio-selected samples (e.g. Willott et al., 2001; Smolčić et al., 2009; McAlpine et al., 2013). One reason for this could be that there is a real difference in the physical properties that generate radio emission in optically-selected quasars compared to the more-general population of radio-selected AGN. We also cannot rule out the possibility of the optical selection creating a bias in the RLF that artificially flattens, or decreases, the bright-end of the RLF below $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]\sim 26$. However, we have been conservative with our optical selection, ensuring that the quasar sample is complete across the full width of all redshift bins. We cannot rule out incompleteness due to the colour selection within the SDSS sample, but we would not expect this to have a significant effect in individual, relatively narrow, redshift bins. A possible explanation could be due to our sample becoming incomplete in terms of the RLF based on the optically-selected sample. This could arise if there is a correlation between the optical emission in these quasars and their radio emission.

Several authors have investigated the link between optical emission and radio emission from quasars (e.g. Serjeant et al., 1998; White et al., 2007; White et al., 2017), finding evidence for a correlation. However, one has to be careful when measuring correlations between flux-limited samples. Therefore, in Fig. 8 we show the radio luminosity where we expect the optical flux limit to start imposing incompleteness on the RLF, based on the absolute magnitude limits shown in Table 1. For this we use the relation between optical luminosity and the star-formation subtracted radio luminosity found by White et al. (2017), from their radio-quiet quasar sample at $z\sim 1$. We also show the radio luminosity limit based on the White et al. (2017) optical luminosity versus total radio-luminosity, for completeness. One can see that the radio luminosity at which the optical selection may lead to incompleteness in the RLF is just above the radio luminosity at which the drop in the RLF occurs. This supports the argument that the drop/turnover at low luminosities is caused by lack of fainter optical quasars. However, there is significant scatter in the White et al. (2017) optical-radio correlation of around 1 order of magnitude in radio luminosity for a given optical luminosity. Therefore, it is certainly possible that some of the flattening could arise from incompleteness in the RLF due to the optical magnitude limit. To test this we increased the optical magnitude limit for our sample in each redshift bin in order to check if the flattening or downturn becomes more prominent. In all bins the turnover (i.e. the value of $\beta_{1}$) became more prominent. We therefore suggest that at least some of the flattening is due to incompleteness introduced by the optical magnitude limit of the parent sample, although we note that the uncertainties increase due using a smaller quasar sample when a higher optical-luminosity threshold is imposed..

6.2 Radio-quiet quasars

This population makes up about 92 per cent (Table 1) of the quasar population in our sample, but the origin of the radio emission is not well understood. Our reconstructed RLFs increase steeply below $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]\sim{24.5}$. This steepening could be attributed to an increasing contribution from SF in the host galaxy (e.g. Terlevich et al. 1987, 1992; Padovani et al. 2011; Kimball et al. 2011; Bonzini et al. 2013; Condon et al. 2013; Kellermann et al. 2016; Stacey et al. 2018; Gürkan et al. 2018) or is AGN-related with a different scaling relation or different emission associated with the AGN (Herrera Ruiz et al. 2016, Zakamska et al. 2016, White et al. 2015, 2017, Hartley et al. 2019) compared to their radio-loud counterparts. Although, we note that the steepening is significantly less pronounced in the two lowest-redshift bins, which may indicate that the optical magnitude limit may play a role in creating an artificially-steepening slope in the observed RLF. In such a case, the distinction between radio-loud and radio-quiet would become more difficult, with evidence that the population has a more continuous distribution (e.g. Lacy et al., 2001; Gürkan et al., 2019).

Kimball et al. (2011) suggested that the ‘bump’ observed at $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]\approx 22.7$ in their low-$z$ volume-limited sample (Fig. 8) corresponds to star-forming galaxies. Kellermann et al. (2016) tested their hypothesis by using mid-infrared data from WISE to search for a correlation between the 22 $\mu$m and 6-GHz flux-densities, which is a characteristic of the radio–far-infrared correlation. However, they found no strong correlation and so suggest that the 22 $\mu$m fluxes do not only measure SF but can also be contaminated by warm dust heated by the AGN (Polletta et al. 2010). Coziol et al. (2017) tested the SF hypothesis by also matching the Kellermann et al. sources against WISE. They found counterparts for all but 7 sources, created a new diagnostic plane based on WISE colours (Coziol et al. 2015), and found that: (i) there is no separation between the radio-quiet and radio-loud quasars in the colour distribution, and (ii) the majority of the Kellermann et al. 2016 quasars (and our lowest $z$ sample) have low star-formation rates.

White et al. (2015) used deep optical and near-infrared data to identify a sample of quasars across a range of redshifts, and conducted a stacking experiment using deep VLA 1.4 GHz data. Comparing their quasar sample with multiple galaxy samples (matched by stellar mass, having assumed the black-hole mass of the quasar), they provided evidence that the radio emission from these quasars – which lie at much higher redshift but cover similar optical luminosities as the Kellermann et al. (2016) sample – predominantly arises from accretion-related activity. Furthermore, by comparing the star-formation rates using far-infrared data of a randomly selected subset of a volume-limited quasar sample at $0.9<z<1.1$, White et al. (2017) showed that the radio emission from star formation is sub-dominant.

The only evidence in our modelled RLFs for star-formation contributing to the radio emission in quasars comes from the observed strong steepening of the RLF towards low luminosities, below the nominal 5$\sigma$ detection threshold at $z>0.7$. However, where our optically-selected quasar sample contains the lowest-luminosity quasars ($z<0.7$), the evidence for this steepening is weaker. On the other hand, comparing the observed upturn in the quasar RLF with the star-forming galaxy RLF at $z>0.8$ from Novak et al. (2018), we find that the steepening occurs at approximately the same radio luminosity that the star-forming galaxies dominate over AGN in radio-selected surveys. This strengthens the suggestion that star formation plays an important role at these low radio luminosities. Indeed, this was used as evidence in favour of the star-formation becoming the dominant contribution to the radio luminosity in this regime by Kimball et al. (2011) and Condon et al. (2013). Nevertheless, it is clear that the RLF is a relatively blunt tool for disentangling the dominant contribution to the low-luminosity radio emission in quasars. A more productive route may be to explore the bivariate optical and radio luminosity function for quasars (e.g. Singal et al., 2011), where the optical selection is naturally accounted for and models that link the optical and radio emission could be incorporated.

A more direct method would be to use high-resolution radio data that can resolve any star formation on the scale of the host galaxy. The VLA has the potential to do this, but would need to move towards a frequency of 6 GHz to achieve the required resolution of $\sim 0.5\arcsec$ for the vast majority of quasars that lie at $z>0.5$. Given the typical spectral index of the synchrotron radiation from both star-formation and AGN-associated emission of $\alpha\sim 0.7$, this would then require longer integration times and may also suffer from contamination from free-free emission, making the results more difficult to interpret. eMERLIN has the potential to carry out similar resolution studies at lower frequencies (e.g. Guidetti et al., 2013; Radcliffe et al., 2018; Jarvis et al., 2019). In the future, the Square Kilometre Array (e.g. Jarvis & Rawlings, 2004; Smolcic et al., 2015; McAlpine et al., 2015) would be able to carry out high-resolution studies to much deeper levels at a range of frequencies and thus help make great strides in our understanding of the dominant radio emission mechanism in radio-quiet quasars.

7 Conclusions

1. [i]

We have built on the work of Roseboom & Best (2014) and Zwart et al. (2015b) by fitting directly quasar radio luminosity functions below the radio detection threshold using a Bayesian stacking approach (bayestack). We tested the technique by fitting three models to mock SKADS simulation catalogues (Wilman et al. 2008; 2010), with random Gaussian noise of 150 $\mu$Jy added. We successfully recovered the SKADS RLF over three orders of magnitude below the $5\sigma$ detection threshold. We ran further tests using mock catalogues with 15 $\mu$Jy Gaussian noise and as expected reconstructed a better-constrained RLF with respect to the true SKADS RLF. 3. 2.

We used FIRST radio flux-densities extracted at the positions of optical quasars from a uniformly-selected (homogeneous) sample of SDSS DR7 divided into seven volume-limited redshift bins. We parameterised the high-luminosity RLF using a double power-law. For our lowest-$z$ sample we found that the $1/V_{\rm{max}}$ and double power-law RLF for luminous sources is in agreement with that from Kellermann et al., but is marginally inconsistent with that of Condon et al. at the luminous and faint ends of the detected RLF. Some of the difference at the faint end is likely due to the different resolution of NVSS and FIRST. In the other redshift bins, we find that each of the bright ends of the RLFs, which broadly represent radio-loud quasars, are well described by a double power-law. This double power-law generally flattens towards low luminosities. A similar drop/flattening is observed for extremely-red quasars (Hwang et al., 2018) and in AGN RLFs (e.g Smolčić et al. 2009). We suggest that some of this flattening could also be due to the optical flux limit of the sample reducing the number of quasars that could contribute radio data to these radio luminosities, although this would need to be tested thoroughly with a deeper optical selection or by considering a bivariate model of the optical and radio luminosity functions. 4. 3.

With bayestack we probe the RLF down to 2 orders of magnitude below the detection threshold of FIRST (1 mJy). The difference in how deep we can probe compared to the simulation is related to the lack of low luminosity radio sources in the sample as compared to SKADS. At low redshift ($z<0.7$) we see a continuous distribution from the bright to the faint end of the RLF, whereas at $z>0.7$, the RLF steeply increases towards fainter luminosities $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]<{24.5}$. This could be due to the source population changing or due to the biased flattening of the RLF because of the optical flux limit described previously. We note, however, that the steep increase coincides with the measured steepening in the RLF from radio-selected samples of star-forming galaxies (e.g. Novak et al., 2018). In order to resolve whether this steepening is indeed due to star formation, higher-resolution radio imaging would be ideal, in order to resolve the radio emission from star formation in the host galaxy. 5. 4.

Finally, the RLF peaks between $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]\sim{22.5}$ and $\log_{10}[L_{1.4}/{\rm W\,Hz}^{-1}]\sim{24}$ (depending on the redshift) and drops rather abruptly after that. This is due to the parent sample containing no quasars that are generating radio emission below this luminosity.

Acknowledgments

We thank the referee for the helpful comments that have contributed to improve this paper. EM, MGS, MJJ and SVW acknowledge support from the South African Radio Astronomy Observatory (SARAO). EM and MGS also acknowledge support from the National Research Foundation (Grant No. 84156). JZ is thankful for a South African Square Kilometre Array Research Fellowship. We are grateful for valuable contributions from Matt Prescott, Imogen Whittam, Margherita Molaro, Jośe Fonseca and Brandon Engelbrecht. We would like to acknowledge the computational resources of the Centre for High Performance Computing.

Appendix

Fig. A0 show the 1-D and 2-D posterior distributions for all of the winning models for each redshift bin. The 1-D posterior distribution is the marginalization of each parameter shown at the end of each row. The parameters have well-defined peaks, except for the boundary parameters ($L_{\rm{min}_{1}}$, $L_{\rm{max}_{1}}$, $L_{\rm{min}_{2}}$ and $L_{\rm{max}_{2}}$) and the second slope $\beta_{2}$ (the faint-end slope for the faint-quasar function in Models A and B), which are not well constrained. The upper limit of the fitted $L_{\rm{max}_{1}}$ is unconstrained or (even) truncated, but this does not significantly affect the fit in our areas of interest.