# The WISE-2MASS Survey: Red Quasars Into the Radio Quiet Regime

E. Glikman Department of Physics, Middlebury College, Middlebury, VT 05753,
USA  M. Lacy National Radio Astronomy Observatory, Charlottesville, VA, USA
S. LaMassa Space Telescope Science Institute, 3700 San Martin Drive, Baltimore
MD, 21218, USA  C. Bradley  Department of Physics, Middlebury College,
Middlebury, VT 05753, USA  S. G. Djorgovski California Institute of
Technology, Pasadena, CA 91125, USA  T. Urrutia Leibniz Institut für
Astrophysik, An der Sternwarte 16, D-14482 Potsdam, Germany  E. L. Gates
UCO/Lick Observatory, P.O. Box 85, Mount Hamilton, CA 95140, USA  M. J. Graham
California Institute of Technology, Pasadena, CA 91125, USA  M. Urry Yale
Center for Astronomy & Astrophysics, Physics Department, P.O. Box 208120, New
Haven, CT 06520, USA  Department of Physics, Yale University, P.O. Box 208121,
New Haven, CT 06520, USA  I. Yoon National Radio Astronomy Observatory,
Charlottesville, VA, USA

###### Abstract

We present a highly complete sample of broad-line (Type 1) QSOs out to $z\sim
3$ selected by their mid-infrared colors, a method that is minimally affected
by dust reddening. We remove host galaxy emission from the spectra and fit for
excess reddening in the residual QSOs, resulting in a Gaussian distribution of
colors for unreddened (blue) QSOs, with a tail extending toward heavily
reddened (red) QSOs, defined as having $E(B-V)>0.25$. This radio-independent
selection method enables us to compare red and blue QSO radio properties in
both the FIRST (1.4 GHz) and VLASS ($2-4$ GHz) surveys. Consistent with recent
results from optically-selected QSOs from SDSS, we find that red QSOs have a
significantly higher detection fraction and a higher fraction of compact radio
morphologies at both frequencies. We employ radio stacking to investigate the
median radio properties of the QSOs including those that are undetected in
FIRST and VLASS, finding that red QSOs have significantly brighter radio
emission and steeper radio spectral slopes compared with blue QSOs. Finally,
we find that the incidence of red QSOs is strongly luminosity dependent, where
red QSOs make up $>40\%$ of all QSOs at the highest luminosities. Overall, red
QSOs comprise $\sim 40\%$ of higher luminosity QSOs, dropping to only a few
percent at lower luminosities. Furthermore, red QSOs make up a larger
percentage of the radio-detected QSO population. We argue that dusty AGN-
driven winds are responsible for both the obscuration as well as excess radio
emission seen in red QSOs.

quasars: general, galaxies: active, infrared: galaxies, surveys

††software: Astropy (Astropy Collaboration et al., 2013, 2018), IRAF (Tody,
1986, 1993), Spextool (Cushing et al., 2004), TOPCAT (Taylor, 2005)

##  1 Introduction

Our incomplete understanding of the relationship between supermassive black
hole (SMBH) growth and the growth of galaxies in the Universe remains an
outstanding problem in astrophysics. There is ample evidence that SMBHs are
linked to their host galaxies through scaling relations, such as the
$M_{BH}-\sigma$ relationship (Gebhardt et al., 2000), which tell us that
galaxies and their nuclear black holes likely grew in tandem. In order for
galaxies to build up stars while growing a nuclear SMBH, a galaxy-scale energy
exchange, or “feedback”, is required to regulate this process and tie the two
systems together. This feedback is still very poorly understood, and may come
in the form of radiation, winds, outflows, and/or jets (c.f., Fabian, 2012).

A population of sources that may help elucidate the nature of feedback are
dust-reddened quasars, which appear to represent an important evolutionary
phase linking galaxy mergers to black hole growth. Luminous quasars are
thought to be triggered by major galaxy mergers (Sanders et al., 1988;
Treister et al., 2012) and simulations of major, gas-rich mergers are able to
reproduce many of the aforementioned correlations and galaxy properties (Di
Matteo et al., 2005; Hopkins et al., 2005). During the merger, some gas loses
angular momentum and feeds the black hole while shocks trigger a starburst. As
the SMBH grows, it starts out in a heavily-obscured state followed by a
relatively brief transitional phase during which the dust is cleared via
feedback mechanisms. Subsequently, an unobscured, blue quasar emerges and
dominates the radiation output for the system. The objects in the brief
transitional phase are moderately reddened (or, red) quasars and can serve as
laboratories for studying how quasar feedback impacts their host galaxies.

Red quasars can be elusive, because their optical and near-infrared colors
resemble those of low-mass stars, which are far more abundant at these
wavelengths. Early work used radio selection to find red quasars and avoid
contamination from red stars, which are weak radio sources (Webster et al.,
1995; White et al., 2003). Results from these studies suggested that red
quasars make up a large fraction (up to $\sim 80\%$) of the overall quasar
population but had been missed by optical selection methods. Glikman et al.
(2004) combined the Faint Images of the Radio Sky at Twenty-centimeters
(FIRST; Becker et al., 1995) radio survey and 2 Micron All Sky Survey (2MASS;
Skrutskie et al., 2006) to develop an efficient selection method for finding
these missed red quasars. Subsequent work identified $\gtrsim 130$ dust-
reddened quasars via the same method (hereafter referred to as F2M quasars
Glikman et al., 2007, 2012; Urrutia et al., 2009) that have broad emission
lines and are moderately obscured by $A_{V}\sim 1-4~{}(0.1<E(B-V)\lesssim
1.5$) across a broad range of redshifts ($0.1<z\lesssim 3$). Follow up studies
of F2M quasars showed that they are accreting with very high Eddington rates
(Urrutia et al., 2012; Kim et al., 2015), are overwhelmingly in merger-
dominated systems (Urrutia et al., 2008; Glikman et al., 2015), and often have
broad absorption lines that are typically associated with outflows and
feedback (LoBALs and FeLoBALs; Urrutia et al., 2009; Farrah et al., 2012;
Glikman et al., 2012). This body of evidence suggests that red quasars are
merger-induced systems, in a transitional phase, emerging from their shrouded
environments, as predicted by the galaxy merger simulations.

Careful comparison with blue quasar samples found that red quasars comprise
$\sim 20-30\%$ of the overall quasar population (Glikman et al., 2012, 2018a),
at least at the highest luminosities. When interpreted as an evolutionary
phase, this fraction implies that the duration of this transitional phase is
$\sim 20-30\%$ as long as the unobscured phase, consistent with theoretical
models of quasar ignition and evolution triggered by a major galaxy merger
(e.g., Hopkins et al., 2005).

However, because the F2M survey used radio selection, those quasars belong to
the rarer radio-loud and radio-intermediate populations that make up $\sim
10\%$ of the overall quasar population. Assuming that radio emission from
these quasars is unrelated to their surrounding dust, we could extend the F2M
results to the entire quasar population. However, if the radio emission and
reddening are not independent, then any conclusions about the red quasar
population derived from the F2M sample, such as the duration of the
transitional phase, could be biased and does not apply to the full quasar
population.

Recent results indeed suggest a correlation between reddening and radio
emission. White et al. (2007) used stacking of FIRST images of known (mostly
radio-quiet) quasars and found that redder quasars have higher median radio
fluxes: objects that are 0.8 mag redder than average have radio fluxes that
are $\sim 3$ times higher than average. Another study of the brightest red
quasars ($K<14.5$) by Georgakakis et al. (2009), using only a $J-K>1.5$ color
selection and no radio constraint, found that 6 out of their 10 objects were
detected in the radio. In a sample of extremely red quasars (ERQs) found in
the Sloan Digital Sky Survey (SDSS) without a radio criterion, all of the mid-
infrared-brightest and reddest sources are detected in FIRST (Ross et al.,
2015; Hamann et al., 2017).

More recently, Klindt et al. (2019) and Fawcett et al. (2020) found distinct
differences between the radio properties of blue and red SDSS quasars. They
find that the redder quasars have a significantly higher detection fraction in
FIRST. When stacked, radio quiet red quasars have higher median radio fluxes
than an unreddened sample. In addition, red quasars’ radio morphologies are
more compact compared with blue quasars. However, the optical selection of
SDSS quasars misses the more heavily reddened sources like those found in the
F2M survey because their optical colors place the sources atop the stellar
locus (Urrutia et al., 2009); most of the red QSOs in the SDSS sample have
$E(B-V)\lesssim 0.2$. To avoid any biases in the SDSS QSO selection algorithm
that misses heavily reddened quasars, a selection method is needed at
wavelengths that are minimally impacted by dust extinction and is also radio
independent. Such a method should have sufficient depth and coverage area to
enable a robust comparison between the red and blue populations.

In this paper, we invoke mid-infrared selection, as it has been shown to
successfully identify broader populations of QSOs111Because the focus of this
paper is on the differences between radio detected and undetected systems, we
adopt the canonical nomenclature that distinguishes quasars – radio-detected
luminous AGN whose radio emission is essential to their selection – from QSOs
– the overall class of luminous AGN.  that were less affected by dust
extinction (e.g., Lacy et al., 2004; Stern et al., 2005; Donley et al., 2012;
Jarrett et al., 2011). The Wide-Field Infrared Space Explorer (WISE; Wright et
al., 2010), scanned the sky at 3.4, 4.6, 12, and 22 $\mu$m down to flux
densities of 0.08, 0.11, 1, and 6 mJy, respectively222The sensitivity limit
varies with latitude due to the scanning cadence of the satellite, with the
deepest fields near the poles., providing the wide-area coverage needed for
identifying large numbers of luminous QSOs. In WISE color-color space QSOs can
be isolated from stars and other extragalactic sources, making the mid-IR an
excellent wavelength region for our purposes (e.g., Mateos et al., 2012;
Jarrett et al., 2011; Stern et al., 2012; Assef et al., 2013, 2018).

In a pilot study, Glikman et al. (2018a) identified a complete sample of QSOs
using near-to-mid-infrared color selection over $\sim 260$ deg2 that overlap
the SDSS Stripe 82 legacy field (Frieman et al., 2008). Here we expand upon
that work over a $7.5\times$ larger area reaching fainter mid-infrared fluxes.
We identify a sample of QSO candidates according to their mid-infrared colors
and obtain spectroscopy of sources missed by SDSS and other optical QSO
surveys. Because mid-infrared selection identifies both blue and red luminous
Type 1 QSOs333However, see LaMassa et al. (2019) and Anicetti et al. (in prep)
for examples of luminous X-ray selected QSOs that are undetected in WISE., we
can compare sources that are drawn from the same mid-infrared selection
criteria.

In this paper, we present a sample of red QSOs (defined as having
$E(B-V)>0.25$; Lacy et al., 2007; Glikman et al., 2018a) without relying on
radio selection and aim to determine whether the fraction of red quasars found
in the F2M survey ($\sim$20%) holds for the full red QSO population, including
radio quiet sources. In addition, we explore the differences between radio-
detected and radio-undetected QSOs in the red and blue populations, as well as
their average radio properties through stacking, to explore possible
differences in the mechanisms giving rise to their radio emission, informing
notions of jet formation, dusty winds, or other physical processes.

The surveys employed in this paper use a mix of AB and Vega for their
photometric zeropoints. Rather than transform to a common system, we adhere to
the native systems presented in each survey. WISE and 2MASS photometry are on
the Vega system, while SDSS uses AB magnitudes. When colors are derived from
catalogs on mixed systems, we provide specificity as to which system we are
using. Throughout this work, we adopt the concordance $\Lambda$CDM cosmology
with $H_{0}=70$ km s-1 Mpc-1, $\Omega_{M}=0.3$, and $\Omega_{\Lambda}=0.7$
when computing cosmology-dependent values (Bennett et al., 2013).

##  2 Sample Selection

Our survey area covers an equatorial region overlapping the SDSS over two
fields: a region spanning a range in right ascension of $\alpha=8^{h}-16^{h}$
and in declination of $\delta=0.5^{\circ}-17^{\circ}$, available for follow-up
in the spring months, and the region over Stripe 82 identical to that in
Glikman et al. (2018a), $\alpha=20^{h}40^{m}-3^{h}56^{m}$ and
$\delta=-1.25^{\circ}$ to $+1.25^{\circ}$ (excluding the region
$00^{h}<\alpha<00^{h}15^{m}$). This covers a total area of 2213 deg2 (1950
deg2 and 2632 in the two regions, respectively).

Our aim is to define a sample of luminous Type 1 QSOs in order to compare the
red and blue sub-populations with minimal reddening and radio biases. We begin
by selecting sources with WISE colors consistent with QSO emission, focusing
on sources brighter than $K=14.7$ mag to enable near-infrared spectroscopy
with 3-m class telescopes (§2.5.2) over the chosen survey area. We include all
spectroscopically confirmed broad-line QSOs from SDSS and identify red QSO
candidates among the sources lacking a spectrum in SDSS, selected by their
optical through near-infrared colors. We perform follow-up spectroscopy of all
such candidates and keep all broad-line QSOs to construct a complete sample of
QSOs that obey uniform mid-to-near infrared selection criteria. To construct
the blue and red QSO subsamples, we fit a reddened QSO template to each
spectrum, subtracting host galaxy emission when necessary, and define red QSOs
as having $E(B-V)>0.25$. Finally, we identify a luminosity-restricted
subsample that excludes the blue quasars that would not have been detected if
they were reddened by $E(B-V)\geq 0.25$ to enable a valid comparison between
the red and blue QSOs.

We note that due to the $K<14.7$ mag limit of our survey, the red QSO
population is incomplete, as there are likely to be red QSOs with
$E(B-V)>0.25$ with lower intrinsic luminosities that when reddened, fall below
the flux limit. On the other hand, our selection does not miss significant
numbers of blue QSOs. Therefore, we can compare the blue and red populations
to arrive at a red quasar fraction, understanding that it is a lower limit.

Figure 1 shows a flowchart following the selection process as described in the
steps below and Table 1 presents an overview of the selection with references
to the sections in the text that elaborate on each step in the process.

Figure 1: Flowchart showing the process for selecting the full QSO sample in this work. Each box reports a major selection stage, with the number of sources that passed each stage shown in boldface. Table 1 complements the flowchart with references to the relevant sections in the text that elaborate on each step in the process. Table 1: Overview of sample selection steps  Step | Definition | Relevant Section Reference | Number of sources
---|---|---|---
I | Mid-Infrared color selection | §2.1 \+ Eqn. 1 | 3808
IIa | SDSS Type 1 QSOs | §2.3 | 1072
IIb | KX selected Type 1 QSOs | §2.4 \+ Eqn. 2 \+ Eqn. 3 \+ §2.5 | 40
III | all W2M QSOs | … | 1112
… | red QSOs | §3 \+ $E(B-V)>0.25$ | 63
… | blue QSOs | §3 \+ $E(B-V)<0.25$ | 1049
IV | luminosity restricted blue QSOs | §4 | 798
… | $z>0.1$ | … | 645

###  2.1 Mid-infrared selection

