The characteristic halo masses of half-a-million WISE-selected quasars

TL;DR

This study refines measurements of dark matter halo masses for WISE-selected quasars, confirming that obscured quasars inhabit significantly more massive haloes than unobscured ones, supporting an environmental or evolutionary link to obscuration.

Contribution

It expands the survey area and improves measurement precision, confirming the halo mass difference between obscured and unobscured quasars using multiple techniques.

Findings

Obscured quasars reside in haloes three times more massive than unobscured quasars.

The halo mass difference is statistically significant at about 5 sigma.

Obscured and unobscured quasars are well-matched in redshift and luminosity, supporting intrinsic differences.

Abstract

Recent work has found evidence for a difference in the bias and dark matter halo masses of WISE-selected obscured and unobscured quasars, implying a distinction between these populations beyond random line-of-sight effects. However, the significance of this difference in the most up-to-date measurements is relatively weak, at $\sim$2$\sigma$ for individual measurements but bolstered by agreement from different techniques, including angular clustering and cross-correlations with cosmic microwave background (CMB) lensing maps. Here, we expand the footprint of previous work, aiming to improve the precision of both methods. In this larger area we correct for position dependent selection effects, in particular fluctuations of the WISE-selected quasar density as a function of Galactic latitude. We also measure the cross-correlation of the obscured and unobscured samples and confirm that they are well-matched in redshift, both centred at $z=1$. Combined with very similar detection fractions and magnitude distributions in the long-wavelength WISE bands, this redshift match strongly supports the fact that IR selection identifies obscured and unobscured quasars of similar bolometric luminosity. Finally, we perform cross-correlations with confirmed spectroscopic quasars, again confirming the results from other methods --- obscured quasars reside in haloes a factor of 3 times more massive than unobscured quasars. This difference is significant at the $\sim$5$\sigma$ level when the measurements are combined, strong support for the idea that obscuration in at least some quasars is tied to the larger environment, and may have an evolutionary component.

Figures

Tables

Table 1. Table 1: Summary of samples and bias and halo mass measurements.

Sample		N	$⟨z⟩$	$b_q$	$\log \left(M_h/M_{\odot }{h}^{-1}\right)$
		Angular Clustering
All		581,728	1.02	1.98$\pm$0.05	${12.79}_{-0.04}^{+0.04}$
All (SDSS footprint)		296,793	1.02	1.94$\pm$0.07	${12.75}_{-0.06}^{+0.06}$
Unobscured		172,537	1.05	1.77$\pm$0.11	${12.56}_{-0.12}^{+0.11}$
Obscured		121,486	0.98	2.15$\pm$0.12	${12.96}_{-0.09}^{+0.08}$
		CMB Lensing Cross-correlation
All		531,545	1.02	1.82$\pm$0.14	${12.65}_{-0.15}^{+0.13}$
All (SDSS)		269,751	1.02	1.85$\pm$0.19	${12.68}_{-0.19}^{+0.16}$
Unobscured		157,172	1.05	1.66$\pm$0.23	${12.45}_{-0.31}^{+0.23}$
Obscured		94,938	0.98	2.14$\pm$0.27	${12.95}_{-0.22}^{+0.18}$
		Combineda
Unobscured			1.05	1.70$\pm$0.10	${12.49}_{-0.11}^{+0.10}$
Obscured			0.98	2.15$\pm$0.11	${12.94}_{-0.08}^{+0.08}$
a Weighted mean of clustering and lensing results, after verifying sufficient independence; see Section 6.

Table 2. Table 2: Summary of cross-correlation results

		$b_1{b}_2$	$b_{\text{unob}}$	$b_{\text{obsc}}$
Obsc $\times$ Unob		4.12$\pm$0.41	1.92$\pm$0.22	2.33$\pm$0.27
Spec $\times$ Unob		3.06$\pm$1.26	1.1$\pm$0.5	…
Spec $\times$ Obsc		5.29$\pm$1.38	…	1.9$\pm$0.6

Equations

$$dP=n(1+{\omega }_{\theta })d\Omega .$$

$$\begin{gathered} {\omega }_{dm}(\theta )=\pi {\int }_{z=0}^{\infty }{\int }_{k=0}^{\infty }\frac{{\Delta }^2(k,z)}{k}{J}_0[k\theta \chi (z)]\times \\ {\left(\frac{dN}{dz}\right)}_1{\left(\frac{dN}{dz}\right)}_2\left(\frac{dz}{d\chi }\right)\frac{dk}{k}dz. \end{gathered}$$

$${\omega }_{\text{qq}}(\theta )=\frac{DD-2DR+RR}{RR},$$

$$\begin{gathered} {\text{C}}_{ij}=\sum _{L=1}^N\sqrt{\frac{R{R}_L({\theta }_i)}{RR({\theta }_i)}}[{\omega }_L({\theta }_i)-\omega ({\theta }_i)]\times \\ \sqrt{\frac{R{R}_L({\theta }_j)}{RR({\theta }_j)}}[{\omega }_L({\theta }_j)-\omega ({\theta }_j)], \end{gathered}$$

$${\sigma }_p^2(\theta )=\frac{f(1+\omega (\theta ){)}^2}{DD}$$

$$b_{\text{unob}}=\frac{b_{\text{cross}}^2({\mu }_2)}{b_{\text{obsc}}({\mu }_2)}.$$

$$W^q(\chi )=\frac{dz}{d\chi }\frac{dN}{dz}{b}_q(\chi ),$$

$$W^{\kappa }(\chi )=\frac{3}{2}{\Omega }_{\mathrm{m}}{\left(\frac{H_0}{c}\right)}^2\frac{\chi }{a(\chi )}\frac{{\chi }_{\mathrm{CMB}}-\chi }{{\chi }_{\mathrm{CMB}}},$$

$$C_l^{\kappa q}=\int dz\frac{d\chi }{dz}\frac{1}{{\chi }^2}{W}^{\kappa }(\chi )W^q(\chi )P\left(\frac{l}{\chi },z\right),$$

$${\delta }_q=\frac{\rho -\bar{\rho }}{\bar{\rho }},$$

$${\text{C}}_{ij}=\frac{1}{N-1}\left[\sum _{k=1}^N(C_{l_i,k}^{\kappa q}-C_{l_i}^{\kappa q})(C_{l_j,k}^{\kappa q}-C_{l_j}^{\kappa q})\right],$$

$$b(\nu )=1-A\frac{{\nu }^a}{{\nu }^a+{\delta }_c^a}+B{\nu }^b+C{\nu }^c,$$

$${\delta }_c=0.15(12\pi {)}^{2/3}{\Omega }_{m,z}^{0.0055},$$

$${\sigma }^2(M)=\frac{1}{2{\pi }^2}\int P(k,z){\hat{W}}^2(k,R_h)k^{2}dk.$$

$$\hat{W}(k,R_h)=\frac{3}{(k{R}_h{)}^3}[\sin (k{R}_h)-k{R}_h\cos (k{R}_h)].$$

Full text

††pagerange: The characteristic halo masses of half-a-million WISE-selected quasars–References††pubyear: 2015

1 INTRODUCTION

