THE HALO MASSES AND GALAXY ENVIRONMENTS OF HYPERLUMINOUS QSOs AT z ≃ 2.7 IN THE KECK BARYONIC STRUCTURE SURVEY^*

An important prerequisite to establishing the galaxy environment of the HLQSOs is an accurate measurement of the HLQSO systemic redshifts. Redshifts for QSOs in the range 2 ≲ z_QSO ≲ 3 are typically measured from the peaks or centroids of broad emission lines of relatively high ionization species in the rest-frame far-UV (e.g., N v λ1240, C iv λ1549, Si iv λ1399, C iii] λ1909). These lines are known to yield redshifts that differ significantly from systemic, and tend to be blueshifted by several hundred to several thousand km s⁻¹ (see, e.g., McIntosh et al. 1999; Richards et al. 2002; Gonçalves et al. 2008). These velocity offsets also tend to increase with QSO luminosity, thus making the present sample of HLQSOs particularly susceptible to this issue. In view of the importance of precise redshifts to locate the HLQSO environments within the survey volume, we obtained near-IR spectra of the entire sample using NIRSPEC on the Keck II 10 m telescope, TripleSpec on the Palomar 200 inch (5 m) telescope, and in some cases, both (see Table 1).

Among the 15 HLQSOs in the sample, narrow forbidden lines ([O iii] λ5007) were detected for only two of them (Q1623+268, Q2343+12), either because no such lines were present in the spectra (common at the highest luminosities), or because the HLQSO redshift was such that the strongest transitions fell in regions between the near-IR atmospheric bands. However, in all cases we were able to measure one or more hydrogen Balmer lines and the Mg ii λ2798 line, which were mutually consistent and are known to be closer to the true systemic redshift than the high ionization lines in the UV (McIntosh et al. 1999; Richards et al. 2002). The redshifts obtained from the Balmer/Mg ii lines (which agree well with that given by [O iii] in the two cases where all were measured) were then subjected to several cross-checks, including: the wavelength at the onset of the Lyα forest measured from the high-resolution HLQSO spectra (a lower limit on the systemic redshift, but one which agrees to within Δz ≃ 0.001 of the Balmer line redshift in all but two cases); the redshift of narrow He ii λ1640 in intermediate-resolution optical spectra of the HLQSOs; and, in several cases, regions exhibiting narrow Lyα emission were discovered with small angular separations from the HLQSO, and we have found that such nebulae lie very close to the systemic redshift of the nearby HLQSO. In two cases (HS1603+3820 and Q1009+29), this last criterion led to a significant modification (Δz ∼ +0.01, or ∼800 km s⁻¹) of the redshift suggested by the near-IR spectroscopy. We adopt an HLQSO redshift uncertainty σ_z = 270 km s⁻¹ (for those HLQSOs without measured [O iii] redshifts) based on the measured dispersion of the Mg ii line with respect to [O iii] by Richards et al. (2002); the broad agreement among our many redshift criteria suggest that this is a conservative estimate of the redshift uncertainties.

Table 1 summarizes the adopted redshifts for all 15 HLQSOs based on these considerations; also given (Column 4; z_old) is the published redshift for each and the redshift and velocity error that would result from adopting the published values (Δz ≡ z_old − z_new). As expected, all but one of the old redshifts are systematically too low (the median shift is ∼ − 1500 km s⁻¹, and the mean ∼ − 2100 km s⁻¹). Failure to account for these large velocity errors would severely compromise our measurements. As we show below, the measured z_QSO values must be quite accurate given the very tight redshift-space correlation between the HLQSOs and the spectroscopically measured, continuum-selected galaxies nearby.

Galaxy redshifts were measured using low-resolution (∼5 Å), rest-frame UV spectra obtained with the LRIS multi-object spectrograph on the Keck I telescope (Oke et al. 1995; Steidel et al. 2004). Candidate galaxies were color-selected using the Lyman break technique and were sorted as BX (z ∼ 2.2), MD (z ∼ 2.6), or CDM (z ∼ 3) galaxies based on the color criteria discussed in Steidel et al. (2003) and Adelberger et al. (2004); the data collection and reduction procedures are described therein. All galaxies in the spectroscopic sample have $\mathcal {R} < 25.5$ (where $\mathcal {R} \equiv m_{{\rm AB}}$(6830 Å)), which corresponds to M_AB(1700 Å) ≲ −19.9 at z ∼ 2.7 (about 1 mag fainter than M_* at this redshift; see Reddy et al. 2008). Redshifts were determined by a combination of Lyα emission or absorption and far-UV interstellar (IS) absorption. Since Lyα emission tends to be redshifted with respect to the systemic redshift of the host galaxy, and IS absorption tends to be blueshifted (see, e.g., Shapley et al. 2003; Adelberger et al. 2003; Steidel et al. 2010), we estimate each galaxy's systemic redshift via the method proposed in Adelberger et al. (2005a) and updated by Steidel et al. (2010). In this method, the average Lyα emission and IS absorption offsets are calculated based on the redshift of the Hα nebular line (NIR spectroscopy is available for a subset of the galaxy sample), which traces ionized gas in star-forming regions of the galaxy, and is thus a more accurate estimate of the systemic redshift. Rakic et al. (2011) derive similar corrections for the same galaxy sample using the expected symmetry of intergalactic medium absorption about the systemic redshift of the galaxy.

We estimate the systemic galaxy redshifts (z_gal) based on a combination of the above results. A more detailed discussion of our correction formulae can be found in Rudie et al. (2012), but the formulae are reproduced below. For galaxies with NIR spectra (e.g., the Hα line), the NIR redshift is used with no correction. For galaxies with measured Lyα emission but without IS absorption, we use the following estimate:

