# Broadband spectral energy distributions of SDSS-selected quasars and of their host galaxies: intense activity at the onset of AGN feedback

Federico Bianchini, Giulio Fabbian, Andrea Lapi, Joaquin Gonzalez-Nuevo, Roberto Gilli, Carlo Baccigalupi

<figure><img src="image_1.png" /></figure>

► To cite this version:

Federico Bianchini, Giulio Fabbian, Andrea Lapi, Joaquin Gonzalez-Nuevo, Roberto Gilli, et al. Broadband spectral energy distributions of SDSS-selected quasars and of their host galaxies: intense activity at the onset of AGN feedback. The Astrophysical Journal, 2019, 871 (2), pp.136. 10.3847/1538-4357/aaf86b. hal-01833720

HAL Id: hal-01833720

https://hal.science/hal-01833720v1

Submitted on 13 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Broadband Spectral Energy Distributions of SDSS-selected Quasars and of Their Host Galaxies: Intense Activity at the Onset of AGN Feedback

Federico Bianchini¹, Giulio Fabbian², Andrea Lapi³,⁴,⁵, Joaquin Gonzalez-Nuevo⁶, Roberto Gilli⁷, and Carlo Baccigalupi³,⁴

¹School of Physics, University of Melbourne, Parkville, VIC 3010, Australia; fbianchini@unimelb.edu.au

²Institut d'Astrophysique Spatiale, CNRS (UMR 8617), Univ. Paris-Sud, Université Paris-Saclay, Bât. 121, F-91405 Orsay, France

³SISSA, Via Bonomea 265, I-34136 Trieste, Italy

⁴INFN-Sezione di Trieste, via Valerio 2, I-34127 Trieste, Italy

⁵INAF-Osservatorio Astronomico di Trieste, via Tiepolo 11, I-34131 Trieste, Italy

⁶Departamento de Física, Universidad de Oviedo, C. Federico García Lorca 18, E-33007 Oviedo, Spain

⁷INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Gobetti I-93/3-40129 Bologna, Italy
Received 2018 June 18; revised 2018 December 10; accepted 2018 December 11; published 2019 January 28

## Abstract

We present the mean spectral energy distribution (SED) of a sample of optically selected quasars (QSOs) at redshifts of $1 \le z \le 5$. To derive it, we exploit photometric information from the Sloan Digital Sky Survey, UKIRT Infrared Deep Sky Survey, and Wide-field Infrared Survey Explorer surveys in combination with a stacking analysis of Herschel, AKARI, and Planck maps at the location of the QSOs. The near-UV and optical parts of the reconstructed mean rest-frame SED are similar to those found in other studies. However, the SED shows an excess at $1-2 \ \mu\text{m}$ (when compared to the aforementioned SEDs normalized in the near UV) and a prominent bump around $4-6 \ \mu\text{m}$, followed by a decrease out to $\sim 20 \ \mu\text{m}$ and a subsequent far-IR bump. From the fitted SEDs, we estimate the average active galactic nucleus (AGN) luminosity $L_{\text{AGN}}$ and star formation rate (SFR) as a function of cosmic time, finding typical $L_{\text{AGN}} \sim 10^{46}-10^{47} \ \text{erg s}^{-1}$ and SFR $\sim 50-1000 \ \text{M}_\odot \text{yr}^{-1}$. We develop mid-IR-based criteria to split the QSO sample, finding that these allow us to move along the average relationship in the SFR versus $L_{\text{AGN}}$ diagram toward increasing AGN luminosities. When interpreted in the context of the in situ coevolution scenario presented by Lapi et al., our results suggest that the detection in the far-IR band is an effective criterion to select objects where the star formation is on the verge of being affected by energy/momentum feedback from the central AGN.

**Key words:** galaxies: active – infrared: galaxies – quasars: general

## 1. Introduction

According to our current understanding of galaxy formation and evolution, the nuclei of active galaxies (AGNs) are powered by accretion of matter onto supermassive black holes (SMBHs; Rees 1984), which makes them among the most luminous extragalactic sources (Soltan 1982). The observed correlations between the black hole mass $M_{\text{BH}}$ and the properties of the host galaxy ($M_{\text{BH}}-\sigma$ relation; Ferrarese & Merritt 2000; Gebhardt et al. 2000; Kormendy & Richstone 1995; Magorrian et al. 1998, $M_{\text{BH}}-M_{\text{bulge}}$) suggest the existence of a coupling between black hole (BH) growth and star formation (e.g., Lapi et al. 2014). A further observation that corroborates this picture is that both the cosmic star formation and the BH accretion histories peak around $z \sim 2$ (Aird et al. 2010; Madau & Dickinson 2014).

In the standard evolutionary framework, the host galaxy and the SMBH grow coevally. The galaxy undergoes phases of intense starburst activity while the central SMBH grows by accretion, until the AGN-driven large-scale outflows expel the gas and dust from the interstellar medium (e.g., King 2003; Di Matteo et al. 2005). The depletion of gas regulates the accretion of matter onto the SMBH and may eventually lead to the suppression of star formation activity, leaving a red and dead galaxy behind (e.g., Granato et al. 2004; Croton et al. 2006).

The specific interplay between AGN feedback and the surrounding material plays an important role in the evolution of the host galaxy, as it is thought to inhibit star formation and reheat the interstellar and circumgalactic medium (e.g.,

Sturm et al. 2011; Fabian 2012; Cicone et al. 2014). AGN feedback also leaves an imprint on the large-scale structure of the universe, as it affects the small-scale clustering of dark matter (van Daalen et al. 2011), potentially biasing the cosmological analysis of clustering and weak lensing data (Semboloni et al. 2011; Eifler et al. 2015). Despite its key role in galaxy evolution, the details of the mechanism and energetics behind AGN feedback are still poorly understood (see Alexander & Hickox 2012; Fabian 2012 for reviews).

The AGN spectral energy distribution (SED), which spans the entire electromagnetic spectrum from the radio band to the X-rays, provides a key tool to unravel the nature of the physical processes that govern the SBMH/host-galaxy coevolution. For example, based on the simple schematic provided by the unified AGN model (Antonucci 1993; Urry & Padovani 1995), the dust produced during the star-forming phase is supposed to form a torus surrounding the accretion disk that reprocesses most of the ultraviolet (UV) and optical photons into the infrared (IR) band. Since the AGN can contribute to the IR luminosity due to the thermal radiation re-emitted by the surrounding torus, it is clear that disentangling the contributions to the total SED from the accretion disk, the dusty torus, and the host galaxy is of pivotal relevance (e.g., Mullaney et al. 2011; Del Moro et al. 2013; Leipski et al. 2013; Delvecchio et al. 2014; Netzer et al. 2016; Symeonidis et al. 2016; Stanley et al. 2017).

The aim of this paper is to provide a multiwavelength picture of AGNs by reconstructing the mean SED of quasars (QSOs)
from the near-UV/optical to the far-IR. To this end, we use the optically selected QSOs from the Sloan Digital Sky Survey (SDSS). By means of template fitting, we derive constraints on the AGN bolometric luminosity and star formation rate (SFR) in the host galaxy as a function of redshift. We also suggest mid-IR-based selection criteria that can probe different parts of the SFR versus AGN luminosity diagram.

The SFR is usually estimated through a number of indicators such as UV, mid-IR, and total IR luminosities, as well as emission lines (e.g., Hα; e.g., Kennicutt 1998). Each of these proxies has its own drawbacks. For example, optical emission lines usually provide reliable SFR estimates (since they probe star formation on timescales shorter than the integrated UV or IR luminosities) but are much more time consuming to obtain, especially for high-z sources. On the other hand, UV and mid-IR luminosities can be heavily contaminated by the intrinsic AGN emission and can return biased estimates of the SFR if the AGN emission has not been properly removed.

In our study, we adopt the integrated far-IR luminosity (which mostly traces the cold dust and provides a less biased SFR measure) as our SFR indicator. This task is best accomplished by exploiting the high-sensitivity far-IR observations from the Herschel satellite (Pilbratt et al. 2010), which has opened a new window for measuring the rest-frame far-IR properties of QSOs at high z (z ≳ 1). Luminous optical QSOs have been extensively studied in the far-IR band, and the broadly accepted view is that they tend to live in galaxies with ongoing star formation, at rates similar to those of the star-forming population (e.g., Hatziminaoglou et al. 2010; Leipski et al. 2013; Rosario et al. 2013; Xu et al. 2015; Netzer et al. 2016; Stanley et al. 2017; Hall et al. 2018). Moreover, several studies argue for a positive correlation between the mean SFRs as a function of the bolometric AGN luminosity (e.g., Bonfield et al. 2011; Rosario et al. 2013; Kalfountzou et al. 2014; Gürkan et al. 2015; Harris et al. 2016; Podigachoski et al. 2016; Stanley et al. 2017). In addition to the far-IR Herschel observations, we also include mid-IR data from the Wide-Field Infrared Survey (WISE) satellite. As we will show, combining these multiwavelength observations can place strong constraints on the SFR of the host galaxy as well as on the AGN power.

This paper is structured as follows. In Section 2, we present the different data sets utilized in the analysis, while the stacking and SED-fitting methodology is discussed in Section 3. We present our results, compare them with the existing literature, and provide an interpretation in terms of galaxy/SMBH coevolution in Section 4. Finally, we draw our conclusions in Section 5.

Throughout this work, we adopt the standard flat cosmology (Planck Collaboration et al. 2018) with rounded parameter values: matter density Ω<sub>M</sub> = 0.32, baryon density Ω<sub>b</sub> = 0.05, Hubble constant H<sub>0</sub> = 100 km s<sup>-1</sup> Mpc<sup>-1</sup> with h = 0.67, and mass variance σ<sub>8</sub> = 0.81 on a scale of 8 h<sup>-1</sup> Mpc. The reported stellar masses and SFRs (or luminosities) of galaxies refer to the Chabrier (2003) IMF.

## 2. Data Sets

The main focus of this paper is the reconstruction of the mean SED of optical-selected QSOs with the aim of constraining the SFR and the AGN luminosity across cosmic time. To this end, we need access to data sets spanning the electromagnetic spectrum from the optical–UV to the

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Unit-normalized redshift distribution of the individual SDSS DR7 (orange) and DR12 (blue) QSO catalogs, as well as the merged sample (green) used in this work (after the radio-loud cut has been applied). The vertical dashed lines represent the redshift bins employed in our analysis. The smooth solid lines represent the kernel density estimates of the underlying redshift distributions.</figcaption></figure>

millimeter band. In this section, we first present the public data sets employed in the analysis and then describe the extracted QSO samples.

### 2.1. Optical SDSS QSO Catalog

We select our QSO sample from the publicly available SDSS-II and SDSS-III Baryon Oscillation Spectroscopic Survey catalogs of spectroscopically confirmed QSOs detected over 9376 deg². Specifically, we create a merged sample by combining the QSO catalogs from the seventh (DR7; Schneider et al. 2010)<fnref n="8" /> and twelfth SDSS data releases (DR12; Pâris et al. 2017).<fnref n="9" /> As can be seen in Figure 1, the DR7 sample mostly contains objects at z ≲ 2.5 while DR12 mainly targets sources at z ≳ 2.15, although a color degeneracy in the target selection of photometric data causes a fraction of low-redshift QSOs to be also observed, resulting in a secondary maximum around z ~ 0.8 (see Ross et al. 2012 for a detailed discussion of the QSO target selection process). Approximately 4% of the DR7 objects at z > 2.15 have been reobserved for DR12: any QSO present in both DR7 and DR12 is included only once.<fnref n="10" /> SDSS QSO catalogs contain several redshift estimates for the sources: in this work, we use the SDSS pipeline redshift as our baseline estimate. The sources in the merged catalog span a broad range of redshifts, out to z ~ 7. We cut the merged catalog to lie in the redshift range 1 ≲ z ≲ 5, discarding QSOs at z < 1 and at z > 5 because they provide a sparse sampling of the bright QSO population. The resulting catalog contains 313,004 sources (or 386,245 if no redshift cuts are applied). As discussed in greater detail in Section 2.5, the main QSO sample, as well as the other subsamples considered in this analysis, is extracted from this merged parent QSO catalog. In order to test the evolution of the SED over cosmic time, we further slice the QSO sample into three redshift bins: (i)

