format: text

text:
# Acknowledgments

## Acknowledgments

A. C. P. is supported by Alexander von Humboldt Foundation and G. S. by the US Department of Energy under grant DE-FG05-91ER40627. A. C. P. would like to thank the University of Tennessee, Knoxville for its kind hospitality.


text:
# High-Redshift Quasars Found in Sloan Digital Sky Survey Commissioning Data IV: Luminosity Function from the Fall Equatorial Stripe Sample1

## 1 Introduction

Soon after quasars were first discovered, Schmidt (1968) found them to evolve strongly with redshift. In the following 25 years, a consensus has been reached that the number density of optical quasars peaks at $`z2.53`$, and declines towards both lower and higher redshifts. The luminosity function of optical quasars at $`z<2.5`$ has been well studied through UV-excess and slitless spectroscopic surveys (e.g., Boyle, Shanks & Peterson (1988), Hewett et al. (1991), Pei (1995), Boyle et al. (2000)). At the high-redshift end, both multicolor surveys (Warren, Hewett & Osmer (1994), WHO hereafter; Kennefick, Djorgovski & de Carvalho (1995), DPOSS hereafter) and grism surveys (Schmidt, Schneider & Gunn (1995), SSG here after) have shown that the number density of quasars declines rapidly from $`z3`$ to $`4.3`$. The sample sizes for these surveys are relatively small, so the rate of this decline and the shape of the high-redshift quasar luminosity function are not nearly as well-constrained as at low redshift. The quasar luminosity function provides important constraints for models of quasar evolution (e.g., Efstathiou & Rees (1988), Turner (1991), Haiman & Loeb (1998), Haiman & Hui (2000)), and the nature of the high-redshift UV ionizing field (Madau, Haardt & Rees (1999)).

In Paper III (Fan et al. 2000a ), we presented a sample of 39 quasars at $`z>3.6`$ and $`i^{}20`$, selected from 182 deg2 of Fall Equatorial Stripe multicolor imaging data taken during the commissioning phase of the Sloan Digital Sky Survey (SDSS, York et al. (2000), see also Fukugita et al. (1996), Gunn et al. (1998)). The selection completeness function of this multicolor sample is calculated based on model quasar colors. In this paper, we derive the high-redshift quasar luminosity function based on the quasar sample and selection function presented in Paper III. In §2, we use three different statistical approaches to find the luminosity function: a $`1/V_a`$ estimator; a maximum likelihood solution; and Lynden-Bell’s non-parametric $`C^{}`$ estimator (Lynden-Bell (1971)). In §3, we discuss the rate of decline of the quasar spatial density towards high redshift and the shape of the quasar luminosity function based on the results of §2, and compare these new results with those of previous studies. In §4, we discuss possible selection biases on the evolution of the luminosity function and theoretical implications of our results.

## 2 Derivation of Luminosity Function

The differential quasar luminosity function, $`\mathrm{\Psi }(L,z)`$, is defined as the number of quasars per unit comoving volume, per unit luminosity as a function of luminosity and redshift. It can also be written as $`\mathrm{\Psi }(M,z)`$, the number of quasars per unit comoving volume, per unit absolute magnitude:

$$\mathrm{\Psi }(M,z)=\mathrm{\Psi }(L,z)\left|\frac{dL}{dM}\right|=0.92L\mathrm{\Psi }(L,z).$$
(1)
Throughout the paper, we use $`\mathrm{\Psi }(L,z)`$ for the mathematical derivations, while we present our final results in terms of the more familiar $`\mathrm{\Psi }(M,z)`$. The quasar luminosity function can also be presented in the cumulative form: $`\mathrm{\Phi }(L,z)_L^{\mathrm{}}\mathrm{\Psi }(L^{},z)𝑑L^{}`$.

In general, the redshift and luminosity dependences of the quasar luminosity function are not separable. In limiting cases, $`\mathrm{\Psi }(L,z)`$ can be represented either in the form of pure density evolution:

$$\mathrm{\Psi }(L,z)=\psi (L)\rho (z),$$
(2)
where the luminosity function remains the same shape, but with different normalization at different redshift; or in the form of pure luminosity evolution

$$\mathrm{\Psi }(L,z)=\mathrm{\Psi }[L/g(z)]/g(z),$$
(3)
where $`g(z)`$, the characteristic luminosity, is a function of redshift while the normalization of the luminosity function is a constant. It might be better presented by a combination of the two cases above (e.g., Koo & Kron (1988)):

$$\mathrm{\Psi }(L,z)=\rho (z)\mathrm{\Psi }[L/g(z)]/g(z).$$
(4)
For $`z3`$, pure luminosity evolution provides a good fit to the quasar luminosity function (e.g. Schmidt & Green (1978), Boyle et al. (2000)). At $`z>3`$, existing quasar samples are still too small to distinguish these models. In general, the evolution can be more complicated. For example, the luminosity function at low and high redshifts may be determined by very different physics and have totally different shapes. Note that if the luminosity function is a power law, $`\mathrm{\Psi }(L,z)L^\beta `$, where $`\beta `$ does not depend on $`z`$, there is no characteristic luminosity, and pure luminosity and pure density evolution cannot be distinguished.

The determination of the high-redshift quasar luminosity function from a flux-limited sample faces several difficulties: (1) the selection function of a multicolor survey is a strong function of luminosity, redshift and the shape of the quasar spectral energy distribution (SED, see Paper III); the selection bias has to be corrected for carefully. (2) The high-redshift quasar samples are typically small (the largest complete sample of quasars at $`z>4`$ prior to this work has 10 objects). The small sample size makes it difficult to parameterize the luminosity function, and the result is also strongly affected by the choice of binning. (3) As the redshift increases, an optical survey samples different parts of the quasar’s intrinsic spectrum. When comparing results from low and high redshift samples, the K-correction is large and quite uncertain, due to the lack of statistics on the continuum shape between the rest-frame UV and optical region of high-redshift quasars.

The luminosity function can be derived using either non-parametric or parametric methods. The most commonly used non-parametric method is the $`1/V_a`$ estimator (e.g. Avni & Bahcall (1980)). We present the $`1/V_a`$ results in §2.1. In §2.2, we generalize the test of correlation between luminosity and redshift distributions described in Efron & Petrosian (1992) and Maloney & Petrosian (1999), by including the selection function. We show that for the limited range of redshift and luminosity our survey covers, the two distributions are not correlated. Thus in §2.3, we assume that the luminosity and redshift dependences are separable, and that the quasar density follows a power law in luminosity and depends on redshift exponentially. The maximum likelihood solutions are derived following the method of Marshall (1985), where each quasar is treated as a $`\delta `$ function in parameter space, and the results are not subject to binning. When the density and luminosity evolution can be separated, one can also use the Lynden-Bell’s $`C^{}`$ estimator (Lynden-Bell (1971)) to calculate the non-parametric marginal distributions of the cumulative luminosity function along the luminosity and redshift directions. In Appendix A, we show how to generalize the $`C^{}`$ estimator to include the survey selection function. The luminosity function estimated using the $`C^{}`$ estimator is presented in §2.4.