$$\begin{equation} z_{\rm gal} \equiv z_{\rm Ly\alpha }+\frac{\Delta v_{{\rm Ly \alpha }}}{c}(1+z_{{\rm Ly\,\alpha }}), \end{equation} \tag{ 1 }$$

where z_Lyα is the redshift of the measured Lyα emission and Δv_Lyα = −300 km s⁻¹ is the velocity shift needed to transform the Lyα redshift to the systemic value, z_gal.

For galaxies with IS absorption, we use an estimate based on the absorption redshift whether or not Lyα emission is present:

$$\begin{equation} z_{{\rm gal}} \equiv z_{{\rm IS}}+\frac{\Delta v_{{\rm IS}}}{c}\left(1+z_{{\rm IS}}\right), \end{equation} \tag{ 2 }$$

where z_IS is the redshift of the measured IS absorption and Δv_Lyα = 160 km s⁻¹ is the velocity shift needed to transform the absorption redshift to the systemic value.

For galaxies with both IS absorption and Lyα emission, we verify that the corrected absorption redshift does not exceed the measured redshift of the Lyα line; that is, we verify that z_IS < z_gal < z_Lyα, where z_gal is calculated using Equation (2) above. If this condition is not satisfied, we recompute the galaxy systemic redshift as the average of the absorption and emission redshifts:

$$\begin{equation} z_{{\rm gal}} \equiv \frac{z_{{\rm IS}}+z_{{\rm Ly\alpha }}}{2}. \end{equation} \tag{ 3 }$$

The residual redshift errors (calculated from the galaxies in the NIR sample) have a standard deviation σ_{v, err} = 125 km s⁻¹, which we adopt as the uncertainty in our galaxy redshift measurements.

Clustering measurements can be grossly misinterpreted when the relevant selection functions are not well understood (Adelberger 2005). While the criteria for selecting galaxies for follow-up spectroscopy were identical for all 15 of the KBSS fields, small differences in image depth and seeing, as well as slight changes in the algorithms used for assigning relative weights in the process of designing slit masks, can lead to field-to-field variations in the redshift selection functions. To at least partially mitigate such variations in the redshift-space sampling between fields, we used the number of successfully observed BX, MD, and CDM galaxies in each field to estimate the form of our field-specific selection functions $\mathcal {N}_{j}(z)$. These estimates of the selection functions were constructed as follows.

First, the redshift distributions of all BX, MD, and CDM galaxies in our sample were arranged in a coarse histogram with bins of width Δz = 0.2. A spline fit was then performed to estimate the smooth distribution functions of each galaxy type—the histograms and spline fits for each type are displayed in Figure 2.

Zoom In Zoom Out Reset image size

For each field, we built a field-specific selection function by combining these galaxy redshift distributions for each color criterion according to the number of those galaxies successfully observed in the field. Thus, for a field with index j, the redshift selection function is given by Equation (5):

$$\begin{eqnarray} \mathcal {N}_j(z) & = & N_{{\rm BX},j}\,\mathcal {N}_{{\rm BX}}(z)+ N_{{\rm MD},j}\,\mathcal {N}_{{\rm MD}}(z)+N_{{\rm CDM},j}\,\mathcal {N}_{{\rm CDM}}(z),\nonumber\\ \end{eqnarray} \tag{ 5 }$$

where N_{BX, j} corresponds to the number of BX-selected galaxies in field j, $\mathcal {N}_{{\rm BX}}$ is the selection function for BX-selected galaxies over all fields, and other variables are defined similarly. We then transform these redshift-space selection functions into units of velocity relative to their corresponding bright HLQSOs using Equation (4). Finally, we combined this set of field-specific velocity-space selection functions (already weighted by the number of galaxies in each field) into a single stacked function:

$$\begin{equation} \mathcal {N}(v)=\sum _{j=1}^{15}\mathcal {N}_j(v). \end{equation} \tag{ 6 }$$

The resulting selection function is fairly flat over the range |δv| < 20, 000 km s⁻¹ with a slight negative slope (yellow shading in Figure 1), indicating our slight bias toward detecting objects "in front" of the HLQSO (that is, at lower redshifts) compared to galaxies slightly "behind" the HLQSO in each field. The selection function is thus a prediction for the observed distribution of galaxies in relative-velocity space in the absence of clustering.

We previously knew one KBSS field (HS1549+1919) to have a large overdensity in the galaxy distribution very close to the redshift of the central HLQSO. The variation in overdensity among fields can be estimated by N₁₅₀₀ in Table 2, which is the number of galaxies within 1500 km s⁻¹ of the HLQSO redshift for that field. In order to ensure that our clustering results are not being dominated by a single field, we repeated our analysis on subsamples of the data consisting of 14 of the 15 fields, removing a different field each time. In each case, the magnitude and scale of the overdensity was consistent with that observed when all 15 fields were included in the analysis, indicating that the observed magnitude and scale of the overdensity are not determined by any single field.

Figure 1 shows the observed galaxy distribution in units of velocity along with the selection function estimate from Section 3.1. The peak in the galaxy distribution near the HLQSO redshifts is clearly visible. Fitting a Gaussian function to the histogram in Figure 1 gives a velocity width σ_{v, fit} = 350 ± 50 km s⁻¹, which includes the effect our σ_{v, err} ≃ 125 km s⁻¹ galaxy redshift errors and the random residual errors in our HLQSO redshifts, assumed to be σ_{v, err} ∼ 270 km s⁻¹. After subtracting the redshift errors in quadrature, we find an intrinsic velocity width of σ_{v, pec} ≃ 200 km s⁻¹ for the galaxy overdensity, which we attribute to peculiar velocities. Note that the residual HLQSO redshift errors are uncertain and likely to be largely systematic (see Section 2.1), so our estimated velocity dispersion is an upper limit on the true peculiar velocity scale if the random component of the HLQSO redshift error is larger than we have assumed.

