Host Dark Matter Halos of Wide-field Infrared Survey Explorer-selected Obscured and Unobscured Quasars: Evidence for Evolution

Though heavily obscured quasars can be challenging to distinguish from normal galaxies in the optical wave band, they can be recovered simply via their red mid-IR colors (Lacy et al. 2004; Stern et al. 2005, 2012; Donley et al. 2012; Assef et al. 2013), which trace the reprocessed emission from the dusty torus. The Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) has thus revolutionized obscured quasar selection by mapping the entire sky in four mid-IR bands centered at (named) 3.6 (W1), 4.5 (W2), 12 (W3), and 22 μm (W4). These data currently provide the only diagnostic capable of producing highly reliable samples of millions of both obscured and unobscured quasars over the entire extragalactic sky.

We therefore elect to use the Assef et al. (2018) 90% reliable (R90) criterion on WISE photometry to select our parent sample of quasars. Rather than adopting the publicly released catalog associated with Assef et al. (2018), which used the criterion to select candidates from the AllWISE catalog (Cutri et al. 2021), we apply the criterion to the newer CatWISE 2020 (Eisenhardt et al. 2020; Marocco et al. 2021) data release. This catalog is generated from WISE observations taken between 2010 and 2018, incorporating six times more exposures in the W1 and W2 channels than were used in generating the AllWISE catalog. This deeper imaging enables more precise measurements of photometric colors and therefore more reliable selection of quasars. Crucially, the deeper imaging will also enable a uniformly complete selection across the sky to fainter fluxes, greatly simplifying the reconstruction of the selection function necessary to perform a clustering measurement. We thus query the IRSA ⁵ database and select all objects in the CatWISE catalog (Marocco et al. 2021) satisfying the R90 criterion, adopting mpro Vega magnitudes. We also apply a magnitude limit in the W2 channel such that the selection is >99% complete across the entire sky, and a bright-end cut to exclude IR stars:

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\rm{W}}1-{\rm{W}}2\gt 0.65\times \exp [0.153\times {\left({\rm{W}}2-13.86\right)}^{2}]\\ 9\lt {\rm{W}}2\lt 16.5.\end{array}\right.\end{eqnarray} \tag{ 1 }$$

As the Assef et al. (2018) criterion was calibrated in the extragalactic Boötes field where Galactic objects are relatively rare, we must apply masks to remove regions which likely suffer from higher contamination rates. For this, we follow the procedure described by Assef et al. (2018), though we create multiorder coverage maps (Fernique et al. 2014) to accomplish this rather than remove sources geometrically. In particular, we mask regions within 10° of the Galactic plane and within 30° of the Galactic center. We also mask regions occupied by Galactic planetary nebulae (Acker et al. 1992), H ii regions (Anderson et al. 2014), star-forming regions (Lynds 1962, 1965), and resolved nearby galaxies in the Catalog and Atlas of the Local Volume Galaxies (Karachentsev et al. 2013) or the Two Micron All Sky Survey Extended Source Catalog (Skrutskie et al. 2006). We refer the reader to Assef et al. (2018) for more detail in this masking procedure. We will refer to this mask as the "contamination mask" throughout the remainder of this work.

One of the advantages of the CatWISE 2020 catalog over previous WISE releases is the order-of-magnitude improvement in detecting proper motions owing to the longer time baseline. This information allows the removal of stellar or solar system object contaminants from our quasar candidate sample. We thus remove sources with measured motions >025 yr⁻¹, which are able to be detected at >5σ across the sky down to our flux limit. This cut removes 1.5% of objects from the masked catalog.

2.1.1. Optical Data: Binning by Obscuration

The distribution of optical to mid-IR colors of IR-selected quasars appears bimodal, as rest-frame UV-optical emission is easily extincted by dust while near-IR emission is less so. Therefore, a simple color cut can be used to classify obscured systems. Quasars with optical-IR colors r−W2 >3.1 [AB] typically show X-ray absorption corresponding to absorbing column densities of N_H > 10²² (Hickox et al. 2007). In order to classify our sources as "obscured" or "unobscured," we sought an optical survey deep enough to detect most of the quasars and as wide as possible to maintain a large sample size. We therefore utilize r-band optical data from the ninth release (DR9) of the Dark Energy Spectroscopic Instrument Legacy Imaging Survey (DESI-LS; Dey et al. 2019). This survey covers ∼20,000 deg² of extragalactic sky to r_AB ∼ 24, and the unprecedented combination of depth and area ensures that we are able to estimate the degree of obscuration for a uniform sample consisting of more than one million quasars for the first time. We match our parent sample of WISE quasars to the DESI-LS catalog with a matching radius of 2'' (Assef et al. 2013) using NOIRLab's Astro Data Lab ⁶ tools. The optical photometry is corrected for Galactic reddening using the map of Schlegel et al. (1998). We retain the ∼8% of WISE quasar candidates within the DESI-LS footprint lacking optical counterparts and assign lower limits to their r−W2 color indices based on the local r-band imaging depth. To estimate the color distribution of sources undetected in DESI-LS, we match the nondetections to the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP) "DEEP" catalog (Aihara et al. 2018) in the areas where the surveys overlap. Ninety-six percent of the DESI-LS nondetections are detected in HSC-DEEP, implying a negligible contamination rate of spurious WISE sources in our sample.

The resulting optical-IR color distribution of our quasar sample is shown in Figure 1. It is clear that the WISE mid-IR selection reveals a population with a bimodal optical-IR color distribution, and thus uncovers both obscured and unobscured sources in roughly equal proportion. We adopt the criterion of Hickox et al. (2007 ; r−W2 [AB] > 3.1) to separate the sample into obscured and unobscured subsets. This results in 52% of the sample being classified as unobscured and 48% as obscured.

Zoom In Zoom Out Reset image size