Since the effects of dust reddening diminish with longer wavelengths, we
expect red QSOs to have nearly the same mid-infrared spectral shape as
unreddened QSOs444At far infrared wavelengths, we expect the spectra to
deviate once again, as the scattered and absorbed UV photons are reprocessed
into thermal radiation from heated dust. However, the QSO continuum dominates
at the WISE wavelengths.. Glikman et al. (2018a, Figure 1) showed that the F2M
red quasars lie in the same region of WISE color space as blue QSOs. We
therefore inform our selection of red QSOs by studying the regions that
normal, blue QSOs occupy in WISE color-color space.

Figure 2: We reproduce Figure 12 from Wright et al. (2010) showing the
location of various astrophysical objects in WISE color-color space. We
overplot a box that outlines our initial selection region (Eqn 1) with thick
purple solid lines and plot within the box sources with SDSS spectra
identified as QSO (blue circles) and GALAXY (orange circles). We also plot our
selected candidates with large outlined circles. Newly confirmed QSOs are
filled red and sources that are not obviously QSOs are filled green. The black
solid line shows our refined selection criterion (Eqn 2) that avoids
significant contamination from non-QSOs (§2.4).

We began by selecting sources from the WISE catalog, which presents its
photometry in the Vega system. We select sources with $0.5<W1-W2<2$ and
$2<W2-W3<4.5$, which is a liberal cut around the location of QSO colors in
these bands, as shown in Figure 12 of Wright et al. (2010), restricting our
sample to sources with $K<14.7$ mag. To generate our sample, we use the IRSA
GATOR catalog selection tool to identify all sources in the AllWISE Source
Catalog obeying the following search criteria via an SQL query:

$\begin{array}[]{l}\tt WHERE~{}(~{}(w1mpro-w2mpro)>0.5~{}and\\\
\tt(w1mpro-w2mpro)<2~{}and\\\ \tt(w2mpro-w3mpro)>2~{}and\\\
\tt(w2mpro-w3mpro)<4.5~{})~{}and\\\ \tt k\\_m\\_2mass<14.7,\\\ \end{array}$
(1)

applied to the two sky regions. The GATOR query returned 3,808 sources (3,330
in the northern region and 478 over Stripe 82).

We match the WISE sources to SDSS DR9 photometric catalog (Ahn et al., 2012)
within 2.0″ using the Centre de Données astronomiques de Strasbourg (CDS)
Upload X-Match service through the TOPCAT tool (Taylor, 2005) and excluded the
handful of sources with $K<10$ mag. This resulted in 3,741 matches. Of these,
2,779 have spectroscopic identifications in SDSS (1398 QSOs, 1372 GALAXYs; 9
STARs). Another seventeen sources were observed with the BOSS spectrograph on
SDSS as part of the DR 14 campaign (Abolfathi et al., 2018). Figure 2 plots
these sources in the WISE $[3.4]-[4.6]$ vs. $[4.6]-[12]$ color-color space
defined by Wright et al. (2010), which shows the location of different
astrophysical populations in this space. The blue circles are SDSS-identified
QSOs and they overlap their expected location (cyan oval) on this diagram. The
orange circles are SDSS-identified GALAXYs and they are concentrated toward
lower $[3.4]-[4.6]$ values (i.e, GALAXYs are bluer at these wavelengths).

###  2.2 Type 1 QSOs

The unification model for AGN (Urry & Padovani, 1995) says that the difference
in viewing angle to the central engine of the accreting-black-hole system
determines the observed spectral shape, including emission line widths. Type 1
(broad-line) sources are understood through this model to be seen at
orientation angles nearer to the pole. Beyond a certain range of viewing
angles, the line-of-sight to the broad line region is blocked by close-in
high-column-density gas and dust (i.e., the so-called ‘torus’). In this study,
we wish to focus on Type 1 sources in both the blue and red samples so as to
compare sources with the same approximate distribution of viewing angles,
knowing that our line-of-sight is not intersecting the dense clouds of the
torus (see Figure 3 for an illustration of this argument). If we include Type
2 (narrow-line) sources in the sample, it becomes more challenging to
determine the location and nature of the obscuring material. Therefore,
focusing exclusively on broad-line sources allows us to directly compare the
blue and red populations. As the canonical definition of a broad-line quasar
requires having line widths $>1000$ km s-1 (Glikman et al., 2004, 2012; Gregg
et al., 1996; Schneider et al., 2003), our sample will focus on all objects
with line widths broader than 1000 km s-1.

Figure 3:  This Figure shows the purported obscuration geometry of high-
Eddington-ratio AGN viewed along an unobscured line-of-sight. In both cases,
broad lines are seen from this viewing angle. Therefore, by focusing on only
broad-line (Type 1) QSOs, under the assumption that red QSOs are not reddened
by the nuclear obscuration (i.e., by the torus, gray regions) but are rather
embedded in a dusty environment (bottom figure) where the dust may arise as
the the consequence of a merger, we can translate the fraction of red QSOs to
the duration of the red QSO phase.

###  2.3 QSOs in SDSS

The bulk of the QSOs in our sample comes from the SDSS, seen among the blue
dots in Figure 2. We selected from the SDSS DR9 spectroscopic catalog all
sources with a class of QSO, including all QSOs in the WISE color selection
box (Eqn. 1; Figure 2; 1398 QSOs).

As noted in Glikman et al. (2018a), among the spectra that SDSS classifies as
QSO there are sources that only show narrow lines. To eliminate the Type 2
QSOs, we utilize the ALPAKA catalog (Mullaney et al., 2013) which provides
detailed line analysis for 25,670 AGN with spectra in SDSS DR7. The spectra of
these sources were fitted with multi-component Gaussians to study their
kinematics and Eddington ratios. Line fitting was performed on [O III] 4959Å,
[O III] 5007Å, [N II] 6548Å, H$\alpha$, and [N II] 6584Å, including a broad
component for forbidden species and an additional broad component for
permitted species, if warranted. The ALPAKA sample is limited to $z<0.4$ in
order not to lose H$\alpha$ beyond the SDSS spectroscopic wavelength limit of
9000Å. We matched the QSO sample to the ALPAKA catalog and found 1016 matches
which we use to examine line-widths and select broad-line sources; 382 sources
are left to be dealt with separately. The ALPAKA catalog classifies sources as
Type 1 when the broad component of their H$\alpha$ line constitutes $\geq
50\%$ of the total line flux and a 600 km s-1 threshold (Mullaney et al.,
2013). However, we impose the additional requirement that the broad component
have a Full Width at Half Maximum (FWHM) velocity, $v_{FWHM}>1000$ km s-1. We
also exclude any source not classified as Type 1 in the ALPAKA catalog; 733
spectra obey these criteria.

For the 382 remaining AGN without line analysis in ALPAKA, we batch-downloaded
the spectra via the Science Archive Server (SAS) web
interface555https://dr15.sdss.org/optical/spectrum/search and used the value-
added measurements provided in their multi-extension FITS headers to further
reduce the size of the SDSS QSO sample by examining the distribution of their
maximum emission-line-width. We examined the emission-line fits performed on
the SDSS spectra through DR9 (Bolton et al., 2012). Single-component Gaussians
were fitted to common UV and optical emission lines for all the spectra, tying
the Balmer line widths to each other. Forbidden line widths are also tied to
each other and fitted with a separate Gaussian. These are provided in the
third extension of the FITS tables (see Bolton et al., 2012, for
details666This link provides details on how to extract the line information
from the FITS tables https://data.sdss.org/datamodel/files/BOSS_SPECTRO_REDUX/
RUN2D/spectra/PLATE4/spec.html). The line widths are reported in terms of the
Gaussian $\sigma$ parameter, which we convert to a FWHM velocity, $v_{\rm
FWHM}=2.355\sigma$. Figure 4 shows the FWHM distribution of the broadest line
component in each spectrum (blue line) where $v_{\rm FWHM}=1000$ km s-1 is
shown with a vertical red line. Thirty-two sources fail this criterion.

We visually examined all the spectra that obeyed the broad-line criterion and
identified several sources with erroneous redshifts. We corrected these using
the catalog of Hewett & Wild (2010) which provides improved redshifts for SDSS
QSOs. We also identified and removed another 11 sources whose spectra were
featureless, with no discernible emission lines, suggesting an error in the
automated line analysis for these spectra. This leaves 339 (382 - 32 -11)
additional broad-line QSOs that we add to the QSO sample.

Therefore, together with the ALPAKA-line-width-selected sources, our QSO
sample contains 1,072 ($733+339$) broad-line QSOs with SDSS spectra. Figure 5
shows a flowchart of the process. As we we show in Section 3, some of these
SDSS-identified QSOs show significant amounts of reddening in their spectra
and will be part of the red QSO subsample.

Figure 4: We show the distribution of the velocity width, $v_{FWHM}$, for the
QSOs derived from both the ALPAKA and Bolton et al. (2012) analyses. The solid
line represents all 1398 SDSS QSOs found within the selection box defined by
Eqn 1. The filled blue histogram shows the line width distribution for QSOs
not analyzed by ALPAKA, whose line widths were derived in Bolton et al.
(2012). We also show the same for the newly identified W2M red QSO spectra
(red filled histogram) which we derived by fitting a single Gaussian profile
to the strongest broad line in each QSO spectrum (§2.5.3). The FWHM of the red
QSO spectra range from 1000 km s-1 to 9000 km s-1. The red vertical line shows
the $v_{\rm FWHM}=1000$ km s-1 cutoff.  Figure 5: Flowchart showing the
process for identifying the sample of broad-line QSOs with SDSS spectra.

###  2.4 KX selection of red QSOs

We are interested in recovering likely reddened QSOs missed by the SDSS and
other QSO selection algorithms by performing spectroscopy on sources lacking a
spectral classification in SDSS. Initially, our follow-up spectroscopy spanned
the full WISE color selection box (Eqn. 1; Figure 2). Although incomplete, we
recovered no new QSOs (red circles) in the lower right region of the diagram.
We therefore aimed to refine our color selection to increase the efficiency of
the follow-up spectroscopy of unidentified sources.

We explore the occurrence of SDSS-identified GALAXYs and QSOs as a function of
$W2-W3$ (i.e., $[4.6]-[12]$) color in Figure 6. On the left, we plot a
histogram of the number of SDSS spectra (gray line), QSOs (blue line), and
GALAXYs (green line) and see a sharp decline of QSO number beyond $W2-W3\sim
3.5$. On the right we show the fraction of QSO-classified spectra in blue,
compared with the fraction of GALAXY-classified spectra in green. The vertical
red dashed line at $W2-W3=3.45$ shows the point where the fraction of GALAXYs
exceeds the QSO fraction. To maximize success of identifying QSOs, we impose
an additional cut:

$\begin{array}[]{l}W1-W2>0.7\\\ \qquad{\rm~{}for~{}}\quad 2<W2-W3<3.45\\\
W1-W2>0.57~{}(W2-W3)-1.27\\\ \qquad{\rm~{}for~{}}\quad
3.45<W2-W3<4.5,\end{array}$ (2)

shown by the thick black line in Figures 2 and 7.

Figure 6: These figures demonstrate the reasoning behind applying a diagonal
color cut at large $[4.6]-[12]$ (i.e., $W2-W3$) colors (Eqn. 2). Left –
Histogram of SDSS-identified objects defined by Eqn. 1 and $[3.4]-[4.6]>0.7$
(grey) and divided by classification: QSOs (blue) and GALAXYs (green). We see
that beyond $[4.6]-[12]\gtrsim 3.5$, GALAXYs begin to overtake QSOs thereby
reducing the efficiency of QSO selection at these colors. Right – Fraction of
SDSS-identified objects classified as QSOs (blue) and GALAXYs (green). The red
line shows the point, $[4.6]-[12]=3.45$, at which the two populations are
equal in proportion.

We note that this cut remains liberal compared to other studies of infrared-
selected AGN (e.g., Stern et al., 2012; Assef et al., 2013, who used
$W1-W2>0.8$) and recovers all but four of the 120 F2M red quasars (as seen in
Figure 1 of Glikman et al., 2018a). The shorter wavelengths of the $W1$ and
$W2$ bands are more affected by reddening than the longer WISE bands, so it is
not surprising that reddened QSOs have slightly redder $W1-W2$ colors.
Applying the color cut in Eqn. 2 leaves 393 objects with no spectrum in SDSS.

Although only 889 Type 1 QSOs with SDSS spectra obey this stricter cut, we
keep the full sample of 1072 QSOs established in §2.3 for our subsequent
study. In Section 5.1.3 we demonstrate that the radio properties, which are
the interest of this paper, are indistinguishable between the full blue QSO
sample and the blue QSOs obeying Eqn. 2 (above the ‘diagonal’ cut).

Figure 7: Here we focus on the sources that obey the selection criteria
defined in Equation 1 with SDSS counterpart. Gray circles are sources without
a spectrum in SDSS. Blue circles are identified as QSO in SDSS, orange circles
are identified as GALAXY, and open white circles are STARs. Red circles are
the objects that meet the selection criteria in Equations 2 and 3 that amount
to the sample that we follow up spectroscopically. Although we defined a
generous selection box around the QSO region, most of our candidates overlap
the SDSS-identified QSOs and the solid black line represents our refined
selection cuts.

The SDSS QSO selection algorithm (Richards et al., 2002; Ross et al., 2012) is
very successful at finding blue, broad-line QSOs whose colors naturally
deviate from the color locus produced by Galactic stars. However, dust-
reddened QSOs often overlap the stellar locus in the optical, making them hard
to find by their SDSS colors alone (Urrutia et al., 2009, Figure 5). As a
remedy, Warren et al. (2000) showed that in optical-to-near-infrared (so-
called KX) color space, reddened QSOs can be cleanly separated from stars. We
apply the following optical-to-near-infrared color cuts (in AB magnitudes),

$J_{AB}-K_{AB}>0.5\quad{\rm and}\quad g_{AB}-J_{AB}>2.0,$ (3)

to the remaining 393 sources using their 2MASS $J$ and $K$ magnitudes combined
with their SDSS $g$ magnitude. These color cuts correspond to
$J_{V}-K_{V}>1.462$ and $g_{AB}-J_{V}>2.938$, which are in line with the cuts
used by Urrutia et al. (2009, $J_{V}-K_{V}>1.3$) and Glikman et al. (2013,
$J_{V}-K_{V}>1.5$).

Figure 8 shows the $J-K$ vs. $g-J$ KX-selection colors, using the same color
scheme as in Figure 7. In this space, the black dashed line separates stars
from QSOs (Maddox et al., 2008). Our candidate sources appear to complete the
cloud of points around $g-J\simeq 2$ but also extend to very red tails in both
$g-J$ and $J-K$. Applying these criteria to the 393 sources results in 91
candidates. Seven of these are in the Stripe 82 region and, of those, three
were originally identified in Glikman et al. (2018a). Another source in Stripe
82 was originally identified in the FIRST-2MASS survey (Glikman et al., 2007)
and was recovered in Glikman et al. (2018a). The 91 selected sources are shown
with large outlined red circles in Figures 7, and 8. They obey the criteria
outlined in Equation 2 at WISE wavelengths and have the optical through near-
infrared colors in Eqn. 3.

