THE RADIO AND OPTICAL LUMINOSITY EVOLUTION OF QUASARS. II. THE SDSS SAMPLE

Quasars are distant active galactic nuclei (AGNs) with emission seen across the electromagnetic spectrum, from radio to X-rays. The optical emission from quasars is thought to be dominated by the radiation from the accretion disk around supermassive black holes, while the radio emission is dominated by the plasma outflowing from the black hole/accretion disk systems. Because of this, different but complementary information can be gathered from both photon energy ranges regarding the cosmological evolution of AGNs. It is therefore important to determine in detail the redshift evolutions of quasars in both radio and optical regimes.

In a previous paper (Singal et al. 2011, hereafter QP1), we explored the luminosity evolutions and radio-loudness distribution of quasars with a data set consisting of the overlap of the FIRST (Faint Images of the Radio Sky at Twenty cm) radio survey with the Automatic Plate Measuring Facility catalog of the Palomar Observatory Sky Survey (POSS-I), as presented by White et al. (2000). In this paper, we present the results from a much larger data set consisting of the overlap of FIRST with the Sloan Digital Sky Survey (SDSS) Data Release 7 (Abazajian et al. 2009) quasar catalog, which contains a factor of ∼10 more objects and spans redshifts from 0.065 to 5.46.

In the literature, the evolution of the quasar luminosity function (LF) has been described not only for optical and radio luminosities but also for X-ray, infrared, and bolometric luminosities (e.g., Ueda et al. 2003; Richards et al. 2006; Matute et al. 2006; Hopkins et al. 2007; Croom et al. 2009). The shape of the LF and its evolution are usually obtained from a flux-limited sample in the a waveband f_a > f_{m, a} with f_{m, a} denoting the flux limit and L_a = 4 π d²_Lf_a/K_a, where d_L(z) is the luminosity distance and K_a(z) stands for the K-correction. For a pure power-law emission spectrum of index ε_a defined as $f_a \propto \nu ^{-\varepsilon _a}$, one has $K_a(z) = (1+z)^{1-\varepsilon _a}$. This simple form may be modified for optical data by the presence of emission lines, as in this work.

In general, the determination of the LF and its evolution requires analysis of the bi-variate distribution Ψ_a(L_a, z). The first step of the process should be the determination of whether the variables of the distributions are correlated or are statistically independent. A correlation between L_a and z is a consequence of luminosity evolution. In the case of quasars with the optical and some other band luminosity, we have at least a tri-variate function. One must determine not only the correlations between the redshift and individual luminosities (i.e., the two luminosity evolutions) but also the possible intrinsic correlation between the two luminosities, before individual distributions can be determined. Knowledge of these correlations and distributions is essential not only for constraining robustly the cosmological evolution of active galaxies but also for interpretation of related observations, such as the extragalactic background radiation (e.g., Singal et al. 2010; Hopkins et al. 2010).

A related question is the distribution of the "radio-loudness parameter," R = L_rad/L_opt, for the quasar population, defined as the ratio of the 5 GHz radio to 2500 Å  optical luminosity spectral densities, and the distinction between so-called radio-loud (R > 10; RL for short) and radio-quiet (R < 10; RQ for short) quasars. Weak hints of a bimodality in the distribution of the radio-loudness parameter described by Kellerman et al. (1989) suggested that log R = 1 could be chosen as the radio-loud/radio-quiet demarcation value. Using this value for the division between RL and RQ quasars, the differences between the two classes have been investigated, including the possibility of distinct cosmological evolution of the RL and RQ populations (e.g., Miller et al. 1990; Goldschmidt et al. 1999; Jiang et al. 2007). Still, the more recent analyses of different samples of objects reported in the literature so far gave rather inconclusive results on whether any bimodality in the distribution of the radio-loudness parameter for quasars is inherent in the population (see Ivezic et al. 2002, 2004; Cirasuolo et al. 2003; Zamfir et al. 2008; Kimball et al. 2011; Mahony et al. 2012). In QP1 we found no evidence for a bimodality in R in the range −1 < log R < 4.

In addition to the above-cited works, there have been many papers dealing with this ratio and RL versus RQ issue, as well as luminosity ratios at other wavelengths: IR/radio, optical/X-ray, etc. However, to the best of our knowledge, except for QP1, none of these works have addressed the correlations between the radio and optical luminosities, which is necessary for such an analysis. Additionally, in general they have not concentrated on obtaining the intrinsic—as opposed to the raw observed—distribution (and/or evolution) of the ratio, which is related to the tri-variate LF Ψ(L_opt, L_rad, z) by⁵

$$\begin{eqnarray} G_R(R,z) &=& \int _0^{\infty } { \Psi (L_{\rm opt}, R \,\, L_{\rm opt}, z) \, L_{\rm opt} \, dL_{\rm opt} } \nonumber \\ &=& \int _0^{\infty } { \Psi \left({L_{\rm rad} \over R}, L_{\rm rad}, z \right) \, L_{\rm rad} \, {{ dL_{\rm rad} } \over {R^2}} }. \end{eqnarray} \tag{ 1 }$$

In Appendix A of QP1 we showed how, even in the simplest cases, the observed distributions can be very different from the intrinsic ones. Thus, for determination of the true distributions the data truncations must be determined and the correlations between all variables must be properly evaluated.

Efron & Petrosian (1992, 1999, EP for short) developed new methods for determination of the existence of correlation or independence of variables from a flux-limited and more generally truncated data set, which were further expanded in QP1. Our aim in this paper is to take all the selection and correlation effects into account in determination of the true evolution of optical and radio luminosities and their ratio and to find their distributions, using the larger SDSS DR7 QSO × FIRST data set.

In Section 2 we describe the radio and optical data used. Section 3 contains a general discussion of luminosity evolution and the sequence of the analysis. In Section 4 we apply the EP method to achieve the luminosity–redshift evolutions and the correlation between the luminosities. We determine the density evolution in Section 5 and the local LFs in Section 6. A discussion of the radio-loudness distribution is presented in Section 7. In Section 8 we investigate some of the assumptions used, and Section 9 contains a discussion of the results. This work assumes the standard ΛCDM cosmology throughout, with H₀ = 71 km s⁻¹ Mpc⁻¹, Ω_Λ = 0.7, and Ω_m = 0.3.

