CLUSTERING OF OBSCURED AND UNOBSCURED QUASARS IN THE BOÖTES FIELD: PLACING RAPIDLY GROWING BLACK HOLES IN THE COSMIC WEB

The quasars¹⁵ are taken from the sample of luminous mid-IR-selected AGNs presented by H07. Quasars are identified on the basis of their colors in the mid-IR as observed by Spitzer IRAC, using the color–color criterion of Stern et al. (2005, Figure 1(a)), and are selected such that their best estimates of redshift are at z>0.7. To the relatively shallow flux limits of the IRAC Shallow Survey, the AGN sample is highly complete and suffers little contamination from star-forming galaxies (as discussed in detail in Section 7 of H07; see also Assef et al. 2010, 2011).

Zoom In Zoom Out Reset image size

H07 showed that at the ISS flux limits, the IR-selected quasars show a bimodal distribution in optical to mid-IR color. The selection boundary at R − [4.5] = 6.1 can be interpreted as dividing quasars into unobscured (optically bright and so "blue" in R − [4.5]) and obscured (optically faint and so "red" in R − [4.5]) subsets (Figure 1(b)). For the purposes of this study these objects will be referred to as "QSO-1s" and "Obs-QSOs," respectively; the reader is reminded that the selection is based not on optical spectroscopy but only on optical to mid-IR color. This selection yields samples of 839 QSO-1s and 640 Obs-QSOs at z>0.7.

A detailed study of the optical colors, morphologies, and average X-ray spectra of these objects is given in H07. To briefly summarize, H07 found that the QSO-1s have blue optical colors, point-like optical morphologies, and soft X-ray spectra characteristic of unobscured quasars, while the Obs-QSOs had redder optical colors, extended optical morphologies, hard X-ray spectra, and high L_X characteristic of obscured quasars. The sample does not include all obscured quasars, as sources with very large extinction may fall below the IR flux limits of the survey or move out of the Stern et al. (2005) selection region (as shown in Figure 1 of H07; see also Gorjian et al. 2008; Assef et al. 2010). The typical absorbing column for the Obs-QSO sample is estimated to be N_H ∼ 10²²–10²³ cm⁻². We expect the Obs-QSOs to suffer little contamination from bright star-forming galaxies. H07 used an X-ray stacking analysis and constraints from deeper surveys to estimate the possible contamination, and concluded that the contamination is at most ≈30%, and likely significantly smaller (<10%).

For our spatial correlation analysis, we limit the IR-selected quasar sample to the redshift range 0.7 < z < 1.8 to maximize overlap with the normal galaxies in the field (Section 3.2). We also include only objects in regions of good optical photometry and away from bright stars. These criteria yield 563 QSO-1s and 361 Obs-QSOs. Finally, we restrict the QSO-1 sample to those spectroscopically identified as broad-line AGNs to ensure that they unambiguously represent a sample of unobscured quasars and to enable clean tests of photo-z errors (see Section 6.3). Of the full sample of QSO-1s all redshifts, the vast majority (80%) have optical spectra from AGES and 96% of these are classified as broad-line AGNs at 0.7 < z < 4.3, supporting their selection as unobscured quasars. We limit the QSO-1 sample to the 445 that have accurate optical spectroscopic redshifts in the range 0.7 < z < 1.8 and clear broad emission line features. (In a sense this is conservative; we verify that including the 20% of objects with only photo-zs has no significant effect on the clustering results.) Based on these selection criteria, our QSO-1 sample is essentially equivalent to other Type 1 quasar samples selected purely on optical photometric colors and/or spectroscopy (e.g., Richards et al. 2001; Croom et al. 2004; Schneider et al. 2007; Richards et al. 2009b), since the vast majority of spectroscopic Type 1 quasars show AGN-like mid-IR colors (Stern et al. 2005; Richards et al. 2009a). The positions on the sky of the final samples of QSO-1s and Obs-QSOs are shown in Figure 2(a), and their distribution in redshift is given in Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The Obs-QSOs are (by definition) optically faint, and so few (only 7%) are bright enough to obtain good redshifts from MMT optical spectroscopy. AGES targeted objects down to a flux limit of I < 20 for sources that are optically extended, which is the case for almost all the Obs-QSOs. Therefore, the vast majority of the Obs-QSO sample have only photometric estimates of redshift, derived using an artificial neural net technique (Brodwin et al. 2006). Uncertainties on photo-zs using this technique for optically bright quasars are typically σ_z = 0.12(1 + z). However, the errors are more difficult to estimate for optically faint Obs-QSOs, for which there are few spectroscopic redshifts for comparison. Photo-z uncertainties were discussed at length by H07, with the conclusion that typical uncertainties are at most σ_z = 0.25(1 + z) and are likely smaller. Figure 4 shows the photo-zs and spec-zs for the handful of Obs-QSOs with spectroscopic redshifts, as well as those for the QSO-1s for comparison. The impact of photo-z errors in the present clustering analysis is addressed in detail in Section 6.3. As discussed in Section 6.3, random errors in the photo-zs can only tend to decrease the observed clustering amplitude, so we expect the present analysis to provide a robust lower limit on the clustering of the Obs-QSOs.

Zoom In Zoom Out Reset image size

Since the primary aim of this analysis is to compare the clustering of quasars with and without obscuration by dust, it is imperative that the samples are otherwise matched in key properties such as redshift and luminosity. We show in Figure 3 that the redshift distributions of the two samples are similar, and we obtain bolometric luminosities (L_bol) for the quasars by scaling from the rest-frame 8 μm luminosity. We compute the flux at rest-frame 8 μm by extrapolating between the fluxes at 8 and 24 μm in the observed frame, and use this flux to obtain the monochromatic luminosity νL_ν at 8 μm. We then multiply by a luminosity-dependent bolometric correction from Hopkins et al. (2007), which ranges from factors of ≈8 to 11, in order to obtain L_bol. Visual inspection of the Spitzer data shows that essentially all of the quasars have broadly power-law SEDs at these wavelengths, indicating that the rest-frame 8 μm emission is indeed dominated by the AGN. We note that 49 quasars lie outside the region covered by the MIPS 24 μm observations, while 26 (≈3%) of those inside the MIPS area are not detected at 24 μm. For these 75 objects, we use the estimates of L_bol derived from the rest-frame 2 μm luminosity as in Section 4.6 of H07¹⁶.