Figure 9 shows a flowchart of the process to select these redder QSO
candidates missed by SDSS. We list these candidates in Table 2.5.2, including
their positions, optical through mid-infrared magnitudes, as well as
classification and redshift based on spectroscopic followup (§2.5).

Figure 8: Optical-to-near-infrared colors used in the KX selection method
defined by Warren et al. (2000) and refined by Maddox et al. (2008). The KX
color cut separating stars from QSOs is shown by the black dashed line. The
red dashed line is the cut we impose (Eqn. 3) in this space to identify red
QSOs. The colors of sources are the same as in Figure 7. Figure 9: Flowchart
showing the process for selecting the red QSO candidate sample in this work.
Each box reports the candidate selection step, with the number of sources that
passed each stage shown in boldface. The final box reports the confirmed QSOs
that are added to the Type 1 QSOs with SDSS spectra.

###  2.5 Spectroscopic Classification

####  2.5.1 Archival Spectroscopy

Among the 91 QSO candidates, we recover twelve objects that were identified in
red QSOs surveys. Four were in the pilot study of Glikman et al. (2018a) over
Stripe 82: W2M J0030$-$0027 is a Type 2 AGN at $z=0.242$ originally identified
in Glikman et al. (2012); W2M J0306+0108 is a Type 2 AGN at $z=0.189$; and W2M
J0349+0054 is a Type 2 AGN at $z=0.109$. F2M J2216$-$0054 is a red quasar at
$z=0.2$ originally identified in Glikman et al. (2007).

The other eight objects are in the spring sky region and were all identified
in the F2M survey (Glikman et al., 2012). Four of them are red quasars. F2M
J1232+1112 is at $z=0.25$, F2M J1248+0531 is at $z=0.749$, F2M J1439+1136 is
at $z=0.296$, and F2M J1554+0714, which was unclassified in Glikman et al.
(2012), has had new spectroscopy (§2.5.2) revealing broad H$\alpha$ and
corresponding, narrower, H$\beta$ at $z=1.64$.

####  2.5.2 Spectroscopic Observations

We obtained spectroscopic classifications of all but four out of the 91
candidates in our sample ($96\%$ spectroscopic completeness). We also obtained
near-infrared spectra for QSOs whose SDSS spectrum revealed strong reddening
(§3) to broaden their wavelength coverage. These observations were conducted
over six observing runs at three different telescope facilities. We used the
SpeX spectrograph (Rayner et al., 2003) at the NASA InfraRed Telescope
Facility (IRTF), TripeSpec (Herter et al., 2008) at the Palomar Observatory’s
200 inch Hale telescope, and TripleSpec (Wilson et al., 2004) on the 3.5 meter
telescope at the Apache Point Observatory. The data were reduced using the
Spextool software package (Cushing et al., 2004), which was originally written
for the SpeX instrument but has been modified to also reduce data from
TripleSpec. We followed the procedure outlined in Vacca et al. (2003) to
correct for telluric absorption using spectra of nearby A0V stars obtained
immediately before or after our target exposures.

We also obtained optical spectroscopy of twenty-seven sources. Nine sources
were observed with the MODS1B Spectrograph on the Large Binocular Telescope
(LBT) observatory, with the red and blue arms simultaneously, with a
0$\farcs$6-wide slit on UT 2013 March 14, covering the wavelength range
$3300-10100$Å. We also obtained 18 optical spectra with the KAST spectrograph
at the Lick Observatory. All data were reduced using standard IRAF routines
and flux calibrated using Feige 34.

Examination of the spectra for the presence of a broad emission line (see
§2.5.3) results in 40 new QSOs; 37 of them have $E(B-V)>0.25$ deeming them as
red QSOs (see §3) and label them with the prefix “W2M”, which is an
abbreviation of WISE-2MASS, consistent with our definition established in
Glikman et al. (2018a). The three newly-discovered blue QSOs are included with
the blue sample.

Figures 10a \- 10e show a spectral atlas of the optical through near-infrared
spectra of the 40 objects classified as QSOs in decreasing redshift order. The
last five columns of Table 2.5.2 provides the details on the objects’
spectroscopic observations, source classification, and assigned redshift.

Figure 10a: Optical through near-infrared spectral atlas of newly identified
QSOs, in decreasing redshift order. Typical QSO emission lines are marked with
vertical dashed lines. Orange circles represent the photometry-based fluxes in
the $g$, $r$, and $i$ bands from SDSS, to which the optical spectrum was
scaled, and $J$, $H$, and $K_{s}$ bands from 2MASS, to which the near-infrared
spectrum has been scaled. The red line is the best-fit reddened QSO template
from which $E(B-V)$ is derived.  Figure 10b: Same as Figure 10a.  Figure 10c:
Same as Figure 10a. Figure 10d: Same as Figure 10a. Figure 10e: Same as Figure
10a. Table 2: WISE-2MASS red QSO candidates

Name | $g$ | $r$ | $i$ | $J$ | $K$ | $W1$ | $W2$ | $W3$ | $W4$ | Redshift | Class | Spectrum | References
---|---|---|---|---|---|---|---|---|---|---|---|---|---
| (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) |  |  |  |
J003009.08$-$002744.2 | 20.35 | 19.18 | 18.44 | 16.50 | 14.20 | 13.86 | 12.81 | 8.90 | 6.48 | 0.242 | NLAGN | LRIS | G18
J004331.44$-$001527.3 | 20.20 | 18.74 | 18.19 | 16.86 | 14.42 | 13.27 | 12.34 | 9.74 | 7.39 | 0.294 | QSO | LRIS/IRTF |
J030654.88+010833.6 | 19.39 | 18.05 | 17.53 | 16.36 | 14.49 | 13.22 | 12.02 | 9.19 | 6.30 | 0.189 | NLAGN | LICK/SpeX | G18
J034902.65+005430.5 | 18.97 | 17.92 | 17.14 | 15.84 | 13.35 | 11.59 | 10.23 | 7.31 | 4.86 | 0.109 | NLAGN | LICK/SpeX | G18
J080121.75+004028.4 | 19.10 | 18.10 | 17.56 | 16.16 | 14.51 | 13.66 | 12.79 | 9.54 | 7.23 | … | ? | IRTF |
J080155.37+132827.7 | 22.73 | 18.14 | 17.55 | 16.36 | 14.62 | 13.47 | 12.72 | 9.63 | 7.17 | 0.110 | Galaxy | LICK |

References. — G18 = Glikman et al. (2018a)

References. — Only a portion of this table is shown here to demonstrate its
form and content. A machine-readable version of the full table is available.

####  2.5.3 Red QSO Line Properties

All the red QSOs for which we obtained spectra have at least one broad
emission line to which we fit a single Gaussian profile. We convert the best-
fit $\sigma$, in wavelength units, to a FWHM in velocity units through the
expression,

$v_{\rm FWHM}=c\times\frac{2.355\sigma}{\lambda_{\rm 0}},$ (4)

where $c$ is the speed of light and $\lambda_{\rm 0}$ is the rest-wavelength
of the line being fit. The red filled histogram in Figure 4 shows the
distribution of $v_{FWHM}$ for the newly-discovered red QSOs.

Our master sample of mid-infrared color selected Type 1 QSOs now contains 1112
sources, with 1072 objects coming from the SDSS spectroscopic sample and 40
newly-identified QSOs selected by their red colors in KX color space. Table 3
lists the full QSO sample – referred to as the W2M sample going forward –
listing their coordinates, SDSS, 2MASS, and WISE magnitudes as well as peak
flux densities in the FIRST (1.4 GHz) and Very Large Array Sky Survey
(VLASS;$2-4$ GHz; Lacy et al., 2020) surveys and the spectroscopic redshift.
The table also indicates whether the object obeyed the stricter selection
criteria of Eqn 2, whether the object qualifies as a red QSO (§3) with its
best-fit $E(B-V)$ value.

Table 3: WISE-2MASS QSOs

Name | $g$ | $r$ | $i$ | $J$ | $K$ | $W1$ | $W2$ | $W3$ | $W4$ |  $F_{pk}$aaFIRST peak flux density. |  $V_{pk}$bbVLASS peak flux density. | Redshift | DiagccBoolean indicating whether the object obeys the stricter ’diagonal’ selection cut of Equation 2. | Class | $E(B-V)$ |  $M_{K}$ddAbsolute $K$-band magnitude, de-reddened by $E(B-V)$ for red QSOs. | LumeeBoolean indicating whether the object is part of the luminosity restricted subsamples. | $\log L_{\rm bol}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
| (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mJy) | (mJy) |  |  |  | (mag) | (mag) | Sel | (erg s-1)
W2M J003147.21$-$004359.2 | 17.35 | 16.47 | 15.95 | 15.50 | 14.24 | 13.22 | 12.68 | 9.88 | 8.20 |  |  | 0.067 | F | QSO |  | $-23.15$ | T | 43.90
W2M J003238.19$-$010035.2 | 17.58 | 17.14 | 16.76 | 15.80 | 14.34 | 13.41 | 12.61 | 9.19 | 6.78 |  |  | 0.092 | T | QSO |  | $-23.77$ | T | 44.36
W2M J003328.05$-$001912.8 | 17.37 | 16.54 | 16.01 | 15.29 | 13.92 | 12.72 | 11.96 | 9.41 | 7.07 | 1.02 | 2.02 | 0.107 | T | QSO |  | $-24.55$ | T | 44.59
W2M J003659.78$+$010544.3 | 18.19 | 17.68 | 17.20 | 16.00 | 14.45 | 13.57 | 12.82 | 10.02 | 8.15 |  |  | 0.121 | T | QSO |  | $-24.31$ | T | 44.42
W2M J003659.82$-$011332.5 | 21.29 | 20.19 | 19.66 | 16.57 | 13.63 | 11.67 | 10.51 | 7.88 | 5.21 | 1.92 | 2.77 | 0.294 | T | red-QSO | 0.79 | $-27.67$ | T | 46.25
W2M J003847.97$+$003457.5 | 17.07 | 16.58 | 16.09 | 15.34 | 14.01 | 13.15 | 12.50 | 9.91 | 7.63 | 1.12 |  | 0.081 | F | QSO |  | $-23.81$ | T | 44.12

Note. — Only a portion of this table is shown here to demonstrate its form and
content. A machine-readable version of the full table is available.

##  3 Defining Red and Blue QSO subsamples

Here, we aim to construct well-defined red and blue QSO subsamples whose
properties can be distinguished and compared. We study the reddening
properties of the QSOs and use that information to determine the sample’s de-
reddened absolute magnitudes to ensure the two samples are intrinsically
similar.

###  3.1 Reddening investigation

We fit a reddened QSO template to each SDSS spectrum, following the formalism
described in Glikman et al. (2007) and plot the distribution of $E(B-V)$ in
the left panel of Figure 11. We use the Gordon & Clayton (1998) SMC dust law
(blue line) and, for comparison, we also used the dust law from Zafar et al.
(2015) that was derived directly from QSOs (orange line). The distributions
are nearly identical and we find no systematic differences between the two
dust laws. Therefore, we choose to use the SMC dust law of Gordon & Clayton
(1998) in order to maintain consistency with previous red quasar studies. The
dashed blue line shows the $E(B-V)$ distribution for the QSOs that obey Eqn 2,
which demonstrates that including sources at lower $W1-W2$ values finds
optically redder objects. This is likely due to contamination from host galaxy
emission, since dust reddening would have the effect of increasing the $W1-W2$
color. We discuss the removal of host galaxy emission to account for this in
Section 3.2, below.

Figure 11:  Left – Distribution of $E(B-V)$ for the SDSS QSOs, which we
determined by fitting a reddened QSO template to each spectrum. The solid line
represents the broad line QSOs obeying Eqn 1 with $E(B-V)$ derived using the
SMC dust law of Gordon & Clayton (1998) shown in blue and the Zafar et al.
(2015) dust law shown in orange. The distributions are nearly identical and
therefore we proceed with the SMC dust law. The dashed line is the $E(B-V)$
distribution for the subsample obeying the diagonal color cut of Eqn 2. The
Gaussian curve shown in gray is centered at $E(B-V)=0$ with excess reddening
likely due to host galaxy contamination. Right – Distribution of $E(B-V)$ for
the SDSS QSOs after subtracting host galaxy light from each spectrum (blue).
Solid lines represent the full sample, while the dashed line is for the
subsample obeying the diagonal color cut. The orange histogram shows, for
comparison, the original values plotted in the left panel. Subtracting a host
galaxy shifts the $E(B-V)$ to lower values and is subsequently well-fit by a
Gaussian distribution centered on $E(B-V)=0$. The vertical red dashed line at
$E(B-V)=0.25$ is our defined limit for a QSO to be considered dust-reddened.

The distributions appear as a Gaussian (gray curve; fitted to the data with
$E(B-V)<0.12$), peaked at $E(B-V)=0$ with a broad tail extending toward redder
colors. The Gaussian distribution is attributed to an intrinsic spread of the
power-law continuum slope of un-reddened quasars. The tail extending to higher
$E(B-V)$ values is due to dust-reddening as well as host galaxy contamination.
This was pointed out earlier by Richards et al. (2002) when examining the
relative $g^{*}-i^{*}$ colors of SDSS QSOs in an early data release of the
SDSS survey.

###  3.2 Removal of host galaxy emission

Because this sample contains a majority of sources at low redshift ($z<0.4$),
lower luminosity AGN may have red colors as a result of added light from a
host galaxy rather than reddening of the AGN continuum. To correct for this,
we used the Gas AND Absorption Line Fitting code (GANDALF; Sarzi et al., 2006)
to fit a model host galaxy simultaneously with Gaussian profiles fitted to
specified emission lines to the $z<0.4$ objects with SDSS spectra. We examined
the fits and subtract the best-fit host galaxy model from the spectra when a
good fit is achieved (i.e., the galaxy continuum traces stellar absorption
features around 4000Å). Figure 12 shows two representative examples of spectra
whose reddened QSO template fits were poor and needed to have a host galaxy
template subtracted. The top panel in each column is the original, poor, fit.
The middle panel shows the fits produced by GANDALF, with the host galaxy
spectrum shown in green. The bottom panel shows the galaxy-subtracted spectrum
and its fitted reddened QSO template. In both cases – as well as in the other
368 sources that needed this treatment – removing the galaxy results in an AGN
spectrum with little to no reddening.

We then re-fit the galaxy-subtracted AGN with the reddened QSO template to
minimize host galaxy effects on our $E(B-V)$ estimates. The right panel of
Figure 11 shows this newer distribution in blue (solid line is the full SDSS
sample, the dashed line shows only sources obeying Eqn 2 with the original
distribution, seen also ion the left, shown for comparison in orange). The
galaxy-subtracted AGN reddening distribution is now symmetric and well-fit by
a slightly broader Gaussian, still centered at $E(B-V)=0$. The fact that host-
galaxy removal shifts the QSO subsample that obeys the diagonal color
selection (dashed blue line; Eqn 2 ) to agree with the full sample provides
reassurance that the excess red color in that subsample was indeed due to host
galaxy light and that the blue QSOs are otherwise similar. The removal of the
host galaxy results in a reddening estimate that is lower by a mean of $0.42$
mag.