The results throughout this paper are presented for two cosmological models: (1) a model with $`\mathrm{\Omega }=1`$ and $`H_0=50`$ km s-1 Mpc-1, which we refer to as the $`\mathrm{\Omega }=1`$ model. Most of the previous studies of the quasar luminosity function have used this cosmology. (2) a $`\mathrm{\Lambda }`$-dominated flat model with $`\mathrm{\Omega }=0.35`$, $`\mathrm{\Lambda }`$=0.65, and $`H_0=65`$ km s-1 Mpc-1 (Ostriker & Steinhardt (1995), Krauss & Turner (1995), Bahcall et al. (1999)), which we refer to as the $`\mathrm{\Lambda }`$-model. The luminosity function is presented in term of $`M_{1450}`$, the absolute AB magnitude of the quasar continuum in the rest-frame at 1450 Å. Assuming a power law $`f_\nu \nu ^\alpha `$, Schmidt et al. (1995) find:

$$M_B=M_{1450}+2.5\alpha \mathrm{log}(4400/1450)+0.12,$$
(5)
where the effective wavelength of the Kron-Cousins B band is 4400 Å, and the factor 0.12 comes from the zero point difference between the AB and Vega-based magnitude systems for quasar-like spectra.

### 2.1 1/$`V_a`$ estimate

We first calculate the space density of quasars in different luminosity and redshift bins using the $`1/V_a`$ method following the discussions in WHO and DPOSS. The available volume $`V_a`$ is defined for a quasar of given luminosity and SED as the comoving volume over which the quasar could have been detected by the survey. In the presence of a selection function, the available volume is weighted by the selection probability. If one is to determine the luminosity in luminosity and redshift bins $`\mathrm{\Delta }z,\mathrm{\Delta }L`$, the available volume for a quasar with given $`L`$, $`z`$ and SED is:

$$V_a=_{\mathrm{\Delta }z}p(L,z,\text{SED})\frac{dV}{dz}𝑑z,$$
(6)
where $`p(L,z,\text{SED})`$ is the selection probability of the quasar as a function of redshift, luminosity and the SED. The luminosity function and its statistical uncertainty can be estimated as:

$$\mathrm{\Psi }(L,z)=\underset{i}{}\frac{1}{V_a^i\mathrm{\Delta }L}$$
(7)
$$\sigma (\mathrm{\Psi })=\left[\underset{i}{}\left(\frac{1}{V_a^i\mathrm{\Delta }L}\right)^2\right]^{1/2}$$
(8)
where $`L`$ and $`z`$ are the average luminosity and redshift over the bin, the sum is over all the objects in the bin, and the available volume $`V_a^i`$ is calculated for each object.

In Paper III, two flavors of selection probability are calculated: (1) the probability for a given quasar SED shape, $`p[L,z,\alpha ,\text{EW (Ly}\alpha \text{+NV)}]`$, where $`\alpha `$ is the slope of the power law continuum and EW(Ly$`\alpha `$+NV) is the rest-frame equivalent width of the Ly$`\alpha `$+NV emission line. These two parameters characterize the shape of the SED. The selection probability of each quasar in the sample is listed in Table 4 of Paper III. (2) The average selection probability $`p(L,z)`$ (§5.3 in Paper III), averaged over the intrinsic distribution of $`\alpha `$ and EW. In Paper III, we assume that this distribution is the product of Gaussians in $`\alpha `$ and $`\mathrm{EW}(\mathrm{Ly}\alpha +\mathrm{NV})`$, and find that their mean and standard deviations are: $`\alpha =0.79\pm 0.34`$ and $`\mathrm{EW}(\mathrm{Ly}\alpha +\mathrm{NV})=69.0\pm 18.3`$ Å. In Figure 1, we present the $`1/V_a`$ results of the quasar luminosity function using both flavors of selection probability under the $`\mathrm{\Omega }=1`$ and the $`\mathrm{\Lambda }`$-model cosmologies. Figure 1(a) and (b) are for the $`\mathrm{\Omega }=1`$ model. In Figure 1(a), the density is calculated using the average selection function $`p(L,z)`$, while in Figure 1(b), the density is calculated using the selection probability of each quasar for its SED type, $`p[L,z,\alpha ,\text{EW(Ly}\alpha \text{+NV)}]`$. Similarly, Figure 1(c) and (d) are for the $`\mathrm{\Lambda }`$-model. In each case, we divide the sample of 39 quasars in 3 redshift bins: $`3.6<z<3.9`$ (18 quasars), $`3.9<z<4.4`$ (14 quasars) and $`4.4<z<5.0`$ (7 quasars). The redshift bins are further divided into two or three luminosity bins. There are 3 – 6 quasars in each bin.

The different methods of treating the selection function give consistent results. This can be understood as follows. For simplicity, we ignore the dependence on EW, $`z`$ and $`L`$, and assume that the quasar distribution only depends on $`\alpha `$. For an intrinsic distribution function $`f(\alpha `$) and a selection function $`p(\alpha `$), the observed distribution will be $`n(\alpha )=p(\alpha )f(\alpha )`$. In the discrete case, $`n(\alpha )=_i\delta (\alpha \alpha _i)`$. The effective density of quasars can be estimated as:

$`\rho _{\mathrm{eff}}`$ $`=`$ $`{\displaystyle \frac{1}{V_a}}{\displaystyle f(\alpha )𝑑\alpha }`$ (9)
$`=`$ $`{\displaystyle \frac{1}{V_a}}{\displaystyle \frac{p(\alpha )f(\alpha )d\alpha }{p(\alpha )}}`$
$`=`$ $`{\displaystyle \frac{1}{V_a}}{\displaystyle \frac{n(\alpha )d\alpha }{p(\alpha )}}`$
$`=`$ $`{\displaystyle \frac{1}{V_a}}{\displaystyle \underset{i}{}}{\displaystyle \frac{1}{p(\alpha _i)}}.`$
Alternatively, if we define the average selection probability:

$$\overline{p}=\frac{p(\alpha )f(\alpha )𝑑\alpha }{f(\alpha )𝑑\alpha },$$
(10)
then the density can also be estimated as:

$$\rho _{\mathrm{eff}}=\frac{p(\alpha )f(\alpha )𝑑\alpha }{\overline{p}V_a}=\frac{N}{\overline{p}V_a},$$
(11)
where $`N`$ is the total number of objects observed. Therefore, if, as we have assumed, the SED shape is not a function of redshift or luminosity, the two approaches only differ in the different weightings applied to the observed $`n(\alpha )`$ distribution; in Eq. (11), each object is given the same weight and normalized by the average selection probability. The two methods should give consistent results if the sample is large enough and all the SED shapes are properly sampled. However, as we showed in Paper III, there are regions in the $`(L,z,\alpha ,\text{EW})`$ parameter space poorly sampled by the selection criteria ($`p1`$). In this case, the average probability $`p(L,z)`$ is much better behaved, and we will use it for most of the analyses which follow.

Due to the small number of quasars in each bin, the appearance of Figure 1 changes somewhat with different binning, but the general trends are obvious: (1) the overall density of quasars declines rapidly with redshift. For the $`\mathrm{\Omega }=1`$ model, the density of quasars at $`M_{1450}26.5`$ decreases from $`2\times 10^8`$ Mpc-3 mag-1 at $`z=3.75`$ to $`6\times 10^9`$ Mpc-3 mag-1 at $`z=4.75`$. (2) the luminosity function rises towards lower luminosity by a factor of $`4`$ per magnitude. These results are very similar for the $`\mathrm{\Lambda }`$-model. At $`z4`$, the behavior of the $`\mathrm{\Omega }=1`$ and $`\mathrm{\Lambda }`$ models are quite similar, as $`\mathrm{\Omega }(z)1`$ when the redshift becomes large.