We also consider the relative overdensity at the HLQSO redshift by comparing the observed density to that predicted by our selection function. The distribution is plotted as a relative overdensity

$$\begin{equation} f_{{\rm ovr}}=(N_{{\rm obs}}-N_{{\rm pred}})/N_{{\rm pred}} \end{equation} \tag{ 7 }$$

in Figure 3, where N_obs is the number of galaxies observed in a given velocity bin and N_pred is the number predicted for that bin by our selection function. The relative overdensity is measured in bins with Δv = 500 km s⁻¹; this scale was chosen to correspond roughly to the transverse scale of our field, since a Hubble flow velocity of 500 km s⁻¹ ∼ 5 h⁻¹ cMpc at these redshifts. Figure 3 shows that the HLQSOs are associated (on average) with a δn/n ∼ 7 overdensity of galaxies when considered on the ∼5 h⁻¹ Mpc scale of our field, with no features of comparable amplitude over a wide range of redshifts (40,000 km s⁻¹ corresponds to Δz ≃ 0.5 at z ≃ 2.7).

Zoom In Zoom Out Reset image size

Repeating this analysis after dividing the galaxies into radial annuli, we find that the redshift association is most pronounced for those galaxies within 25''  of the HLQSO line of sight (∼200 pkpc), though a lower level of redshift clustering does extend to larger projected distances (see Figure 4). If this distance is taken as an isotropic spatial scale of the galaxy overdensity, then the line-of-sight velocity dispersion due to the Hubble flow would be only ∼65 km s^-1. However, a less-significant overdensity does extend to larger radii, and thus likely includes many galaxies that are clustered around the HLQSO but move with the Hubble flow. In order to ensure that the measured velocity width is not inflated by these non-virialized galaxies, we directly measure the velocity dispersion among the 15 galaxies within 1500 km s⁻¹ and 0.5 h⁻¹ cMpc (200 pkpc) of the HLQSOs; as discussed in Section 5, these galaxies are likely to be virialized and associated with the HLQSO, and our selection functions predict only 1.5 galaxies in this volume in the absence of clustering. These 15 galaxy velocities have a sample standard deviation of 335 km s⁻¹, consistent with the velocity width measured for the entire overdensity. The observed velocity spread is thus presumably set by peculiar velocities of σ_{v, pec} ≃ 200 km s^-1 among the HLQSO-associated galaxies.

Zoom In Zoom Out Reset image size

A comparison of Figures 3 and 4 demonstrates that the relative overdensity is highly scale dependent. If we assume that the width of the overdensity in velocity space is entirely due to peculiar velocities, and hence that all 15 of the galaxies observed with R < 0.5 h⁻¹ cMpc and |δv| < 500 km s⁻¹ are physically located within a three-dimensional distance r < 0.5 h⁻¹ cMpc from their nearest HLQSO, then the number of galaxies in this composite volume is ∼50× the number predicted by our redshift and angular selection functions (described in Sections 3.1 and 4.1.1, respectively).

Much of the recent work on QSO clustering relies on large-scale two-point correlation functions, particularly the QSO autocorrelation function (see, e.g., Shen et al. 2007). The galaxy–HLQSO cross-correlation function ξ_Q can provide a complementary estimate of HLQSO host halo mass.

The correlation function is defined as the excess conditional probability of finding a galaxy in a volume dV at a distance r = |r₁ − r₂|, given that there is an HLQSO at point r₁, such that P(r₂|r₁)dV = P₀[1 + ξ(r)]dV, where P₀dV is the probability of finding a galaxy at an average place in the universe. Here we assume a power-law form for the correlation function: ξ_GQ = (r/r^GQ₀)^−γ, where γ is the slope parameter and r₀ corresponds to the comoving distance at which the local number density of galaxies is twice that of an average place in the universe.

Many recent analyses of the two-point correlation function have dispensed with power-law fits in favor of directly modeling the halo-occupation distribution (HOD; see, e.g., Seljak 2000; Berlind & Weinberg 2002; Zehavi et al. 2004) based on the theory of Press & Schechter (1974) and a statistical method of populating DM halos with galaxies. A general feature of these HOD models is a deviation from a single power law at distances near 1 h⁻¹ cMpc due to a transition from the single-halo regime (the clustering of galaxies/QSOs within a single DM halo) to the two-halo regime (the clustering of galaxies/QSOs hosted by distinct halos).

In this paper, we implement the simpler power-law fitting technique for the following reasons. First, our smaller sample (with respect to the large surveys at low redshift) does not allow us to detect a deviation from a power-law fit with any significance, particularly for the galaxy–HLQSO cross-correlation. Second, our choice to fix the power-law slope γ (see below) desensitizes our result to the precise shape of the correlation function, leaving the clustering length r^GQ₀ to primarily reflect the integrated pair-probability excess over the range of projected distances in our sample.

In practice, the three-dimensional correlation function ξ(r) is not directly measurable: line-of-sight velocities are an imperfect proxy for radial distance due to peculiar velocities and redshift errors. As such, it is more useful to consider the reduced angular correlation function, w_p(R|Δz) by integrating over a redshift or velocity window:

$$\begin{equation} P(R)d\Omega =P_{0}^\prime d \Omega [1+w_p(R)]= d \Omega \int _{\Delta z} P(r) dz, \end{equation} \tag{ 8 }$$