Figure 12: Examples of two QSO spectra whose host galaxies were fitted with
GANDALF and subtracted to improve the QSO template fitting process and better-
estimate the reddening, parametrized by $E(B-V)$, experienced by the QSO. The
top panel in each column is the original, poor, fit. The black line is the
original spectrum. The red line is the best-fit, reddened QSO template and the
cyan line is the unreddened QSO template. The middle row shows the fits
produced by GANDALF, with the host galaxy spectrum shown in green, atop the
original spectrum in black. The bottom panel shows the galaxy-subtracted
spectrum in black and the newly best-fitted reddened QSO template in red. In
both examples, and in the vast majority of spectra in our sample, host galaxy
emission results in artificially large extinctions. Subtracting the host
galaxy exposes predominantly unreddened QSOs.

###  3.3 Reddening in the W2M sample

The 40 newly identified W2M QSOs do not have the uniform optical spectroscopy
that SDSS provides, and thus our estimates of $E(B-V)$ for these sources
requires an object-by-object approach, depending on what data exist for each
source. Glikman et al. (2007, 2012, 2013, 2018a) discuss the challenges of
estimating $E(B-V)$ depending on the available spectroscopy and photometry.
Here we use wisdom gained from those studies to determine the reddening for
the W2M QSO sample.

Twelve sources have only a near-infrared spectrum and, when fit with a
reddened composite template, six of those yield $E(B-V)\geq 0.25$. The
remaining six are fit with a lower $E(B-V)$ but when the reddened template is
compared with the SDSS photometry, there are significant mis-matches,
suggesting that the optical emission is far more reddened than is constrained
by the the near-infrared spectrum. Some of this is due to the near-infrared
spectra becoming noisier toward shorter wavelength, where the shape of the
spectrum is more sensitive to $E(B-V)$. For the six sources whose reddened
fits do not match their optical photometry, we fit a reddened template to the
effective flux in eight photometric bands, $u$,$g$,$r$,$i$,$z$ from SDSS plus
$J$, $H$, $K_{s}$ from 2MASS, following the procedure that we describe in
Section 6.1 of Glikman et al. (2013). In this process, three QSOs were well
fit with $E(B-V)<0.25$ and we do not classify them as red QSOs.

The two highest redshift QSOs in the W2M sample, W2M J1542+1259 and W2M
J1042+1641777W2M J1042+1641 is also a quadruply lensed quasar with significant
flux anomaly. The system is thoroughly analyzed in Glikman et al. (2018b) and
Glikman et al. (2022, submitted). , both at $z=2.52$, are heavily absorbed in
their optical (rest-frame UV) spectrum due to outflowing gas that gives rise
to broad absorption lines (BAL). For these sources, we conducted the template
fitting to the spectrum considering only wavelengths above $\lambda>8000$Å and
$\lambda>10000$Å, respectively. Details on this fitting approach are provided
in Glikman et al. (2018b)

The remaining sources that had both an optical and near-infrared spectrum were
combined in one of two ways. Objects that had overlapping optical and near-
infrared spectral regions with sufficiently high signal-to-noise were scaled
to match, and the combined spectrum was fitted to a reddened QSO template. In
other cases, where the spectral scaling was either impossible due to no
overlap between the optical and near-infrared spectrum or where the
overlapping regions were so noisy that small differences in the region used
for comparison yielded wildly different results, we scaled the spectra to
their SDSS and 2MASS photometry, combined them into a single spectrum and fit
that spectrum with a reddened QSO template.

We note that in many of the low-redshift QSOs ($z\lesssim 0.3$) the optical
spectrum does not match well the reddened template, often exposing excess
light from the host galaxy. We attempted to remove the galaxy as described in
Section 3.2, but were then unable to reliably combine the galaxy-subtracted
optical QSO spectrum with the infrared spectrum. In one case, W2M J1250+1318
at $z=0.3$, the galaxy-subtracted optical spectrum had sufficient signal in
the observed $g$-band to allow a scaling to that part of the spectrum,
yielding an excellent fit to the reddened template. However, the galaxy-
subtracted spectra for the other sources were not easily scalable to the near-
infrared spectrum. Therefore, since the QSO dominates at longer wavelengths,
we use the near-infrared spectra alone to determine $E(B-V)$.

We find 37 W2M QSOs with $E(B-V)>0.25$ which we classify as red QSOs, and 3
W2M QSOs with $E(B-V)<0.25$ which we classify as blue QSOs. This information
is recorded in a column in Table 3.

###  3.4 Red QSOs in the SDSS sample

Red QSOs also exist among the sources with SDSS spectra and we wish to
identify them and combine them with the 37 newly-identified red QSOs. As seen
in Figure 8, there exist SDSS-identified QSOs in the KX-selection box along
with the newly identified W2M red QSOs. All but three of them obey the
diagonal color cut (Eqn 2) and 28 of them meet the broad-line criterion. One
source, J131327.46+145338.5, whose SDSS spectrum failed the broad-line
criterion is a known F2M red quasar at $z=0.584$ (F2M1313+1453 Glikman et al.,
2012) where strong broad H$\alpha$ and Pa$\beta$ lines are seen in its near-
infrared spectrum. The object was assigned a redshift of $z=5.513$ by SDSS
based on a mis-identification of [O III] $\lambda 5007$ as Ly$\alpha$, which
is a common occurrence in automated redshift assignments of red QSOs (see,
e.g., §3.1 of Glikman et al., 2018a).

We plot the $E(B-V)$ values for these KX-obeying sources as filled blue bars
in the histogram shown in Figure 11. After correcting for the host galaxy,
many of them recover their blue colors, but 19 of them have $E(B-V)>0.25$,
meeting our red QSO definition. Another 35 sources have $E(B-V)>0.25$, based
on their spectral fits, but their KX colors are outside the selection criteria
of Eqn 3. Most of these sources are at low redshift ($z<0.3$) and have
$E(B-V)<0.3$, lying in the tail of the reddening distribution. Visual
inspection of their spectra confirm that the host-galaxy subtracted AGN
spectrum does not resemble the AGN-dominated red systems we are pursuing and
they are not included in the red QSO counts. However, five of these sources
are at higher redshifts ($z>0.4$) and their spectra are well fit by a reddened
QSO template. Their colors are outside of the KX-selection box either due to
strong emission lines affecting their photometry or variability between the
2MASS and SDSS epochs yielding artificially blue colors.

We also matched the SDSS QSOs to the F2M red quasar sample of Glikman et al.
(2012) and identified 11 additional sources. Two of these had $E(B-V)<0.25$
based on the SDSS spectrum alone, but when combined with their near-infrared
spectrum from Glikman et al. (2012), the fit yielded $E(B-V)>0.25$ which
supersedes the SDSS-derived value. The other nine F2M sources were already
recognized as red quasars through our various methods above. Therefore, we add
the 19 SDSS QSOs in the KX selection box plus the five $z>0.4$ SDSS QSOs with
$E(B-V)>0.25$ and two F2M quasars with SDSS spectra to our sample of 37 red
QSOs, as they obey the same selection criteria, including the infrared color
conditions of Eqn 2, as the W2M sources listed in Table 2.5.2. The final
number of red QSOs that we identify in this survey is thus 63 (37+19+5+2),
leaving 1049 blue QSOs.

We note that this number of red QSOs is a lower limit, since near infrared
spectroscopy can reveal broad lines in objects that show narrow emission their
optical spectra (e.g., F2M1313+1453). We also miss red QSOs whose flux is
dimmed below our $K<14.7$ mag limit when reddened with $E(B-V)>0.25$, but
whose unreddened counterparts are present in the blue QSO sample.

##  4 Luminosity-Restricted sample

The host-galaxy subtraction (§3.2) allowed us to compute the reddening
experienced by the QSOs without contamination from stellar light. We used the
distribution of those reddenings (Figure 12, right) to define a sample of
reddened QSOs with $E(B-V)>0.25$. The resultant $E(B-V)$ values can then be
used compute de-reddened absolute $K$-band magnitudes, $M_{K}$. At the same
time, since all the QSOs have a measured WISE $W4$ (22$\mu$m) magnitude – a
wavelength at which extinction is negligible – we can check whether $W4$
absolute magnitude and de-reddened absolute $K$-band magnitude are consistent.

In the left panel of Figure 13, we plot the de-reddened, absolute $K$
magnitude versus absolute $W4$ (22$\mu$m) magnitude, in the observed frame,
with no reddening correction. The unreddened QSOs are shown with blue circles,
with the QSOs that had host-galaxy subtraction overplotted with orange
squares. The orange squares lie atop the blue points with no apparent offset,
which provides reassurance that the reddenings derived after host galaxy
removal are reliable (i.e., had we correccted their $K$-band magnitudes by the
erroneous $E(B-V)$ values assumed from the poor spectral fits, the orange
points would have been shifted toward higher luminosities). The red circles
are the red QSOs whose $K$-magnitudes are de-reddened and they lie along the
same relation as all the other QSOs, providing further assurance that the
reddening analysis of our sample is reliable.

In the right panel of Figure 13, we plot the $W4$ and de-reddened $K$-band
absolute magnitudes versus the bolometric luminosity, which we compute by
interpolating between the infrared fluxes in the WISE bands, corrected for
reddening in the red QSOs, to measure the rest-frame $6\mu$m luminosity and
apply a bolometric correction factor of 7.82 using the QSO spectral energy
distribution (SED) from Richards et al. (2006). The $W4$ magnitudes are
plotted in the upper relation and the de-reddened $K$-band magnitudes are
shown in the lower relation. Both relations correlate well, though there is
some scatter in the high luminosity $K$-band measurements.

Figure 13:  Left – De-reddened absolute $K$-band magnitude from 2MASS versus
absolute 22$\mu$m ($W4$) magnitude, computed in the observed frame, for the
QSO sample. A direct relationship is seen across the luminosity range, except
at the highest luminosities where the $K$-band magnitudes are slightly
enhanced compared with $W4$. Overplotted are the sources whose host galaxy
light was subtracted to determine the QSO reddnenings (§3.2; orange squares).
The red circles are the red QSOs whose absolute $K$-band magnitude have been
de-reddened. The red QSOs lie on the same relation with no apparent offset.
Right – Absolute magnitude in $W4$ (upper locus) and de-reddened $K$ (lower
locus) magnitudes versus the QSO bolometric luminosity. Both relations
correlate well, although there is more scatter in the $K$-band relation at
higher luminosities. These agreements provide reassurance that the reddening
estimates in Section 3 are reliable.

Although our use of WISE colors to select QSOs avoids the color bias
experienced by optically selected QSO samples that miss most reddened QSOs,
the imposition of the $K<14.7$ mag limit introduces a luminosity bias because,
even at 2.2$\mu$m, moderately reddened QSOs will be dimmed below the flux
limit and therefore missed unless they are intrinsically more luminous. The
bias has a strong redshift dependence because, at higher redshifts, 2.2$\mu$m
represents rest-frame optical emission, which is more sensitive to the effects
of dust-extinction. Therefore, to conduct a valid comparison between blue and
red populations, we must compare objects with similar intrinsic absolute
magnitude thresholds at all redshifts.

We plot the absolute $K$-band magnitudes for the sample in Figure 14, where
the red QSOs are de-reddened and colored by their $E(B-V)$ values versus
redshift. The dashed lines indicate the flux limit of $K<14.7$ mag with no
reddening and with increasing amounts of extinction. Here, the luminosity bias
becomes obvious with an absence of red QSOs near the flux limit especially
toward increasing redshifts. This trend is well-explained by the fact that
higher-redshift sources with even moderate amounts of reddening must be more
luminous in order to pass the selection criteria. At lower redshifts, the
effect of reddening is weaker and we identify heavily reddened sources with
intrinsic luminosities consistent with their blue QSO counterparts.

A similar analysis on the F2M red quasars revealed that they are more luminous
than their unobscured counterparts, when corrected for extinction (Glikman et
al., 2007, 2012). In the W2M sample, we see that at lower redshifts ($z<1$)
the red and blue QSOs occupy similar luminosities, while at higher redshifts
$z>1$ there appear to be more high luminosity red QSOs compared with blue
QSOs.

Figure 14:  Absolute $K$-band magnitude from 2MASS, where the red QSOs have
been de-reddened using $E(B-V)$ determined in §3 and color coded to reflect
the original value, defined in the legend. The dashed lines indicate the
survey limit of $K<14.7$, and for increasing amounts of extinction. The small
dots are the blue QSOs.

In order to define a blue QSO subsample that has a similar intrinsic
luminosity limit, we plot in Figure 15 the difference between the flux limit
(dashed line labeled $E(B-V)=0$ in Figure 14) and the de-reddened $K$-band
magnitude. Open circles are blue QSOs while filled red circles are red QSOs.
We also plot the flux limit for sources with $E(B-V)=0.25$ with a dotted line
and use this as a cut to separate out the lower luminosity blue QSOs that have
no red QSO counterparts. Imposing this luminosity cut leaves 798 blue QSOs
which we plot as blue open circles in Figure 15 (for a total of 861 QSOs, when
including the 63 red QSOs). We use this sample to compare the fraction of red
and blue QSOs and their radio properties in the sections that follow. We list
in the last three columns of Table 3 the de-reddened absolute $K$-band
magnitude, $M_{K}$, the bolometric luminosity, $L_{\rm bol}$, and a boolean
indication of whether the object is part of the luminosity restricted
subsample.

Figure 15: The magnitude difference between the $K=14.7$ magnitude limit of
the W2M survey minus the de-reddened absolute $K$ magnitude for red and blue
QSOs showing a luminosity bias imposed by the flux limit on reddened QSOs. De-
reddened red QSOs are plotted with red filled circles. Blue QSOs are shown
with open symbols. The dotted black line indicates the magnitude limit when
$E(B-V)=0.25$ such that blue QSOs above the line (blue circles) have
comparable absolute magnitudes to red QSOs.

We note that two of the most luminous infrared sources seen at $z\simeq 2.5$
are both gravitationally lensed QSOs of the rare quad variety. The blue QSO is
the well-studied Cloverleaf Quasar H1413+117 (Kayser et al., 1990) while the
red QSO was newly-discovered in the current survey and is fully analyzed in
(Glikman et al. 2018b). We exclude these sources from the statistical analyses
that follow.