<footnote n="8">⁸ Available at http://classic.sdss.org/dr7/products/value_additional/qsocat_dr7.html.</footnote>

<footnote n="9">⁹ Available at http://www.sdss.org/dr12/algorithms/boss-dr12-quasar-catalog/.</footnote>

<footnote n="10">¹⁰ We retain the redundant DR12 information because of the extensive multiwavelength information available for this catalog.</footnote>
1 ≦ z < 2.15, (ii) 2.15 ≦ z < 2.5, and (iii) 2.5 ≦ z ≦ 5, such that they approximately contain the same number of sources.

In this work, we make use of the SDSS photometric information; specifically, we use the Galactic-extinction-corrected point-spread function (PSF) magnitudes in the *u*, *g*, *r*, *i*, and *z* bands, which are then converted to fluxes. The DR12 QSO catalog also provides multiwavelength matching with external surveys. Namely, we use cross-matching with (i) the latest Faint Images of the Radio Sky at Twenty-Centimeters (FIRST; Becker et al. 1995) radio survey at 1.4 GHz that has almost full overlap with SDSS with a sensitivity of about 1 mJy, (ii) the mid-IR data from *WISE*, specifically to the newly released AllWISE Source Catalog (Cutri et al. 2013), and (iii) to the sources detected in the near-IR images from the UKIRT Infrared Deep Sky Survey (UKIDSS; Lawrence et al. 2007). The radio-band information from FIRST is exploited to enforce a radio-loud cut in order to minimize the synchrotron emission from the QSOs.

## 2.2. Near-IR and Mid-IR Photometry

The mid-IR properties of our sample have been defined via cross-matching with the AllWISE sky survey as provided in the SDSS catalogs. *WISE* has scanned the whole sky in four bands centered at wavelengths of 3.4, 4.6, 12, and 22 µm (*W1*, *W2*, *W3* and *W4* channels). The *WISE* All-Sky Source catalog provides PSF-fitting magnitudes and uncertainties in the Vega system for the four mid-IR bands. We convert from Vega magnitudes to flux densities by using zero points of 309.540, 171.787, 31.674, and 8.363 Jy for the *W1*, *W2*, *W3*, and *W4* bands, respectively. As recommended in the *WISE* “Cautionary Notes,” we only use those matches for which the contamination and confusion flag *cc_flags* is equal to 0000, in order to exclude sources that are flagged as spurious in any band.<fnref n="11" /> Sources in the catalog with less than a 2σ significance in a given band are assigned an upper limit given by the integrated flux density measurement plus two times the estimate uncertainty. If the flux density estimate is negative, then the upper limit is defined as two times the measurement uncertainty.<fnref n="12" />

Similar to the case of *WISE*, the near-IR information is included via UKIDSS cross-matching with SDSS QSOs. UKIDSS observes the sky in the four *Y*, *J*, *H* and *K* bands (from ≈1 µm to ≈2.2 µm), and fluxes are provided in Jy.

## 2.3. Herschel Far-IR Photometry and Maps

For our work, we make use of data from the first public release of the *Herschel*-Astrophysical Terahertz Large Area Survey (H-ATLAS; Eales et al. 2010), the largest extragalactic key project carried out in open time with the *Herschel Space Observatory*. It was allocated 600 hr of observing time and covered about 600 deg² of sky in five photometric bands with the Photodetector Array Camera and Spectrometer (PACS; 100 and 160 µm) and the Spectral and Photometric Imaging Receiver (SPIRE; 250, 350 and 500 µm) instruments. Specifically, we use a subset of the PACS and SPIRE image products provided by the H-ATLAS team<fnref n="13" /> (see Valiante et al. 2016 for

a detailed discussion of the map processing and properties), covering the three equatorial fields surveyed by the GAMA spectroscopic survey, GAMA09 (G09), GAMA12 (G12), and GAMA15 (G15), each one of about ≈54 deg². All of the maps utilized in our analysis have had the large-scale background subtracted through the application of the Nebuliser algorithm, so that extended emission from the Galaxy or from any other extragalactic signal on scales larger than ≈3' is removed. In addition, SPIRE maps have also been convolved with a matched-filter derived from the PSF, such that the value in each pixel represents the best estimate of the measured flux that a point source would have in that position. We recall that the SPIRE maps are provided in Jy/beam units while PACS maps are released in Jy/pixel. The pixel sizes are ≈3", 4", 6", 8", and 12" at 100, 160, 250, 350, and 500 µm respectively, while the angular resolutions at the different frequencies are provided in Table 1. The three GAMA fields, as well as the SDSS footprint, are shown in Figure 2. In the case of our baseline QSO sample outlined in Section 2.1, approximately 2300 QSOs fall within the H-ATLAS footprint in each redshift bin.

## 2.4. Additional Maps: Planck and AKARI

We complement the far-IR information from *Herschel* with *AKARI* data at higher frequencies and with *Planck* observations at similar and lower frequencies.

