© ESO 2021

1. Introduction

Since the discovery of quasars (QSOs, Schmidt 1963), the most energetic kind of active galactic nuclei (AGN), an in-depth census of QSOs and characterisation of their properties has been carried out. These have revealed millions of QSOs (e.g. Flesch 2019) with a diversity of black hole (BH) masses (10⁷ − 10¹⁰ M_⊙) and luminosities ranging up to ∼10⁴⁸ ergs s⁻¹, pervading the Universe to its earliest epochs (z = 7.5, Bañados et al. 2018). Large optical spectroscopic surveys such as the Sloan Digital Sky Survey (SDSS; e.g. Pâris et al. 2018) have played a key role in the identification of QSOs, delivering statistical samples that allow us to achieve a more significant understanding of QSO properties. Optical spectroscopy is an efficient tracer of QSOs, since their primary emission, the accretion disk, is most luminous in the ultraviolet (UV) and blue optical regime of the spectrum. The accretion disk emission in QSOs arises from the accretion of the surrounding material onto the BH gravitational potential, heating gas clouds in the immediate surroundings to produce the broad and narrow emission lines characteristic of QSO spectra. Unobscured QSOs, which constitute the majority of the optically-selected QSO population, are thus characterised by a luminous blue thermal continuum. However, a small fraction of QSOs exhibit a diversity of redder optical and near-infrared colours. Although the origin of these red colours have been the subject of a long list of studies (e.g. Whiting et al. 2001; Wilkes et al. 2002; Urrutia et al. 2008; Glikman et al. 2007; Rose et al. 2013; Kim & Im 2018; Klindt et al. 2019), the nature of red QSOs remains an intriguing question in AGN astronomy.

Observationally, several scenarios have been investigated to explain the colours in red QSOs. These include red QSOs having an intrinsically red continuum due to particular accretion mechanisms such as higher Eddington ratios (e.g. Richards et al. 2003; Kim & Im 2018), or the red colours being just an observational effect due to the luminous contamination from other physical components of the AGN (e.g. synchrotron emission, Serjeant 1996; Whiting et al. 2001) or the stellar emission from the host galaxy. The most widely accepted explanation is, however, that the red colours are due to attenuation by dust along the line of sight (Webster et al. 1995; Glikman et al. 2007; Klindt et al. 2019). Despite extensive investigations on this topic, the origin of this attenuation and whether the dusty medium is located at nuclear or galaxy scales remains uncertain (e.g. Hickox & Alexander 2018).

Generally, obscuration in AGN is mostly linked to the geometrically thick, most likely clumpy, structure, known as the dusty ‘torus’ (e.g. Silva et al. 2004; Fritz et al. 2006; Nenkova et al. 2008; Hatziminaoglou et al. 2009; Mor et al. 2009; Alonso-Herrero et al. 2011; Stalevski et al. 2016). This structure is believed to explain the bulk of the mid- and near-infrared emission, as well as the angle-dependent obscuration of the accretion disk emission and broad line region. Indeed, the latter leads to the well-established unification model that explains the emission line variations of Type 1 (broad and narrow lines) and Type 2 QSOs (narrow lines only) as an effect of the viewing angle (Antonucci 1993; Urry & Padovani 1995; Ramos Almeida & Ricci 2017; Hickox & Alexander 2018), which constitutes one of the pillars of AGN astronomy. Within this framework, a straightforward explanation for the optical colour of red Type 1 QSOs is that it is caused solely by orientation, while red QSOs are intrinsically the same object as their unobscured counterparts. In this scenario, viewing angles of red QSOs would constitute an intermediate inclination angle between those from unobscured Type 1 QSOs and Type 2 QSOs.

A list of observational evidence, however, suggests that red QSOs present other fundamental differences which cannot be explained by orientation. Although with a lack of consensus, differences in their host galaxies have been suggested, with host galaxies of red QSOs being more likely to be undergoing mergers (e.g. Urrutia et al. 2008; Glikman et al. 2015; cf. Zakamska et al. 2019) or have higher star-formation rates than unobscured QSOs (Georgakakis et al. 2009; cf. Wethers et al. 2020). Banerji et al. (2015) reported that the red QSO luminosity functions are flatter than those of unobscured ‘typical’ QSOs, and a higher incidence of low-ionization broad absorption lines (LoBALs) has been found for red QSOs (Urrutia et al. 2009), indicating the presence of strong outflows. This set of evidence, although sometimes associated with large uncertainties based on small sample sizes and heterogeneous selections, suggests that red QSOs are a special population, potentially linked to a key phase of galaxy and black-hole evolution.

Indeed, QSOs are expected to play a key role in galaxy evolution as they are invoked as regulators of the growth of their host galaxies through a process known as AGN feedback (e.g. Harrison 2017, and references therein). Galaxy evolution models (e.g. Sanders et al. 1988; Hopkins et al. 2008; Narayanan et al. 2010; Alexander & Hickox 2012) suggest unobscured Type 1 QSOs constitute an evolutionary stage posterior to dusty star-forming galaxies (such as sub-millimetre galaxies). Red Type 1 QSOs are often invoked as potential transitioning populations in between these stages. Although a number of previous observational studies have supported this scenario, statistical samples are required to confirm that these findings are representative (e.g. Hickox & Alexander 2018).

Fundamental differences in the radio properties of red QSOs compared to the overall QSO population have recently been discovered (Klindt et al. 2019; Rosario et al. 2020; Fawcett et al. 2020). By matching the comparison samples in luminosity and redshifts, the incidence of radio detection for red QSOs have been found to be up to three times higher than for control QSOs. Using wide-field radio surveys at high radio frequencies (Klindt et al. 2019; Fawcett et al. 2020) and low radio frequencies (Rosario et al. 2020), at low and high-resolution, these studies have found that the enhanced radio detection fractions in red QSOs are connected to mostly ‘radio-quiet’ emission (as defined by the radio-to-optical luminosity ratios) and to preferentially compact morphologies. These observations cannot be explained by the orientation model which would predict the opposite result: an enhancement in compact emission from face-on blue QSOs (due to relativistic beaming) in comparison to red QSOs. These observations are, consequently, in good support of the evolutionary model whereby red QSOs are an early phase in the evolution of QSOs.

Motivated by the key findings achieved by this comparative approach, in this work we aim to investigate the origin of QSO reddening by studying their photometric and spectroscopic properties across the electromagnetic spectrum. Indeed, QSOs are prominent examples of multiwavelength phenomena, with characteristic spectral energy distributions (SEDs) from the radio to the gamma rays, which encompass a wealth of information on the physical processes that constitute their nature. While previous works have studied the SEDs of red QSOs in the optical-near-infrared regime (e.g. Glikman et al. 2012; LaMassa et al. 2017), or from the FIR to the UV or X-rays (e.g. Georgakakis et al. 2009; Urrutia et al. 2012; Mehdipour & Costantini 2018), these have focused their studies only on a handful red QSOs (1−13 sources), and without control samples to compare these with. Building on the robust selection and comparison strategy developed in Klindt et al. (2019), we use statistical samples of red and control QSOs to characterise their spectral energy distributions. In particular, we apply the Bayesian MCMC-based code for SED-modelling (AGNFITTER; Calistro Rivera et al. 2016), which combines multiwavelength-models of host galaxy and AGN physical components specifically tailored for SED studies of QSOs and AGN. A clear advantage of AGNFITTER over other SED-fitting codes is the inclusion of an independent accretion disk emission component or the ‘big blue bump’, parametrised by a free reddening parameter E(B − V)_BBB which is key for the study of red QSOs (e.g. Calistro Rivera et al. 2016; LaMassa et al. 2017).

Combining the modelled broad-band photometry with spectroscopic properties inferred from SDSS spectra (Rakshit et al. 2020), we aim to characterise the nature of red QSOs from a multiwavelength perspective and understand the origin of their red colours. In Sect. 2 we present our selection strategy and the multiwavelength broad-band photometry and spectroscopic catalogues used in this work. Additionally, we describe the settings used within the SED-fitting code AGNFITTER and other statistical methods used for this analysis. In Sect. 3.1 we present the SED-fitting results, and combine our measurements with spectroscopic data to infer intrinsic accretion properties of QSOs in Sect. 3.2. In Sect. 3.3 we present our results on physical parameters related to the obscuration in QSOs and in Sect. 3.4 we compare the host galaxy properties inferred from the SED-modelling for red and control QSOs. Finally, in Sect. 4 we discuss our results and explore the connection of our findings with the incidence of AGN dusty winds in red and control QSOs. Throughout this work, we adopt a cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3 and Ω_Λ = 0.7.

2. Data sets and methods

2.1. Main sample selection

The starting sample of our investigation are optically-selected sources defined as QSOs by the SDSS DR14 Quasar catalogue (Pâris et al. 2018). The SDSS DR14 Quasar catalogue is the largest compilation of QSOs with detailed emission line measurements and black-hole mass constraints (Rakshit et al. 2020) to date. It consists of 526 356 spectroscopically confirmed QSOs (emission line full widths at half maximum FWHM > 500 km s⁻¹) out to z = 7, comprising different sets of all archival SDSS spectroscopic campaigns. These include QSOs spectra selected from SDSS I/II Legacy surveys, targeted QSOs from the Baryon Oscillation Spectroscopic Survey (BOSS) and the extended Baryon Oscillation Spectroscopic Survey (eBOSS). The selection procedure for our main sample is shown graphically in Fig. 1 and follows the approach presented by Klindt et al. (2019), Rosario et al. (2020), Fawcett et al. (2020).

Fig. 1. Diagram to summarize our selection process. The shaded area corresponds to the general definition of red and blue QSO samples from the overall SDSS DR14 (Pâris et al. 2018), following the strategy previously presented by Fawcett et al. (2020), Rosario et al. (2020). The purple shaded area correspond to the field selection, determined by the availability of deep multiwavelength catalogues provided by the LoTSS data release in three northern fields and the public compilation of the HELP collaboration in the three equatorial fields of the GAMA survey.

One key property of the SDSS QSO catalogue is the availability of reliable spectroscopic redshift information for all QSOs in the parent sample. To provide a robust colour selection in the following steps, the first cut restricts our sample to QSOs at redshifts z < 2.5, to avoid the Lyman-break feature artificially reddening the optical spectra short-ward of λ ∼ 3200 Å. An unbiased multiwavelength comparison of QSOs requires properly accounting for the main fundamental properties that may dominate over other parameters, such as the redshift and the intrinsic luminosities of the QSOs. A challenge lies, however, in the calculation of the intrinsic bolometric quasar luminosity, which can be very different from the observed luminosity due to dust attenuation. Following Klindt et al. (2019), we adopt the rest-frame luminosity at 6 μm (L_6 μm), as the best proxy for the intrinsic luminosity. The motivation for this wavelength choice relies on the rest-frame MIR emission arising from the dust structure heated by the accretion disk, and not suffering any further extinction. The 6 μm luminosities L_6 μm are calculated for each quasar based on the extrapolation of near- and mid-IR measurements from the WISE All Sky survey (Wright et al. 2010). To perform this calculation consistently for all QSOs we thus require a second cut, selecting only sources with WISE detections at SNR > 3 (W1 − W3). The implications of this selection have been investigated by Rosario et al. (2020), showing this cut only focuses our study on the most powerful AGN, without any significant bias on colour. Recently, a study of quasar bolometric luminosities from z = 0 − 6 (Shen et al. 2020) indeed shows that the infrared luminosity is ideal for such a study as the required correction factor from the monochromatic to bolometric luminosities remains almost constant across the whole redshift range, in particular for QSOs of luminosities log L_bol > 44 ergs s⁻¹. The number of optically (SDSS) selected QSOs at z < 2.5, with WISE detections in at least three bands is 218 747 in total.

The SDSS QSO catalogues also provides photometric information from the SDSS ugriz bands. Based on the g^* (4770 Å) and i^* (7625 Å) band photometry, our third selection criterion uses the g^* − i^* colours, to define the populations of red QSOs and a control sample of ‘normal’ QSOs, following the same strategy presented by Rosario et al. (2020) and Fawcett et al. (2020). We define red QSOs as the 10% of the total QSO population with reddest g^* − i^* colours, and control QSOs as the 50% of the total QSO population closest to the median value of the total population colours. In particular, to account for the redshift dependence of a colour selection, the red and control samples were defined by binning the g^* − i^* distribution as a function of redshift, in contiguous redshift bins of 1000 sources each. Following this classification scheme, from the 218 747 optically (SDSS) selected QSOs at z < 2.5, 21 800 are classified as red QSOs, 109 000 are classified as control QSOs.

While the entire QSO catalogue covers a region of 9376 deg², for the purposes of this investigation we select QSOs located within six survey fields which cover ∼175 deg² in total and include three surveys in the northern sky: the Boötes field, the Lockman-Hole Field, and the ELAIS-N1 Field, and three areas in the equatorial sky covered by the GAMA surveys (Driver et al. 2011): G9, G12, and G15. These fields were selected for two reasons. First, they have rich multiwavelength data from multiple surveys from the FIR to the UV, including deep Herschel coverage, some of which were recently compiled by the Herschel Extragalactic Legacy Program (HELP; see more details in Sect. 2.3), and second, all these fields have deep LOFAR surveys, including the LoTSS Deep Fields surveys (Tasse et al. 2021; Sabater et al. 2021) and the LOFAR GAMA surveys (Williams, Hardcastle et al., in prep.). The total number of selected QSOs defined as red or control QSOs in each field are summarised in Fig. 1, comprising 1483 control QSOs and 306 red QSOs across all six fields, making a total of 1789 sources included in this study.

