# Host Dark Matter Halos of WISE-selected Obscured & Unobscured Quasars:
Evidence for Evolution

Grayson C. Petter Department of Physics and Astronomy, Dartmouth College, 6127
Wilder Laboratory, Hanover, NH 03755, USA  Ryan C. Hickox Department of
Physics and Astronomy, Dartmouth College, 6127 Wilder Laboratory, Hanover, NH
03755, USA  David M. Alexander Centre for Extragalactic Astronomy, Department
of Physics, Durham University, South Road, Durham, DH1 3LE, UK  Adam D. Myers
Department of Physics and Astronomy, University of Wyoming, Laramie, WY 82071,
USA  James E. Geach Centre for Astrophysics Research, School of Physics,
Astronomy & Mathematics, University of Hertfordshire, Hatfield, AL10 9AB, UK
Kelly E. Whalen Department of Physics and Astronomy, Dartmouth College, 6127
Wilder Laboratory, Hanover, NH 03755, USA  Carolina P. Andonie Centre for
Extragalactic Astronomy, Department of Physics, Durham University, South Road,
Durham, DH1 3LE, UK

(Accepted 2023 January 30 to ApJ)

###### Abstract

Obscuration in quasars may arise from steep viewing angles along the dusty
torus, or instead may represent a distinct phase of supermassive black hole
growth. We test these scenarios by probing the host dark matter halo
environments of $\sim 1.4$ million WISE-selected obscured and unobscured
quasars using angular clustering measurements as well as cross-correlation
measurements of quasar positions with the gravitational lensing of the cosmic
microwave background (CMB). We interpret these signals within a halo
occupation distribution (HOD) framework to conclude that obscured systems
reside in more massive effective halos ($\sim 10^{12.9}h^{-1}M_{\odot}$) than
their unobscured counterparts ($\sim 10^{12.6}h^{-1}M_{\odot}$), though we do
not detect a difference in the satellite fraction. We find excellent agreement
between the clustering and lensing analyses and show that this implies the
observed difference is robust to uncertainties in the obscured quasar redshift
distribution, highlighting the power of combining angular clustering and weak
lensing measurements. This finding appears in tension with models that ascribe
obscuration exclusively to orientation of the dusty torus along the line-of-
sight, and instead may be consistent with the notion that some obscured
quasars are attenuated by galaxy-scale or circumnuclear material during an
evolutionary phase.

Quasars — Gravitational Lensing

††facilities: WISE, DESI, IRSA, Planck, Spitzer††software: astropy: Astropy
Collaboration et al. (2018), CAMB: Lewis et al. (2000), colossus: Diemer
(2018), Corrfunc: Sinha & Garrison (2020), emcee: Foreman-Mackey et al.
(2013), healpy/HEALPix: Zonca et al. (2019); Górski et al. (2005)

##  1 Introduction

Quasars are the most luminous class of Active Galactic Nuclei (AGNs),
manifestations of accretion onto supermassive black holes (SMBHs) at the
centers of galaxies. Since their discovery (Schmidt, 1963), AGNs and quasars
have been classified into a taxonomy according to their observed
multiwavelength properties (Padovani et al., 2017). However, it remains
unclear whether or not many of the observed differences between quasar sub-
populations reflect intrinsic or cosmologically relevant features of black
hole growth.

The class of quasars which are “obscured” has been increasingly recognized as
important to characterize. These systems often lack the optical and soft X-ray
signatures of unobscured quasars, implying that significant columns of dust
and gas lie between the sub-pc region surrounding the black hole and our line
of sight. A majority of the AGN activity in the Universe is obscured (e.g.,
Ueda et al., 2014; Aird et al., 2015; Buchner et al., 2015; Ananna et al.,
2019), which implies that traditional AGN surveys in the optical or X-rays
have potentially led to a biased picture of SMBH growth. Some of the most
pressing questions in extragalactic astrophysics require a complete census of
AGN activity extending to even the most heavily obscured systems, including
understanding the physical structure of AGN systems, probing the cosmic
history of black hole growth, and uncovering the role that AGNs play in
galactic evolution (Hickox & Alexander, 2018). Crucial to these endeavors is
understanding the nature of the obscuring material and whether obscured
quasars represent intrinsically different systems to their unobscured
counterparts.

The canonical view of the AGN structure is termed the “unified model”
(Antonucci, 1993; Urry & Padovani, 1995; Elvis, 2000; Netzer, 2015; Ramos
Almeida & Ricci, 2017), which posits that most AGN systems consist of a
similar axisymmetric structure of material surrounding the SMBH, but appear to
differ because of our chance line-of-sight toward this structure. This model
implies that a quasar will appear obscured when our viewing angle happens to
nearly coincide with the plane of the dusty “torus”. In this view, obscured
and unobscured quasars represent intrinsically similar objects and are thus
expected to occupy similar environments, even when accounting for potential
relationships between accretion rate and torus covering factor (Whalen et al.,
2020).

Alternatively, evolutionary models suggest that AGNs vary in their observed
properties over the course of their lifetimes through interaction with their
broader environments. In this scenario, obscured quasar activity could
represent a specific evolutionary stage of SMBH growth. One such model posits
that rapid star formation and SMBH growth are triggered during major galaxy
mergers which funnel gas to the nuclear region of the merger remnant. This
process is expected to generate quasar activity obscured by galactic or
circumnuclear star-forming gas, before this activity heats and expels the gas
to unveil an unobscured quasar (e.g., Sanders et al., 1988; Hopkins et al.,
2005, 2006, 2008; Alexander & Hickox, 2012; Hickox & Alexander, 2018). A
number of recent studies have found evidence for AGN obscuration taking place
on galactic scales or correlated with star formation activity (e.g., Chen et
al., 2015; Buchner et al., 2017; Ricci et al., 2017a; Circosta et al., 2019;
Li et al., 2021; Yan et al., 2021; Andonie et al., 2022; Gilli et al., 2022;
Juneau et al., 2022), which may be connected to this evolutionary scenario.

A natural test of these models lies in measuring how quasars populate the
large-scale structure (LSS) of the Universe, or the manner in which they
occupy dark matter halos. This is because the properties of a quasar-hosting
halo could not feasibly be connected with the obscuring torus’ orientation
along our particular line-of-sight, while a connection between halo properties
and galaxy evolution is expected. Thus, this work aims to probe the nature of
infrared-selected obscured and unobscured quasars by estimating the host dark
matter halo properties of each class.

The connection between halo properties and obscuration is contested. Some
studies find that obscured systems occupy more massive halos (Hickox et al.,
2011; Elyiv et al., 2012; Donoso et al., 2014; DiPompeo et al., 2014, 2015,
2016a, 2017a; Powell et al., 2018), some find the opposite (Cappelluti et al.,
2010; Allevato et al., 2014), and others find no trend (Coil et al., 2009;
Gilli et al., 2009; Ebrero et al., 2009; Mountrichas & Georgakakis, 2012;
Geach et al., 2013; Mendez et al., 2016; Jiang et al., 2016; Koutoulidis et
al., 2018; Krumpe et al., 2018). These studies, though, vary in their
selection, statistical power, analysis technique, obscuration definition, and
sample distributions across luminosity and redshift.

In this work, we utilize the largest sample of mid-infrared-selected quasars
to date in order to put precise constraints on the host halo properties of
obscured and unobscured quasars . We probe these properties first by measuring
angular autocorrelation functions, and interpreting these signals in a halo
occupation distribution (HOD) framework. We find that obscured quasars occupy
significantly more massive halos than unobscured quasars on average, but we do
not detect a difference in the fraction that are satellites. Next, we test the
halo properties with an entirely independent method, the cross-correlation of
quasar positions with Planck’s map of the gravitational lensing of the CMB. We
interpret the lensing signals with a linearly-biased model, and again find
that obscured quasars occupy significantly more massive halos. The implied
effective halo masses from the clustering and lensing analyses are in
excellent agreement, which we show implies our results are robust against
uncertainties in the obscured quasar redshift distribution. We interpret these
results as favoring an evolutionary explanation for the obscuration of at
least some quasars.

Throughout this work, we adopt a “Planck 2018” CMB+BAO $\Lambda$-CDM
concordance cosmology (Planck Collaboration et al., 2020a), with $h=H_{0}/100\
\mathrm{km\ s}^{-1}\mathrm{Mpc}^{-1}=0.6766$, $\Omega_{m}=0.3111$,
$\Omega_{\Lambda}=0.6888$, $\sigma_{8}=0.8102$, and $n_{s}=0.9665$.

##  2 Data

###  2.1 WISE Quasar Sample