In order to robustly interpret clustering measurements of quasar candidates, it is important to optimize and quantify the purity of the quasar selection. Mid-IR quasar selection can be contaminated by low-redshift star-forming galaxies (SFGs) as well as by luminous high-redshift galaxies. Barrows et al. (2021) showed that the most significant source of sample contamination in the Assef et al. (2018) R90 catalog is from z = 0.2 to 0.3 galaxies with high specific star formation rates. Messias et al. (2012) developed a criterion to discriminate between low-redshift SFGs and AGNs using a near-IR to mid-IR color, K − [4.5]. In this work, we show that an alternative three-band diagnostic can effectively isolate low-redshift star-forming contaminants, as these objects are redder at r−W2 versus z−W2 than quasars. This was discovered by observing that optically bright candidates appear bimodal in this color space, with objects redder at r−W2 compared to z−W2 having photometric redshifts z < 0.7 (Barrows et al. 2021) and appearing resolved in DESI-LS imaging.

We display our candidates in this color space in Figure 2. Using a separate color, we show candidates which are highly resolved in the optical imaging, defined as being significantly better fit with a round exponential profile than as a point source (${\chi }_{\mathrm{REX}}^{2}-{\chi }_{\mathrm{PSF}}^{2}\gt 1000$). We also show predicted colors of optical quasars using the type-1 template of Hickox et al. (2017), and of the prototypical starburst galaxy M82 from Polletta et al. (2007). It is clear that the optically bright end of the R90 quasar selection criterion includes low-redshift galaxy-dominated contaminants. We use a simple intersection of two lines to excise sources in the SFG region of color space, which removes 12% of our full sample. We note that the quasar color distribution shown in Figure 1 was produced after culling these galaxy contaminants.

Zoom In Zoom Out Reset image size

We briefly consider the reliability of this new catalog by matching it with the photometric redshift catalog produced in the Boötes field (Duncan et al. 2021). This catalog contains AGN diagnostic flags for every source denoting whether it appears in the "Milliquas" complilation of spectroscopically confirmed quasars from the literature (Flesch 2019), exhibits mid-IR Spitzer-IRAC colors indicative of AGN activity (Donley et al. 2012), or coincides with a luminous X-ray source (Kenter et al. 2005). Ninety-five percent of our sample sources pass at least one of these criteria, while only 75% of sources in the Assef et al. (2018) R90 catalog do. We attribute this reliability increase to the more precise CatWISE photometric colors as well as our removal of galaxy contaminants. Our catalog is thus expected to be highly reliable and any systematics introduced into our clustering measurements from nonquasar contamination will be subdominant.

Interpreting the angular clustering or CMB lensing signal for a sample of quasars requires knowledge of the sample's redshift distribution. Targeted quasar redshift surveys are currently limited to optically bright relatively unobscured quasars (e.g., Lyke et al. 2020; Alexander et al. 2023; Chaussidon et al. 2023), and individual blind redshift surveys lack either the breadth or depth to constrain the redshift distribution of our full sample. In order to obtain a nearly complete and statistically representative estimate of the redshifts of our sample quasars, we therefore look to the well-characterized Boötes and COSMOS (Scoville et al. 2007) fields.

We select the portion of our sample falling within the footprint of the AGN and Galaxy Evolution Survey (AGES; Kochanek et al. 2012) in Boötes or the zCOSMOS (Lilly et al. 2007) footprint in the COSMOS field. We then performed a literature search for photometric and spectroscopic redshifts of the objects in these fields by matching them against various publicly available redshift catalogs with a 2'' search radius. We prioritize spectroscopic over photometric measurements, and adopt the spectroscopic estimate with the highest-quality flag when available. We also reject photometric redshift estimates of z > 3 as unrealistic and likely the result of poor fits. In Boötes, we take photometric redshifts from Chung et al. (2014) and then Duncan et al. (2021), with the latter taking priority. In COSMOS, we adopt photometric redshifts from Delvecchio et al. (2017) and then Marchesi et al. (2016), which in turn depend on analyses from Salvato et al. (2009, 2011), Ilbert et al. (2009), and Laigle et al. (2016). The surveys providing the spectroscopic redshifts include AGES (Kochanek et al. 2012), the Sloan Digital Sky Survey Data Release 16 (Ahumada et al. 2020), zCOSMOS (Scoville et al. 2007), PRIMUS (Coil et al. 2011; Cool et al. 2013), IDEOS (Hernan-Caballero et al. 2016), FMOS-COSMOS (Kashino et al. 2019), DEIMOS10K (Hasinger et al. 2018), and the studies of Prescott et al. (2006), Trump et al. (2009), Bussmann et al. (2012), Casey et al. (2012), Lacy et al. (2013), Assef et al. (2013), Comparat et al. (2015), Kartaltepe et al. (2015), Onodera et al. (2015), and Schulze et al. (2018). Overall, we find that 100% and 97% of unobscured and obscured quasars of the 899 total sources in these fields have a redshift measurement, respectively. We thus treat the redshift distributions as representative of their respective samples.

We show the resulting redshift distributions of the obscured and unobscured samples in Figure 3. Incorporating both photometric and spectroscopic redshifts, the distributions appear similar between the samples, implying their clustering properties can be appropriately compared. We verify this by measuring a p-value of p = 0.15 in a Kolmogorov–Smirnov (KS) test, failing to reject the null hypothesis that the two samples are drawn from the same distribution. However, only photometric redshifts are available for approximately half of the obscured sample. Therefore, we will test the robustness of our clustering results against possible systematics in the obscured population's redshift distribution in Section 6.

Zoom In Zoom Out Reset image size

The CMB radiation from z ≈ 1090 has been gravitationally lensed by the intervening LSS, and thus encodes information about the dark matter halos that host galaxies across cosmic time. Therefore, measuring the cross-correlation between the lensing convergence and the angular overdensity of quasars independently constrains the sample's host halo properties. In this work, we utilize the 2018 release of the CMB lensing convergence map (Planck Collaboration et al. 2020b) provided by data from the Planck satellite to probe the halo properties of quasars as a function of obscuration. We adopt the map generated using the minimum-variance estimator, which used both temperature and polarization data to reconstruct the lensing convergence. We also make use of simulated noise maps provided with the data release to estimate uncertainties. To generate the lensing and noise maps, we transform the provided κ_{ℓ m} coefficients at ℓ < 2048 into NSIDE = 1024 resolution HEALPix maps.