2.2. Lz-matched subsample selection

In Fig. 2 we show the redshift and L_6 μm distributions of red and control QSOs that result from the sample selection described above. In the central panel red QSOs are represented by red circles, and as red dashed histograms in the side panels. The control QSOs are shown as sky-blue circles and sky-blue filled histograms in the side panels. The single parameter histograms in the side panels show that our selection has produced samples of slightly different redshift and L_6 μm distributions. In particular, the difference is significant at redshifts z > 1.5, where red QSOs appear to peak, and in connection with this, a larger fraction of the selected red QSOs have 6 μm luminosities larger than log L_6 μm = 45.5.

Fig. 2. Monochromatic 6 μm luminosities (L6 μm) of red QSOs (red filled circles), control QSOs (sky-blue filled circles) and Lz-matched control QSOs (dark-blue circles) as a function of redshift are plotted in the central panel. Side panels show normalised histograms of the L6 μm and redshift distributions for each population, shown as solid lines (red QSO), dashed lines (L − z matched samples) and shaded areas (entire control sample) of respective colours. By comparing the distributions of the red QSOs and the Lz-matched control QSO sample, it can be clearly seen that equivalent QSO properties are being sampled.

In order to overcome these disparities and guarantee comparable quasar properties between red and control samples, that is to day equivalent distributions of luminosities and redshift, we assume a conservative approach and define a sub-sample of control QSOs that is directly matched in redshift and luminosity to our selected red QSO sample. Given the larger number of control QSOs compared to red QSOs (5 times larger), we select the closest control QSOs in L_6 μm and z for each red quasar, matched using the Nearest Neighbour method defining a two-dimensional metric on the L_6 μm − z space. Using this approach we define the L_6 μm − z-matched sample, which contain 612 sources, composed of the 306 red QSO, and 306 L − z matched control QSOs. In Fig. 2 we overplot the L_6 μm − z-matched subsample as dark blue circles in the central plot and as a solid dark blue histogram in both side panels. The similarity achieved in the distributions of red and Lz-matched control QSOs is now compelling, providing a sub-sample for our further comparative investigation between red and control QSOs.

We note that the entire analysis included in this paper was performed using both the complete and the L_6 μm − z-matched samples, finding that our main results do not change irrespectively of the sample choice (in agreement to what was previously found in the radio properties by Klindt et al. 2019). Thus, for the purpose of clarity, throughout this paper we use ‘control QSOs’ or ‘cQSOs’ and dark blue markers to refer to the control sample matched in z and L_6 μm to the sample of red QSOs (or ‘rQSOs’), which will be shown as red markers.

2.3. Multiwavelength photometry

For consistency with the sample definition and colour selection, we constructed the multiwavelength SEDs based on the SDSS ugriz bands for the optical regime, and the four WISE bands in the infrared (Pâris et al. 2018). For the three northern Deep Fields we complement these basic SEDs with the aperture-matched photometry compiled for the LoTSS data release (Kondapally et al. 2021). For the equatorial GAMA fields we complement this with aperture-matched photometry from the Herschel Extragalactic Legacy Project HELP project (Shirley et al. 2019), utilising their merged catalogues dmu32¹ to construct the respective multiwavelength SEDs.

We use optical position coordinates to find the multiwavelength counterparts to the SDSS QSOs in our sample, using a search radius of 1 arcsec with respect to the SDSS positions. Given the point-source nature and large luminosities of the optical QSOs, this process is straightforward and overall we find counterparts for 99.9% of our sample. In the three northern fields, 97% of the covered sources (469) have counterparts in the multiwavelength bands. For the remaining sources in these three fields (17) we construct the SEDs based on SDSS and WISE data only. We note that no biases exist towards red or control QSOs for these 17 sources. Given the small number of those without counterparts in the LoTSS catalogue, the impact of the inclusion of these sources on our results is negligible. In the three GAMA fields, we find catalogued counterparts for all QSOs from our main sample.

In the three northern fields, we complement the SDSS photometry with the aperture-matched UV-optical and infrared data available in each field, which includes FUV and NUV from GALEX surveys, and a wide suite of ground-based optical and near-infrared filters reaching out to the K-band. For a detailed description of the photometry in each of the three fields please refer to Kondapally et al. (2021). In the equatorial fields, the HELP catalogues compile calibrated data from a large number of surveys, including the GAMA survey data² in the optical and the NUV and/or FUV GALEX measurements when available. In the near-infrared regime, these are complemented with Y, J, H, Ks bands from both the VISTA and UKIDSS survey.

To cover the MIR regime, the WISE bands are complemented by four IRAC bands at 3.6, 4.5, 5.8 and 8.0 μm and the 24 μm MIPS band in the three northern fields. Since the Spitzer IRAC and MIPS bands are not available in the equatorial fields, here we only use the WISE data. Given that our sample selection already requires strong detections for all our sources in at least three WISE bands, and the IRAC and WISE bands cover similar wavelength regimes, the different sensitivities and lack of IRAC coverage in these field do not significantly affect our conclusions.

Finally, the infrared SEDs of the QSOs in our selection are well sampled by measurements in the 100, 160, 250, 350, 500 μm bands from deep Herschel surveys, including the HerMES project (Oliver et al. 2012) for the three northern fields and the H-ATLAS project (e.g. Valiante et al. 2016) for the GAMA fields. In the northern fields, these have been probabilistically deblended using the XID+ tool (Hurley et al. 2017) as presented by McCheyne et al. (in prep.). In a similar manner, in the equatorial fields HELP computes the deblended Herschel emission for all sources in their catalogues based on optical and MIR priors using the same XID+ method, recovering deblended fluxes in the Herschel bands. While the Herschel flux measurements were deconvolved and robustly assessed, a significant fraction of these are measurements with large uncertainties and in a few cases consistent with zero. In the three northern fields, which were observed with the Herschel PACS and SPIRE instruments with up to 10 times longer exposures per square degree than the H-ATLAS fields (e.g. Lutz 2014), 96.6% have robust Herschel measurements with uncertainties lower than 70%, inconsistent with zero. In the three equatorial GAMA fields, which are part of the H-ATLAS survey, 87% of the sources have Herschel estimates. From these, a significant fraction of sources have large photometric uncertainties and consequently fluxes consistent with zero: ∼23% of the total sample in the case of the SPIRE 250 μm band, and up to ∼59% in the case of the SPIRE 500 μm. In order to properly account for this probabilistic information, and robustly infer star formation properties from such diverse statistical properties a Bayesian method is imperative, and we expand on this in Sect. 2.5.

In total the SEDs for this study consist of 25−30 photometric data points, from the FIR to the UV. We note that we can reliably discard a possible contribution from synchrotron emission to the total SEDs based on available deep low-frequency radio imaging of these fields (LoTSS, Tasse et al. 2021; Sabater et al. 2021; Kondapally et al. 2021). We expand on the radio properties of our sample in a companion paper (Calistro-Rivera, in prep.). To summarise the results relevant to this paper, we find that most QSOs are in the radio quiet regime, in line with the findings by Klindt et al. (2019), Fawcett et al. (2020), Rosario et al. (2020). In particular, we find the synchrotron emission contribution to the FIR-to-UV SED is negligible for the large majority (85%) of the sources in our sample, relatively minor for 11% of the sources (contribution to the infrared luminosities < 10%), and is significant only for the remaining 4%.

2.4. Spectroscopic properties

In addition to the photometry, we compiled spectroscopic information for our QSO sample. The spectra of all SDSS DR14 QSOs were analysed, decomposed and modelled by Rakshit et al. (2020), providing a catalogue of spectral properties which we use for this study. The estimated spectral properties in this catalogue include line fluxes, FWHM estimates for broad and narrow components of several emission lines. In particular in this work we use broad emission line properties from the Hβ, Mg II, and C IV lines, including their FWHMs estimates and associated uncertainties. In Rakshit et al. (2020), the Hβ and Mg II lines were modelled with broad and narrow components, where multiple Gaussians were used to fit the broad component, and a single Gaussian was used to fit the narrow component. During the fitting, the narrow component was allowed to have a maximum FWHM of 900 km s⁻¹ while each of the Gaussians in the broad component had FWHM > 900 km s⁻¹. While the spectral quantities reported for Hβ and Mg II were estimated from the broad component, in the case of C IV the total profile was used without the subtraction of a narrow component, because of the ambiguity in the presence of narrow components in these lines (Shen et al. 2011, 2020). A detailed description of the spectral fitting can be found in Rakshit et al. (2020).

Additionally, we also use the spectral fitting output by Rakshit et al. (2020) on the [O III]λ5007 emission line (through private communication), covered by the SDSS spectral band only for sources with z < 1. A double Gaussian model was used to represent the [O III] emission line; one for the core narrow component and another for the wing component, to characterise potential outflows. More details are discussed in Sect. 4.

Finally, we inspect continuum measurements reported by Rakshit et al. (2020) for comparison to our own estimates from SED fitting, such as bolometric luminosities and monochromatic luminosities. For control QSOs we find a general agreement, as expected. Since the work by Rakshit et al. (2020) does not include corrections for QSO reddening, the reported bolometric luminosities are underestimated for red QSOs when compared to the reddening-corrected continuum results. In Sect. 3.2 we discuss these differences and use the reported emission line properties together with reddening-corrected continuum information from our SED-fitting analysis, to estimate virial M_BH and Eddington ratios for control and red QSOs.

2.5. Fitting spectral energy distributions with AGNfitter

To infer physical properties of the active nuclei and their host galaxies we model the panchromatic spectral energy distributions constructed in Sect. 2.3 using an updated version of the SED-fitting code AGNFITTER (Calistro Rivera et al. 2016). The total AGN SED model in AGNFITTER consists of four physical components: the QSO accretion disk, the dusty torus, the stellar populations, and the galaxy cold dust component. The accretion disk emission model consists of a modified version of the empirical template by Richards et al. (2006) attenuated by a free reddening parameter E(B − V) which is applied following the Small Magellanic Cloud reddening law by Prevot et al. (1984). The dusty torus emission is modelled using the model by Silva et al. (2004), which includes a library of 80 torus templates parametrised by column density in the range of 10²¹ cm⁻² < N_H < 10²⁵ cm⁻², which is associated with the inclination angle. While the fitting machinery and models employed for the inclusion of the host galaxy contribution are mostly the same as those described by Calistro Rivera et al. (2016), we made a few changes to the publicly available version which we describe below.

The new version of AGNFITTER (AGNFITTER2.0, in prep.) was used which increases the flexibility for the easy inclusion of new models. One new model included in this version is the cold dust SED library developed by Schreiber et al. (2018). These models were constructed using the dust continuum models from Galliano et al. (2011) and have been calibrated to include high redshift (z = 2 − 4) galaxies. This model depends on three free parameters: the dust mass (M_dust), the dust temperature (T_dust), and the mid-to-total infrared colour (IR8 = LIR/L8), which measures the relative contribution of polycyclic aromatic hydrocarbon (PAH) molecules to the total infrared luminosity. The main advantage of these models compared to previously used models lies especially in the MIR regimes, where the Schreiber et al. (2018) library allows for higher flexibility in the PAH contributions, which is necessary for star forming galaxies at higher redshift. A second advantage of these models lies in the possibility of exploring a different parameter space than our previous low-z calibrated models, where the dust temperature is of high relevance as it has been observed to increase at higher redshifts (Magnelli et al. 2014; Béthermin et al. 2015). A second variation in the model libraries is the inclusion of a larger and more detailed template library of the stellar population synthesis model. For this we use Bruzual & Charlot (2003), assuming a Chabrier (2003) initial mass function and an exponentially declining star-formation history parametrised by the timescale τ as ψ(age)∝e^−age/τ. The Bruzual & Charlot (2003) model library in AGNFITTER has now been extended to include metallicity as a free parameter, in contrast to the previous version of the code where our library assumed a solar metallicity. More details on this model can be found in Calistro Rivera et al. (2016).

The new version of AGNFITTER also includes the addition of three informative priors, specifically tailored for this QSO population. A prior on the maximum expected contribution of the host galaxy at rest-frame UV has been updated, based on the expected unobscured UV galaxy luminosity functions presented by Parsa et al. (2016), which are calibrated using data from low to high redshift. This prior is applied on the ratio of accretion disk to galaxy contribution ratios at 1500 Å, in order to avoid unphysically large stellar masses, giving preference to a higher AGN contribution when the rest-frame data at 1500 Å is 10 times higher than the maximum expected by the Parsa et al. (2016) galaxy luminosity function at the given redshift. To achieve a conservative and smooth prior we use a wide Gaussian function, which applies this condition unless the data strongly prefers the galaxy contribution.