where R = D_A(z)θ(1 + z) is the projected comoving distance from the HLQSO, and D_A(z) is the angular diameter distance to the HLQSO. In the limit Δz → ∞, and assuming a power-law form of the three-dimensional correlation function, it can be shown that the reduced angular correlation function has an equally simple power-law form:

$$\begin{equation} w_p(R)=A R^{-\eta }. \end{equation} \tag{ 9 }$$

with η = γ − 1. However, we would like to restrict our analysis to small redshift/velocity scales, choosing a value of Δz that includes the entire clustering signal while eliminating the noise contribution of uncorrelated structure at large line-of-sight separations from the HLQSO, and this priority precludes the assumption of Δz → ∞. In the case of a truncated redshift range, the reduced angular correlation function does not simplify to a power law, and instead takes the form of a Gaussian hypergeometric function, denoted as ₂F₁(a, b; c; z). In our particular case, the reduced angular cross-correlation function w^GQ_p is expressed by the following:

$$\begin{eqnarray} w_p^{{\rm GQ}}(R) & = & \int _{-z_0}^{z_0} (r/r_{0}^{{\rm GQ}})^{-\gamma } dz = \int _{-z_0}^{z_0} \big(\sqrt{R^2+z^2}/r_{0}^{{\rm GQ}}\big)^{-\gamma } dz \nonumber \\ & = & \left(\frac{r_0^{{\rm GQ}}}{R}\right)^\gamma {}_2F_1\left(\frac{1}{2},\frac{\gamma }{2};\frac{3}{2};\frac{-z_0^2}{R^2}\right), \end{eqnarray} \tag{ 10 }$$

where z₀ is the half-width of the redshift window over which we compute the clustering strength. We choose a value z₀ = 1500 km s⁻¹ ≃ 14 h⁻¹ cMpc in order to encompass the entire observed overdensity (see Figure 1) and the range of projected distances we are able to probe (R < 4.2 h⁻¹ cMpc) without including excess noise. We then fit the reduced angular correlation function w^GQ_p(R|r₀^GQ, γ) to the data by variation of the correlation length r^GQ₀. We fix γ = 1.5 for simplicity in matching our data to halo populations (see Section 4.3); this value of γ was chosen as it is a reasonably good fit to both the galaxy autocorrelation function and galaxy–HLQSO cross-correlation function, as well as the correlation functions among the simulated halo populations. Increasing the value of γ causes the best-fit value of r^GQ₀ to decrease, but the corresponding halo mass is very insensitive to the choice of γ, so long as the same value is used for both the galaxy and the simulated halo populations.

In order to estimate w^GQ_p(R) from our data, we separate our fields into projected circular annuli of varying widths, constructed so that each annulus has a roughly similar signal-to-noise ratio, and with our largest annulus having its outer edge ∼200''  from the HLQSO, a projected distance of R = 4.2 h⁻¹ cMpc. As noted above, we wish to restrict our analysis to those galaxies that are closely associated with an HLQSO in redshift as well as projected position, so we also separate our galaxy sample into two velocity groups: one with |δv| ⩽ 1500 km s⁻¹ and one with |δv| > 1500 km s⁻¹. In this way, we define N_v(R_k) as the number of velocity-associated galaxies in the kth annular bin and N₀(R_k) as the number of non-associated galaxies in the bin.¹ The solid angle covered by the kth annular bin is given by $\mathcal {A}_k=\pi (R_{{\rm outer},k}^2-R_{{\rm inner},k}^2)$; thus, we likewise define the area densities of galaxies $\Sigma _{v,0} (R_k) = N_{v,0}(R_k)/\mathcal {A}_k$.

We then used the selection function constructed in Section 3.1 to estimate the expected number of galaxies in each velocity group (i.e., $\mathcal {N}_v=\mathcal {N}(|\delta v|<1500)$ and $\mathcal {N}_0=\mathcal {N}(|\delta v|>1500)$), which we convert to expected average area densities for each velocity group ($\Sigma _{v,0}^{{\rm pred}}=\mathcal {N}_{v,0}/\mathcal {A}_{{\rm field}}$). Finally, we divide the measured area densities by the average predicted value to define the relative overdensity of each annulus. If the overdensity Σ_v(R_k)/Σ^pred_v is purely due to the clustering signal, then the reduced angular cross-correlation function is given by w^GQ_p = Σ_{v, obs}^GQ(R_k)/Σ^GQ_{v, pred} − 1; however, this assumption is invalid if the angular sampling of the field is not uniform, which we explore below.

4.1.1. Angular Selection Function

The use of the redshift selection function in N_pred ensures that large-scale variations and sampling biases in redshift are taken into account in our analysis. Selection biases can also occur in the plane of the sky; because our fields are centered on their HLQSOs, any bias that varies with distance from the center of the field will mimic a change in the correlation function. To account for this effect, we recall that galaxies with |δv| > 1500 km s⁻¹ show no association with the HLQSO (see Figure 1), and therefore should be uniformly distributed on average. Therefore, if the function Σ₀(R) is not a constant, it must describe a non-uniform angular selection function, which encapsulates variations in optical selection sensitivity (e.g., due to non-uniform extinction or field coverage) as well as any biases in slit positions on our masks. We assume these biases are independent of redshift, and that they produce the same fractional excess of galaxy counts in all velocity bins. Therefore, the measured values of Σ_v(R)/Σ^pred_v correspond to the true reduced correlation function 1 + w_p(R) multiplied by a transverse (angular) selection function, which we estimate by Σ₀(R_k).

$$\begin{equation} \frac{\Sigma _v(R_k)}{\Sigma _v^{{\rm pred}}} = \frac{\Sigma _0(R)}{\Sigma _0^{{\rm pred}}}[1+w_p (R_k)]. \end{equation} \tag{ 11 }$$