Though heavily obscured quasars can be challenging to distinguish from normal
galaxies in the optical waveband, they can be recovered simply via their red
mid-infrared 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-infrared bands centered at (named) 3.6 (W1), 4.5 (W2), 12
(W3), and 22 $\mu$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., 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
IRSA111https://irsa.ipac.caltech.edu database and select all objects in the
CatWISE catalog 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
infrared stars:

$\left\\{\begin{aligned}
&\mathrm{W}1-\mathrm{W}2>0.65\times\mathrm{exp}[0.153\times(\mathrm{W}2-13.86)^{2}],\\\
&9<\mathrm{W}2<16.5\end{aligned}\right.$ (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 multi-order coverage (MOC) 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),
$\mathrm{H}_{\mathrm{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 (LVG; Karachentsev et al., 2013) or the
2MASS extended source catalog (XSC; 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 $>0.25\arcsec/\mathrm{yr}$, which are able to be
detected at $>5\sigma$ 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-infrared colors of infrared-selected
quasars appears bimodal, as rest-frame UV-optical emission is easily extincted
by dust, while near-infrared emission is less so. Therefore, a simple color
cut can be used to classify obscured systems. Quasars with optical-infrared
colors $r-$W2 $>3.1$ [AB] typically show X-ray absorption corresponding to
absorbing column densities of $N_{\mathrm{H}}>10^{22}$ (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 $\sim 20,000\ \mathrm{deg}^{2}$ of extragalactic sky to
$r_{\mathrm{AB}}\sim 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 greater 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\arcsec$ (Assef et al., 2013) using NOIRLab’s Astro Data
Lab222https://datalab.noirlab.edu tools. The optical photometry is corrected
for Galactic reddening using the map of Schlegel et al. (1998). We retain the
$\sim 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 non-detections to the Hyper
Suprime-Cam Subaru Strategic Program (HSC-SSP) “DEEP” catalog (Aihara et al.,
2018) in the areas where the surveys overlap. 96$\%$ of the DESI-LS non-
detections are detected in HSC-DEEP, implying a negligible contamination rate
of spurious WISE sources in our sample.

The resulting optical-infrared color distribution of our quasar sample is
shown in Figure 1. It is clear that WISE mid-infrared selection reveals a
population with a bimodal optical-infrared color distribution, and thus
uncovers both obscured and unobscured sources in roughly equal proportion. We
adopt the criterion of Hickox et al. (2007) ($r-\mathrm{W}2$ [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.

Figure 1: The optical-infrared ($r-$W2) color distribution of 1.4 million
WISE-selected quasar candidates. The distribution is clearly bimodal,
reflecting that infrared quasar selection recovers both obscured and
unobscured populations. We classify quasars with $r-$W2 ¿ 3.1 [AB] as
obscured, displaying them in red, and show unobscured quasars in blue. We note
that all masking described in $\S$2.1 and contaminant removal ($\S$2.2) has
been applied before producing this distribution. We also show the color
distribution of the $8\%$ of sources not detected in DESI-LS using black
hatched boxes, obtained by matching against the HSC-SSP DEEP catalog.

###  2.2 Sample Purity

In order to robustly interpret clustering measurements of quasar candidates,
it is important to optimize and quantify the purity of the quasar selection.
Mid-infrared 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-0.3$ galaxies with high
specific star formation rates. Messias et al. (2012) developed a criterion to
discriminate between low redshift star-forming galaxies and AGN using a near-
infrared to mid-infrared 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
(REX) than as a point source
($\chi^{2}_{\mathrm{REX}}-\chi^{2}_{\mathrm{PSF}}>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.

Figure 2: The distribution of Assef et al. (2018) R90 infrared quasar
candidates in an optical ($r$), near-infrared ($z$), mid-infrared (W2) color
space. We show the density of candidates in greyscale, and in green we show
the density of highly resolved sources in the optical band representing low-
redshift galaxy contaminants. Unobscured candidates ($r$-W2 ¡ 3.1) appear
bimodal roughly along a line of constant $r-z$, which differentiates low-
redshift SFG contaminants from quasars. Template colors of Type-1 quasars
(Hickox et al., 2017) at $0<z<2$ and the starburst M82 at $0<z<0.7$ (Polletta
et al., 2007) are shown with blue open circles and green stars, respectively,
and increasing marker size represents linear steps towards increasing
redshift. We show our criterion to remove low-$z$ star-forming galaxy
contaminants with a black dashed line.

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-infrared
Spitzer-IRAC colors indicative of AGN activity (Donley et al., 2012), or
coincides with a luminous X-ray source (Kenter et al., 2005). 95% 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 non-quasar
contamination will be subdominant.

###  2.3 Redshift Information

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., 2022;
Chaussidon et al., 2022), 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\arcsec$ 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), SDSS DR16
(Ahumada et al., 2020), zCOSMOS (Scoville et al., 2007), PRIMUS (Coil et al.,
2011; Cool et al., 2013), IDEOS (Hernán-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. . 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 $\S$6.

Figure 3: The redshift distributions $(dN/dz)$ of the unobscured (top panel in
blue) and obscured (bottom panel in red) quasar samples obtained by cross-
matching to surveys in the Boötes and COSMOS fields. The samples are similarly
distributed in redshift , implying that we can appropriately compare their
clustering properties. Redshift distribution information stemming from
spectroscopic redshifts is shown with a dashed line, a dotted line shows
information from photometric redshifts, and a solid line shows the
combinations. The unobscured distribution is derived mostly from spectroscopic
redshifts, while we rely on photometric redshifts for half of the obscured
population. We also display the shape of the CMB lensing kernel (Eq. 16) with
a dashed black line, showing that the LSS surrounding our samples should
efficiently lens CMB photons.

###  2.4 Planck Lensing Convergence Map

The CMB radiation from $z\approx 1090$ has been gravitationally lensed by the
intervening large-scale structure (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 (MV)
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 $\kappa_{\ell m}$ coefficients at
$\ell<2048$ into NSIDE=1024 resolution HEALPix maps.

###  2.5 Ancillary Photometry

Finally, we estimate the bolometric luminosity distributions of obscured and
unobscured quasars to compute space densities and occupation statistics in
$\S$6.3. We use the rest-frame $6\ \mu$m band as a bolometric luminosity
tracer which is minimally affected by obscuration. As $90\%$ 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. For checking the location of these sources in Spitzer color-color
space, we utilize mid-infrared 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 (HELP, Shirley et al., 2019, 2021). .

Figure 4:

##  3 Measurements

###  3.1 Masking

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.,
2022)333https://portal.nersc.gov/cfs/cosmo/data/legacysurvey/dr9/randoms/
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-$\sigma$ 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 3.4$\arcmin$ 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
($\S$ 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$ degree)
diffraction spikes. As our candidate quasars are selected on their very red
mid-infrared colors, these regions appear to be affected by the reddest and
brightest infrared 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 W$3\ [\mathrm{Vega}]<-1.5$ (of
which there are only $\sim 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 apply a mask within 10∘ of each Ecliptic pole.

Figure 5:

We apply all 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.

###  3.2 Weighting Scheme

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.

###  3.3 Angular Correlation Functions

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(\theta)$ is defined as the excess probability —
above that expected from an unclustered Poisson distribution — of finding a
pair of objects at angular separation $\theta$ (Peebles, 1980). We estimate
the angular correlation function in this work using the Landy & Szalay (1993)
estimator:

$w(\theta)=\frac{DD(\theta)-2DR(\theta)+RR(\theta)}{DD(\theta)},$ (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}<\theta<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
$\sim 6\arcsec$ 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 $\gtrsim 1^{\circ}$ 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., 2012).
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).

###  3.4 Lensing Cross-Correlations

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
$\S$3.2. We correct this density field for effects of sub-pixel masking by
dividing the density by the fractional-area-lost map. Finally, we convert this
into a fractional overdensity map:

$\delta_{q}=\frac{\rho-\langle\rho\rangle}{\langle\rho\rangle}$ (3)

We estimate the cross spectrum between the two maps using the
psuedo-$C_{\ell}$ (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 ($\ell$) from
$100<\ell<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.

##  4 Modeling

###  4.1 Modeling Correlation Functions

We model the correlation functions in a halo occupation distribution (HOD)
framework (e.g., Berlind & Weinberg, 2002). To first order, the HOD $\langle
N(M)\rangle$ is the mean number of quasars belonging to halos of mass $M$,
decomposed into contributions from quasars at the centers of halos, $\langle
N_{c}(M)\rangle$, and secondary or ‘satellite” quasars belonging to the same
halos $\langle N_{s}(M)\rangle$:

$\langle N\rangle=\langle N_{c}\rangle+\langle N_{s}\rangle$ (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):

$n_{q}(z)=\int dM\langle N\rangle\frac{dn}{dM}$ (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:

$\langle
N_{c}\rangle=\frac{1}{2}\left[1+\mathrm{erf}\left(\frac{\mathrm{log}_{10}(M/M_{\mathrm{min}})}{\sigma_{\mathrm{log}_{10}M}}\right)\right]$
(6)

where erf is the error function, $M_{\mathrm{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:

$\langle
N_{s}\rangle=\Theta(M-M_{0})\left(\frac{M-M_{0}}{M_{1}}\right)^{\alpha}$ (7)

where $\Theta$ is the Heaviside step function, $M_{0}$ is the minimum mass to
host a satellite quasar, and $M_{1}$ is the mass at which the term transitions
to the power-law form.

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_{\mathrm{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_{\mathrm{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 AGN. We note
that fixing $M_{1}$ to be the same for obscured and unobscured quasars may not
be valid. However, $M_{1}$ and $\alpha$ 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
$\alpha$ and $M_{1}$ is not possible with present data. Importantly however,
breaking this degeneracy in Eq. 7 is not critical to our analysis or
interpretation, as we are primarily concerned with deriving the satellite
fraction via $\langle N_{s}\rangle$. Thus, the two free parameters in our
model are $M_{\mathrm{min}}$ and $\alpha$. These parameters roughly govern the
two-halo term and the one-halo term, respectively.

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

$b_{\mathrm{eff}}=\int dz\frac{1}{n_{q}}\frac{dN}{dz}\int
dM\frac{dn}{dM}\langle N\rangle b_{h}(M,z)$ (8)

and the satellite fraction:

$f_{\mathrm{sat}}=\int dz\frac{1}{n_{q}}\frac{dN}{dz}\int
dM\frac{dn}{dM}\langle N_{s}\rangle$ (9)

The effective halo mass $M_{\mathrm{eff}}$ is computed by solving for the mass
which would result in the effective bias (Eq. 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
“NFW” (Navarro et al., 1997) model. We use the analytic solution of this
profile’s transform given by Scoccimarro et al. (2001):

$P_{1h}=\frac{1}{n_{q}^{2}}\int dM\frac{dn}{dM}[2\tilde{u}\langle
N_{c}N_{s}\rangle+\tilde{u}^{2}\langle N_{s}(N_{s}-1)\rangle]$ (10)

We assume the satellites follow Poisson statistics such that $\langle
N_{s}(N_{s}-1)\rangle=\langle N_{s}\rangle^{2}$. Setting
$M_{0}=M_{\mathrm{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 $\langle
N_{c}N_{s}\rangle=\langle N_{s}\rangle$.

The halo density profile transform $\tilde{u}$ approaches unity as
$k\rightarrow 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^{-2}\ h$/Mpc:

$P_{1h}\rightarrow P_{1h}\times\frac{(k/k^{*})^{4}}{1+(k/k^{*})^{4}}$ (11)

The two-halo term is given by the matter power spectrum $P_{m}(k,z)$
multiplied by a bias factor:

$P_{2h}=P_{m}\left[\int dM\frac{dn}{dM}b_{h}(M,z)\left(\langle
N_{c}\rangle+\tilde{u}\langle N_{s}\rangle\right)\right]^{2}$ (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):

$P_{qq}=\left[(P_{1h})^{\beta}+(P_{2h})^{\beta}\right]^{1/\beta},$ (13)

which remedies this effect by smoothing the power spectrum in the transition
region. We adopt the best-fit value of $\beta=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):

$w(\theta)=\int dz\frac{H(z)}{c}\left(\frac{dN}{dz}\right)^{2}\int\frac{dk\
k}{2\pi}P_{qq}\ J_{0}[k\theta\chi(z)]$ (14)

where $J_{0}$ is the zeroth-order Bessel function of the first kind. The
$k$-space integral of Eq. 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_{\mathrm{min}}$ and $\alpha$ along
with a source redshift distribution, a model angular correlation function may
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).

###  4.2 Modeling the Lensing Cross-Spectrum

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\rightarrow l+1/2$;
LoVerde & Afshordi, 2008):

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

Here, $c$ is the speed of light, $\chi(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\chi/dz=c/H(z)$)
as a function of redshift is given by (e.g., Cooray & Hu, 2000; Song et al.,
2003):

$W^{\kappa}(z)=\frac{3}{2}\Omega_{m,0}\left(\frac{H_{0}}{c}\right)^{2}(1+z)\chi(z)\frac{\chi_{\mathrm{CMB}}-\chi(z)}{\chi_{\mathrm{CMB}}}$
(16)

The quasar overdensity kernel is in turn given by:

$W^{q}(z)=b_{h}\frac{H(z)}{c}\frac{dN}{dz}$ (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 $\ell$-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 Eq. 17.

##  5 Results

The results of the angular autocorrelation measurements are displayed in
Figure 6. It is apparent that obscured quasars cluster significantly more
strongly than their unobscured counterparts (across a similar redshift range),
implying that obscured quasars are more biased tracers of matter and occupy
more massive dark matter halos. The data for both samples are very well fit by
our two parameter HOD model. However, the relative strength of the one-halo
term is similar for the two populations, implying similar satellite fractions.

Figure 6: The angular autocorrelation functions of obscured and unobscured
quasars are shown with red and blue markers, respectively. The model dark
matter correlation functions for the corresponding redshift distributions are
shown with solid lines. Finally, the best HOD model fits are shown with
dashed/dotted lines. The dotted line indicates where the one-halo term
dominates, while the dashed line indicates two-halo term domination. Obscured
quasars cluster more strongly than their unobscured counterparts, implying
they are a more biased tracer of matter and occupy more massive halos. Figure
7: The posterior distributions of our two-parameter HOD model from MCMC fits
to the angular autocorrelation functions of obscured and unobscured quasars,
shown with red and blue contours, respectively. It is clear that the host halo
properties differ significantly, with obscured quasars requiring a larger
minimum host halo mass, implying that the effective mass and bias are also
larger. However, we do not detect a difference in the one-halo term index
$\alpha$, or the satellite fraction $f_{\mathrm{sat}}$.

The posterior distributions of the HOD parameters from the MCMC fits to the
correlation functions are displayed in Figure 7. It is clear that obscured
quasars occupy their host halos in a signifcantly different manner than
unobscured quasars. The minimum mass required to host an obscured quasar is
higher, implying both larger effective biases and host halo masses. In
particular, obscured quasars appear to occupy effective halos of $\sim
10^{12.9}h^{-1}M_{\odot}$, while unobscured quasars occupy halos of $\sim
10^{12.6}h^{-1}M_{\odot}$. However, we do not detect a difference in the
satellite power law index $\alpha$, nor the derived satellite fraction,
finding that $\sim 5-20\%$ of obscured and unobscured quasars are satellites
within their halos. A non-detection of a difference in the satellite fraction
is in contrast to the study of Mitra et al. (2018). We thus find from angular
clustering measurements that obscured quasars occupy more massive halos than
unobscured systems, but are not more likely to be satellites within their
halos.

We also investigate the host halo properties of quasars as a function of
obscuration using an independent technique, by calculating cross-correlations
of quasar overdensities with the CMB lensing convergence map measured by
Planck. We display the result of this in Figure 8. In agreement with the
clustering analysis, we find that obscured quasars are a significantly more
biased tracer of dark matter, again implying more massive host halo
environments. We are unable to constrain the one-halo term, and thus satellite
fraction with present data, as the Planck CMB lensing map is noisiest at small
scales, as discussed in $\S$4.2. Therefore, we fit the observed spectra with
linearly-biased dark matter models, which provide good fits to the data.

Figure 8: The cross-power spectra (multiplied by $\ell$ to reduce dynamic
range) between obscured/unobscured quasar overdensity and Planck CMB lensing
convergence $\kappa$ are shown with red and blue markers. The model spectra of
linear dark matter for the corresponding redshift distributions are shown
using solid lines, while the best model fits are shown with dashed lines.
Obscured quasar density correlates more strongly with lensing convergence,
implying that obscured quasars are more biased tracers of matter and occupy
more massive halos.

We summarize our results in Figure 9 and tabulate them in Table 1. This shows
most strikingly that the implied effective halo masses from clustering and CMB
lensing analyses are in excellent agreement, with both showing that obscured
quasars occupy significantly more massive effective halos. We show that
infrared-selected unobscured quasars occupy similar effective halos as
optically-selected spectroscopic Type-1 quasars from eBOSS (e.g., Laurent et
al., 2017). We also show that our results are consistent with those of
DiPompeo et al. (2017a), demonstrating that the enhanced clustering of WISE-
selected obscured quasars persists both with improved precision and at higher
redshifts.

Figure 9: The bias (top panel) and effective halo mass (bottom panel) of obscured and unobscured quasars. We show the results for obscured quasars in red, and unobscured quasars in blue. Our angular clustering results are shown with closed circles, while the CMB lensing results are displayed with open circles, offset horizontally for visual clarity. We also show clustering measurement results of infrared quasars at similar effective redshifts from (Hickox et al., 2011) and (DiPompeo et al., 2017a), as well as for SDSS Type-1 quasars (Laurent et al., 2017).  Analysis | Sample | $b_{q}$ | log${}_{10}(M_{h}/h^{-1}M_{\odot}$
---|---|---|---
Clustering |  |  |
| Unobscured | $2.27\pm 0.06$ | $12.58\pm_{0.05}^{0.04}$
| Obscured | $2.96\pm 0.07$ | $12.93\pm 0.03$
Lensing |  |  |
| Unobscured | $2.29\pm 0.18$ | $12.60\pm 0.09$
| Obscured | $3.02\pm 0.19$ | $12.96\pm 0.06$
Table 1: A tabulation of the measured halo biases and the corresponding
effective halo masses for obscured and unobscured quasars. We display results
from angular clustering analyses as well as from CMB lensing, which are in
excellent agreement and show obscured quasars occupy more massive halos. The
contents of this table are visually represented in Figure 9.

##  6 Discussion

In $\S$5, we have shown that WISE-selected obscured quasars occupy
significantly more massive halos than unobscured quasars do. In this section,
we will discuss possible systematics in these measurements and speculate on
interpretations of this observed halo occupation difference.

###  6.1 Redshift Systematics

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
$\mu_{\mathrm{obs}}=1.4$, $\sigma_{\mathrm{obs}}=0.5$. We then vary these
parameters, shifting the mean by $-1<\mu_{\mathrm{sim}}-\mu_{\mathrm{obs}}<3$
and varying the dispersion between
$0.1<\sigma_{\mathrm{sim}}/\sigma_{\mathrm{obs}}<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.

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\approx 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 $\sim 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\approx 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.

Figure 10: A corner plot demonstrating a test for systematics in our effective
halo mass results due to photometric redshift distribution uncertainties. We
have approximated the measured photometric redshift distribution of obscured
quasars as a Gaussian, then shifted the mean $\mu$ and dispersion $\sigma$
from their observed values until the implied effective halo mass for obscured
quasars matches that of unobscured quasars. We show the 68$\%$ and 95$\%$
confidence contours for the redshift distribution properties required to
resolve the halo mass discrepancy for the angular clustering and CMB lensing
measurements in teal and purple, respectively. In order to resolve the halo
mass tension while reproducing both the clustering and lensing signals, the
half of our obscured quasar sample in deep fields without spectroscopic
redshifts would need to be distributed around $z\approx 3.5$ instead of the
value derived from photometric redshifts of $z=1.4$, which we regard as
implausible. Our result is thus robust against redshift systematics.

###  6.2 Interpretation of Clustering Difference

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 non-
evolutionary 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 non-
evolutionary 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 nor 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 (SFR). 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
infrared-selected quasars in the COSMOS field, finding no difference in the
stellar mass distribution for obscured and unobscured quasars.

###  6.3 Non-torus-obscured Population

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 exorcise 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_{\mathrm{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:

$b_{\mathrm{NTO}}=\frac{b_{\mathrm{ob}}-b_{\mathrm{unob}}}{f_{\mathrm{NTO}}}+b_{\mathrm{unob}}$
(18)

We convert between this NTO bias and a minimum halo mass required to host a
quasar as a function of the NTO fraction $f_{\mathrm{NTO}}$ using:

$b(M>M_{\mathrm{min}})=\frac{\int_{M_{\mathrm{min}}}^{\infty}dM\frac{dn}{dM}b(M)}{\int_{M_{\mathrm{min}}}^{\infty}dM\frac{dn}{dM}}.$
(19)

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. This appears to be a compelling line of evidence
that a non-negligible $(\gtrsim 15\%)$ fraction of obscured quasars are
obscured by a mechanism aside from orientation, such as galactic or
circumnuclear dust. .

Figure 11:  the implied minimum host halo mass for “non-torus-obscured” (NTO)
and unobscured quasars as a function of the NTO fraction of our sample. If
only a small fraction of obscured quasars belong to the NTO population, the
bias and host halo masses of these systems must be very large in order to
drive the observed clustering of the entire obscured population. Conversely,
this constrains the NTO fraction to be $\gtrsim 15\%$ in order to avoid
unphysically massive host halos, where the quasar occupation fraction exceeds
unity (see $\S$6.3).

###  6.4 Possible Connection to Dust-obscured Galaxies

By matching our quasar catalog to SDWFS 24 $\mu$m photometry, we find that
$\sim 30\%$ of the objects in the the obscured quasar sample would be
classified as dust-obscured galaxies (DOGs, Dey et al., 2008) with
$F_{24}/F_{r}>1000$. DOGs are ultra-luminous infrared galaxies (ULIRGs) 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 post-merger 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 infrared colors satisfy the
Donley et al. (2012) criterion of power-law infrared spectra. We observe that
the DOG fraction increases in our obscured quasar sample towards 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
$\sim 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.

Figure 12: A further exploration of the effective host halo masses of quasars
as a function of obscuration. We have now split the obscured sample into two
subsets using a color cut ($r$-W2 [AB] = 4.5). We show the effective halo
masses of the samples split by optical-infrared color using gradients, where
the vertical span represents the 68$\%$ confidence interval, and the opacity
of the gradient represents the relative number of sources of a given color.
The effective halo mass appears to continue increasing for obscured quasars
towards redder colors, suggesting the enhanced clustering may be dominantly
driven by the reddest and most obscured subset.

Interestingly, the effective halo mass of WISE quasars appears to continue
increasing towards more extreme red optical-infrared 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 $\sim
10^{13.0}h^{-1}M_{\odot}$. 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 $\sim
10^{13.6}h^{-1}M_{\odot}$. Given a majority of the reddest sample would be
selected as DOGs, we speculate that the observed enhanced clustering of
obscured infrared-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.

###  6.5 Comparison to Literature Results

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 infrared
(Hickox et al., 2011; Donoso et al., 2014; DiPompeo et al., 2014, 2015, 2016a,
2017a, this work) or ultra-hard 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 infrared quasar selection may uncover
obscured quasars which 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 wavebands. 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 infrared 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 towards redder optical-infrared 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; Del Moro et al., 2016; Assef 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 infrared-selected obscured quasars is
driven by deeply-buried systems missed by soft X-ray surveys.

##  7 Conclusions

In this paper, we have selected $\sim 1.4$ million quasars over $\sim 15,000\
\mathrm{deg}^{2}$ of sky in the mid-infrared with WISE and split them into
obscured and unobscured samples using optical-infrared colors. We have then
probed their host halo properties by measuring their angular clustering
statistics as well as the typical gravitational deflection of CMB photons
passing through their surrounding large-scale structures. We have interpreted
the clustering signals within an HOD framework, finding that the minimum and
effective halo masses of obscured quasars are significantly higher than their
unobscured counterparts, with obscured quasars occupying effective halos of
$\sim 10^{12.9}h^{-1}M_{\odot}$, while unobscured quasars occupy effective
halos of $\sim 10^{12.6}h^{-1}M_{\odot}$. We however do not detect a
difference in the one-halo term, finding that $\sim 5-20\%$ of both obscured
and unobscured quasars are satellites within their halos. Interpreting the CMB
lensing signals with a linearly biased halo model, we find excellent agreement
with the clustering results. We showed that this agreement confirms the halo
mass difference is robust against uncertainty in the obscured sample’s
redshift distribution. We discuss interpretations of the observed clustering
difference, and favor an evolutionary explanation for the obscuration of at
least some quasars. Finally, we detect a hint that the enhanced clustering is
driven by the systems with the most extremely red optical-infrared colors,
which have a significant overlap (¿$50\%$) with the Dust-obscured Galaxy (DOG)
population.

With this work, we have now estimated the halo properties of infrared quasars
across the majority of the extragalactic sky. We are therefore likely
approaching the limit of the statistical power available to probe the halos of
purely WISE-selected quasars through angular autocorrelations. Consequently,
further investigation of the halo properties of these quasars must utilize
other techniques. While the Planck map of CMB lensing will remain the only
all-sky dataset for the near future, ground-based experiments such as the
Atacama Cosmology Telescope, South Pole Telescope, and soon the Simons
Observatory and CMB-S4 will improve the lensing precision on smaller scales
(e.g., Omori et al., 2017; Darwish et al., 2021), and may allow direct
measurements of the quasar one-halo term. The weak lensing of background
galaxies can also be used to study AGN halo properties (e.g., Leauthaud et
al., 2015; Luo et al., 2022). Upcoming missions including Euclid and the Rubin
Observatory will measure the weak gravitational lensing of galaxies at $z>2$
over a wide area, probing the large-scale structure around the peak of quasar
activity at $z=1-2$. These missions will also detect and measure photometric
redshifts for galaxies at $z=1-2$, enabling cross-clustering measurements to
increase the signal-to-noise ratio for the inherently rare population of
bright quasars. Finally, spectroscopic surveys of quasars would allow spatial
clustering measurements and more precise obscuration criteria, but extending
spectroscopy to the most heavily obscured systems will prove challenging due
to their inherent optical faintness.

In order to interpret the clustering difference between the obscured and
unobscured quasars observed in this work, a comparison of the host galaxy
properties through full panchromatic SED modeling (e.g., Andonie et al., 2022)
is required in future work. Such a characterization, combined with the results
of this work, will provide a probe of the full halo-galaxy-SMBH connection,
testing models of AGN structure and of galaxy-SMBH coevolution.

We thank the anonymous reviewer for comments which significantly improved this
work. GCP acknowledges support from the Dartmouth Fellowship. DMA thanks the
Science Technology Facilities Council (STFC) for support from the Durham
consolidated grant (ST/T000244/1). CA acknowledges support from EU H2020-MSCA-
ITN-2019 Project 860744 “BiD4BESt: Big Data applications for black hole
Evolution STudies. This research has made use of the NASA/IPAC Infrared
Science Archive, which is funded by the National Aeronautics and Space
Administration and operated by the California Institute of Technology. This
research uses services or data provided by the Astro Data Lab at NSF’s
National Optical-Infrared Astronomy Research Laboratory. NOIRLab is operated
by the Association of Universities for Research in Astronomy (AURA), Inc.
under a cooperative agreement with the National Science Foundation. The Legacy
Surveys consist of three individual and complementary projects: the Dark
Energy Camera Legacy Survey (DECaLS; Proposal ID 2014B-0404; PIs: David
Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop. ID
2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey
(MzLS; Prop. ID 2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS together
include data obtained, respectively, at the Blanco telescope, Cerro Tololo
Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward
Observatory, University of Arizona; and the Mayall telescope, Kitt Peak
National Observatory, NOIRLab. The Legacy Surveys project is honored to be
permitted to conduct astronomical research on Iolkam Du’ag (Kitt Peak), a
mountain with particular significance to the Tohono O’odham Nation. This
project used data obtained with the Dark Energy Camera (DECam), which was
constructed by the Dark Energy Survey (DES) collaboration. , MANGLE: (Hamilton
& Tegmark, 2004; Swanson et al., 2008), TOPCAT: (Taylor, 2005).

## References

  * Acker et al. (1992) Acker, A., Marcout, J., Ochsenbein, F., et al. 1992, The Strasbourg-ESO Catalogue of Galactic Planetary Nebulae. Parts I, II.
  * Ahumada et al. (2020) Ahumada, R., Prieto, C. A., Almeida, A., et al. 2020, ApJS, 249, 3, doi: 10.3847/1538-4365/ab929e
  * Aihara et al. (2018) Aihara, H., Armstrong, R., Bickerton, S., et al. 2018, PASJ, 70, S8, doi: 10.1093/pasj/psx081
  * Aird & Coil (2021) Aird, J., & Coil, A. L. 2021, MNRAS, 502, 5962, doi: 10.1093/mnras/stab312
  * Aird et al. (2018) Aird, J., Coil, A. L., & Georgakakis, A. 2018, MNRAS, 474, 1225, doi: 10.1093/mnras/stx2700
  * Aird et al. (2015) Aird, J., Coil, A. L., Georgakakis, A., et al. 2015, MNRAS, 451, 1892, doi: 10.1093/mnras/stv1062
  * Alam et al. (2020) Alam, S., Peacock, J. A., Kraljic, K., Ross, A. J., & Comparat, J. 2020, MNRAS, 497, 581, doi: 10.1093/mnras/staa1956
  * Alexander & Hickox (2012) Alexander, D. M., & Hickox, R. C. 2012, New A Rev., 56, 93, doi: 10.1016/j.newar.2011.11.003
  * Alexander et al. (2008) Alexander, D. M., Brandt, W. N., Smail, I., et al. 2008, AJ, 135, 1968, doi: 10.1088/0004-6256/135/5/1968
  * Alexander et al. (2022) Alexander, D. M., Davis, T. M., Chaussidon, E., et al. 2022, arXiv e-prints, arXiv:2208.08517.  https://arxiv.org/abs/2208.08517
  * Allevato et al. (2014) Allevato, V., Finoguenov, A., Civano, F., et al. 2014, ApJ, 796, 4, doi: 10.1088/0004-637X/796/1/4
  * Alonso et al. (2019) Alonso, D., Sanchez, J., Slosar, A., & LSST Dark Energy Science Collaboration. 2019, MNRAS, 484, 4127, doi: 10.1093/mnras/stz093
  * Ananna et al. (2019) Ananna, T. T., Treister, E., Urry, C. M., et al. 2019, ApJ, 871, 240, doi: 10.3847/1538-4357/aafb77
  * Ananna et al. (2022) Ananna, T. T., Weigel, A. K., Trakhtenbrot, B., et al. 2022, ApJS, 261, 9, doi: 10.3847/1538-4365/ac5b64
  * Anderson et al. (2014) Anderson, L. D., Bania, T. M., Balser, D. S., et al. 2014, ApJS, 212, 1, doi: 10.1088/0067-0049/212/1/1
  * Andonie et al. (2022) Andonie, C., Alexander, D. M., Rosario, D., et al. 2022, Monthly Notices of the Royal Astronomical Society, doi: 10.1093/mnras/stac2800
  * Antonucci (1993) Antonucci, R. 1993, ARA&A, 31, 473, doi: 10.1146/annurev.aa.31.090193.002353
  * Ashby et al. (2009) Ashby, M. L. N., Stern, D., Brodwin, M., et al. 2009, ApJ, 701, 428, doi: 10.1088/0004-637X/701/1/428
  * Assef et al. (2018) Assef, R. J., Stern, D., Noirot, G., et al. 2018, ApJS, 234, 23, doi: 10.3847/1538-4365/aaa00a
  * Assef et al. (2013) Assef, R. J., Stern, D., Kochanek, C. S., et al. 2013, ApJ, 772, 26, doi: 10.1088/0004-637X/772/1/26
  * Assef et al. (2016) Assef, R. J., Walton, D. J., Brightman, M., et al. 2016, ApJ, 819, 111, doi: 10.3847/0004-637X/819/2/111
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Barrows et al. (2021) Barrows, R. S., Comerford, J. M., Stern, D., & Assef, R. J. 2021, ApJ, 922, 179, doi: 10.3847/1538-4357/ac1352
  * Berlind & Weinberg (2002) Berlind, A. A., & Weinberg, D. H. 2002, ApJ, 575, 587, doi: 10.1086/341469
  * Buchner et al. (2017) Buchner, J., Schulze, S., & Bauer, F. E. 2017, MNRAS, 464, 4545, doi: 10.1093/mnras/stw2423
  * Buchner et al. (2015) Buchner, J., Georgakakis, A., Nandra, K., et al. 2015, ApJ, 802, 89, doi: 10.1088/0004-637X/802/2/89
  * Bussmann et al. (2012) Bussmann, R. S., Dey, A., Armus, L., et al. 2012, ApJ, 744, 150, doi: 10.1088/0004-637X/744/2/150
  * Cappelluti et al. (2010) Cappelluti, N., Ajello, M., Burlon, D., et al. 2010, ApJ, 716, L209, doi: 10.1088/2041-8205/716/2/L209
  * Casey et al. (2012) Casey, C. M., Berta, S., Béthermin, M., et al. 2012, ApJ, 761, 140, doi: 10.1088/0004-637X/761/2/140
  * Chatterjee et al. (2012) Chatterjee, S., Degraf, C., Richardson, J., et al. 2012, MNRAS, 419, 2657, doi: 10.1111/j.1365-2966.2011.19917.x
  * Chaussidon et al. (2022) Chaussidon, E., Yèche, C., Palanque-Delabrouille, N., et al. 2022, arXiv e-prints, arXiv:2208.08511.  https://arxiv.org/abs/2208.08511
  * Chen et al. (2015) Chen, C.-T. J., Hickox, R. C., Alberts, S., et al. 2015, ApJ, 802, 50, doi: 10.1088/0004-637X/802/1/50
  * Chung et al. (2014) Chung, S. M., Kochanek, C. S., Assef, R., et al. 2014, ApJ, 790, 54, doi: 10.1088/0004-637X/790/1/54
  * Circosta et al. (2019) Circosta, C., Vignali, C., Gilli, R., et al. 2019, A&A, 623, A172, doi: 10.1051/0004-6361/201834426
  * Coil (2013) Coil, A. L. 2013, in Planets, Stars and Stellar Systems. Volume 6: Extragalactic Astronomy and Cosmology, ed. T. D. Oswalt & W. C. Keel, Vol. 6, 387, doi: 10.1007/978-94-007-5609-0_8
  * Coil et al. (2009) Coil, A. L., Georgakakis, A., Newman, J. A., et al. 2009, ApJ, 701, 1484, doi: 10.1088/0004-637X/701/2/1484
  * Coil et al. (2011) Coil, A. L., Blanton, M. R., Burles, S. M., et al. 2011, ApJ, 741, 8, doi: 10.1088/0004-637X/741/1/8
  * Comparat et al. (2015) Comparat, J., Richard, J., Kneib, J.-P., et al. 2015, A&A, 575, A40, doi: 10.1051/0004-6361/201424767
  * Cool et al. (2013) Cool, R. J., Moustakas, J., Blanton, M. R., et al. 2013, ApJ, 767, 118, doi: 10.1088/0004-637X/767/2/118
  * Cooray & Hu (2000) Cooray, A., & Hu, W. 2000, ApJ, 534, 533, doi: 10.1086/308799
  * Cutri et al. (2021) Cutri, R. M., Wright, E. L., Conrow, T., et al. 2021, VizieR Online Data Catalog, II/328
  * Darwish et al. (2021) Darwish, O., Madhavacheril, M. S., Sherwin, B. D., et al. 2021, MNRAS, 500, 2250, doi: 10.1093/mnras/staa3438
  * Del Moro et al. (2016) Del Moro, A., Alexander, D. M., Bauer, F. E., et al. 2016, MNRAS, 456, 2105, doi: 10.1093/mnras/stv2748
  * Delvecchio et al. (2017) Delvecchio, I., Smolčić, V., Zamorani, G., et al. 2017, A&A, 602, A3, doi: 10.1051/0004-6361/201629367
  * Dey et al. (2008) Dey, A., Soifer, B. T., Desai, V., et al. 2008, ApJ, 677, 943, doi: 10.1086/529516
  * Dey et al. (2019) Dey, A., Schlegel, D. J., Lang, D., et al. 2019, AJ, 157, 168, doi: 10.3847/1538-3881/ab089d
  * Diemer (2018) Diemer, B. 2018, ApJS, 239, 35, doi: 10.3847/1538-4365/aaee8c
  * DiPompeo et al. (2017a) DiPompeo, M. A., Hickox, R. C., Eftekharzadeh, S., & Myers, A. D. 2017a, MNRAS, 469, 4630, doi: 10.1093/mnras/stx1215
  * DiPompeo et al. (2016a) DiPompeo, M. A., Hickox, R. C., & Myers, A. D. 2016a, MNRAS, 456, 924, doi: 10.1093/mnras/stv2681
  * DiPompeo et al. (2017b) DiPompeo, M. A., Hickox, R. C., Myers, A. D., & Geach, J. E. 2017b, MNRAS, 464, 3526, doi: 10.1093/mnras/stw2589
  * DiPompeo et al. (2014) DiPompeo, M. A., Myers, A. D., Hickox, R. C., Geach, J. E., & Hainline, K. N. 2014, MNRAS, 442, 3443, doi: 10.1093/mnras/stu1115
  * DiPompeo et al. (2015) DiPompeo, M. A., Myers, A. D., Hickox, R. C., et al. 2015, MNRAS, 446, 3492, doi: 10.1093/mnras/stu2341
  * DiPompeo et al. (2016b) DiPompeo, M. A., Runnoe, J. C., Hickox, R. C., Myers, A. D., & Geach, J. E. 2016b, MNRAS, 460, 175, doi: 10.1093/mnras/stw986
  * Donley et al. (2012) Donley, J. L., Koekemoer, A. M., Brusa, M., et al. 2012, ApJ, 748, 142, doi: 10.1088/0004-637X/748/2/142
  * Donoso et al. (2014) Donoso, E., Yan, L., Stern, D., & Assef, R. J. 2014, ApJ, 789, 44, doi: 10.1088/0004-637X/789/1/44
  * Duncan et al. (2021) Duncan, K. J., Kondapally, R., Brown, M. J. I., et al. 2021, A&A, 648, A4, doi: 10.1051/0004-6361/202038809
  * Ebrero et al. (2009) Ebrero, J., Mateos, S., Stewart, G. C., Carrera, F. J., & Watson, M. G. 2009, A&A, 500, 749, doi: 10.1051/0004-6361/200911670
  * Efron (1982) Efron, B. 1982, The Jackknife, the Bootstrap and other resampling plans
  * Eisenhardt et al. (2020) Eisenhardt, P. R. M., Marocco, F., Fowler, J. W., et al. 2020, ApJS, 247, 69, doi: 10.3847/1538-4365/ab7f2a
  * Eisenstein & Hu (1998) Eisenstein, D. J., & Hu, W. 1998, ApJ, 496, 605, doi: 10.1086/305424
  * Elvis (2000) Elvis, M. 2000, ApJ, 545, 63, doi: 10.1086/317778
  * Elyiv et al. (2012) Elyiv, A., Clerc, N., Plionis, M., et al. 2012, A&A, 537, A131, doi: 10.1051/0004-6361/201117983
  * Fedeli et al. (2014) Fedeli, C., Semboloni, E., Velliscig, M., et al. 2014, J. Cosmology Astropart. Phys, 2014, 028, doi: 10.1088/1475-7516/2014/08/028
  * Fernique et al. (2014) Fernique, P., Boch, T., Donaldson, T., et al. 2014, MOC - HEALPix Multi-Order Coverage map Version 1.0, IVOA Recommendation 02 June 2014, doi: 10.5479/ADS/bib/2014ivoa.spec.0602F
  * Fiore et al. (2008) Fiore, F., Grazian, A., Santini, P., et al. 2008, ApJ, 672, 94, doi: 10.1086/523348
  * Flesch (2019) Flesch, E. W. 2019, arXiv e-prints, arXiv:1912.05614.  https://arxiv.org/abs/1912.05614
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Geach et al. (2013) Geach, J. E., Hickox, R. C., Bleem, L. E., et al. 2013, ApJ, 776, L41, doi: 10.1088/2041-8205/776/2/L41
  * Georgakakis et al. (2017) Georgakakis, A., Aird, J., Schulze, A., et al. 2017, MNRAS, 471, 1976, doi: 10.1093/mnras/stx1602
  * Georgakakis et al. (2019) Georgakakis, A., Comparat, J., Merloni, A., et al. 2019, MNRAS, 487, 275, doi: 10.1093/mnras/sty3454
  * Gilli et al. (2009) Gilli, R., Zamorani, G., Miyaji, T., et al. 2009, A&A, 494, 33, doi: 10.1051/0004-6361:200810821
  * Gilli et al. (2022) Gilli, R., Norman, C., Calura, F., et al. 2022, arXiv e-prints, arXiv:2206.03508.  https://arxiv.org/abs/2206.03508
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759, doi: 10.1086/427976
  * Haiman & Hui (2001) Haiman, Z., & Hui, L. 2001, ApJ, 547, 27, doi: 10.1086/318330
  * Hamilton (2000) Hamilton, A. J. S. 2000, MNRAS, 312, 257, doi: 10.1046/j.1365-8711.2000.03071.x
  * Hamilton & Tegmark (2004) Hamilton, A. J. S., & Tegmark, M. 2004, MNRAS, 349, 115, doi: 10.1111/j.1365-2966.2004.07490.x
  * Hasinger et al. (2018) Hasinger, G., Capak, P., Salvato, M., et al. 2018, ApJ, 858, 77, doi: 10.3847/1538-4357/aabacf
  * Hernán-Caballero et al. (2016) Hernán-Caballero, A., Spoon, H. W. W., Lebouteiller, V., Rupke, D. S. N., & Barry, D. P. 2016, MNRAS, 455, 1796, doi: 10.1093/mnras/stv2464
  * Hickox & Alexander (2018) Hickox, R. C., & Alexander, D. M. 2018, ARA&A, 56, 625, doi: 10.1146/annurev-astro-081817-051803
  * Hickox et al. (2017) Hickox, R. C., Myers, A. D., Greene, J. E., et al. 2017, ApJ, 849, 53, doi: 10.3847/1538-4357/aa8c77
  * Hickox et al. (2007) Hickox, R. C., Jones, C., Forman, W. R., et al. 2007, ApJ, 671, 1365, doi: 10.1086/523082
  * Hickox et al. (2011) Hickox, R. C., Myers, A. D., Brodwin, M., et al. 2011, ApJ, 731, 117, doi: 10.1088/0004-637X/731/2/117
  * Hivon et al. (2002) Hivon, E., Górski, K. M., Netterfield, C. B., et al. 2002, ApJ, 567, 2, doi: 10.1086/338126
  * Hopkins et al. (2005) Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2005, ApJ, 630, 705, doi: 10.1086/432438
  * Hopkins et al. (2006) —. 2006, ApJS, 163, 1, doi: 10.1086/499298
  * Hopkins et al. (2008) Hopkins, P. F., Hernquist, L., Cox, T. J., & Kereš, D. 2008, ApJS, 175, 356, doi: 10.1086/524362
  * Ilbert et al. (2009) Ilbert, O., Capak, P., Salvato, M., et al. 2009, ApJ, 690, 1236, doi: 10.1088/0004-637X/690/2/1236
  * Jiang et al. (2016) Jiang, N., Wang, H., Mo, H., et al. 2016, ApJ, 832, 111, doi: 10.3847/0004-637X/832/2/111
  * Juneau et al. (2022) Juneau, S., Goulding, A. D., Banfield, J., et al. 2022, ApJ, 925, 203, doi: 10.3847/1538-4357/ac425f
  * Karachentsev et al. (2013) Karachentsev, I. D., Makarov, D. I., & Kaisina, E. I. 2013, AJ, 145, 101, doi: 10.1088/0004-6256/145/4/101
  * Kartaltepe et al. (2015) Kartaltepe, J. S., Sanders, D. B., Silverman, J. D., et al. 2015, ApJ, 806, L35, doi: 10.1088/2041-8205/806/2/L35
  * Kashino et al. (2019) Kashino, D., Silverman, J. D., Sanders, D., et al. 2019, ApJS, 241, 10, doi: 10.3847/1538-4365/ab06c4
  * Kenter et al. (2005) Kenter, A., Murray, S. S., Forman, W. R., et al. 2005, ApJS, 161, 9, doi: 10.1086/444379
  * Kochanek et al. (2012) Kochanek, C. S., Eisenstein, D. J., Cool, R. J., et al. 2012, ApJS, 200, 8, doi: 10.1088/0067-0049/200/1/8
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511, doi: 10.1146/annurev-astro-082708-101811
  * Koutoulidis et al. (2018) Koutoulidis, L., Georgantopoulos, I., Mountrichas, G., et al. 2018, MNRAS, 481, 3063, doi: 10.1093/mnras/sty2429
  * Krishnan et al. (2020) Krishnan, C., Almaini, O., Hatch, N. A., et al. 2020, MNRAS, 494, 1693, doi: 10.1093/mnras/staa815
  * Krolewski et al. (2020) Krolewski, A., Ferraro, S., Schlafly, E. F., & White, M. 2020, J. Cosmology Astropart. Phys, 2020, 047, doi: 10.1088/1475-7516/2020/05/047
  * Krumpe et al. (2018) Krumpe, M., Miyaji, T., Coil, A. L., & Aceves, H. 2018, MNRAS, 474, 1773, doi: 10.1093/mnras/stx2705
  * Lacy et al. (2004) Lacy, M., Storrie-Lombardi, L. J., Sajina, A., et al. 2004, ApJS, 154, 166, doi: 10.1086/422816
  * Lacy et al. (2013) Lacy, M., Ridgway, S. E., Gates, E. L., et al. 2013, ApJS, 208, 24, doi: 10.1088/0067-0049/208/2/24
  * Laigle et al. (2016) Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24, doi: 10.3847/0067-0049/224/2/24
  * Landy & Szalay (1993) Landy, S. D., & Szalay, A. S. 1993, ApJ, 412, 64, doi: 10.1086/172900
  * Lanzuisi et al. (2009) Lanzuisi, G., Piconcelli, E., Fiore, F., et al. 2009, A&A, 498, 67, doi: 10.1051/0004-6361/200811282
  * Laurent et al. (2017) Laurent, P., Eftekharzadeh, S., Le Goff, J.-M., et al. 2017, J. Cosmology Astropart. Phys, 2017, 017, doi: 10.1088/1475-7516/2017/07/017
  * Leauthaud et al. (2015) Leauthaud, A., J. Benson, A., Civano, F., et al. 2015, MNRAS, 446, 1874, doi: 10.1093/mnras/stu2210
  * Lewis et al. (2000) Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473, doi: 10.1086/309179
  * Li et al. (2021) Li, J., Silverman, J. D., Ding, X., et al. 2021, ApJ, 918, 22, doi: 10.3847/1538-4357/ac06a8
  * Lilly et al. (2007) Lilly, S. J., Le Fèvre, O., Renzini, A., et al. 2007, ApJS, 172, 70, doi: 10.1086/516589
  * Limber (1953) Limber, D. N. 1953, ApJ, 117, 134, doi: 10.1086/145672
  * LoVerde & Afshordi (2008) LoVerde, M., & Afshordi, N. 2008, Phys. Rev. D, 78, 123506, doi: 10.1103/PhysRevD.78.123506
  * Luo et al. (2022) Luo, W., Silverman, J. D., More, S., et al. 2022, arXiv e-prints, arXiv:2204.03817.  https://arxiv.org/abs/2204.03817
  * Lyke et al. (2020) Lyke, B. W., Higley, A. N., McLane, J. N., et al. 2020, ApJS, 250, 8, doi: 10.3847/1538-4365/aba623
  * Lynds (1962) Lynds, B. T. 1962, ApJS, 7, 1, doi: 10.1086/190072
  * Lynds (1965) —. 1965, ApJS, 12, 163, doi: 10.1086/190123
  * Marchesi et al. (2016) Marchesi, S., Civano, F., Elvis, M., et al. 2016, ApJ, 817, 34, doi: 10.3847/0004-637X/817/1/34
  * Marocco et al. (2021) Marocco, F., Eisenhardt, P. R. M., Fowler, J. W., et al. 2021, ApJS, 253, 8, doi: 10.3847/1538-4365/abd805
  * Marocco et al. (2020) Marocco et al. 2020, CatWISE2020 Catalog, IPAC, doi: 10.26131/IRSA551
  * Martini & Weinberg (2001) Martini, P., & Weinberg, D. H. 2001, ApJ, 547, 12, doi: 10.1086/318331
  * Mead et al. (2021) Mead, A. J., Brieden, S., Tröster, T., & Heymans, C. 2021, MNRAS, 502, 1401, doi: 10.1093/mnras/stab082
  * Mead et al. (2015) Mead, A. J., Peacock, J. A., Heymans, C., Joudaki, S., & Heavens, A. F. 2015, MNRAS, 454, 1958, doi: 10.1093/mnras/stv2036
  * Meisner et al. (2019) Meisner, A. M., Lang, D., Schlafly, E. F., & Schlegel, D. J. 2019, PASP, 131, 124504, doi: 10.1088/1538-3873/ab3df4
  * Mendez et al. (2016) Mendez, A. J., Coil, A. L., Aird, J., et al. 2016, ApJ, 821, 55, doi: 10.3847/0004-637X/821/1/55
  * Messias et al. (2012) Messias, H., Afonso, J., Salvato, M., Mobasher, B., & Hopkins, A. M. 2012, ApJ, 754, 120, doi: 10.1088/0004-637X/754/2/120
  * Mitra et al. (2018) Mitra, K., Chatterjee, S., DiPompeo, M. A., Myers, A. D., & Zheng, Z. 2018, MNRAS, 477, 45, doi: 10.1093/mnras/sty556
  * Mountrichas & Georgakakis (2012) Mountrichas, G., & Georgakakis, A. 2012, MNRAS, 420, 514, doi: 10.1111/j.1365-2966.2011.20059.x
  * Murray et al. (2021) Murray, S. G., Diemer, B., Chen, Z., et al. 2021, Astronomy and Computing, 36, 100487, doi: 10.1016/j.ascom.2021.100487
  * Myers et al. (2006) Myers, A. D., Brunner, R. J., Richards, G. T., et al. 2006, ApJ, 638, 622, doi: 10.1086/499093
  * Myers et al. (2022) Myers, A. D., Moustakas, J., Bailey, S., et al. 2022, arXiv e-prints, arXiv:2208.08518.  https://arxiv.org/abs/2208.08518
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493, doi: 10.1086/304888
  * Netzer (2015) Netzer, H. 2015, ARA&A, 53, 365, doi: 10.1146/annurev-astro-082214-122302
  * Norberg et al. (2009) Norberg, P., Baugh, C. M., Gaztañaga, E., & Croton, D. J. 2009, MNRAS, 396, 19, doi: 10.1111/j.1365-2966.2009.14389.x
  * Omori et al. (2017) Omori, Y., Chown, R., Simard, G., et al. 2017, ApJ, 849, 124, doi: 10.3847/1538-4357/aa8d1d
  * Onodera et al. (2015) Onodera, M., Carollo, C. M., Renzini, A., et al. 2015, ApJ, 808, 161, doi: 10.1088/0004-637X/808/2/161
  * Padovani et al. (2017) Padovani, P., Alexander, D. M., Assef, R. J., et al. 2017, A&A Rev., 25, 2, doi: 10.1007/s00159-017-0102-9
  * Pannella et al. (2009) Pannella, M., Carilli, C. L., Daddi, E., et al. 2009, ApJ, 698, L116, doi: 10.1088/0004-637X/698/2/L116
  * Peacock (1991) Peacock, J. A. 1991, MNRAS, 253, 1P, doi: 10.1093/mnras/253.1.1P
  * Pedregosa et al. (2012) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2012, arXiv e-prints, arXiv:1201.0490.  https://arxiv.org/abs/1201.0490
  * Peebles (1973) Peebles, P. J. E. 1973, ApJ, 185, 413, doi: 10.1086/152431
  * Peebles (1980) —. 1980, The large-scale structure of the universe
  * Petter et al. (2022) Petter, G. C., Hickox, R. C., Alexander, D. M., et al. 2022, ApJ, 927, 16, doi: 10.3847/1538-4357/ac4d31
  * Piconcelli et al. (2015) Piconcelli, E., Vignali, C., Bianchi, S., et al. 2015, A&A, 574, L9, doi: 10.1051/0004-6361/201425324
  * Planck Collaboration et al. (2020a) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020a, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Planck Collaboration et al. (2020b) —. 2020b, A&A, 641, A8, doi: 10.1051/0004-6361/201833886
  * Polletta et al. (2007) Polletta, M., Tajer, M., Maraschi, L., et al. 2007, ApJ, 663, 81, doi: 10.1086/518113
  * Powell et al. (2020) Powell, M. C., Urry, C. M., Cappelluti, N., et al. 2020, ApJ, 891, 41, doi: 10.3847/1538-4357/ab6e65
  * Powell et al. (2018) Powell, M. C., Cappelluti, N., Urry, C. M., et al. 2018, ApJ, 858, 110, doi: 10.3847/1538-4357/aabd7f
  * Prescott et al. (2006) Prescott, M. K. M., Impey, C. D., Cool, R. J., & Scoville, N. Z. 2006, ApJ, 644, 100, doi: 10.1086/503325
  * Ramos Almeida & Ricci (2017) Ramos Almeida, C., & Ricci, C. 2017, Nature Astronomy, 1, 679, doi: 10.1038/s41550-017-0232-z
  * Ricci et al. (2017a) Ricci, C., Bauer, F. E., Treister, E., et al. 2017a, MNRAS, 468, 1273, doi: 10.1093/mnras/stx173
  * Ricci et al. (2017b) Ricci, C., Trakhtenbrot, B., Koss, M. J., et al. 2017b, Nature, 549, 488, doi: 10.1038/nature23906
  * Richardson et al. (2012) Richardson, J., Zheng, Z., Chatterjee, S., Nagai, D., & Shen, Y. 2012, ApJ, 755, 30, doi: 10.1088/0004-637X/755/1/30
  * Salvato et al. (2009) Salvato, M., Hasinger, G., Ilbert, O., et al. 2009, ApJ, 690, 1250, doi: 10.1088/0004-637X/690/2/1250
  * Salvato et al. (2011) Salvato, M., Ilbert, O., Hasinger, G., et al. 2011, ApJ, 742, 61, doi: 10.1088/0004-637X/742/2/61
  * Sanders et al. (1988) Sanders, D. B., Soifer, B. T., Elias, J. H., et al. 1988, ApJ, 325, 74, doi: 10.1086/165983
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525, doi: 10.1086/305772
  * Schmidt (1963) Schmidt, M. 1963, Nature, 197, 1040, doi: 10.1038/1971040a0
  * Schulze et al. (2018) Schulze, A., Silverman, J. D., Kashino, D., et al. 2018, ApJS, 239, 22, doi: 10.3847/1538-4365/aae82f
  * Scoccimarro et al. (2001) Scoccimarro, R., Sheth, R. K., Hui, L., & Jain, B. 2001, ApJ, 546, 20, doi: 10.1086/318261
  * Scoville et al. (2007) Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1, doi: 10.1086/516585
  * Shen et al. (2020) Shen, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2020, MNRAS, 495, 3252, doi: 10.1093/mnras/staa1381
  * Shirley et al. (2019) Shirley, R., Roehlly, Y., Hurley, P. D., et al. 2019, MNRAS, 490, 634, doi: 10.1093/mnras/stz2509
  * Shirley et al. (2021) Shirley, R., Duncan, K., Campos Varillas, M. C., et al. 2021, MNRAS, 507, 129, doi: 10.1093/mnras/stab1526
  * Simon (2007) Simon, P. 2007, A&A, 473, 711, doi: 10.1051/0004-6361:20066352
  * Sinha & Garrison (2020) Sinha, M., & Garrison, L. H. 2020, MNRAS, 491, 3022, doi: 10.1093/mnras/stz3157
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163, doi: 10.1086/498708
  * Song et al. (2003) Song, Y.-S., Cooray, A., Knox, L., & Zaldarriaga, M. 2003, ApJ, 590, 664, doi: 10.1086/375188
  * Stern et al. (2005) Stern, D., Eisenhardt, P., Gorjian, V., et al. 2005, ApJ, 631, 163, doi: 10.1086/432523
  * Stern et al. (2012) Stern, D., Assef, R. J., Benford, D. J., et al. 2012, ApJ, 753, 30, doi: 10.1088/0004-637X/753/1/30
  * Stern et al. (2014) Stern, D., Lansbury, G. B., Assef, R. J., et al. 2014, ApJ, 794, 102, doi: 10.1088/0004-637X/794/2/102
  * Swanson et al. (2008) Swanson, M. E. C., Tegmark, M., Hamilton, A. J. S., & Hill, J. C. 2008, MNRAS, 387, 1391, doi: 10.1111/j.1365-2966.2008.13296.x
  * Talman (1978) Talman, J. D. 1978, Journal of Computational Physics, 29, 35, doi: 10.1016/0021-9991(78)90107-9
  * Taylor (2005) Taylor, M. B. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton, & R. Ebert, 29
  * Tinker et al. (2008) Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709, doi: 10.1086/591439
  * Tinker et al. (2010) Tinker, J. L., Robertson, B. E., Kravtsov, A. V., et al. 2010, ApJ, 724, 878, doi: 10.1088/0004-637X/724/2/878
  * Toba et al. (2017) Toba, Y., Nagao, T., Kajisawa, M., et al. 2017, ApJ, 835, 36, doi: 10.3847/1538-4357/835/1/36
  * Trump et al. (2009) Trump, J. R., Impey, C. D., Elvis, M., et al. 2009, ApJ, 696, 1195, doi: 10.1088/0004-637X/696/2/1195
  * Ueda et al. (2014) Ueda, Y., Akiyama, M., Hasinger, G., Miyaji, T., & Watson, M. G. 2014, ApJ, 786, 104, doi: 10.1088/0004-637X/786/2/104
  * Urry & Padovani (1995) Urry, C. M., & Padovani, P. 1995, PASP, 107, 803, doi: 10.1086/133630
  * Vito et al. (2018) Vito, F., Brandt, W. N., Stern, D., et al. 2018, MNRAS, 474, 4528, doi: 10.1093/mnras/stx3120
  * Whalen et al. (2020) Whalen, K. E., Hickox, R. C., DiPompeo, M. A., Richards, G. T., & Myers, A. D. 2020, ApJ, 888, 71, doi: 10.3847/1538-4357/ab5d31
  * Whitaker et al. (2017) Whitaker, K. E., Pope, A., Cybulski, R., et al. 2017, ApJ, 850, 208, doi: 10.3847/1538-4357/aa94ce
  * Woo et al. (2008) Woo, J.-H., Treu, T., Malkan, M. A., & Blandford, R. D. 2008, ApJ, 681, 925, doi: 10.1086/588804
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868, doi: 10.1088/0004-6256/140/6/1868
  * Yan et al. (2021) Yan, W., Hickox, R. C., Chen, C.-T. J., et al. 2021, ApJ, 914, 83, doi: 10.3847/1538-4357/abfaa0
  * Yang et al. (2018) Yang, G., Brandt, W. N., Darvish, B., et al. 2018, MNRAS, 480, 1022, doi: 10.1093/mnras/sty1910
  * Zehavi et al. (2011) Zehavi, I., Zheng, Z., Weinberg, D. H., et al. 2011, ApJ, 736, 59, doi: 10.1088/0004-637X/736/1/59
  * Zheng et al. (2007) Zheng, Z., Coil, A. L., & Zehavi, I. 2007, ApJ, 667, 760, doi: 10.1086/521074
  * Zonca et al. (2019) Zonca, A., Singer, L., Lenz, D., et al. 2019, Journal of Open Source Software, 4, 1298, doi: 10.21105/joss.01298