The probability that the accretion disk emission in QSOs outshines the host galaxy contribution in the optical is large, thus limiting the inference of the stellar mass parameters in the case of the most powerful QSOs. Due to this we observed that the fitting of a fraction of the QSO SEDs in our sample resulted in stellar masses that are unrealistically low for the host galaxies of such powerful SMBHs. In order to use information on a lower limit of stellar masses in a conservative manner, we applied a prior of a wide Gaussian shape on the normalization parameter of the host galaxy (GA) with a mean value of GA = 4.5 and a standard deviation of σ = 1.5. We choose these values by taking into account that the normalization parameter GA correlates strongly with the output on total host-galaxy stellar-masses, where GA < 3 corresponds approximately to total stellar masses of M_* < 10⁹ M_⊙, an unlikely scenario for our sample, with a scatter that depends on the variety of SED shapes. We tested the robustness of this prior finding clearly different prior and posterior distributions for the parameters. This shows that while supporting the estimation of physical stellar masses, the prior does not entirely determine the posterior probability distribution but there is significant constraining power from the data.

Finally we include a conservative prior to consider an energy balance between the dust-attenuated emission from the stellar populations and the reprocessed cold dust emission. In particular, we assume that the energy budget emitted by the cold dust in the infrared is at least as high as the energy comprised by the attenuated stellar emission. In contrast to other methods used in the literature, we adopt this conservative approach, so that our model is flexible enough to account for uncertainties evident in dusty star forming galaxies; in these sources, large uncertainties on light to mass ratios, accuracy of dust attenuation laws, and the increasing number of observed uncorrelated distributions of dust and stars, may challenge more stringent energy-balance assumptions (e.g. Calistro Rivera et al. 2018; Buat et al. 2019).

We apply the SED-fitting routine on the entire sample of 1789 QSO SEDs, achieving high quality fits (ln_like < 50) for 92% of the sources, and fair quality fits for another 5%. Although we fit the entire sample, as stated in Sect. 2.2 we consider only the red and control QSO samples matched in redshift and luminosity for the analysis and interpretation of the results in Sects. 3 and 4. Examples of the fitting results can be found in Fig. 3, for a few cases of control and red QSOs. We note that an important assumption in the physical modelling in AGNFITTER is that any variations in the optical-UV shapes of the accretion disk emission SED are explained by reddening through dust attenuation. This assumption would be wrong if the red colours would be due to intrinsic accretion properties instead. However, extensive evidence exists in the literature (e.g. Richards et al. 2003; Glikman et al. 2012, 2013; Kim & Im 2018; Wethers et al. 2018; Klindt et al. 2019; Temple et al. 2021) which supports the dust-reddening scenario as the most plausible description of the data. In particular, the QSO reddening model in AGNFITTER follows a Small Magellanic Cloud (SMC) attenuation law, which is assumed to be a good description for the reddening in Type 1 SDSS QSO spectra (e.g. Glikman et al. 2012). In addition to the mentioned studies, in a forthcoming investigation (Fawcett et al., in prep.) we apply detailed spectral fitting using different reddening models to X-shooter UV–near-IR spectroscopy of a sample selected in the same manner as this study, demonstrating that the dust-reddening scenario is a robust description for our QSO sample. Finally, using SDSS spectroscopic data and reddening-corrected fluxes, we find no significant differences in the intrinsic properties of the red and control QSOs (see Sect. 3.2 for more details), validating our dust-attenuation approach.

Fig. 3. Examples of the SED decomposition with AGNfitter. In the main panels the photometric data points are shown with black error bars. The best fit and ten realisations of the posterior distributions of the fit are drawn as coloured solid, and transparent lines, respectively. The realisations display the uncertainties associated with the fits. The different colours represent the total SED (red lines), the accretion disk emission (blue lines), the host galaxy emission (yellow lines), the torus emission (purple lines) and the galactic cold dust emission (green lines). Lower panels: residuals expressed in terms of significance given the data noise level for the best fit and the 10 realisations. Through this selection of SED fits we exemplify the diversity of the obscuration levels found. Two examples of red QSO SEDs which clearly present an infrared excess emission at ∼3 μm, as discussed in Sect. 3.3.2, are included in the two lower left panels.

2.6. Statistical comparative methods

The Bayesian method inherent to the AGNFITTER code allows us to sample the probability density functions for all of the output physical parameters. The AGNFITTER output therefore includes the parameter posterior PDFs for each source sampled using 100 random realisations. To achieve robust comparisons of the properties of red and control QSOs in this study, we need to account for the uncertainties and shapes of the posterior PDFs of each source. To this end we produce composite probability distributions for the red QSO and control QSO populations, as a superposition of the PDFs for the single sources within each group.

For purposes of visualisation, we apply a Gaussian kernel-density estimator (KDE) with a given bandwidth (specified in each plot). To assess the significance of the differences between the distributions we apply a bootstrap version of the non-parametric two-sample Kolmogorov–Smirnov (KS) test. This method consists in defining bootstrapped pairs of red and control QSO samples and calculating the KS probability through the application of the two-sample KS test. We note that, throughout Sects. 3 and 4, we consistently use red and control samples matched in redshift and luminosity for the comparisons, if not otherwise stated. To quantify the distribution of KS p-values that arise from the whole composite distributions, we perform 3000 bootstrap iterations, after testing for convergence of the results at this iteration number. Finally, the median and percentiles values of the KS p-values represent the probability that the parameter distributions for the two samples are drawn for the same distribution, i.e. the distributions have no significant differences. In particular, median KS p-values of p < 0.05 would imply that we find statistically significant differences in properties of red and control QSOs (Sects. 3 and 4).

3. Results

3.1. The far-infrared–UV SEDs of blue and red QSOs

We show average SEDs for the different QSO populations in Fig. 4. All have been renormalised to have the same L_6 μm for purposes of visualisation. Although there is a significant scatter, the median SEDs depicted as solid lines for both red and control QSOs show remarkably similar SEDs across the electromagnetic spectrum, starting to show differences only in the optical/UV regime consistent with the shape expected for red QSOs (Sect. 2).

Fig. 4. Ten modelled SEDs are plotted for each source (reconstructed from 10 randomly-picked realisations of the posterior PDFs) as transparent red (red QSOs) and blue lines (control QSOs) to display the uncertainties. Solid red and dotted blue lines depict the composite median SED for the red and control QSO population respectively, matched in redshift and L6 μm. At frequencies log ν < 14.5 Hz, strikingly similar composite SEDs can be observed for the red and control QSO populations. At frequencies log ν > 14.5 Hz a clear difference in the composite SEDs can be observed which is the signature of dust reddening consistent with the optical colour selection.

To understand in detail the origin of the SED diversity, in Fig. 5 we present average SEDs for each individual physical component, normalised at rest-frame 6 μm, and binned in groups of redshift and rest-frame 6 μm luminosities. A quick visual inspection of the overall SEDs (grey solid and dashed lines) reveals the prevalence of flatter SEDs at lower redshift, where lower QSO luminosities are sampled, whereas SEDs at higher redshifts present ‘deeper’ features in the near-infrared. These deeper features are seen at high redshifts since the host galaxy contribution is no longer significant and does not fill the near-IR emission bands. Indeed, a clear advantage of our SED decomposition approach is the ability to disentangle the host galaxy contamination which is more prominent at low-redshifts, and in particular in the low-luminosity bins. Our study shows that for redshifts of z < 1.6, the host galaxy emission can have a similar contribution to the near-infrared SED as the dusty torus and accretion disk components, in particular at luminosities of log L_6 μm ≤ 45.

Fig. 5. Median composite SEDs are shown for the red QSO sample (solid lines of different colours) and matched control QSO samples (dotted lines of different colours), now split into different bins of L6 μm (left-right) and z (top-bottom). All composite SEDs have been renormalised to the same L6 μm for purposes of visualisation. The grey solid and dotted lines represent the total SEDs for the red and control samples, respectively. The remaining colours codify the different physical components as described in the legend of Fig. 3. The luminosity and redshift bins considered and the number of sources included in each bin are annotated. The contribution of the stellar emission of the host galaxies has a clear contribution to the total SEDs in the lower redshift bins and to luminosities of L6 μm < 45. A remarkable similarity in the composite torus SEDs of red and control QSOs across all redshift bins can be seen.

A significant decrease of the relative contribution of the stellar population component with respect to the accretion disk emission is found as a function of redshift and luminosity. In line with the average result in Fig. 4, the median SEDs within different redshift and luminosity bins also show overall similar SEDs across the electromagnetic spectrum apart from the optical/UV regime.

A key result from the SED decomposition is the remarkable similarity between the torus SEDs for red and control QSOs, shown by the solid and dashed purple lines in Fig. 5, across all redshift and luminosity bins. This result hints towards no differences with respect to the dusty torus component between the control and red QSOs, suggesting the origin of the reddening does not arise from inclination. A more detailed analysis of obscuration properties is carried out in Sect. 3.3.

3.2. Intrinsic accretion properties

To further test whether the optical-colour diversity used for the selection of our samples arise from intrinsic black hole accretion properties, rather than dust attenuation, we estimate black hole masses, bolometric luminosities, and Eddington ratios. While estimates of these properties were provided for the entire SDSS catalogue by Rakshit et al. (2020), these estimates have not been corrected for accretion disk reddening. Our SED-fitting approach allow us to include these corrections and produce more reliable comparisons.

We estimate black hole masses for our sample adopting the same approach employed by Rakshit et al. (2020) and Shen et al. (2011), with

$$\begin{array}{c}\hfill log\left(\frac{M_{\mathrm{BH}}}{M_{\odot }}\right)=A+Blog\left(\frac{\lambda L_{\lambda }}{{10}^{44}\mathrm{erg}{\mathrm{s}}^{-1}}\right)+2log\frac{\mathit{FWHM}}{{\mathrm{km}\mathrm{s}}^{-1}},\end{array}$$(1)

where λL_λ are the monochromatic luminosities L₁₃₅₀, L₃₀₀₀, L₅₁₀₀, estimated from the reddening and host galaxy corrected SEDs. The FWHMs are obtained from [C IV], [Mg II] and Hβ lines, reported by Rakshit et al. (2020) for the SDSS QSOs at z < 0.8, 0.8 < z < 1.9, and 1.9 < z < 2.5 respectively. The constants A and B³ are given by Rakshit et al. (2020) as well, where they use corrections based on virial estimates to the corresponding lines inferred by Vestergaard & Peterson (2006) and Vestergaard & Osmer (2009). We note that the monochromatic luminosities from our SEDs agree with those reported by Rakshit et al. (2020) for the control and red QSOs before the reddening correction is applied. After the reddening correction, the monochromatic luminosities of red QSOs are significantly larger than those reported by Rakshit et al. (2020), as expected. While red QSOs would have lower median M_BH values than the control sample based on black hole masses reported by Rakshit et al. (2020) (Δ log M_BH ∼ 0.25), this effect disappears after applying the reddening correction (Δ log M_BH < 0.1). In Fig. 6 both populations have similar distributions with almost identical median M_BH values at log M_BH ∼ 9 M_⊙ and a slight difference in the width of the distributions. These results are in line with the emission line width and luminosity analysis presented by Klindt et al. (2019), where they investigate the intrinsic accretion properties by comparing line width distributions for blue, control and red QSO, finding no significant differences.

Fig. 6. Distributions of black hole mass estimates (MBH), reddening-corrected bolometric luminosities Lbol and Eddington ratios λEdd for the control and red QSO samples. Only MBH estimates from spectra with quality tag 0 are considered, which comprise 90% of our entire sample. A Gaussian kernel-density estimator (KDE) with a bandwidth of 0.3 is used to display these distributions. With exception of the differences in the bolometric luminosity distributions, where significance is weak, no significant differences are found between the red and the control QSOs in these parameters. The significance of the difference between the two distributions (KS-test p-value) is reported below. The p-values were computed using the estimates reported by Rakshit et al. (2020) (MBH and λEdd), and using the bootstrapped KS-test for the Lbol, where the posterior distributions from the MCMC-based SED fitting are compared and therefore uncertainties could be calculated.

Bolometric luminosities are also estimated from the integrated emission from the accretion disk component (BBB) over the wavelength range of 500 Å−1 μm, with a small additive correction factor of δlog L_bol = 0.3 to account for X-ray emission not included in the SEDs. The bolometric luminosities reported by Rakshit et al. (2020), on the other hand, are computed using a bolometric correction to the monochromatic luminosity at 3000 Å, following L_bol = 5.15 × L₃₀₀₀ (Shen et al. 2011). We find that our SED-based bolometric luminosity distribution is equivalent to that reported by Rakshit et al. (2020) for the entire sample before the reddening correction is applied. After the application of the reddening correction, however, our estimates for the red QSOs are significantly larger, as expected. The second panel of Fig. 6 shows the L_bol distributions, where slightly larger bolometric luminosities for red QSOs can be recognised. A bootstrapped KS-test (Sect. 2.6) shows these differences are small though statistically significant (p-value = $0.{01}_{-0.01}^{+0.09}$). We explore these implications in Sect. 3.3.

Based on our estimates of M_BH and L_bol we calculate Eddington luminosities and ratios for our comparative samples of red and control QSOs in the right panel of Fig. 6 following,

$$\begin{array}{c}\hfill L_{\mathrm{Edd}}=3.28\times {10}^4\left(\frac{M_{\mathrm{BH}}}{M_{\odot }}\right)L_{\odot }{\lambda }_{\mathrm{Edd}}=\frac{L_{\mathrm{bol}}}{L_{\mathrm{Edd}}}·\end{array}$$(2)