Figure 16 shows the redshift distribution of the Type 1 SDSS QSOs (in blue)
along with the newly added W2M red QSOs (in red). The left panel shows a
histogram of all sources on a logarithmic scale to better view both the blue
and red populations. On the right we normalize the blue and red histograms by
the total number of QSOs in each subsample so that their redshift
distributions could be better compared. There is a dearth of red QSOs at low
redshifts ($z<0.1$). A similar observation was made in Glikman et al. (2018a)
for the pilot sample of luminous infrared selected red and blue QSOs in Stripe
82 where it was shown that Type 2 (narrow line) QSOs dominated the obscured
QSO population at low redshifts, while the numbers of red QSOs increased at
higher redshifts. This behavior is also seen in X-ray selected red QSOs,
suggesting an evolutionary explanation (LaMassa et al., 2017). On the other
hand, The absence of $z<0.1$ red QSOs may be a selection effect, as these will
be lower luminosity AGN such that, when further reddened, will have their mid-
infrared fluxes contaminated by host galaxy light, potentially shifting them
out of the selection box. With this in mind, we also consider comparisons
between red and blue QSOs restricted to $z>0.1$, which amounts to 645 blue
sources (for a total of 708 QSOs, when including the 63 red QSOs). Although
the red QSO sample is otherwise predominantly low redshift ($z\lesssim 1.0$),
which is largely due to the shallow $K$-band flux limit imposed on our
selection, higher-redshift red QSOs appear to be more represented among the
red QSOs than in the blue population. This too could point to an evolutionary
explanation.

Figure 16: Left – QSO redshift distribution plotted on a logarithmic scale.
The solid blue line represents the full set of 1112 QSOs, while the dashed
blue line shows the 861 luminosity-restricted sources. The filled red
histogram overplots the 63 red QSOs while the dark-red shading represents the
37 newly identified W2M red QSOs (§3.3). Right – Redshift distribution of the
QSO subsamples normalized by the total number of objects in each sample,
plotted on a linear scale. Here, the solid and dashed blue lines are the blue
QSOs and luminosity-restricted blue QSOs, respectively. The red QSOs are shown
in red. The redshift distributions of the blue QSOs are similar for the full
and luminosity-restricted samples, while the red QSOs are shifted such that
they have a higher relative incidence at higher redshifts.

Figure 17: Bolometric luminosity distribution of W2M QSOs. Left – The full QSO
sample. Right –  The FIRST-detected subsample. The gray filled histogram shows
all W2M sources. The light blue histogram represents the luminosity restricted
subsample of blue QSOs, and the dark blue shows blue QSOs with $z>0.1$. The
red histogram represents the red QSOs.

Figure 17 shows a histogram of the bolometric luminosities for the W2M QSOs,
divided into various subsamples. The left panel shows all the the W2M QSOs,
while the right panel shows only the FIRST-detected sources. In both panels,
the gray shaded histogram represents all the QSOs within that category, and
reflects the full luminosity range accessible to our selection criteria. The
blue histogram shows the full luminosity restricted blue QSO sample, and the
dark blue shaded subset represents only QSOs with $z>0.1$. We use these
luminosity-restricted subsamples to compare red and blue QSOs in the sections
that follow. However, we caution that because we are unable to create
perfectly luminosity-matched samples, owing to the fact that we do not know
how many red QSOs we are missing as they drop below our selection threshold,
some of our results could be explained by differences in the luminosity
functions of the two samples.

We note that this luminosity restricting effort is the minimum needed to
establish comparable populations. The blue QSOs excluded by this process
removes objects that would fall below the flux limit if reddened by
$E(B-V)=0.25$. However, the red QSO sample is still incomplete in ways that
cannot be easily corrected without knowing the distribution of $E(B-V)$ as a
function of luminosity and redshift. The lower luminosity bins of the red QSO
histogram are therefore incomplete and those sources can only be recovered
with a deeper survey.

##  5 Results

We have constructed a carefully-selected, luminosity-restricted sample of 63
red and 798 blue Type 1 QSOs based on their mid and near-infrared properties.
All of these sources are spectroscopically confirmed either through publicly
available spectra from SDSS or from supplementary spectroscopy. With this
sample of QSOs in hand, we study their radio properties and demographics in
the sections below.

###  5.1 Radio Properties

A recent study by Klindt et al. (2019) investigated the fraction of FIRST-
detections among blue and red QSOs in SDSS, defined according to how their
observed $g-i$ colors compared to the median QSO $g-i$ at a given redshift
(this is effectively the ‘relative color’ defined in Richards et al., 2003)
and found that red QSOs have a FIRST-detection fraction that is $\sim 3$ times
higher than blue QSOs, across redshifts. A follow-up study by Fawcett et al.
(2020) used stacking of the radio images of red and blue QSOs, identified in a
similar way but that lacked a radio detection, and found that the integrated
flux density of the median red QSO is 30% higher than the median blue QSO’s
integrated flux density.

The SDSS spectroscopic survey is not sensitive to heavily reddened QSOs,
therefore the sample defined as ‘red’ by Klindt et al. (2019) is dominated by
objects with $E(B-V)\sim 0.05-0.2$ and contains very few sources with
$E(B-V)>0.25$. With the W2M sample, we have taken a more conservative approach
in separating reddened QSOs from blue QSOs, having shown that most QSOs have a
natural distribution of spectral slopes that, because of their power-law
shape, can mimic the presence of dust (or ‘negative dust’ for very blue
objects) when a template QSO spectrum is fitted to a host-galaxy-subtracted
AGN spectrum.

If the differences in radio properties seen in Klindt et al. (2019) and
Fawcett et al. (2020) are best-explained by dust-reddened QSOs being
fundamentally distinct from ‘normal’ QSOs then those differences should be
more pronounced when comparing radio properties between blue and red QSOs in
the W2M sample.

####  5.1.1 radio detected fraction

We matched the luminosity-restricted W2M sample to two radio surveys that
overlap our fields. There were 249 matches within 1$\farcs$5 to the FIRST
survey which reaches a 5-$\sigma$ sensitivity of 1 mJy at 1.4 GHz and has an
angular resolution of 5″. We found 186 matches to VLASS, which reaches a
median rms of 120 $\mu$Jy at $2-4$ GHz with a spatial resolution of
2$\farcs$5, using the catalog constructed by Gordon et al. (2021).

Figure 18, left, shows the fraction of FIRST-detected sources among the red
and blue QSO subsets as a function of redshift. We divided the redshift range
of our sample into four bins each representing 2.85 Gyrs in lookback time
corresponding to redshift limits of $z=0.248,0.612,1.267,3.271$. We confirm
the finding from Klindt et al. (2019) that red QSOs have a significantly
higher fraction of FIRST-detections and that the fraction increases toward
lower redshifts. However, we find a higher overall fraction for both red and
blue QSOs; Klindt et al. (2019) find $\sim 7\%$ and $\sim 17\%$ for blue and
red QSOs, compared to $\sim 28\%$ (226/798) and $\sim 56\%$ (35/62) for blue
and red QSOs in the W2M sample, respectively. The right panel shows the same
for VLASS, largely corroborating the FIRST behavior. We attribute this overall
increase in the detection fraction to be due to the different flux limits of
the two samples, as it is known that more luminous QSOs in general have a
higher radio-loudness fraction (Lacy et al., 2001). While the red QSOs have a
higher detection fraction in the lower redshift bins, at higher redshifts the
difference between the red and blue fractions cannot be distinguished due to
the small numbers of sources in those bins. The breakdown of FIRST- and VLASS-
detected sources for the blue and red subsamples is provided in Table 4.

Figure 18: The fraction of FIRST-detected (left) and VLASS-detected (right)
QSOs as a function of redshift for the red QSOs (red symbols) and blue QSOs
(dark blue symbols). Each redshift bin spans 2.85 Gyr in lookback time. At low
redshifts, the red QSOs have a significantly higher fraction compared to all
the sources we deem unreddened ($E(B-V)<0.25$). However, at higher redshifts,
the radio fraction of red QSOs declines and cannot be distinguished from the
fraction of blue QSOs due to the small sample sizes in those bins. These
trends are similar in both VLASS and FIRST. The error bars here are computed
using binomial proportion confidence interval using a Wilson (1927) interval
for small number counts.  Table 4: Radio detection fraction in redshift bins
for blue and red QSOs, as shown in Figure 18.

Sample | $N_{\rm all}$ | $N(0.023\leq z<0.248)$ | $N(0.248\leq z<0.612)$ | $N(0.612\leq z<1.267)$ | $N(1.267\leq z<3.271)$
---|---|---|---|---|---
Parent QSO | 858aaWe remove one gravitational lens and two sources lacking coverage in FIRST from the total sample of 861. | 577 | 185 | 37 | 79
Color-selected QSO samples
blue QSOs | 798 | 577 | 169 | 36 | 14
red QSOs | 62 | 14 | 34 | 7 | 7
FIRST blue | 226 | 155 | 50 | 15 | 6
FIRST red | 35 | 7 | 21 | 4 | 3
VLASS blue | 167 | 107 | 42 | 13 | 5
VLASS red | 29 | 5 | 19 | 3 | 2
| Radio-detected fraction (%)
|  | FIRST | VLASS
Redshift |  | blue | red | blue | red
$0.023\leq z<0.248$ |  | $26.9^{+1.9}_{-1.8}$ | $50.0\pm 12.9$ | $18.5^{+1.7}_{-1.6}$ | $35.7^{+13.4}_{-11.5}$
$0.248\leq z<0.612$ |  | $29.6^{+3.6}_{-3.4}$ | $61.7^{+7.9}_{-8.6}$ | $24.8^{+3.5}_{-3.2}$ | $55.9^{+8.2}_{-8.6}$
$0.612\leq z<1.267$ |  | $41.7^{+8.3}_{-7.9}$ | $57.1^{+16.6}_{-18.4}$ | $36.1^{+8.3}_{-7.5}$ | $42.9^{+18.4}_{-16.6}$
$1.267\leq z<3.271$ |  | $42.9^{+13.3}_{-12.3}$ | $42.9^{+18.4}_{-16.6}$ | $35.7^{+13.4}_{-11.5}$ | $28.6^{+18.9}_{-13.5}$

####  5.1.2 radio morphologies

A major result from Klindt et al. (2019) was that when categorized by their
FIRST morphology among FIRST-detected sources, a much larger fraction of red
QSOs had a compact appearance compared to blue QSOs (7% vs 2% of the entire
subsample, respectively) while their fractions of extended sources are
approximately the same. A similar finding was reported in Fawcett et al.
(2020).

We examined the FIRST morphologies of our red and blue samples following the
same morphological classifications that Klindt et al. (2019) employed (Faint,
Compact, Extended, FRII-like). The left panel of Figure 19 shows the fraction
of red and blue QSOs with FIRST detections in the different morphological
categories. We see a similar trend with a higher detection fraction of faint
compact red QSOs compared to blue QSOs (32% and 16%, respectively).

We note that Klindt et al. (2019) find that red QSOs have a compact morphology
fraction that is $\sim 3.5$ times higher than for blue QSOs, whereas the W2M
red QSOs have a compact morphology fraction that is only $\sim 2$ times higher
than for blue QSOs. This may be due to the fact that the red QSOs in our
sample are significantly redder, and that some of the sources that we define
as blue, with $E(B-V)\sim 0.05-0.2$, would be categorized as red in the Klindt
et al. (2019) sample. In addition, there may be luminosity-based effects as
there are more radio-loud blue QSOs overall in the W2M sample.

We repeat this exercise in the right panel of Figure 19 with the VLASS images,
which sample the sources at a higher frequency than FIRST. With both surveys,
we separate the ‘Faint’ from ‘Compact’ categories using a flux cutoff of
$S_{\rm pk}<3$ mJy to maintain consistency with Klindt et al. (2019). Table 5
lists the percentages for each morphology class in the red and blue QSO
samples for FIRST and VLASS.

Table 5: Percentage of QSOs by morphology  Sample | Faint | Compact | Extended | FR II-like
---|---|---|---|---
| (%) | (%) | (%) | (%)
FIRST
blue | $16.4_{-1.3}^{+1.4}$ | $10.8_{-1.0}^{+1.1}$ | $1.0_{-0.3}^{+0.4}$ | $0.9_{-0.3}^{+0.4}$
red | $32.3_{-5.6}^{+6.1}$ | $24.2_{-5.0}^{+5.8}$ | $<1.6$ | $<1.6$
VLASS
blue | $11.3_{-1.1}^{+1.2}$ | $10.8_{-1.0}^{+1.1}$ | $0.9_{-0.3}^{+0.4}$ | $0.6_{-0.2}^{+0.3}$
red | $25.8_{-5.1}^{+5.9}$ | $19.4_{-4.5}^{+5.5}$ | $1.6_{-1.0}^{+2.5}$ | $<1.6$

Figure 19: The percentage of red and blue QSO samples with different radio
morphologies, as defined by Klindt et al. (2019). Left – FIRST-detected QSOs
in our survey. Right – VLASS-detected QSOs. The four morphological categories
are illustrated by examples at the top of the figure, where the Faint category
is distinguished from the Compact category by $F_{\rm pk,FIRST}<3$ mJy. The
fraction of red QSOs in the ‘Faint’ and ‘Compact’ categories is approximately
twice as high as for the blue QSOs in the FIRST-detected sample and more than
ten times higher in the VLASS detected sample. Due to the overall small number
of red QSOs (62) we detect no (or few) ‘Extended’, or ‘FRII-like’ sources.

####  5.1.3 median radio properties via stacking

The majority of the QSOs in our sample are undetected in neither FIRST nor
VLASS. We therefore employ the method of image stacking to study the median
radio properties of the red and blue QSOs in our sample following White et al.
(2007). Image stacking involves creating a three-dimensional image cube made
up of image cutouts centered on the individual QSO positions and collapsing it
onto a two-dimensional image with each pixel having the median value along the
stacked axis. The stacking process involves adding signal from sources well
within the noise while reducing the background rms. The brightness of the
resultant stacked image represents the median flux density of the input
sample. The measured flux density of stacked images of sub-threshold FIRST
sources experience a flux density bias (dubbed “snapshot bias” in White et
al., 2007) whereby the stacked source has a $\sim 30\%$ deficit in its
measured flux. Sources that are detected with a peak flux density above 0.75
mJy also have a deficit of 0.25 mJy beam-1, known as the CLEAN bias. As a
result, White et al. (2007) carefully calibrated the peak flux densities
derived from stacking and established a bias correction formula: $S_{p,{\rm
corr}}={\rm min}(1.40S_{p},~{}S_{p}+0.25~{}{\rm mJy})$.

When performing the stacking with the VLASS data, we include only a single
epoch observation for each QSO so that each source has a uniform weighting. We
use the first epoch data if available and only use a second epoch observation
if not. Because the broad frequency range of the VLASS data ($2-4$ GHz)
results in minimal side-lobes, which are the target of the CLEAN algorithm,
they do not suffer from CLEAN bias the way FIRST sources do (Rau et al.,
2016). To investigate whether a “snapshot” bias correction is needed, we
stacked the positions of 12 unresolved COSMOS sources, observed with the same
$2-4$ GHz band with the VLA, whose fluxes were between 0.2 and 0.5 mJy, with
an average of 0.306 mJy (Smolčić et al., 2017). We measure a peak flux density
to their VLASS image stack of 0.289$\pm$0.40 mJy, which is consistent within
the errors. However, given that this test was only conducted on a single stack
of a small number of sources, there may be a correction that is up to $\sim
20\%$ of the stacked value if the uncertainty to the stacked flux is
considered.