The distributions in L_bol are almost identical for the QSO-1 and Obs-QSO samples, as shown in the top panel of Figure 1(b). The median and dispersion in log L_bol (erg s⁻¹) are (45.86, 0.37) and (45.83, 0.39) for QSO-1s and Obs-QSOs, respectively, indicating that the two samples are very well matched in bolometric luminosity. For completeness, we note that if we use the L_bol estimates derived from rest-frame 2 μm in H07 and restrict our analysis to QSO-1 and Obs-QSO samples that are matched in L_bol, this has a negligible effect on the clustering results.

The sample of 256,124 galaxies is selected from the deeper SDWFS IRAC observations, with a flux limit [4.5] < 18.6. The galaxies are selected to have best estimates of photometric redshift between 0.5 and 2, with an average photo-z of 〈z〉 = 1.09. The sample includes an optical magnitude cutoff of I < 24 to restrict it to optical fluxes for which the photo-zs are well calibrated. To eliminate powerful AGNs, we have also excluded any object detected in 5 ks Chandra X-ray observations (Kenter et al. 2005) or with 5σ SDWFS detections in all four IRAC bands and colors in the Stern et al. (2005) AGN selection region. The exclusion of AGNs from the galaxy sample removes only 6979 objects and has negligible effect on the results.

The distribution on the sky of the 256,124 galaxies is shown in Figure 2(b), and their distribution in photometric redshift is shown in Figure 3. Photometric redshifts are obtained using an updated version of the Brodwin et al. (2006) algorithm, which is based on template fitting to the optical-IR SEDs. The SED fitting produces a redshift probability density function (PDF) for each object, where P(z) represents the probability that the object lies at redshift z. (Note that the neural net used for the quasar photo-zs does not produce an equivalent estimate of the PDF. Thus, for the quasars we use the best value for the redshift, as discussed in Section 4.) P(z) is normalized such that ∫P(z)dz = 1. For most galaxies, the PDF is roughly Gaussian in shape, although often with a broader tail toward higher redshift. The typical redshift uncertainties are σ_z ∼ 0.1(1 + z), and only a small fraction (≈0.6%) of galaxies show multiple significant peaks in the PDF at different redshifts. Typical galaxy PDFs are shown in Figure 5.

Zoom In Zoom Out Reset image size

In addition to the observed galaxy catalog, the correlation analysis requires a reference sample of objects with random sky positions, in order to compare the observed quasar–galaxy pair counts with the number expected for an uncorrelated distribution. We use a catalog of 8 × 10⁶ random "galaxies" that is assigned to random positions in the regions of good photometry, reflecting the spatial selection function for the SDWFS galaxies.

To measure the spatial clustering of quasars, we can in principle derive the autocorrelation of the quasars themselves, or measure their cross-correlation with a sample of other objects (specifically, normal galaxies) at the same redshifts. Our quasar sample is too small to obtain sufficiently good measurements of their autocorrelation function. However, cross-correlation with galaxies (of which there are ≃300 times as many objects in the Boötes data set) allows far greater statistical power. Furthermore, cross-correlation requires knowledge only of the selection function for the galaxies, which is generally better constrained than that for AGNs. Cross-correlations of AGNs with galaxies have proved an effective technique in a number of previous studies (e.g., Croom et al. 2004; Adelberger & Steidel 2005; Serber et al. 2006; Li et al. 2006; Coil et al. 2007; Wake et al. 2008; Coil et al. 2009; Hickox et al. 2009; Padmanabhan et al. 2009; Mandelbaum et al. 2009; Mountrichas et al. 2009; Donoso et al. 2010; Krumpe et al. 2010).

For the present analysis, the uncertainties in the galaxy photo-zs restrict our ability to perform a full three-dimensional clustering analysis. However, making use of the quasar redshifts and the galaxy photo-z information, we can derive a projected spatial correlation function (w_p(R), with R in comoving h⁻¹ Mpc) that has both higher signal to noise, and a more straightforward physical interpretation than, for example, the purely angular correlation function ω(θ).

The two-point correlation function ξ(r) is defined as the probability above Poisson of finding a galaxy in a volume element dV at a physical separation r from another randomly chosen galaxy, such that

$$\begin{equation} dP=n[1+\xi (r)]dV, \end{equation} \tag{ 1 }$$

where n is the mean space density of the galaxies in the sample. The projected correlation function w_p(R) is defined as the integral of ξ(r) along the line of sight,

$$\begin{equation} w_p(R)=2\int _{0}^{\pi _{\rm max}}\xi (R,\pi)d\pi, \end{equation} \tag{ 2 }$$

where R and π are the projected comoving separations between galaxies in the directions perpendicular and parallel, respectively, to the mean line of sight from the observer to the two galaxies. By integrating along the line of sight, we eliminate redshift-space distortions owing to the peculiar motions of galaxies, which distort the line-of-sight distances measured from redshifts. w_p(R) has been used to measure correlations in a number of surveys, for example, Sloan Digital Sky Survey (SDSS; Zehavi et al. 2005; Li et al. 2006; Myers et al. 2009; Krumpe et al. 2010), 2SLAQ (Wake et al. 2008), DEEP2 (Coil et al. 2007, 2008, 2009), Boötes (Hickox et al. 2009; Starikova et al. 2010), COSMOS (Gilli et al. 2009), and GOODS (Gilli et al. 2007b).

In the range of separations 0.3 ≲ r ≲ 50 h⁻¹ Mpc, ξ(r) for galaxies and quasars is roughly observed to be a power law

$$\begin{equation} \xi (r)=(r/r_0)^{-\gamma }. \end{equation} \tag{ 3 }$$

For sufficiently large π_max such that we average over all line-of-sight peculiar velocities, w_p(R) can be directly related to ξ(r) (for a power-law parameterization) by

$$\begin{equation} w_p(R)=R\left(\frac{r_0}{R}\right)^\gamma \frac{\Gamma (1/2)\Gamma [(\gamma -1)/2]}{\Gamma (\gamma /2)}. \end{equation} \tag{ 4 }$$