The *Planck* satellite has scanned the entire sky in nine frequency bands from 30 to 857 GHz. Its broad spectral coverage allows it to characterize and separate the Galactic and extragalactic emissions, making it particularly suited for astrophysical studies. For a detailed description of the public data release and the scientific results of the *Planck* team, we refer the reader to Planck Collaboration et al. (2016a). In this work, we make use of the publicly available<fnref n="14" /> 2015 single-frequency intensity *Planck* maps from 30 to 857 GHz (see Table 1 for a summary of the maps' properties). Maps are provided in the HEALPix<fnref n="15" /> (Górski et al. 2005) format at a resolution of $N_{side} = 2048$ size. Maps up to 353 GHz are converted from units of $K_{CMB}$ to MJy/sr with the unit conversion factors given in the *Planck* 2015 release (Planck Collaboration et al. 2016b). In addition, we utilize *Planck* Cosmic Infrared Background (CIB) maps obtained with the Generalized Needlet Internal Linear Combination method (GNILC; Planck Collaboration et al. 2016c) to estimate the contamination of CIB in our measurements of QSO fluxes at 353, 545, and 857 GHz. GNILC is a generalization of the internal linear combination method capable of extracting multiple correlated sources of emission in a multifrequency data set, taking into account the properties of such signals as a function of both direction in the sky and angular scale using needlets. We refer the reader to Remazeilles et al. (2011) for further details on the GNILC method.

In order to bridge the observational gap between the mid-IR band covered by *WISE* and the far-IR portion of the spectrum sampled by *Herschel*, we also make use of the all-sky far-IR data from the *AKARI* satellite (Doi et al. 2015). *AKARI* carried out a photometric survey over more than 99% of the sky in four bands centered at 65, 90, 140, and 160 µm, with spatial resolutions ranging from 1' to 1'5. The absolute flux calibration was performed through a comparison of *AKARI* images with

<footnote n="11">¹¹ See http://wise2.ipac.caltech.edu/docs/release/allwise/expsup/sec2_2.html#comp.</footnote>

<footnote n="12">¹² See http://wise2.ipac.caltech.edu/docs/release/allwise/expsup/sec2_1a.html.</footnote>

<footnote n="13">¹³ Available at http://www.h-atlas.org/public_data/.</footnote>

<footnote n="14">¹⁴ http://pla.esac.esa.int/pla/</footnote>

<footnote n="15">¹⁵ http://healpix.jpl.nasa.gov</footnote>
<table><caption>**Table 1**

Overview of the Properties of the Maps Used in the Stacking Analysis</caption><thead><tr><th></th><th colspan="10">Planck</th><th colspan="4">Herschel</th><th>AKARI</th></tr></thead><tbody><tr><th>Frequency (GHz)</th><td>30</td><td>44</td><td>70</td><td>100</td><td>143</td><td>217</td><td>353</td><td>545</td><td>857</td><td>600</td><td>857</td><td>1200</td><td>1870</td><td>3000</td><td>3300</td></tr><tr><th>Wavelength (μm)</th><td>10600</td><td>6800</td><td>4300</td><td>3000</td><td>2100</td><td>1350</td><td>830</td><td>540</td><td>350</td><td>500</td><td>350</td><td>250</td><td>160</td><td>100</td><td>90</td></tr><tr><th>Beam FWHM</th><td>32′.3</td><td>27′.12</td><td>13′.31</td><td>9′.68</td><td>7′.30</td><td>5′.02</td><td>4′.94</td><td>4′.83</td><td>4′.64</td><td>35″.2</td><td>24″.2</td><td>17″.8</td><td>13″.7</td><td>11″.4</td><td>1′.30</td></tr></tbody></table>

Note. The observational wavelengths of the *Planck* satellite have been computed from the effective frequencies of the observational bands (Planck Collaboration et al. 2014a, 2014b).

<figure><img src="image_5.png" /><figcaption>**Figure 2.** SDSS (purple shaded area) and Herschel-ATLAS (pink shaded area) footprints in galactic coordinates. The background image shows the *Planck* thermal dust map from Commander, one of the component separation algorithms adopted by the *Planck* collaboration.</figcaption></figure>

COBE/DIRBE data sets. A reliable flux calibration down to low intensities is provided only for the 90 µm channel, which also has the lowest noise level. For this reason, we only consider the 90 µm AKARI map for the present analysis. The map<fnref n="16" /> is provided in the HEALPix format at a resolution of $N_{side} = 4096$, corresponding to 0.86 per pixel, in MJy/sr units.

## 2.5. Sample Selections

In this subsection, we describe the construction of our main QSO sample and subsamples that are used in our analysis.

All of the QSO samples analyzed in this work are extracted from the merged SDSS DR10 and DR12 QSO catalog introduced in Section 2.1. From this parent QSO sample, which contains 386,245 sources, we create three samples by gradually enforcing additional selection criteria:

1. AllQSO sample: Our main QSO sample includes 3,13,004 sources within the redshift range of $1 \le z \le 5$, which we divide into three redshift bins. As mentioned in Section 2.1, we additionally use the radio-band information from FIRST to impose a radio-lout cut that minimizes the emission due to synchrotron. Specifically, we reject any object in the parent SDSS QSO catalog that has a counterpart in FIRST, reducing the overall sample size by approximately 4%. These cuts define our main QSO sample, also referred to as AllQSO in the following figures and tables. The angular distribution of these QSOs is shown as the purple shaded area in Figure 2. Note that

Herschel observations only cover a fraction of the entire SDSS footprint. However, unless otherwise stated, when estimating the average QSO SED, we combine all of the available data sets, i.e., we do not only consider the objects that fall in the H-ATLAS footprint. By doing so, we assume the different patches of the sky to be statistically independent. We test the representativity of the H-ATLAS subsample in Appendix A.

2. W4 sample: Most of the optical SDSS QSOs are undetected at an observed wavelength of 22 µm. In creating this subsample, we first start from our main QSO sample (AllQSO) and require a detection in the W4 channel of WISE, i.e., $S/N \ge 2$. As discussed later, this selection implies a stronger mid-IR emission and, as a consequence, higher emission in the far-IR and UV-optical bands. The number of QSOs selected in this case is approximately 10% of the main QSO sample and depends on the redshift bin as well as on the specific data sets. For example, in the case of SDSS, we find ~13,000 objects in the first redshift bin and ~8000 in the last one, while for Herschel the numbers are roughly ~130/200 for the first/last bin, respectively.

3. W4 + W4/r > 100 sample: Since we are also interested in studying the QSOs that live in obscured environments, and considering that the W4 selection also increases the mean UV-optical emission (see Figure 5), we add a further constraint to the W4 subsample by requesting individual QSOs with the lowest UV-optical emission. In order to achieve this, we impose a ratio between the mid-IR (W4) emission and the optical one (r band) greater than 100 in flux density. The number of SDSS QSOs that fulfills both criteria is around 6100, while the ones falling in the H-ATLAS footprint are approximately 100. If we raise the $S_{W4}/S_r$ threshold to higher values, e.g., $S_{W4}/S_r > 500$, we are left with only a handful of sources; therefore, we set the threshold ratio to 100.

The three samples outlined above include QSOs that are both undetected (the majority) and detected in the far-IR bands covered by Herschel. To assess the differences between the reconstructed SED of submillimeter-detected and -undetected QSOs, we proceed in the following way. We start from the H-ATLAS catalogs (Valiante et al. 2016; Bourne et al. 2016) and cross-match them with each one of the three QSO samples described above, retaining only those objects that are detected at Herschel frequencies. By doing so, we create three more restrictive QSO samples with the additional property of containing only submillimeter-detected QSOs. Note that, differently from the previous case, we limit the analysis of the sole QSOs that fall in the H-ATLAS footprint. Moreover, in

<footnote n="16">¹⁶ Available at http://cade.irap.omp.eu/dokuwiki/doku.php?id=akari.</footnote>
<table><caption>**Table 2**
The Optically Selected Quasar Samples that We Create for This Work, along with Their Redshift Coverage</caption><thead><tr><th>Sample</th><th>Redshift Bin</th><th>N<sub>SDSS</sub></th><th>N<sub>WISE</sub></th><th>N<sub>Herschel</sub></th></tr></thead><tbody><tr><td rowspan="3">AllQSO</td><td>1 &lt; z &lt; 2.15</td><td>62,955</td><td>22868</td><td>2460</td></tr><tr><td>2.15 &lt; z &lt; 2.5</td><td>87,644</td><td>23475</td><td>2345</td></tr><tr><td>2.5 &lt; z &lt; 5</td><td>92,042</td><td>23391</td><td>2263</td></tr><tr><td>AllQSO + H-ATLAS detected</td><td>1 &lt; z &lt; 5</td><td>147</td><td>64</td><td>147</td></tr><tr><td rowspan="3">W4</td><td>1 &lt; z &lt; 2.15</td><td>13,285</td><td>6322</td><td>125</td></tr><tr><td>2.15 &lt; z &lt; 2.5</td><td>8277</td><td>3950</td><td>198</td></tr><tr><td>2.5 &lt; z &lt; 5</td><td>8141</td><td>3885</td><td>194</td></tr><tr><td>W4 + H-ATLAS detected</td><td>1 &lt; z &lt; 5</td><td>36</td><td>34</td><td>36</td></tr><tr><td>W4 + W4/r &gt; 100</td><td>1 &lt; z &lt; 5</td><td>6171</td><td>2974</td><td>103</td></tr><tr><td>W4 + W4/r &gt; 100 + H-ATLAS detected</td><td>1 &lt; z &lt; 5</td><td>7</td><td>6</td><td>7</td></tr></tbody></table>

**Note.** Note that the number of stacked sources depends on the specific survey and wavelength. We report here the number of sources stacked in SDSS, WISE, and Herschel as representative of the surveys most sensitive to the optical, mid-IR, and starburst components of the SED, respectively.

this case, there is no need to perform the stacking analysis on the H-ATLAS maps to infer the submillimeter density flux as the best-fit PACS/SPIRE fluxes are already provided in the H-ATLAS catalogs. Requiring submillimeter detection and limiting the analysis over the H-ATLAS footprint severely reduce the number of selected QSOs. For this reason, we decided to consider only one redshift bin between 1 < z < 5. The number of sources included in these cross-matched samples is 147 for the main QSO sample, 36 for the “W4” selection, and 7 in the case of the “W4 W4/r > 100” subsample. Unfortunately, the number of sources is not sufficient to carry out the stacking analysis with AKARI or Planck. Table 2 summarizes the number of objects available in the SDSS, WISE, and Herschel surveys for all of the QSO samples we create.

### 3. Methods

The main focus of this work is the estimation of the mean SED of optically selected QSOs from the near-UV to the millimeter band and, through that, constrain the average star formation and AGN luminosity as a function of cosmic time. Since the typical flux density of optically selected QSOs is fainter than Herschel, Planck, and AKARI detection limits, we rely on a stacking analysis to recover photometric information for the undetected sources. Stacking is a powerful technique that allows the extraction of signal hidden in noisy data sets, albeit at the price of losing information about the individual sources. In particular, stacking can be shown to be an unbiased maximum likelihood estimator of the average flux density of a catalog, provided that background noise is Gaussian or that the data are confusion-limited (Marsden et al. 2009). Note that we assume the different patches of the sky to be statistically independent for the purpose of stacking.

#### 3.1. Herschel Stacking

We perform the stacking analysis in each redshift bin by cutting out 5' × 5' stamps at the locations of all QSOs in our

sample, thus creating a pixel cube. The number *N* of pixels in each stamp then depends on the pixel resolution, and as such, it varies between different bands. The mean of the Herschel maps is not zero because of the application of the Nebuliser algorithm. As recommended by Valiante et al. (2016), we set to zero the mean of each map before performing the stacking analysis. After extracting the cutouts at the different redshift bins and wavelengths, we combine them into a single image, a process known as data cube reduction. In literature, this is usually achieved by computing the mean or median flux of all the stamps in a given pixel. In this work, we use the mean in order to represent the typical flux density of the sample. To this end, we use the noise maps extracted at the same QSO locations to perform a weighted average of each pixels in the stacked image:

$$ \bar{S}_{\nu} = \frac{\sum_{i}^{N_{\text{QSO}}} w_{\nu}^{i} S_{\nu}^{i}}{\sum_{i}^{N_{\text{QSO}}} w_{\nu}^{i}}. \qquad (1) $$

Here, $S_{\nu}^{i}$ is the pixel flux density in the stamp of the QSO *i* at a band *ν*, $w_{\nu}^{i}$ is the inverted noise flux density in the same pixel, and the sum runs over all $N_{\text{QSO}}$ sources in a given redshift bin. We adopt two different approaches for measuring the stacked flux from the PACS and SPIRE cutouts. For the PACS stacks, which are in units of Jy/pixel, we perform aperture photometry and integrate the flux density out to 5''6 and 10''2 for the 100 and 160 µm channels, respectively. The sky background is calculated in a ring of inner and outer radii of 35'' and 45'' and subtracted from the previous estimate. In the case of the SPIRE observations, which are in the units of Jy/beam, we fit a 2D Gaussian (plus a constant background) to the image stack and take its peak value as the flux density estimate.

We employ a bootstrap approach to estimate the $1\sigma$ errors on the stacked flux in both the PACS and SPIRE bands for each redshift bin. To this end, we randomly resample (with replacement) the sources and perform the stacking procedure 100 times for each redshift bin and frequency band. The uncertainty in the flux is then computed as the standard deviation of the measured stacked mean flux in these 100 trials. An example of stacked cutouts in different redshift bins and PACS/SPIRE frequencies at the positions of our main QSO sample is shown in Figure 3.

#### 3.2. Planck and AKARI Stacking

Planck and AKARI data sets differ from the Herschel maps in that they are full sky and have a lower angular resolution. For these reasons, as well as to avoid projection effects associated with the conversion from the HEALPix full-sky to the flat-sky geometry and pixelization, we follow Soergel et al. (2017) and rely on a different approach based on matched filtering to extract the QSO emission in the Planck and AKARI maps. Matched filtering, first proposed by Haehnelt & Tegmark (1996) in the context of kinetic Sunyaev-Zel'dovich studies, is a powerful technique that allows the construction of an optimal spatial filter to recover the signal of interest (the QSO flux density in our case) in the presence of a noisy background (the primary cosmic microwave background (CMB), foregrounds, and instrumental noise). One of the assumptions at the base of the matched-filter approach is that the spatial signature of the sought-after signal is known. Given that the QSO host halos are unresolved at all Planck and AKARI frequencies (100 kpc at
<figure><img src="image_6.png" /><figcaption>Figure 3. Stacked thumbnail images from H-ATLAS PACS and SPIRE data at the locations of the sources included the main QSO sample (AllQSO). The number of objects stacked in each redshift bin is approximately 2300. Each box is 5' × 5', and the units are mJy/pixel. The orange circle in each panel represents the instrument PSF at the various bands. Its diameter corresponds to the FWHM size.</figcaption></figure>

$z \sim 2$ subtend an angle of about 0.2, much smaller than their respective angular resolutions; see Table 1), we can safely assume the QSO spatial profiles to be approximated by the instrumental beam profile. Note that in this way the extracted fluxes are model-independent.

The filtered map harmonic coefficients $a_{\ell m}^{\text{filt}}$ are related to the unfiltered ones through $a_{\ell m}^{\text{filt}} = F_{\ell} a_{\ell m}$ (Schäfer et al. 2006), where the filter function reads as