We found that the measured values of Σ₀(R_k) are well matched by a power law in R, and therefore, rather than using Equation (11) directly to estimate w_p(R_k), we found best-fit parameters α and β for the following model:

$$\begin{equation} \frac{\Sigma _0(R)}{\Sigma _0^{{\rm pred}}} = \alpha R^\beta. \end{equation} \tag{ 12 }$$

The best-fit parameters for this model are α = 1.59, β = −0.58 with R in h⁻¹ cMpc; the fit selection function is displayed in Figure 5. Combining Equations (10), (11), and (12), we arrive at an explicit model for Σ_v(R) in terms of the galaxy–HLQSO cross-correlation length r₀ and correlation slope γ:

$$\begin{equation} \frac{\Sigma _v(R_k)}{\Sigma _v^{{\rm pred}}} = \alpha R^\beta \left[1+\left(\frac{r_0}{R}\right)^\gamma {}_{2}F_1\left(\frac{1}{2},\frac{\gamma }{2};\frac{3}{2};\frac{-z_0^2}{R^2}\right)\right], \end{equation} \tag{ 13 }$$

where z₀ = (1500 km s⁻¹)H⁻¹₀(1 + z)⁻¹ is the half-width of the redshift window in physical units, and again α and β are set by fitting Σ₀(R) to Σ₀(R_k). We then adjust the free parameter r₀ corresponding to the cross-correlation function to fit the measured values of Σ_v(R_k).

Zoom In Zoom Out Reset image size

The fit to Σ_v was performed via a simple χ²-minimization using an error vector constructed assuming Poisson uncertainties in the galaxy counts. The binned data cover a range of projected distances 0.22–3.57 h⁻¹ cMpc.

Our empirical estimate for w_p(R) obtained via the above methods is displayed in Figure 6. We find a best-fit correlation length r₀ = (7.3 ± 1.3) h⁻¹ Mpc after fixing γ = 1.5, where the error is a 1σ uncertainty computed via a bootstrap estimate. This procedure consisted of repeating the entire analysis 100 times (computation of selection functions, counting of pairs, and χ²-fitting of w_p) using a random bootstrap sample of 15 of the 15 independent fields selected with replacement. The quoted uncertainty is the standard deviation of parameter values derived from these 100 bootstrap samples. The results of this procedure were consistent with the results of jackknifing estimates performed using 14 of the 15 fields (i.e., an n − 1 jackknife estimate) or using 8 of the 15 fields (i.e., an approximately ∼n/2 jackknife estimate). The χ² value for the fit is 3.7 on 4 degrees of freedom.

Zoom In Zoom Out Reset image size

As noted in Section 1, Adelberger & Steidel (2005) performed a cross-correlation measurement with a similar sample of color-selected galaxies to compare BH and galaxy masses over a large range of AGN luminosities (−20 ≳ M_AB(1350 Å) ≳ −30) at a similar range of redshifts to our galaxy sample (1.5 ≲ z ≲ 3.6). That study separated the AGN sample into two bins of BH mass, obtaining galaxy–AGN cross-correlation lengths r₀ = 5.27^+1.59_−1.36 h⁻¹ cMpc for AGNs with 10^5.8 < M_BH/M_☉ < 10⁸ and r₀ = 5.20^+1.85_−1.16 h⁻¹ cMpc for AGNs with 10⁸ < M_BH/M_☉ < 10^10.5. These measurements are fairly consistent with our own measurement of r₀ for the galaxy–HLQSO cross-correlation, given the size of the uncertainties, and Adelberger & Steidel (2005) assume a correlation function slope of 1.6, rather than the slope of 1.5 used in this study. As noted above, the best-fit value of r₀ varies inversely with the chosen slope for the range of separations our measurements include, and this effect likely accounts for the slight discrepancy between these two estimates.

The strength of the clustering signal corresponds to the mass scale of the HLQSO host halos, which we are interested in comparing to the average halo mass scale of non-active galaxies. As such, the relative strengths of the galaxy–galaxy (GG) and galaxy–HLQSO (GQ) clustering reveal the relative mass scales of their respective host halos, and therefore illuminate any halo-mass requirements for the formation of HLQSOs at z ≃ 2.7.

Our estimate of the galaxy–galaxy correlation function is based on the same technique as our galaxy–HLQSO estimates but is modified by centering on each galaxy, rather than the HLQSO, in turn. In addition, we restrict our GG analysis to those galaxies at redshifts z_gal > 2.25 so that the GG autocorrelation function probes a similar redshift range to that of the GQ cross-correlation; the 909 galaxies with z > 2.25 have a median redshift z^med_gal = 2.63, while the median HLQSO redshift is z^med_QSO = 2.66. For each galaxy in our sample, we consider the number density (per unit solid angle) of galaxies as a function of projected distance from our fiducial galaxy, separating between redshift-associated galaxies (those within 1500 km s⁻¹ and in the same field as the fiducial galaxy) and non-associated galaxies (those outside the velocity range or in a different field). We then integrate over the redshift selection function (Section 3.1) to find the expected number of galaxies in each interval, from which we can define an angular correlation function for each interval (by analogy to Equation (11)):

$$\begin{equation} \frac{\Sigma _{ v,{\rm obs}}^{{\rm GG}}(R)}{\Sigma _{ v,{\rm pred}}^{{\rm GG}}} = \frac{\Sigma _{ 0,{\rm obs}}^{{\rm GG}}(R)}{\Sigma _{0,{\rm pred}}^{{\rm GG}}}\big[1+w_{{p}}^{{\rm GG}} (R)\big], \end{equation} \tag{ 14 }$$