We use Equation (4) to obtain power-law parameters for the observed correlation functions, to facilitate straightforward comparisons to other works. However, we note that a number of recent studies have shown evidence for separate terms in the correlation function owing to pairs of galaxies found within a single DM halo (the "one-halo" term), and from pairs in which each galaxy is in a different halo (the "two-halo" term; e.g., Zehavi et al. 2004; Zheng et al. 2007; Coil et al. 2008; Brown et al. 2008; Zheng et al. 2009). A halo occupation distribution (HOD) analysis accounting for both the one- and two-halo terms can provide valuable constraints on the distribution of objects within their DM halos; however, a full HOD calculation is beyond the scope of the present analysis.

To measure w_p(R) for the quasar–galaxy cross-correlation, we employ the method developed by M09. This technique makes use of the full photo-z PDF for every galaxy, to weight quasar–galaxy pairs based on the probability of their being associated in redshift space. We describe the formalism briefly here, and refer the reader to M09 for further details.

For a set of spectroscopic quasars all at the same comoving distance χ_* from the observer, the angular cross-correlation between the (spectroscopic) quasars and (photometric) galaxies can be expressed in terms of the physical transverse comoving distance by (e.g., Shanks et al. 1983)

$$\begin{equation} w_\theta (R) = \frac{N_R}{N_G}\frac{D_Q D_G (R)}{D_Q R_G (R)} - 1, \end{equation} \tag{ 5 }$$

where R is the projected comoving distance for a given angular separation θ, such that R = χ_*θ. N_G and N_R are the total numbers of photometric galaxies and random galaxies, respectively, and D_QD_G and D_QR_G are the number of quasar–galaxy and quasar–random pairs, respectively, in each bin of R.

Defining the radial distribution function for the full galaxy sample as f(χ), where ∫f(χ)dχ = 1, and assuming that f(χ) varies slowly at the redshifts of interest, then the angular correlation function w_θ(R) is related to the projected real-space correlation function w_p(R) by

$$\begin{equation} w_\theta (R) = f(\chi _*)w_p(R) \end{equation} \tag{ 6 }$$

(for a derivation see Section 3.2 of Padmanabhan et al. 2009). As discussed in detail in M09, we can generalize the analysis such that the contribution to w_p(R) is calculated individually for each quasar–galaxy pair, with f_i,j defined as the average value of the radial PDF f(χ) for each photometric object i, in a window of size Δχ around the comoving distance to each spectroscopic object j. We use Δχ = 100 h⁻¹ Mpc to effectively eliminate redshift-space distortions, although the results are insensitive to the details of this choice.

In this case of weighting by pairs, we obtain, as in Equation (13) of M09:

$$\begin{equation} w_p(R) = N_R N_Q \sum _{i,j} c_{i,j} \frac{D_Q D_G (R)}{D_Q R_G (R)} - \sum _{i,j} c_{i,j}, \end{equation} \tag{ 7 }$$

where

$$\begin{equation} c_{i,j} = f_{i,j}{\bigg/}\sum _{i,j}f_{i,j}^2. \end{equation} \tag{ 8 }$$

We refer the reader to Section 2 of M09 for a detailed derivation and discussion of these equations.

M09 use Equation (7) to compute the cross-correlation between spectroscopic and photometric quasars from SDSS in the relatively narrow redshift bin 1.8 < z < 2.2, corresponding to comoving distances 3400 ≲ χ ≲ 3800 h⁻¹ Mpc; they obtain a cross-correlation length r₀ = 4.56 ± 0.48 h⁻¹ Mpc, assuming γ = 1.5. (Note that M09 do not derive γ = 1.5, they assume it purely to describe their method. Higher values of gamma are typically obtained in the recent literature, and we obtain γ ≈ 1.8 in the present work.) Our quasar sample spans a comparatively larger range in redshift (0.7 < z < 1.8, corresponding to 1750 ≲ χ ≲ 3400 h⁻¹ Mpc).

We evaluate Equation (7) by calculating the D_QD_G/D_QR_G term individually for each quasar. That is, for each quasar and each bin in separation R, we sum the redshift weights c_i,j for galaxies in the given range of distance from the quasar, and divide by the number of random galaxies in the same distance range (note that this implies N_Q = 1 in Equation (7)). The advantage of this procedure is that it consists of a simple sum and accounts exactly for the comoving distance to each quasar. However, the calculation is limited by shot noise on small scales where we have small numbers of quasar–galaxy and quasar–random pairs. To check that this does not significantly affect the results, we also divide the quasar sample into bins of width Δz = 0.1 (over which the comoving distance variations are small enough that there is little mixing between bins in R), and calculate the D_QD_G/D_QR_G term for all the N_Q quasars in the bin. We then average the w_p(R) values for the different bins to obtain a mean w_p(R) over the redshift range of interest. The resulting clustering amplitude differs by ≲10% (and in the majority of cases, < a few percent) compared to evaluating Equation (7) treating each quasar separately. The choice of method does not affect any of our conclusions, but to account for these differences we conservatively include an additional 10% systematic uncertainty on the measurement of the clustering amplitude. Finally, we emphasize that we are averaging w_p(R) over the whole redshift range of 0.7 < z < 1.8. The validity of this procedure depends on the fact that the observed w_p(R) varies slowly in the redshift range of interest, which we verify explicitly in Section 6.2.

To estimate DM halo masses for the quasars, we calculate the relative bias between quasars and galaxies from which we derive the absolute bias of the quasars relative to DM. As discussed below, calculation of absolute bias (and thus halo mass) requires a measurement of the autocorrelation function of the SDWFS galaxies. The large sample size enables us to derive the clustering of the galaxies accurately from the angular autocorrelation function ω(θ) alone. Although we expect the photometric redshifts for the SDWFS galaxies to be well constrained (as discussed in Section 3.2), by using the angular correlation function we minimize any uncertainties relating to individual galaxy photo-zs for this part of the analysis. The resulting clustering measured for the galaxies has much smaller uncertainties than that for the quasar–galaxy cross-correlation. To save computation time, for the galaxy autocorrelation analysis we use a significantly smaller random catalog with only 5 × 10⁵ random "galaxies." This likely introduces some additional shot noise into the calculation of ω(θ); however, since the resulting uncertainties are still far smaller than those for the quasar–galaxy cross-correlation, they are more than sufficient for the present analysis.