Finally, we estimate the bolometric luminosity distributions of obscured and unobscured quasars to compute space densities and occupation statistics in Section 6.3. We use the rest-frame 6 μm band as a bolometric luminosity tracer that is minimally affected by obscuration. As 90% of our sample falls within the Donley et al. (2012) AGN wedge, we suggest this wavelength is also primarily tracing the torus luminosity and is uncontaminated by star formation activity. To check the location of these sources in Spitzer color–color space, we utilize mid-IR Spitzer photometry for our sources from the Spitzer Deep, Wide-field Survey (SDWFS; Ashby et al. 2009) as collated by the Herschel Extragalactic Legacy Project (Shirley et al. 2019, 2021). We display the 6 μm luminosity distributions for obscured and unobscured quasars in Figure 4.

Zoom In Zoom Out Reset image size

Clustering measurements are inherently sensitive to systematics arising from a nonuniform selection function across the sky. Therefore, in addition to the contamination mask generated following Assef et al. (2018), we create several more masks to excise regions likely to introduce systematics into our measurements. As we permit optically undetected quasars in our sample, we must carefully characterize the imaging footprint of DESI-LS to avoid classifying quasars as obscured simply because they lack deep optical imaging. Fortunately, DESI-LS DR9 provides random catalogs (Myers et al. 2023) ⁷ populated within the imaging footprint, with an associated local depth for each random point and each band. We take the median of the r-band depth in HEALPix cells, and mask cells where the 5σ depth is r < 23. Next, we create a "reddening mask," excising regions with a high degree of Galactic reddening, E(B − V) > 0.1 (Schlegel et al. 1998). All of these masks are generated with HEALPix at NSIDE = 1024 resolution.

We must also remove regions with severe WISE imaging artifacts, such as diffraction spikes, saturated pixels, latents, and ghosts. However, these regions are often localized to scales of a few arcseconds, so masking the parent NSIDE = 1024 pixel (corresponding to a 34 resolution) of every flagged region would cause us to remove an unacceptably large fraction of the sky from our analysis. Instead, we use the UnWISE bitmasks (Meisner et al. 2019) to remove sources in the data and random catalogs overlapping with flagged pixels. This process, as well as the contamination-masking procedure (Section 2.1), effectively masks a fraction of each parent pixel, which we must correct for when generating the quasar overdensity map to cross-correlate with the CMB lensing map. Thus, we compute the fractional area lost within each NSIDE = 1024 pixel, naming this the "lost-area" mask. We further apply a Boolean mask to remove pixels in which the lost area is greater than 20% (e.g., Krolewski et al. 2020).

After applying the aforementioned masks, we visually inspect the sky distribution of our sample using TOPCAT (Taylor 2005) to check for remnant artifacts not captured by the masks. We find a small number of remaining regions which suffer from egregious bright star artifacts, where the quasar candidate sky distribution mirrors the pattern of large-scale (>1°) diffraction spikes. As our candidate quasars are selected on their very red mid-IR colors, these regions appear to be affected by the reddest and brightest IR stars in the sky. We therefore find that using stars selected in the W3 band allows us to effectively mask these regions. We query the AllWISE catalog for sources brighter than W3 [Vega] < −1.5 (of which there are only ∼200 sources in the extragalactic sky), and mask regions around them. We find that using a mask diameter in degrees equal to the absolute value of the W3 Vega magnitude effectively removes the affected regions, as brighter stars tend to affect larger areas. We note that this choice represents a convenient rather than optimal solution. Finally, we observe an increase of the density of quasar candidates without optical counterparts within a few degrees of the ecliptic poles, where the WISE imaging is deepest. We thus conservatively apply a mask within 10° of each ecliptic pole. We summarize the sample selection and masking as follows. The Assef et al. (2018) R90 selection of Equation (1) curates a sample of ∼2.3 million WISE-detected sources within the ∼20,000 deg² DESI-LS footprint. We report the percent reduction for each masking step in reference to this original sample as many sources are masked by multiple masking steps. Five percent of the original sources are excised by the contamination mask generated following Assef et al. (2018; Section 2.1). Ten percent of the original DESI-LS footprint is masked by our optical imaging depth requirement, while 6% of the area is vetoed by our reddening mask. The UnWISE bitmask-flagging procedure removes 6% of sources, and the lost-area mask removes 3%. The bright-IR star mask removes 1% of the footprint, while the ecliptic-pole mask removes 2%. After removing the 12% of likely contaminants with galaxy-like colors (Section 2.2) and the 1% of stellar interlopers with large proper motions (Section 2.1), we are left with a sample of ∼1.4 million quasar candidates over an effective area of ∼15,000 deg². We display the final combined mask atop the sky density of quasar candidates as a visual representation of this summary in Figure 5.

Zoom In Zoom Out Reset image size

We apply all the above masks to both the data and random catalogs required to perform a clustering measurement. We also mask the Planck lensing map with the mask provided alongside the Planck Collaboration et al. (2020b) data release.

After applying the aforementioned masks to the data and random catalogs, we develop weights to further correct for large-scale angular systematics in the density of quasars (e.g., DiPompeo et al. 2017a). These weights will be used in both the lensing and the clustering analyses.

We might expect our sample density to vary with Galactic latitude due to increasing stellar density near the Galactic plane or ecliptic coordinates reflecting the WISE scanning pattern. Furthermore, by using optical data to bin the sample into obscured and unobscured objects, we might expect each subsample's density to vary with optical imaging depth. We measure the quasar density as a function of each of these variables and find the variation is largest for ecliptic latitude and logarithmic optical imaging depth. We thus generate two-dimensional histograms of each sample along these dimensions and assign weights as the ratio of data points to random points in each bin. We have subsampled the random catalog using the assigned weights such that the randoms match the data, and use these random catalogs in the correlation function measurements. Alternatively, for the CMB lensing cross-correlation, we assign the inverse of the weights to the data points rather than the randoms, using these weights when constructing the quasar overdensity field.