As part of the analysis of FIRST image stacking, White et al. (2007) found
that redder QSOs, as parametrized by their $g-i$ relative color, have higher
median radio fluxes. Compared to the SDSS composite, objects redder by 0.8 mag
have radio fluxes that are $\gtrsim 3$ times higher. The median radio flux
densities of QSOs bluer than the mean are flat and do not change with color.

We stacked the FIRST images of the 62 red and 796 blue QSOs888Excluding the
two gravitational lenses, as noted in §4. and examined their median radio
image properties. The top row of Figure 20 shows the resultant blue (left) and
red (right) median stacked images. The scaling of the color bar is matched for
both images, reflecting the result that red QSOs have a higher median flux
density at 1.4 GHz. After fitting a two dimensional Gaussian profile to the
image and applying the bias correction to the peak flux density, we find that
the blue QSOs have $S_{p,{\rm blue}}=0.390\pm 0.008$ mJy and red QSOs have
$S_{p,{\rm red}}=1.183\pm 0.024$ mJy. This is a significant difference in
median flux density, and is consistent with the results from the previous
sections showing a higher incidence of FIRST detections among the red QSOs.
Fawcett et al. (2020) also found a higher flux density for red QSOs but with
significantly lower luminosities and dominated by sources with much smaller
$E(B-V)$ values compared with the W2M red QSOs. They find that the red QSO
stack is 35% brighter than the blue QSO stack.

Similarly, we performed image stacking using cutouts from the VLASS survey and
show the stacked red and blue QSO images in the bottom row of Figure 20. With
a median peak flux density of $S_{p,{\rm red}}=0.698\pm 0.022$ mJy, the red
QSOs have a higher flux density than the blue QSOs, $S_{p,{\rm blue}}=0.301\pm
0.006$ mJy, at higher frequencies as well.

Figure 20: Stacked radio images of blue (left) and red (right) QSOs in our
survey. Top –  FIRST. Images are 33 pixels on a side (59.4″) and are displayed
with a linear color scale. Bottom –  VLASS. Images are 61 pixels on a side
(61″) and are displayed with a asinh color scale. In both cases, the
brightness limits are fixed to the min/max of the blue QSO stack so that the
brighter average flux density of the red QSO is apparent.

To test for biases among the blue QSOs in our sample, we conducted additional
median stacks of several subsets. Given that the ‘box’ selection contains
sources with significant host galaxy emission in their optical spectra, we
constructed a median stack of only QSOs that obey the ‘diagonal’ selection to
compare with the blue QSO stack in Figure 20. Figure 21 shows the resultant
stacks, FIRST in the top row and VLASS in the bottom row, with the blue QSOs
from the ‘box’ selection (Eqn. 1) shown on the left (same as in Figure 20) and
the ‘diagonal’ selection (Eqn. 2) shown on the right. The bias-corrected peak
flux densities in FIRST are $S_{p,{\rm blue}}=0.390\pm 0.008$ mJy and
$S_{p,{\rm blue,diag}}=0.402\pm 0.008$ mJy, and for VLASS are $S_{p,{\rm
blue}}=0.301\pm 0.006$ mJy and $S_{p,{\rm blue,diag}}=0.323\pm 0.007$ mJy,
both of which are near identical. This provides reassurance that we can
compare the red QSOs, which are shifted toward higher $W1-W2$ colors, to blue
QSOs that extend to lower $W1-W2$ colors without biasing our analysis of their
radio properties.

Figure 21: Stacked radio images of blue QSOs obeying the box (left; Eqn. 1)
and diagonal (right; Eqn. 2) selected QSOs in our survey. Top –  FIRST. The
images are 33 pixels on a side (59.4″) and are displayed with a linear color
scale. Bottom –  VLASS. Images are 61 pixels on a side (61″) and are displayed
with a asinh color scale. The stacked images are near-identical. In both
cases, the brightness limits are fixed to the min/max of the box blue QSO
stack on the left.

Since all of the red QSOs are at $z>0.1$, we also stacked the luminosity-
restricted blue QSOs with $z>0.1$ to ensure that the 153 low-redshift and low-
luminosity blue QSOs are not skewing the median radio flux density to lower
values. We find that the median radio flux densities of the $z>0.1$ blue QSOs
have $S_{p,{\rm blue,highz}}=0.373\pm 0.007$ mJy in FIRST and $S_{p,{\rm
blue,highz}}=0.293\pm 0.007$ mJy in VLASS which are both consistent with – and
even slightly fainter than – the full stacked sample.

Considering that the red QSOs are skewed toward higher luminosities we also
stacked only blue QSOs with $L_{\rm bol}>10^{45}$ erg s-1 which overlaps the
histogram of dereddened $K$-band magnitudes for the red QSOs (see Figure 17)
to ensure that the lower-lumniosity blue QSOs are not skewing the median radio
flux densities. This subsample contains 424 such sources (considering only
those that are luminosity-restricted). Here too, we find that the median radio
flux densities are largely consistent with the full stacked sample. The median
FIRST flux density of the high luminosity blue QSOs is $S_{p,{\rm
blue,highlum}}=0.435\pm 0.009$ mJy. The median VLASS flux density is
$S_{p,{\rm blue,highlum}}=0.336\pm 0.009$ mJy. Both values are slightly higher
than the full blue QSO sample stack, but still significantly lower than the
red QSO value.

Finally, we broke the blue sample into three bins of $E(B-V)$ to explore
whether QSOs that are not dust-reddened, but whose intrinsic continua vary to
give slightly flatter and steeper slopes, might exhibit different median radio
properties. To test this, we divided all 796 blue QSOs into three equal sized
bins999Here, as noted in §3, $E(B-V)$ is merely a proxy for the intrinsic
spectral shape of the QSO continuum.. The bins span $E(B-V)<-0.05$ for the
first bin (Bin 1), $-0.05\leq E(B-V)<0.07$ for the second bin (Bin 2),
encompassing the peak of the distribution, and $0.070\leq E(B-V)$ (Bin 3).
Figure 22 shows the median stacked images of these three subsets, with FIRST
at the top and VLASS on the bottom. The bias-corrected FIRST peak flux
densities for the subsets are $S_{p,{\rm Bin1}}=0.404\pm 0.012$ mJy,
$S_{p,{\rm Bin2}}=0.375\pm 0.011$ mJy, $S_{p,{\rm Bin3}}=0.398\pm 0.012$ mJy,
and $S_{p,{\rm Bin1}}=0.268\pm 0.011$ mJy, $S_{p,{\rm Bin2}}=0.316\pm 0.011$
mJy, $S_{p,{\rm Bin3}}=0.314\pm 0.011$ mJy in VLASS. Regardless of intrinsic
color, the overall the median radio flux densities are largely very similar
and significantly lower than the average red QSO flux density.

Figure 22: Stacked radio images of blue QSOs broken into three bins spanning
the blue side, peak, and red side of the $E(B-V)$ distribution seen in Figure
11, right, which we interpret as a proxy for the intrinsic distribution of QSO
continuum slopes rather than dust reddening. Top –  FIRST. The images are 33
pixels on a side (59.4″) and are displayed with a linear color scale. Bottom –
VLASS. Images are 61 pixels on a side (61″) and are displayed with a asinh
color scale. In both cases, the brightness limits are fixed to the min/max of
the blue QSO stack in the middle image. The stacked images are very similar
suggesting that unreddened QSOs do not have significant intrinsic variations
in their radio properties with the shape of their optical spectrum.

In addition, there are differences among the average morphologies of the
various QSO subsets based on the widths along the x and y axes derived from
the two-dimensional Gaussian fits to the blue and red stacked images. We find
that the red QSOs are more compact than the blue QSOs, with FWHMred = 3.1
pixels and FWHM${}_{\rm blue}\sim$4.1 pixels101010The size of a pixel in the
FIRST images is 1$\farcs$8.. Similarly, in VLASS, the red quasars have
FWHM${}_{\rm red}=3.4$ pixels versus FWHM${}_{\rm blue}\sim$ 4.0
pixels111111The size of a pixel in the VLASS images is 1″.. In both cases, the
sources are slightly resolved (FIRST resolution is 2.8 pixels, VLASS
resolution is 3 pixels). We note, that stacking will cause some morphological
distortion by spreading out the PSF depending on where the source’s true peak
exists within the central image pixel. Nonetheless, we find consistency with
the earlier findings (§5.1.2) that, among the radio-detected sources, red
quasars have more compact morphologies.

Table 6 summarizes all the stacking results from the various sub-samples and
their derived quantities. Our main takeaway is that red QSOs display
significantly enhanced radio emission compared with the un-reddened sample,
particularly when red QSOs are defined as having $E(B-V)>0.25$, beyond the
normal spread of spectral shapes.

Table 6: Statistics of stacking results

QSO | Num | $F_{pk}$ | FWHMaaA single FWHM is reported as the average FWHM along the x and y directions. | Num | $F_{pk}$ | FWHMaaA single FWHM is reported as the average FWHM along the x and y directions.
---|---|---|---|---|---|---
Subset |  | (mJy) | (pixels) |  | (mJy) | (pixels)
| FIRST | VLASSbbAn additional $20\%$ uncertainty may be added to the VLASS stacks due to possible “snapshot bias” (see §5.1.3).
Red | 62 | $1.183\pm 0.024$ | 3.1 | 62 | $0.698\pm 0.022$ | 3.4
Blue | 796 | $0.390\pm 0.008$ | 4.1 | 790 | $0.301\pm 0.006$ | 4.1
Blue Diag | 640 | $0.402\pm 0.008$ | 4.0 | 636 | $0.323\pm 0.007$ | 3.9
Blue, $z>0.1$ | 643 | $0.373\pm 0.007$ | 4.1 | 639 | $0.293\pm 0.007$ | 4.0
Blue, $L_{\rm bol}>10^{45}$ | 424 | $0.435\pm 0.009$ | 4.1 | 421 | $0.336\pm 0.009$ | 3.9
Bin 1 | 259 | $0.404\pm 0.012$ | 4.0 | 258 | $0.268\pm 0.011$ | 4.6
Bin 2 | 278 | $0.375\pm 0.011$ | 4.2 | 275 | $0.316\pm 0.011$ | 3.8
Bin 3 | 259 | $0.398\pm 0.012$ | 4.1 | 257 | $0.314\pm 0.011$ | 4.0

Note. — The color-separated bins in the bottom three rows are defined as Bin
1: $E(B-V)<-0.05$, Bin 2: $-0.05\leq E(B-V)<0.07$, Bin 3: $0.070\leq E(B-V)$.

The differences between the reported FIRST and VLASS numbers in columns (2)
and (6) are due to some objects not having coverage in the respective radio
imaging survey.

###  5.2 a note on incompleteness

The red W2M QSOs obey the KX color cuts (i.e., red box in Figure 8) and are
nearly spectroscopically complete with either a spectrum from SDSS or spectra
that we obtained in the near-IR and optical. The blue QSO sample is not
spectroscopically complete. Among the sources obeying the WISE color selection
outside the KX color cuts, 747 sources lacked a spectrum in SDSS. These
sources could be a mix of QSOs, galaxies, and perhaps a few stars. Might the
exclusion of the QSOs missed from this set of objects affect the different
mean radio properties that we observe?

Given that the radio properties of blue QSOs do not change with $W1-W2$ color,
and that an object is more likely to be a galaxy below the diagonal line cut
(Eqn. 2), we consider the 262 objects with no SDSS spectrum above that line.
We found that 69 of these objects had a FIRST match within 2″, which is 26%,
consistent with the 28% found for confirmed blue QSOs. Even under the most
extreme assumption that all 69 sources are QSOs, while the remaining objects
are not, the FIRST-detected fraction would rise to 31%, which is still far
below the 52% found for red QSOs. A visual examination of the SDSS images of
these spectrum-free sources reveals that a majority appear to be extended
galaxies. We therefore conclude that the spectroscopic incompleteness of the
blue QSO sample does not bias our conclusions about the differences between
the radio properties of red and blue QSOs.

##  6 Discussion

We have defined a sample of QSOs selected in a radio-independent manner such
that difference in their radio properties ought to reflect intrinsic
differences between the two populations. We took care to correct for host
galaxy light that reddens the spectral shape of otherwise unobscured QSOs to
define a clear distinction between dust-reddened QSOs and unreddened, blue
QSOs with an intrinsically redder optical continuum. We then identified a
luminosity-restricted subsample to further minimize the effects of reddening
combined with the survey’s flux limit. Below, we discuss the implications of
the differences seen in the radio properties of these red and blue QSOs. We
also estimate the fraction of red QSOs in this radio-independent population,
given the differences in the radio properties of red and blue QSOs uncovered
in this work.

###  6.1 Radio emission in red quasars

Although most of the QSOs in our sample are undetected in FIRST and VLASS, we
investigate the radio emission of the sources with detections to understand
the radio luminosities and radio-loudnesses of the blue and red populations.
Because radio-loudness is defined as the ratio of radio to optical emission,
the presence of reddening and extinction at optical wavelengths will result in
red QSOs being artificially apparent as radio-loud. This was addressed in
Glikman et al. (2012) by de-reddening the optical flux based on $E(B-V)$
(e.g., §3). However, Klindt et al. (2019) and Fawcett et al. (2020) define
radio-loudness using the QSO luminosity at 6 $\mu$m,

$\mathcal{R}=\log_{10}\frac{1.4\times 10^{16}L_{1.4{\rm GHz}}}{L_{6\mu{\rm
m}}},$ (5)

which is less sensitive to dust extinction but still probes the QSO continuum.
In this formulation, the radio-quiet/radio-loud divide occurs at
$\mathcal{R}=-4.6$. Figure 23 shows the 6 $\mu$m luminosity versus the 1.4 GHz
radio luminosity for the QSOs with FIRST detections that are part of the
luminosity-restricted subsample and at $z>0.1$. The solid black line
represents $\mathcal{R}=-4.6$ and demonstrates that all but four of the red
QSOs are radio quiet. A similar investigation is performed in Fawcett et al.
(2020) who find a larger representation of radio-loud red QSOs.

Figure 23: Rest-frame infrared luminosity at 6 $\mu$m vs. radio luminosity at
1.4 GHz for the luminosity-restricted blue and red QSOs at $z>0.1$ with FIRST
detections. The diagonal black line represents $\mathcal{R}=-4.6$, which is
the curtoff for radio-loud QSOs, seen to the right. The dashed black line is
$\mathcal{R}=-3.5$ which includes radio-intermediate sources.

The median radio properties of red and blue QSOs are significantly distinct,
with red QSOs emitting $\sim 2-3$ times more flux, at both FIRST (1.4 GHz) and
VLASS frequencies (3.0 GHz). Assuming a power-law shape to the radio spectrum,
$S_{\nu}\propto\nu^{\alpha}$, we can compare the ratio of flux densities for
the stacked red and blue QSOs to investigate differences in their spectral
slopes. Using the median flux values in FIRST and VLASS, reported in Table 6,
we compute an estimated spectral index,