Once we have accounted for host galaxy and dust reddening corrections, we find strikingly similar distributions of black hole masses and Eddington ratios for red and control QSOs.

Both QSO populations have median Eddington ratios around λ_Edd = 0.1, suggesting that the intrinsic accretion properties of the SDSS QSOs included in this study are very similar for both red and control QSOs. Based on these observations, we conclude that we find no differences in the intrinsic accretion properties of red and control QSOs.

3.3. Obscuration properties of blue and red QSO

To investigate the origin of the red colours in our QSO sample we focus on two physical components of our SED fitting model: the reddened accretion disk emission itself, and the reprocessed emission by hot dust at nuclear scales commonly modelled with a torus structure. Two output quantities that parametrise these components are the dust extinction parameter that attenuates the blue part of the accretion disk emission, E(B − V)_BBB, and the column density of the obscuring medium inferred from the torus infrared SED, N_H. In Fig. 7 we show the posterior distributions of the dust attenuation of the accretion disk SED E(B − V)_BBB for red QSOs and control QSOs. The posterior distributions are reconstructed by the superposition of 100 draws from the PDF of each source.

Fig. 7. Probability distributions of the accretion disk reddening parameter E(B − V), and the column density NH. The E(B − V) parameter was inferred assuming the SMC attenuation law by Prevot et al. (1984) on the accretion disk emission, while the NH is a quantity that parametrizes the (X-ray-normalised) infrared torus templates by Silva et al. (2004). These distributions are constructed as the superposition of 100 random draws from each source. A gaussian kernel-density estimator (KDE) with a bandwidth of 0.1 is used to display these distributions.

As expected by their selection, we find a significant difference in the distributions of E(B − V)_BBB for red and control QSOs (bootstrapped KS-test p-value ∼ 10⁻⁶⁵), where red QSOs have $E{(B-V)}_{\mathrm{BBB}}=0.{12}_{-0.08}^{+0.21}$ mag ($A_V=0.{32}_{-0.21}^{+0.57}$ mag), and control QSOs have $E{(B-V)}_{\mathrm{BBB}}=0.{02}_{-0.01}^{+0.04}$ mag ($A_V=0.{05}_{-0.03}^{+0.11}$ mag).

In contrast to what is expected from the dust attenuation, we find that both red and control QSOs have strikingly similar distributions of dust column density N_H as measured by the torus emission component (Fig. 7). The bulk of both populations (∼83%) show low column densities of log N_H < 22., with no significant differences (bootstrapped KS-test p-value = 0.17). This is in line with the similar shapes of the dusty torus SED components found among red and control QSOs shown in Fig. 5, suggesting that the main infrared emitting component is similar for the two populations. We remark that the log N_H values quoted here are not inferred from the X-rays but from the infrared alone, using the log N_H term to parametrize the torus templates by Silva et al. (2004). These estimates of N_H are, however, overall consistent with measures from X-rays spectroscopic studies of SDSS QSOs (Wilkes et al. 2002; Urrutia et al. 2005), although these are found to be significantly higher for red QSOs selected in the near-infrared (Goulding et al. 2018; Lansbury et al. 2020).

The similarity in the shape of the torus SED suggests that the nuclear dust structure, which is responsible for the AGN unification scheme, has no direct connection with the dust attenuation that characterises the red QSO population. Indeed, would the reddening arise from a viewing angle effect, in particular from viewing angles closer to an edge-on view of the dust structure (Type 2 QSO), this would imprint the torus SED shape due to the higher observed column density. The lack of diversity of the torus SEDs for all red and control QSOs, as inferred by our model, suggests that the torus viewing angle effect is not causing the reddening. One systematic limitation that would synthetically produce this observation is the possibility that the torus emission templates used for the SED fitting do not recover the effect of torus inclination properly. However, our previous work has shown that this is not the case. Using the same machinery and semi-empirical templates (Silva et al. 2004) on spectroscopically confirmed Type 1 and Type 2 QSO, in Calistro Rivera et al. (2016) we demonstrated that these models successfully recover a diversity of torus structures. In particular, these were able to recover the spectroscopic Type 1 and Type 2 QSO classification based solely on the shape of the SEDs, using the diversity found in the N_H and E(B − V) parameters as the drivers of the classification. Additionally, the column density distributions and torus SED shapes observed in Type 1 QSOs in the COSMOS field by Calistro Rivera et al. (2016), are fully consistent with the distributions observed by the SDSS QSOs in this work, independently of colour.

In the following sections we explore two scenarios that explain the lack of differences in the infrared emission of red and control QSO, despite the clear difference in the dust-attenuated optical and UV emission. Firstly, we investigate the dust reprocessing efficiency in Sect. 3.3.1. Secondly, in Sect. 3.3.2 we investigate whether there is emission that is not recovered by our models by computing the infrared residuals from the SED fitting.

3.3.1. Reprocessing efficiency

In Fig. 8 we investigate the fraction of the intrinsic bolometric AGN luminosities that is reprocessed by the dusty torus and re-emitted in the infrared, i.e the torus reprocessing efficiency, which is used as a proxy for the covering factor of the obscuring medium (Maiolino et al. 2007; Treister et al. 2008; Lusso et al. 2013; Stalevski et al. 2016). We compute the AGN bolometric luminosities L_bbb and torus luminosities L_tor by integrating over the fitted templates of the accretion disk and torus physical components, covering the range between 0.05 and 1 μm, and 1 μm, and 1000 μm, respectively. We note that our results do not change significantly if the range 0.1−1 μm is used for the accretion disk emission. We also stress that the SED decomposition method we use to compute the L_bbb is more robust than bolometric luminosities scaled up from one or a few wavelength measurements using a single bolometric correction, since more information/data is included in the estimation.

Fig. 8. Composite two-dimensional posterior distributions of the torus reprocessing efficiency, defined as the ratio of torus luminosities and bolometric luminosities, plotted as a function of bolometric luminosities. These are constructed by sampling 100 realisations from the posterior of each source to properly account for the uncertainties. The red contour lines (for red QSOs) and the blue contour areas (for control QSOs) range between 20 and 100%, in intervals of 20%. Upper panels: this ratio for the observed bolometric luminosities (left panel) and for the reddening-corrected (intrinsic) bolometric luminosities (right panel) for red and control QSOs. Lower panel: torus reprocessing efficiency corrected for the torus anisotropy following Stalevski et al. (2016), to provide a better description of the covering factor.

We compute the torus reprocessing efficiency as the ratio of the accretion disk emission (bolometric luminosity) and torus luminosities, considering both the observed accretion disk luminosity (L_tor/L_{bbb − obs} Fig. 8, left) and the reddening-corrected one (L_tor/L_{bbb − dered}, Fig. 8, right). This ratio is related to the covering factor, which is in turn connected to the probability to observe an AGN as obscured or unobscured (Lusso et al. 2013; Stalevski et al. 2016). We remind the reader that, in contrast to other SED fitting models, AGNFITTER models the accretion disk and torus components independently, allowing to investigate the variations in this ratio.

In the left upper panel of Fig. 8, we show this ratio as a function of the intrinsic accretion disk luminosities (reddening corrected) for our QSO samples. A clear difference between red and control QSOs can be recognised, which is expected from the dust attenuation in red QSOs. We next test whether this observation holds after correcting the bolometric luminosities for reddening of the accretion disk. In the right panel of Fig. 8, we show this ratio as a function of the intrinsic reddening-corrected accretion disk luminosities. We find that the distribution of the log L_tor/L_bol has significantly changed after the reddening correction, now around physical values below unity, and that red QSOs are more similarly distributed to the control QSOs. We find that the median value of the torus reprocessing efficiency is 0.437 ± 0.002 for red QSOs, which is lower than that observed for the control sample, at values of 0.550 ± 0.002, where uncertainties are bootstrap errors. The linear fits reveal slope values of −0.35 for red QSOs and −0.37 for the control sample. These values reflect an apparent negative evolution of the obscuration in both QSO populations as a function of bolometric luminosity, consistent with the receding torus scenario, where the increasing radiation pressure of higher QSO luminosities may push the obscuring medium away from the luminous accretion disk (e.g. Simpson 2005, but also note caveats and different interpretations in Netzer et al. 2016; Stalevski et al. 2016). We perform a non-parametric two-sample KS statistic described in Sect. 2.6 to test the hypothesis that the two distributions of red and control QSOs are drawn from the same distribution of reprocessing efficiency. We find we can reject this hypothesis with a confidence level of 99.9% (p-value = 2 × 10⁻⁵) and therefore conclude this difference is statistically significant. The tentative observation, supported by the KS-test, that red QSOs have lower reprocessing efficiency would hint towards a slight difference in the torus composition, possibly related to larger gas to dust fractions and a more diffuse obscuring medium (Lansbury et al. 2020).

Finally, we estimate the torus covering factor by accounting for the torus anisotropy and its effect on our estimates of the reprocessing efficiency. Based on a 3D radiative transfer model, Stalevski et al. (2016) quantified this effect, showing that its influence can be non-linear and strongly dependent on the assumed torus optical depth. Motivated by the output on torus column densities from the SEDs, we assume a torus optical depth of τ = 3 for both red and control QSO samples, and apply the corrections reported by Stalevski et al. (2016) to account for the torus anisotropy. We show the resulting relation in the lower panel of Fig. 8. The differences between distributions of red and control QSOs decrease after this correction, however they remain significant (p-value = 2 × 10⁻⁵). Nonetheless, whether it is indeed the case that red QSOs have less efficient torus obscuration, or if we adopt the more conservative interpretation that both populations have equivalent efficiencies based on the weak differences, none of these scenarios can explain the origin of the red QSO colours, reinforcing our previous hypothesis that no connection exists between QSO reddening and their dusty tori.

Possible systematic effects, such as a uncertainties in choosing appropriate extinction curves to describe QSO reddening might affect our results, especially where no data points are available to constrain the curve. However, the E(B − V) distributions from the Prevot et al. (1984) model are consistent with what has been found in other studies, and the quality of the fits show SMC laws are a good description of the reddening in these sources (e.g. Hopkins et al. 2004; Glikman et al. 2012; Fawcett et al., in prep.). Finally, as mentioned above, we verified that a choice of a different wavelength range (0.1−1 μm, instead of 0.05−1 μm) does not change our results significantly. Another possible systematic effect is the choice of intrinsic accretion disk emission model, which could over-predict the intrinsic luminosities for red QSOs. We expect that the effect of such a scenario is not significant, since the intrinsic accretion model of Richards et al. (2006) performs well for the control sample, and no differences in the intrinsic accretion properties are expected based on our spectral analysis in Sect. 3.2 (Fig. 6). Overall, we find the infrared reprocessing efficiency of red QSOs is slightly lower or at least comparable to the control QSO sample, suggesting the torus structure is disconnected to the reddening properties of red QSOs.

3.3.2. Residual dust distributions

We now investigate the existence of infrared features that could be linked to the dust attenuation in red QSOs but are not recovered by torus models. To quantify such emission we compute the residuals of the best fits for the QSO SEDs in the infrared, expressed in terms of significance given the data uncertainties. We plot these as a function of rest-frame log frequency in Fig. 9. Although the residuals are overall small and not significant in the majority of the frequencies covered, some structure can be recognised in the mid-infrared regime. As expected, the SEDs do not have coverage in the rest-frame wavelength regime between λ ∼ 20 − 30 μm. Interestingly, two peaks are recognised in the regions highlighted with colours in Fig. 9 where significant positive residuals are found (σ > 3). Positive residuals mean that there is emission from the data which cannot be recovered by the models.

Fig. 9. Residuals from the best fits for each source as a function of rest-frame frequency. Two small ‘bumps’ are identified in the MIR regime, where the torus models included in our model show the largest discrepancy with the data. The main discrepancies have positive residuals meaning the model cannot reproduce a part of the flux at these frequencies. Right panels: cumulative distributions of the incidence of residuals in each sample. We investigate whether these residuals are connected to the optical colour of the QSO, despite the fact these bands are independent in the fitting. We find that red QSOs have larger MIR emission excess compared to the model in the regime of 13.75 < log ν < 14.2 (2−5 μm), showing significant residuals ( > 2σ) in 20% of the cases, in contrast to the 10% in the case of control QSOs.

We investigate the origin of these peaks and compute the incidence of residual luminosities for each QSO sample in the right side panels of Fig. 9. We find that there is a larger incidence of residuals in the SEDs of red QSOs as compared to control QSOs, which is particularly present in the infrared regime of rest-frame log ν ∼ 14 Hz (2−5 μm; panel A). A KS test on this difference reports a p-value of 0.0006, suggesting the enhanced flux excess in red QSOs is highly significant. We further evaluate the fraction of sources with residuals above this value. We find that 20% (54/306) of the red QSO sample have residuals of σ > 2, whereas the fraction reduces to 11% (35/306) for the control QSO population. A few examples of red QSO SEDs which clearly present these residuals were also included in Fig. 3 (two lower left SEDs).

