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

We present a highly complete sample of broad-line (Type 1) QSOs out to z ~ 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 ~ 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.


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 −σ 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 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 ∼ 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 130 dustreddened quasars via the same method (hereafter referred to as F2M quasars Glikman et al. 2007Glikman et al. , 2012Urrutia et al. 2009) that have broad emission lines and are moderately obscured by A V ∼ 1 − 4 (0.1 < E(B − V ) 1.5) across a broad range of redshifts (0.1 < z 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 mergerdominated systems (Urrutia et al. 2008;Glikman et al. 2015), and often have broad absorption lines that are typically associated with outflows and feedback (LoB-ALs 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 ∼ 20 − 30% of the overall quasar population (Glikman et al. , 2018a, at least at the highest luminosities. When interpreted as an evolutionary phase, this fraction implies that the duration of this transitional phase is ∼ 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 ∼ 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 ∼ 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 ) 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 QSOs 1 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 µm down to flux densities of 0.08, 0.11, 1, and 6 mJy, respectively 2 , 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. 2013Assef et al. , 2018. In a pilot study, Glikman et al. (2018a) identified a complete sample of QSOs using near-to-mid-infrared color selection over ∼ 260 deg 2 that overlap the SDSS Stripe 82 legacy field (Frieman et al. 2008). Here we expand upon that work over a 7.5× 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 QSOs 3 , 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 (∼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 ΛCDM cosmology with H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, and Ω Λ = 0.7 when computing cosmology-dependent values (Bennett et al. 2013).

SAMPLE SELECTION
Our survey area covers an equatorial region overlapping the SDSS over two fields: a region spanning a range in right ascension of α = 8 h − 16 h and in declination of δ = 0.5 • − 17 • , available for follow-up in the spring months, and the region over Stripe 82 identical to that in Glikman et al. (2018a), α = 20 h 40 m − 3 h 56 m and δ = −1.25 • to +1.25 • (excluding the region 00 h < α < 00 h 15 m ). This covers a total area of 2213 deg 2 (1950 deg 2 and 263 2 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 subpopulations 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 broadline 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 ) ≥ 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.

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  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. We began by selecting sources from the WISE catalog, which presents its photometry in the Vega system. We select sources with 0.5 < W 1 − W 2 < 2 and 2 < W 2 − W 3 < 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: 4 At 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.
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 Wright et al. (2010), which shows the location of different astrophysical populations in this space. The blue circles are SDSSidentified 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).

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  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).
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. 2004Gregg et al. 1996;Schneider et al. 2003), our sample will focus on all objects with line widths broader than 1000 km s −1 .

QSOs in SDSS
Blue Type-1 Red Type-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.
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 Hα, 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α 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α line constitutes ≥ 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 F W HM > 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 interface 5 and used the value-added measurements provided in their multiextension 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 details 6 ). The line widths are reported in terms of the Gaussian σ parameter, which we convert to a FWHM velocity, v FWHM = 2.355σ. Figure 4 shows the FWHM distribution of the broadest line component in each spectrum (blue line) where v 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-widthselected 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, vF W HM , 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 vFWHM = 1000 km s −1 cutoff. 0.5 < W1 -W2 < 2 2 < W2 -W3 < 4.5 K < 14.7 SDSS QSO Spectrum 1398

Broad line 339
QSOs with SDSS spectrum: 1072 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 followup 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 W 2 − W 3 (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 W 2 − W 3 ∼ 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 W 2 − W 3 = 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: shown by the thick black line in Figures 2 and 7. 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 W 1 − W 2 > 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 W 1 and W 2 bands are more affected by reddening than the longer WISE bands, so it is not surprising that reddened QSOs have slightly redder W 1 − W 2 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).
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 rem-edy, Warren et al. (2000) showed that in optical-to-nearinfrared (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), 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 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 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  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, including their positions, optical through mid-infrared magnitudes, as well as classification and redshift based on spectroscopic followup ( §2.5).

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 . 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, 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.

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  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. 6-wide slit on UT 2013 March 14, covering Blue 3 Red 37 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.
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 provides the details on the objects' spectroscopic observations, source classification, and assigned redshift.

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 σ, in wavelength units, to a FWHM in velocity units through the expression, where c is the speed of light and λ 0 is the rest-wavelength of the line being fit. The red filled histogram in Figure 4 shows the distribution of v F W HM for the newlydiscovered 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 se-lection criteria of Eqn 2, whether the object qualifies as a red QSO ( §3) with its best-fit E(B − V ) value.

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.

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 W 1 − W 2 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 W 1 − W 2 color. We discuss the removal of host galaxy emission to account for this in Section 3.2, below.
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 unreddened 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. 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 Ks 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.        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.