With the development of simple infrared (IR) quasar selection criteria (e.g. Lacy et al., 2004; Stern et al., 2005, 2012; Mateos et al., 2012; Assef et al., 2013), large samples can now be assembled across the entire sky from the Wide-field Infrared Survey Explorer (WISE; Wright et al., 2010). Critically, selecting quasars in the mid-IR identifies obscured quasars, which have been previously poorly represented in statistical samples, while still also selecting the much better studied unobscured population. While the bolometric luminosities of the two classes selected this way are quite similar (Hickox et al., 2007), large obscuring columns of dust make obscured quasars difficult to identify photometrically in optical surveys, such as the Sloan Digital Sky Survey (SDSS; York et al., 2000), despite comprising roughly half of the quasar population (Assef et al., 2013). The ability of mid-IR radiation to penetrate this dust provides a unique tool for exploring these “hidden monsters.”

With large populations of both obscured and unobscured quasars selected in a similar manner, we can now deeply explore questions regarding the nature of obscured quasars in the context of both AGN unification models (e.g. Antonucci, 1993) and evolutionary models of quasars and galaxies (e.g. Hopkins et al., 2008; Croton, 2009). In particular, the location and distribution of the obscuring material in the most luminous quasars, such as those selected with simple colour-cuts from relatively shallow WISE data, is still a matter of debate. The simplest unification models ascribe all of the obscuration to random orientations and an axis-symmetric structure just beyond the accretion disk commonly known as the “dusty torus”, the properties of which may depend on other factors such as luminosity (e.g. Lawrence, 1991; Netzer, 2015). This explanation implies that the bulk properties of the unobscured and obscured populations on larger scales should be the same, including those of their host galaxies, dark matter haloes, and environment more generally.

However, obscuration can also occur due to dust in the host galaxy (e.g. Page et al., 2004; Goulding et al., 2012; Buchner & Bauer, 2017). While this can also be due to the randomness of our line of sight, models that invoke galaxy mergers as a significant fueling mechanism for black hole growth suggest that large-scale, high covering fraction obscuration can occur naturally in an evolutionary quasar sequence (e.g. Sanders et al., 1988; Hopkins et al., 2008). In such a scenario, where obscured quasars are a phase in the overall quasar life cycle, differences in the host galaxies (e.g. Chen et al., 2015) and larger scale environments may be present (e.g. Ellison, Patton & Hickox, 2015)

Over the last several years, much work has been done to test these models by measuring the characteristic dark matter halo masses of IR-selected obscured and unobscured quasars (see DiPompeo, Hickox & Myers, 2016, hereafter D16, for a detailed review of the IR-selection of quasars and the use of such samples to probe large-scale structure). While the results have varied depending on the methodology and precise sample selection, most have agreed that IR-selected obscured quasars reside in haloes at least a factor of a few times more massive than their unobscured counterparts (Hickox et al., 2011; Donoso et al., 2014; DiPompeo et al., 2014, 2015b; DiPompeo, Hickox & Myers, 2016). Using simple relationships between halo mass, galaxy mass, and black hole mass, DiPompeo et al. (2017) showed that a simple evolutionary trend from obscured to unobscured quasar, with black hole growth lagging behind that of the rest of the system, could broadly reproduce the measured halo mass difference, as well as explain some discrepancies in the literature (i.e. the similar IR-selected obscured/unobscured quasar halo masses at lower redshift and luminosity; Mendez et al., 2016).

Though the halo mass difference was still present in the most recent measurements of D16, the significance was reduced to $\sim$2$\sigma$ after employing the latest photometric samples and CMB lensing maps. While this difference is likely a lower limit on the potentially evolutionarily distinct obscured population (due to a subset being unified by orientation; DiPompeo et al., 2016), the measurements thus far have relied on a portion of the total available area, where systematics in the survey data are less prevalent. Given the important implications of any potential halo mass difference between the populations and the current state of the measurements, it is prudent to expand the footprint, increasing the sample sizes to their full potential for the most precise possible measurement. That, along with exploring cross-correlations, is our goal in this work.

All models and measured properties use a cosmology of $H_{0}=70.2$ km s*-1* Mpc*-1*, $\Omega_{\textrm{M}}=\Omega_{\textrm{CDM}}+\Omega_{\textrm{b}}=0.229+0.046=0.275$, $\Omega_{\Lambda}=0.725$, and $\sigma_{8}=0.82$ (Komatsu et al., 2011). Magnitudes are given in the Vega system unless otherwise specified.

2 Quasar Sample

2.1 WISE Quasar Selection

We make our initial sample selection using both the original ALLSKY and the updated ALLWISE catalogs from WISE. The WISE mission mapped the entire sky in four bands: 3.4, 4.6, 12, and 22 $\mu$m, designated as $W1$, $W2$, $W3$, and $W4$, respectively. The two shortest wavelength bands are the most sensitive and highest resolution ($\sim$0.1 mJy 5$\sigma$ sensitivity compared to $\sim$1 mJy, and $\sim$6 arcsec versus $\sim$10 arcsec resolution), and are the most efficient for selecting luminous quasars (Stern et al., 2012). The ALLSKY catalog was the first full source catalogue release in March 2012, and was supplanted by the ALLWISE release in November 2013, which included new imaging in the $W1$ and $W2$ bands and an improved data calibration pipeline.

D16 provided a detailed analysis of quasar samples selected from these two catalogs in the context of the measurements performed here (angular clustering in Section 3 and CMB lensing cross-correlations in Section 4). While they found that the ALLWISE survey provides a slightly brighter and bluer but less-contaminated sample, there were potentially problematic systematic correlations in the source distribution with Moon contamination that were not present in the ALLSKY catalog. In order to provide the benefits of both and a conservatively generated sample, we adopt their approach and select objects that satisfy the Stern et al. (2012) selection criteria of $W1-W2>0.8$ and $W2<15.05$ in both catalogs.

After selecting objects from the catalog data, we apply the same “mask” to the data as D16, which discards regions of likely contamination by defining holes in the survey footprint. We refer the reader to D16 for complete details, but in general we remove the Galactic plane within $|b|\leq 20$ degrees, regions of high Galactic extinction where $A_{g}>0.18$, regions of high Moon contamination, regions around clustered data flagged as spurious in the WISE catalogs (using cc_flags; see DiPompeo et al., 2014), regions of poor WISE photometric quality (based on the ph_qual flag), and areas around bright stars as catalogued by SDSS. These masks are generated for both the ALLSKY and ALLWISE catalogs and both are applied — we provide mangle (Swanson et al., 2008) polygon files of the mask at http://faraday.uwyo.edu/~admyers/wisemask2017/wisemask.html.

After applying the mask, we are left with 581,728 WISE-selected quasars over 12,745 deg2. The distribution of this sample on the sky is shown in red in Figure 1.

2.2 Obscured and Unobscured Quasars

While the IR data alone can select quasars, additional optical information is needed to divide the sample into obscured and unobscured populations. Here we follow Hickox et al. (2007) and exploit the bimodal distribution of quasar optical-IR colours, using $r-W2=6$ to divide the samples111We stress that this definition is different from the more traditional type 1 (unobscured) and type 2 (obscured) quasar definitions, which are based on the presence or absence of broad emission lines in spectra. (see also Hickox et al. 2017, in prep).