$$F_{\ell} = \left[ \sum_{\ell'} \frac{(2\ell + 1)\tilde{B}_{\ell}^2}{4\pi C_{\ell}^{\text{tot}}} \right]^{-1} \frac{\tilde{B}_{\ell}}{C_{\ell}^{\text{tot}}}. \quad (2)$$

In the above equation, $\tilde{B}_{\ell} = 2\pi\sigma B_{\ell} p_{\ell}$, $B_{\ell}$ is the instrumental beam window function, $p_{\ell}$ is the pixel window function, and $\sigma = \theta_{\text{FWHM}}/\sqrt{8 \ln 2}$ is a factor that accounts for the beam normalization to unity in real space. Instrumental beams are assumed to be Gaussian with the FWHM shown in Table 1. We checked that this approximation was sufficiently accurate by repeating the analysis for the Planck frequencies while adopting the $B_{\ell}$ of the Planck Reduced Instrument Model<fnref n="17" /> and finding no significant difference with our baseline results. For the GNILC maps, we assumed a Gaussian beam with 5'

<footnote n="17">¹⁷ https://wiki.cosmos.esa.int/planckpla2015/index.php/The_RIMO</footnote>
FWHM following Planck Collaboration et al. (2016c). The power spectrum describing all of the noise sources is taken to be the total power spectrum of the unfiltered maps, $C_{\ell}^{\text{tot}}$, since the QSOs' emission signal is negligible with respect to the primary CMB, foreground, and instrumental noise. The total power spectrum is extracted from the full-sky maps using the XPOL code<fnref n="18" /> (Tristram et al. 2005) after applying an apodized (2°) 40% Galactic mask to mitigate the effect of Galactic foregrounds and a point-source mask obtained by the combination of the Planck 2015 Catalog of Compact Sources and the AKARI Far-Infrared Surveyor Bright Source Catalog. As a conservative choice to mitigate any residual large-scale foreground and following Soergel et al. (2017), we further apply an additional smooth high-pass filter to remove large angular scales at $\ell \lesssim 300$. After filtering the full-sky maps, we can estimate the flux density for the QSO at position i by multiplying the pixel value $I_{\nu}^{\text{filt}}(\hat{n}_i)$ of the filtered map by the beam area,

$$S_{\nu}^{i} = I_{\nu}^{\text{filt}}(\hat{n}_{i}) \int d\Omega B_{\nu}(\theta). \quad (3)$$

Then, the mean flux density is computed as

$$\bar{S}_{\nu} = \frac{1}{N_{\text{QSO}}} \sum_{i} S_{\nu}^{i}, \quad (4)$$

while the frequency-frequency covariance reads as

$$C_{\nu\nu'} = \frac{1}{N_{\text{QSO}}(N_{\text{QSO}} - 1)} \sum_{i} (S_{\nu}^{i} - \bar{S}_{\nu})(S_{\nu'}^{i} - \bar{S}_{\nu'}), \quad (5)$$

from which we extract the $1\sigma$ error bars as the square root of the diagonal.

### 3.3. Mean Flux Densities from Cataloged Sources

We also estimate the QSO fluxes at the catalog level using the available SDSS, UKIDSS, and WISE photometry. In order to be consistent with the stacking analysis discussed above, we report the mean flux density as the representative quantity for the QSO sample at the different bands. For the SDSS and UKIDSS fluxes, we estimate the $1\sigma$ uncertainties as the standard error deviation of the flux density distribution of the individual sources divided by the square root of the number of objects. In the case of WISE photometry, the simple mean of the QSOs' flux distribution would yields a biased value because the fraction of sources with upper limits in the W3 and W4 channels range from 32% to 44% and from 70% to 81% in the three redshift bins, respectively (for our baseline QSO sample). Given that the upper limits show a random-enough distribution, we exploit the non-parametric Kaplan-Meier estimator (Feigelson & Nelson 1985) to take into account the upper limits and calculate the WISE mean flux density. To assess the WISE uncertainties, we conservatively consider two extreme cases, one where we calculate the mean flux assuming all sources to be detected (upper error bar), and one where we compute the mean flux assigning a null flux to the sources for which only upper limits are provided (lower error bar).

### 3.4. Modeling the Total SED

In order to study the average physical properties of our selected QSOs, we need to perform an SED decomposition. We are mainly interested in deriving the bolometric luminosity for each component and comparing their variation with redshift and for the different subsamples described in the next subsection. Therefore, a detailed analysis of the 0.5–500 µm SEDs is beyond the purpose of this paper. We simplified the SED emission in three components: AGN emission from the central engine of the QSOs, a starburst emission from the star formation activity in the galaxy host, and a mid-IR emission from the torus surrounding the central engine.

1. UV-optical component: To represent the UV-optical emission (residual after obscuration) from normal quasars, we selected the AGN intrinsic template based on Elvis et al. (1994), developed by Xu et al. (2015, see their Appendix C). We assume that this SED remains unchanged for all considered subsamples and redshifts.

2. Starburst component: We used the Pearson et al. (2013) two-graybody template that was found to best describe the intermediate-high-redshift starburst galaxies observed by Herschel. Although simplistic, it provides enough information to perform a rough estimation of a photometric redshift and to fit the stacked far-IR emission (see Section 2.5).

3. Torus component: Classical SED decomposition approaches to mid-IR observations, i.e., torus emission, use rigid templates that do not provide a good fit to our stacked QSO SEDs (e.g., Elvis et al. 1994; Richards et al. 2006; Lani et al. 2017). Following Lyu et al. (2017), we decided to fit the mid-IR emission by taking into account three graybody components in the following temperature ranges: 1000–2000 K, 700–1000 K, and 200–700 K. For all of them, we fixed the dust emissivity to $\beta = 1.5$ as found in studies of high-z AGNs (e.g., Beelen et al. 2006; Lyu et al. 2017).

Taking into account these three components, we perform a two-step SED decomposition. (i) We exclude the mid-IR data (2 µm < $\lambda$ < 90 µm), focusing only on the optical ($\lambda$ < 2 µm) and far-IR data ($\lambda$ > 90 µm). We fit at the same time the UV-optical and starburst components, i.e., we determine the normalization (vertical shift) for each component and the common redshift (horizontal shift or a kind of mean photometric redshift). The curve_fit, a nonlinear least-squares algorithm from the scipy package, allows us to derive the best-fit parameters for the sum of both components to the measured data. The estimated photometric redshifts are very close to the spectroscopic mean redshift in each bin, as expected. (ii) We add the mid-IR measurements, and we attempt to fit them with the addition of a torus component described as three graybody functions. We start with the UV-optical and starburst components, considering the best-fit parameters determined in the previous step. The best-fit parameters from the previous step remains fixed at this step. Then, we sum the three graybody functions at the same best-fit redshift from the previous step. Finally, using again curve_fit, we determine the best-fit parameters for the three graybody functions.

Obviously, the mid-IR data fit is degenerate: six free parameters (three temperatures plus three normalization values) versus four photometric points. However, as the UV-optical

<footnote n="18">¹⁸ https://gitlab.in2p3.fr/tristram/Xpol</footnote>
<figure><img src="image_7.png" /></figure>

**Figure 4.** Averaged SEDs obtained for the main sample of QSOs in the three redshift bins with the following photometry: SDSS (dark blue), UKIDSS (light blue), WISE (green), AKARI (yellow), Herschel (SPIRE + PACS, red), and Planck (dark red). The CIB contamination estimated with the GNILC map has been subtracted from the average flux density at 540 and 830 µm. Unfilled markers at those frequencies represent the values of the measurements prior to this correction. CIB-subtracted measurements at 350 µm are not shown as they are affected by systematics in the GNILC maps. In the case of AKARI and Planck, lower triangles indicate upper limits. As described in detail in the text, the SEDs are modeled using three components: UV-optical emission (dashed blue; Xu et al. 2015), torus emission (dotted green; three graybody components) and a starburst emission (dashed magenta; Pearson et al. 2013). Note that the Planck data (empty/filled dark red squares) were not used in the fitting process.

and starburst components are already fixed and the temperatures are restricted to a narrow range, the algorithm is able to converge to a reasonable fit, at least, good enough to estimate the mid-IR bolometric luminosity. Consequently, the estimated temperatures are reasonable but very uncertain and, for this reason, they are neither presented nor discussed here. We plan to study in detail the mid-IR range in a future work. See Figure 4 for an example of the SED fitting for the main QSO sample in the three redshift bins. Note that the *Planck* data (empty/filled dark red squares) are just plotted for a visual comparison but not used in the fitting process. As commented later on, we defer the analysis of CMB and radio frequencies to future work.

Once the fit is obtained, we estimate three bolometric luminosities, one for each of the components: for the AGN component, we consider only the UV-optical SED, for the mid-IR the sum of the three graybody fitted functions, and for the far-IR the starburst SED. This process is repeated for each sample and redshift bin, as discussed in Section 4.1.

As a robustness test, we check that our derived averaged SEDs are independent of the exact applied methodology (see Appendix C for more details).

# 4. Results and Discussion

In this section, we discuss the main findings of our analysis. We start by showing the reconstructed SED for the different QSO samples and reporting the inferred bolometric luminosities of the various QSO/host-galaxy components. We then compare our mean SED with those present in the literature and

conclude by discussing the results in terms of the galaxy–SMBH coevolutionary framework.

## 4.1. Reconstructed SEDs

In Figure 4, we present the measured average SEDs for our main QSO sample in the three redshift bins (see Table 3 for their estimated photometric redshifts), along with the best-fit SEDs and their decomposition in terms of the single components described in Section 3.4. In all three cases, we can easily identify the UV-optical emission associated with the AGN, the starburst component at submillimeter wavelengths, as well as the mid-IR bump sourced by the heated torus emission.

Note also that although there is a qualitatively good agreement between the extrapolated SEDs and *Planck's* measured fluxes at millimeter and radio wavelengths, we do not make use of *Planck* data in the fitting process as the flux derived from the stacking appears to generally overestimate the flux density. The overestimation of the flux density becomes more important at the shorter wavelengths of *Planck*, in particular at those overlapping with the *Herschel* bands, i.e., 540 and 830 µm. Residual miscalibration between *Herschel* and *Planck* cannot account for the difference in the measured flux density. Thus, we suspect that the coarser resolution and higher noise levels might lead to a measured flux density that systematically picks up background graybody-like signal not associated with the QSOs, e.g., CIB or residual Galactic dust. If we subtract the estimated CIB contamination obtained from the stacking of the CIB GNILC maps at 540, 830 µm, we recover
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Averaged SEDs obtained for the different samples of QSOs in the redshift range $2.15 < z < 2.5$ (or $1 < z < 5$ for the H-ATLAS detected cases; see the text for more details). Filled symbols correspond to SEDs where only sources detected in H-ATLAS catalogs are considered (no map-based stacking was needed in this case). On the contrary, the SEDs estimated considering also undetected H-ATLAS sources (stacking applied) are represented with empty symbols. Lower triangles indicate upper limits.</figcaption></figure>

<caption>**Table 3**</caption>

Inferred Bolometric Luminosities of the Optical, Mid-IR, and Starburst Components and Estimated Photometric Redshift for the Different Samples and Redshift Bins

<table><thead><tr><th>Sample</th><th>Redshift Bin</th><th>Photometric z</th><th>log<sub>10</sub> L<sub>opt</sub> (L<sub>☉</sub>)</th><th>log<sub>10</sub> L<sub>milR</sub> (L<sub>☉</sub>)</th><th>log<sub>10</sub> L<sub>SB</sub> (L<sub>☉</sub>)</th></tr></thead><tbody><tr><td rowspan="3">AllQSO</td><td>1 &lt; z &lt; 2.15</td><td>1.79<sup>+0.01</sup><sub>-0.01</sub></td><td>12.32<sup>+0.01</sup><sub>-0.01</sub></td><td>12.30<sup>+0.17</sup><sub>-0.14</sub></td><td>11.84<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>2.15 &lt; z &lt; 2.5</td><td>2.15<sup>+0.01</sup><sub>-0.01</sub></td><td>12.36<sup>+0.01</sup><sub>-0.01</sub></td><td>12.34<sup>+0.22</sup><sub>-0.15</sub></td><td>11.79<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>2.5 &lt; z &lt; 5</td><td>2.65<sup>+0.01</sup><sub>-0.01</sub></td><td>12.57<sup>+0.01</sup><sub>-0.01</sub></td><td>12.52<sup>+0.24</sup><sub>-0.20</sub></td><td>11.96<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>AllQSO + H-ATLAS detected</td><td>1 &lt; z &lt; 5</td><td>1.59<sup>+0.07</sup><sub>-0.07</sub></td><td>12.52<sup>+0.05</sup><sub>-0.06</sub></td><td>12.40<sup>+0.17</sup><sub>-0.07</sub></td><td>12.60<sup>+0.05</sup><sub>-0.05</sub></td></tr><tr><td rowspan="3">W4</td><td>1 &lt; z &lt; 2.15</td><td>1.68<sup>+0.02</sup><sub>-0.02</sub></td><td>12.61<sup>+0.01</sup><sub>-0.01</sub></td><td>12.57<sup>+0.01</sup><sub>-0.02</sub></td><td>12.02<sup>+0.05</sup><sub>-0.06</sub></td></tr><tr><td>2.15 &lt; z &lt; 2.5</td><td>2.15<sup>+0.03</sup><sub>-0.03</sub></td><td>12.79<sup>+0.01</sup><sub>-0.01</sub></td><td>12.79<sup>+0.01</sup><sub>-0.01</sub></td><td>12.11<sup>+0.05</sup><sub>-0.05</sub></td></tr><tr><td>2.5 &lt; z &lt; 5</td><td>2.65<sup>+0.03</sup><sub>-0.03</sub></td><td>12.99<sup>+0.01</sup><sub>-0.01</sub></td><td>13.03<sup>+0.01</sup><sub>-0.01</sub></td><td>12.32<sup>+0.04</sup><sub>-0.05</sub></td></tr><tr><td>W4 + H-ATLAS detected</td><td>1 &lt; z &lt; 5</td><td>1.53<sup>+0.14</sup><sub>-0.14</sub></td><td>12.72<sup>+0.10</sup><sub>-0.11</sub></td><td>12.64<sup>+0.11</sup><sub>-0.12</sub></td><td>12.62<sup>+0.09</sup><sub>-0.10</sub></td></tr><tr><td>W4 + W4/r &gt; 100</td><td>1 &lt; z &lt; 5</td><td>2.21<sup>+0.27</sup><sub>-0.27</sub></td><td>12.83<sup>+0.01</sup><sub>-0.01</sub></td><td>13.29<sup>+0.01</sup><sub>-0.01</sub></td><td>12.51<sup>+0.08</sup><sub>-0.10</sub></td></tr><tr><td>W4 + W4/r &gt; 100 + H-ATLAS detected</td><td>1 &lt; z &lt; 5</td><td>3.80<sup>+0.04</sup><sub>-0.04</sub></td><td>12.27<sup>+0.12</sup><sub>-0.14</sub></td><td>13.49<sup>+0.19</sup><sub>-0.23</sub></td><td>12.95<sup>+0.12</sup><sub>-0.13</sub></td></tr><tr><td>HotDOG</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>SDSS 022052.11+013711.1</td><td>z = 3.138</td><td>...</td><td>11.02<sup>+0.08</sup><sub>-0.09</sub></td><td>12.42<sup>+0.16</sup><sub>-0.19</sub></td><td>12.63<sup>+0.07</sup><sub>-0.08</sub></td></tr></tbody></table>

values of the flux densities compatible with the expected values of our best-fit SEDs. However, the GNILC maps at 340 µm are highly dominated by spurious systematic effects. The estimated flux density of this map, in fact, exceeds the value extracted from the *Planck* single wavelength map, which is an unphysical situation. Hints that the GNILC CIB flux density might be overestimated can also be observed in the $1 < z < 2.25$ redshift bin at 540 µm. Despite templates of the Galactic dust emission having been published by the *Planck* collaboration using several component separation methods, we did not

perform stacking using these maps as on small angular scales they are heavily contaminated by CIB emission (Planck Collaboration et al. 2016b, 2016c). For all the reasons just discussed, we decided to postpone an accurate inclusion of *Planck* data to a future work. We stress that our findings also suggest that any use of the *Planck* data for analysis similar to ours should be treated with caution, particularly in the highest frequency channels.

The estimated mean SEDs for all of the samples described in Section 2.5 are shown and compared in Figure 5. We remind
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Comparison of the mean SED reconstructed in this work in the redshift range $2.15 < z < 2.5$ (solid black line) with classic templates from the literature. We show the SED from Richards et al. (2006) obtained using SDSS-selected QSOs in combination with *Spitzer* data (dashed red line), the template by Elvis et al. (1994) as revisited by Xu et al. (2015) from star formation in the host galaxy (dotted–dashed magenta), and the WDD template by Lyu et al. (2017). Moreover, we have also plotted the most recent intrinsic AGN SED derived by Lani et al. (2017; green line).</figcaption></figure>

that the SEDs for the $W4+W4/r > 100$ sample and all of the H-ATLAS detected samples are calculated over the redshift range $1 < z < 5$, while the remaining SEDs are obtained for the intermediate-redshift bin $2.15 < z < 2.5$. When looking at the samples comprising both undetected and detected Herschel (empty symbols), we notice that requiring an $S/N \ge 2$ in the W4 band results in an overall shift to higher flux densities with respect to the values found for the main QSO sample. Instead, imposing a further constraint on the mid-IR to optical flux ratio $S_{W4}/S_r > 100$ translates to similar or slightly smaller UV-optical flux densities and higher flux values at mid- and far-IR wavelengths when compared to the baseline QSO sample, while pushing down the measured fluxes at $\lambda \lesssim 100$ µm with respect to the “W4” case. As expected, requiring a detection at far-IR frequencies affects the reconstructed SED by boosting the far-IR emission from a few mJy to a few tens of mJy, almost independently of the specific selection criteria (filled symbols). A similar trend is seen at wavelengths $\lambda \lesssim 100$ µm, where the average fluxes for the submillimeter-detected sources are enhanced, although in a negligible way for $W4+S_{W4}/S_r > 100$ selected sources in the UV-optical band range. In Figure 5, we also plot the SED of the only hot-dust-obscured galaxy (HotDOG) from Fan et al. (2016), which is included in our SDSS-based QSO catalog (SDSS 022052.11 +013711.1). This object, situated at a redshift of $z = 3.138$, is a violently star-forming galaxy featuring a high-temperature dust component with a strong submillimeter emission and is characterized by a W4 to r ratio of about 1000.

In Table 3, we report the estimated photometric redshift and the inferred bolometric luminosities obtained by applying the fitting method outlined in Section 3.4 to each SEDs. As can be seen, the luminosities associated with the three QSO components tend to increase over cosmic time, with similar

UV-optical and mid-IR luminosities (of course when $S_{W4}/S_r$ criterion is not enforced). Once again, as expected, the submillimeter detection criterion results in starburst luminosities above $\log_{10}(L_{SB}/L_\odot) \gtrsim 12.60$.

## 4.2. Comparison with SEDs from Literature

In Figure 6, we illustrate the stacked rest-frame SED for the full QSO sample from our analysis (in the redshift range $2.15 < z < 2.5$) and compare it to classic AGN templates from the literature. Specifically, we consider the template by Richards et al. (2006; red dashed line) obtained from a sample of 259 SDSS-selected quasar with *Spitzer* data; the template by Elvis et al. (1994; magenta dotted–dashed line) based on a Palomar–Green (PG) sample of radio-quiet QSOs, as revisited by Xu et al. (2015) from star formation in the host. Moreover, we have also plotted the most recent warm-dust-deficient template by Lyu et al. (2017; cyan triple-dotted–dashed line) applying to high-luminosity PG QSOs and the intrinsic AGN SED derived by Lani et al. (2017; green line). This latest SED corresponds to the median SED after subtracting the star formation contribution using a PAH-based method to the PG QSO sample. All of the templates have been normalized to unity at a rest-frame wavelength of 0.3 µm.

As expected, our SED is very similar in shape to the other templates in the optical/UV part of the spectrum. At around 1–2 µm, it shows a little excess with respect to the other templates, likely due to the combination of the emissions from old stars (e.g., RGB/AGB) in the host galaxy and from dust near the sublimation temperature ($T \sim 2000$ K). Around 4–6 µm, our SED features a prominent bump, followed by an abrupt decrease out to 10 µm. This articulated behavior is not present in the Richards template, while it is just hinted at in the
<figure><img src="image_10.png" /></figure>

**Figure 7.** SFR vs. $L_{\text{AGN}}$ plane between the bolometric AGN luminosity and star formation rate in the host galaxy at $z \gtrsim 2$. Dashed lines illustrate three typical evolutionary tracks (forward time direction indicated by arrows) from the in situ coevolution scenario by Lapi et al. (2014, 2017; see also Mancuso et al. 2016, 2017) corresponding to SFR ~ 400, 1200, and 4000 $M_{\odot}$ yr$^{-1}$. The dotted line represents the locus where the bolometric luminosities from AGN emission and from star formation are equal; galaxies to the right of such a line are prone to be affected by substantial AGN activity (see text). The thick solid line shows the average relationship and the shaded areas render its $1\sigma$ (light gray), $2\sigma$ (gray), and $3\sigma$ (dark gray) variance, computed as in Mancuso et al. (2016) by taking into account the number density of galaxies and AGNs, the relative time spent by individual objects in different portions of the tracks, and a far-IR detection threshold around SFR ~ 100 $M_{\odot}$ yr$^{-1}$. Literature data are from Stanley et al. (2015, 2017; circles), Netzer et al. (2016; squares), Venemans et al. (2018; pentagons), Fan et al. (2016; inverse triangles), Xu et al. (2015; triangles), and Page et al. (2012; diamonds). Red symbols refer to X-ray-selected AGNs, blue to optically selected AGNs, and green to mid-IR-selected AGNs; the SFR measurements are obtained from far-IR observations with Herschel, both via detections (filled symbols) and via stacking (empty symbols); the sample by Venemans et al. (2018) is based on detections from ALMA millimeter data. Results from this work are illustrated as magenta stars for the overall sample at $\langle z \rangle \sim 1.5-2$ (tiny stars), for WISE W4-detected galaxies at $\langle z \rangle \sim 1.5-2$ (small stars), for WISE W4-detected galaxies with W4/r-band luminosity ratio exceeding 100 at $\langle z \rangle \sim 2-3$ (medium-size stars), and for the unique HotDOG falling in the Herschel-ATLAS survey area at $z \sim 3$ (big stars); filled stars refer to far-IR detections with Herschel and unfilled symbols to stacked measurements (see Table 3 for more detailed information on photometric redshift). The thin solid black line is the average relationship obtained according to the AGN variability model (Hickox et al. 2014; Stanley et al. 2015, 2017; see text for details).

Elvis No-SF WDD templates. In fact, the stacked data from WISE and AKARI are fundamental in revealing the bump and the decrease in the mid-IR; similar evidences have also been found in individual hyperluminous QSOs and HotDOGs observed with WISE by Tsai et al. (2015, cf. their Figure 3). Finally, above 20 µm, our stacked SED increases again and behaves similarly to the Richards template in the far-IR.

We stress that our SED has been obtained by stacking an unprecedented number of QSOs with multiband data from the UV to the far-IR band (including SDSS, UKIDSS, WISE, and Herschel). This has enabled us to reveal detailed features of the average AGN SED that were missed by previous studies based on more limited statistics and spectral coverage. In particular, the presence of the bump around 4–6 µm and of the steep decrease in the mid-IR constitute the most relevant finding from our stacked analysis.

A qualitative physical interpretation of these features may be related to the impact of the QSO emission on the surrounding galactic environment. Recall that the stacked signal in our analysis is dominated by individual QSO hosts with far-IR emission undetected in Herschel; the stacked sample includes both galaxies with intrinsic SFR below the Herschel detection limit and galaxies with large-scale star formation already substantially reduced by the energy/momentum feedback from the QSO. Furthermore, a determination of the stellar mass for these sources is fundamental in interpreting the data.

On smaller scales, the dusty torus around the central SMBH, in general responsible for most of the mid-IR emission, may be affected by the QSO feedback. Specifically, the torus may be partly peeled off, the associated warm-dust covering factor may be reduced substantially, and the residual dust component may be instead excited to high temperatures. A partially eroded warm-dust torus with a residual hot-dust component may indeed be at the origin of the mid-IR behavior characterizing our measured SED (see the WDD template in Figure 6). However, a detailed radiative transfer model would be required to test such an hypothesis; this is beyond the scope of the present work and will be performed in a forthcoming paper.

### 4.3. Interpretation in Terms of Galaxy–BH Coevolution Scenario

We now turn to physically interpret our results in terms of an in situ coevolution scenario between the host galaxy and the central SMBH, developed by Lapi et al. (2014, 2017; see also Mancuso et al. 2016, 2017). For such purposes, we focus on the SFR versus $L_{\text{AGN}}$ plane, namely, a diagram where the SFR in the host galaxy is plotted against the bolometric AGN luminosity $L_{\text{AGN}}$; see Figure 7.

The dotted line represents the locus where the bolometric AGN luminosity $L_{\text{AGN}}$ equals the bolometric luminosity $L_{\text{SFR}}$ [erg s$^{-1}$] $\sim 3 \times 10^{43}$ SFR [$M_{\odot}$ yr$^{-1}$] associated with star
formation in the host. In our scenario (see Lapi et al. 2014; Mancuso et al. 2016, for details), galaxies located to the right of such a line have accumulated appreciable stellar masses, host a growing SMBH that has attained considerable luminosities, and so are prone to be strongly affected by the resulting AGN energy/momentum feedback. The dashed lines illustrate three typical evolutionary tracks for individual galaxies based on our coevolution scenario, which correspond to values of SFR ~ 400, 1200, and 4000 $M_{\odot}$ yr$^{-1}$. Note that these tracks refer to the time-average behavior of the AGN luminosity over 100 Myr or more, but do not exclude the presence of variability over much shorter timescales.

During the early stages of a massive galaxy's evolution, the SFR is approximately constant with the galactic age, while the AGN luminosity is smaller than that associated with star formation, but increases much more rapidly; in the SFR versus $L_{\text{AGN}}$ plane, this generates an evolutionary track that is closely parallel to the $L_{\text{AGN}}$ axis. After some $10^8$ yr, the nuclear power progressively increases to values overwhelming that from star formation in the host, which is affected and reduced on a rapid timescale, and soon after even the AGN luminosity decreases because fuel in the central region gets exhausted; in the SFR versus $L_{\text{AGN}}$ plane, this generates a track receding toward the bottom-left portion of the diagram.

The average statistical relationship between SFR and AGN luminosity in the SFR versus $L_{\text{AGN}}$ plane can be determined on the basis of these tracks for individual galaxies by taking into account three additional ingredients: (i) the number density of galaxies with a given value of the SFR, implying that galaxies with higher SFRs are rarer, (ii) the duty cycle, specifying that a galaxy spends different relative time intervals in the different portions of the track, and (iii) the detection threshold for star formation in the host by current far-IR measurements, amounting to SFR ~ $10^2 M_{\odot}$ yr$^{-1}$. The thick solid lines show the resulting average relationship, and the shaded areas render its $1\sigma$ (light gray), $2\sigma$ (gray), and $3\sigma$ (dark gray) variance (see Mancuso et al. 2016 for details on the computation). The rise of the average SFR for increasing values of $L_{\text{AGN}}$ just occurs because, statistically, to achieve a higher AGN luminosity, the BH must reside in a more massive galaxy with a higher initial SFR.

For comparison, we report in the coevolution plane the literature data from Stanley et al. (2015, 2017, circles) at $z \sim 2$, Netzer et al. (2016, squares) at $z \sim 3$, Venemans et al. (2018, pentagons) at $z \sim 6$, Fan et al. (2016, inverse triangles) at $z \sim 3$, Xu et al. (2015, triangles) at $z \sim 2$, and Page et al. (2012, diamonds) at $z \sim 2$; the SFRs have been converted to a Chabrier IMF (see Table 2 in Bernardi et al. 2010, for conversion factors). Red symbols refer to X-ray-selected AGNs, blue to optically selected AGNs, and green to mid-IR-selected AGNs; the SFR measurements are obtained from far-IR observations with Herschel, both from detections (filled symbols) and from stacking techniques (unfilled symbols); the sample by Venemans et al. (2018) is based on detections from ALMA millimeter data. For example, in the upper right portion of the diagram, the mid-IR-selected objects by Fan et al. (2016) are hyperluminous HotDOGs; these are violently star-forming galaxies featuring a high-temperature dust component likely heated up by the powerful energy/momentum feedback from the central AGN.

According to our scenario outlined above, the position on the diagram of these data points can be understood as follows.

Optical data tend to pick up objects close to the peak of AGN luminosity when dust has been partially removed and the SFR is being progressively affected; X-ray data can pick up objects before or after the AGN peak, hence with dusty SFR in the host still sustained or reduced, respectively. Mid-IR data with the current observational limits strike an intermediate course between the former two. On top of that, detections in the far-IR with Herschel tend to probe galaxies when the SFR is still sustained to high values, while stacked measurements tend to probe the epochs after the SFR has been partly affected by AGN feedback. Very interestingly, the ALMA millimeter data by Venemans et al. (2018) in powerful quasars at $z \sim 6$ start to detect star formation in the host even after its suppression has already started/occurred.

Finally, results from this work are illustrated as magenta stars. In order of increasing symbol size, we show the outcome for the overall sample at $\langle z \rangle \sim 1.5-2$ (tiny stars), for WISE W4-detected galaxies at $\langle z \rangle \sim 1.5-2$ (small stars), for WISE W4-detected galaxies with $S_{W4}/S_r$ ratio exceeding 100 at $\langle z \rangle \sim 2-3$ (medium-size stars), and for the unique HotDOG falling in the Herschel-ATLAS survey area (big stars) with $S_{W4}/S_r$ ratio exceeding 1000, which is located at $z \sim 3$; filled stars refer to far-IR detections with Herschel and unfilled symbols to the corresponding stacked measurements.

We can draw the following picture. First, our results follow well the average SFR versus $L_{\text{AGN}}$ relationship within its $1\sigma$ dispersion; actually, we can identify two parallel sequences, one for detected sources that lie above the mean relationship and one for stacked sources that lie below. It is remarkable that, despite our sample being primarily pre-selected in the optical band from the SDSS survey, our additional selection involving mid-IR bands (i.e., requiring detection in the WISE W4 band, or $S_{W4}/S_r \gtrsim 100$, or prominent emission in $S_{W4}/S_r \gtrsim 1000$ as in the HotDOG) constitutes an effective way to move along the average relationship toward increasing AGN luminosities.

Second, in terms of evolutionary track for individual galaxies, the results of this work are located just before (average for far-IR-/(-sub-)millimeter-detected sources) or just after (average for undetected sources) the epoch when the star formation starts to be strongly affected by AGN feedback. Note that, given the rarity of such hyperluminous objects and the relatively small Herschel-ATLAS survey area, it is problematic to pinpoint the stacked measurements corresponding to the HotDOG. Stacey et al. (2018), albeit with noticeable uncertainties, has tentatively detected suppressed star formation in powerful quasar hosts at $z \sim 3$ with Herschel, thanks to gravitational lensing by foreground objects. At much higher redshift $z \sim 6$, Venemans et al. (2018) has found similar suppressed star formation in powerful quasar hosts, thanks to millimeter ALMA observations. In the future, it would be extremely interesting to complement these data with information from X-ray surveys that could allow individual evolutionary tracks to be probed both in the early stages of obscured BH accretion and in the late stage after the onset of feedback.

Finally, it is worth reporting that another popular interpretation of the SFR versus $L_{\text{AGN}}$ diagram calls into play AGN variability (see Aird et al. 2013; Hickox et al. 2014; Stanley et al. 2015, 2017), as inspired by merger-driven models of galaxy formation and numerical simulations (e.g., Novak et al. 2011; Veale et al. 2014; Volonteri et al. 2015; McAlpine et al. 2017). Such a scenario envisages that, when averaged over timescales longer than 100 Myr, the AGN luminosity $\langle L_{\text{AGN}} \rangle$ is
correlated to the SFR to keep a constant ratio of the BH accretion rate to the SFR of 1/3000, as observed in the local universe (e.g., Chen et al. 2013). However, the luminosities $L_{AGN}$ of individual AGNs are allowed to vary on shorter timescales according to an assumed AGN luminosity distribution $dp/d \log L_{AGN} \propto L_{AGN}^{0.2} \exp(-L_{AGN}/100 \langle L_{AGN} \rangle)$, which features a broken power-law form at the faint end and an exponential cutoff. Coupling such a distribution with the statistics of galaxies of a given SFR and redshift yields an average relationship in the SFR versus $L_{AGN}$ plane. We perform the computation with the variability model by using the same SFR galaxy statistics and detection threshold employed above for our scenario, and show in Figure 7 the outcome as a thin solid line. With the present data, it is challenging to differentiate between the in situ coevolution and variability scenarios. A marked difference can be recognized in the statistics and the interpretation of the objects in the bottom right region of the SFR versus $L_{AGN}$ plane, i.e., at high AGN luminosity $L_{AGN} \gtrsim$ some $10^{47}$ erg s⁻¹ and relatively low SFRs $\lesssim$ a few $10^2 M_{\odot}$ yr⁻¹. In the in situ coevolution picture, these are galaxies close to the end of their evolution, when strong AGN activity powered by a BH with mass $\gtrsim 10^9 M_{\odot}$ is removing dust and suppressing the SFR with respect to the high values $\gtrsim 10^3 M_{\odot}$ yr⁻¹ applying to the early dust-enshrouded phase; as such, these systems should feature large stellar masses $M_{\star} \gtrsim 10^{11} M_{\odot}$ already accumulated in their central regions $\lesssim$ a few kiloparsecs. On the contrary, in the AGN variability framework, these galaxies are selected by chance at AGN luminosities substantially higher than the average because of short-time variability; the adopted AGN luminosity distribution will limit their number density and will imply independence from the stellar mass. In the near future, accurate estimates of the number density of the galaxies populating the bottom right region of the SFR versus $L_{AGN}$ plane, and precise measurements of the stellar masses in these systems will contribute to distinguishing between the two competing scenarios; these observations will likely become feasible, even at substantial redshifts, with the advent of the *James Webb Space Telescope*.

## 4.4. Future Outlook

We now move on to discuss natural extensions of the analysis presented in this paper.

A promising avenue that has recently emerged for investigating AGN activity—and in particular the feedback energetics—is its thermal Sunyaev–Zel'dovich (tSZ) signature (Sunyaev & Zeldovich 1970, 1972) at millimeter wavelengths. The energy released by the AGN to the surrounding diffuse baryons affects their temperature and distribution. As a consequence, the heated electron gas can in principle Compton-upscatter the CMB photons, leaving an imprint in the CMB blackbody spectrum by tilting it to higher frequencies (Voit 1994; Natarajan & Sigurdsson 1999; Lapi et al. 2003; Chatterjee & Kosowsky 2007; Scannapieco et al. 2008). However, the AGN feedback-induced tSZ is expected to be too feeble to be detected at the single source level; therefore, one must resort to statistical techniques such as cross-correlation and stacking. The first tentative detection of tSZ due to the AGN feedback effect was reported by Chatterjee et al. (2010), while subsequent works by Ruan et al. (2015), Crichton et al. (2016), Verdier et al. (2016), Spacek et al. (2016, 2017), and Soergel et al. (2017) found evidence of this effect at different

significance levels. However, there is no general consensus as the inferred thermal energies differ substantially among the various studies. One of the main obstacle in this analysis is represented by the strong dust signal that dominates the SED at $\nu \gtrsim 100$ GHz (Crichton et al. 2016; Verdier et al. 2016; Soergel et al. 2017) and creates a degeneracy with the tSZ signal. An additional complication involves distinguishing between the AGN-induced signal and that of the virialized gas in the host dark matter halo; see, e.g., Dutta Chowdhury & Chatterjee (2017).

Following Soergel et al. (2017), in this paper we have performed the first steps in addressing this aspect by stacking maps at *Planck* wavelengths (above and below the SZ null frequency at $\simeq 217$ GHz). One of the outcomes of this work also informs us about the best approach and most suited AGN samples for investigating the AGN feedback-induced tSZ imprint, which is ultimately limited by the dust emission. A possible way forward to look for this SZ signature would be to use high-resolution and high-sensitivity ground-based telescopes with multifrequency coverage.

## 5. Conclusions

In this paper, we have provided a multiwavelength picture of SDSS optically selected QSOs in the redshift range $1 \le z \le 5$ from the near-UV to the submillimeter bands. To this end, we have combined diverse photometric information from the SDSS, UKIDSS, and WISE surveys, as well as performing a stacking analysis on Herschel, AKARI, and *Planck* maps.

The shape of the reconstructed mean SED is similar to other templates in the UV-optical part of the spectrum, while it features a little excess at $1-2 \ \mu$m (likely due to old stars and heated dust), a prominent bump around $4-6 \ \mu$m followed by a decrease out to $\sim 20 \ \mu$m and a subsequent far-IR bump comparable to that measured by other studies. By decomposing the SED into an AGN, starburst, and torus components, we have estimated the AGN luminosity and SFR, finding typical $L_{AGN} \sim 10^{46}-10^{47}$ erg s⁻¹ and SFR $\sim 50-1000 M_{\odot}$ yr⁻¹.

When interpreted in the context of the in situ coevolution scenario for galaxies and hosted SMBHs presented by Lapi et al. (2014), our results suggest that the mid-IR-based selection criteria developed in this work (detection in the WISE W4 band and $S_{W4}/S_r > 100$) have been proven effective in moving along the SFR versus $L_{AGN}$ diagram toward higher AGN luminosities; on the other hand, the far-IR/submillimeter detection criterion, if (not) enforced, allows us to probe the epoch just (after) before the onset of feedback from the central AGN.

Natural extensions of the current work involve complementing these data sets with X-ray observations, in order to shed light on the early stages of the BH accretion and the evolution after the onset of AGN feedback, as well as with millimeter and radio information to investigate the alleged AGN feedback-induced Sunyaev–Zel’dovich emission and constrain the non-gravitational heating of the intergalactic medium.

We thank the referee for constructive comments that helped in improving the paper. We are indebted to Luigi Danese for his support during the early stages of this work. We acknowledge very fruitful discussions with S. Mateos and F. J. Carrera. We also thank Simone Aiola, Srinivasan Ragunathan, Christian L. Reichardt, and Reinier Janssen for providing comments on an earlier version of the manuscript. F.B. acknowledges support from an Australian Research Council Future
Fellowship (FT150100074). G.F. acknowledges the support of the CNES postdoctoral program. J.G.N. acknowledges financial support from the I+D 2015 project AYA2015-65887-P (MINECO/FEDER) and from the Spanish MINECO for a "Ramon y Cajal" fellowship (RYC-2013-13256). A.L. acknowledges the RADIOFOREGROUNDS grant (COM-PET-05-2015, agreement No. 687312) of the European Union Horizon 2020 research and innovation program and the MIUR grant "Finanziamento annuale individuale attività base di ricerca." R.G. acknowledges support from the agreements ASI-INAF n.2017-14-H.O and ASI-INAF I/037/12/0.

Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III website is http://www.sdss3.org/. SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofísica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.

*Software:* Astropy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018), IPython (Pérez & Granger 2007), matplotlib (Hunter 2007), scipy (Jones et al. 2001), HEALPix (Górski et al. 2005).

## Appendix A
### Testing the Far-IR Subsample Representativity

In this appendix, we assess to what extent the fluxes of sources lying in the H-ATLAS footprint (i.e., the H-ATLAS

subsample) are representative of the entire QSO sample. Adopting representative subsamples is in fact a crucial point when studying the demographics of a given sources population.

As discussed in Section 2.5, when no submillimeter detection criterion is enforced, the average photometry for the optical and mid-IR part of the SED is calculated over the whole QSO sample, while the far-IR/submillimeter portion of the spectrum is estimated using only sources included in the H-ATLAS footprint. The main reason for this strategy is simply to increase the sample size statistics, an issue that mostly affects low-resolution and relatively high detection-limit experiments such as AKARI and Planck. However, one might wonder whether the fluxes measured for the H-ATLAS subsample are representative of the full sample and therefore provide a reliable estimate of the average flux. In order to test this assumption, we proceed in the following way. From the main QSO sample,<fnref n="19" /> we randomly draw $N_{\text{MC}} = 100$ sub-samples containing a number of sources comparable to that of the H-ATLAS one (about 7000 sources, or 3% of the whole QSO sample with SDSS fluxes). Then, from each of these 100 random subsamples, we estimate the average flux $\bar{S}_{\nu}$ distribution in each photometric band and compare it to the average of fluxes measured first in the H-ATLAS subsample and second over the whole SDSS footprint. The results of this test are shown in Figure 8. In each panel, the gray line shows the distribution of the mean fluxes across the random subsamples at a given photometric band, while the solid red and dashed black lines denote the average flux measured for the whole QSO sample and the H-ATLAS subsample, respectively. The darker and lighter shaded regions indicate the confidence levels equivalent to the $\pm 1\sigma$ and $\pm 2\sigma$ of a Gaussian distribution for the average flux distribution.<fnref n="20" /> The average fluxes measured within the H-ATLAS fields do not appear to be outliers of the distribution, as they fall within the $\pm 1\sigma$ band. The mean flux in 2 out of the 13 photometric bands is within $\pm 2\sigma$, but approximately one-third of the measurements are expected to be scattered more than $1\sigma$. Hence, we conclude that we can safely consider the far-IR/submillimeter average fluxes measured in the H-ATLAS fields to be representative of the whole QSO sample.

<footnote n="19">¹⁹ Note that, for the purposes of this check, we consider all of the sources within the AllQSO sample in a single redshift bin $1 \le z \le 5$.</footnote>

<footnote n="20">²⁰ Evaluated as the (13.6, 86.4) and (2.1, 97.9) percentiles of the distribution respectively.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Assessment of the H-ATLAS subsample representativity. In each panel, the gray line shows a kernel density estimate of the distribution of the mean fluxes in a certain photometric band estimated from a set of 100 random subsamples. The darker and lighter shaded areas indicate the ±1σ and ±2σ of the average flux distribution, evaluated as the (13.6, 86.4) and (2.1, 97.9) percentiles of the distribution, respectively. The solid red and dashed black lines denote the average flux measured for whole QSO sample and the H-ATLAS subsample, respectively.</figcaption></figure>
## Appendix B
SED Templates

In this appendix, we report the best-fit SEDs for the main QSO sample in the three redshift bins between $1 \le z \le 5$. The tabulated values as a function of the rest-frame wavelength are presented in Table 4.

<table><caption>**Table 4**
Tabulated Values of the Overall Average SED for the Baseline Optical-selected QSO Sample as a Function of Redshift and in Terms of the Rest-frame Wavelength</caption><thead><tr><th rowspan="2">log<sub>10</sub> λ (μm)</th><th colspan="3">S<sub>ν</sub> (mJy)</th></tr><tr><th>1 &lt; z &lt; 2.15</th><th>2.15 &lt; z &lt; 2.5</th><th>2.5 &lt; z &lt; 5</th></tr></thead><tbody><tr><td>-1.0</td><td>1.67 × 10<sup>-2</sup></td><td>1.33 × 10<sup>-2</sup></td><td>1.49 × 10<sup>-2</sup></td></tr><tr><td>-0.9</td><td>2.25 × 10<sup>-2</sup></td><td>1.79 × 10<sup>-2</sup></td><td>2.00 × 10<sup>-2</sup></td></tr><tr><td>-0.8</td><td>2.80 × 10<sup>-2</sup></td><td>2.23 × 10<sup>-2</sup></td><td>2.50 × 10<sup>-2</sup></td></tr><tr><td>-0.7</td><td>3.38 × 10<sup>-2</sup></td><td>2.69 × 10<sup>-2</sup></td><td>3.01 × 10<sup>-2</sup></td></tr><tr><td>-0.6</td><td>4.12 × 10<sup>-2</sup></td><td>3.28 × 10<sup>-2</sup></td><td>3.66 × 10<sup>-2</sup></td></tr><tr><td>-0.5</td><td>4.40 × 10<sup>-2</sup></td><td>3.50 × 10<sup>-2</sup></td><td>3.92 × 10<sup>-2</sup></td></tr><tr><td>-0.4</td><td>4.05 × 10<sup>-2</sup></td><td>3.21 × 10<sup>-2</sup></td><td>3.60 × 10<sup>-2</sup></td></tr><tr><td>-0.3</td><td>3.78 × 10<sup>-2</sup></td><td>3.03 × 10<sup>-2</sup></td><td>3.40 × 10<sup>-2</sup></td></tr><tr><td>-0.2</td><td>4.23 × 10<sup>-2</sup></td><td>3.49 × 10<sup>-2</sup></td><td>3.96 × 10<sup>-2</sup></td></tr><tr><td>-0.1</td><td>5.40 × 10<sup>-2</sup></td><td>4.56 × 10<sup>-2</sup></td><td>5.29 × 10<sup>-2</sup></td></tr><tr><td>0.0</td><td>7.85 × 10<sup>-2</sup></td><td>6.29 × 10<sup>-2</sup></td><td>7.45 × 10<sup>-2</sup></td></tr><tr><td>0.1</td><td>1.18 × 10<sup>-1</sup></td><td>8.35 × 10<sup>-2</sup></td><td>9.88 × 10<sup>-2</sup></td></tr><tr><td>0.2</td><td>1.76 × 10<sup>-1</sup></td><td>1.09 × 10<sup>-1</sup></td><td>1.21 × 10<sup>-1</sup></td></tr><tr><td>0.3</td><td>2.64 × 10<sup>-1</sup></td><td>1.68 × 10<sup>-1</sup></td><td>1.61 × 10<sup>-1</sup></td></tr><tr><td>0.4</td><td>3.73 × 10<sup>-1</sup></td><td>2.83 × 10<sup>-1</sup></td><td>2.68 × 10<sup>-1</sup></td></tr><tr><td>0.5</td><td>4.70 × 10<sup>-1</sup></td><td>4.27 × 10<sup>-1</sup></td><td>4.82 × 10<sup>-1</sup></td></tr><tr><td>0.6</td><td>6.02 × 10<sup>-1</sup></td><td>5.61 × 10<sup>-1</sup></td><td>7.97 × 10<sup>-1</sup></td></tr><tr><td>0.7</td><td>9.05 × 10<sup>-1</sup></td><td>7.45 × 10<sup>-1</sup></td><td>1.13 × 10<sup>0</sup></td></tr><tr><td>0.8</td><td>1.35 × 10<sup>0</sup></td><td>1.00 × 10<sup>0</sup></td><td>1.30 × 10<sup>0</sup></td></tr><tr><td>0.9</td><td>1.69 × 10<sup>0</sup></td><td>1.23 × 10<sup>0</sup></td><td>1.23 × 10<sup>0</sup></td></tr><tr><td>1.0</td><td>1.70 × 10<sup>0</sup></td><td>1.27 × 10<sup>0</sup></td><td>9.86 × 10<sup>-1</sup></td></tr><tr><td>1.1</td><td>1.42 × 10<sup>0</sup></td><td>1.10 × 10<sup>0</sup></td><td>6.90 × 10<sup>-1</sup></td></tr><tr><td>1.2</td><td>1.03 × 10<sup>0</sup></td><td>8.18 × 10<sup>-1</sup></td><td>4.37 × 10<sup>-1</sup></td></tr><tr><td>1.3</td><td>6.97 × 10<sup>-1</sup></td><td>5.60 × 10<sup>-1</sup></td><td>2.77 × 10<sup>-1</sup></td></tr><tr><td>1.4</td><td>6.40 × 10<sup>-1</sup></td><td>4.85 × 10<sup>-1</sup></td><td>3.05 × 10<sup>-1</sup></td></tr><tr><td>1.5</td><td>1.19 × 10<sup>0</sup></td><td>8.06 × 10<sup>-1</sup></td><td>7.18 × 10<sup>-1</sup></td></tr><tr><td>1.6</td><td>2.41 × 10<sup>0</sup></td><td>1.57 × 10<sup>0</sup></td><td>1.56 × 10<sup>0</sup></td></tr><tr><td>1.7</td><td>3.83 × 10<sup>0</sup></td><td>2.57 × 10<sup>0</sup></td><td>2.52 × 10<sup>0</sup></td></tr><tr><td>1.8</td><td>5.01 × 10<sup>0</sup></td><td>3.22 × 10<sup>0</sup></td><td>3.31 × 10<sup>0</sup></td></tr><tr><td>1.9</td><td>5.84 × 10<sup>0</sup></td><td>3.75 × 10<sup>0</sup></td><td>3.87 × 10<sup>0</sup></td></tr><tr><td>2.0</td><td>6.11 × 10<sup>0</sup></td><td>3.92 × 10<sup>0</sup></td><td>4.04 × 10<sup>0</sup></td></tr><tr><td>2.1</td><td>5.53 × 10<sup>0</sup></td><td>3.56 × 10<sup>0</sup></td><td>3.67 × 10<sup>0</sup></td></tr><tr><td>2.2</td><td>4.29 × 10<sup>0</sup></td><td>2.76 × 10<sup>0</sup></td><td>2.84 × 10<sup>0</sup></td></tr><tr><td>2.3</td><td>2.88 × 10<sup>0</sup></td><td>1.85 × 10<sup>0</sup></td><td>1.91 × 10<sup>0</sup></td></tr><tr><td>2.4</td><td>1.71 × 10<sup>0</sup></td><td>1.10 × 10<sup>0</sup></td><td>1.14 × 10<sup>0</sup></td></tr></tbody></table>

## Appendix C
SED-fitting Methodology Test

When dealing with stacked sources (in our case their spectroscopic redshift is known but not their individual flux densities), there is always an assumption to make in order to extract useful information from the analyzed data. The methodology used for the SED fitting in this work assumes that the different adopted SED templates are representative of the whole population inside the redshift bin. Therefore, they are simply used to estimate the main redshift of the sample and the flux density normalization to eventually derive the intrinsic luminosities (see Section 3.4 for more details).

However, we note the caveat that high-redshift objects ($z > 1$) can suffer strong K-corrections. One may wonder whether our results are stable when taking into account the individual redshifts of the sources, so as to avoid mixing information from different wavelengths. For example, this is the approach followed by Crichton et al. (2016) in the IR band. These authors considered a graybody emission with emissivity index $\beta$ and the dust temperature as fitting parameters; on the other hand, they were forced to assume the same luminosity for all of the objects inside the redshift bin. Then they estimated the stacked redshifted SEDs following the same procedure applied to the data and used it to get the best SED normalization (in terms of bolometric luminosity).

In order to check the robustness of our results, we perform a simple test comparing both approaches. For the model SED, we use the averaged SED estimated for the AllQSO sample in the first redshift bin (the most critical one since the redshift distribution there is rather flat). We use the full SED in order to show the main differences between the methodologies, which are only visible in the most prominent SED features.

We consider two simulated populations of 1000 QSOs with two different redshift distributions (see inset in Figure 9): a uniform distribution (red histogram) and a Gaussian distribution (blue histogram). Taking into account the different redshifts for each QSO, we then estimate the main redshifted SED for both redshift distributions. We finally compare them with the initial SED at the mean redshift $z = 1.6$.

As can be observed in Figure 9, the three estimated SEDs are almost identical. The main effect of considering a redshift distribution is a negligible smoothing, mostly observed in the valley feature around 60 µm at $z = 1.6$. Therefore, we can safely conclude that our results and conclusions are robust and independent of the adopted SED-fitting methodology.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Results of the simple test using two different methodologies to estimate the average SED (see the text for more details): main sample of QSOs in the first redshift bin (1.0 < z < 2.15) derived at the mean redshift (z = 1.6, black solid line); rest-frame averaged SED considering uniform (dashed red line) and Gaussian (dotted–dashed blue line) redshift distributions inside the redshift bin.</figcaption></figure>

## ORCID iDs

Federico Bianchini https://orcid.org/0000-0003-4847-3483

Giulio Fabbian https://orcid.org/0000-0002-3255-4695

Andrea Lapi https://orcid.org/0000-0002-4882-1735

Joaquin Gonzalez-Nuevo https://orcid.org/0000-0003-1354-6822

Roberto Gilli https://orcid.org/0000-0001-8121-6177

## References

* Aird, J., Coil, A. L., Moustakas, J., et al. 2013, ApJ, 775, 41

* Aird, J., Nandra, K., Laird, E. S., et al. 2010, MNRAS, 401, 2531

* Alexander, D. M., & Hickox, R. C. 2012, NewAR, 56, 93

* Antonucci, R. 1993, ARA&A, 31, 473

* Becker, R. H., White, R. L., & Helfand, D. J. 1995, ApJ, 450, 559

* Beelen, A., Cox, P., Benford, D. J., et al. 2006, ApJ, 642, 694

* Bernardi, M., Shankar, F., Hyde, J. B., et al. 2010, MNRAS, 404, 2087

* Bonfield, D. G., Jarvis, M. J., Hardcastle, M. J., et al. 2011, MNRAS, 416, 13

* Bourne, N., Dunne, L., Maddox, S. J., et al. 2016, MNRAS, 462, 1714

* Chabrier, G. 2003, PASP, 115, 763

* Chatterjee, S., Ho, S., Newman, J. A., & Kosowsky, A. 2010, ApJ, 720, 299

* Chatterjee, S., & Kosowsky, A. 2007, ApJ, 661, L113

* Chen, C.-T. J., Hickox, R. C., Alberts, S., et al. 2013, ApJ, 773, 3

* Cicone, C., Maiolino, R., Sturm, E., et al. 2014, A&A, 562, A21

* Crichton, D., Gralla, M. B., Hall, K., et al. 2016, MNRAS, 458, 1478

* Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11

* Cutri, R. M., Wright, E. L., Conrow, T., et al. 2013, yCat, 2328, 0

* Del Moro, A., Alexander, D. M., Mullaney, J. R., et al. 2013, A&A, 549, A59

* Delvecchio, I., Gruppioni, C., Pozzi, F., et al. 2014, MNRAS, 439, 2736

* Di Matteo, T., Springel, V., & Hernquist, L. 2005, Natur, 433, 604

* Doi, Y., Takita, S., Ootsubo, T., et al. 2015, PASJ, 67, 50

* Dutta Chowdhury, D., & Chatterjee, S. 2017, ApJ, 839, 34

* Eales, S., Dunne, L., Clements, D., et al. 2010, PASP, 122, 499

* Eifler, T., Krause, E., Dodelson, S., et al. 2015, MNRAS, 454, 2451

* Elvis, M., Wilkes, B. J., McDowell, J. C., et al. 1994, ApJS, 95, 1

* Fabian, A. C. 2012, ARA&A, 50, 455

* Fan, L., Han, Y., Nikutta, R., Drouart, G., & Knudsen, K. K. 2016, ApJ, 823, 107

* Feigelson, E. D., & Nelson, P. I. 1985, ApJ, 293, 192

* Ferrarese, L., & Merritt, D. 2000, ApJL, 539, L9

* Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJL, 539, L13

* Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

* Granato, G. L., De Zotti, G., Silva, L., Bressan, A., & Danese, L. 2004, ApJ, 600, 580

* Gürkan, G., Hardcastle, M. J., Jarvis, M. J., et al. 2015, MNRAS, 452, 3776

* Haehnelt, M. G., & Tegmark, M. 1996, MNRAS, 279, 545

* Hall, K. R., Crichton, D., Marriage, T., Zakamska, N. L., & Mandelbaum, R. 2018, MNRAS, 480, 149

* Harris, K., Farrah, D., Schulz, B., et al. 2016, MNRAS, 457, 4179

* Hatziminaoglou, E., Omont, A., Stevens, J. A., et al. 2010, A&A, 518, L33

* Hickox, R. C., Mullaney, J. R., Alexander, D. M., et al. 2014, ApJ, 782, 9

* Hunter, J. D. 2007, CSE, 9, 90

* Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python, www.scipy.org/

* Kalfountzou, E., Stevens, J. A., Jarvis, M. J., et al. 2014, MNRAS, 442, 1181

* Kennicutt, R. C., Jr. 1998, ARA&A, 36, 189

* King, A. 2003, ApJL, 596, L27

* Kormendy, J., & Richstone, D. 1995, ARA&A, 33, 581

* Lani, C., Netzer, H., & Lutz, D. 2017, MNRAS, 471, 59

* Lapi, A., Cavaliere, A., & De Zotti, G. 2003, ApJL, 597, L93

* Lapi, A., Mancuso, C., Bressan, A., & Danese, L. 2017, ApJ, 847, 13

* Lapi, A., Raimundo, S., Aversa, R., et al. 2014, ApJ, 782, 69

* Lawrence, A., Warren, S. J., Almaini, O., et al. 2007, MNRAS, 379, 1599

* Leipski, C., Meisenheimer, K., Walter, F., et al. 2013, ApJ, 772, 103

* Lyu, J., Rieke, G. H., & Shi, Y. 2017, ApJ, 835, 257

* Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

* Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285

* Mancuso, C., Lapi, A., Prandoni, I., et al. 2017, ApJ, 842, 95

* Mancuso, C., Lapi, A., Shi, J., et al. 2016, ApJ, 833, 152

* Marsden, G., Ade, P. A. R., Bock, J. J., et al. 2009, ApJ, 707, 1729

* McAlpine, S., Bower, R. G., Harrison, C. M., et al. 2017, MNRAS, 468, 3395

* Mullaney, J. R., Alexander, D. M., Goulding, A. D., & Hickox, R. C. 2011, MNRAS, 414, 1082

* Natarajan, P., & Sigurdsson, S. 1999, MNRAS, 302, 288

* Netzer, H., Lani, C., Nordon, R., et al. 2016, ApJ, 819, 123

* Novak, G. S., Ostriker, J. P., & Ciotti, L. 2011, ApJ, 737, 26
Page, M. J., Brindle, C., Talavera, A., et al. 2012, MNRAS, 426, 903

Pâris, I., Petitjean, P., Ross, N. P., et al. 2017, A&A, 597, A79

Pearson, E. A., Eales, S., Dunne, L., et al. 2013, MNRAS, 435, 2753

Pérez, F., & Granger, B. E. 2007, CSE, 9, 21

Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1

Planck Collaboration, Adam, R., Ade, P. A. R., et al. 2016a, A&A, 594, A1

Planck Collaboration, Adam, R., Ade, P. A. R., et al. 2016b, A&A, 594, A10

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014b, A&A, 571, A9

Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2018, arXiv:1807.06209

Planck Collaboration, Aghanim, N., Armitage-Caplan, C., et al. 2014a, A&A, 571, A2

Planck Collaboration, Aghanim, N., Ashdown, M., et al. 2016c, A&A, 596, A109

Podigachoski, P., Rocca-Volmerange, B., Barthel, P., Drouart, G., & Fioc, M. 2016, MNRAS, 462, 4183

Price-Whelan, A. M., Sip'ocz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123

Rees, M. J. 1984, ARA&A, 22, 471

Remazeilles, M., Delabrouille, J., & Cardoso, J.-F. 2011, MNRAS, 418, 467

Richards, G. T., Lacy, M., Storrie-Lombardi, L. J., et al. 2006, ApJS, 166, 470

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Rosario, D. J., Trakhtenbrot, B., Lutz, D., et al. 2013, A&A, 560, A72

Ross, N. P., Myers, A. D., Sheldon, E. S., et al. 2012, ApJS, 199, 3

Ruan, J. J., McQuinn, M., & Anderson, S. F. 2015, ApJ, 802, 135

Scannapieco, E., Thacker, R. J., & Couchman, H. M. P. 2008, ApJ, 678, 674

Schäfer, P. M., Pfrommer, C., Hell, R. M., & Bartelmann, M. 2006, MNRAS, 370, 1713

Schneider, D. P., Richards, G. T., Hall, P. B., et al. 2010, AJ, 139, 2360

Semboloni, E., Hoekstra, H., Schaye, J., van Daalen, M. P., & McCarthy, I. G. 2011, MNRAS, 417, 2020

Soergel, B., Giannantonio, T., Efstathiou, G., Puchwein, E., & Sijacki, D. 2017, MNRAS, 468, 577

Soltan, A. 1982, MNRAS, 200, 115

Spacek, A., Scannapieco, E., Cohen, S., Joshi, B., & Mauskopf, P. 2016, ApJ, 819, 128

Spacek, A., Scannapieco, E., Cohen, S., Joshi, B., & Mauskopf, P. 2017, ApJ, 834, 102

Stacey, H. R., McKean, J. P., Robertson, N. C., et al. 2018, MNRAS, 476, 5075

Stanley, F., Alexander, D. M., Harrison, C. M., et al. 2017, MNRAS, 472, 2221

Stanley, F., Harrison, C. M., Alexander, D. M., et al. 2015, MNRAS, 453, 591

Sturm, E., González-Alfonso, E., Veilleux, S., et al. 2011, ApJL, 733, L16

Sunyaev, R. A., & Zeldovich, Y. B. 1970, CoASP, 2, 66

Sunyaev, R. A., & Zeldovich, Y. B. 1972, CoASP, 4, 173

Symeonidis, M., Giblin, B. M., Page, M. J., et al. 2016, MNRAS, 459, 257

Tristram, M., Macías-Pérez, J. F., Renault, C., & Santos, D. 2005, MNRAS, 358, 833

Tsai, C.-W., Eisenhardt, P. R. M., Wu, J., et al. 2015, ApJ, 805, 90

Urry, C. M., & Padovani, P. 1995, PASP, 107, 803

Valiante, E., Smith, M. W. L., Eales, S., et al. 2016, MNRAS, 462, 3146

van Daalen, M. P., Schaye, J., Booth, C. M., & Dalla Vecchia, C. 2011, MNRAS, 415, 3649

Veale, M., White, M., & Conroy, C. 2014, MNRAS, 445, 1144

Venemans, B., Decarli, R., Walter, F., et al. 2018, ApJ, 886, 159

Verdier, L., Melin, J.-B., Bartlett, J. G., et al. 2016, A&A, 588, A61

Voit, G. M. 1994, ApJL, 432, L19

Volonteri, M., Capelo, P. R., Netzer, H., et al. 2015, MNRAS, 452, L6

Xu, L., Rieke, G. H., Egami, E., et al. 2015, ApJS, 219, 18