In order to evaluate the luminosity evolution in both radio and optical and to separate and compare these effects, we require a data set that has both radio and optical fluxes to reasonable limits across a broad range of redshifts that contains a significant number of both RL and RQ objects. The overlap of the FIRST bright quasar radio survey (Becker et al. 1995; White et al. 1997) with the SDSS DR7 quasar catalog (DR7Q; Schneider et al. 2010) can form such a data set.

The SDSS DR7 quasar catalog has a limiting i-band magnitude of 21, or f_{m, i} = 0.015 mJy, and contains over 105,000 objects with redshifts ranging from 0.065 to 5.46 over 9380 deg². We work with a more restrictive i-band optical magnitude limit of 19.1 (f_{m, i} = 0.083 mJy) to form a parent optical catalog that has a smoother redshift distribution with a reduced bias toward objects with z > 2 and with a uniform limiting flux for every redshift (see Section 6 of DR7Q). The i < 19.1 catalog consists of 63,942 objects, with a maximum redshift of 4.98, shown in Figure 1.

Zoom In Zoom Out Reset image size

In the SDSS DR7 catalog, objects are identified as quasar candidates either if they have the requisite optical colors or if they have a radio match within 2''. Only 1% of the quasars in the catalog were selected from the candidate list for spectroscopic follow-up based on the latter criterion. We also note the presence of an upper limit i-band magnitude of 15.0 for inclusion in the catalog. However, this criterion is not completely rigorous, as mentioned in DR7Q. Further properties of the DR7 quasar catalog, such as black hole masses obtained from line emission measurements, are presented in Shen et al. (2011).

The FIRST survey, carried out with the Very Large Array (VLA) in B configuration between 1995 and 2004, has a limiting peak pixel 1.4 GHz flux of 1 mJy and contains 816,000 sources over 9500 deg². Different criteria can be used to determine radio and optical matches to form a joint SDSS DR7 QSO × FIRST catalog. Jiang et al. (2007, hereafter J07), in their analysis of SDSS quasars, used a matching radius of 5'' for single radio matches and 30'' for multiple matches to an optical source. We construct a fiducial catalog using these radius matching criteria, but, unlike J07, we do not include quasars with no radio detection (see the discussion in Section  9.1).

This fiducial catalog contains 5677 objects ranging in redshift from 0.064 to 4.82 and with log R ranging from −0.47 to 4.44. A total of 4327 of the objects in the fiducial catalog are single matches (which J07 calls "Fanaroff–Riley type I sources," FR1s, perhaps incorrectly; see the discussion in Section 9.1), and 1350 are multiple matches (which J07 calls "Fanaroff–Riley type II objects," or FR2s for short). This sample seems to be skewed toward RL as opposed to RQ quasars, with 4543 (80%) having R > 10. It should be noted that J07 used the SDSS Data Release 3, so our realization with the same matching criteria, even with the additional i < 19.1 magnitude cut and no objects without radio detections, results in more objects (J07 have 2566 objects).

We also construct a catalog with a universal 5'' matching criterion, which we feel is more appropriate, as discussed in Section 9.1. In this catalog there are 5445 objects with only 25 multiple radio matches to an optical source. We will refer to this as the canonical sample, and base our analysis primarily on it, while comparing to results we obtain with the fiducial sample, which as discussed below are quite similar. In the canonical sample 4466 (82%) have R > 10. The major changes from the fiducial catalog to the canonical one are that some optical objects with multiple radio matches in the fiducial sample are reduced to having only single matches in the canonical sample, with a corresponding reduction in radio flux, while a small number of objects drop out entirely, having had multiple matches within 30'' but none within 5''. The fiducial sample and canonical sample are nearly identical in their redshift distribution, and the average radio luminosity differs only slightly, being 10% lower in the canonical sample.

As shown in Richards et al. (2006, hereafter R06; e.g., Figure 8 of that work), the SDSS optical quasar sample contains significant biases in the redshift distribution due to emission-line effects and the differing flux limits at z > 2. We have reduced the latter effect by using only the universally flux-limited i < 19.1 sample, and the former by adopting the full ε_opt = 0.5 power-law continuum plus emission-line K-corrections presented in R06 and discussed in Section 5 and Table 4 of that work. The methods of this work can then account for any bias resulting from emission-line effects, as long as they are included in the conversions from luminosity to flux, i.e., in the K-corrections. For the radio data we assume a power law with ε_rad = 0.6. Figures 2 and 3 show the radio and optical luminosities versus redshifts of the quasars in the canonical constructed SDSS × FIRST catalog. Figure 4 shows the radio-loudness parameter R versus redshift for the canonical data set, while Figure 5 shows the radio luminosity versus the optical luminosity for different redshift bins.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1. Luminosity and Density Evolution

The LF gives the number of objects per unit comoving volume V per unit source luminosity, so that the number density is dN/dV = ∫dL_aΨ_a(L_a, z) and the total number is N = ∫dL_a ∫dz (dV/dz) Ψ_a(L_a, z). To examine luminosity evolution, without loss of generality, we can write an LF in some waveband a as

$$\begin{equation} \Psi _{a}(L_a,z) = \rho (z)\,\psi _a\big(L_{a}/g_{a}(z), \eta _a^j\big)/g_a(z), \end{equation} \tag{ 2 }$$

where g_a(z) and ρ(z) describe the luminosity evolution and comoving density evolution with redshift, respectively, and η^j_a stands for parameters that describe the shape (e.g., power-law indices and break values) of the a-band LF. In what follows we assume a non-evolving shape for the LF (i.e., η^j_a = const, independent of L and z), which is a good approximation for determining the global evolutions. Once the luminosity evolution g_a(z) is determined using the EP method, we can obtain the mono-variate distributions of the independent variables L'_a = L_a/g_a(z) and z, namely, the density evolution ρ(z) and "local" LF ψ_a. The total number of observed objects is then

$$\begin{equation} N_{{\rm tot}} = \int _0^{z_{{\rm max}}} dz \int _{L_{\rm min}(z)}^\infty { dL_{a} \, \rho (z) \, {dV \over dz} \, { {\psi _a\!\left(L_{a}/g_{a}(z)\right)} \over {g_a(z)} } }. \end{equation} \tag{ 3 }$$