Removal of host galaxy emission
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 galaxysubtracted 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.

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. (2007Glikman et al. ( , 2012Glikman et al. ( , 2013Glikman et al. ( , 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 ) ≥ 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+1641 7 , both at z = 2.52, are heavily absorbed in their optical (restframe 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 λ > 8000Å and λ > 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 nearinfrared spectral regions with sufficiently high signal-tonoise 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 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 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.
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.

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α and Paβ 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] λ5007 as Lyα, 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. 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.

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 W 4 (22µm) magnitude -a wavelength at which extinction is negligible -we can check whether W 4 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 W 4 (22µ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 W 4 and dereddened K-band absolute magnitudes versus the bolo-metric 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µm luminosity and apply a bolometric correction factor of 7.82 using the QSO spectral energy distribution (SED) from Richards et al. (2006). The W 4 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.
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µ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µ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 wellexplained 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 . 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.
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) Figure 13. Left -De-reddened absolute K-band magnitude from 2MASS versus absolute 22µm (W 4) 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 W 4. 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 W 4 (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.

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.
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 bol , and a boolean indication of whether the object is part of the luminosity restricted subsample.
We note that two of the most luminous infrared sources seen at z 2.5 are both gravitationally lensed QSOs of the rare quad variety. The blue QSO is the wellstudied 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 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 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.

RESULTS
We have constructed a carefully-selected, luminosityrestricted 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.

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 ∼ 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 ) ∼ 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)    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.

radio detected fraction
We matched the luminosity-restricted W2M sample to two radio surveys that overlap our fields. There were 249 matches within 1. 5 to the FIRST survey which reaches a 5-σ 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 µJy at 2−4 GHz with a spatial resolution of 2. 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 ∼ 7% and ∼ 17% for blue and red QSOs, compared to ∼ 28% (226/798) and ∼ 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.

radio morphologies
A major result from Klindt et al. (2019) was that when categorized by their FIRST morphology among FIRSTdetected 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 ∼ 3.5 times higher than for blue QSOs, whereas the W2M red QSOs have a compact morphology fraction that is only ∼ 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 ) ∼ 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 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.

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 subthreshold FIRST sources experience a flux density bias (dubbed "snapshot bias" in White et al. 2007) whereby the stacked source has a ∼ 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,corr = min(1.40S p , S p + 0.25 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 sidelobes, 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±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 ∼ 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 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 QSOs 8 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 im-  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. 28.6 +18.9 −13.5 a We remove one gravitational lens and two sources lacking coverage in FIRST from the total sample of 861.
ages, 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,blue = 0.390 ± 0.008 mJy and red QSOs have S p,red = 1.183 ± 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,red = 0.698±0.022 mJy, the red QSOs have a higher flux density than the blue QSOs, S p,blue = 0.301±0.006 mJy, at higher frequencies as well.
To test for biases among the blue QSOs in our sample, we conducted additional median stacks of several sub-   sets. 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 biascorrected peak flux densities in FIRST are S p,blue = 0.390 ± 0.008 mJy and S p,blue,diag = 0.402 ± 0.008 mJy, and for VLASS are S p,blue = 0.301 ± 0.006 mJy and S p,blue,diag = 0.323 ± 0.007 mJy, both of which are near identical. This provides reassurance that we can com- pare the red QSOs, which are shifted toward higher W 1 − W 2 colors, to blue QSOs that extend to lower W 1 − W 2 colors without biasing our analysis of their radio properties. 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 lowluminosity 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,blue,highz = 0.373 ± 0.007 mJy in FIRST and S p,blue,highz = 0.293 ± 0.007 mJy in VLASS which are both consistent with -and even slightly fainter thanthe full stacked sample.
Considering that the red QSOs are skewed toward higher luminosities we also stacked only blue QSOs with L 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,blue,highlum = 0.435 ± 0.009 mJy. The median VLASS flux density is S p,blue,highlum = 0.336 ± 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 dustreddened, 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 bins 9 . The bins span E(B − V ) < −0.05 for the first bin (Bin 1), −0.05 ≤ E(B − V ) < 0.07 for the second bin (Bin 2), encompassing the peak of the distribution, and 0.070 ≤ 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 biascorrected FIRST peak flux densities for the subsets are S p,Bin1 = 0.404±0.012 mJy, S p,Bin2 = 0.375±0.011 mJy, S p,Bin3 = 0.398 ± 0.012 mJy, and S p,Bin1 = 0.268 ± 0.011 mJy, S p,Bin2 = 0.316±0.011 mJy, S p,Bin3 = 0.314±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.
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 twodimensional Gaussian fits to the blue and red stacked images. We find that the red QSOs are more compact than the blue QSOs, with FWHM red = 3.1 pixels and FWHM blue ∼4.1 pixels 10 . Similarly, in VLASS, the red quasars have FWHM red = 3.4 pixels versus FWHM blue ∼ 4.0 pixels 11 . 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 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. 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.

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 W 1 − W 2 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.

DISCUSSION
We have defined a sample of QSOs selected in a radioindependent 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.

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 radioloud. 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 µm, which is less sensitive to dust extinction but still probes the QSO continuum. In this formulation, the radioquiet/radio-loud divide occurs at R = −4.6. Figure  23 shows the 6 µm luminosity versus the 1.4 GHz radio luminosity for the QSOs with FIRST detections that b An additional 20% uncertainty may be added to the VLASS stacks due to possible "snapshot bias" (see §5.1.3).
Note-The color-separated bins in the bottom three rows are defined as Bin 1: 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. . Rest-frame infrared luminosity at 6 µ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 R = −4.6, which is the curtoff for radioloud QSOs, seen to the right. The dashed black line is R = −3.5 which includes radio-intermediate sources.
are part of the luminosity-restricted subsample and at z > 0.1. The solid black line represents 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 radioloud red QSOs. The median radio properties of red and blue QSOs are significantly distinct, with red QSOs emitting ∼ 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 ν ∝ ν α , 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, with an uncertainty, α err , derived using standard error propagation. We find that the red quasars have a median index of α red = −0.70 ± 0.05 while the blue quasars have significantly flatter spectra, with α blue = −0.34 ± 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 pk,FIRST > 10 mJy) and faint (S 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 α = −0.68 and most sources having α between −0.5 and −1.1 12 . Given that only 7 W2M red QSOs have S 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 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 R ≤ −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 α red = −0.69 ± 0.05 and α blue = −0.36 ± 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 radioluminosity 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 6GHz [W Hz −1 ]) 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.4GHz [W Hz −1 ]) ∼ 23. Since these detected red QSOs account for 50% of the red QSO sample, we cannot attribute most of the radio emission to starformation processes. Other studies of radio emission from radio quiet QSOs have also argued against starformation as the dominant source of emission in radioquiet QSOs (e.g., White et al. 2017;Laor et al. 2019) and 12 Among the 11 F2M red quasars recovered here, only two had spectral indices measured in Glikman et al. (2007). F2M1004+1229 has α = 0.03 and, with a FIRST flux density of 12.3 mJy, belongs in the bright sample. F2M2216−0054 has α = −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. 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. 2 at ∼ 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 α −0.9, which is even steeper. Laor et al. (2019) find a strong correlation between α and Eddington ratio (L/L 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).

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  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 With this assumption, we can compute the fraction of red quasars to be, 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 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 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 ∼ 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 ∼ 30−40% of the overall QSO population at the highest luminosities, above νL 5µ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 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 (∼ 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 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.

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 deg 2 . 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 luminosityrestricted 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 popula-tion at log(L bol [erg s −1 ]) > 46 with the radio-detected red QSOs having a ∼ 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 dustreddened 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.