where again the v subscript denotes quantities corresponding to the redshift-associated sample, and the 0 subscript denotes those corresponding to the non-associated sample. The non-associated sample is not expected to cluster about the arbitrary line of sight defined by the position of our fiducial galaxy, so we interpret the quantity Σ_{obs, 0}(R)/Σ_{pred, 0}(R) as an estimate of the relative completeness of our angular sampling. We found the completeness (i.e., the angular selection function of the galaxy–galaxy pairs) of our sample to be well described by a linear model in R with negative slope: Σ₀ = aR + b. This shape reflects the fact that we are able to measure the power on small scales for essentially all galaxies, while we can see the maximum separation ≃ 2R_max only for the small fraction of galaxies at the edge of our fields (and even then we see only the subset of pairs that lie entirely within the field). The best-fit parameter values for the GG angular selection function are a = −0.160, b = 1.04 with R in h⁻¹ cMpc (Figure 5).

As in the case of the GQ cross-correlation function, we then fit a model to Σ^GG_v that is a combination of the underlying clustering signal described by w^GG_p(R) and the selection function described by Σ^GG₀. The combined model is given by Equation (15):

$$\begin{eqnarray} &&\frac{\Sigma _{v}^{{\rm GG}}(R)}{\Sigma _v^{{\rm pred}}} = (a R+b)\left[1+\left(\frac{r_0^{{\rm GG}}}{R}\right)^\gamma {}_2F_1\left(\frac{1}{2},\frac{\gamma }{2};\frac{3}{2};\frac{-z_0^2}{R^2}\right)\right].\nonumber\\ \end{eqnarray} \tag{ 15 }$$

Fitting this model to the measured values of Σ^GG_v(R_k), we find the best-fit galaxy autocorrelation length to be r^GG₀ = (6.0 ± 0.5) h⁻¹ Mpc, again fixing the slope γ = 1.5. In this case, the errors in Σ^GG_v cannot be considered Poissonian because each galaxy is counted in several pairs and the galaxy counts are correlated between bins. However, the 15 HLQSO fields each provide an independent estimate of Σ^GG_v, and the error used for the χ²-fitting is based on the scaled scatter among these values. The quoted error on r^GG₀ is the 1σ uncertainty from the same bootstrap and jackknife procedures described for r^GQ₀ (Section 4.1.1). The χ² value for the fit is 9.3 on 5 degrees of freedom.

This autocorrelation length is significantly larger than that found by Adelberger et al. (2005b) (r₀ = (4.0 ± 0.6) h⁻¹ cMpc at z = 2.9), despite both studies relying on a similarly selected set of galaxies. However, this study is restricted to the spectroscopically observed galaxies, which have a higher mean luminosity than the galaxies in the photometric sample used in Adelberger et al. (2005b), and are more comparable to the higher-luminosity subsample of galaxies used in that paper, for which the authors estimated r₀ = (5.2 ± 0.6) h⁻¹ cMpc. In addition, the much larger set of spectroscopic redshifts used in our sample allows us to characterize the redshift selection function with much greater accuracy, as well as to restrict our analysis to those galaxies associated with the HLQSOs in three-dimensional space. For example, an injudicious choice of z₀ in Equation (10) would lower the estimated cross-correlation length, either by failing to count HLQSO-associated galaxies (z₀ < 1500 km s⁻¹) or by diluting the clustering signal by the inclusion of the voids adjacent to the HLQSO in redshift space (z₀ > 1500 km s⁻¹).

The measured clustering of the galaxies in our sample is primarily useful in its connection to the mass scale of the galaxy host halos. Because the clustering strength of DM halos is a function of halo mass, we can invert this relation to obtain the halo mass for a population of objects with a given autocorrelation length. In practice, we perform this inversion numerically, finding the population of simulated halos (for which the mass is known) that match the clustering strength of our galaxy sample.

Using halo catalogs from the MultiDark MDR1 simulation (Prada et al. 2012; accessed via the MultiDark Database of Riebe et al. 2012), we measured the correlation length r₀ as a function of minimum halo mass M_h using the Landy–Szalay estimator (Landy & Szalay 1993) and assuming the same power-law slope (γ = 1.5) used in our fit to the galaxy autocorrelation function. The correlation function was measured for halo populations of differing masses by varying the minimum M_h in steps of 0.05 dex; the correlation lengths for a subset of the halo samples are listed in Table 3. A power-law correlation function is a poor fit to the halo clustering at small scales due to the effect of halo exclusion; therefore, we restricted our fit to pairs with separations 1 ⩽ d/(h⁻¹ Mpc) ⩽ 5, a range that avoids the halo-exclusion zone while still closely matching the range of projected distances in our observed sample. In this manner, we find that our galaxy sample is most consistent with having a minimum halo mass log (M_h/M_☉) > 11.7 ± 0.1, the fit to which is displayed in Figure 6. The halos in this mass range have a median mass log (M_h/M_☉) = 11.9  ±  0.1. The statistical error in the mass estimate is entirely due to the propagated error in the autocorrelation function, as the uncertainty in the autocorrelation function among simulated halos is negligible by comparison.

In addition to matching clustering strengths, we can also attempt to match the abundances of observed galaxies and simulated halos. Although our spectroscopic sample of galaxies is incomplete, we can compare this halo population to the galaxy luminosity function (GLF) of Reddy et al. (2008), which corrects for incompleteness in both the spectroscopic and photometric samples. Luminosity functions are measured separately for galaxies with 1.9 ⩽ z < 2.7 and 2.7 ⩽ z < 3.4, while our sample straddles these two redshift intervals, but the GLF evolves very little over this redshift range, and the predictions of either model are quite similar. Using the Schechter (1976) GLF parameters listed in Table 7 of Reddy et al. (2008), and taking our magnitude limit $\mathcal {R}<25.5$ to correspond to M_AB(1700 Å) ≲ −19.9 at z ≃ 2.7, the Reddy et al. models predict a galaxy number density ϕ_gal = (2.4–7.0) × 10⁻³ h³ Mpc⁻³ (including the 1σ limits on ϕ*). The number density of log (M_h/M_☉) > 11.7 halos in the MultiDark MDR1 simulation is ϕ_sim = 4.4 × 10⁻³ h³ Mpc⁻³, entirely consistent with the measured value of ϕ_gal.