$\alpha=\frac{\log{(S_{\rm FIRST}/S_{\rm
VLASS})}}{\log{(1.4{\rm~{}GHz}/3.0{\rm~{}GHz})}},$ (6)

with an uncertainty, $\alpha_{\rm err}$, derived using standard error
propagation. We find that the red quasars have a median index of $\alpha_{\rm
red}=-0.70\pm 0.05$ while the blue quasars have significantly flatter spectra,
with $\alpha_{\rm blue}=-0.34\pm 0.04$.

It is interesting to note that the median spectral index for the red QSOs is
the spectral index typically used for applying $k$-corrections to radio
luminosities (e.g. Kimball et al., 2011). However, given the lack of robust
calibration of a ‘snapshot bias’ correction for VLASS, these values are best
used in comparison between the two populations rather than as a true value.
However, the ‘snapshot bias’ correction involves scaling the flux upward,
which in both the red and blue cases would result in a more negative (steeper)
slope.

Indeed, a similar trend was seen in the F2M quasars with contemporaneous VLA
observations of 44 F2M sources at 1.4 GHz and 8.3 GHz (i.e., 20 cm and 3.6 cm;
Glikman et al., 2007). These were compared to 214 FIRST Bright Quasar Survey
(FBQS; Gregg et al., 1996) sources, which is the blue quasar equivalent of
F2M, observed with the same VLA configuration (Lacy et al., 2001). When broken
up into bright ($S_{\rm pk,FIRST}>10$ mJy) and faint ($S_{\rm pk,FIRST}<10$
mJy) subsamples, the F2M red quasars were found to have a similar spread of
spectral indices in the bright sample, but the faint red quasars had steeper
spectral slopes with the a median spectral index of $\alpha=-0.68$ and most
sources having $\alpha$ between $-0.5$ and $-1.1$121212Among the 11 F2M red
quasars recovered here, only two had spectral indices measured in Glikman et
al. (2007). F2M1004+1229 has $\alpha=0.03$ and, with a FIRST flux density of
12.3 mJy, belongs in the bright sample. F2M2216$-$0054 has $\alpha=-1.08$ and,
with a FIRST flux density of 1.3 mJy, is faint. With only two sources that
span the range of spectral indices, we are unable to generalize more broadly.
. Given that only 7 W2M red QSOs have $S_{\rm pk,FIRST}>10$ mJy and that the
spectral index estimated from the stacked images is similar to the median
spectral index for the F2M red quasars with $S_{\rm pk,FIRST}<10$ mJy, we are
likely probing a similar population of faint but enhanced radio sources that
are somehow associated with a dusty environment that is reddening their
optical to near-infrared spectra.

This apparent connection between enhanced radio emission and reddening
suggests a different physical mechanism driving the radio emission of the red
and blue populations. It may be tempting to explain the flatter median
spectral index for the blue QSOs as due to orientation effects from
relativistic jets. To check for the influence of core-dominated emission from
radio-loud QSOs, whose beamed emission may flatten the slope, we stacked only
QSOs with $\mathcal{R}\leq-3.5$. These objects appear to the left of the
dashed line in Figure 23 and excludes the radio-loudest systems (18 blue and 1
red). We found $\alpha_{\rm red}=-0.69\pm 0.05$ and $\alpha_{\rm
blue}=-0.36\pm 0.04$, which are unchanged from the full sample, within the
uncertainties. Furthermore, given that the vast majority of the QSOs are blue,
it is unlikely that they are all viewed within the small angle needed to
witness relativistic beaming effects.

Given that our sample is dominated by relatively low redshift sources, we
consider the possibility that star formation is a significant contributor to
the radio emission. Kimball et al. (2011) constructed a radio-luminosity
function for luminous blue QSOs, including radio-quiet sources, which has a
shape characterized by two-components; star formation dominates in sources
with $\log(L_{6\rm GHz}{\rm[W~{}Hz^{-1}]})\lesssim 22.5$. In Figure 23, we see
that among the FIRST-detected sources, there are very few QSOs near this
luminosity threshold at 1.4 GHz. Furthermore, the red QSOs are not found below
$\log(L_{1.4\rm GHz}{\rm[W~{}Hz^{-1}]})\sim 23$. Since these detected red QSOs
account for $\gtrsim 50\%$ of the red QSO sample, we cannot attribute most of
the radio emission to star-formation processes. Other studies of radio
emission from radio quiet QSOs have also argued against star-formation as the
dominant source of emission in radio-quiet QSOs (e.g., White et al., 2017;
Laor et al., 2019) and even in red QSOs, Fawcett et al. (2020) argue that the
enhanced radio emission is likely due to AGN activity. We note that these
conclusions are based on indirect analyses and more targeted studies, such as
high spatial resolution, multi-frequency radio imaging of the red QSOs would
directly test whether the radio emission is AGN dominated.

Recently, employing high-resolution radio imaging (0$\farcs$2 at $\sim 1.5$
GHz) of a sample of red and blue QSOs chosen from the Klindt et al. (2019)
study, Rosario et al. (2021) found unresolved radio cores in the majority of
both groups arising from regions smaller than 2 kpc in size. However, the red
QSOs did have a significantly higher fraction of extended or multi-component
radio emission compared with the blue QSOs. The authors propose that dusty
winds are both reddening the QSOs and driving shocks that generate radio
emission. This interpretation is corroborated by Calistro Rivera et al. (2021)
who find an excess of near-infrared emission in the mean SEDs of red QSOs also
drawn from the Klindt et al. (2019) study. The authors interpret this excess
emission as arising from hot outflowing dust, i.e., a dusty wind (see also
Zakamska & Greene, 2014).

We note that if we consider just the excess radio emission in the red QSO
population, we find a spectral index of $\alpha\simeq-0.9$, which is even
steeper. Laor et al. (2019) find a strong correlation between $\alpha$ and
Eddington ratio ($L/L_{\rm Edd}$) for a sample of luminous radio-quiet QSOs
such that the steeper the slope, the higher the accretion rate. The F2M red
quasars are known to have Eddington ratios that are significantly higher than
comparable blue quasars (e.g., Urrutia et al., 2012; Kim et al., 2015) and the
same may be true for the W2M red QSOs. Laor et al. (2019) interpret the steep
spectral slope in radio-quiet high-accretion QSOs as possibly due to AGN-
driven winds generating outflows and associated shocked gas that results in
synchrotron radiation emitted in the radio. This is consistent with the
presence of broad absorption line systems seen in the F2M red quasars (Urrutia
et al., 2009; Glikman et al., 2012). The highest redshift W2M QSOs also show
evidence for outflows, either in absorption or emission, when an optical
(i.e., rest-frame UV) spectrum exists. These interpretations are also
consistent with the dusty wind scenario proposed by Klindt et al. (2019),
Rosario et al. (2021), and Zakamska & Greene (2014).

###  6.2 The fraction of red QSOs

If the difference between blue and red QSOs is not due to orientation with
respect to our viewing angle, as is suggested by the radio results, then we
can directly compare the two populations to find a true fraction of each sub
group within the Type 1 QSO population (i.e., we can control for viewing
angle, as illustrated in Figure 3, while exploring dust-reddening as the
variable). With this assumption, we can compute the fraction of red quasars to
be,

$\frac{N_{red~{}QSO}}{N_{red~{}QSO}+N_{blue~{}QSO}}.$ (7)

However, as Figure 17 shows, the fraction of red QSOs appears to have a strong
luminosity dependence which must be considered when determining a red QSO
fraction. In both panels, red QSOs make up a large percentage of all QSOs at
high luminosities. On the lower luminosity end blue QSOs dominate.

Figure 24: Red QSO fraction as a function of bolometric luminosity, determined
from the histograms in Figure 17 and obeying Equation 7. The data were
smoothed to 1 dex bins and uncertainty intervals (shaded gray areas) were
computed using binomial proportion confidence interval for small number counts
(Wilson, 1927). The black line and dark shaded area represent the full W2M
sample, while the gray line and light gray shading represent the FIRST-
detected subsample. At the highest luminosities, we see that red QSOs dominate
the overall QSO population, regardless of their radio properties. At
decreasing luminosities, the radio-detected sources have a higher fraction
compared the full sample of QSOs. These results are consistent with previous
findings that red QSOs make up $\sim 20\%$ (and up to $40\%$) of all QSOs, at
least at the highest luminosities.

Figure 24 shows the fraction of red QSOs found in our survey as a function of
bolometric luminosity. We determine this fraction by using the red QSO
histogram and the luminosity-restricted blue QSO histogram, with a binning of
1 dex to smooth out fluctuations, and taking a ratio following Equation 7.
There is a strong luminosity effect, showing that red QSOs dominate the
overall QSO populations regardless of radio properties. The radio detected red
QSO fraction is higher at decreasing luminosities. However, in both cases, red
QSOs make up at least 20% and up to 40% of the overall QSO population at the
highest luminosities. These fractions are lower limits, since we miss more
heavily reddened QSOs, while accounting for the blue QSOs, at the same
luminosities.

We note that these fractions are consistent with the fraction of red quasars
estimated in Glikman et al. (2012) for the F2M sample, which extends to
fainter magnitudes ($K\lesssim 15.5$ mag, the 2MASS limit). Figure 17 of that
paper shows that red quasars dominate at the highest de-reddened absolute
$K$-band magnitudes of a similar range. LaMassa et al. (2017) see similar
behavior for X-ray-selected red QSOs, finding that red QSOs make up a larger
fraction of all X-ray-selected QSOs when corrected for absorption; above
$L_{X}=10^{44}$ erg s-1, red QSOs make up $\sim 20\%$ of all quasars in that
luminosity regime. In Glikman et al. (2018a), the red QSO sample identified
over Stripe 82 is combined with deeper, mid-infrared red QSO surveys conducted
over smaller areas to enable a luminosity function calculation and the
comparison of red QSO space density versus blue and Type 2 QSOs (which are
more heavily obscured, likely due to orientation). That study finds that red
QSOs make up $\sim 30-40\%$ of the overall QSO population at the highest
luminosities, above $\nu L_{5\mu{\rm m}}=10^{45.5}$ erg s-1, which this work
corroborates.

However, the limitations of the relatively shallow $K<14.7$ mag limit of the
W2M survey limits our ability to disentangle the effects of redshift,
luminosity, reddening, and radio emission and provides motivation for a
deeper, mid-infrared-selected QSO study. In fact, Glikman et al. (2013) showed
that the space density of FIRST-selected quasars rises more steeply than for
blue quasars when approaching a deeper flux limit of $K=17$ mag. Work is
currently underway to expand the W2M QSO sample to a fainter $K$-band
magnitude limit.

###  6.3 Mergers, radio emission, and red quasars as an evolutionary phase

The distinct differences between the radio properties of red and blue QSOs
complicates our ability to address the role of mergers in the co-evolutionary
picture for QSOs and their host galaxies. It is unclear whether the dusty
winds proposed to explain the enhanced radio and reddening properties of red
QSOs are associated with the high merger fraction ($>80\%$) seen in the host
galaxies of F2M red quasars (Urrutia et al., 2008; Glikman et al., 2015),
which are radio selected. Many of the F2M red quasar properties are consistent
with a ‘blowout’ phase. The existence of broad lines rules out the source of
reddening being due to orientation along a line-of-sight that intersects with
a dusty torus. Their spectra show an unusually high fraction of low-ionization
broad absorption line sources (LoBALs; Urrutia et al., 2009; Glikman et al.,
2012), high accretion rates compared to unreddened quasars (Urrutia et al.,
2012; Kim et al., 2015), and absorption-corrected bolometric luminosities that
are higher than blue quasars at similar redshifts (Glikman et al., 2012;
Treister et al., 2012).

However, the results for red and obscured QSOs at lower radio luminosities are
a less clear-cut. Zakamska et al. (2019) analyzed HST images of ERQs, which
are hyper-luminous broad-lined QSOs at $z=2-3$. These objects are infrared
bright, heavily obscured in X-rays, exhibit outflows in their emission line
profiles, and are possibly accreting at super-Eddington rates. The radio
properties of ERQs are interpreted as also arising from accretion-driven
winds. However, these objects, which are similar to the F2M red quasars in
many ways, do not show a significant merger fraction ($\sim 20\%$). Their main
distinction is in their radio properties, with the F2M red quasars being an
order of magnitude more luminous at rest-frame 1.4 GHz ($10^{41.9}$ erg s-1
versus $10^{40.9}$ erg s-1; Zakamska et al., 2019).

Interestingly, Chiaberge et al. (2015) studied a large sample of radio
galaxies from the 3C catalog between $1<z<2.5$ with HST and found a remarkably
high merger rate ($>$90%), while radio-quiet Type 2 analogs at the same
redshifts have a merger fraction consistent with inactive galaxies. They
conclude that major mergers not only trigger star formation and SMBH growth as
the models predict, but are also responsible for launching jets. It is
possible, then, that the enhanced radio emission in the W2M red QSOs can be
due to a combination of phenomena, where winds are a ubiquitous main driver of
low-level radio emission and jets, possibly associated with mergers,
contribute to the radio emission in the more radio-luminous sources.

The main driver for dust-driven winds may then be bolometric luminosity or
accretion rate. Dust-driven winds are invoked to explain the properties of the
most luminous Type 2 QSOs at $z\lesssim 1$ (Zakamska & Greene, 2014). Gas-rich
mergers are particularly effective at providing an abundant fuel source to the
SMBH, resulting in a higher luminosity AGN possibly explaining them being a
significant presence among the radio-selected red quasars. Given that the
‘blowout’ phase derives its meaning in a merger scenario, it is important to
distinguish between different feedback phenomena before interpreting the red
QSO fraction (e.g., Figure 24) as a phase duration. One way to test this would
be to study the host morphologies of W2M QSOs with HST to determine whether
mergers are a universal phenomenon among all red QSOs or whether radio-
selection is biased toward merging systems. This analysis is currently
underway with an HST imaging program.

##  7 Conclusions

We have identified a sample of 1,112 QSOs selected according to their mid-
infrared colors with a near-infrared flux limit of $K<14.7$ mag over 2,213
deg2. This selection method identifies blue and red QSOs with minimal
contamination and reddening bias and without the need for a radio selection
criterion. We performed a careful analysis of the QSO spectra, removing host
galaxy light which artificially mimics reddening by dust in some otherwise
blue QSOs. We also defined a luminosity-restricted subsample in the $K$-band,
after correcting for reddening. This enabled us to create intrinsically blue
(798) and red (63) QSO subsamples whose properties we studied and compared.