We consider this form of the LF for luminosities in different bands, allowing for separate (optical and radio) luminosity evolutions.

In the previous work (QP1) with the White et al. data set, we assumed a simple power law for the luminosity evolution

$$\begin{equation} g_{\rm a}(z)=(1+z)^{k_{\rm a}}. \end{equation} \tag{ 4 }$$

However, with the SDSS data set, which contains objects past z = 4 and many more objects past z = 3, we need to use a more complicated parameterization:

$$\begin{equation} g_{\rm a}(z)={ {(1+z)^{k_{\rm a}}} \over {1+\big({{1+z} \over {z_{cr}}}\big)^{k_{\rm a}}} }. \end{equation} \tag{ 5 }$$

With either definition of g_a(z) that gives g_a(0) = 1, the luminosities L'_a refer to the de-evolved values at z = 0, hence the name "local."

We have verified the form of Equation (5) as a good fit to the evolution by considering the optical evolution in the large parent optical-only sample in narrow redshift bins with the appropriate simple form of Equation (4) and then verifying the form of the cumulative evolution over a larger range. We have determined the appropriate value for z_cr and confirmed that the same exponent is appropriate in the numerator and denominator, using the large set of optical-only data and the methods of Section 4.1. With the exponents fixed to be the same, the optimal value of z_cr is determined to be 3.7 ± 0.3, by considering the 1σ range of the best-fit evolution while letting that numerical factor vary. We would expect that if the truncations and correlations have been properly accounted for by the methods of this work, g_opt(z) as determined from the simultaneous radio and optical data set should match that as determined from the parent optical-only data set, which is shown to be the case (see Section 4).