### 2.2 Test of Correlation between the Luminosity and Redshift Distributions

In this subsection, we test whether the intrinsic luminosity and redshift distributions of the quasar sample are correlated or whether they can be treated separately (as in the case of pure density evolution, Eq. 2). At low redshift, the luminosity function can be fit with a double power law (e.g., Boyle, Shanks & Peterson (1988)). The sample in this paper covers a limited range of absolute magnitude ($`27.5<M_{1450}<25.5`$) and a limited redshift range ($`3.6<z<5.0`$). If, in this case, the luminosity function can be fit with a single power law whose slope is independent of the redshift, the luminosity and redshift distributions of our sample are not expected to be correlated. This correlation test is also crucial for the non-parametric $`C^{}`$ estimator in §2.4, which requires that the data be expressed in terms of two uncorrelated variables.

Correlation tests for truncated data sets (such as flux-limited samples) have been developed by Efron & Petrosian (1992) and Maloney & Petrosian (1999). In our survey, there is essentially no bright-end truncation in the selection (the saturation limit is $`i^{}14`$, four magnitudes brighter than the brightest quasar in the sample). Thus, our data set $`\{L_i,z_i\}`$ is one-side truncated in luminosity. We first consider the case in which the selection function is a step function and the selection boundary is sharp:

$$p(L,z)=\{\begin{array}{cc}1\hfill & \text{if }L>L^{}(z)\hfill \\ 0\hfill & \text{if }L<L^{}(z)\hfill \end{array}$$
(12)
where $`L^{}(z)`$ is the limiting luminosity at a given redshift. In this case, following the discussion in Maloney & Petrosian (1999), one first defines the comparable or associated set for each object $`i`$:

$$J_i=\{j:L_j>L_i,L_j^{}<L_i\}$$
(13)
consisting of all objects which are brighter than the object in question and which would be selected by the survey if they had the same luminosity as the object in question. This is the same set that will be used in the $`C^{}`$ estimator below. It is the luminosity and volume limited subset of the sample that can be constructed for each object. If $`z`$ and $`L`$ are independent, the rank $`R_i`$ of $`z_i`$ in the comparable set, defined as

$$R_i=\text{number of }j\text{, such that }\{jJ_i;z_j<z_i\},$$
(14)
should be distributed uniformly between 0 and $`N_i`$, where $`N_i`$ is the number of points in the comparable set for each object (not including the object in question). The expectation value of $`R_i`$, is $`E_i=(1/2)N_i`$ and the variance $`V_i=(1/12)N_i^2`$. Therefore, one can construct Kendell’s $`\tau `$ statistic:

$$\tau =\frac{(R_iE_i)}{\sqrt{V_i}}.$$
(15)
For $`|\tau |1`$, the luminosity and redshift would not be correlated at $`1\sigma `$ level and could be treated independently.

However, the discussion above only applies to a sample with a sharp boundary. For the quasar survey we discuss in this paper, even though the selection criteria have a sharp cut at the limiting apparent magnitude ($`i^{}<20.05`$ for $`gri`$ selected quasars and $`i^{}<20.2`$ for $`riz`$ selected quasars, see Paper III), the selection function smoothly approaches zero at the limiting luminosity (Figure 8 of Paper III). There is no one-to-one correspondence between the apparent $`i^{}`$ magnitude and the luminosity of the quasar at a given redshift, due to scatter in the quasar SED shape and photometric errors in the SDSS measurements. We thus generalize the correlation test to include the selection function $`p(L,z)`$. First, we define a generalized comparable set:

$$J_i=\{j:L_j>L_i\},$$
(16)
which includes all the objects more luminous than the object in question, since all might have a selection probability between 0 and 1. We then define the total number in the set by weighting each point $`j`$ in the comparable set $`J_i`$. This weight is proportional to $`p(L_i,z_j)`$, the selection probability if the object $`j`$ had the same luminosity $`L_i`$ as the object in question (object $`i`$), and it is proportional to the inverse of its own selection probability $`p(L_j,z_j)`$. It is useful to define a quantity $`T_i`$:

$$T_i=\underset{j=1}{\overset{N_i}{}}\frac{p(L_i,z_j)}{p(L_j,z_j)}.$$
(17)
The definition of the rank changes accordingly:

$$R_i=\underset{j=1}{\overset{N_i}{}}\frac{p(L_i,z_j)}{p(L_j,z_j)},\text{if}z_j<z_i.$$
(18)
The expectation value of the distribution $`R_i`$ is $`E_i=(1/2)T_i`$, and the variance $`V_i=(1/12)(T_i)^2`$. The test statistic $`\tau `$ is defined in the same way as in Eq. (15). Note that these definitions return to its normal form if the selection function has a sharp boundary, since $`p(L_i,z_j)=0`$ when $`L_i<L_j^{}`$. In the case where the selection boundary is smooth, the ranking is weighted by $`p(L_i,z_j)/p(L_j,z_j)`$. The selection probability always increases with larger luminosity and smaller photometric error. Therefore, since $`p(L_i,z_j)p(L_j,z_j)`$ for $`L_i>L_j`$, the weight is always between 0 and 1. In Appendix A, we prove that the comparable set defined in this way gives the correct $`C^{}`$ estimator for the luminosity function in the presence of a selection function.

We calculate the ranking statistic $`\tau `$ of the SDSS Equatorial Stripe quasar sample using the average selection function in Paper III. For the $`\mathrm{\Omega }=1`$ model, $`\tau =0.20`$; for the $`\mathrm{\Lambda }`$ model, $`\tau =0.15`$. Therefore, for the limited redshift and luminosity range this sample covers, the two distributions can be treated independently. For the full quasar sample of the main SDSS survey which will cover a much larger redshift and luminosity range, the possibility of luminosity evolution will need to be introduced. The ranking technique presented here can then be used to determine the functional form of the luminosity evolution ($`g(z)`$ in Eq. 3) of the quasar population (e.g. Maloney & Petrosian 1999).

### 2.3 Maximum Likelihood Fits

In this subsection, we model the luminosity function using the result above that the redshift and luminosity dependences are separable. Imagine dividing up the sample into small bins in $`L`$ and $`z`$, $`L_i`$ and $`z_j`$, with sizes $`\mathrm{\Delta }L`$ and $`\mathrm{\Delta }z`$. The likelihood that $`n_{ij}`$ quasars are found in the ($`L_i`$, $`z_j`$) bin can be written as :

$$=\underset{i,j}{}\frac{e^{\mu _{ij}}\mu _{ij}^{n_{ij}}}{n_{ij}!},$$
(19)
where $`\mu _{ij}`$ is the average number of quasars expected in the bin:

$$\mu _{ij}=_{\mathrm{\Delta }L,\mathrm{\Delta }z}\mathrm{\Psi }(L,z)p(L,z)\frac{dV}{dz}𝑑L𝑑z,$$
(20)
where $`\mathrm{\Psi }(L,z)`$ is the differential luminosity function, and the integral is over the bin. The maximum likelihood solution is obtained by minimizing the function

$$S=2\mathrm{ln},$$
(21)
with respect to the parameters describing the luminosity function.

One could find the maximum likehood solution by fitting the $`1/V_a`$ results in §2.1, but the result would be strongly affected by the choice of binning for a small sample. However, in the limit of an infinitesimal bin size, $`n_{ij}`$ is either 0 and 1, and the result is not dependent on binning. Following Marshall (1985), the likelihood function can then be written as:

$$S=2\underset{i}{\overset{N}{}}\mathrm{ln}[\mathrm{\Psi }(L_i,z_i)p(L_i,z_i)]+2\mathrm{\Psi }(L,z)p(L,z)\frac{dV}{dz}𝑑L𝑑z,$$
(22)
where the sum is over all quasars in the sample. The integral is over the entire $`(L,z)`$ range of the survey, and is equals to the number of quasars expected in the survey for a given luminosity function. This term provides the normalization for the likelihood function.

The likelihood function can be expressed in a slightly different way (see also §5.2 in Paper III):

$$^{}=\underset{i}{}\frac{p_i\mathrm{\Psi }(L_i,z_i)}{\mathrm{\Psi }(L,z)p(M,z)\frac{dV}{dz}𝑑L𝑑z}$$
(23)
$$S^{}=2\mathrm{ln}^{}=2\underset{i}{\overset{N}{}}\mathrm{ln}[\mathrm{\Psi }(L_i,z_i)p(M_i,z_i)]+2N\mathrm{ln}\left(\mathrm{\Psi }(L,z)p(L,z)\frac{dV}{dz}𝑑L𝑑z\right),$$
(24)
where $`N`$ is the sample size. As in Eq. (22), the integral term provides a normalization. These two expressions are closely related. If we define $`N^{}=\mathrm{\Psi }(L,z)p(L,z)\frac{dV}{dz}𝑑L𝑑z`$, then for any parameter $`x`$, it can be shown that

$$\frac{S^{}}{x}=\frac{S}{x}+2\frac{N^{}}{x}\left(\frac{N}{N^{}}1\right).$$
(25)
Thus if we constrain the total number of objects expected based on the maximum likelihood solution to be equal to the total number actually observed, $`N^{}=N`$, Eqs. (22) and (24) are maximized at the same parameters. In practice, they give almost identical results.

The calculation in §2.2 shows that the luminosity and redshift distributions of the quasar sample are not correlated. We assume that the quasar density is a power law function of the quasar luminosity : $`\psi (L)L^\beta `$, or, expressed in terms of absolute magnitude, $`\psi (M)10^{0.4(\beta +1)M}`$. Following SSG, we further assume that the density at a given luminosity declines exponentially with redshift $`\rho (z)10^{Bz}`$. SSG express their results in the form of a cumulative luminosity function

$$\mathrm{log}\mathrm{\Phi }(z,<M_{1450})=AB(z3)+C(M_{1450}+26),$$
(26)
This is equivalent to the differential luminosity function

$$\mathrm{\Psi }(z,M_{1450})=\frac{\mathrm{\Psi }^{}}{10^{0.4[M_{1450}+26\alpha (z3)](\beta +1)}},$$
(27)
where $`\mathrm{\Psi }^{}=10^{0.4\mathrm{ln}10(\beta +1)A}`$, $`\beta =2.5C1`$, and $`\alpha =2.5B/(\beta +1)`$.

As discussed in §5.2 of Paper III, the average probability $`p(M,z)`$ depends on the underlying distribution of quasar SED parameters $`\alpha `$ and EW(Ly$`\alpha `$+NV). The likelihood function (Eqs. 22 and 24) has to be maximized not only with regard to the luminosity function, but also to the intrinsic distribution of $`\alpha `$ and EW(Ly$`\alpha `$+NV) at the same time. In practice, however, the maximum likelihood solutions of the ($`\alpha `$, EW) distribution and of the luminosity function are only weakly correlated. The solutions can be found by first assuming the $`\alpha `$ and EW distributions based on their weighted average (albeit biased by selection effects) to derive the first-order luminosity function and then iterating the results.

Using the likelihood function in Eq. (22), we solve for the parameters of the cumulative luminosity function (Eq. 26). For the $`\mathrm{\Omega }=1`$ model,

$$\mathrm{log}\mathrm{\Phi }(z,<M_{1450})=(7.24\pm 0.19)(0.48\pm 0.15)(z3)+(0.63\pm 0.10)(M_{1450}+26),$$
(28)
where $`\mathrm{\Phi }`$ is in units of Mpc-3. From the maximum likelihood fit, we also find $`\alpha =0.79\pm 0.34`$ and EW(Ly$`\alpha `$+NV) = $`69.3\pm 18.0`$ Å, as reported in Paper III. The best-fit parameters in both cosmological models and in both cumulative and differential forms are given in Table 1, and these models are plotted as dashed lines on Figures 1 and 2. The error bars on the parameters are estimated both by finding the parameters that yield $`S=S_{min}+1`$ (which gives the correct 1$`\sigma `$ errors for a Gaussian distribution), and by using the bootstrap technique. The results from the two methods are consistent with each other within 10%; the values given in Table 1 are from the $`S_{min}+1`$ method. We find that the best-fit luminosity function predicts a total number of $`\mathrm{\Psi }(M,z)p(M,z)\frac{dV}{dz}𝑑M𝑑z=38.7(38.6)`$ quasars observed in the survey area for the $`\mathrm{\Omega }=1`$ ($`\mathrm{\Lambda }`$) model, compared to the 39 quasars actually observed in the survey.

### 2.4 Non-Parametric Determination : Lynden-Bell’s C- Estimator

If the luminosity and redshift distributions of objects in a sample are uncorrelated (§2.2), the bivariate luminosity function $`\mathrm{\Psi }(L,z)`$ can be separated: $`\mathrm{\Psi }(L,z)=\psi (L)\rho (z)`$, where $`\psi (L)`$ and $`\rho (z)`$ are the marginal distributions of the luminosity function in the redshift and luminosity directions. In this case, the cumulative marginal distribution in the luminosity direction, $`\varphi (L)\psi (L)𝑑L`$ can be estimated using Lynden-Bell’s $`C^{}`$ estimator (Lynden-Bell (1971)):

$$\varphi (L_j)=\varphi (L_1)\underset{k=2}{\overset{j}{}}(1+1/N_k),$$
(29)
where $`N_k`$ is the number of objects in the comparable set of the object $`k`$ as defined in Eq. (13), and the objects are sorted according to their luminosities: $`L_1>\mathrm{}>L_{i1}>L_i>\mathrm{}>L_N`$. Similarly, the cumulative marginal distribution in redshift $`\sigma (z)=\rho (z)𝑑z`$ can be written as:

$$\sigma (z_i)=\sigma (z_1)\underset{k=2}{\overset{i}{}}(1+1/M_k),$$
(30)
where $`M_k`$ is the number of objects in the box $`L>L_k`$ and $`z<z_{max}(L_k)`$, and $`z_{max}(L_k)`$ is the maximum redshift the survey can detect at the luminosity $`L_k`$. Therefore, $`M_k`$ is the size of the comparable set in the redshift direction, defined in a similar fashion to $`N_k`$. Here the objects are sorted in the order of their redshifts: $`z_1<\mathrm{}<z_{i1}<z_i<\mathrm{}<z_N`$. If the $`L`$ and $`z`$ distributions are correlated, one should find a set of parameters $`x=x(L,z)`$ and $`y=y(L,z)`$, whose distributions are uncorrelated, and use the $`C^{}`$ estimator to obtain the marginal distributions in $`x`$ and $`y`$ (e.g., Maloney & Petrosian 1999). Note that these equations only give the shape of the distribution, but not the overall normalization of the function. The distributions can be normalized by requiring tha the total predicted number of objects equal that observed:

$$N_{obs}=_0^{\mathrm{}}\psi (L)\sigma [z_{max}(L)]𝑑L=_0^{\mathrm{}}\rho (z)\varphi [L_{min}(z)]𝑑z$$
(31)
As shown in Petrosian (1992), a variety of non-parametric methods can be reduced to Lynden-Bell’s $`C^{}`$ method (Lynden-Bell (1971)) in the limiting case of one object per bin. The $`C^{}`$ estimator makes efficient use of the data by deriving the marginal distribution for the uncorrelated variables, and does not require binning when calculating the cumulative distribution. Obviously, binning is still needed when deriving the differential distribution. But in this case, rather than binning the data in both $`L`$ and $`z`$ as for the $`1/V_a`$ estimator, we only need to bin the data one axis at a time.

For survey with a complicated selection function $`p(L,z)`$, the definition of the comparable set should be changed (Appendix A). For the luminosity direction:

$$N_i=\underset{j}{}\frac{p(L_i,z_j)}{p(L_j,z_j)},\text{where the sum extends over}L_j>L_i,$$
(32)
and for the redshift axis:

$$M_i=\underset{k}{}\frac{p(L_k,z_i)}{p(L_k,z_k)},\text{where the sum extends over}z_k<z_i.$$
(33)
With these definitions, Eqs. (29) and (30) remain valid.

We use the equations above to derive the marginal distributions for the SDSS sample. For the shape of the luminosity function, we calculate the quantity $`\psi (M)/\varphi (M<25.5)`$. the number density of quasars as function of redshift, normalized by the total number at $`M<25.5`$. For the redshift evolution, we calculate $`\rho (z,M<25.5)`$, the total spatial density of quasars with $`M_{1450}<25.5`$ as a function of redshift normalized via Eq. (31). Note that $`\rho (z,M<25.5)=\rho (z)\varphi (M<25.5)`$, therefore the product of the two quantities gives the differential luminosity function: these distributions are shown in Figure 2 for both the $`\mathrm{\Omega }=1`$ model and the $`\mathrm{\Lambda }`$ model, and compared with the maximum likelihood results. The data are binned to have roughly the same number of objects per bin. From Figure 2, we find $`\psi (M_{1450})10^{0.55M_{1450}}`$, and $`\rho (z)10^{0.51z}`$ for the $`\mathrm{\Omega }=1`$ model. These two coefficients are $`0.63\pm 0.10`$ and $`0.48\pm 0.15`$ in the maximum likelihood fits (the dashed lines in Figure 2). Comparing Figures 1 and 2, the $`C^{}`$ results benefit from the smaller bin size and have higher signal-to-noise ratio because they use marginal distributions.

## 3 Evolution of the Quasar Luminosity Function at High Redshift

In §2, we calculated the quasar luminosity function from the SDSS Equatorial Stripe sample using both parametric and non-parametric methods. This quasar sample includes 39 quasars at $`i^{}20`$ and covers 182 deg2. Among them, 18 objects are at $`z>4.0`$ and six at $`z>4.5`$. This is the largest complete quasar sample at $`z>3.6`$ to date. We derive the luminosity function over the range $`3.6<z<5.0`$ and $`27.5<M_{1450}<25.5`$ (for the $`\mathrm{\Omega }=1`$ model). The three methods (the $`1/V_a`$ estimator, the maximum likelihood solution and Lynden-Bell’s $`C^{}`$ estimator) yield consistent results:

1. Over the redshift and luminosity range considered, the distribution of $`z`$ and $`L`$ are uncorrelated, therefore, the luminosity function can be separated: $`\mathrm{\Psi }(L,z)=\psi (L)\rho (z)`$.
2. The quasar number density declines rapidly towards high redshift. It can be fitted by an exponential decline: $`\rho (z)10^{0.5z}e^{1.15z}`$; the spatial density drops by a factor of $`3`$ per unit redshift.
3. The quasar number density rises towards fainter luminosity as $`\psi (M)10^{0.6}M`$, or $`\psi (L)L^{2.5}`$. The density increases by a factor of $`4`$ per magnitude.

In this section, we first compare these results with previous samples at high redshift, and then discuss the shape of the high-redshift luminosity function.

### 3.1 Comparison with Other Surveys at $`z>3.6`$

The quasar luminosity function at $`z>3.6`$ has been calculated for several previous optical quasar surveys (WHO, SSG, DPOSS). The selection of all these surveys is based upon the observations in the rest-frame UV, as in our survey. However, some of the results are expressed in the absolute $`B`$-band magnitude $`M_B`$ in the rest-frame. The UV magnitudes are converted to B-band magnitudes assuming a power law continuum $`f_\nu \nu ^{0.5}`$ (SSG), which we adopt for consistency. From Eq. (5), we have $`M_B=M_{1450}0.48`$. We re-write our maximum likelihood results in the $`\mathrm{\Omega }=1`$ model:

$$\mathrm{log}\mathrm{\Phi }(z,<M_B)=(6.91\pm 0.19)(0.48\pm 0.15)(z3)+(0.63\pm 0.10)(M_B+26)$$
(34)
where $`\mathrm{\Phi }`$ has units of Mpc-3. Since all the surveys to which we will compare are selected from the rest-frame UV fluxes at similar redshifts, uncertainties in the K-correction will not affect the comparison (see also §4.1).

SSG derive the quasar luminosity function in the redshift range $`2.75<z<4.75`$ for 90 quasars at $`r21`$, selected by their Ly$`\alpha `$ emission in the Palomar Transit Grism Survey. The survey covers 61 deg2, and includes 20 quasars at $`z>3.6`$ and 9 quasars at $`z>4.0`$. For the $`\mathrm{\Omega }=1`$ model, SSG find:

$$\mathrm{log}\mathrm{\Phi }(z,<M_B)=6.840.43(z3)+0.75(M_B+26).$$
(35)
In Figure 3, we compare the quasar evolution from the SDSS and from the WHO and SSG surveys. We also show the low-redshift results from the 2dF survey (Boyle et al. (2000)). The new SDSS result is consistent with previous surveys. It shows very similar redshift evolution to that of SSG, and this trend continues towards $`z5`$. The SDSS luminosity function is somewhat flatter than that of SSG, although the difference is only at the 1-$`\sigma `$ level, and is based on rather small sample sizes. Note, however, that the median redshift of the SDSS sample ($`z4`$) is considerably higher than that of SSG’s ($`z3.3`$), which could have an effect on the luminosity function slope (see §3.2).

The consistency of these two surveys is important because they used totally different selection methods and have very different selection functions: SSG is a slitless spectroscopic survey whose selection is based on the detection of emission lines in grism spectra. It is therefore is biased towards quasars with strong emission lines, but is not selected on the continuum shape of the quasars. On the other hand, the SDSS sample selection is based on broad-band colors. It depends on the strength of the emission lines rather weakly, but is biased towards objects with stronger continuum breaks and/or bluer continua, and is not sensitive to red quasars ($`\alpha 1.6`$, see the discussion in Paper III).

The WHO survey (Warren, Hewett & Osmer (1994)) and the DPOSS survey (Kennefick, Djorgovski & de Carvalho (1995)) used broad-band colors to select quasar candidates. These studies are based on photographic photometry from Schmidt plates, which typically have photometric errors $`0.10`$ mag. They use a very similar simulation technique to Paper III to correct for selection effects, although they suffer from a larger selection incompleteness ($`50\%`$). The WHO survey covers an effective area of 43 deg2 down to $`m_{or}=20`$. It contains 86 quasars at $`z>2.2`$, including 8 at $`z>3.5`$ and 2 at $`z>4`$. They found that the quasar number density drops by a factor of $`6`$ from $`z=3.3`$ to $`z=4.0`$, although given the small number of objects in their sample, the rate of evolution is poorly constrained. The DPOSS survey includes 10 quasars at $`z>4`$ and $`16.5<r<19.6`$, covering 681 deg2. They derive the quasar luminosity function at their median redshift $`z=4.35`$. In Figure 4, we compare the cumulative quasar luminosity function from the SDSS sample at $`z4.3`$ with the results from these three previous surveys. The results are all consistent within the error bars. Note that although the two previous multicolor surveys needed a much larger selection correction than the SDSS sample, the final corrected results agree with one another quite well.

The comparisons above show that the decline of the number density of luminous quasars from optical surveys is unlikely to be strongly affected by either the selection technique or the correction for the selection functions. In §4.1, we discuss the effect of K-correction and extinction on the evolution of high-redshift quasar densities.

### 3.2 The Shape of the High-redshift Quasar Luminosity Function

Using the SDSS sample, we find that at $`z4`$ and $`M_B<26`$ , the shape of the quasar luminosity function can be fitted with a power law

$$\psi (L)L^{2.5\pm 0.25},$$
(36)
where the power law index is the average of the results from maximum likelihood and $`C^{}`$ estimators.

At $`z<3`$, the quasar luminosity function is often fitted with a double power law:

$$\mathrm{\Psi }(M_B,z)=\frac{\mathrm{\Psi }^{}(M_B)}{10^{0.4[(\beta _1+1)(M_BM_B^{}(z))]}+10^{0.4[(\beta _2+1)(M_BM_B^{}(z))]}},$$
(37)
In the expression above, pure luminosity evolution is assumed, with a redshift dependence given by the evolution of the characteristic luminosity $`M_B^{}(z)`$. Boyle et al. (2000) find that for the 2dF sample, the normalization of the luminosity function $`\mathrm{\Psi }^{}(M_B)=1.1\times 10^6`$ Mpc-3 mag-1, the faint end slope $`\beta _1=1.58`$, and the bright end slope $`\beta _2=3.43`$ (the $`\mathrm{\Omega }=1`$ model). The error bar on the slopes $`\beta _1`$ and $`\beta _2`$ are $`0.1`$.

If we assume that the luminosity function in Eq. (37) applies to $`z4`$, i.e., the shape of the luminosity function remains the same while only the characteristic luminosity $`M_B^{}(z)`$ varies, we find that in order to match our results at $`M_B26`$, we must have $`M_B^{}24.7`$, more than one magnitude fainter than our survey limit. Therefore, the SDSS results in this paper clearly probe the bright end of the quasar luminosity function.

The results on the bright-end slope of the low-redshift quasar luminosity function from various surveys appear to be quite robust: $`\varphi (L)L^\beta `$, with $`\beta 3.43.5`$ (Boyle et al. 1988, Hewett et al. 1991, Pei 1995, Boyle et al. 2000), much steeper than our high-redshift results ($`\beta 2.5`$). At low redshift, the quasar number density drops by a factor of $`10`$ per magnitude at the bright end, while at $`z4`$, it only drops by a factor of 4.

We can further test the significance of this discrepancy. We derive a new maximum likelihood solution from the SDSS sample by forcing the slope of the luminosity $`\beta =3.43`$ but allowing the normalization $`\mathrm{\Psi }_0`$ and the redshift evolution $`\alpha `$ to vary in the luminosity function (Eq. 27). This is shown as the dashed line in Figure 5. The luminosity function from the SDSS sample using both the $`C^{}`$ estimator and the maximum likelihood estimator from Figure 2 are shown as the points and solid line in the figure. It is evident that this steep slope is not consistent with the data. The statistic $`S`$ (Eq. 22) of the maximum likelihood solution follows a $`\chi ^2`$ distribution. In this case, the likelihood ratio for a specific $`\beta _0`$ can be estimated by calculating $`\mathrm{\Delta }S=S(\beta _0)S_{min}`$, where $`S_{min}`$ is the minimum $`S`$ value when all three parameters are allowed to vary, and $`S(\beta _0)`$ is the $`S`$ value if $`\beta `$ is fixed to be $`\beta _0`$ and the other two parameters are allowed to vary. We find that by forcing $`\beta =3.43`$, $`\mathrm{\Delta }S=8.9`$ with two degrees of freedom. This corresponds to a probability $`p=1\%`$, so the two slopes are inconsistent at $``$ 2.5-$`\sigma `$ level.

We can use another method to estimate the significance of the inconsistency between low and high-redshift slopes: there are six quasars at $`M_{1450}<27`$ in the sample. Using the model with $`\beta =3.43`$, only $`1.7`$ quasars are predicted over the whole survey. The probability of $`6`$ quasars being observed in our survey is only $`0.8\%`$ from Possion statistics.

It is unlikely that this flattening arises because we are missing a large fraction of low luminosity quasars. The photometric errors are quite small and the selection function is very uniform as a function of magnitude (Figure 8 of Paper III) for all objects a few tenths of a magnitude brighter than the survey limit. Figure 5 shows that the slope is much flatter than $`\beta =3.43`$ even over the three most luminous magnitude bins. Thus the bright end slope of the quasar luminosity function is considerably flatter at $`z4`$ than at $`z2`$. Note that SSG find $`\psi (L)L^{2.9}`$, while their median redshift is somewhat lower. The flattening of the bright end luminosity function at high-redshift, if confirmed using a larger sample, indicates that the evolution of quasar luminosity function at $`z>3`$ cannot be pure luminosity evolution, unlike the case for $`z<3`$. Instead, the shape of the quasar luminosity function evolves with redshift as well.

## 4 Discussion

Figure 3 shows that the comoving quasar number density peaks at $`z2.5`$, and drops by a factor of $`20`$ from $`z2.5`$ to 5. This decline is evident in all optical high-redshift quasar surveys to date. However, these results might be biased by three factors: the survey incompleteness; uncertainties in the K-correction when comparing low and high redshift results; and the possibility that a large fraction of quasars are not detectable in optical survey due to dust extinction.

In Paper III, we show that the survey selection criteria are only sensitive to a range of the quasar SED shapes, and that the incompleteness is a function of redshift. When calculating the luminosity function, we first correct this incompleteness by assuming that the distribution of the continuum power-law index $`\alpha `$ is a Gaussian with mean $`0.79`$ and standard deviation $`0.34`$. At $`4.5<z<5.2`$, the survey is sensitive for quasars even with $`\alpha <2.0`$, although no quasars with $`\alpha <1.6`$ are detected in the sample. At $`z<4.5`$, the survey is not sensitive to quasars with $`\alpha 1.6`$. Therefore, the density at $`z4`$ could be underestimated if there was a significant contribution from quasars at $`\alpha 2.0`$ that we failed to detect. However, since we are not biased against these red quasars at the highest redshift bin ($`z>4.5`$), the strong decline in the quasar density seen at $`z5`$ will not be affected by them. Since our selection is not sensitive to very red quasars ($`\alpha <2.5`$), if there were a large population of such objects only at high redshift, the inferred density decline at high redshift from our color-selected sample might be an overestimate. But as we pointed out in §3.1, the same decline is observed in a spectroscopically selected sample (SSG), which argues that neither survey suffers from large incompleteness.

Low-redshift quasar luminosity functions usually are expressed in the rest-frame B band magnitude $`M_B`$. When comparing with high-redshift results, which are measured in the rest-frame UV (e.g. $`M_{1450}`$), a K-correction is needed. It is customary to assume a $`\alpha =0.5`$ power law continuum for the K-correction. An incorrect K-correction might result in a large difference in density. In fact, we find that the UV continua of quasars at $`z4`$ is better fit with $`\alpha 0.8`$, a value that was also found in the similar sized sample of $`z>3`$ quasars (Schneider, Schmidt & Gunn 1991). A change of continuum slope by 0.4 will result in a 0.2 mag difference in the K-correction, or a 30% difference in quasar density. However, when we compare the observed results at $`z4`$ with $`z2`$ the situation is actually better than this calculated difference would imply. The low-redshift quasar surveys (such as 2dF) are selected based on the observed $`B`$-band magnitude ($`\lambda _{eff}4400`$Å). At $`z2`$, the $`B`$-band corresponds to 1470 Å in the rest-frame, almost the same as the case in which we select $`z4`$ quasars based on their $`i^{}`$ band measurements ($`\lambda _{eff}7600`$Å). Therefore, both results are based on the rest-frame UV flux; the comparison is not strongly affected by the K-correction, as long as the two results are transferred to the rest-frame B-band in the same way.

The presence of dust might bias the results on quasar evolution. Dust could arise along the line of sight in the damped Ly$`\alpha `$ systems, or in the quasar environment. Fall & Pei (1993) studied the effect of dust in damped Ly$`\alpha `$ systems on the determination of the quasar luminosity function. They found that 10% - 70% of quasars at $`z3`$ could be missed from optical surveys whose selection was based on B band magnitudes. This fraction cannot be calculated precisely, due to uncertainties in the dust content and dust distribution in damped Ly$`\alpha `$ systems. But the effect of intervening dust must increase quickly with redshift, both because the total number of damped systems increases, and because the extinction curve rises sharply towards the rest-frame UV. The number density of quasars at $`z>4`$ could be significantly affected by dust along the line of sight, but it is not clear whether it is enough to change the general trend of the decline.

High-redshift quasars are known to reside in dusty environments (Omont et al. (1996), Carilli et al. (2000)), and the new millimeter and submillimeter surveys have revealed a population of high-redshift dusty galaxies (e.g. Barger et al. (1998)). The existence of dust extinction changes the bolometric correction from the rest-frame UV luminosity observed in the optical surveys. If the bolometric correction is a strong function of redshift, the evolution of the bolometric luminosity function of quasars could be very different from that of optically selected quasars. However, we’ve seen that the $`\alpha `$ distribution of $`z4`$ quasars is quite similar to that of lower-redshift quasars, implying that the objects in our sample are no more reddened than those at low redshift. Thus in order to invoke dust to explain the decline of optical quasars with redshift, one would require a completely separate, very extincted population of objects whose numbers increase with redshift.

Radio and X-ray selected samples are not strongly affected by dust extinction. Shaver et al. (1996) found that the space density of high-redshift radio quasars also declines rapidly at $`z>3`$. On the other hand, Miyaji, Hasinger & Schmidt (2000) found that there is no evidence for a decrease in the number density of ROSAT-selected X-ray quasars at $`z>3`$, although the statistical significance of this result is only marginal. One difficulty is how to normalize the results from different wavelength samples. At $`z3`$, the spatial density of the quasars in the Miyaji, Hasinger & Schmidt (2000) sample is considerably higher than in the SSG sample, which implies that these X-ray quasars are intrinsically fainter objects than those in SSG. The redshift evolution is expected to depend on luminosity, and there are theoretical reasons to believe that the evolution is flatter at low luminosities (see below).

Studies of the quasar luminosity function put important constraints on models of quasar evolution. Haiman & Loeb (1998) show that the observed quasar luminosity function can be fit with a theoretical model based on the Press-Schechter (1974) approximation, with reasonable assumptions about the relation between quasar luminosity and dark matter halo mass, and the duty cycle of quasar activity. Haiman & Hui (2000) and Martini & Weinberg (2000) further demonstrate that the duty cycle can be better constrained by fitting both the luminosity function and two-point correlation function of high-redshift quasars. Faint quasars reside in less massive dark matter halos in this model. They represent less rare peaks in the density field, and are expected to evolve more slowly than are the bright quasars. The luminosity function in this paper covers only two magnitudes in luminosity, which is not enough to make a detailed comparison with these models. A deep high-redshift quasar survey is needed in order to constrain the faint end slope of the quasar luminosity function and to compare with theoretical models.

Observations of a $`z=5.8`$ quasar (Fan et al. 2000b ) demonstrates that the universe is already highly ionized at $`z6`$, by either the UV ionizing radiation from quasars or from star-forming galaxies. Assuming that the high-redshift quasar luminosity function has the same shape as at low redshift, Madau et al. (1999) estimated that the UV ionizing photons from high-redshift quasars is not sufficient to keep the universe ionized at $`z4`$. In this paper, we have shown that the shape of the high-redshift quasar luminosity function is shallower than at low redshift for bright quasars, implying even fewer ionizing photons from the these quasars than assumed in Madau et al. (1999). However, the faint quasars might contribute more UV ionizing photons than do bright quasars at high redshift. We currently have no knowledge of the faint end slope of the quasar luminosity function at $`z>3`$. Therefore, a deep quasar survey is also needed in order to determine the nature of high-redshift ionizing background.

The SDSS southern survey (see Paper III) will image the Fall Equatorial Stripe (the same area of the sky as the sample used in this paper) 35 – 40 times in the five-year survey period. By co-adding these data, we will be able to select high-redshift quasar candidates down to $`i^{}22.5`$ at $`z4`$, and even fainter for low-redshift quasars. In Paper III, we present 18 quasars down to $`i^{}21`$ selected from areas that have been observed twice during the SDSS commissioning. The high-redshift quasar luminosity function from the southern survey will reach $`M_B24`$, close to the traditional boundary between quasars and Seyfert galaxies ($`M_B=23`$). This sample will give a more complete description of the statistical properties of high-redshift quasars.

The Sloan Digital Sky Survey (SDSS) is a joint project of the University of Chicago, Fermilab, the Institute for Advanced Study, the Japan Participation Group, The Johns Hopkins University, the Max-Planck-Institute for Astronomy, Princeton University, the United States Naval Observatory, and the University of Washington. Apache Point Observatory, site of the SDSS, is operated by the Astrophysical Research Consortium. Funding for the project has been provided by the Alfred P. Sloan Foundation, the SDSS member institutions, the National Aeronautics and Space Administration, the National Science Foundation, the U.S. Department of Energy, and Monbusho, Japan. The SDSS Web site is http://www.sdss.org/. The Hobby-Eberly Telescope (HET) is a joint project of the University of Texas at Austin, the Pennsylvania State University, Stanford University, Ludwig-Maximillians-Universität München, and Georg-August-Universität Göttingen. The HET is named in honor of its principal benefactors, William P. Hobby and Robert E. Eberly. XF and MAS acknowledge support from Research Corporation, NSF grant AST96-16901, the Princeton University Research Board, and a Porter O. Jacobus Fellowship. DPS acknowledges support from NSF grant AST99-00703.

## Appendix A Including the Selection Function in Lynden-Bell’s $`C^{}`$ Estimator

Lynden-Bell’s (1971) $`C^{}`$ estimator calculates the marginal distribution of quasars in luminosity or redshift by constructing a comparable set, the luminosity and volume limited subsample for each object in the sample. The comparable set in luminosity is defined as:

$$J_i=\{j:L_j>L_i,L_j^{}<L_i\},$$
(A1)
where $`L_j^{}`$ is the limiting luminosity at redshift $`z_j`$. For a survey with a complicated selection function such as ours, the limiting luminosity is not well-defined: the selection probability goes to zero smoothly at the survey limit; elsewhere, the selection probability is a complicated function of $`z`$, $`L`$ and the spectral energy distribution of the object. In this case, the comparable set has to be redefined to include the contribution from this selection function. We follow the derivation of Petrosian (1992) for the $`C^{}`$ estimator, now including selection probabilities.

Here we assume that the distribution of $`L`$ and $`z`$ are not correlated, thus the luminosity function can be written as:

$$\mathrm{\Psi }(L,z)=\psi (L)\rho (z).$$
(A2)
The cumulative luminosity function is defined as $`\mathrm{\Phi }(L,z)=_L^{\mathrm{}}\mathrm{\Psi }(L^{},z)𝑑L^{}`$,

$`{\displaystyle \frac{d\mathrm{\Phi }(L,z)}{\mathrm{\Phi }(L,z)}}`$ $`=`$ $`{\displaystyle \frac{\frac{n(L,z)}{p(L,z)}dL}{_L^{\mathrm{}}\frac{n(L^{},z)}{p(L^{},z)}𝑑L^{}}}`$ (A3)
$`=`$ $`{\displaystyle \frac{n(L,z)dL}{_L^{\mathrm{}}\frac{n(L^{},z)}{p(L^{},z)}p(L,z)𝑑L^{}}}`$
where $`n(L,z)`$ is the observed distribution and $`p(L,z)`$ is the selection function.

In order to get the marginal cumulative luminosity function, $`\varphi (L)=_L^{\mathrm{}}\psi (L^{})𝑑L^{}`$, we need to integrate over $`z`$. Since the quantity $`d\varphi (L)/\varphi (L)=d\mathrm{\Phi }(L,z)/\mathrm{\Phi }(L,z)`$ does not depend on $`z`$, it is easy to show that

$`{\displaystyle \frac{d\varphi (L)}{\varphi (L)}}`$ $`=`$ $`{\displaystyle \frac{_0^{\mathrm{}}𝑑zn(L,z)𝑑L}{_0^{\mathrm{}}𝑑z_L^{\mathrm{}}𝑑L^{}\frac{n(L^{},z)}{p(L^{},z)}p(L,z)}}.`$ (A4)
Integrating it gives:

$$\varphi (L)=A\mathrm{exp}\left\{_L^{\mathrm{}}𝑑L\frac{_0^{\mathrm{}}𝑑zn(L,z)𝑑L}{_0^{\mathrm{}}𝑑z_L^{\mathrm{}}𝑑L^{}\frac{n(L^{},z)}{p(L^{},z)}p(L,z)}\right\},$$
(A5)
where $`A`$ is a constant.

In the discrete case, the observed distribution is the sum of a series of $`\delta `$ functions over individual objects:

$$n(L,z)=\underset{i}{}\delta (LL_i)\delta (zz_i).$$
(A6)
Let us assume that the objects in the sample are sorted by their luminosities: $`L_1>L_2>\mathrm{}>L_{i1}>L_i>\mathrm{}>L_N`$. Note that for the discrete case, $`\varphi (L)`$ is the step function: $`\varphi (L_jϵ)=\varphi (L_{j1})`$, and $`\varphi (L_j+ϵ)=\varphi (L_j)`$, when $`ϵ1`$. Thus, one re-writes the integral from $`L_jϵ`$ to $`L_j+ϵ`$, and regards $`p(L,z)`$ as a constant over the integrated range:

$$\varphi (L_j)=\varphi (L_{j1})\mathrm{exp}_{L_jϵ}^{L_j+ϵ}\frac{\delta (LL_j)dL}{_{L_i>L_j}\frac{p(L_j,z_i)}{p(L_i,z_i)}+\mathrm{\Theta }(LL_j)},$$
(A7)
where $`\mathrm{\Theta }(L)`$ is a step function. We can define the comparable set for object $`j`$ as:

$$J_j=\{i:L_i>L_j\},$$
(A8)
and

$$N(L_j)=\underset{iJ_j}{}\frac{p(L_j,z_i)}{p(L_i,z_i)}.$$
(A9)
Note that $`d\mathrm{\Theta }(x)=\delta (x)`$. We have:

$`\varphi (L_j)`$ $`=`$ $`\varphi (L_{j1})\mathrm{exp}\left\{{\displaystyle _{L_jϵ}^{L_j+ϵ}}{\displaystyle \frac{d\mathrm{\Theta }(LL_j)}{N(L_j)+\mathrm{\Theta }(LL_j)}}\right\}`$ (A10)
$`=`$ $`\varphi (L_{j1}){\displaystyle \frac{N(L_j)+1}{N(L_j)}}`$ (A11)
Thus, we have the $`C^{}`$ estimator for the cumulative luminosity function:

$$\varphi (L_j)=\varphi (L_1)\underset{i=2}{\overset{j}{}}(1+1/N_i)$$
(A12)
for the case of an arbitrary selection function.

Similarly, if we define:

$$M_i=\underset{j}{}\frac{p(L_j,z_i)}{p(L_j,z_j)},\text{where the sum extends over}z_j<z_i,$$
(A13)
and the sample is sorted by redshift, $`z_1<\mathrm{}<z_{j1}<z_j<\mathrm{}<z_N`$, the cumulative redshift distribution can be estimated as:

$$\sigma (z_k)=\sigma (z_1)\underset{i=2}{\overset{k}{}}(1+1/M_i).$$
(A14)