A slight enhancement of infrared residuals can be recognised as well in the rest-frame log ν ∼ 13.5 Hz regime (λ ∼ 10 μm), albeit with lower significance. With black body temperatures at ∼150 − 500 K, the residuals in this longer-wavelength regime suggest that warm dust emission contributions might be present for Type 1 QSO, irrespective of colour. For the purposes of this investigation we restrict our further analysis to the infrared residuals in the 2−5 μm regime.

We further investigate a direct relation between the incidence of infrared residuals in the 2−5 μm regime as a function of QSO reddening in Fig. 10. A clear trend can be observed, where the infrared residuals increase with larger reddening values for all samples, in spite of the fact that both parameters (E(B − V)_BBB and residuals) are determined by independent regions of the SED (rest-frame near-infrared and optical/UV, respectively). This observation offers additional support to a potential link between this excess hot dust emission and the colours in red QSOs. It is important to note that the fact that the infrared residuals are characterised by low levels of σ is mainly determined by the low signal-to-noise of the infrared photometric measurements which is our only means to parametrize such excess.

Fig. 10. Median infrared residuals are shown as a function of the QSO reddening parameter, E(B − V)BBB. Black circles represent the median residual values for each E(B − V)BBB bin for all QSOs irrespective of colour classification. The error bars represent bootstrap 1-σ errors. The same relation is shown for the red and control QSOs as red and blue shaded areas, respectively. A clear trend can be observed, where the infrared residuals increase with larger reddening values for all samples.

Excess emission at this wavelength regime (2 − 5 μm) corresponds to dust temperatures of 550−1700 K, assuming a black body and applying Wien’s law. In the near-infrared regime, we should also consider the potential contribution from the host galaxy. However, we can rule out this possibility for the following reasons. First, based on the general SED of galaxies it is highly unlikely that the host galaxy SED presents a peak around 2 − 5 μm, since the reddest possible peaks for the most massive, older and cooler stellar populations are located at wavelengths bluer than 1.5 μm. Instead, galaxy SEDs show overall flat SEDs (Bruzual & Charlot 2003) in the wavelength regime relevant here, being unable to explain the red excess bumps around λ ∼ 4 μm. Most crucially, the overall galaxy contribution is negligible in the majority of our sources due to their high-luminosities, as seen in Fig. 5.

We conduct a few exercises to test the robustness of this finding. First we investigate whether the higher infrared excess in red QSOs is driven by a fraction of sources at a specific redshift. We find this is not the case, finding a clear enhancement of infrared residuals for red QSOs compared to the control samples across most redshift ranges covered by the study, with p-values (for the respective redshift ranges) of 0.01 (0.2 < z < 0.9), 0.16 (0.9 < z < 1.6), 0.02 (1.6 < z < 2.2), and 0.001( 2.2 ≤ z < 2.5). We further test whether this excess is driven by the observations of a specific instrument (WISE or IRAC). Although we find a significant enhancement in both samples which include WISE+IRAC (northern fields) and those with only WISE data available (GAMA fields), it is clear that those with deeper measurements and a better sampling of the SED around rest-frame 3 μm have excess detections of higher significance. In particular, since the bulk of the QSO population is located at z > 1 the coverage of the IRAC channels ch3 (5.8 μm) and ch4 (8.0 μm) is crucial for detecting this excess. We further test whether this emission could be recovered by more complex state-of-the-art torus models which include other free parameters such as opening angle of the torus, and different dust grain sizes. In particular, we perform a SED-fitting test run on a small subsample of our data with significant residuals using the torus model by Stalevski et al. (2016), finding no improvement in the fit of those features.

Given that these infrared residuals cannot be reproduced by torus SED models we investigate whether these have originated from additional hot dust emission components, such as the presence of dusty winds at radii close to the sublimation radius. We discuss this scenario in Sect. 4. In particular, hot dust models exist which include the emission from polar cones of dusty winds lifted up from the inner parts of the torus by radiation pressure (Hönig & Kishimoto 2017; Costa et al. 2018a; Venanzi et al. 2020). We note that including these models in the SED fitting implies sampling highly complex parameter spaces (11 parameters, Hönig & Kishimoto 2017). This is out of the scope of this investigation given the limitations of the available infrared photometry and the lack of MIR spectroscopic data.

3.4. Host galaxy properties of red and control QSOs

Finally, we also investigate a possible link between the properties of the host galaxies and the reddening of QSOs. In particular, we focus on the star-formation rates, stellar masses, and molecular gas masses estimated from the SED decomposition technique. We do not include parameters such as the galaxy age or galaxy dust attenuation in our analysis since these are constrained based on the bluer parts of the SED which is largely dominated by the QSO emission, making these estimates highly unreliable. Following the same argument, we note that the SFR estimates presented in this work are obscured star-formation rates, i.e. these have been estimated solely based on the integrated luminosities of the infrared regime between 8 and 1000 μm, by fitting the photometric data with the Schreiber et al. (2018) cold dust model. Although estimating stellar masses from photometry for Type 1 QSOs can be a challenging endeavour, the near-infared data has significant constraining power especially for moderate QSO contributions and in the lower redshifts bins. Based on our knowledge of the QSO accretion disk emission constrained by the blue optical bands and UV, and supporting the inference with prior information on expected host galaxy luminosities (see Sect. 2.5), we are able to disentangle the host galaxy component and derive estimates of the stellar masses and their associated uncertainties. We note, however, that in the higher QSO luminosity and high redshift bins (log L_6 μm > 44.6 ergs s⁻¹, z > 1.), where the host galaxy emission starts becoming less important compared to the QSO luminosities, the uncertainties are expected to increase, and we rely on our Bayesian technique to recover these as demonstrated by the statistical tests presented in Calistro Rivera et al. (2016). Molecular masses are calculated using the empirical calibration of the Rayleigh-Jeans luminosity-to-mass ratio by Scoville et al. (2016), based on the SED-estimated rest-frame 850 μm luminosities and associated uncertainties.

In Fig. 11 we examine the composite posterior probability density functions of the relevant parameters: stellar masses, infrared star-formation rates, and molecular gas masses. We construct these composite PDFs, again by combining 100 draws from each source’s PDF and apply a Gaussian KDE for display. The overall stellar mass distribution for red QSOs has a median value of $log{M}_{\ast }(\mathrm{rQSO})=10.{48}_{-1.04}^{+0.69}{M}_{\odot }$ (where uncertainties correspond to the 14th and 86th percentiles) whereas control QSOs have $log{M}_{\ast }(\mathrm{cQSO})=10.{32}_{-0.72}^{+0.72}{M}_{\odot }$. To assess the significance of the difference between these two distributions we explore the probability that these have been drawn from the same distributions, using the KS-test described in Sect. 2.6. We find a median p-value = $0.{01}_{-0.01}^{+0.09}$, suggesting the differences in stellar masses are significant, although marginally given the uncertainties in the p-value. In order to further test the SED-based stellar masses, we compare these to those estimated based on the black-hole-mass-stellar-mass relation, empirically calibrated by Bennert et al. (2011) as a function of redshift. The M_BH-based stellar mass distribution for the QSO samples have median values of $log{M}_{\ast }(\mathrm{rQSO})=10.{82}_{-0.57}^{+0.52}{M}_{\odot }$ and $log{M}_{\ast }(\mathrm{cQSO})=10.{89}_{-0.47}^{+0.32}{M}_{\odot }$. This shows that, although the SED-inferred stellar masses are overall consistent with the M_BH-based estimates within the uncertainties, the M_BH-based estimates do not recover the slightly enhanced stellar masses for red QSOs reported above. We therefore conclude that the overall stellar mass distributions show no difference between red and control QSOs.

Fig. 11. Composite posterior distributions of host galaxy parameters inferred through SED-fitting. These are constructed by sampling 100 realisations from the posterior distributions of each source, and we apply a KDE with a bandwidth for 0.3 for display as explained in detail in Sect. 2. Stellar masses are estimated based on the modelling of near-infrared and optical emission with stellar population models (Bruzual & Charlot 2003). Star formation rates are computed based on total infrared luminosities (LIR(8 − 1000 μm)) from the modelling of the infrared photometry with cold dust models (Schreiber et al. 2018). Molecular gas masses are estimated from monochromatic luminosities of the dust continuum SED at 850 μm and the empirical relation by Scoville et al. (2016). No significant differences are found between the distributions of red and control QSOs for these host-galaxy parameters.

The SFR distributions are also remarkably similar, where the QSO samples have median values of SFR${}_{\mathrm{IR}}(\mathrm{rQSO})={95}_{-85}^{+364}{M}_{\odot }$ yr⁻¹ and SFR${}_{\mathrm{IR}}(\mathrm{cQSO})={70}_{-62}^{+309}{M}_{\odot }$ yr⁻¹. Here we find that the distributions are not significantly different with a median p-value = $0.{05}_{-0.05}^{+0.25}$. Similarly, the molecular masses have median values of $log{M}_{\mathrm{mol}}(\mathrm{rQSO})=10.5_{-1.1}^{+0.9}{M}_{\odot }$ and $log{M}_{\mathrm{mol}}(\mathrm{cQSO})=10.4_{-1.0}^{+0.8}{M}_{\odot }$. We find that the distributions of molecular gas masses are not significantly different with a median p-value = $0.{05}_{-0.05}^{+0.21}$.

Despite this overall agreement of properties, we carry out a more focused analysis for the low-z sample since we can assume the host galaxy properties are more reliably determined in these sources through the SED fitting. In Fig. 12 we plot the M_*–SFR relation for the red QSO sample in comparison to the luminosity-and-redshift-matched control sample for the redshift range 0.2 < z < 0.8, where the overall SEDs had significant contribution from the host galaxy emission. While both distributions are consistent with the main sequence and have similar sSFRs, the red QSOs are distributed in a more massive region ($log{M}_{\ast }(\mathrm{rQSO})=10.{49}_{-1.11}^{+0.43}{M}_{\odot }$) than control QSOs ($log{M}_{\ast }(\mathrm{cQSO})=10.{14}_{-0.70}^{+0.60}{M}_{\odot }$) although with considerable overlap. The region occupied by red QSOs is consistent with that populated by ULIRGs at z < 1 as shown by the median M_*–SFR relation of the GAMA Field ULIRGs (Driver et al. 2018). However, this difference disappears if we use the M_BH-based stellar mass estimates finding values of $log{M}_{\ast }(\mathrm{rQSO})=10.{47}_{-0.59}^{+0.45}{M}_{\odot }$ and $log{M}_{\ast }(\mathrm{cQSO})=10.{44}_{-0.58}^{+0.54}{M}_{\odot }$. While this tentative finding is interesting, it is limited to this redshift bin (only 16% of the red QSO sample) and no conclusions can be drawn on the stellar masses of the general red QSO population based on this observation.

Fig. 12. Composite two-dimensional posterior distributions of star-formation rates versus stellar mass for red QSOs (red shaded contour areas) and control QSOs (blue contour lines) for the low-z sources. These are constructed by sampling 100 realization from the posterior of each source to properly account for uncertainties. Stellar masses and star-formation rates are estimated based on SED modelling. The distributions are consistent with the main sequence (grey shaded region), although as well with ULIRGS (GAMA survey; Driver et al. 2018) within the uncertainties.

Overall, we find no strong differences in the host galaxy stellar populations and infrared properties between the red and control QSO populations. In particular, the infrared luminosities from which star-formation rates were estimated and the molecular gas masses, as measured using a conversion from the 850 μm luminosities, show equivalent distributions, suggesting no differences in the potential impact of AGN feedback on these properties in red or control QSOs. These results suggest that there is no clear connection between the reddening properties of red QSOs and integrated properties at galaxy scales. We discuss these points in more detail in Sect. 4.

4. Discussion

We have studied the multiwavelength physical properties of red QSOs in comparison to a control QSO sample, carefully selected based on the observed g^* − i^* colours from the SDSS survey photometry. Following the selection approach presented by Klindt et al. (2019), we define our red and control QSO samples as the reddest 10%, and the central 50% of the colour distributions, respectively. Motivated by fundamental differences found in the radio properties among these samples reported by Klindt et al. (2019), Fawcett et al. (2020), Rosario et al. (2020), we undertook a comprehensive investigation of the FIR-to-UV emission to elucidate the nature of red QSOs. We utilised the Bayesian SED-fitting code AGNFITTER to disentangle multiple physical components (i.e. the AGN accretion disk and dusty torus; the host-galaxy stellar and cold dust emission) responsible for the overall panchromatic emission over the FIR–UV waveband. Overall, we found remarkable similarities in the SEDs between red and control QSOs which provide us with interesting insights on the nature of the reddening in red QSOs. We discuss different implications of these results in the next sections.

Combining our SED estimates and the spectroscopic properties from SDSS QSO spectra, in Sect. 3.2 we used the virial approximation to infer the distributions of black hole masses and Eddington ratios for our red and control QSO samples. We found remarkably similar distributions and no significant differences in these properties. Although this is in tension with a few results which relate red colours to different accretion properties (e.g. Richards et al. 2003; Kim & Im 2018) these literature results are based on small samples (e.g. few to handfuls of near-infrared-selected red QSO), and usually their comparison samples are subject to other selection strategies (e.g. optical selected) making direct comparisons challenging. On the other hand, the result that the reddening is linked to dust attenuation rather than intrinsic red QSO continua is in overall agreement with the majority of QSO investigations in the literature (e.g. Richards et al. 2003; Glikman et al. 2012, 2013; Kim & Im 2018) and consistent with our detailed fitting to X-shooter UV–near-IR spectroscopy for red and control QSOs selected in the same manner as this study (Fawcett et al., in prep.). However, although dust attenuation is the most straightforward explanation for the reddening in QSOs, the origin or distribution of the obscuring material is still unclear. We investigate this question in the following sections.