We gather optical $r$-band magnitudes from data release 8 of the SDSS (DR8222Note that the astrometric error in SDSS DR8 that was subsequently corrected in DR9 (Ahn et al., 2012) does not impact our matching, given the comparatively low resolution of the WISE imaging. All of our objects are matched to the same optical counterpart in DR8 and DR9., adopting psf_mags though the use of model_mags does not impact the results strongly; DiPompeo et al., 2014), matching to the ALLWISE catalog positions with a radius of 2 arcseconds. We limit the optical footprint to that of the BOSS survey333http://www.sdss3.org/dr9/algorithms/boss_tiling.php#footprint, which covers 10,269 deg2. The imaging in many of the ancillary regions outside of the BOSS footprint is less uniform, which complicates the splitting of the sample and can lead to problems with the relative obscured and unobscured density as a function of position. We include mangle polygon files of the BOSS footprint, as well as a pixelated version of the DR8 footprint (which is also used in order to remove small regions where imaging is not complete) at online with the other mask files above. We also apply the SDSS “bad fields” mask (e.g. White et al., 2011), which is included online with the other masking files for completeness.

Within the available SDSS$+$WISE footprint (5314 deg2), there are 294,023 WISE-selected quasars, overlaid in blue on the full WISE-selected sample in red in Figure 1. Of these, 245,800 (83.5%) have optical matches with $15<r<25$. The $r-W2$ colours are calculated using the ALLWISE $W2$ magnitudes. Sources without optical detections are considered obscured, and the final samples consist of 172,537 (58.6%) unobscured and 121,486 (41.4%) obscured quasars, roughly a factor of two larger than the samples in D16. Due to more strict masking (Section 4), these numbers are reduced somewhat for the CMB lensing correlations. We will consider the impacts of including or neglecting sources without optical matches in later sections.

2.3 Redshift distributions and luminosities

In order for fair comparisons between the unobscured and obscured populations, they must be well-matched in redshift and luminosity. However, without individual source redshifts, we are limited to characterizing the overall properties of the sample more generally by studying subsets in small fields with deep multi-band imaging and follow up spectroscopy. As in D16, we apply our selection criteria and mask to sources in the Boötes field (Brodwin et al., 2006; Hickox et al., 2007), and find 368 WISE-selected quasars in the $\sim$7 deg2 portion with deep spectroscopy (Kochanek et al., 2012). Nearly all of these sources (361, or 98%) have spectroscopic redshifts, which allows us to generate accurate redshift distributions ($dN/dz$) for both the obscured and unobscured samples, shown in Figure 2. Both samples cover a similar range ($0<z<3$) and have mean redshifts of $\sim$1 (0.98 and 1.05 for obscured and unobscured, respectively). We will utilize cross-correlations between samples to further probe the similarity in redshift distributions in Section 3.2.3.

Because there are trends between black hole mass, galaxy mass, and halo mass (e.g. Häring & Rix, 2004; Booth & Schaye, 2010; Guo et al., 2010; Kormendy & Ho, 2013, and references therein), it is expected that there should also be correlations between quasar luminosity and halo mass. However, observational evidence of this relationship has been elusive (Shen et al., 2009; Eftekharzadeh et al., 2015), potentially due to the breadth of the Eddington ratio distribution and variability timescales (Hopkins & Hernquist, 2009; Hickox et al., 2014; Jones et al., 2016). Regardless, because of this potential link, comparisons of obscured and unobscured quasar halo masses should only be made for similar luminosities. Hickox et al. (2007, 2011) used detailed SED fitting and bolometric corrections to the Spitzer rest-frame 8 $\mu$m luminosities of IR-selected quasars in the Boötes field to illustrate that simple IR colour-cuts select high luminosity quasars (peaking at $L_{\textrm{bol}}\sim 10^{46}$ erg/s), and that the obscured and unobscured populations selected using optical-IR colours have very similar luminosity distributions.

An additional complication in calculating intrinsic luminosities for WISE obscured quasars is that in the most sensitive (shorter wavelength) bands, the IR emission may still be significantly anisotropic (e.g. Hönig et al., 2011), as well as affected by reddening in sources with the highest column densities. A simple check on the importance of these effects in our samples is to look at the detection fractions and magnitude distributions in the long-wavelength WISE bands, as well as colours across the WISE filters.

By definition, all of our targets have $W1$ and $W2$ measurements. In the $W3$ and $W4$ bands, we find obscured detection fractions ($>5\sigma$) of 96% and 64%, respectively, while for the unobscured sample we find detection fractions of 97%, and 62%, respectively. In the top panel of Figure 3 we show the colour distributions of each sample at shorter ($W1-W2$) and longer ($W3-W4$) wavelengths. It is clear that at shorter wavelengths reddening (and potentially orientation) effects are present, as the obscured sample is redder on average due to a significant red tail in the distribution. However, this difference completely disappears at the longer wavelengths, highlighting the utility of these bands for direct comparisons.

In the second panel of Figure 3 we show the magnitude distributions in $W3$ and $W4$ of the detected sources, and again find excellent agreement between the samples. The median magnitudes in the detected obscured sources are 11.3 ($W3$) and 8.5 ($W4$), and 11.4 and 8.5 for the unobscured sources. Considering the well-matched redshift distributions, these similarities suggest excellent correspondence in luminosity.

3 ANGULAR CLUSTERING

3.1 Methodology

Without individual source redshifts, we are limited to projected angular correlation functions and rely on overall redshift distributions (see Section 2.3). D16 presented a detailed description of angular correlation functions, including the underlying theory, methods for calculation, and model generation444IDL codes used to produce our measurements are available at https://github.com/mdipompe/angular_clustering; see D16.. We refer the reader to their section 3.1 for complete details, and provide a brief summary here.

The angular correlation function $\omega_{\theta}$ quantifies the probability that a pair of objects separated by angular distance $\theta$ and with mean number density $n$ will be found within the solid angle $d\Omega$ (Totsuji & Kihara, 1969; Peebles, 1980):

[TABLE]

Massive galaxies that host quasars are found in the largest potential wells of the dark matter distribution, and are therefore biased relative to the underlying density field (Kaiser, 1984; Bardeen et al., 1986). We quantify this via the quasar bias $b_{q}$, which is tied to the underlying dark matter angular correlation function: $\omega_{q}(\theta)=\omega_{\textrm{dm}}(\theta)b_{q}^{2}$, under the limiting assumption that the bias is not a strong function of redshift or spatial scale. In the case of a cross-correlation between two samples, rather than the autocorrelation of one, this becomes $\omega_{\textrm{cross}}(\theta)=\omega_{\textrm{dm}}(\theta)b_{1}b_{2}$, where $b_{1}$ and $b_{2}$ are the biases of the two populations of interest (we will refer to the term $b_{1}b_{2}$ as $b_{\textrm{cross}}^{2}$ below).