Without redshift information for each of the objects in our quasar sample, we are limited to calculating angular correlation functions, which measure the clustering projected onto the surface of the sky. The two-point angular autocorrelation function, w(θ), is defined as the excess probability—above that expected from an unclustered Poisson distribution—of finding a pair of objects at angular separation, θ (Peebles 1980). We estimate the angular correlation function in this work using the Landy & Szalay (1993) estimator:

$$\begin{eqnarray}&&w(\theta )=\displaystyle \frac{{DD}(\theta )-2{DR}(\theta )+{RR}(\theta )}{{DD}(\theta )},\end{eqnarray} \tag{ 2 }$$

where DD, DR, and RR are weighted counts, normalized by number density, as a function of separation, for data–data pairs, data–random pairs, and random–random pairs, respectively.

We measure the angular correlation function using Corrfunc (Sinha & Garrison 2020) in 15 logarithmically spaced bins between scales of 10^−2.5 < θ < 10^−0.25 degrees. This range probes both regimes in which the one-halo term dominates, as well as the two-halo term, while avoiding the ∼6'' resolution limit of WISE and also the regime at large scales where the Limber approximation begins to break down at the >10% level, which is ≳1° for our sample's redshift distribution (Simon 2007).

To estimate uncertainties on the correlation functions, we perform a bootstrap resampling of the data (e.g., Efron 1982; Norberg et al. 2009). We divide the sample footprint into 30 equally sized patches using the k-means clustering algorithm within the scikit-learn package (Pedregosa et al. 2011). We then randomly draw from these patches with replacement, and recalculate a correlation function with the data and random points within these patches. This process is repeated 500 times and the variance across realizations is taken to be the variance of our measurement. We find that this error estimation agrees well with analytic Poisson errors in all cases on these scales (e.g., Myers et al. 2006).

We also calculate the cross-power spectra, ${C}_{{\ell }}^{\kappa q}$, of the WISE quasar overdensity fields with the CMB lensing convergence field measured by the Planck satellite. First, a quasar overdensity map is produced at the same resolution as the lensing map by performing a weighted count of sources in each cell, where the weights correct for the angular systematics discussed in Section 3.2. We correct this density field for effects of subpixel masking by dividing the density by the fractional area-lost map. Finally, we convert this into a fractional overdensity map:

$$\begin{eqnarray}&&{\delta }_{q}=\displaystyle \frac{\rho -\langle \rho \rangle }{\langle \rho \rangle }.\end{eqnarray} \tag{ 3 }$$

We estimate the cross-spectrum between the two maps using the psuedo-C_ℓ (e.g., Peebles 1973) algorithm MASTER (Hivon et al. 2002), as implemented in the NaMaster package (Alonso et al. 2019). This algorithm allows for fast and nearly unbiased estimation of angular power spectra in the presence of sky masks. We measure the quasar-lensing cross-spectrum in 10 logarithmically spaced bins of angular multipole moment (ℓ) from 100 < ℓ < 1000. When fitting models to the observed spectrum, we bin the theoretical curves using the same mode-coupling matrix and bandpower binning scheme as used for the data, though in all figures we show the unbinned model spectra for simplicity.

To estimate uncertainties on the cross-spectra, we utilize simulated noise maps released as part of the Planck lensing data product. After generating 60 noise maps in the same manner as the data, we measure the cross-spectrum between the quasar overdensity field and each of the noise maps, taking the variance as our uncertainty estimate.

We model the correlation functions in a HOD framework (e.g., Berlind & Weinberg 2002). To first order, the HOD 〈N(M)〉 is the mean number of quasars belonging to halos of mass M, decomposed into contributions from quasars at the centers of halos, 〈N_c (M)〉, and secondary or "satellite" quasars belonging to the same halos, 〈N_s (M)〉:

$$\begin{eqnarray}&&\langle N\rangle =\langle {N}_{c}\rangle +\langle {N}_{s}\rangle .\end{eqnarray} \tag{ 4 }$$

The average number density of quasars is then an integral of the HOD over the halo mass function ${dn}/{dM}(M,z)$, for which we adopt the model of Tinker et al. (2008):

$$\begin{eqnarray}&&{n}_{q}(z)=\int {dM}\langle N\rangle \displaystyle \frac{{dn}}{{dM}}.\end{eqnarray} \tag{ 5 }$$

We adopt the HOD model developed by Zheng et al. (2007) and Zehavi et al. (2011), which models the occupation of central objects as a softened step function:

$$\begin{eqnarray}&&\langle {N}_{c}\rangle =\displaystyle \frac{1}{2}\left[1+\mathrm{erf}\left(\displaystyle \frac{{\mathrm{log}}_{10}(M/{M}_{\min })}{{\sigma }_{{\mathrm{log}}_{10}M}}\right)\right],\end{eqnarray} \tag{ 6 }$$

where "erf" is the error function, ${M}_{\min }$ is the characteristic minimum halo mass required to host a quasar, and ${\sigma }_{{\mathrm{log}}_{10}M}$ is the softening parameter. The satellite HOD is given by

$$\begin{eqnarray}&&\langle {N}_{s}\rangle ={\rm{\Theta }}(M-{M}_{0}){\left(\displaystyle \frac{M-{M}_{0}}{{M}_{1}}\right)}^{\alpha },\end{eqnarray} \tag{ 7 }$$

where Θ is the Heaviside step function, M₀ is the minimum mass to host a satellite quasar, and M₁ is the mass at which the term transitions to the power-law form. We note that a HOD model of this form has often been used in halo occupation studies of AGNs and quasars (e.g., Chatterjee et al. 2012; Richardson et al. 2012; Georgakakis et al. 2019; Alam et al. 2020), and thus we adopt this model for the sake of comparison to results from the literature.