4.1. The origin of QSO reddening: No link with orientation and host galaxy ISM

In this section we use our SED-fitting results to search for connections between the dust attenuation in red QSOs and the source of the attenuation; i.e. whether this occurs at nuclear or host-galaxy scales. Firstly, we investigate whether the reddening in QSOs arise from specific configurations of the torus structure, which is in turn linked to viewing-angle and parametrised by a column density parameter in our dusty torus model. We find the dusty ‘torus’ components have almost identical shapes for both populations, with modelled column densities of log N_H ∼ 21.5 cm⁻² in both cases (Fig. 7) implying that no direct connection exists between the dust reddening in red QSOs and the emission from the dusty torus, responsible for the unification scheme. On the contrary, in the context of unification, the tori modelled for the majority of red QSOs are consistent with a face-on unobscured view of prototypical Type 1 QSOs. We believe systematics from our modelling technique are not affecting these results, given that this model has been shown to characterise inclination properties reliably, being able to robustly recover spectroscopic classifications of Type 1 and Type 2 AGN based solely on photometry (Calistro Rivera et al. 2016). This observation argues against the possible link between QSO reddening and inclination, and suggests the reddening must arise from other structures, potentially at different scales.

Indeed, these results are supported by the observations presented by Klindt et al. (2019), Rosario et al. (2020), and Fawcett et al. (2020), which find fundamental differences in the radio properties of red and control QSOs that are inconsistent with the scenario where the red QSOs are more inclined, on average, in comparison to the control QSOs. This finding is also in agreement with several recent studies of red QSOs that support an evolutionary scenario (Urrutia et al. 2012; Glikman et al. 2012; Banerji et al. 2012, 2017), although in tension with the interpretation by Wilkes et al. (2002), which used X-ray spectroscopy in red QSOs and interpret their values of log N_H(X-rays)∼21.5 as a signature of intermediate inclination between Type 1 and Type 2. We argue X-ray-based column densities are not necessarily linked to the torus structure, but might be connected to any gas structure on the line of sight, or different dust-to-gas ratios in the composition of the obscuring medium (Circosta et al. 2019; Lansbury et al. 2020).

We further studied a possible link between the reddening and the torus reprocessing efficiency (Fig. 8), and found the tori in red QSOs are at most equally efficient in reprocessing the intrinsic QSO emission, i.e. no additional dust-attenuation can occur in the torus structure. This reinforces our first conclusion that no clear link is found between the nature of the tori in red QSOs and their colour. In agreement with our results, Kim & Im (2018) present analysis of near-infrared spectra of 20 red QSOs at z < 0.9 finding no differences in their covering factors as compared to a larger samples of unobscured SDSS QSOs. In conclusion, this set of evidence puts forward the scenario where no clear connection between reddening in QSOs and inclination exists, suggesting QSO reddening may be linked to a different source of obscuration, possibly related with evolution.

Our SED-fiting results have also allowed us to explore the connection between QSO reddening and the stellar masses, star-formation rates and interstellar medium of their host galaxies (Fig. 11), finding no significant differences in these properties between the red and control QSO populations. In particular, we find that the molecular gas masses inferred from the rest-frame infrared 850 μm emission, together with total infrared luminosities (used to estimate SFR) are equivalent for red and control QSOs. We note that there are large uncertainties associated with the estimation of the rest-frame 850 μm emission, since these are extrapolated from the inferred SEDs given our lack of photometric data coverage at those wavelengths (our reddest band is SPIRE 500 μm). However, these uncertainties are well constrained by our Bayesian inference approach. We aim to investigate this point more in detail in a forthcoming publication based on ALMA 1 μm observations for a subset of red and control QSOs. Based on our current results on the inferred stellar population and interstellar medium properties, we argue that no clear signature is found that the reddening might be linked to host galaxy scales.

4.2. The origin of QSO reddening: the dusty wind scenario

To identify infrared features that are not fully recovered by the torus main structure, nor host galaxy emission, we searched for residuals from our SED-fitting of the red and control QSOs (Fig. 9). We find significant dust residuals evident as infrared bumps at 2 − 5 μm (consistent with hot dust at 500−1500 K), which are on average more present in red QSOs than in QSOs from the control sample. This trend is observed to be present across all redshift bins covered in our study, precluding any significant host galaxy contribution as the latter is negligible at higher redshifts and QSO luminosities. Additionally, we find a general correlation between the prominence of these residuals and the QSOs reddening parameter E(B − V)_BBB (Fig. 10). We interpret this observed infrared excess as an indicator of the nature of the obscuring medium, which causes the reddening in red QSOs.

Indeed, the temperatures indicated by these wavelengths can provide information on the scales over which the obscuring dust is distributed. Following the prescription of Barvainis (1987) (similar to Kawaguchi & Mori 2010; Stalevski et al. 2016), we assume a spatial distribution of dust as a flared disk, and we can estimate the radius where the dust needs to be located in order to produce this emission. We estimate this following:

$$\begin{array}{c}\hfill \left({\displaystyle \frac{R}{\mathrm{pc}}}\right)\simeq 1.3{\left({\displaystyle \frac{L_{\mathrm{bol}}}{{10}^{46}\mathrm{erg}{\mathrm{s}}^{-1}}}\right)}^{0.5}{\left({\displaystyle \frac{T_{\mathrm{dust}}}{1500\mathrm{K}}}\right)}^{-2.8}{\left({\displaystyle \frac{a}{0.05\mathrm{\mu }\mathrm{m}}}\right)}^{-0.5},\end{array}$$(3)

where we used the median bolometric luminosity values for the red QSO sample log L_bol = 46.5, T_dust is the black body temperature corresponding the peak of the residual peak, which lies approximately at 1000 K, and a is the dust grain size, for which we assume the Galactic ISM value of ∼0.05 μm. We estimate that the dust producing the infrared excess emission is located at a radius of R ∼ 7 pc from the central engine, which corresponds to a factor of 3−14 of the sublimation radius for average Galactic grains and large graphites, respectively. Although we note these estimates are based on a very approximative approach, these values are consistent with the inner region of the torus. Observationally, similar infrared excess or hot dust bumps around the same wavelengths have been previously found based on accurate infrared spectroscopy in nearby QSOs and Seyfert galaxies (Mor et al. 2009; Riffel et al. 2009; Burtscher et al. 2015), predominantly in unobscured, Type 1 AGN, local equivalents to those in our sample. In particular, Mor et al. (2009) analysed detailed infrared spectra of Palomar-Green (PG) QSOs finding that a hot dust black body component at temperatures 900−1800 K is required to fit the observed residuals. Despite not including optical data points, in contrast to our study, they are also able to reject any possible host galaxy contribution (see also Hernán-Caballero et al. 2016). A similar infrared excess has been also reported in luminous Type 1 QSOs at intermediate redshifts (z > 1), connecting it to emission arising from the sublimation radii of the torus (Temple et al. 2021).

While our infrared modelling describes only the torus emission, an additional emission component could be responsible for causing the reddening of the red QSOs. One scenario which is being increasingly discussed in the observational and theoretical literature is the presence of dusty polar winds in AGN (Mor et al. 2009; Hönig & Kishimoto 2017; DiPompeo et al. 2018; Baron & Netzer 2019; Stalevski et al. 2019; Venanzi et al. 2020; Harrison et al. 2018; Costa et al. 2018a). Indeed, the presence of winds at the scales estimated above, and the corresponding infrared emission at the wavelength where the infrared excess was found (2−5 μm) is predicted by theoretical work (e.g. Hönig 2019), as a consequence of AGN and infrared radiation pressure leading to a puffed-up disk/wind-launching region (e.g. Roth et al. 2012). The scenario that the obscuring component is distributed in such a wind seems plausible according to recent work which suggests that dust ‘forms preferentially’ in outflowing gas at these scales (Sarangi et al. 2019). A widely-used tracer of winds and outflows in AGN is the [O III]λ5007 emission line (from now on referred as [O III]), which traces low-density warm ionised gas in the narrow-line region of AGN. The presence of high velocities and strongly kinematically disturbed [O III]-emitting gas is a sign-post for outflows on scales of few pc to few kpc (e.g. Baskin & Laor 2005; Harrison et al. 2014; Woo et al. 2016; Zakamska et al. 2016). Since these scales are consistent with our estimation from the infrared residuals, we explore the potential presence of such winds and outflows in our data.

To examine the distribution of FWHMs of [O III] wings, we use the spectral fitting output from the SDSS spectral catalogue by Rakshit et al. (2020) (partially through priv. comm.). The [O III] line emission was covered by the SDSS spectral band only for sources at z < 1. A double Gaussian model was used to represent the [O III] emission line; one for the core and another for the wing components. During the fit, the FWHM of the core component was allowed to vary between 100 and 900 km s⁻¹ and the velocity offset was restricted to ±1000 km s⁻¹, while for the wing component, the FWHM was varied between 150 and 2100 km s⁻¹ and the velocity offset was restricted to ±3000 km s⁻¹. To avoid possible contamination from host galaxy features, we make a further cut for QSO luminosities with log L_6 μm > 44.5 ergs s⁻¹. We note that the luminosity cut and the redshift cut applied to cover the [O III] emission reduce our original sample by 87% (42 rQSOs and 31 cQSOs). In the left panel of Fig. 13 we plot the distributions of [O III] FWHM for the wing components of the red and control QSO samples for all sources with [O III] measurements and luminosities of log L_6 μm > 44.5 ergs s⁻¹.

Fig. 13. Cumulative distributions of the FWHMs (in the units of km s−1) of the wing component of the [O III]λ5007 emission line for the red and control QSOs. Two cuts were required from our parent sample due to observability of the [O III] line (z < 1) and to avoid contamination from the host galaxy (log L6 μm > 44.5). Left panel: entire population, without any further quality requirement. Right panel: we restrict our study to high quality detections of the wing components, based on a SNR criterion defined in the plot, which is consistent with our visual inspection of the lines. The high-quality detection cut increases the evidence for a prevalence of high velocity winds in red QSOs as compared to control QSOs.

We further test for reliability of the detection of the [O III] wing components both through visual inspection and by comparing the significance of the emission line peak compared to the continuum noise at 5100 Å, defining the SNR${}_{\mathrm{OIII}}=F_{\mathrm{peak}}^{[\mathrm{O}\mathrm{I}\mathrm{I}\mathrm{I}]}/{\mathrm{rms}}_{5100}$ (following Rakshit & Woo 2018). Through visual inspection, we find this criterion is a good representation of the reliability of the wing. After applying a reliability cut at SNR_OIII > 3 we remain with a sample of 23 sources (11 cQSOs, and 12 rQSOs). The right panel of Fig. 13 clearly shows the prevalence of the [O III] wings in red QSOs as compared to the control QSOs. In order to overcome the limitation of the small sample size and increase the statistical significance of our results, in Appendix A we explore the [O III] line emission properties for the red and control QSOs in the entire SDSS sky. Using samples of > 1000 QSO spectra we confirm the prevalence of higher [O III] wing velocities in the overall red QSO populations as compared to the control sample with high significance (p-value = 10⁻⁴). The most striking difference lies in the fractions of red and control QSOs which fulfil the requirement of having [O III] detections at all, or [O III] measurements of a given SNR_OIII. Indeed, the comparison of the distribution of SNR_OIII for red and control QSO, matched in redshift and luminosity, is highly significant as well (p-value = 10⁻²⁷), as shown in Fig. A.1.

To make a more direct connection between the [O III] high-velocity winds and the infrared residuals, in Fig. 14 we show the cumulative distributions of the infrared residuals for all our QSO samples, divided into two bins of [O III] line velocities in each panel. In the high-velocity bin (right panel) a clear predominance of infrared residuals can be recognised for red QSOs compared to control QSOs, whereas in the low-velocity bin (left panel) the difference is negligible. This result supports the scenario where high-velocity nuclear winds, which are more prevalent in red QSOs (Fig. 13), are the main transport of obscuring medium into the line of sight of the QSO, both producing the hot dust excess and the colours of red QSOs.

Fig. 14. Cumulative distributions of the infrared residuals defined in Fig. 9 for the red and control QSO samples, split into two bins of low-velocity (left panel) and high-velocity [O III] wings (right panel). A clear predominance of residuals can be recognised for red QSOs as compared to control QSOs in the right panel, whereas this difference disappears in the left panel. This result supports the scenario where high-velocity winds, which are more prevalent in red QSO, are the main transport of obscuring medium in the QSO line of sight, producing the hot dust reprocessed emission and the colours of red QSOs.