We calculate the angular autocorrelation function ω(θ) using the Landy & Szalay (1993) estimator

$$\begin{equation} \omega (\theta)=\frac{1}{{\it RR}}({\it DD}-2{\it DR}+{\it RR}), \end{equation} \tag{ 9 }$$

where DD, DR, and RR are the number of data–data, data–random, and random–random galaxy pairs, respectively, at a separation θ, where each term is scaled according to the total numbers of quasars, galaxies, and randoms.

The galaxy autocorrelation varies with redshift, owing to the evolution of large-scale structure, and because the use of a flux-limited sample means that we select more luminous galaxies at higher z. This will affect the measurements of relative bias between quasars and galaxies, since the redshift distribution of the quasars peaks at higher z than that for the galaxies and so relatively higher z galaxies dominate the cross-correlation signal. To account for this in our measurement of galaxy autocorrelation, we randomly select galaxies based on the overlap of the PDFs with the quasars in comoving distance (in the formalism of Section 4.2 this is f_i,j for each galaxy, averaged for all quasars). We select the galaxies so that their distribution in redshift is equivalent to the weighted distribution for all galaxies (weighted by 〈f_i,j〉). The redshift distribution of this galaxy sample is shown in Figure 6. We use this smaller galaxy sample to calculate the angular autocorrelation of SDWFS galaxies.

Zoom In Zoom Out Reset image size

In fields of finite size, estimators of the correlation function based on pair counts are subject to the integral constraint, which can be expressed as (Groth & Peebles 1977)

$$\begin{equation} \int\!\!\int \omega (\theta _{12}) {d}\Omega _1 {d}\Omega _2 \simeq 0, \end{equation} \tag{ 10 }$$

where θ₁₂ is the angle between the solid angle elements dΩ₁ and dΩ₂ and the integrals are over the survey area. If the number density fluctuations in the volume are small, and the angular correlations are smaller than the variance within the volume, then to first order the correlation function is simply biased low by a constant equal to the fractional variance of the number counts. A straightforward way to remove this bias is to add to the observed ω(θ) the term

$$\begin{equation} \omega _\Omega = \frac{1}{\Omega ^2}\int\!\!\int \omega (\theta _{12}) {d}\Omega _1 {d}\Omega _2, \end{equation} \tag{ 11 }$$

where Ω is the area of the survey region. The value of $\bar{n}^2 \omega _\Omega$, where $\bar{n}$ is an estimate of the mean number of galaxies per unit area, is the contribution of clustering to the variance of the galaxy number counts (Groth & Peebles 1977; Efstathiou et al. 1991). Evaluating Equation (11) for the Boötes survey area and the typical slope of the ω(θ) for the objects considered here, we obtain ω_Ω ≈ 0.03ω(1'). We estimate ω(1') by interpolating the observed ω(θ), then add ω_Ω to ω(θ) before performing model fits. For the projected real-space correlation functions w_p(R) (which is ultimately derived from individual estimates of ω(θ), as in Equation (5)), we perform an approximate correction for the integral constraint. We determine the value of w_p(R) at the physical scale (typically 0.5–1 h⁻¹ Mpc) corresponding to 1' for each quasar, and add the average of these estimates (multiplied by 0.03) to the observed w_p(R). These corrections increase the observed clustering amplitude by ≈10%, but have little effect on our overall conclusions.

Ideally, uncertainties in w_p(R) and ω(θ) would be determined by calculating the correlation function for various random realizations of mock IR-selected quasar and galaxy samples, for example by populating DM N-body simulations. In the absence of such mock catalog, we instead determine uncertainties in w_p(R) directly from the data through bootstrap resampling.

In a standard bootstrap analysis, the survey volume is divided into N_sub subvolumes, and these subvolumes are drawn randomly (with replacement) for inclusion in the calculation of the correlation function. Owing to the relatively small size of the field compared to large surveys such as SDSS or 2dF, we are only able to divide the field into a small number of subvolumes (we choose N_sub = 8). The width of one subvolume corresponds to ≈50 h⁻¹ Mpc at z = 1.2, so that correlations between the subvolumes should be relatively weak. (We verify explicitly that using a larger N_sub = 22 has no significant effect on the results.) For each bootstrap sample draw a total of three N_sub subvolumes (with replacement), which has been shown to best approximate the intrinsic uncertainties in the clustering amplitude (Norberg et al. 2009). We then recalculate w_p(R) including only the subvolumes in the bootstrap sample. For the calculations of w_p(R), we use 10,000 bootstrap samples for which the uncertainties at each scale converge to better than 1%. (To save computing time, we limit the analysis to 2000 bootstrap samples for the angular correlation analyses, for which the uncertainties converge to within ≈1.5%.)

This bootstrap technique works well for the galaxy autocorrelation, for which we have a large number of objects and the uncertainty is dominated by the clustering of the sample rather than by counting statistics. However, for the quasar–galaxy cross-correlation, the bootstrap analysis results in very small errors that are significantly smaller than the observed scatter between points. This appears to be caused by the fact that, owing to the small quasar samples of only a few hundred objects, the uncertainties are dominated by shot noise that is not fully characterized by randomly selecting entire subvolumes. To account for the shot noise, we therefore take the sets of three N_sub bootstrap subvolumes and randomly draw from them (with replacement) a sample of objects (quasars or galaxies) equal in size to the parent sample; only pairs including these objects are used in resulting cross-correlation calculation. This procedure yields a good estimate of the shot noise (the resulting χ²_ν ∼ 1) while also accounting for covariance due to the large-scale structure.

When fitting power-law models to the observed correlation functions, we compute parameters by minimizing χ², taking into account covariance between different bins in R. From the bootstrap analysis, we can estimate the covariance matrix C_ij by