This HOD model has five free parameters, which our present data will not be able to simultaneously constrain. We thus simplify the model as follows. We set ${M}_{0}={M}_{\min }$, such that the minimum halo mass required to host a satellite is the same as that required to host a central quasar. We fix the relatively unimportant parameter governing the softness of the step function to ${\sigma }_{{\mathrm{log}}_{10}M}$ = 0.4 dex. Finally, we fix ${M}_{1}=12\times {M}_{\min }$, motivated by the AGN HOD simulation results of Georgakakis et al. (2019) relying on the empirical accretion rate distributions of Georgakakis et al. (2017) and Aird et al. (2018) for luminous AGNs. We note that fixing M₁ to be the same for obscured and unobscured quasars may not be valid. However, M₁ and α are degenerate with one another, roughly driving the relative strength of the one-halo term and thus the fraction of quasars which are satellites. Breaking the degeneracy between α and M₁ is not possible with present data. Importantly, however, breaking this degeneracy in Equation (7) is not critical to our analysis or interpretation, as we are primarily concerned with deriving the satellite fraction via 〈N_s 〉. Thus, the two free parameters in our model are ${M}_{\min }$ and α. These parameters roughly govern the two-halo term and the one-halo term, respectively.

With the HOD parameters specified, the derived halo properties of interest can be expressed. These include the effective bias:

$$\begin{eqnarray}&&{b}_{\mathrm{eff}}=\int {dz}\displaystyle \frac{1}{{n}_{q}}\displaystyle \frac{{dN}}{{dz}}\int {dM}\displaystyle \frac{{dn}}{{dM}}\langle N\rangle {b}_{h}(M,z),\end{eqnarray} \tag{ 8 }$$

and the satellite fraction:

$$\begin{eqnarray}&&{f}_{\mathrm{sat}}=\int {dz}\displaystyle \frac{1}{{n}_{q}}\displaystyle \frac{{dN}}{{dz}}\int {dM}\displaystyle \frac{{dn}}{{dM}}\langle {N}_{s}\rangle .\end{eqnarray} \tag{ 9 }$$

The effective halo mass, M_eff, is computed by solving for the mass which would result in the effective bias (Equation (8)) when integrated over a sample's redshift distribution.

For a given HOD, the power spectrum of the quasar overdensity can be written as the sum of power contributed by pairs of quasars within the same halos (the one-halo term, P_1h (k, z)) and pairs of quasars between distinct halos (the two-halo term, P_2h (k, z)).

The one-halo term can be decomposed into pairs between satellites and pairs between satellites and central quasars. The satellites are prescribed to follow the density profile of halos. The function $\tilde{u}(k,M,z)$ is the Fourier transform of the dark matter density profile, for which we adopt the Navarro–Frenk–White (Navarro et al. 1997) model. We use the analytic solution of this profile's transform given by Scoccimarro et al. (2001):

$$\begin{eqnarray}&&{P}_{1h}=\displaystyle \frac{1}{{n}_{q}^{2}}\int {dM}\displaystyle \frac{{dn}}{{dM}}[2\tilde{u}\langle {N}_{c}{N}_{s}\rangle +{\tilde{u}}^{2}\langle {N}_{s}({N}_{s}-1)\rangle ].\end{eqnarray} \tag{ 10 }$$

We assume the satellites follow Poisson statistics such that 〈N_s (N_s − 1)〉 = 〈N_s 〉². Setting ${M}_{0}={M}_{\min }$ (see above) also imposes the "central condition" (Murray et al. 2021), which stipulates that a halo can host satellite quasars only if it hosts a central quasar, such that 〈N_c N_s 〉 = 〈N_s 〉.

The halo density profile transform, $\tilde{u}$, approaches unity as k → 0, such that the one-halo power spectrum approaches a constant at large scales, adding spurious power which can dominate over the two-halo term. We suppress this unphysical power by modifying the one-halo power spectrum with the ad hoc correction of Mead et al. (2021), with a characteristic damping scale k* = 10⁻² h Mpc⁻¹:

$$\begin{eqnarray}&&{P}_{1h}\to {P}_{1h}\times \displaystyle \frac{{\left(k/{k}^{* }\right)}^{4}}{1+{\left(k/{k}^{* }\right)}^{4}}.\end{eqnarray} \tag{ 11 }$$

The two-halo term is given by the matter power spectrum, P_m (k, z), multiplied by a bias factor:

$$\begin{eqnarray}&&{P}_{2h}={P}_{m}{\left[\int {dM}\displaystyle \frac{{dn}}{{dM}}{b}_{h}(M,z)\left(\langle {N}_{c}\rangle +\tilde{u}\langle {N}_{s}\rangle \right)\right]}^{2}.\end{eqnarray} \tag{ 12 }$$

A known problem in the halo model framework is the underprediction of power compared to N-body simulations in the "quasi-linear regime," the transition scales between which the one- and two-halo terms dominate (e.g., Fedeli et al. 2014; Mead et al. 2015). This arises due to a natural breakdown of linear perturbation theory as well as halo exclusion effects, though self-consistently modeling these effects is an active area of research. Here, we adopt the empirical function of Mead et al. (2015):

$$\begin{eqnarray}&&{P}_{{qq}}={\left[{\left({P}_{1h}\right)}^{\beta }+{\left({P}_{2h}\right)}^{\beta }\right]}^{1/\beta },\end{eqnarray} \tag{ 13 }$$

which remedies this effect by smoothing the power spectrum in the transition region. We adopt the best-fit value of β = 0.719 from Mead et al. (2021).

With a HOD power spectrum specified, we can compute an angular correlation function over the redshift distribution of our sample using the Limber (1953) approximation (e.g., Peebles 1980; Peacock 1991; DiPompeo et al. 2017a):

$$\begin{eqnarray}&&w(\theta )=\int {dz}\displaystyle \frac{H(z)}{c}{\left(\displaystyle \frac{{dN}}{{dz}}\right)}^{2}\int \displaystyle \frac{{dk}\ k}{2\pi }{P}_{{qq}}\ {J}_{0}[k\theta \chi (z)],\end{eqnarray} \tag{ 14 }$$

where J₀ is the zeroth-order Bessel function of the first kind. The k-space integral of Equation (14) is simply a Hankel transform of the power spectrum, for which we utilize the FFTLog algorithm (Talman 1978; Hamilton 2000).

Thus, given two free HOD parameters of ${M}_{\min }$ and α along with a source redshift distribution, a model angular correlation function can be calculated. We fit the observed correlation functions using a Markov Chain Monte Carlo (MCMC) method, as implemented in the emcee package (Foreman-Mackey et al. 2013).