Although these results provide significant evidence for the prevalence of high velocity winds in red QSOs, our conclusions are limited to low redshifts (z < 1). To rule out a redshift dependence of this result, we further investigate the presence of high-velocity winds in our red and control QSOs at 1.4 < z < 2.5 by analysing the spectral properties of the C IV emission line, reported by Rakshit et al. (2020). In particular, the C IV line profile is commonly used as a tracer of outflows that are driven off the accretion disk, as it often presents asymmetric shapes associated with high-velocity components. Such components are able to shift the observed velocities from the rest-frame velocity of the line by up to 5000 km s⁻¹ (e.g. Richards et al. 2011; Coatman et al. 2019; Temple et al. 2021; Vietri et al. 2020). We calculate the velocity shifts for our matched sample of red and control QSOs as |Δv_C IV| = |c × (λ_rest − λ_peak)/λ_rest|, where λ_rest is the rest-frame wavelength of the C IV doublet at 1549.48 Å, and λ_peak represent the peak wavelength measured and reported by Rakshit et al. (2020). We remind the reader that the C IV emission line parameters in Rakshit et al. (2020) correspond to the entire C IV profile without subtracting a narrow component. This choice was made due to the ambiguity in the presence of narrow components in these lines (Shen et al. 2011, 2020). A detailed description of the spectral fitting can be found in Rakshit et al. (2020). In addition to the redshift limits imposed to cover the C IV line in SDSS (1.4 < z < 2.5), we apply a cut in luminosity (log L_6 μm > 44.5 ergs s⁻¹) for consistency with our [O III] study and require a C IV line equivalent width with C IV EW > 10 Å to ensure robust velocity shift estimates.

In Fig. 15, we show the cumulative |Δv| distributions for our luminosity-redshift matched sample. We find that there is a significant difference in the distributions of |Δv| of red and control QSOs, where red QSO show higher C IV velocity shifts. In line with the results presented by Coatman et al. (2019) and Vietri et al. (2020), where they find a correlation between the outflow signatures in the [O III] and C IV emission lines, we conclude that this result is consistent with our findings for red and control QSOs at low z. This finding demonstrates that the enhanced velocities in the broad and narrow line region in red QSOs is a property which holds for both the low and high-z population. Finally, in Fig. 16 we also investigate a direct relation between the high-velocity wind signature and the infrared residuals. We find that red QSOs with higher C IV velocity shifts exhibit stronger hot dust excess, while no differences are found for the control QSOs irrespective of C IV spectral properties. These results are also consistent with what we observed using the [O III] spectral properties, providing evidence of a connection between high-velocity dusty winds and QSO reddening across the entire redshift range probed (0.2 < z < 1., and 1.4 < z < 2.5).

Fig. 15. Cumulative distributions of the velocity shifts of the C IV emission line (in units of km s−1) for our samples of red and control QSOs, matched in redshift and luminosity. These samples cover the redshift range 1.4 < z < 2.5, complementing the results presented using the [O III] lines at low redshift (z < 1). A significant enhancement of C IV velocity shifts is found for red QSOs as compared to control QSOs, which is consistent with our finding in Fig. 13.

Fig. 16. Cumulative distributions of the infrared residuals defined in Fig. 9 for the red and control QSO samples, split into two bins of low (left panel) and high C IV velocity shifts (right panel). An enhancement of residuals can be seen for red QSOs with high C IV velocity shifts, a signature of potential winds or outflows. This difference is not seen for control QSOs, which show overall a low incidence of residuals. These results are consistent with our [O III] study at low redshifts (z < 1), and support the proposed scenario of a connection between QSO reddening and the presence of dusty winds.

This result is the first study presenting a potential systematic connection between QSO optical reddening and the presence of dusty outflows. However, it is in agreement with a large number of recent studies on luminous optical QSO samples (e.g. Mullaney et al. 2013; Bischetti et al. 2017), infrared (e.g. DiPompeo et al. 2018), as well as lower luminosity local AGN samples (e.g. Rojas et al. 2020). Additionally, our result provides an interesting connection to extremely red QSOs (ERQs), where very high [O III] velocities of FWHM > 2000 km s⁻¹ have been observed (Zakamska et al. 2016; Hamann et al. 2017; Mehdipour & Costantini 2018; Perrotta et al. 2019; Temple et al. 2019; Villar Martín et al. 2020). In particular, Zakamska et al. (2016) studied 4 ERQs at z = 2.5 and log L_6 μm > 47 ergs s⁻¹ finding luminous, broad (FWHM = 2600 − 5000 km s⁻¹) and blue-shifted lines, estimating ∼3% of the bolometric luminosity of the QSOs was converted to outflow kinetic luminosity extended to a few kpcs. In line with our study, they found [O III] kinematics positively correlated with infrared luminosities, although these four sources are extreme counterparts of our sample in terms of most properties presented here. A similar connection has been found by Temple et al. (2021), reporting enhanced rest-frame 2 μm emission correlated with the blueshift of the [C IV]λ1550 line for a sample selected in the optical and near-infrared bands. In agreement with our findings, they find that objects with stronger signatures of nuclear outflows tend to have a larger covering fraction of dust around the sublimation-temperature, although they find the sources with higher excess have bluer continua and argue that no connection exists to QSO reddening. While several differences exist between our sample selection as they do not focus on a colour selection, this result is intriguing and underlines the need for more detailed spectroscopic investigations both in the optical and mid-infrared regimes to confirm our interpretation.

4.3. On the role of red QSOs in galaxy evolution

The host galaxies of red QSOs have been explored for just a few objects, given the difficulty in disentangling the stellar emission from the dominant accretion disk component and the lack of homogeneous deep Herschel or ALMA imaging. Indeed, due to these challenges existing studies find diverging evidence regarding several questions around the host galaxies, such as whether the host galaxies of red QSOs are undergoing mergers in a similar manner (Zakamska et al. 2019) or more often than blue QSOs (Urrutia et al. 2008, 2012), and whether red QSOs host galaxies are more strongly star-forming than typical QSOs (Georgakakis et al. 2009) or whether there are no significant differences (Urrutia et al. 2012; Wethers et al. 2020). Two main challenges in the task of comparing these different studies is that they often involve low statistics and inhomogeneous selection effects.

In our work we exploit the availability of rich multiwavelength UV-to-FIR ancillary data and the capability of AGNFITTER to robustly infer underlying uncertainties, in order to model the host galaxy properties of our controlled samples of red and control QSOs. Applying our careful comparative exercise on a statistical sample, we find no significant differences in the stellar mass, star-formation rates, and molecular gas estimates of red and control QSOs (Fig. 11). Taking into account the uncertainties associated with our measurements, we find both populations are overall consistent with the massive end of the main sequence of star forming galaxies, in line with previous studies of the star formation in general samples of optical QSOs (Rosario et al. 2013; Stanley et al. 2017; Wethers et al. 2018). Interestingly, we infer no significant differences in the molecular gas masses of the red and control QSOs.

Since our SED-based results argue against a connection between the reddening in QSOs and inclination of the dusty torus, we further evaluate our observations in the context of the galaxy evolution models (Sanders et al. 1988; Hopkins et al. 2008; Narayanan et al. 2010) that propose red QSOs as intermediate stages between dusty star forming galaxies, and UV-bright unobscured QSOs. Based on the similarities found in the host-galaxy properties of red and control QSOs, we conclude that while an evolutionary link might exist between red and control QSOs, the impact it has on the host galaxies is not large enough for it to be apparent in our statistical study. This could occur if the timescales separating these phases are not long enough to produce such an effect. To put this into context, recent simulations by Costa et al. (2020) indicate small-scale AGN winds would require timescales of ∼100 Myr−1 Gyr to significantly reduce star formation in their host galaxies, through the removal of high-density gas or suppression of halo gas accretion when coupled to large-scale outflows. By comparison, the typical estimated lifetime for a QSO is ∼10 Myr (e.g. Martini & Weinberg 2001). Assuming all QSOs undergo a red QSO phase then the lifetime for this phase on the basis of our sample would be 1 Myr (i.e. red QSOs comprise 10% of our overall QSO sample). We therefore conclude that the QSO lifetime is too brief to be able to expect significant differences in the host galaxies between red and control QSOs. We note though that given our uncertainties, subtle changes in these properties might be still present.

In the context of the proposed evolutionary sequence between red and control QSO, we speculate that the presence of these winds could be a key ingredient for the evolution. In such a scenario, these small-scale dusty winds would redden the QSO continuum, while expelling a fraction of the gas and dust content from the central regions, becoming progressively less energetic with less available material. This would be in line with models (e.g. Roth et al. 2012) which suggest that powerful small-scale winds would accelerate through radiation pressure, which within a given timescale, would be efficient in unobscuring the AGN, with limited impact on the host galaxy (e.g. Costa et al. 2018a,b). The evidence that the infrared residuals correlate on average with the E(B − V) parameter (estimated independently from the optical SED), and the fact that a fraction of control QSOs also present high-velocity winds, is consistent with the scenario of a potential transition between these populations.

5. Conclusions and summary

Motivated by recent evidence of fundamental differences in the radio properties of red and control QSO, we presented a multiwavelength comparative analysis of ∼1800 red and control SDSS QSOs, where red were defined as the 10% reddest, and control as the central 50% of the population in the (g^* − i^*) colour distribution. Using the Bayesian MCMC-based AGNFITTER code, we performed a detailed AGN-tailored SED-fitting analysis on a rich compilation of multiwalength photometry, covering up to 29 photometric bands from the UV to the far-infrared. The SED fitting results together with spectroscopic information extracted from the QSO spectra led us to the following conclusions:

• (i)

  AGNFITTER models the dust attenuation of the accretion disk emission using an SMC reddening law, finding distributions of reddening of $E{(B-V)}_{\mathrm{BBB}}=0.{02}_{-0.01}^{+0.04}$ mag ($A_V=0.{05}_{-0.03}^{+0.11}$ mag) for our control sample of ‘average’ optically-selected QSOs. Our red QSO selection has median reddening values of $E{(B-V)}_{\mathrm{BBB}}=0.{12}_{-0.08}^{+0.21}$ mag ($A_V=0.{32}_{-0.21}^{+0.57}$ mag), demonstrating that we are probing mild attenuation in the most reddened optically-selected QSO. Already with this mild attenuation, fundamental differences have been previously found in the radio properties of these sources (Klindt et al. 2019; Rosario et al. 2020; Fawcett et al. 2020).

• (ii)

  Apart from the diversity in attenuation of the accretion disk emission, we find no strong variations in the multiwavelength SEDs of red and control QSOs across all optical-to-FIR wavelengths. Based on spectral properties and reddening-corrected continuum measurements, we evaluate intrinsic properties of our QSO samples, such as black hole masses and Eddington ratios, finding no significant differences in the accretion properties of the two populations.

• (iii)

  We find the bulk of the infrared emission of both red and control QSOs are well described with the same torus emission structure, associated with low column densities of log N_H < 22 and face-on viewing angle. Additionally, we compare the torus reprocessing efficiency (L_tor/L_bol) and covering factors, finding that the red QSO population has covering factors equivalent to those of control QSOs. The similarity of the infrared emission of the torus component suggests that there is no connection between the QSO reddening and the emission expected from the AGN torus, disproving unification models as the origin of the red QSO nature.

• (iv)

  We investigate the presence of residual emission from infrared components other than the torus. Interestingly, we find a higher incidence of residual flux in the rest-frame 2 − 5 μm SEDs of red QSO, associated with temperatures of 500 − 1700 K and locations of ∼5 − 10 pc from the nucleus. We suggest the origin of this emission lies in nuclear high-velocity winds, for which we find tentative evidence through high-velocity ionised winds (traced by broad [O III] line wings of FWHM > 1000 km s⁻¹ at low-z, and C IV velocity shifts at high-z). In particular, we find a higher incidence of high velocity winds for red QSOs as compared to control QSOs, and more prominently when residual fluxes are present.

• (v)

  We measure the host galaxy properties of red and control QSOs. We find no significant differences in the stellar masses or star-formation rates of red and control QSOs. Interestingly, we find no significant differences in the masses of the molecular gas reservoirs of the two populations, as inferred from extrapolations of the cold dust SEDs. We conclude no clear link can be recognised between the host galaxy properties and the reddening of QSOs, suggesting that the driver of the differences between red and control QSOs must be on nuclear scales.

Based on our statistical comparative study, we find tentative evidence for a connection between the reddening in red QSOs and the presence of dusty winds in the QSO nuclear region. Within the QSO evolutionary framework, high-velocity winds can indeed play an important role in ejecting the material away from the nuclei, consistent with observations that red QSOs are usually more abundant at higher bolometric luminosities. In a forthcoming paper we connect our UV-to-FIR findings with deep low-frequency radio observations to expand our understanding of the nature of red QSOs and the origin of their enhanced radio detections. More detailed observations of the infrared spectra of a controlled QSO sample are required in order to confirm the prevalence of the hot dust emission component around ∼2 − 5 μm in red QSOs and its link to dusty winds, as the origin of QSO reddening.

Acknowledgments