$$\begin{eqnarray} C_{ij} &=& \frac{1}{1-N} \left[\sum _{k=1}^{N} \big (w_p^k(R_i)-w_p(R_i)\big)\right. \nonumber \\ &&\times \big (w_p^k(R_j)-w_p(R_j)\big)\bigg ] \,, \end{eqnarray} \tag{ 12 }$$

where w^k_p(R_i) and w^k_p(R_j) are the projected correlation function derived for the k-th bootstrap samples, N is the total number of bootstrap samples, and w_p(R) is the correlation function for the full sample. This formalism is equally valid for bins of angular separation θ in calculations of ω(θ). The 1σ uncertainty in each bin in R is the square root of the diagonal component of this matrix (σ_i = (C_ii)^1/2).

Taking into account covariance, χ² is defined as

$$\begin{eqnarray} &&\chi ^2 = \sum _{i=1}^{N_{\rm bins}}\sum _{j=1}^{N_{\rm bins}}\big (w_p(R_i)-w_p^{\rm model}(R_i)\big) \nonumber \\ &&\times \,C^{-1}_{ij} \big (w_p(R_j)-w_p^{\rm model}(R_j)\big), \end{eqnarray} \tag{ 13 }$$

where C⁻¹_ij is the inverse of the covariance matrix C_ij. We determine best-fit parameters by minimizing χ², and derive 1σ errors in each parameter by the range for which Δχ² = 1. As a check, we also estimate parameter uncertainties by calculating best-fit parameters for each of the bootstrap samples and calculating the variance between them; this obtains almost identical estimates of the errors. Furthermore, we note that if we use only the diagonal terms in the covariance matrix in determining χ², the variation in the best-fit parameters is significantly smaller than the statistical uncertainties, indicating that the precise details of the covariance matrix are relatively unimportant.

We also note that while in principle the SDWFS field is large enough to enable measurements of clustering up to ∼50 h⁻¹ Mpc at z ∼ 1, we limit the analysis to scales <12 h⁻¹ Mpc, because of edge effects that skew the correlation function on large scales but have minimal effect on smaller scales. An investigation of this effect is given in the Appendix.

For the projected real-space quasar–galaxy cross-correlation analysis, we fit power-law models of w_p(R) using Equation (4). We also fit power laws to the angular correlation functions (both galaxy autocorrelations and quasar–galaxy cross-correlations), using the simple expression

$$\begin{equation} \omega (\theta) = A\theta ^{-\delta }. \end{equation} \tag{ 14 }$$

For meaningful comparison to other clustering measurements obtained using samples with different distributions in redshift, we wish to convert the observed parameters A and δ to the real-space r₀ and γ as defined in Equation (3). Inverting Limber's equation, the conversion between these parameters can be computed analytically (here we follow Section 4.2 of Myers et al. 2006; for the full derivation see Peebles 1980)

$$\begin{eqnarray} \gamma &=& \delta +1 \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} A &=& H_\gamma \frac{\int ^{\infty }_0({d}N_1/{d}z)({d}N_2/{d}z) E_z \chi ^{1-\gamma }{d}z}{\big[\int ^{\infty }_0({d}N_1/{d}z){d}z\big]\big[\int ^{\infty }_0({d}N_2/{d}z){d}z\big]} r_{0}^\gamma, \end{eqnarray} \tag{ 16 }$$

where H_γ = Γ(0.5)Γ(0.5[γ − 1])/Γ(0.5γ), Γ is the gamma function, χ is the radial comoving distance, dN_1,2/dz are the redshift distributions of the samples (for an autocorrelation dN₁/dz = dN₂/dz), and E_z = H_z/c = dz/dχ. The Hubble parameter H_z can be found via

$$\begin{equation} H_z^2 = H_0^2[\Omega _m(1+z)^3+\Omega _\Lambda ]. \end{equation} \tag{ 17 }$$

Equation (16) assumes no evolution with redshift in the clustering of the sample (equivalent to the implicit assumption made in fitting w_p(R) with Equation (4)). For each angular correlation analysis, we derive A and δ from the observed ω(θ) and then obtain the corresponding γ and r₀ from Equations (15) and (16).

The masses of the DM halos in which galaxies and quasars reside are reflected in their absolute clustering bias relative to the DM distribution. To determine absolute bias (following e.g., Myers et al. 2007; Coil et al. 2008; Hickox et al. 2009) we first calculate the two-point autocorrelation of DM as a function of redshift. We use the HALOFIT code of Smith et al. (2003) to determine the nonlinear-dimensionless power spectrum Δ²_NL(k, z) of the DM assuming our standard cosmology, and the slope of the initial fluctuation power spectrum, Γ = Ω_mh = 0.21. The Fourier transform of the Δ²_NL(k, z) gives us the real-space correlation function ξ(r), which we then integrate to π = 100 h⁻¹ Mpc following Equation (2) to obtain the DM projected correlation function w^DM_p(R, z). The uncertainty in the DM power spectrum obtained from HALOFIT is ∼5%; this corresponds to a systematic uncertainty ∼0.05 dex in M_halo, but does not impact the relative halo masses of the different subsamples.

To derive quasar absolute bias from the projected real-space correlation function, we average the w^DM_p(R) over the redshift distribution of the sample, weighted by the overlap with the galaxy PDFs. The overlap of each quasar with the galaxy PDFs is given by

$$\begin{equation} W_i = \sum _{j} f_{i,j} \end{equation} \tag{ 18 }$$

and the corresponding w_p(R) for the DM is given by

$$\begin{equation} w_p^{\rm DM}(R) = \sum _i W_i w_p^{\rm DM}(R,z_i) \big/ \sum _i W_i, \end{equation} \tag{ 19 }$$

where z_i is the quasar redshift.

The redshift distributions for the QSO-1s and Obs-QSOs are essentially identical (the resulting w^DM_p(R) values for the two samples differ by <2% on all scales), so for simplicity we use the same w^DM_p(R) (defined for the QSO-1s) for both sets of quasars. We obtain the bias by calculating the average ratio between the best-fit power-law model and w^DM_p(R) over the range of scales of 1–10 h⁻¹ Mpc, for which w^DM_p(R) corresponds closely to a power law and is dominated by the two-halo term. The observed clustering amplitude relative to the DM corresponds to b_Qb_G, where b_Q and b_G are the absolute linear biases of the quasars and SDWFS galaxies, respectively.