Taking the population of log (M_h/M_☉) > 11.7 halos to represent the host halos of the galaxies in our sample, we then estimate the mass of the HLQSO hosts by finding the population of simulated halos whose cross-correlation with the representative galaxy halos is equal to our measured galaxy–HLQSO cross-correlation function. Again varying the minimum M_h (of fiducial HLQSO hosts) in 0.05 dex increments, and using a Landy–Szalay variant for a cross-correlation with γ = 1.5, we find that an HLQSO host halo mass log (M_h/M_☉) > 12.1 ± 0.5 (a median halo mass of log (N_h/M_☉) = 12.3 ± 0.5) is most consistent with our galaxy–HLQSO cross-correlation measurement. The error in the mass is due to the error on both r^GQ₀ and the propagated error on the galaxy host halo mass, since the strength of the cross-correlation function depends on the mass of both the HLQSO host and galaxy host halo populations. The fit to the corresponding simulated cross-correlation function is shown in Figure 6. The MultiDark halos of log (M_h/M_☉) > 12.1 ± 0.5 have an autocorrelation length of 6–15 h⁻¹ cMpc, which we consider to be an estimate of the HLQSO autocorrelation length, and such halos have an abundance ϕ_sim = (0.17–5.9) × 10⁻³ h³ Mpc⁻³ at z ∼ 2.5 in the simulation.

The MultiDark suite of simulations used cosmological parameters based on the WMAP 5 year results, {Ω_m,Ω_Λ,σ₈,h} = {0.27,0.73,0.82,0.70}, which are consistent with the most recent WMAP 7 year results from Larson et al. (2011): {0.276 ± 0.029,0.734 ± 0.029,0.801 ± 0.030,0.710 ± 0.025}. In comparison, the older (and widely utilized) Millennium Simulation (Springel et al. 2005) used cosmological parameters based on the WMAP 1 year results, {Ω_m,Ω_Λ,σ₈,h} = {0.25,0.75,0.9,0.73}. We here consider how such a variation in these cosmological parameters affects the halo-matching process employed in this study.

The parameters Ω_m and σ₈ both affect halo abundances, and thus affect the halo bias and the mapping from clustering strengths to halo masses. Zehavi et al. (2011) conduct an HOD analysis on a large sample of galaxies and find that varying the matter density over the range 0.25 ⩽ Ω_m ⩽ 0.3 produces only a ∼2% variation in their clustering measurements, which is quite small compared to the statistical uncertainty in our measurements.

However, the amplitude of the linear DM fluctuations, σ₈, is tied to the clustering in a more pronounced and complicated manner. The clustering of galaxies in linear theory is given by the galaxy bias and the DM clustering: ξ_GG(M) = b²(M)ξ_DM. Decreasing σ₈ decreases the value of ξ_DM but also greatly decreases the number density of high-mass halos, which causes the bias at a given halo mass b(M) to increase. For high-mass halos, the overall effect is to increase the clustering strength at a given mass when σ₈ is decreased, which suggests that mapping halo masses to clustering strengths using the Millennium Simulation would result in a shift toward larger halo masses. Repeating our halo-matching analysis on Millennium halo catalogs, we find that the best-matched halo population of galaxies has a minimum halo mass log (M_{h, Mill}/M_☉) > 12.0 ± 0.1 and a median mass log (M_{h, Mill}/M_☉) = 12.2 ± 0.1 in that simulation. The discrepancy between the two simulations is ∼3× the statistical uncertainty in the galaxy halo-mass measurements, confirming that the clustering of these massive halos is quite sensitive to the chosen cosmological parameters.

The mass scales of the host halos for galaxies and HLQSOs map to halo abundances, as described in Section 4.3. A galaxy host halo abundance of ϕ_{sim, 11.7} = 4.4 × 10⁻³ h³ Mpc⁻³ (the MultiDark simulation abundance of halos with log (M_h/M_☉) > 11.7) and an HLQSO host halo abundance of ϕ_{sim, 12.1} = 1.2 × 10⁻³ h³ Mpc⁻³ suggests that halos massive enough to host an HLQSO are only ∼4× less abundant than those massive enough to host the average galaxy in our sample; the fact that far fewer than one quarter of the galaxies in our sample host an HLQSO is a strong constraint on the duty cycle of these objects. However, the precise value of the HLQSO duty cycle depends on the number density of HLQSOs, which in turn depends on the choice of QSO population.

All of our HLQSOs have luminosities at rest-frame 1450 Å  of log (νL_ν/L_☉) ∼ 14, or an absolute magnitude M(1450 Å) ∼ -30.² This is brighter than the luminosity range for which large-sample statistics are available in surveys such as SDSS and SLAQ (see, e.g., Croom et al. 2009, whose M_g(z = 2) is comparable to M(1450 Å)), but we can obtain an estimate of the z ∼ 2.7 quasar luminosity function by extrapolating the results of the highest-redshift bins of Croom et al. (2009) to slightly higher redshifts and luminosities; in this way we roughly estimate the number density of M(1450 Å) ≳ −30 QSOs to be ϕ_QSO ∼ 10^−9.5 h³ Mpc⁻³. Integrating this density over the total comoving volume between redshifts 2 ⩽ z ⩽ 3 predicts ∼25 QSOs in this luminosity range over the entire sky, suggesting that a large fraction of the comparably bright QSOs at these redshifts are already in our sample.