We also model the cross-spectrum between matter overdensity and CMB lensing convergence. The quasar overdensity and lensing fields are both projections of three-dimensional density fields onto the plane of the sky, and thus the cross-spectrum is given by a line-of-sight integral of the matter power spectrum P(k, z) over the two respective projection kernels under the Limber (1953) approximation with the first-order correction (l → l + 1/2; LoVerde & Afshordi 2008):

$$\begin{eqnarray}&&{C}_{{\ell }}^{\kappa q}=c\int {dz}\ \displaystyle \frac{{W}^{\kappa }(z){W}^{q}(z)}{{\chi }^{2}(z)H(z)}P\left(k=\displaystyle \frac{{\ell }+1/2}{\chi (z)},z\right).\end{eqnarray} \tag{ 15 }$$

Here, c is the speed of light, χ(z) is the comoving distance, and H(z) is the Hubble parameter. We generate the matter power spectra using the analytic form of Eisenstein & Hu (1998). The CMB lensing kernel, a measure of the efficiency of lensing by structure (when multiplied by d χ/dz = c/H(z)) as a function of redshift is given by (e.g., Cooray & Hu 2000; Song et al. 2003)

$$\begin{eqnarray}&&{W}^{\kappa }(z)=\displaystyle \frac{3}{2}{{\rm{\Omega }}}_{m,0}{\left(\displaystyle \frac{{H}_{0}}{c}\right)}^{2}(1+z)\chi (z)\displaystyle \frac{{\chi }_{\mathrm{CMB}}-\chi (z)}{{\chi }_{\mathrm{CMB}}}.\end{eqnarray} \tag{ 16 }$$

The quasar overdensity kernel is in turn given by

$$\begin{eqnarray}&&{W}^{q}(z)={b}_{h}\displaystyle \frac{H(z)}{c}\displaystyle \frac{{dN}}{{dz}},\end{eqnarray} \tag{ 17 }$$

where dN/dz is the normalized redshift distribution of the lenses and b_h is the linear halo bias.

We have elected to fit the lensing spectra as linearly biased with respect to dark matter as we have found that our present data are unable to constrain the two-parameter HOD model introduced in the previous section. In particular, the Planck lensing map is noisiest at the small scales (large ℓ-modes) where the one-halo term dominates. Therefore, we instead fit the lensing spectra by assuming an effective halo mass as the single free parameter and inserting the mass–bias relation of Tinker et al. (2010) into Equation (17).

Inferring the host halo properties of a sample from angular statistics is inherently sensitive to systematics in the sample's redshift distribution, with the uncertainty in the latter oftentimes dominating over the statistical uncertainty of the clustering measurement (e.g., Coil 2013). As approximately half of our estimation of the obscured quasar redshift distribution relies on photometric redshifts (Figure 3), it is important to quantify whether our result could be explained by photometric redshift systematics, as was suggested by the authors of Mendez et al. (2016) to explain the similar results of Donoso et al. (2014) and DiPompeo et al. (2014).

We test this possibility by exploring the potential obscured quasar redshift distributions which would resolve the observed effective halo mass difference between obscured and unobscured quasars while reproducing our measurements. This was achieved by replacing the observed photometric redshift distribution of obscured quasars with many simulated distributions and refitting the data. First, the obscured quasar photometric redshift distribution of Figure 3 was fit as a normal distribution, finding a mean and dispersion of μ_obs = 1.4, σ_obs = 0.5. We then vary these parameters, shifting the mean by $-1\lt {\mu }_{\mathrm{sim}}-{\mu }_{\mathrm{obs}}\lt 3$ and varying the dispersion between $0.1\lt {\sigma }_{\mathrm{sim}}/{\sigma }_{\mathrm{obs}}\lt 1.5$. Redshifts are randomly drawn from this new distribution, and only values z > 0 are retained. Distributions with >10% of the cumulative distribution function below zero are rejected outright. These simulated redshifts are finally recombined with the observed spectroscopic redshifts. We then refit the obscured quasar clustering and lensing signals (on scales dominated by the two-halo term) and determine which combinations of the redshift distribution parameters imply an effective halo mass consistent with that measured for unobscured quasars. The combinations which resolve the tension are explored in an MCMC analysis, for which we display the results in Figure 10.

Zoom In Zoom Out Reset image size

It is clear that the dependence of clustering and CMB lensing measurements on redshift systematics are independent, with clustering being primarily sensitive to the distribution width and the lensing sensitive to the mean. We observe that the only possible configuration which would resolve the halo mass tension while reproducing both the clustering and lensing measurements would be if the obscured quasars without spectroscopic redshifts were distributed around z ≈ 3.5, which we regard as infeasible. The photometric redshifts stem mostly from Duncan et al. (2021), which reports a typical photometric redshift scatter of 7% and an outlier fraction of ∼20% for AGNs. In order to resolve the observed clustering differences, the true outlier fraction for obscured quasars would need to be nearly 100% if these systems were truly narrowly distributed at z ≈ 3.5. Therefore, we argue that our measurements of halo mass differences are robust against systematics from uncertainties in the redshift distribution. This test highlights the power of combining clustering and lensing measurements to mitigate redshift systematics.

A significant observed difference in the host halo properties of matched samples of obscured and unobscured quasars appears to negate the simplest unified model, which attributes obscuration solely to viewing angle. An interesting alternative in the context of galaxy evolution is an evolutionary model which posits that obscured and unobscured quasars represent different phases in a coevolutionary scheme between galaxies and their nuclear black holes. DiPompeo et al. (2017b) proposed a simple evolutionary model in which black hole growth lags behind dark matter halo growth (e.g., Alexander et al. 2008; Woo et al. 2008; Kormendy & Ho 2013). If obscured quasars represent an early evolutionary phase of rapid SMBH growth and thus their black holes are systematically undermassive with respect to their halos (they have not yet "caught up" to the SMBH halo relation), then the obscured subset of a luminosity-limited quasar sample (which is in turn black hole mass-limited for a fixed Eddington ratio distribution) will systematically occupy more massive halos.