To measure b_G from the galaxy autocorrelation function (or b_Qb_G from the quasar–galaxy angular cross-correlation, described in Section 6.1), we require an estimate of the corresponding ω(θ) of the DM. To obtain ω(θ), we use Limber's equation to project the power spectrum Δ²_NL(k, z) into the angular correlation (Limber 1953; Peebles 1980; Peacock 1991; Baugh & Efstathiou 1993). Specifically, we perform a Monte Carlo integration of Equation (A6) of Myers et al. (2007) to obtain ω(θ) for the DM. The key parameter in this equation is (dN_G/dz)² where dN_G/dz is the redshift distribution of the galaxies. We calculate dN_G/dz from the sum of the PDFs of the galaxies for which we perform the autocorrelation. In deriving the DM ω(θ) for the quasar–galaxy cross-correlation, we replace (dN/dz)² with (dN_Q/dz)(dN_G/dz), where dN_Q/dz is the distribution of quasar redshifts. For each angular correlation analysis, we compute the average ratio between the best-fit power-law model and the DM ω(θ) on scales 1'–10', where ω(θ) is dominated by the two-halo term. This ratio yields b²_G for galaxy autocorrelations or b_Qb_G for quasar–galaxy cross-correlations.

Finally, we use b_Q and b_G to estimate the characteristic mass of the DM halos hosting each subset of galaxies or quasars. Sheth et al. (2001) derive a relation between dark matter halo mass and large-scale bias that agrees well with the results of cosmological simulations. We use Equation (8) of Sheth et al. (2001) to convert b_abs to M_halo for the mean redshift of each subset of objects. If we use a different relation between b_abs and M_halo (Tinker et al. 2005), we obtain estimates for M_halo that are similar, although slightly larger by 0.2–0.3 dex; these differences do not significantly affect our conclusions.

We note that to estimate M_halo, we have performed fits to the observed w_p(R) on scales of 0.3–12 h⁻¹ Mpc. In principle, the DM and galaxy correlation functions can have somewhat different shapes such that the bias depends on the range of scales considered. If we limit the fits on scales 1–12 h⁻¹ Mpc, the results change by ≲5%, but with slightly larger uncertainties. We also note that our estimates of M_halo are relatively insensitive to our choice of σ₈. If we change σ₈ from 0.84 to 0.8 (as favored by the more recent Wilkinson Microwave Anisotropy Probe cosmology, e.g., Spergel et al. 2007) our M_halo estimates for quasars and galaxies increase by ≈0.1 dex.

We first calculate the cross-correlation of quasars and SDWFS galaxies using a simple angular clustering analysis, and check whether the corresponding absolute bias is consistent with that derived from the more sophisticated w_p(R) calculation. To calculate the ω(θ), we use an estimator corresponding to Equation (9) but for cross-correlations:

$$\begin{equation} \omega (\theta)=\frac{1}{RR}(D_QD_G-D_QR-D_GR+RR), \end{equation} \tag{ 20 }$$

where each term is scaled according to the total numbers of galaxies and randoms. To maximize the signal-to-noise ratio by cross-correlating objects associated in redshift space, the galaxies include only the redshift-matched SDWFS sample of 151,256 objects described in 4.3. Uncertainties are estimated using bootstrap resampling as described in Section 4.5. We fit the observed cross-correlation with a power law as described in Section 4.6. Owing to the limited statistics which provide only very weak constraints on the power-law slope, we fix δ = 0.8 (corresponding to real-space γ = 1.8).

The resulting cross-correlations and scaled DM fits are shown in Figure 9, and fit parameters are given in Table 1. The estimates of r₀ and b_Qb_G are in broad agreement between the two estimators, although as may be expected, the statistical uncertainties for the angular correlation analysis are larger (by ∼50%) than for the real-space analysis with fixed γ. Given that the absolute bias derived from the projected correlation function corresponds broadly to the bias from the noisier, but simpler angular cross-correlation, we conclude that there are no significant systematic effects that skew our estimate of w_p(R).

Zoom In Zoom Out Reset image size

Our calculation of the real-space quasar–galaxy correlation function over the redshift range 0.7 < z < 1.8 requires that w_p(R) varies only slowly between these redshifts, as discussed in Section 4.2. If the objects reside in similar halos at all redshifts, then we may expect r₀ to change slowly; simulations suggest that the typical r₀ for the autocorrelation of DM halos of mass ∼10¹²–10¹³ h⁻¹ M_☉ should change by ≲0.2 h⁻¹ Mpc between z = 0.5 and 2 (see Figure 10 of Starikova et al. 2010). To test explicitly the redshift variation for the clustering in our sample, we re-derived w_p(R) using the method outlined in Section 4 but selecting quasars over smaller redshift bins of 0.7 < z < 1.25 and 1.25 < z < 1.8. Uncertainties are calculated using the bootstrap method as for the full quasar samples, and DM and power-law fits are again performed over the range of separations 0.3 < R < 12 h⁻¹ Mpc. We evaluate w^DM_p(R) as in Section 4.7, but only including the quasars in the redshift ranges of interest. Owing to larger statistical errors and for simple comparison to the results over the full redshift range, in the power-law fits we fix γ to 1.8.

The resulting w_p(R) for the separate redshift bins and the power-law fits are shown in Figure 10. For the QSO-1s, we obtain r₀ = 5.1 ± 0.8 h⁻¹ Mpc and 6.3 ± 0.7 h⁻¹ Mpc for the low- and high-redshift bins, respectively, and for the Obs-QSOs we correspondingly obtain r₀ = 6.2 ± 0.7 h⁻¹ Mpc and 5.8 ± 1.0 h⁻¹ Mpc. The measured quasar–galaxy cross-correlation should be largely independent between the two redshift bins. Although the quasars are cross-correlated against the same galaxy sample in each bin, the galaxy samples will be weighted toward higher and lower redshifts in the high- and low-z bins, respectively. For the high- and low-redshift bins, the best-fit r₀ values bracket those for the full redshift samples, are broadly consistent within the uncertainties. Interestingly, the best-fit clustering amplitude for the Obs-QSOs increases with redshift while it decreases for the Obs-QSOs; however, given the uncertainties we decline to speculate on any possible difference in redshift evolution between the two subsets. Overall, the results in the different redshift bins confirm that any variation in the observed w_p(R) is sufficiently weak over the redshift range of interest, so that the method outlined in Section 4.2 should provide a reasonable estimate of the average clustering amplitude over the full redshift range.