Given this number density, we can extract the duty cycle of HLQSOs from the ratio ϕ_QSO/ϕ_{sim, 12.1} ≃ 10⁻⁶ to 10⁻⁷, defining the duty cycle as the fraction of halos massive enough to host an HLQSO [log (L/L_☉) ≳ 14] that actually do host such a QSO. This extreme rarity with respect to the number of potential host halos indicates that the formation of the HLQSO must rely on a correspondingly rare event occurring on scales much smaller than those probed by our analysis, perhaps related to an extremely atypical merger or galaxy interaction scenario.

It is interesting to compare the host halo masses of the HLQSOs to the minimum BH masses allowed by their luminosities under the assumption of Eddington-limited accretion; we will refer to this minimum mass as M_BH. The minimum BH masses for each HLQSO are listed in Table 2. We calculate the value of M_BH directly from L₁₄₅₀ (the value of νL_ν at a rest-frame wavelength of 1450 Å):

$$\begin{equation} M_{{\rm BH}} = \frac{\sigma _{T} L_{{\rm 1450}}}{4 \pi G m_{{\rm p}} c} = 3.1 \times 10^{-5} \left(\frac{L_{{\rm 1450}}}{L_{\odot }}\right) M_{\odot }\,. \end{equation} \tag{ 16 }$$

We use L₁₄₅₀ in place of the bolometric luminosity L_bol in order to avoid the additional uncertainty in the bolometric correction. The true bolometric correction is likely to be small: Nemmen & Brotherton (2010) use a thin accretion disk model to predict a correction factor L_bol/L₁₄₅₀ ∼ 3, which is consistent with the empirical correction estimated by Netzer & Trakhtenbrot (2007) for L₅₁₀₀ and scaled by the L₅₁₀₀/L₁₄₅₀ relationship of Netzer et al. (2007). However, there is substantial scatter in these corrections, and it may be expected that L_bol/L₁₄₅₀ approaches unity for QSOs selected by the most extreme rest-UV luminosities, so L₁₄₅₀ is a useful estimate (and likely a lower limit) on L_bol.

The HLQSOs in our sample span a range of ∼5× in M_BH (with the possible exception of Q0142−10), and it is notable that the field (HS1549+1919) with the largest value of M_BH is also associated with the largest redshift overdensity in the galaxy distribution (see column N₁₅₀₀ in Table 2). However, there is no clear relation between N₁₅₀₀ and M_BH among the other fields, and our galaxy samples are not large enough to comment on the variation of M_BH with halo mass.

The median value of M_BH for our sample is log (M_BH/M_☉) ≃ 9.7. This indicates that the HLQSO host DM halos are only ∼300–2000× more massive than their associated supermassive BHs, even assuming accretion at the Eddington limit. The relationship between BH mass and DM halo mass is uncertain even at z ≃ 0 (compare, e.g., Ferrarese 2002, Booth & Schaye 2010, and Kormendy & Bender 2011), but these BHs lie well above the predictions of the M_BH–σ or M_BH–v_c relations for any reasonable mapping of the DM halo mass to the bulge velocity dispersion σ or circular velocity v_c, as demonstrated below.

Three such mappings are considered in Ferrarese (2002), in which the halo virial velocity v_vir is related to the circular velocity by considering v_c = v_vir (a zeroth-order approximation), v_c = 1.8v_vir (based on observational constraints on DM halo-mass profiles by Seljak 2002), or v_c/v_vir given by a function of halo mass extracted from the N-body simulations of Bullock et al. (2001). These three different assumptions predict BH masses of log (M_BH/M_☉) = 7.0, 8.4, and 7.5, respectively, for a halo of mass log (M_h/M_☉) = 12.3, corresponding to M_DM/M_BH ≃ 2 × 10⁵, 8 × 10³, and 6 × 10⁴. In any of these cases, the minimum BH masses for the HLQSOs in our sample are 1–2 orders of magnitude higher than the predictions of the low-redshift associations, implying that the host halos must "catch up" with the BHs in order to fall on the established relations by z ≃ 0.

Though estimates of BH masses at high redshift are highly uncertain, as are the stellar masses of their host galaxies, this result agrees qualitatively with several observational studies that find BH host galaxies at high redshift (1 ≲ z ≲ 4) of a given stellar mass have systematically higher BH masses than in the local universe (e.g., Peng et al. 2006; Decarli et al. 2010; Merloni et al. 2010; Greene et al. 2010); Booth & Schaye (2010, 2011) describe an interpretation of this evolution in terms of the compactness of DM halos, which are more tightly bound at high redshifts. These studies generally find a smaller deviation from the z ≃ 0 relations than is present in our sample, but the extreme luminosities of the HLQSOs in our sample force us to probe the highest-mass end of the BH mass distribution, so our measurements are very sensitive to the scatter in the M_h–M_BH relation as well as evolution in the mean.

If we consider the possibility that the dynamical mass discussed in Section 5 includes matter external to the HLQSO host halo at z ≃ 2.7, which may merge into a single more massive halo by z ≃ 0, we can calculate where such a halo would fall in the M_h–M_BH relations of Ferrarese (2002). Taking a halo mass log (M_h/M_☉) = 13, the above prescriptions predict BH masses of log (M_BH/M_☉) = 8.3, 9.6, and 8.7, respectively, which lie much closer to the range of BH masses seen in our sample. However, it seems clear that the extremely high BH masses indicate that the HLQSOs are atypical (with respect to the general population of QSOs) at the smallest scales.