We investigated the fraction of sources detected in two radio surveys, FIRST
(1.4 GHz) and VLASS ($2-4$ GHz), and we employed radio stacking to study the
flux densities of sources undetected in these surveys. We found that red QSOs
are significantly more likely to be detected at both 1.4 GHz and 3.0 GHz and
are more likely to appear compact in morphology. We also found that red QSOs
have brighter median radio flux densities compared with blue QSOs. These
results are consistent with recent work by Klindt et al. (2019) and Fawcett et
al. (2020) who find similar radio enhancement at 1.4 GHz for SDSS quasars that
have redder colors compared with their blue counterparts. We note that,
compared with red QSOs in the SDSS sample, the W2M red QSOs reach higher
$E(B-V)$ values and show a more pronounced distinction in the stacked radio
flux ratios. Considering both frequencies, we find that red QSOs have steeper
median radio spectra compared with blue QSOs (i.e., red QSOs have higher FIRST
to VLASS flux ratios than blue QSOs). We speculate that a dusty AGN-driven
wind can account for both the unique radio and reddening properties of red
QSOs, as has been noted elsewhere (Zakamska & Greene, 2014; Zakamska et al.,
2019; Klindt et al., 2019; Fawcett et al., 2020; Calistro Rivera et al., 2021;
Rosario et al., 2021).

The red QSOs in this study are among the more luminous QSOs, especially at
high redshift ($z>1.5$), though our survey is not sensitive to the same
luminosity distributions sampled for the red and blue QSOs. We also note an
absence of red QSOs at $z<0.1$, which is consistent with evolutionary behavior
seen in previous work (e.g., Glikman et al., 2018a; LaMassa et al., 2017). We
therefore investigated the fraction of red and blue QSOs as a function of de-
reddened absolute $K$-band magnitude in a de-reddened-luminosity-restricted
and redshift-matched subsample. We find that red QSOs dominate the QSO
population at the highest luminosities, remaining a significant fraction of
the QSO population at $\log(L_{\rm bol}[{\rm erg~{}s^{-1}}])>46$ with the
radio-detected red QSOs having a $\sim 40\%$ higher fraction.

The results of this study suggest that previous conclusions about the fraction
of red quasars, determined from radio-selected samples, is too high to be
extended to the overall QSO population, and a radio-independent selection is
essential for understanding the nature of dust-reddened QSOs. The fact that
red QSOs appear to be a predominantly high luminosity phenomenon with distinct
radio properties showing enhanced emission, strongly implies that red QSOs are
not simply an apparent orientation effect but are rather a distinct population
that can shed light on supermassive black-hole growth and the quasar
phenomenon.

We thank the anonymous referee for a careful reading of the text and for
suggesting revisions and analyses that greatly improved the presentation of
the paper. EG acknowledges the generous support of the Cottrell Scholar Award
through the Research Corporation for Science Advancement. SGD and MJG
acknowledge a partial support from the NASA grant 16-ADAP16-0232 and the NSF
grants AST-1749235, AST-1815034, and AST-2108402. This research made use of
the cross-match service provided by CDS, Strasbourg. Research at Lick
Observatory is partially supported by a generous gift from Google. This
publication makes use of data products from the Wide-field Infrared Survey
Explorer, which is a joint project of the University of California, Los
Angeles, and the Jet Propulsion Laboratory/California Institute of Technology,
funded by the National Aeronautics and Space Administration. This publication
makes use of data products from the Two Micron All Sky Survey, which is a
joint project of the University of Massachusetts and the Infrared Processing
and Analysis Center/California Institute of Technology, funded by the National
Aeronautics and Space Administration and the National Science Foundation.
Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the
Participating Institutions, the National Science Foundation, and the U.S.
Department of Energy Office of Science. The SDSS-III web site is
http://www.sdss3.org/. SDSS-III is managed by the Astrophysical Research
Consortium for the Participating Institutions of the SDSS-III Collaboration
including the University of Arizona, the Brazilian Participation Group,
Brookhaven National Laboratory, Carnegie Mellon University, University of
Florida, the French Participation Group, the German Participation Group,
Harvard University, the Instituto de Astrofisica de Canarias, the Michigan
State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence
Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max
Planck Institute for Extraterrestrial Physics, New Mexico State University,
New York University, Ohio State University, Pennsylvania State University,
University of Portsmouth, Princeton University, the Spanish Participation
Group, University of Tokyo, University of Utah, Vanderbilt University,
University of Virginia, University of Washington, and Yale University.
Palomar, IRTF, APO, LBT, Lick

## References

  * Abolfathi et al. (2018) Abolfathi, B. et al. 2018, ApJS, 235, 42
  * Ahn et al. (2012) Ahn, C. P. et al. 2012, ApJS, 203, 21
  * Assef et al. (2013) Assef, R. J. et al. 2013, ApJ, 772, 26
  * Assef et al. (2018) Assef, R. J., Stern, D., Noirot, G., Jun, H. D., Cutri, R. M., & Eisenhardt, P. R. M. 2018, ApJS, 234, 23
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al. 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2013) —. 2013, A&A, 558, A33
  * Becker et al. (1995) Becker, R. H., White, R. L., & Helfand, D. J. 1995, ApJ, 450, 559
  * Bennett et al. (2013) Bennett, C. L. et al. 2013, ApJS, 208, 20
  * Bolton et al. (2012) Bolton, A. S. et al. 2012, AJ, 144, 144
  * Calistro Rivera et al. (2021) Calistro Rivera, G. et al. 2021, A&A, 649, A102
  * Chiaberge et al. (2015) Chiaberge, M., Gilli, R., Lotz, J. M., & Norman, C. 2015, ApJ, 806, 147
  * Cushing et al. (2004) Cushing, M. C., Vacca, W. D., & Rayner, J. T. 2004, PASP, 116, 362
  * Di Matteo et al. (2005) Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604
  * Donley et al. (2012) Donley, J. L. et al. 2012, ApJ, 748, 142
  * Fabian (2012) Fabian, A. C. 2012, ARA&A, 50, 455
  * Farrah et al. (2012) Farrah, D. et al. 2012, ApJ, 745, 178
  * Fawcett et al. (2020) Fawcett, V. A., Alexander, D. M., Rosario, D. J., Klindt, L., Fotopoulou, S., Lusso, E., Morabito, L. K., & Calistro Rivera, G. 2020, MNRAS, 494, 4802
  * Frieman et al. (2008) Frieman, J. A. et al. 2008, AJ, 135, 338
  * Gebhardt et al. (2000) Gebhardt, K. et al. 2000, ApJ, 539, L13
  * Georgakakis et al. (2009) Georgakakis, A., Clements, D. L., Bendo, G., Rowan-Robinson, M., Nandra, K., & Brotherton, M. S. 2009, MNRAS, 394, 533
  * Glikman et al. (2004) Glikman, E., Gregg, M. D., Lacy, M., Helfand, D. J., Becker, R. H., & White, R. L. 2004, ApJ, 607, 60
  * Glikman et al. (2007) Glikman, E., Helfand, D. J., White, R. L., Becker, R. H., Gregg, M. D., & Lacy, M. 2007, ApJ, 667, 673
  * Glikman et al. (2018a) Glikman, E. et al. 2018a, ApJ, 861, 37
  * Glikman et al. (2018b) Glikman, E., Rusu, C. E., Djorgovski, S. G., Graham, M. J., Stern, D., Urrutia, T., Lacy, M., & O’Meara, J. M. 2018b, arXiv e-prints, arXiv:1807.05434
  * Glikman et al. (2015) Glikman, E., Simmons, B., Mailly, M., Schawinski, K., Urry, C. M., & Lacy, M. 2015, ApJ, 806, 218
  * Glikman et al. (2012) Glikman, E. et al. 2012, ApJ, 757, 51
  * Glikman et al. (2013) —. 2013, ApJ, 778, 127
  * Gordon & Clayton (1998) Gordon, K. D., & Clayton, G. C. 1998, ApJ, 500, 816
  * Gordon et al. (2021) Gordon, Y. A. et al. 2021, arXiv e-prints, arXiv:2102.11753
  * Gregg et al. (1996) Gregg, M. D., Becker, R. H., White, R. L., Helfand, D. J., McMahon, R. G., & Hook, I. M. 1996, AJ, 112, 407
  * Hamann et al. (2017) Hamann, F. et al. 2017, MNRAS, 464, 3431
  * Herter et al. (2008) Herter, T. L. et al. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7014, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series
  * Hewett & Wild (2010) Hewett, P. C., & Wild, V. 2010, MNRAS, 405, 2302
  * Hopkins et al. (2005) Hopkins, P. F., Hernquist, L., Martini, P., Cox, T. J., Robertson, B., Di Matteo, T., & Springel, V. 2005, ApJ, 625, L71
  * Jarrett et al. (2011) Jarrett, T. H. et al. 2011, ApJ, 735, 112
  * Kayser et al. (1990) Kayser, R., Surdej, J., Condon, J. J., Kellermann, K. I., Magain, P., Remy, M., & Smette, A. 1990, ApJ, 364, 15
  * Kim et al. (2015) Kim, D., Im, M., Glikman, E., Woo, J.-H., & Urrutia, T. 2015, ApJ, 812, 66
  * Kimball et al. (2011) Kimball, A. E., Kellermann, K. I., Condon, J. J., Ivezić, Ž., & Perley, R. A. 2011, ApJ, 739, L29
  * Klindt et al. (2019) Klindt, L., Alexander, D. M., Rosario, D. J., Lusso, E., & Fotopoulou, S. 2019, MNRAS, 488, 3109
  * Lacy et al. (2020) Lacy, M. et al. 2020, PASP, 132, 035001
  * Lacy et al. (2001) Lacy, M., Laurent-Muehleisen, S. A., Ridgway, S. E., Becker, R. H., & White, R. L. 2001, ApJ, 551, L17
  * Lacy et al. (2007) Lacy, M., Petric, A. O., Sajina, A., Canalizo, G., Storrie-Lombardi, L. J., Armus, L., Fadda, D., & Marleau, F. R. 2007, AJ, 133, 186
  * Lacy et al. (2004) Lacy, M. et al. 2004, ApJS, 154, 166
  * LaMassa et al. (2019) LaMassa, S. M., Georgakakis, A., Vivek, M., Salvato, M., Ananna, T. T., Urry, C. M., MacLeod, C., & Ross, N. 2019, ApJ, 876, 50
  * LaMassa et al. (2017) LaMassa, S. M. et al. 2017, ApJ, 847, 100
  * Laor et al. (2019) Laor, A., Baldi, R. D., & Behar, E. 2019, MNRAS, 482, 5513
  * Maddox et al. (2008) Maddox, N., Hewett, P. C., Warren, S. J., & Croom, S. M. 2008, MNRAS, 386, 1605
  * Mateos et al. (2012) Mateos, S. et al. 2012, MNRAS, 426, 3271
  * Mullaney et al. (2013) Mullaney, J. R., Alexander, D. M., Fine, S., Goulding, A. D., Harrison, C. M., & Hickox, R. C. 2013, MNRAS, 433, 622
  * Rau et al. (2016) Rau, U., Bhatnagar, S., & Owen, F. N. 2016, AJ, 152, 124
  * Rayner et al. (2003) Rayner, J. T., Toomey, D. W., Onaka, P. M., Denault, A. J., Stahlberger, W. E., Vacca, W. D., Cushing, M. C., & Wang, S. 2003, PASP, 115, 362
  * Richards et al. (2002) Richards, G. T. et al. 2002, AJ, 123, 2945
  * Richards et al. (2003) —. 2003, AJ, 126, 1131
  * Richards et al. (2006) —. 2006, ApJS, 166, 470
  * Rosario et al. (2021) Rosario, D. J., Alexander, D. M., Moldon, J., Klindt, L., Thomson, A. P., Morabito, L., Fawcett, V. A., & Harrison, C. M. 2021, MNRAS, 505, 5283
  * Ross et al. (2015) Ross, N. P. et al. 2015, MNRAS, 453, 3932
  * Ross et al. (2012) —. 2012, ApJS, 199, 3
  * Sanders et al. (1988) Sanders, D. B., Soifer, B. T., Elias, J. H., Madore, B. F., Matthews, K., Neugebauer, G., & Scoville, N. Z. 1988, ApJ, 325, 74
  * Sarzi et al. (2006) Sarzi, M. et al. 2006, MNRAS, 366, 1151
  * Schneider et al. (2003) Schneider, D. P. et al. 2003, AJ, 126, 2579
  * Skrutskie et al. (2006) Skrutskie, M. F. et al. 2006, AJ, 131, 1163
  * Smolčić et al. (2017) Smolčić, V. et al. 2017, A&A, 602, A1
  * Stern et al. (2012) Stern, D. et al. 2012, ApJ, 753, 30
  * Stern et al. (2005) —. 2005, ApJ, 631, 163
  * 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
  * Tody (1986) Tody, D. 1986, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 627, Instrumentation in astronomy VI, ed. D. L. Crawford, 733
  * Tody (1993) Tody, D. 1993, in Astronomical Society of the Pacific Conference Series, Vol. 52, Astronomical Data Analysis Software and Systems II, ed. R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, 173
  * Treister et al. (2012) Treister, E., Schawinski, K., Urry, C. M., & Simmons, B. D. 2012, ApJ, 758, L39
  * Urrutia et al. (2009) Urrutia, T., Becker, R. H., White, R. L., Glikman, E., Lacy, M., Hodge, J., & Gregg, M. D. 2009, ApJ, 698, 1095
  * Urrutia et al. (2008) Urrutia, T., Lacy, M., & Becker, R. H. 2008, ApJ, 674, 80
  * Urrutia et al. (2012) Urrutia, T., Lacy, M., Spoon, H., Glikman, E., Petric, A., & Schulz, B. 2012, ApJ, 757, 125
  * Urry & Padovani (1995) Urry, C. M., & Padovani, P. 1995, PASP, 107, 803
  * Vacca et al. (2003) Vacca, W. D., Cushing, M. C., & Rayner, J. T. 2003, PASP, 115, 389
  * Warren et al. (2000) Warren, S. J., Hewett, P. C., & Foltz, C. B. 2000, MNRAS, 312, 827
  * Webster et al. (1995) Webster, R. L., Francis, P. J., Peterson, B. A., Drinkwater, M. J., & Masci, F. J. 1995, Nature, 375, 469
  * White et al. (2007) White, R. L., Helfand, D. J., Becker, R. H., Glikman, E., & de Vries, W. 2007, ApJ, 654, 99
  * White et al. (2003) White, R. L., Helfand, D. J., Becker, R. H., Gregg, M. D., Postman, M., Lauer, T. R., & Oegerle, W. 2003, AJ, 126, 706
  * White et al. (2017) White, S. V., Jarvis, M. J., Kalfountzou, E., Hardcastle, M. J., Verma, A., Cao Orjales, J. M., & Stevens, J. 2017, MNRAS, 468, 217
  * Wilson (1927) Wilson, E. B. 1927, Journal of the American Statistical Association, 22, 209
  * Wilson et al. (2004) Wilson, J. C. et al. 2004, in Proc. SPIE, Vol. 5492, Ground-based Instrumentation for Astronomy, ed. A. F. M. Moorwood & M. Iye, 1295–1305
  * Wright et al. (2010) Wright, E. L. et al. 2010, AJ, 140, 1868
  * Zafar et al. (2015) Zafar, T. et al. 2015, A&A, 584, A100
  * Zakamska & Greene (2014) Zakamska, N. L., & Greene, J. E. 2014, MNRAS, 442, 784
  * Zakamska et al. (2019) Zakamska, N. L. et al. 2019, MNRAS, 489, 497