Zoom In Zoom Out Reset image size

The primary uncertainty in our estimate of w_p(R) for the Obs-QSOs is the lack of accurate (that is, spectroscopic) redshifts and difficulty in estimating the photo-z uncertainties from the neural net calculations. As described in Section 4, in calculating w_p(R) for the Obs-QSO–galaxy cross-correlation, we simply assume that Obs-QSOs lie exactly at the best redshifts output by the neural net estimator. Any uncertainties in the photo-zs or systematic offsets from the true redshifts could therefore affect the resulting clustering measurement. The fact that we obtain very similar estimates of b_Q for the Obs-QSOs from a simple angular cross-correlation analysis as from our w_p(R) calculation suggests that uncertainties on individual photo-zs do not strongly affect our estimates of the quasar bias, as long as the overall distribution in redshifts for the Obs-QSOs is correct. However, it is possible that very large discrepancies from the true photo-zs, or any systematic shift in the redshift distribution, could affect both the estimate of w_p(R) and the real-space clustering parameters derived from ω(θ).

To precisely explore the effect of these errors, we take advantage of the fact that we have an equivalent sample of objects (the QSO-1s) that have a similar redshift distribution and for which the redshifts are known precisely from spectroscopy. We can therefore adjust the redshifts of the QSO-1s and re-calculate w_p(R) to determine how uncertainties or systematic shifts affect the observed correlation amplitude. As a simple first test, we calculate the w_p(R) for the QSO-1s using the photo-zs (as shown in Figure 4) rather than using spectroscopic redshifts. Figure 11 shows that the resulting w_p(R) differs little from that obtained using spectroscopic redshifts; the clustering amplitude for a power-law fit with fixed γ = 1.8 is lower by 12%. Note that if we allow γ to float, the average bias for the photo-z sample is actually larger by ≈10% (owing to a slightly flatter slope) but well within the statistical uncertainties. We conclude that for the QSO-1s, photo-z errors do not have a significant effect on our measurements of the clustering.

Zoom In Zoom Out Reset image size

6.3.1. Random Errors

To explore photo-z errors in more detail, we systematically test the effects of Gaussian random errors in the quasar redshifts. For each quasar, we shift the best estimates of redshift (spec-zs for QSO-1s and photo-zs for Obs-QSOs) by offsets Δz/(1 + z) selected from a Gaussian random distribution with dispersion σ_z/(1 + z). Using these new redshifts, we recalculate w_p(R), using the full formalism described in Section 4. We perform the calculation 10 times for each of several values of σ_z/(1 + z) from 0.02 up to 0.2 (which smears the redshifts across most of the redshift range of interest). To ensure that this step does not artificially smear out the redshift distribution beyond the range probed by the galaxies, we require that the random redshifts lie between 0.7 < z < 1.8; any random redshift that lies outside this range is discarded and a new redshift is selected from the random distribution. For each trial, we obtain the relative bias by calculating the mean ratio of w_p(R), on scales 1–10 h⁻¹ Mpc, relative to the w_p(R) for the best estimates of redshift. We then average the 10 trials at each σ_z to obtain a relation between relative bias and σ_z, shown in Figure 12(a).

Zoom In Zoom Out Reset image size

As may be expected, Figure 12(a) shows that shifting the QSO-1 redshifts from their true values causes a decrease in the cross-correlation amplitude, as the quasars are preferentially correlated with galaxies that are not actually associated in redshift space. We find a monotonic decrease in relative bias with σ_z, from ≈0.95 for σ_z/(1 + z) = 0.05 to ≈0.8 for σ_z/(1 + z) = 0.2. Repeating this calculation for the Obs-QSOs reveals a very similar trend. The decrease in bias with σ_z/(1 + z) shown in Figure 12(a) indicates that such errors would affect the measurements of the clustering amplitude by at most ∼20%.

6.3.2. Systematic Shifts

While the above analysis suggests that random errors in the Obs-QSO photo-zs do not strongly affect the observed clustering amplitude, it is also possible that systematic uncertainties in the photo-z (consistent over- or under-estimates of the redshift) could significantly alter the observed bias. To test this, we shift the redshifts of the quasars as discussed in Section 6.3.1, but in place of random shifts, we compress all redshifts toward one end or the other of the 0.7 < z < 1.8 range. (This procedure allows us to test the effects of systematic shifts in redshift while keeping the same overall redshift range.) The shift in redshift is defined by a redshift scaling parameter S_z, such that

$$\begin{eqnarray} z_{\rm new} = \left\lbrace \begin{array}{@{}l l}z+S_z(z-0.7) & \quad S_z < 0\\ z+S_z(1.8-z) & \quad S_z \ge 0. \\ \end{array} \right. \end{eqnarray} \tag{ 21 }$$

As an additional check we also perform a simple linear offset Δz of the redshifts, allowing the redshifts to move outside the selection range of 0.7 < z < 1.8. As in Section 6.3.1, we use these new redshifts to recalculate w_p(R) via the full formalism described in Section 4, and determine the relative bias on scales 1–10 h⁻¹ Mpc. Relative bias versus S_z and Δz are shown in Figures 12(b) and (c). For the QSO-1s, the peak of the observed clustering amplitude is very close to S_z = 0, while shifting the redshifts down or up systematically decreases the bias. The Obs-QSOs show a similar peak near S_z = 0, indicating that the Obs-QSO photo-zs are not systematically offset higher or lower than the true redshifts by a large factor. We note that for the Obs-QSOs, a slight shift to higher redshifts (Δz = 0.1) increases the clustering by a small amount (≈8%).

Finally, we emphasize that any possible low-redshift contaminants (such as star-forming galaxies, as discussed in Section 7 of H07), will serve only to decrease the observed clustering signal, as they will be completely uncorrelated in angular space with the higher redshift SDWFS galaxies that lie in entirely separate large-scale structures. Therefore, the observed w_p(R) represents a robust lower limit to the clustering amplitude for the Obs-QSOs.