To model $\omega_{\textrm{dm}}(\theta)$, we produce the non-linear matter power spectrum $P(k,z)$ using CAMB555Code for Anisotropies in the Microwave Background (http://lambda.gsfc.nasa.gov/toolbox/tb_camb_ov.cfm), via our IDL wrapper camb4idl which can be found at https://github.com/mdipompe/camb4idl (Lewis, Challinor & Lasenby, 2000). On scales $<<1$ radian, Limber’s approximation can be used to project the matter power spectrum onto the sky (in a flat Universe; Limber, 1953; Peebles, 1980; Peacock, 1991; Myers et al., 2007):

[TABLE]

Here, $\Delta^{2}(k,z)=(k^{3}/2\pi^{2})P(k,z)$ is the dimensionless power spectrum at redshift $z$ and wavenumber $k$, $J_{0}$ is the zeroth-order Bessel function of the first kind, $\chi$ is the comoving distance along the line of sight, and $dz/d\chi=H_{z}/c=(H_{0}/c)[\Omega_{\textrm{m}}(1+z)^{3}+\Omega_{\Lambda}]^{1/2}$. The $dN/dz$ terms represent the redshift distributions of the samples of interest (Section 2.3), and in the case of an autocorrelation are identical.

The actual angular autocorrelation measurements are made by comparing the number counts of quasar pairs on different scales with what is expected for a random distribution (Landy & Szalay, 1993):

[TABLE]

where $DD$, $DR$, and $RR$ are the normalized numbers of data-data, data-random, and random-random pairs in each bin of $\theta$ (i.e. $DD=DD(\theta)=N_{\textrm{data pairs}}/(N_{D}N_{D})$). In the case of a cross-correlation, the $2DR$ term is split into $D_{1}R_{2}+D_{2}R_{1}$ and the $RR$ terms become $R_{1}R_{2}$, where the subscripts indicate the distinct samples in their corresponding random catalogs.

The random data must have the same angular distribution as the real data for an accurate comparison in the presence of selection effects. In previous work over a limited footprint the quasar samples were consistent with being uniformly distributed over the field (after masking), simplifying the generation of the random catalog using mangle (Hamilton & Tegmark, 2004; Swanson et al., 2008). In the expanded footprint used here, there are position-dependent fluctuations in the quasar density, and we apply corrections to an initially uniform random sample (see Section 3.2.1). In all cases, the final random catalog is always at least 50 times larger than the real data set so that Poisson noise in the random counts do not limit the precision of the measurement. In general, the angular correlation function is calculated in four bins per dex, from $\sim$7 arcseconds to 1.1 degrees, unless otherwise noted.

We use inverse-variance weighted jackknife resampling to estimate uncertainty on the angular correlations (e.g. Scranton et al., 2002; Myers et al., 2005, 2007, D16), with the footprint broken into $N=100$ equal-area regions. We iteratively remove one region at a time and repeat the measurement in the remaining area. This provides the covariance matrix $\textrm{{{C}}}_{ij}=\textrm{{{C}}}(\theta_{i},\theta_{j})$ ($i$ and $j$ denote angular size bins):

[TABLE]

where $\omega$ is the angular correlation function for the entire sample and $\omega_{L}$ is the same for subset $L$. Here, the $RR$ terms are not normalized by the size of the random sample, allowing the ratio to compensate for the different number counts in each region in the case where the areas are not precisely equal. The square-root of the diagonal elements of the covariance matrix are adopted as the 1-$\sigma$ uncertainty on the measurements.

In this work we introduce an updated fitting technique beyond the simple $\chi^{2}$ minimization used in prior work. The motivation for this update is that the uncertainties on these types of clustering measurements are not necessarily Gaussian, which can impact the results of the fit. We note that in the case of our measurements here, the distribution of $\omega(\theta)$ values at each angular scale are roughly Gaussian — however, we introduce this fitting procedure here in case it is useful to the community and so that our procedure is more general for future work. It also allows us to generate a posterior distribution of bias values, which may not be purely Gaussian and therefore not fully parametrized with a single value and error bar.

The jackknife resampling provides 100 measurements of $\omega(\theta)$ at each angular scale. We draw randomly from these values to create a new autocorrelation measurement, and assign Poisson errors to each based on the number of pairs contributing to each scale:

[TABLE]

where $DD$ is the unnormalized number of pairs, and $f=2$ in the case of an autocorrelation (because pairs are not independent) and $f=1$ in the case of a cross-correlation. We then fit our dark matter autocorrelation model using standard $\chi^{2}$ minimization, and repeat the process $n=10,000$ times to build a distribution of bias values. When single values and uncertainties are reported, they represent the median and symmetric 67% confidence interval of this distribution666We note that this technique is a pseudo-Bayesian approach. We have also implemented a fully Bayesian MCMC approach (based on e.g. Hogg, Bovy & Lang, 2010) — however fully incorporating non-Gaussian uncertainties into this method is computationally expensive, whereas the method used here produces very similar results while still allowing for arbitrary uncertainty distributions..

3.2 Results

3.2.1 All WISE Quasars

We first measure the angular autocorrelation of the largest quasar sample to date — the 581,728 selected with WISE from the entire sky (after masking). However, expanding beyond the initial pilot region of D16 presents a complication in that the quasar density is not uniform across the sky, as shown in the top panel of Figure 4. There is a decrease in the number of quasars selected closer to the Galactic plane, even after removing $|b|<20$ degrees, likely because the stellar density increases and it becomes more likely that a Galactic star falls within a WISE resolution element of a distant quasar (this was also noted in e.g. Krolewski & Eisenstein, 2015). In the lower two panels of Figure 4, we normalize histograms of the Galactic coordinates of the quasar sample and a uniformly distributed random sample and take their ratio (dashed lines) — it is clear that the systematic variation is with Galactic latitude, while variations in Galactic longitude appear mostly dominated by noise.

In order to correct the random sample such that it follows the data distribution, we first pixelate the sky into equal area regions using HEALPix (Górski et al., 2005), with $n_{\textrm{side}}=8$, which has $\sim$7 deg2 pixels. We explored different pixel sizes and settled on the smallest pixel size possible without obvious “over-fitting” such that intrinsic variance was being traced. We calculate the quasar density in each pixel (taking care to trace the available area of each pixel in the presence of our mask), and take the mean in pixels with the same central latitude. Normalizing these to a maximum value of unity provides weights for each pixel as a function of $b$, which are then used to remove an appropriate number of random sources from each pixel. This provides a resampled random catalog that better matches the distribution of the data, as can be seen in the solid lines of the lower panels of Figure 4.

The angular autocorrelation for the full WISE-selected sample is shown in the top panel of Figure 5. With such a large sample, we find an excellent clustering signal with error bars smaller than the data points in many cases. The data match the shape of the model remarkably well on scales where the two-halo term dominates (where pairs are largely quasars in distinct haloes rather than subhaloes, on scales of a few tenths of a Mpc — note that at $z=1$, 0.01 degrees corresponds to $\sim$0.4 Mpc/$h$ comoving), including the turnover on scales approaching one degree. On small scales, there appears to be significant signal from the one-halo term, which can be used to probe the quasar halo occupation distribution and satellite fraction (though see the discussion in the next section regarding these scales for the obscured and unobscured populations; Zheng et al., 2005; Zheng & Weinberg, 2007; Chatterjee et al., 2013).

Fitting the dark matter model to these data (on scales $>0.025$ degrees, which will be adopted throughout unless otherwise noted) results in a WISE-selected quasar bias of $1.98\pm 0.05$ (Table 1).

3.2.2 Obscured and Unobscured Quasars

We now turn to the angular autocorrelation of WISE quasars split into obscured and unobscured subtypes — these results are shown in the middle panel of Figure 5. We apply a similar correction to the random sample in the SDSS footprint as for the full WISE sample, though the fact that the SDSS footprint does not extend as far in to the Galactic plane means that the correction is not as significant. Again, we overall see excellent agreement with the shape of the models, including the turnover on large scales.

Fitting the models to these data results in $b_{\textrm{unob}}=1.77\pm 0.11$ and $b_{\textrm{obsc}}=2.15\pm 0.12$ (Table 1). These are consistent with D16, but the reduced uncertainties due to the increase in sample size means that this difference is $\sim 3.3\sigma$. Using our new procedure for fitting the bias, we also extract posterior probability distributions for the bias of each sample, shown in the lower panel of Figure 5. We note that these distributions are not purely Gaussian, and tend to have a slight asymmetry with a larger tail toward lower bias values. There is only a small amount of overlap between the obscured and unobscured distributions — integrating the obscured PDF where it overlaps the unobscured distribution yields a total overlap probability of 0.08, while the unobscured overlap with the obscured distribution is 0.05.

The obscured autocorrelation does show an anomalous increase on scales of $\sim 0.1$ degrees, with this point falling roughly 3$\sigma$ from the best-fit line. While this could just be a statistical fluctuation, as a check that this is not due to contaminants that specifically affect this scale (e.g. haloes around bright stars, globular clusters, etc.), we visually inspected 100 random SDSS images of objects contributing to this angular bin, as well as another larger bin for comparison. While this is a small fraction of the roughly 100,000 pairs that contribute to this bin, we found no evidence of additional contaminants on this scale. To test whether this single point is a strong driver of the enhanced obscured bias, we re-perform the fit with this scale excluded. This results in an obscured bias reduced slightly to $b_{\textrm{obsc}}=2.06\pm 0.16$, which is consistent with the full measurement and still larger than the unobscured bias.

Additional signal on small scales, again possibly due to the one-halo term and providing valuable information on the satellite fraction (e.g. Mitra et al. in prep), is present in both of these subsamples as well. However, we do note that in D14 it was found that after additional masking of the WISE quasar sample (as compared to Donoso et al., 2014) the autocorrelation signal was reduced most prominently in the obscured sample on small scales. This suggests that we exercise caution when analyzing these scales, as they are sensitive to artifacts in the data. The fact that the unobscured sample shows this excess signal, and was not affected strongly by the additional masking in D14, however suggests that it may be real. In addition, D16 found that incorporating the updated ALLWISE data release, which we also use here, removed a significant amount of contamination, based on comparison of the positional dependence of the quasar density and the signal on increasingly large scales.

It is possible that the non-optically detected subset of the obscured sample (roughly 40% of the full obscured sample) suffers more contamination and affects our measurement of the bias. We repeat the autocorrelation measurement using only the optically-detected obscured sample, and find that they are nearly identical. The increased signal at 0.1 degrees is largely smoothed over, and the resulting bias measurement is $b_{\textrm{obsc,opt}}=2.18\pm 0.19$.

3.2.3 WISE Obscured-Unobscured Cross-Correlation

Because we are only able to estimate the redshift distributions of the samples based on data in the Boötes field, it is informative to explore the cross-correlation of the obscured and unobscured populations, which carries extra information about the overlap of the samples in redshift space. We show the cross-correlation compared to the autocorrelations in Figure 6. Qualitatively, the cross-correlation behaves as we would expect for samples well-matched in redshift, with the clustering amplitude falling generally in between the individual autocorrelations.

As the model curves in Figure 6 (dotted lines) show, the cross-correlation model is slightly lower in amplitude than the autocorrelation models due to small differences in the redshift distributions (Figure 2). Thus, when fitting the model to the cross-correlation, a slightly larger rescaling is needed. Since the cross-correlation signal is proportional to $b_{\textrm{obsc}}b_{\textrm{unob}}$, we can fit the model to the data, and divide by the biases inferred from the autocorrelations to determine the bias of one sample or the other. These results are listed in Table 2, and while the decreased model amplitude for the cross-correlations produces slightly larger bias values for each sample, these results are quite consistent with the results from the autocorrelation alone.

If we assume that the unobscured redshift distribution is better constrained than that of the obscured sample, which is not unreasonable given the optical brightness of the unobscured objects, we can use the cross-correlation and models generated using simulated obscured redshift distributions to inform us about the true obscured $dN/dz$. To this end, we generate Gaussian $dN/dz$ distributions with standard deviation 0.55 (similar to the standard deviation from sources in the Boötes field) and means $\mu_{2}$ from 0.1 to 2 in steps of 0.1. Note that negative values of $z$ are not possible, so these are discarded and the positive portions of the distributions are properly re-normalized such that $\int\frac{dN}{dz}dz=1$ in all cases.

In the top panel of Figure 7, we show model dark matter cross-correlations using the observed unobscured $dN/dz$ and the mock obscured distributions, normalized to unity at $\sim$0.035 degrees. Note that the amplitudes of the models vary, but we show the normalized versions here to highlight that there is also variation in the shapes depending on the $z$ overlap.

In the second panel of Figure 7 we show the $\chi^{2}$ values obtained by fitting each of these models to the cross-correlation, which provides information on how well the shape (not amplitude) of the data matches the various models. There is a clear minimum when $\mu_{2}\sim 1$, suggesting that the obscured sample is not at significantly different mean redshift than the analysis of the Boötes field implies.

Taking this one step farther, and assuming that the unobscured autocorrelation (and resulting bias) is robust, we can use the above models to estimate the unobscured bias from the cross-correlation for each $\mu_{2}$, and search for the model that provides the best agreement with the unobscured autocorrelation. In addition to the model cross-correlations, we also generate model obscured autocorrelations for each $\mu_{2}$. We fit these model autocorrelations to the obscured data, and fit the corresponding model cross-correlation to the cross-correlation data, which provides the unobscured bias as a function of $\mu_{2}$:

[TABLE]

We compare these results with the unobscured bias derived from the unobscured autocorrelation alone in the last panel of Figure 7. Again, we see that the best agreement from the two approaches is when the mean redshift of the obscured sample is at $\sim$1. This is additional confirmation that the obscured $dN/dz$ we infer from the Boötes field is not significantly skewed by observational biases.

Figure 8 shows the obscured bias inferred from the autocorrelation with each mock $dN/dz$, with the values from the autocorrelations using the Boötes redshift distributions overlaid. In order for the obscured bias to be consistent with that of the unobscured sample, the mean redshift would have to be less than $\sim$0.8, which seems unlikely based on the observational data and the modeling of the cross-correlations in Figure 7. To further quantify this, we explored artificially inflating the obscured $dN/dz$ values (Figure 2) at $z<0.8$ until the overall mean was less than 0.8 — in order for the true obscured mean to be this low, we would need to be missing roughly half of the population at low redshift.

Finally, we note that we have also performed these analyses using the obscured redshift distribution inferred from the Boötes field and shifting its mean while keeping the same shape, rather than using mock Gaussians. The overall conclusions do not change.

3.2.4 Cross-correlation with SDSS Spectroscopic Quasars

In order to study the obscured and unobscured biases in an independent way, we also cross-correlate each WISE sample with spectroscopically confirmed quasars from the SDSS. We begin with the sample of Shen et al. (2009), which has a well-studied bias and complementary random catalogue that mimics the data distribution (Shen, private communication). We first confirm that we can reproduce the bias results of Shen et al. (2009) with an angular measurement, as the original work was done in three dimensions since the spectroscopic sample has individual source redshifts. This is shown in the top panel of Figure 9. Despite the additional scatter in the angular measurement, we find good agreement with Shen et al. (2009) — our inferred bias is $b_{\textrm{spec}}=2.7\pm 0.5$, compared to $\sim 2.4\pm 0.3$.

We limit the spectroscopic sample to the footprint available for the WISE data, discarding the Southern Galactic Cap where there is very little overlap between the samples. We also only include sources with $z_{\textrm{spec}}<3$, where there is overlap with the WISE quasars. We confirm that we find a consistent bias using this subset ($b_{\textrm{spec}}=2.8\pm 0.6$), and we adopt this value for the remainder of the analysis.

The cross-correlations between the spectroscopic and WISE-selected obscured and unobscured quasars are shown in the lower panel of Figure 9. Since the density of spectroscopic quasars is fairly low, the cross-correlation has larger uncertainty than the autocorrelations. However, fitting models to these data and using the spectroscopic autocorrelation results, we find bias values lower than but generally consistent with those of the autocorrelations, as shown in Table 2.

4 CMB Lensing Correlations

4.1 Methodology

As CMB photons travel to us from the surface of last scattering at $z\sim 1100$, they are deflected by the matter distribution along the line of sight. By cross-correlating maps of the quasar density with maps of the CMB lensing convergence $\kappa$, the quasar bias $b_{q}$ can be determined in yet another independent way. As with the angular correlation functions, the details of the theory behind these measurements (Bleem et al., 2012; Sherwin et al., 2012), the calculations, and generation of models are discussed in detail in D16. We refer the reader there for complete details, and provide a basic summary here777See also our code library for these measurements at https://github.com/mdipompe/lensing_xcorr..

For a quasar sample with bias $b_{q}$ (which may be a function of redshift and therefore distance) their host dark matter halo distribution kernel is:

[TABLE]

where $\chi$ is the comoving distance. The CMB lensing kernel, which describes the strength of lensing as a function of distance to the lens, since the CMB is at a fixed distance, is given by (Cooray & Hu, 2000; Song et al., 2003):

[TABLE]

where $a(\chi)$ is the scale factor and $\chi_{\textrm{CMB}}$ is the co-moving distance to the CMB. The cross-correlation in Fourier space at mode $\ell$ of these two components is given by:

[TABLE]

where $P(k=\ell/\chi,z)$ is the matter power spectrum, again generated using CAMB. Setting the bias in the distribution kernel equal to unity provides the model cross-power spectrum for the underlying dark matter distribution.

We utilize the all-sky CMB lensing potential map from the Planck second data release (Planck Collaboration et al., 2016a, b) to perform the cross-correlation. Using tools within HEALpix, the native format of the Planck data, we convert the provided coefficients ($A_{\ell m}$) of the spherical harmonic transform of the lensing potential $\phi$ into an $n_{\textrm{side}}=2048$ (resolution of $\sim$1.7 arcmin) map of the lensing convergence ($\kappa=1/2\nabla^{2}\phi$).

The quasar positions are converted into a HEALPix map (with the same resolution) of the relative quasar density:

[TABLE]

where $\rho$ is the density in a given HEALPix pixel and $\bar{\rho}$ is the mean density. To account for the fact that mean density changes as a function of position (see the top panel of Figure 4), rather than use the overall mean density over the full footprint, the relative density is calculated with respect to the local mean density in larger ($n_{\textrm{side}}=8$ or $\sim$7 deg) pixels.

Finally, we take the Fourier transform of each map and multiply them, keeping only the real portion, using the HEALpix routine anafast. We bin the cross-power in $\ell$, using 4 bins per dex, as with the angular correlation functions. Uncertainties on the cross-correlation are obtained by repeating the measurement with several rotated Planck lensing maps, which should have null cross-correlation signal in the absence of noise and systematics. We use 17 rotations in Galactic longitude in steps of 20 degrees, and 17 more with an additional 180 degree rotation in latitude, for a total of 34 rotated maps (see DiPompeo et al., 2015b, and D16 for more details on this procedure to generate uncertainties). Cross-correlating with these maps provides the covariance matrix:

[TABLE]

where the terms in parentheses are the cross-correlations with a given map rotation $k$ and $N=34$.

To fit the models generated with equation 9 to the data, we use a similar approach as for the clustering measurements. We draw randomly a value of $C_{l}$ for each $\ell$ from the cross-correlations with the rotated maps, and add it to the actual measurement at that $\ell$. We then use typical $\chi^{2}$ minimization to fit the model to this iteration, using the full covariance matrix in the fit (we are not able to generate Poisson errors as in the clustering measurements because it is not possible to separate out how many quasars contribute to the measurement at a given $\ell$). Repeating the process $n=10,000$ times allows us to build a posterior distribution of bias values.

4.2 Results

While we defer a detailed discussion to Section 6, we have included an additional cut on the data for the CMB lensing cross-correlation, limiting our analysis to Galactic latitude $|b|>40$ degrees. This is due to apparent increasing systematic effects with the expanded footprint that results in no significant improvement in the precision, as well as systematically lower bias values, particularly in the unobscured sample ($b_{\textrm{unob}}=1.5$, quite low for quasars selected in any way). Our “best” final result could be considered to be the result with both the best accuracy and the best precision. Precision might be expected to be driven by the total number of objects included in our analysis. Accuracy, however, will be a strong function of the systematics that contaminate our measurements. Our main sources of systematics (likely Moon levels, stellar density, Galactic dust, increasing contamination, etc.) tend to be most prevalent at low Galactic latitudes. Making this additional cut, while reducing the number of objects contributing to the measurement, leads to similar precision as the full footprint (and still slightly worse than D16). We note that the clustering results are not strongly impacted by applying a similar cut.

Figure 10 highlights the high level of correlation between the quasar density maps and CMB lensing convergence. We identify pixels within bins of $\delta_{q}$, and calculate the mean value of $\kappa$ in these pixels — clearly, as the quasar relative density increases, the lensing convergence increases as well, showing that the quasar host haloes are indeed lensing the CMB photons. Differences in the shape of this relationship for different samples will be discussed in Section 6.

The formal quasar and CMB lensing cross-correlations in this high Galactic latitude measurement are shown in the top panel of Figure 11, where the grey points show the results for the complete WISE sample not limited to the SDSS footprint, and the blue and red points are the unobscured and obscured results over the more limited footprint, respectively. A few important things are evident in this figure. First, the obscured quasar cross-correlation signal is higher on nearly all scales. Second, there is an increase in the obscured cross-correlation signal on scales of $\sim$0.2 degrees, which is similar to the scales at which additional signal was present in the clustering measurement (the offset between the two could simply be due to the way the data are binned). Third, it appears that for all data sets the data diverge from the model on scales of a few degrees, such that the fit to the data on smaller scales over-predicts the measurement on large scales. This discrepancy appears to begin on slightly smaller scales in the unobscured sample. This was not seen in previous measurements (e.g. D16).

The biases resulting from fitting the data on scales of $100<\ell<1000$ are given in Table 1, and the posterior distributions are shown in the lower panel of Figure 11. While the distributions are wider, the two samples clearly have distinct peaks, with a significance of $\sim$2$\sigma$. The obscured bias distribution overlaps the unobscured by 24%, and the unobscured overlaps the obscured by 15%.

Due to the discrepancy on different scales discussed above, the regions over which we perform the fit to the data impact the final results more so than in previous work. For example, fitting the same region as D16 ($40<\ell<400$), the bias for the obscured sample is reduced to $1.87\pm 0.18$ and the unobscured bias is reduced to $1.31\pm 0.15$ — note that the unobscured sample is affected more strongly. These measurements are consistent with the others considering the uncertainties, but particularly for the unobscured sample are quite low compared to any previous measurements, including samples selected using other multiwavelength or spectroscopic methods.

Finally, we also note that unlike for the clustering measurements, expanding the footprint of the measurement does not reduce the uncertainty on the final bias measurement, though this does depend on what scales are fit. While the error bars on individual points are reduced, roughly by a factor of $\sqrt{N}$, where $N$ is the difference in sample sizes, over the $100<\ell<1000$ range the final $\chi^{2}$ values are not significantly changed as the quality of the fit is not improved. However, no matter which region is fit, the general result that the obscured sample has a higher bias does not change.

5 Dark matter halo masses

Using model fits to cosmological simulations, the large-scale bias of a sample can be converted into the typical or “effective” dark matter halo mass of the population. For a more complete description of this methodology, see D16888And the code at https://github.com/mdipompe/halomasses.. We use the bias parameterization of Tinker et al. (2010, see their Table 2 for the values of the coefficients), using overdensity parameter $\Delta=200$:

[TABLE]

where $\nu=\delta_{c}/\sigma(M)$. The critical density for collapse of a halo $\delta_{c}$ is:

[TABLE]

where $\Omega_{m,z}$ is the matter density parameter at $z$ (Navarro, Frenk & White, 1997). The linear matter variance on the size scale $R_{h}$ of the halo $\sigma(M)$ is

[TABLE]

The matter power spectrum $P(k,z)$ is calculated with CAMB, and $\hat{W}(k,R_{h})$ is the Fourier transform of a top-hat window function of radius $R_{h}$:

[TABLE]

The halo masses calculated from the bias measurements are listed in Table 1. Note that due to the shape of $b(M)$, which flattens at low halo mass and rises steeply at high mass, the uncertainties on $M_{h}$ are generally asymmetric, and are larger for low halo masses.

6 DISCUSSION

Figure 12 shows the results of this analysis (filled symbols) along with other recent work on obscured quasars, including Hickox et al. (2011) (triangles), and D16 (light open symbols). Our final bias and halo masses for the obscured and unobscured samples are consistent with previous work (e.g. D16) , but the significance of the difference between them is increased. While the measurements all rely on the same underlying data set (the WISE-selected quasars), they are made in independent ways, with different systematics that allow us to combine them for a robust comparison. Before combining the results, we repeat the CMB lensing cross-correlation analysis in the 100 sub-regions used for the jackknife clustering uncertainties (for the full WISE sample), and fit the bias for both measurements in each iteration. This provides a direct test of the independence of the clustering and CMB lensing cross-correlations. Interestingly, the final values are not correlated, with $r=0.11$ and $p=0.26$ (using a simple Spearman rank correlation). This is likely because the error budget in each case is dominated by different sources — based on the analysis above, it is possible that the clustering result is dominated by Poisson noise, while the CMB lensing is dominated by systematics in either or both the CMB lensing and quasar density maps. If we combine the measurements here, we find weighted mean values of $b_{\textrm{obsc}}=2.15\pm 0.11$ and $b_{\textrm{unob}}=1.70\pm 0.10$, or $M_{h\textrm{,obsc}}=12.94\pm 0.08$ $M_{\odot}$ $h^{-1}$ and $M_{h\textrm{,unob}}=12.49\pm 0.11$ $M_{\odot}$ $h^{-1}$, a difference of $\sim$5$\sigma$. These values are listed in Table 1.

As mentioned in Section 4.2, it appears that additional systematics in WISE quasar selection may become more important over large areas, affecting in particular the CMB lensing measurements. This manifests in bias uncertainties consistent with or larger than previous measurements and systematically lower bias values, which we mitigate somewhat with the more stringent cut in Galactic latitude. Interestingly, the additional systematics appear to impact the unobscured sample more strongly than the obscured.

While we have corrected for the fluctuating WISE-selected quasar density as a function of Galactic latitude, it appears that contamination may become an increasing problem over wider areas. While Stern et al. (2012) find a reliability of 95% for their selection criteria in the COSMOS field, it is possible that this is not true across the entire sky. In addition, it is also the case that simply expanding the area allows more rare contaminants to appear in the sample. While the amount and type of contamination is difficult to determine from these measurements alone, there are some indications, particularly in Figure 10, that this could be the case. In general, the trend between $\delta_{q}$ and $\kappa$ will flatten for smaller samples, as Poisson noise dilutes the signal, particularly in the more extreme bins (this is clearly seen by simply randomly down-sampling one of the samples and comparing it to the full measurement). However, this is not the case at the low-$\delta_{q}$ end in Figure 10, where the larger unobscured sample flattens more rapidly than the obscured. This is in contrast to the high $\delta_{q}$ side, where they are more consistent. In a similar figure in DiPompeo et al. (2015b), which covered a smaller area, this difference is not apparent. This suggests that there are more contaminants in the unobscured sample over the expanded footprint, which are scattered randomly across the area so that they are likely to fall in low quasar density pixels, and these are uncorrelated with the CMB.

The most likely such contamination comes from stars, though we find this unlikely to be the case, at least for typical stars. As a simple check, we match our unobscured quasar sample with the catalog of DiPompeo et al. (2015a), which provides quasar and star probabilities for all point sources in the SDSS, based on their mid-IR through UV fluxes. Since resolved sources are not included in the catalog, we find 79% of our unobscured sample in this catalog. However, comparing the distributions of quasar probabilities $P_{\textrm{QSO}}$ over the full footprint and the region of D16 shows that they are nearly identical, with $\sim$97% of the sources in both regions having $P_{\textrm{QSO}}>0.9$. Any stellar contamination is thus likely to be from rare stars with an AGN-like spectral energy distribution (e.g. AGB stars), and therefore difficult to disentangle without extensive spectroscopy. Finally, we also create HEALPix maps at several resolutions of the unobscured to obscured ratio, and find no strongly inconsistent outliers, nor any trend with sky position, suggestive that potential stellar contamination affects both samples equally.

We also explored a more conservative cut of $W1-W2>1.0$, for which Stern et al. (2012) report a reliability of nearly 100%. This reduces all of the samples by $\sim$40%. While such a cut also impacts the redshift distributions of the samples, changing the interpretation of the bias results, we use it simply to check if it improves any of the signs of contamination discussed so far. The binned $\kappa$-$\delta_{q}$ relationship does appear to behave more as expected, and the outlier in the unobscured sample in the lowest $\delta_{q}$ bin is consistent with expectations based just on sample size. However the CMB cross-correlation bias results remain sensitive to the choice of fitting region, with a seemingly more rapid turnover on large scales, particularly for the unobscured sample.

Regardless of these issues, we stress that the difference in bias and halo mass for the two populations continues to persist, even if the final adopted values fluctuate. This is not only bolstered here by the reduced uncertainties on the clustering measurement, but the additional analysis of angular cross-correlations between the WISE-selected samples and with SDSS spectroscopic quasars. The fact that the WISE obscured-unobscured cross-correlation signal lies generally within the individual autocorrelations is strong support for the well-matched redshift distributions based on the Boötes field. If the samples were mismatched in redshift, we would expect to see a rapid decrease in the cross-correlation signal, as evidenced by the lower amplitude model generated for the cross-correlation with the very similar but not perfectly overlapping Boötes distributions. Finally, though the uncertainties are large, the cross-correlation with spectroscopic quasars supports the matching in redshift, as well as the final bias results. All of these taken together are strong evidence that the bias and halo mass difference we find is real.

The most common concerns with these measurements for WISE quasars are that the redshift and/or luminosity distributions of the obscured and unobscured populations are not well matched. The cross-correlations discussed above strongly favor a well-matched redshift distribution, as does the modeling we performed using mock obscured $dN/dz$ distributions at different mean redshifts (Figures 7 and 8). This, combined with the SED analysis of Hickox et al. (2007) and the analysis of the $W3$ and $W4$ magnitude distributions here (Figure 3) strongly suggest a well-matched bolometric luminosity distribution. In addition, as discussed previously, there is little evidence for a strong dependence of the bias on luminosity, implying that even if the samples are at slightly different mean $L_{\textrm{bol}}$, this likely isn’t enough to explain the observed difference in halo mass. However, this may not be the case on small scales (e.g. Degraf & Sijacki, 2017), so care will need to be taken in future analysis of the difference in clustering signal on small scales.

There is general consistency between D16 and this work, especially for the obscured quasars. The unobscured bias, particularly from the CMB lensing analysis, is somewhat low, both compared to D16 and other work on optically selected unobscured quasars more generally. Compared to Hickox et al. (2011), we again find general agreement. We also note that our measurement of the overall WISE-selected quasar bias (not shown in Figure 12 for clarity, but given in Table 1) agrees well with a cross-correlation with a CMB lensing map from the South Pole Telescope (Geach et al., 2013). The interpretation in that work was that the similarity of the overall IR-selected quasar bias to the well-studied optically-selected quasar bias implied that the obscured and unobscured populations were similar. However, we see here that while the overall WISE quasar bias agrees with the mean optical quasar bias, the IR-selected unobscured bias is slightly lower and the obscured bias is slightly higher.

Finally, we point out that this difference in bias and halo mass, supported from many lines of evidence, is likely a lower limit on the value for a truly evolutionarily distinct population (DiPompeo et al., 2016). This is simply because the hot dust that makes up the torus is likely to obscure some quasars that are otherwise intrinsically the same as the unobscured population — i.e. the obscured population is made of of some objects unified by orientation and some unified by evolution or other factors.

Given that extending the survey area seems to lead to an increase in contamination that prevents us from improving the precision of bias and halo mass measurements of quasars selected with simple IR criteria, at least for CMB lensing cross-correlations, even in the era of expanding optical surveys (e.g. Pan-STARRS, LSST) that would allow an obscured/unobscured split over larger areas there may not be much more room for improvement simply be increasing sample sizes. Of course, such surveys also provide additional depth and time-domain information that will certainly be useful for understanding contamination. Instead, for now we will look towards smaller samples with redshift information, either spectroscopic or photometric, as well as cross-correlations with increasingly large galaxy and spectroscopic quasar catalogs.

7 Summary

Our goal in this work is to improve the precision in the bias and halo mass measurements of WISE-selected quasars by expanding the footprint of previous measurements. This is to provide the most accurate WISE-selected quasar measurement to date, as well as due to the important but not highly significant difference in halo mass between obscured and unobscured populations identified in previous work. After correcting for observational systematics that lead to a decrease in the WISE-selected quasar density closer to the Galactic plane, we measure the angular clustering and cross-correlations with CMB lensing maps for both samples.

The precision of the angular clustering measurement improves as expected. However, the CMB lensing cross-correlations do not improve in precision, possibly due to increased contamination (particularly in the unobscured sample), though it is difficult to say for certain. Combining results from the measurements leads to a difference in halo mass of a factor of 3 at a significance of $\sim$5$\sigma$. We also perform a cross-correlation of both samples with a sample of spectroscopic quasars, also confirming the halo mass difference. This implies that the simplest unification models that rely only on orientation are not sufficient to explain the full quasar population.

We also measure the cross-correlation of the obscured and unobscured populations, and find that this gives consistent halo mass measurements. Using models of the redshift distribution of the obscured population at different mean $z$, we show that the cross-correlations strongly suggest that the obscured and unobscured populations have very similar redshift distributions, both centred around $z\sim 1$. We also compare the 12 and 22 $\mu$m WISE detection fractions and flux distributions and find them to be very well matched, suggesting good agreement in the bolometric luminosities of the obscured and unobscured populations selected in this way. These analyses highlight that it is fair to directly compare the bias and halo mass measurements of the two samples, even if these parameters evolve with redshift and/or luminosity.

The highly significant difference in obscured and unobscured quasar halo mass reported here, using several lines of evidence along with careful testing of the match in fundamental properties, suggests very strongly that obscuration in at least a subset of IR-selected quasars is tied to the host galaxy. The simplest explanation for this is a link between the evolutionary state of the galaxy and the availability of material to obscure the central engine.

Acknowledgements

We thank Andy Goulding for his valuable input, in particular with regards to the new fitting methods used here. MAD, RCH, and ADM were partially supported by NASA through ADAP award NNX12AE38G and by the National Science Foundation through grant number 1211096. MAD and ADM were also supported by NSF grant number 1211112, and MAD, ADM, and SE were supported by NSF grant 1515404. MAD and RCH were also supported by NSF grant number 1515364 and NASA ADAP award NNX16AN48G. RCH also acknowledges support from an Alfred P. Sloan Research Fellowship, and a Dartmouth Class of 1962 Faculty Fellowship.