We would like to thank the anonymous referee for their insightful comments and suggestions which have considerably improved the manuscript. We would like to thank Dominika Wylezalek and George Lansbury for discussions that contributed to the analysis and interpretation of our results. G.C.R. acknowledges the ESO Fellowship Program and a Gruber Foundation Fellowship grant sponsored by the Gruber Foundation and the International Astronomical Union. D.M.A and D.J.R. gratefully acknowledge support from the Science Technology Facilities Council (STFC) through grant ST/000244/1. M.S. acknowledges support by the Ministry of Education, Science and Technological Development of the Republic of Serbia through the contract no. 451-03-68/2020-14/200002 and the Science Fund of the Republic of Serbia, PROMIS 6060916, BOWIE. V.A.F. acknowledges a quota studentship funded by the STFC. P.N.B. is grateful for support from the UK STFC via grant ST/R000972/1. M.B. acknowledges support from INAF under PRIN SKA/CTA FORECaST and from the Ministero degli Affari Esteri della Cooperazione Internazionale – Direzione Generale per la Promozione del Sistema Paese Progetto di Grande Rilevanza ZA18GR02. R.K. acknowledges support from the STFC through an STFC studentship via grant ST/R504737/1. The analysis in this work made use of AstroPy (Astropy Collaboration 2013), Matplotlib (Hunter 2007), SciPy (Virtanen et al. 2020) and Pandas (Wes 2010). Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the US Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS website is http://www.sdss.org/. The SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions. The Participating Institutions are the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, University of Cambridge, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington. 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. The Herschel Extragalactic Legacy Project, (HELP), is a European Commission Research Executive Agency funded project under the SP1-Cooperation, Collaborative project, Small or medium-scale focused research project, FP7-SPACE-2013-1 scheme, Grant Agreement Number 607254.

References

Appendix A: [O III] wing properties of red QSOs and control QSOs in the SDSS survey

To complement the analysis of the [O III] wing properties of our red and control QSO sample we reapply our sample selection to the entire sky coverage of the SDSS survey. In this way we also overcome the limitation of the small sample statistics discussed in Sect. 4 after the consideration of signal-to-noise ratio cuts for the [O III] measurements. After applying our selection (WISE detected sources, classified based on their g^* − i^* colours) and a redshift cut to z < 1 to ensure the [O III] is covered by the SDSS spectra, overall we find a total of 6753 red QSOs and 6753 control QSOs, matched in redshift and luminosity. We then compare the signal-to-noise ratio of the [O III] measurements of this sample as shown in the lower panel of Fig. A.1, finding a significant difference in the signal to noise ratio of the [O III] emission in red QSOs. Interestingly, a larger fraction of control QSOs was discarded since these did not fulfil the requirement of having [O III] detections. We also investigate the difference in the FWHM of the [O III] wings, finding a significant enhancement of high-velocity ionised outflows in red QSOs, confirming the measurements discussed in Sect. 4 for a subset of these.

Fig. A.1. [O III] emission line properties for the luminosity-and-redshift-matched red and control QSO samples from the entire SDSS survey. Upper panel: comparison of all red and control QSOs with [O III] line measurements with Fpeak > 10 × rms5100 Å. It is important to note the number of red QSOs satisfying the condition is ∼50% larger than the control QSOs. This is shown more clearly in the lower panel, where the distribution of signal-to-noise ratio of the wing emission-line component, defined as = Fpeak(wing)/rms5100 Å, is shown for red and control QSOs. While here the only condition for a source to be included was to have [O III] detections, the number of red QSOs satisfying this condition is 5 times larger than that of the control sample.

All Figures

Fig. 1. Diagram to summarize our selection process. The shaded area corresponds to the general definition of red and blue QSO samples from the overall SDSS DR14 (Pâris et al. 2018), following the strategy previously presented by Fawcett et al. (2020), Rosario et al. (2020). The purple shaded area correspond to the field selection, determined by the availability of deep multiwavelength catalogues provided by the LoTSS data release in three northern fields and the public compilation of the HELP collaboration in the three equatorial fields of the GAMA survey.

In the text

Fig. 2. Monochromatic 6 μm luminosities (L6 μm) of red QSOs (red filled circles), control QSOs (sky-blue filled circles) and Lz-matched control QSOs (dark-blue circles) as a function of redshift are plotted in the central panel. Side panels show normalised histograms of the L6 μm and redshift distributions for each population, shown as solid lines (red QSO), dashed lines (L − z matched samples) and shaded areas (entire control sample) of respective colours. By comparing the distributions of the red QSOs and the Lz-matched control QSO sample, it can be clearly seen that equivalent QSO properties are being sampled.

In the text

Fig. 3. Examples of the SED decomposition with AGNfitter. In the main panels the photometric data points are shown with black error bars. The best fit and ten realisations of the posterior distributions of the fit are drawn as coloured solid, and transparent lines, respectively. The realisations display the uncertainties associated with the fits. The different colours represent the total SED (red lines), the accretion disk emission (blue lines), the host galaxy emission (yellow lines), the torus emission (purple lines) and the galactic cold dust emission (green lines). Lower panels: residuals expressed in terms of significance given the data noise level for the best fit and the 10 realisations. Through this selection of SED fits we exemplify the diversity of the obscuration levels found. Two examples of red QSO SEDs which clearly present an infrared excess emission at ∼3 μm, as discussed in Sect. 3.3.2, are included in the two lower left panels.

In the text

Fig. 4. Ten modelled SEDs are plotted for each source (reconstructed from 10 randomly-picked realisations of the posterior PDFs) as transparent red (red QSOs) and blue lines (control QSOs) to display the uncertainties. Solid red and dotted blue lines depict the composite median SED for the red and control QSO population respectively, matched in redshift and L6 μm. At frequencies log ν < 14.5 Hz, strikingly similar composite SEDs can be observed for the red and control QSO populations. At frequencies log ν > 14.5 Hz a clear difference in the composite SEDs can be observed which is the signature of dust reddening consistent with the optical colour selection.

In the text

Fig. 5. Median composite SEDs are shown for the red QSO sample (solid lines of different colours) and matched control QSO samples (dotted lines of different colours), now split into different bins of L6 μm (left-right) and z (top-bottom). All composite SEDs have been renormalised to the same L6 μm for purposes of visualisation. The grey solid and dotted lines represent the total SEDs for the red and control samples, respectively. The remaining colours codify the different physical components as described in the legend of Fig. 3. The luminosity and redshift bins considered and the number of sources included in each bin are annotated. The contribution of the stellar emission of the host galaxies has a clear contribution to the total SEDs in the lower redshift bins and to luminosities of L6 μm < 45. A remarkable similarity in the composite torus SEDs of red and control QSOs across all redshift bins can be seen.

In the text

Fig. 6. Distributions of black hole mass estimates (MBH), reddening-corrected bolometric luminosities Lbol and Eddington ratios λEdd for the control and red QSO samples. Only MBH estimates from spectra with quality tag 0 are considered, which comprise 90% of our entire sample. A Gaussian kernel-density estimator (KDE) with a bandwidth of 0.3 is used to display these distributions. With exception of the differences in the bolometric luminosity distributions, where significance is weak, no significant differences are found between the red and the control QSOs in these parameters. The significance of the difference between the two distributions (KS-test p-value) is reported below. The p-values were computed using the estimates reported by Rakshit et al. (2020) (MBH and λEdd), and using the bootstrapped KS-test for the Lbol, where the posterior distributions from the MCMC-based SED fitting are compared and therefore uncertainties could be calculated.

In the text

Fig. 7. Probability distributions of the accretion disk reddening parameter E(B − V), and the column density NH. The E(B − V) parameter was inferred assuming the SMC attenuation law by Prevot et al. (1984) on the accretion disk emission, while the NH is a quantity that parametrizes the (X-ray-normalised) infrared torus templates by Silva et al. (2004). These distributions are constructed as the superposition of 100 random draws from each source. A gaussian kernel-density estimator (KDE) with a bandwidth of 0.1 is used to display these distributions.

In the text

Fig. 8. Composite two-dimensional posterior distributions of the torus reprocessing efficiency, defined as the ratio of torus luminosities and bolometric luminosities, plotted as a function of bolometric luminosities. These are constructed by sampling 100 realisations from the posterior of each source to properly account for the uncertainties. The red contour lines (for red QSOs) and the blue contour areas (for control QSOs) range between 20 and 100%, in intervals of 20%. Upper panels: this ratio for the observed bolometric luminosities (left panel) and for the reddening-corrected (intrinsic) bolometric luminosities (right panel) for red and control QSOs. Lower panel: torus reprocessing efficiency corrected for the torus anisotropy following Stalevski et al. (2016), to provide a better description of the covering factor.

In the text

Fig. 9. Residuals from the best fits for each source as a function of rest-frame frequency. Two small ‘bumps’ are identified in the MIR regime, where the torus models included in our model show the largest discrepancy with the data. The main discrepancies have positive residuals meaning the model cannot reproduce a part of the flux at these frequencies. Right panels: cumulative distributions of the incidence of residuals in each sample. We investigate whether these residuals are connected to the optical colour of the QSO, despite the fact these bands are independent in the fitting. We find that red QSOs have larger MIR emission excess compared to the model in the regime of 13.75 < log ν < 14.2 (2−5 μm), showing significant residuals ( > 2σ) in 20% of the cases, in contrast to the 10% in the case of control QSOs.

In the text

Fig. 10. Median infrared residuals are shown as a function of the QSO reddening parameter, E(B − V)BBB. Black circles represent the median residual values for each E(B − V)BBB bin for all QSOs irrespective of colour classification. The error bars represent bootstrap 1-σ errors. The same relation is shown for the red and control QSOs as red and blue shaded areas, respectively. A clear trend can be observed, where the infrared residuals increase with larger reddening values for all samples.

In the text

Fig. 11. Composite posterior distributions of host galaxy parameters inferred through SED-fitting. These are constructed by sampling 100 realisations from the posterior distributions of each source, and we apply a KDE with a bandwidth for 0.3 for display as explained in detail in Sect. 2. Stellar masses are estimated based on the modelling of near-infrared and optical emission with stellar population models (Bruzual & Charlot 2003). Star formation rates are computed based on total infrared luminosities (LIR(8 − 1000 μm)) from the modelling of the infrared photometry with cold dust models (Schreiber et al. 2018). Molecular gas masses are estimated from monochromatic luminosities of the dust continuum SED at 850 μm and the empirical relation by Scoville et al. (2016). No significant differences are found between the distributions of red and control QSOs for these host-galaxy parameters.

In the text

Fig. 12. Composite two-dimensional posterior distributions of star-formation rates versus stellar mass for red QSOs (red shaded contour areas) and control QSOs (blue contour lines) for the low-z sources. These are constructed by sampling 100 realization from the posterior of each source to properly account for uncertainties. Stellar masses and star-formation rates are estimated based on SED modelling. The distributions are consistent with the main sequence (grey shaded region), although as well with ULIRGS (GAMA survey; Driver et al. 2018) within the uncertainties.

In the text

Fig. 13. Cumulative distributions of the FWHMs (in the units of km s−1) of the wing component of the [O III]λ5007 emission line for the red and control QSOs. Two cuts were required from our parent sample due to observability of the [O III] line (z < 1) and to avoid contamination from the host galaxy (log L6 μm > 44.5). Left panel: entire population, without any further quality requirement. Right panel: we restrict our study to high quality detections of the wing components, based on a SNR criterion defined in the plot, which is consistent with our visual inspection of the lines. The high-quality detection cut increases the evidence for a prevalence of high velocity winds in red QSOs as compared to control QSOs.

In the text

Fig. 14. Cumulative distributions of the infrared residuals defined in Fig. 9 for the red and control QSO samples, split into two bins of low-velocity (left panel) and high-velocity [O III] wings (right panel). A clear predominance of residuals can be recognised for red QSOs as compared to control QSOs in the right panel, whereas this difference disappears in the left panel. This result supports the scenario where high-velocity winds, which are more prevalent in red QSO, are the main transport of obscuring medium in the QSO line of sight, producing the hot dust reprocessed emission and the colours of red QSOs.

In the text

Fig. 15. Cumulative distributions of the velocity shifts of the C IV emission line (in units of km s−1) for our samples of red and control QSOs, matched in redshift and luminosity. These samples cover the redshift range 1.4 < z < 2.5, complementing the results presented using the [O III] lines at low redshift (z < 1). A significant enhancement of C IV velocity shifts is found for red QSOs as compared to control QSOs, which is consistent with our finding in Fig. 13.

In the text

Fig. 16. Cumulative distributions of the infrared residuals defined in Fig. 9 for the red and control QSO samples, split into two bins of low (left panel) and high C IV velocity shifts (right panel). An enhancement of residuals can be seen for red QSOs with high C IV velocity shifts, a signature of potential winds or outflows. This difference is not seen for control QSOs, which show overall a low incidence of residuals. These results are consistent with our [O III] study at low redshifts (z < 1), and support the proposed scenario of a connection between QSO reddening and the presence of dusty winds.

In the text

Fig. A.1. [O III] emission line properties for the luminosity-and-redshift-matched red and control QSO samples from the entire SDSS survey. Upper panel: comparison of all red and control QSOs with [O III] line measurements with Fpeak > 10 × rms5100 Å. It is important to note the number of red QSOs satisfying the condition is ∼50% larger than the control QSOs. This is shown more clearly in the lower panel, where the distribution of signal-to-noise ratio of the wing emission-line component, defined as = Fpeak(wing)/rms5100 Å, is shown for red and control QSOs. While here the only condition for a source to be included was to have [O III] detections, the number of red QSOs satisfying this condition is 5 times larger than that of the control sample.

In the text