We discuss the determination of the evolution factors g_a(z) with the EP method, which in this parameterization becomes a determination of k_a, in Section 4. The density evolution function ρ(z) is determined by the method shown in Section 5. Once these are determined, we construct the local (de-evolved) LF ψ_a(L'_a), shown in Section 6.

3.2. Joint LFs

In general, determination of the evolution of the LF of extragalactic sources for any wavelength band except optical involves a tri-variate distribution because redshift determination requires optical observations, which introduces additional observational selection bias and data truncation. Thus, unless redshifts are known for all sources in a radio survey, we need to determine the combined LF Ψ(L_opt, L_rad, z) from a tri-variate distribution of z and the fluxes in the optical and radio bands. If the optical and radio luminosities were statistically independent variables, then this LF would be separable, Ψ(L_opt, L_rad, z) = Ψ_opt(L_opt, z) × Ψ_rad(L_rad, z), and we would be dealing with two bi-variate distributions.

However, the luminosities in the different wavebands may be correlated. The degree and form of the correlation must be determined from the data, and as described below, the EP method allows us to determine whether any pair of variables are independent or correlated. Once it is determined that the luminosities are correlated (see Section 4.2), the question must be asked whether this luminosity correlation is intrinsic to the population or induced in the data by flux limits and/or similar luminosity evolutions with redshift. Determination of which is the case is quite intricate, as discussed in Appendix B of QP1, and has not been explored much in the literature. This will be the subject of a future work. Here we will simply consider both possibilities.

At one extreme, if the luminosity correlation is intrinsic and not induced, one should seek a coordinate transformation to define a new pair of variables that are independent. This requires a parametric form for the transformation. One can define a new luminosity that is a combination of the two, which we call a "correlation-reduced radio luminosity," L_crr = L_rad/F(L_opt/L_fid), where the function F describes the correlation between L_rad and L_opt and L_fid is a fiducial luminosity taken here⁶ to be 10²⁸ erg s⁻¹ Hz⁻¹. For the correlation function we will assume a simple power law

$$\begin{equation} L_{\rm crr} = {{L_{\rm rad}} \over {(L_{\rm opt}/L_{\rm fid})^{\alpha }}}, \end{equation} \tag{ 6 }$$

where α is a bulk power-law correlation index to be determined from the data. This is essentially a coordinate rotation in the log–log luminosity space. As shown in Section 4 below, EP also prescribe a method to determine a best-fit value for the index α that orthogonalizes the new luminosities. Given the correlation function, we can then transform the data (and their truncation) into the new independent pair of luminosities (L_opt and L_crr).

At the other extreme, if the correlation between the luminosities is entirely induced, then the LFs are separable as described above and the analysis can proceed from there. As mentioned, we will consider both possibilities here. It turns out that the results obtained in both cases are very similar, implying that the luminosity–redshift correlations in the analyzed sample are much stronger than any intrinsic luminosity–luminosity correlations (if present). Further mechanics of the procedure for intrinsically correlated luminosities is discussed in the Appendix.

In this work we are also interested in the distribution of radio-loudness parameter R, specifically its de-evolved value R' = L'_rad/L'_opt. In the case of uncorrelated radio and optical luminosities or an induced, non-intrinsic radio–optical luminosity correlation, the local distribution of R', denoted hereafter as G(R'), can be constructed from the local optical and radio LFs as in Equation (1), with the redshift evolution function

$$\begin{equation} g_{\rm R} = { {g_{\rm rad}} \over {g_{\rm opt}} }. \end{equation} \tag{ 7 }$$

In the case of an intrinsic radio–optical luminosity correlation, $G_{\rm R^{\prime }}$ and the evolution g_R can be constructed by Equations (A3) and (A4) given in the Appendix.

We now describe results obtained from the use of the procedures described in Section 3 on the data described in Section 2. Here we first give a brief summary of the algebra involved in the EP method. We generally follow the procedures described in more detail in QP1. This method uses the Kendall tau test to determine the best-fit values of parameters describing the correlation functions using the test statistic

$$\begin{equation} \tau = {{\sum _{j}{(\mathcal {R}_j-\mathcal {E}_j)}} \over {\sqrt{\sum _j{\mathcal {V}_j}}}} \end{equation} \tag{ 8 }$$

to test the independence of two variables in a data set, say, (x_j, y_j) for j = 1, ..., n. Here $\mathcal {R}_j$ is the dependent variable (y) rank of the data point j in a set associated with it. For untruncated data (i.e., data truncated parallel to the axes) the set associated with point j includes all of the points with a lower independent variable value (x_k < x_j). If the data are truncated, one must form the associated set consisting only of those points of lower independent variable (x) value that would have been observed if they were at the x value of point j given the truncation. As an example, if we have one-sided truncations as in Figures 2 and 3, ignoring the upper optical flux limit for the moment, then the associated set A_j = {k: y_k > y_j, y⁻_k < y_j}, where y⁻_k is the limiting y value of data point j (see EP for a full discussion of this method).

If (x_j, y_j) were independent, then the rank $\mathcal {R}_j$ should be distributed uniformly between 0 and 1 with the expectation value and variance $\mathcal {E}_j=(1/2)(n+1)$ and $\mathcal {V}_j=(1/12)(n^{2}+1)$, respectively, where n is the number of objects in object j's associated set. Independence is rejected at the m σ level if | τ | > m. To find the best-fit correlation, the y data are then adjusted by defining y'_j = y_j/F(x_j) and the rank test is repeated, with different values of parameters of the function F.

4.1. Optical-only Data Set Luminosity–Redshift Correlation

With the optical-only data set (here we use the i < 19.1 set described in Section 2), determination of the luminosity–redshift correlation function g_opt(z) by Equation (5) reduces to determination of the value of the index k_opt. As evident from Figure 1, the L_opt − z data are heavily truncated due to the flux limits of SDSS. The associated set for object j includes only those objects that are sufficiently luminous to exceed the optical flux minimum for inclusion in the survey if they were located at the redshift of the object in question, i.e., L_k ⩾ L_min(z_j). As the values of k_opt are adjusted and L_opt is scaled by g_opt(k, z), the luminosity cutoff limits for a given redshift are also scaled by g_opt(k, z).

Figure 6 shows the value of the test statistic τ for a range of values of k_opt given the parameterization of Equation (5). It is seen that the best-fit value of k_opt is 3.3 with a 1σ range of 0.1. This indicates that quasars have strong optical luminosity evolution, i.e., the LF shifts to higher luminosities at larger redshifts.

Zoom In Zoom Out Reset image size

We have performed an optimization procedure to determine the best-fit value of z_cr, the numerical constant in the denominator of Equation (5), as well as whether the same exponent is appropriate for the numerator and denominator. This was done by varying these parameters and determining the best-fit range of values in the same way as described above. We find z_cr = 3.7 ± 0.3.

In this analysis we have ignored the upper optical flux limit of SDSS quasars corresponding roughly to an i-band magnitude of 15.0 as discussed in Section 2. The reason for this is that data truncations are only consequential in this analysis if the truncation is actually denying the sample data points that would otherwise be there. As can be seen in Figure 1, there are very few objects approaching this upper truncation limit, indicating that the truncation is not consequential, i.e., it does not appreciably alter the sample from the underlying population.

4.2. Radio–Optical Luminosity Correlation

Turning to the combined radio and optical data set, we will determine the correlation between the radio and optical luminosities. The radio and optical luminosities are obtained from radio and optical fluxes from a two-flux-limited sample so that the data points in the two-dimensional flux space are truncated parallel to the axes that we consider to be untruncated. Since the two luminosities have essentially the same relationship with their respective fluxes, except for a minor difference in the K-correction terms, we can consider the luminosity data to also be untruncated in the L_opt–L_rad space. In that case the determination of the associated set is trivial and one is dealing with the standard Kendall tau test. Assuming the correlation function between the luminosities F(x) = x^α, we calculate the test statistic τ as a function of α. Figure 7 shows the absolute value of the test statistic τ versus α, from which we get the best-fit value of α = 1.2 with 1σ range ±0.1. The result is similar whether the canonical or fiducial combined data set is used. As expected, α is near unity, and this result is compatible with that obtained in QP1 (α = 1.3 ± 0.2).

Zoom In Zoom Out Reset image size

As discussed in Section 3.2, this correlation may be inherent in the population or may be an artifact of the flux limits and wide range of redshifts. The general question of determining whether an observed luminosity correlation is intrinsic or induced will be explored in a future work. For this analysis going forward, we will consider both possibilities. We see that it does not make a significant difference for the major conclusions of this work.

4.3. Joint Data Set Luminosity–Redshift Correlations

The basic method for determining simultaneously the best-fit k_opt and k_rad, given the evolution forms in Equation (5), is the same as described in Section 4.1, but in this case the procedure is more complicated because we now are dealing with a three-dimensional distribution (L_rad, L_opt, z) and two correlation functions (g_rad(z) and g_opt(z)).

Since we have two criteria for truncation, the associated set for each object includes only those objects that are sufficiently luminous in both bands to exceed both flux minima for inclusion in the survey if they were located at the redshift of the object in question. Consequently, we have a two-dimensional minimization problem, because both the optical and radio evolution factors, g_opt(z) and g_rad(z), come into play, as the luminosity cutoff limits for a given redshift are also adjusted by powers of k_opt and k_rad.

We form a test statistic $\tau _{\rm comb} = \sqrt{\vphantom{A^A}\smash{\hbox{$\tau _{\rm opt}^2 + \tau _{\rm rad}^2$}}}$, where τ_opt and τ_rad are those evaluated considering the objects' optical and radio luminosities, respectively. The favored values of k_opt and k_rad are those that simultaneously give the lowest τ_comb, and, again, we take the 1σ limits as those in which τ_comb < 1. For visualization, Figure 8 shows a surface plot of τ_comb. Figure 9 shows the best-fit values of k_opt and k_rad taking the 1 and 2σ contours, for the canonical radio–optical data set. We have verified this method with a simulated Monte Carlo data set as discussed in QP1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We see that positive evolution in both radio and optical wavebands is favored. The minimum value of τ_comb favors an optical evolution of k_opt = 3.5 and a radio evolution of k_rad = 5.5, with a very small uncertainty at the 1σ level. The fact that k_opt as determined here from the combined data set is equal to that determined from the much larger optical-only data set (see Section 4.1 and Figure 6) indicates that the truncations have been properly handled.

In the above analysis we have assumed sharp truncation boundaries and that the data are complete above the boundaries. As discussed in Section 8, this may not be the case for the FIRST radio data. If we restrict the sample to only data with f_rad > 2 mJy, the favored optical evolution range for the canonical data set lowers slightly to k_opt = 3.0 ± 0.5 and the favored radio evolution range lowers slightly to k_rad = 5.0 ± 0.25. These overlap with the results considering the entire combined canonical sample at the 1σ level.

These results indicate that quasars have undergone a significantly greater radio evolution relative to optical evolution, in general agreement with the result in QP1, although the best-fit differential evolution was even more dramatic there, given the form of g(z) used.

The results for the best-fit k_opt and k_rad are identical if the different fiducial radio–optical data set matching radius criteria are used, indicating that the addition of extra radio flux and a few additional sources do not matter for this determination. If we consider that the radio–optical luminosity correlation is inherent in the data, then the orthogonal luminosities are L_opt and L_crr (see Section 3.2) and we can determine the best-fit evolutions g_crr(z) and g_opt(z). These results favor k_opt = 3.5 and k_crr = 2. In this case the best-fit radio evolution can be recovered by Equation (A2) in the Appendix, and it is k_rad = 5.5, in perfect agreement with the results obtained from considering k_opt and k_rad as orthogonal, indicating the robustness of the result.

Next we determine the density evolution ρ(z). One can define the cumulative density function

$$\begin{equation} \sigma (z) = \int _0^z { {{dV} \over {dz}} \, \rho (z) \, dz}, \end{equation} \tag{ 9 }$$

which, following Petrosian (1992) based on Lynden-Bell (1971), can be calculated by

$$\begin{equation} \sigma (z) = \prod _{j}{\left(1 + {1 \over m(j)}\right)}, \end{equation} \tag{ 10 }$$

where j runs over all objects with a redshift lower than or equal to z, and m(j) is the number of objects with a redshift lower than the redshift of object j that are in object j's associated set. In this case, the associated set is again those objects with sufficient optical and radio luminosity that they would be seen if they were at object j's redshift. The use of only the associated set for each object removes the biases introduced by the data truncation. Then the density evolution ρ(z) is

$$\begin{equation} \rho (z) = {d \sigma (z) \over dz} \times {1 \over dV/dz}. \end{equation} \tag{ 11 }$$

However, to determine the density evolution, the previously determined luminosity evolution must be taken out. Thus, the objects' optical and radio luminosities, as well as the optical and radio luminosity limits for inclusion in the associated set for given redshifts in the calculation of σ by Equation (10), are scaled by taking out factors of g_opt(z) and g_rad(z), determined as above. Figures 10 and 11 show the cumulative and differential density evolutions, respectively.⁷ The number density of quasars seems to peak at between redshifts 1 and 1.5, earlier than generally thought for the most luminous quasars (e.g., Shaver et al. 1996), and earlier than that found in R06, but more similar to the peak found for less luminous quasars by Hopkins et al. (2007), and in agreement with Maloney & Petrosian (1999). This is slightly earlier than we found in QP1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The normalization of ρ(z) is determined by Equation (3), with the customary choice of $\int _{L_{\rm min}^{\prime }}^{\infty } {\psi (L^{\prime }) \, dL^{\prime }}=1$. As stated in R06, the main sources of bias in the redshift distribution of the SDSS quasar sample are (1) the differing magnitude limits for z > 2, (2) the effects of emission lines on i-band flux at different redshifts, and, at a somewhat less important level, (3) the inclusion of extended sources at the lowest redshifts (z ∼ 0.3). As discussed in Section 2, we have, following R06, dealt with issue 1 by restricting the sample to a universal i < 19.1 magnitude limit and issue 2 by adopting the R06 K-corrections, which include the effect of emission lines and the continuum spectrum. Since we do not see any hint of a peak or deviation in ρ(z) or σ(z) at z ∼ 0.3, we do not deem issue 3 to significantly affect our conclusions. The density evolution ρ(z) can be fit with a broken polynomial in z, with coefficients shown in Table 1.

Knowing both the luminosity evolutions g_a(z) and the density evolution ρ(z), one can form the luminosity density functions ${\it \unicode{xa3}}_{\rm a}(z)$, which are the total rate of production of energy of quasars as a function of redshift, by

$$\begin{equation} {\it \unicode{xa3}}_{\rm a}(z) = \int dL_{\rm a} \, {L_{\rm a} \, \rho (z) \, {dV \over dz}} \, dz \, = \, \langle L^{\prime }_{\rm a} \rangle \, g_a(z) \, \rho (z) \, {dV \over dz}. \end{equation} \tag{ 12 }$$

We show both the 2500 Å  and optical and 1.4 GHz radio luminosity density functions in Figure 12. As expected, because of the greater luminosity evolution in the radio band, the radio luminosity density peaks earlier in redshift than the optical luminosity density.

Zoom In Zoom Out Reset image size

6.1. General Considerations

In a parallel procedure we can use the local (redshift evolution taken out, or "de-evolved") luminosity distributions (and de-evolved luminosity thresholds) to determine the "local" LFs ψ_a(L'_a), where again the subscript a denotes the waveband and the prime indicates that the luminosity evolution has been taken out. We first obtain the cumulative LF

$$\begin{equation} \Phi _a(L_{a}^{\prime }) = \int _{L_{a}^{\prime }}^{\infty } {\psi _a(L_{a}^{\prime \prime }) \, dL_{a}^{\prime \prime }}. \end{equation} \tag{ 13 }$$

Following Petrosian (1992), Φ_a(L'_a) can then be calculated by

$$\begin{equation} \Phi _a(L_{a}^{\prime }) = \prod _{k}{\left(1 + {1 \over n(k)}\right)}, \end{equation} \tag{ 14 }$$

where k runs over all objects with a luminosity greater than or equal to L_a, and n(k) is the number of objects with a luminosity higher than the luminosity of object k that are in object k's associated set, determined in the same manner as above. The luminosity function ψ_a(L'_a) is

$$\begin{equation} \psi _a(L_{a}^{\prime }) = - {d \Phi _a(L_{a}^{\prime }) \over dL_{a}^{\prime }}. \end{equation} \tag{ 15 }$$

In Section 4 we determined the luminosity evolution for the optical luminosity L_opt and the radio luminosity L_rad. We can form the local optical ψ_opt(L'_opt) and radio ψ_rad(L'_rad) LFs straightforwardly, by taking the evolutions out. As before, the objects' luminosities, as well as the luminosity limits for inclusion in the associated set for given redshifts, are scaled by taking out factors of g_rad(z) and g_opt(z), with k_rad and k_opt determined in Section 4. We use the notation L → L' ≡ L/g(z). For the local LFs, we use the customary normalization $\int _{L_{\rm min}^{\prime }}^{\infty } {\psi (L^{\prime }) \, dL^{\prime }}=1$. This normalization may be biased by around 8% due to quasar variability as discussed in Section 8.

6.2. Local Optical Luminosity Function

Figures 13 and 14 show the local cumulative Φ_opt(L'_opt) and differential ψ_opt(L'_opt) local optical LFs of the quasars in the canonical combined data set and that of the optical-only data set.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The optical LF shows possible evidence of a break at ∼10³⁰ erg s⁻¹ Hz⁻¹, which was present already in data analyzed in Petrosian (1973) and was also seen in QP1. With the optical-only data set, fitting a broken power law yields values −2.8  ±  0.2 and −4.1 ± 0.4 below and above the break, respectively. With the canonical combined radio–optical data set, fitting a broken power law yields values −2.8 ± 0.3 and −3.8 ± 0.5. If we allow for the possibility of additional uncertainty resulting from the consideration of possible radio incompleteness at faint fluxes (see discussion in Section 8), the range on the power law below the break changes to −2.1 ± 0.1, and above the break the results are not affected.

As the optical LF has been studied extensively in various AGN surveys, we can compare the slope of ψ_opt(L'_opt) obtained here to values reported in the literature. In QP1, we found power-law values of −2.0 ± 0.2 and −3.2 ± 0.2 below and above the break, respectively. Boyle et al. (2000), using the 2dF optical data set (but with no radio overlap criteria), use a customary broken power-law form for the LF, with values ranging from −1.39 to −3.95 for different realizations, showing reasonable agreement.⁸

6.3. Local Radio Luminosity Function

Figure 15 shows the local radio LF ψ_rad(L'_rad). It is seen that the local radio LF contains a possible break around 2 × 10³¹ erg s⁻¹ Hz⁻¹, with a power-law slope of −1.5 ± 0.1 below the break and −2.6 ± 0.1 above the break. The range for the power law above the break is increased slightly to −2.5  ±  0.3 if the effects of possible radio incompleteness are included, as in Section 8, while below the break it is unchanged, indicating that this is a small effect. It is interesting to note that the value of the radio-loudness parameter R corresponding to the break in the radio and optical LFs is very close to the critical value R ∼ 10 widely discussed in the context of the RL/RQ dichotomy.

Zoom In Zoom Out Reset image size

We can also construct the local radio LF from ψ_opt(L'_opt) and ψ_crr(L'_crr) with Equation (A1) of the Appendix, under the assumption that the radio and optical luminosity correlation is intrinsic. This is shown by the small black points in Figure 15. It is seen that ψ_rad(L'_rad) constructed in this way is nearly identical to that determined by considering the radio and optical luminosity correlation to be induced, indicating that the radio luminosity function result is robust to this consideration.

In QP1 for the local radio luminosity function we found a break at 10³¹ erg s⁻¹ Hz⁻¹, with a power-law slope of −1.7 ± 0.1 below the break and −2.4 ± 0.1 above it, in excellent agreement with the results here. The slope above the break seen here is similar to earlier results of Schmidt (1972) and Petrosian (1973), which probed only high luminosities. A more complete comparison can be made with Mauch & Sadler (2007), who form radio LFs of local sources in the Second Incremental Data Release of the 6 degree Field Galaxy Survey (6dFGS) radio catalog. For the sources they identify as AGNs, they find a break at 3.1 × 10³¹ erg s⁻¹ Hz⁻¹, with slopes of −2.27 ± 0.18 and −1.49 ± 0.04 above and below the break, respectively, in good agreement with our results.

In Figure 15 we also show the radio LFs for the same range of L'_rad, as determined by Kimball et al. (2011) for two samples; one for a sample of SDSS optical and NRAO Very Large Array Sky Survey (NVSS) radio data, and another for deep radio observations with the Extended Very Large Array in which almost every SDSS quasar in a narrow redshift range in a field was detected in radio. The results presented there are the luminosity function only for sources in the redshift range 0.2 < z < 0.3, which we have scaled by means of Equation (5) and converted from 6 GHz to 1.4 GHz with a radio spectral index of 0.6 to obtain the local luminosity function at 1.4 GHz for direct comparison with our results. It is seen that our LFs agree on the faint end of the range considered here but potentially disagree only on the very bright end.

As stated in the introduction, naively one may expect that because the ratio R is independent of cosmological model and nearly independent of redshift, the raw observed distribution of R would provide a good representation of the true distribution of this ratio. In Appendix A of QP1 we show that this assumption would likely be wrong. In Figure 16 we show this raw distribution of R by the triangles, arrived at by using the raw values of R from the data and forming a distribution in the manner of Equations (14) and (15) with no data truncations accounted for. As discussed in Section 3, we can reconstruct the local distribution of $G_{R^{\prime }}(R^{\prime })$, as in Equation (1), which provides for a more proper accounting of the biases and truncations and is an estimate of the true, intrinsic distribution. The results of this calculation are also shown in Figure 16 by the filled circles. We include a range of error for the extremal possibilities that the correlation between the radio and optical luminosities is entirely induced or entirely intrinsic (open circles). In the latter case, $G_{R^{\prime }}(R^{\prime })$ is given by Equation (A3) of the Appendix.

Zoom In Zoom Out Reset image size

The distribution $G_{R^{\prime }}(R^{\prime })$ calculated from the data taking into account the evolutions and truncations is clearly different from the raw distribution and shows no evidence of bimodality in the range of R considered. The raw distribution is seen to be weighted toward much higher values of R. Using the naive, raw distribution instead of the reconstructed intrinsic one would result in estimates of quasar properties that would be significantly and systematically biased.

We also know the best-fit redshift evolution of the ratio, given Equation (A4). The change in the distribution of R with increasing redshift is also shown in Figure 16.⁹ Another way to look at this is that we have found that the radio luminosity evolves at a different rate than the optical luminosity, with the consequence that their ratio is a function of redshift. The radio loudness of the population increases by a factor of 3 by redshift 1, a factor of 8 by redshift 3, and a factor of 11 by redshift 5. This is a less dramatic evolution of the ratio than we found in QP1, where we used a simpler parameterization (see Section 3.1), although still quite significant. The general trend toward increasing R with redshift can be validated in a simple way by examining the median value of R versus redshift in bins of redshift, which indeed shows a steady increase.

This differential evolution is in disagreement with the result presented by J07, who show a decrease in fraction of RL sources with increasing redshift, which could be the case if the radio luminosity were to evolve more slowly than the optical luminosity. They, however, do not determine individual evolutions or LFs. On the other hand, Miller et al. (1990) have noted that the fraction of RL quasars may increase with redshift, which they attribute to a difference in the evolutions of the two populations (RL and RQ). Donoso et al. (2009) compute radio and optical LFs at different redshifts and reach the same conclusion. Cirasuolo et al. (2006) also find that the radio-loud fraction may modestly increase at high redshift. Although not directly comparable, LaFranca et al. (2010) show a similar evolution for R_x, the ratio of radio to X-ray luminosity, as we show here for R. As discussed in Section 9.1, we believe that our differences with J07 stem from our inclusion of only radio-detected quasars in the joint radio–optical sample.

We note that our results favor one population in the range of R considered here, in the sense that the distribution of G(R), recovered from considering the data truncations inherent in the survey and correlations between the luminosities, is continuous and without any bimodality. This is in agreement with the conclusion we reached in QP1. One may hypothesize about a bimodality that is not centered on the commonly assumed value of R = 10 or above that value, but at a lower value of R below those values probed in this analysis. Our results disfavor that as well. As seen in Figures 14 and 15, the lack of a significant change in the shapes of the computed local optical or radio LFs when the radio flux limits for inclusion in the analysis are changed, or the parent optical-only data set is considered, argues against the presence of an additional population of quasars that would become more prominent below the lowest values of R considered here. The results of Kimball et al. (2011), which probe even lower radio fluxes and extend the radio luminosity function to lower values, do not show a large enough population of low-R quasars that would produce a bimodality in the distribution of R.

In Figure 17 we plot our results along with two determinations of the raw observed R distribution available in the literature, that of Ivezic et al. (2004) and that of Cirasuolo et al. (2006). These determinations from the literature are scaled to the ratio of 5 GHz radio to 2500 Å  optical luminosity assuming optical and radio spectral indices of 0.5 and 0.6, respectively. In the former case (red and orange squares in the figure), the data set is from an SDSS × FIRST sample that is flux limited in radio and optical as in this work. The raw R distribution seen there is similar to the one seen here, although with a more pronounced dip at R ∼ 10. In the latter case (blue stars in the figure), the data are from a small patch of sky with deep VLA observations where the sample of SDSS quasars has 100% completeness down to S_1.4 GHz = 60 μJy, and with such radio completeness we would expect the observed distribution to be consistent with our reconstructed intrinsic distribution at low values of R, which indeed it is, but not necessarily at high values of R, where there may be sources missing, either due to cosmic variance or because they have significant radio flux but insufficient optical flux to be included in the SDSS sample (see the discussion in Section 9).

Zoom In Zoom Out Reset image size

Luminosity-dependent density evolution. One may be concerned that luminosity-dependent density evolution (LDDE), which is not considered in the functional forms for the LF used here, may be necessary to represent the evolution of the LF. As a test of this scenario, we divide the data into high and low halves of de-evolved luminosity L' (cutting on optical luminosity) and check the similarity of the computed density evolutions for the two sets. The density evolutions computed for both cuts are similar to each other, with the high-redshift half peaking in ρ(z) at z = 1.3 and the low-redshift half peaking in ρ(z) at z = 1.1. Given the similarity of these distributions to each other and to that computed from the data set as a whole, we conclude that we are justified in neglecting LDDE for the purposes of this analysis.

Optical measurement errors. There is the possibility that errors in the optical magnitudes could lead to a bias that could affect the results (Eddington 1940). The bias introduced by these errors would be negligible if the source counts as a function of flux, i.e., the log N–log S distribution, were flat (i.e., N(> S)∝S⁰ and dN/dS∝S⁻¹). Since the number density of sources increases with decreasing flux, it is more likely that a source will be included than excluded. However, the magnitude of this effect will depend on the faint-end source counts slope, and the shallower the slope, the smaller the effect. For constant fractional flux measurement errors, an error will be introduced on the normalization of the source counts, and therefore the LFs, and can be approximated by 1/2 δ² m_below (m_below + 1) (Teerikorpi 2004), where δ is the fractional error in flux and m_below is the faint-end differential source counts power-law slope. To be conservative, we will adopt errors of 0.2 mag in i-band, to account for the intrinsic rms scatter due to source variability, although the typical reported SDSS DR7Q measurement errors are lower, on the order of a few hundredths of a magnitude. For the faint-end magnitudes of 19.1 and this error, δ is around 0.16. For m_below ∼ 2 the bias on the luminosity function normalization will be around 7.6%. On the other hand, there will be an effect on the reconstructed slope of the luminosity function only if the fractional measurement errors change with luminosity. The data show only a modest dependence on i-band reported error with magnitude at magnitudes below 19.1, with magnitude errors at most a factor of 1.5 higher at the highest magnitude end of that range than for the lowest magnitudes, which corresponds to a small fractional flux error. Because data at a wide range of luminosities correspond to a given flux, we consider this effect to be negligible, even if scatter due to variability has some flux dependence.

Radio incompleteness. Lastly, the selection function for the FIRST objects used here might not be a sharp Heaviside function at a peak pixel flux density of 1 mJy, but rather smeared out. According to Figure 1 of J07, the selection function of FIRST for SDSS optically identified quasars in SDSS data release 3 is such that at an integrated flux density of 1 mJy only about 55% of sources are seen, and this number rises to 75% at 1.5 mJy and about 85% at 2 mJy.¹⁰ We have considered the sample to be limited by the peak pixel flux (i.e., surface brightness limited) in the radio rather than being limited by the integrated flux, in accordance with the criteria set forward in Becker et al. (1995).

The way to test the effects on our analysis of possible radio incompleteness is to repeat the analysis limiting the sample to a higher radio flux, where the sample would presumably be more complete, and determine the extent to which the calculated parameters change in a systematic way. We have done so with a lower radio flux limit of 2 mJy (4251 objects in the canonical sample), as opposed to the original 1 mJy (5445 objects in the canonical sample). The effect, propagated through the analysis, is quite minor. The 1σ uncertainties on k_opt and k_rad are extended slightly. There is no discernible effect on the correlation parameter α. There is a negligible effect on the density evolutions, a small effect on the local LFs that are reported in Section 6, and a small effect on the shape of the G_R(R) distribution that is sub-dominant to the effect of considering the radio–optical luminosity correlation to be intrinsic or induced. To the extent that faint flux radio incompleteness is present in the sample considered here, it does not seem to have a large systematic effect on the determination of the parameters in, and conclusions of, this analysis.

We have used a general and robust method to determine the radio and optical luminosity evolutions simultaneously for the quasars in the SDSS DR7 QSO × FIRST data set, which combines 1.4 GHz radio and i-band optical data for thousands of quasars ranging in redshift from 0.64 to 4.82 and over five orders of magnitude in radio loudness.

As can be seen, the results for the SDSS DR7 QSO × FIRST data set employed here are similar in bulk to the results for the White et al. data set that we presented in QP1. In this work we show that, importantly, the major conclusions are not sensitive to whether the correlation between the radio and optical luminosities is considered to be induced by similar redshift evolution or intrinsic. In the previous work, we assumed the correlation to be intrinsic. Additionally, the luminosity evolution found for the optical-only data set matches that found for the combined radio–optical data set, providing a test that the technique has properly handled the truncations and correlations inherent in the data.

We have seen that the results of this analysis are nearly identical for the two different radio–optical matching radius criteria considered: (1) that adopted by J07 with a 5'' radius for single matches and a 30'' radius for multiple matches, and (2) a universal 5'' matching radius, which we feel is more appropriate. The reason for this preference is that 30'' radius corresponds to a physical scale >100 kpc projected, at any redshift z > 0.2. Extended lobes in radio galaxies are known to reach similar or even larger sizes. These are, however, objects viewed at large inclinations, unlike quasars observed, by definition, at much smaller viewing angles (Barthel 1989). With such small viewing angles, projection effects are expected to limit the projected sizes of the extended structures in radio-loud quasars. In addition, evolutionary effects may play a role as well, with distant quasars being typically younger, and therefore more compact in radio, than nearby evolved radio galaxies. This is supported by a detailed investigation of radio morphologies of SDSS (DR3) × FIRST quasars by De Vries et al. (2006), who noted that an overwhelming majority of the selected sources are characterized by compact radio morphologies. For all these reasons, we believe that applying a universal 5'' matching radius criterion in assigning FIRST counterparts to the SDSS bona fide quasar sources is more appropriate, while the 30'' criterion for multiple matches has a high probability of including radio flux from objects unrelated to the optical quasar that they are supposedly linked to.

9.1. Differential Evolutions toward Radio Loudness with Increasing Redshift

9.2. Lack of Bimodality in Radio Loudness

Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the U.S. Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS Web site is http://www.sdss.org/.

As discussed in Section 3.2, in the case of intrinsically correlated luminosities, the LFs in different wavebands are not a priori separable, and one must do a coordinate rotation of the form in Equation (6) to a correlation-reduced luminosity to achieve a separable luminosity function. In that case, one must determine the two now independent luminosity–redshift correlation functions g_opt(z) and g_crr(z) that describe the luminosity evolutions. The full procedure is detailed in Section 4, and this result can be compared to the evolutions determined from assuming no intrinsic (i.e., only redshift induced) correlation between the radio and optical luminosities, to determine whether the effect is significant.

The local LFs of uncorrelated luminosities L'_opt and L'_crr can then be used to recover the local radio LF by a straightforward integration over L'_crr and the true local optical LF as

$$\begin{eqnarray} &&\psi _{\rm rad}(L_{\rm rad}^{\prime }) = \nonumber \\ &&\quad \int _0^{\infty } { \psi _{\rm opt}(L_{\rm opt}^{\prime }) \, \psi _{\rm crr}\left({{ L_{\rm rad}^{\prime } } \over ({{L_{\rm opt}^{\prime }/L_{\rm fid}})^{\alpha } } } \right) \, {{dL_{\rm opt}^{\prime }} \over ({{L_{\rm opt}^{\prime }/L_{\rm fid}})^{\alpha }} } \,}.\qquad \end{eqnarray} \tag{ A1 }$$

As stated above, this procedure can be used for the determination of the radio LF at any redshift, from which one can deduce that the radio luminosities also undergo luminosity evolution with

$$\begin{equation} g_{\rm rad}(z) = g_{\rm crr}(z) \, \times \, [g_{\rm opt}(z)]^{\alpha } \end{equation} \tag{ A2 }$$

(cf. Equation (6)).

Similarly, we can determine the local distribution of the radio-to-optical luminosity ratio, R' = L'_rad/L'_opt = L'_crr × L'_opt^{α − 1} × L_fid^−α, as

$$\begin{equation} G_{R^{\prime }} = \int _0^{\infty } { \psi _{\rm opt}(L_{\rm opt}^{\prime }) \, \psi _{\rm crr}\left({{ R^{\prime } \, L_{\rm fid}} \over ({{L_{\rm opt}^{\prime }/L_{\rm fid}})^{\alpha -1} } } \right) \, {{dL_{\rm opt}^{\prime }} \over {{L_{\rm opt}^{\prime }}^{\alpha - 1} \, L_{\rm fid}} } \,} \end{equation} \tag{ A3 }$$

and its evolution

$$\begin{equation} g_{R}(z) = g_{\rm crr}(z) \, \times \, [g_{\rm opt}(z)]^{\alpha -1} = {g_{\rm rad} \over g_{\rm opt}}. \end{equation} \tag{ A4 }$$

These can then be compared with ψ_rad(L'_rad) and $G_{R^{\prime }}$ as determined from assuming no intrinsic correlation between the radio and optical luminosities, which is done in Figures 15 and 16. We see that the main conclusions of this work are not dependent on whether the observed correlation between L_rad and L_opt is intrinsic or induced.