However, it is also possible that such a difference could arise from nonevolutionary effects. The modeling work of Whalen et al. (2020) explored whether the observed clustering difference between obscured and unobscured quasars reported by DiPompeo et al. (2017a) could be explained by such nonevolutionary scaling relations. In particular, they studied whether the results could be understood within the "radiation-regulated unification" schema of Ricci et al. (2017b), where an observed relationship between AGN Eddington ratio and torus covering factor is interpreted to imply that the radiation pressure from black holes accreting closer to their limits can expel toroidal dust and expose more unobscured sight lines to a potential observer. Therefore, one would expect a reduced incidence of obscured quasars represented at high Eddington ratios, which has been confirmed at low redshift by Ananna et al. (2022). This implies that at fixed luminosity obscured quasar activity will preferentially be generated by more massive black holes accreting at lower Eddington ratios, which will thus lie in more massive halos according to black hole–halo mass relations (e.g., Kormendy & Ho 2013). Such a model could potentially explain the observed host halo difference without necessitating an evolutionary scheme.

Whalen et al. (2020) also studied whether the observed clustering difference could be explained by scaling relations between host-galaxy stellar mass and ISM column density (Pannella et al. 2009; Buchner et al. 2017; Whitaker et al. 2017). As more massive galaxies are observed to contain larger columns of obscuring gas, obscured quasars may be preferentially selected in massive galaxies and thus halos without regard to an evolutionary process.

The modeling work of Whalen et al. (2020) showed that the observed halo mass difference of DiPompeo et al. (2017a) cannot be explained either through radiation-regulated unification or through scaling relations with galaxy mass while simultaneously reproducing the observed fraction of obscured quasars. With the updated effective host halo masses presented in this work, which are consistent with but roughly twice as precise as the results of DiPompeo et al. (2017a), these arguments become even more stringent. We thus favor an evolutionary explanation of the host halo mass difference observed in this work.

However, accurate comparison of the clustering properties between populations of quasars requires controlling across properties such as redshift, quasar luminosity, and host-galaxy properties such as stellar mass and star formation rate. Indeed, several works suggest that AGN clustering properties depend only on the host-galaxy property distributions revealed by different AGN selection techniques (e.g., Mendez et al. 2016; Yang et al. 2018; Powell et al. 2018, 2020; Krishnan et al. 2020; Aird & Coil 2021). In this work, we have shown that WISE obscured and unobscured quasars exhibit similar redshift and bolometric luminosity distributions. However, a potential concern is that obscured quasars may appear to occupy more massive halos because of a selection effect in which they are detected in host galaxies differing from the hosts of unobscured systems. Future work on modeling the full spectral energy distributions of WISE quasars is required to test whether the observed clustering difference may be understood in terms of host-galaxy properties. However, Andonie et al. (2022) recently performed panchromatic SED modeling of IR-selected quasars in the COSMOS field, finding no difference in the stellar mass distribution for obscured and unobscured quasars.

Although a robust trend between quasar obscuration and halo properties rules out the possibility that obscuration is always driven by viewing angle with the torus, it does not preclude the role of nuclear obscuration. Instead, we expect a fraction of obscured quasars to appear obscured due to orientation effects. As the clustering of torus-obscured objects should match that of unobscured systems, the measured clustering of the obscured population in fact represents a lower limit for the population of systems obscured by evolutionary processes, as noted by DiPompeo et al. (2016b). The true clustering of this "non-torus-obscured" (NTO) population then depends on the NTO fraction, f_NTO, the proportion of obscured quasars which are obscured by evolutionary effects as opposed to orientation.

We thus investigate the halo properties required of NTO quasars in order to produce the observed clustering as a function of NTO fraction. We adopt the method of DiPompeo et al. (2016b) to compute the implied host halo properties for NTO quasars using our updated measurements. Assuming similar redshift distributions, the bias of the NTO population is given by

$$\begin{eqnarray}&&{b}_{\mathrm{NTO}}=\displaystyle \frac{{b}_{\mathrm{ob}}-{b}_{\mathrm{unob}}}{{f}_{\mathrm{NTO}}}+{b}_{\mathrm{unob}}.\end{eqnarray} \tag{ 18 }$$

Therefore, the smaller the fraction of obscured quasars belonging to the NTO class, the more biased and rare halos NTO quasars must occupy to drive the observed clustering enhancement of all obscured quasars. For a physical interpretation, we use this NTO bias to estimate the minimum host halo mass, halo occupation fraction, and characteristic lifetime of NTO quasars as a function of NTO fraction. We convert between this NTO bias and a minimum halo mass required to host a quasar as a function of the NTO fraction f_NTO using

$$\begin{eqnarray}&&b(M\gt {M}_{\min })=\displaystyle \frac{{\int }_{{M}_{\min }}^{\infty }{dM}\tfrac{{dn}}{{dM}}b(M)}{{\int }_{{M}_{\min }}^{\infty }{dM}\tfrac{{dn}}{{dM}}}.\end{eqnarray} \tag{ 19 }$$

We note that this minimum mass differs slightly from that defined in our HOD modeling, which includes a softening of the step function. Given a minimum halo mass, we are able to compute an occupation fraction, f_occ, the number fraction of halos more massive than a threshold mass which host a given type of quasar, by comparing the space density of quasars with the density of halos more massive than the minimum mass. To estimate the space density of our quasar samples, we integrate the Shen et al. (2020) quasar luminosity function over the observed 6 μm luminosity distributions (Section 2.5), also integrating over the sample's redshift distribution. Next, we divide the resulting density by 2 in order to isolate the density of the obscured population (which comprises roughly half of the total; see Section 2). For NTO quasars we finally multiply the density by the NTO fraction. Next, we compute the space density of halos as a function of minimum mass by performing mass and redshift integrals over the halo mass function of Tinker et al. (2008).