We compare our observed absolute bias for the mid-IR-selected quasars with the clustering of optically selected (Type 1) quasars, which has been well established by a number of works. Among the most precise measurements to date are studies that have used data from the 2dF and SDSS surveys. Recently, Ross et al. (2009) measured the evolution of the quasar three-dimensional autocorrelation function based on spectroscopic quasars from SDSS Data Release 5, and compared to previous results from spectroscopic samples from the 2QZ (Croom et al. 2005) and 2SLAQ (da Ângela et al. 2008), as well as the clustering of photometrically selected quasars from SDSS (Myers et al. 2006). Figure 13(a) shows the redshift evolution of linear bias for Type 1 quasars from these studies, taken from Figure 12 of Ross et al. (2009). Where appropriate, we have converted the bias to our adopted cosmology using the formalism in the Appendix of Starikova et al. (2010), assuming a shape factor s = 1 for simplicity. Figure 13(b) shows the corresponding estimates of characteristic M_halo derived from the linear bias using the prescription of Sheth et al. (2001). The linear bias and halo masses for the QSO-1s and Obs-QSOs are shown for comparison.

Zoom In Zoom Out Reset image size

It is readily apparent from Figure 13(a) that the observed bias of Type 1 quasars increases with redshift, as discussed in Section 1. (These results are also consistent with a number of other quasar clustering studies using other data see, e.g., Figure 15 of Hopkins et al. 2008.) The dashed curves in Figure 13(a) show the increase in bias with redshift for halos of constant mass, clearly showing that at all redshifts the QSO-1s reside in DM halos of roughly a few ×10¹² h⁻¹ M_☉. Our measurement for the QSO-1s is in excellent agreement with the evolution in linear bias and roughly constant M_halo observed in previous measurements of Type 1 quasar clustering.

For the Obs-QSOs, the best-fit bias is marginally larger (≈1σ), corresponding to a factor of roughly four difference in the characteristic M_halo between the QSO-1s and Obs-QSOs (Figure 13(b)). As discussed above, random errors in the photo-zs can only decrease the observed clustering amplitude and the inferred M_halo. Thus, our measurement of Obs-QSO clustering represents a lower limit, and it is possible that the true Obs-QSO bias is somewhat higher (although the results of Section 6 suggest that the true bias may be larger by at most ∼20%). Based on this analysis we can make the robust conclusion that the Obs-QSOs are at least as strongly clustered as their QSO-1 counterparts.

Stronger clustering for obscured quasars would have significant implications for physical models of the obscured quasar population. In terms of unified models, a difference in clustering between obscured and unobscured quasars would rule out the simplest picture in which obscuration is purely an orientation effect, but may be consistent with more complicated scenarios where the effective covering fraction changes with environment. Alternatively, if obscuration is caused by large (∼kpc scale) structures, then the processes that drive these asymmetries (e.g., mergers, disk instabilities, accretion of cold gas from the surrounding halo) may be more common in halos of larger mass. Indeed, given that some fraction of quasars might naturally be expected to be obscured by a "unified"-model torus, any observed differences in clustering may reflect even stronger intrinsic dependence of large-scale obscuration and environment.

An intriguing scenario for obscured quasars is that they represent an early evolutionary phase of rapid black hole growth before a "blowout" of the obscuring material from the central regions of the galaxy and the emergence of an unobscured quasar (e.g., Figure 1 of Hopkins et al. 2008). Quasars tend to radiate at large fraction of the Eddington rate (McLure & Dunlop 2004; Kollmeier et al. 2006; although see Kelly et al. 2010), so that the similar L_bol for QSO-1s and Obs-QSOs would imply that they host black holes of similar masses. Any correlation (e.g., Ferrarese 2002; Booth & Schaye 2010), even if indirect (e.g., Kormendy & Bender 2011), between the final masses of black holes and those of their host halos would thus suggest that our obscured and unobscured quasars would have the same M_halo, as long as their black holes are near their final masses. However, if obscured quasars are in an earlier phase of rapid growth and so are in the process of "catching up" to their final mass (e.g., King 2010), then they would have a larger M_halo compared to unobscured quasars with the same M_BH. In light of recent debate as to whether black holes generally grow before or after their hosts (e.g., Peng et al. 2006; Alexander et al. 2008; Woo et al. 2008; Decarli et al. 2010), this scenario would imply that black hole growth lags behind that of the host halo.

In any physical picture, a significant difference in clustering between obscured and unobscured quasars would also imply a difference in accretion duty cycles (or equivalently, lifetime). QSO-1s and Obs-QSOs are found in roughly equal numbers, but the abundance of DM halos drops rapidly with mass (e.g., Jenkins et al. 2001), thus implying that if one type of quasars are found in larger halos then they must be longer-lived. With our current results, we are able to rule out any model in which obscured quasars are substantially less strongly clustered or have shorter lifetimes than their unobscured counterparts. With more accurate future measurements, detailed studies of halo masses and lifetimes for obscured and unobscured quasars could place powerful constraints on evolutionary scenarios such as those described above.

Our results demonstrate the potential for studying the clustering of obscured quasars in extragalactic multiwavelength surveys, and the marginally significant difference in clustering we observe for obscured and unobscured quasars provides strong motivation for more precise measurements in the future. The two main avenues for progress are improvements in redshift accuracy and selection of larger samples for better statistical accuracy. Upcoming sensitive, wide-field multi-object spectrographs will enable efficient measurements of redshift for large numbers of optically faint sources and so improve calibrations of obscured quasar photo-zs, or with large enough samples enable fully three-dimensional clustering studies. In addition, we will soon have the capability to detect many thousands of obscured quasars based on very wide-field observations in the mid-IR with the Wide-Field Infrared Survey Explorer (Wright 2010) and in X-rays with eROSITA (Predehl et al. 2007) or the Wide-Field X-ray Telescope (Murray et al. 2010). These data sets will allow us to measure obscured quasar clustering with statistical precision that is comparable to current measurements of unobscured quasars.