Finally, we provide estimates of the quasar lifetime, the characteristic time a SMBH is expected to be observable as a quasar of a given phase. Assuming that every halo hosts one SMBH, the quasar lifetime is simply proportional to the halo occupation fraction (Haiman & Hui 2001; Martini & Weinberg 2001). A majority of our sample quasars lie at z ≈ 0.5–2, corresponding to a span of ∼5 Gyr of cosmic time. An approximate estimate of the quasar lifetime is thus given by this factor multiplied by the occupation fraction.

We show the results of these calculations in Figure 11. It is apparent that if NTO quasars make up a small proportion of the obscured population, their host halos must be very massive in order to drive the clustering of the full obscured population to the observed value. Conversely, we can use this relation to put constraints on the minimum NTO fraction required by our measurements in order to avoid unphysically large host halo masses for these objects. At sufficiently small NTO fractions (≲15%), the implied host halos of NTO quasars become massive and thus rare enough that the corresponding occupation fraction exceeds unity. This forbidden region is shown in Figure 11 with a gray hatched region. This appears to be a compelling line of evidence that a nonnegligible (≳15%) fraction of obscured quasars are obscured by a mechanism aside from orientation, such as galactic or circumnuclear dust. We also note that our quasar lifetime analysis suggests that if NTO quasars represent a distinct evolutionary phase in the evolution of SMBH growth, the duration of the obscured phase appears greater than the unobscured phase by a factor ≳3.

Zoom In Zoom Out Reset image size

By matching our quasar catalog to SDWFS 24 μm photometry, we find that ∼30% of the objects in the the obscured quasar sample would be classified as dust-obscured galaxies (DOGs; Dey et al. 2008) with F₂₄/F_r > 1000. DOGs are ultraluminous IR galaxies often accompanied by quasar activity at z = 1–2, where the optical emission from both stars and SMBH accretion appears heavily attenuated, implying strong galactic-scale absorption. These objects thus are expected to be good candidates for galaxies caught in the violent postmerger evolutionary stage expected to produce obscured quasar activity and star formation in the models of Sanders et al. (1988) and Hopkins et al. (2008). The DOGs in our sample are expected to be quasar-dominated given that their IR colors satisfy the Donley et al. (2012) criterion of power-law IR spectra. We observe that the DOG fraction increases in our obscured quasar sample toward systems with redder r−W2 colors. To test whether the enhanced clustering of obscured quasars observed in this work could be driven by the reddest tail of sources with significant overlap with the DOG population, we split the obscured sample into two further subsets and again perform each clustering and lensing analysis as previously presented. We split the obscured population in two using a cut of r−W2 = 4.5, which corresponds to the peak of the obscured quasar distribution (Figure 1). The subset redder than this cut consists of ∼50% DOGs. We do not split the unobscured sample, as Petter et al. (2022) showed that unobscured spectroscopic quasars' clustering is not connected with their r−W2 colors. We display the effective halo masses from these analyses as a function of r−W2 color in Figure 12.

Zoom In Zoom Out Reset image size

Interestingly, the effective halo mass of WISE quasars appears to continue increasing toward more extreme red optical-IR colors. We find that the results from clustering and lensing again agree, which indicates this trend is robust. The reddest sample appears to occupy halos of effective mass ∼10^13.0 h⁻¹ M_⊙. Meanwhile, the study of Toba et al. (2017) measured the clustering of WISE-selected bright DOGs (which tend to be quasar-dominated), finding that they occupy very massive effective halos of ∼10^13.6 h⁻¹ M_⊙. Given a majority of the reddest sample would be selected as DOGs, we speculate that the observed enhanced clustering of obscured IR-selected quasars may be driven by the reddest and most heavily obscured subset, which are often DOGs and expected to be obscured by host-galaxy material. However, we note that we perform a full HOD fit of each subsample but do not detect an excess of small-scale clustering in the reddest bin, as might be expected for a sample dominated by merger systems.

It is important to understand the origin of the wide array of seemingly conflicting results on the clustering of obscured and unobscured AGNs in the literature. Assessing the literature, we broadly find that studies detecting enhanced clustering of obscured quasars tend to be selected in the IR (Hickox et al. 2011; Donoso et al. 2014; DiPompeo et al. 2014, 2015,2016a, 2017a; this work) or ultrahard X-ray (Powell et al. 2018), while studies which do not detect a difference tend to use softer X-ray selection (e.g., with Chandra or XMM-Newton) over relatively smaller fields (Coil et al. 2009; Gilli et al. 2009; Ebrero et al. 2009; Mountrichas & Georgakakis 2012; Koutoulidis et al. 2018), though there are exceptions (Geach et al. 2013; Mendez et al. 2016; Krumpe et al. 2018).

We posit that the diversity of results may be understood by one of two selection effects. In one scenario, quasar selection at certain wavelengths may reveal obscured and unobscured populations which are not matched in host-galaxy properties that correlate with dark matter halo mass. In this scenario, it is possible, for instance, that IR quasar selection may uncover obscured quasars that occupy more massive galaxies than unobscured quasars, which would produce a difference in observed halo mass given that galaxy stellar mass correlates with halo mass. However, the modeling work of Whalen et al. (2020) appears to point against this scenario.

Alternatively, it is possible that the effective halo mass of a sample of quasars is driven by an extreme subset of that sample which are not equally recoverable in all wave bands. In this context, we suggest that the enhanced clustering of obscured quasars may be driven by the most heavily obscured subset, which are able to be selected at IR and hard X-ray wavelengths but not with soft X-rays. We have shown in Figure 12 that the clustering of obscured quasars continues to increase toward redder optical-IR colors and greater overlap with the DOG population, possibly supporting this view. Quasar-dominated DOGs have been shown to often be obscured by Compton-thin or Compton-thick columns (e.g., Fiore et al. 2008; Lanzuisi et al. 2009; Stern et al. 2014; Piconcelli et al. 2015; Assef et al. 2016; Del Moro et al. 2016; Vito et al. 2018), which would prevent their detection in typical soft X-ray surveys. We thus speculate that the seemingly conflicting results in the literature regarding the host halos of obscured and unobscured quasars may be explained if the enhanced clustering of IR-selected obscured quasars is driven by deeply buried systems missed by soft X-ray surveys.