SPECTRAL ENERGY DISTRIBUTIONS OF QSOs AT z > 5: COMMON ACTIVE GALACTIC NUCLEUS-HEATED DUST AND OCCASIONALLY STRONG STAR-FORMATION

The parent sample for this study consisted of all quasars with redshift z > 5 that were known at the time of submission of the original Herschel proposal (early 2007). Due to various factors (e.g., revised sensitivity estimates after the launch of Herschel, sparse supplemental data coverage, revised redshifts, uncertain identifications), a small number of sources was subsequently removed from the target list. The final sample includes 69 quasars at z > 5, all of which have been observed with Herschel in five bands. For 68 of them, we also present Spitzer photometry in 5 bands.

Most of the quasars in our sample come from the Sloan Digital Sky Survey (SDSS), either from the main survey or from the deeper Stripe 82. A small complement of objects consists of serendipitously discovered high-redshift quasars (Sharp et al. 2001; Romani et al. 2004; Mahabal et al. 2005; McGreer et al. 2006). The final sample is presented in Table 1, along with an observation log for the Herschel data. In this table we give the full name of the source. For the following tables, figures and in the text we use abbreviated source names in the format of Jhhmm ± ddmm.

The UV continuum brightness of high-redshift quasars is typically indicated by their monochromatic flux at a rest frame wavelength of 1450 Å, often expressed in terms of apparent AB magnitudes. We have compiled these values for our quasars from the literature and report them in Table 1. Two objects (J0841+2905 and J2245+0024) had only their absolute magnitude at 1450 Å (rest frame) given (Sharp et al. 2001; Goto 2006). For these cases we calculated the apparent magnitudes using the world models cited in the respective papers.

For objects that did not have mag(1450 Å) available in the literature, we retrieved the spectrum from the SDSS data base, corrected it for galactic foreground extinction using the map of Schlegel et al. (1998) and determined mag(1450 Å) from the corrected spectrum following the procedure of Fan et al. (2004). This approach has been adopted for 31 objects with z ⩽ 5.41 (see Table 1).

Where no values for the 1450 Å flux were provided in the literature and no spectra where available in electronic form (six sources, see Table 1), we scaled a redshifted version of the SDSS quasar template spectrum (Vanden Berk et al. 2001) to match the (extinction corrected) z-band magnitude (taking into account the filter curve). From the redshifted and scaled template spectrum we then determined mag(1450 Å) as in Fan et al. (2004).

We also compiled z-band and y-band photometry for the majority of the sample (see Table 2). For most of the quasars, the z-band photometry is taken from SDSS or from the discovery papers, which sometimes presented deeper observations.

The y-band photometry was mainly provided by Pan-STARRS (Kaiser et al. 2010), complemented by data from UKIDSS (Lawrence et al. 2007). For objects where the NIR photometry (see below) was taken from UKIDSS, we also used the y-band flux from this survey, for consistency. In most cases where y-band data exist from both surveys, they agree within the combined errors.

An important source of photometry in the J, H, and/or K bands were the discovery papers (or follow-up work on those). In the majority of cases, magnitudes were given in a Vega-based system and were obtained with a multitude of instruments across our sample. We here consistently use the values given in Hewett et al. (2006) to convert all the Vega-based magnitudes into the AB system. The NIR photometry from the literature was complemented by photometry from UKIDSS for a sizable fraction of our sample. For additional nine objects we obtained J-band photometry using Omega2000 at the 3.5 m telescope of the Calar Alto observatory. For the data reduction and photometry we followed standard procedures. Magnitudes and the corresponding references are reported in Table 2.

Mid-infrared imaging from Spitzer exists for all Herschel targets, with the exception of J2054−0005. These data consist of observations at 3.6, 4.5, 5.8, and 8.0 μm with IRAC (Fazio et al. 2004) as well as at 24 μm with MIPS (Rieke et al. 2004). A small number of objects was also observed at 16 μm with the peak-up array of the InfraRed Spectrograph (IRS) on board Spitzer (Houck et al. 2004).

The Spitzer data were processed in a standard manner using procedures within the MOPEX software package provided by the Spitzer Science Center (SSC). The resulting maps from IRAC and MIPS are presented in the Appendix in Figure 15. Aperture photometry on the final images was performed in IDL. In most cases we used apertures with a radius of 3.6, 5.4, and 7 arcsec in IRAC, IRS, and MIPS, respectively. For some objects the aperture size was reduced to avoid contamination from nearby objects. Appropriate aperture corrections were taken from the respective instrument handbooks (also available from the SSC Web site).

Errors on the photometry were determined by measuring the fluxes in 500 apertures (with sizes identical to the science target aperture) which were randomly placed on source-free regions of the background, avoiding area of low coverage. The distribution of these 500 fluxes was fit by a Gaussian. The sigma of this Gaussian was taken as the 1σ uncertainty on the photometry. The measured fluxes and uncertainties are presented in Table 5 for IRAC and MIPS. The additional IRS photometry for a small subset of objects is presented in Table 3. We note that some of the Spitzer data have been published previously (e.g., Hines et al. 2006; Jiang et al. 2006, 2010) and in these cases our photometry is consistent with the earlier results.

For a few objects, direct aperture photometry (even with smaller apertures) was difficult to obtain due to severe blending with neighboring sources. In such cases we used the point source extraction tool APEX in the MOPEX software package to subtract the confusing source from the science image. We then performed aperture photometry as described above on the residual image (i.e., where the confusing source has been removed) for consistency with the rest of the sample.

2.6.1. PACS

2.6.2. SPIRE

2.6.3. Millimeter Regime

Ten quasars in our sample have been detected in at least four of the five Herschel bands (Table 5). In combination with their Spitzer fluxes and using supplemental NIR data, the combined photometry provides SEDs covering the rest frame wavelengths from 0.1 to ∼80 μm. To assess some basic physical properties of these objects, we perform SED fitting, following the approach presented in Leipski et al. (2013). To summarize briefly, the SEDs are fitted with four components: a power-law in the UV/optical mainly representing emission from the accretion disk, a blackbody from hot (∼1200 K) dust, a torus model from the library of Hönig & Kishimoto (2010), and an additional cool dust component in the form of a modified blackbody (β = 1.6). We illustrate this approach and the arrangement of the fitted components schematically in Figure 2.

Zoom In Zoom Out Reset image size

In Leipski et al. (2013), we presented 5 of the 10 FIR-detected quasars, all of which had millimeter detections. The five additional sources presented here do not have millimeter detections and submillimeter/millimeter upper limits exist only for two of the five newly presented objects. In the case of J1204−0021 (Carilli et al. 2001; Priddey et al. 2003), those data points can be used to provide additional constraints on the temperature which is consequently treated as a free parameter. For J1602+4228 (Wang et al. 2008a) the 250 GHz upper limit does not strongly constrain the temperature of the fitted FIR component and we fix T_FIR to a value of 47 K for this object (Beelen et al. 2006; Wang et al. 2007; Leipski et al. 2013). For the remaining three objects, the temperature of the FIR component was also fixed to 47 K. Due to the lack of millimeter data which would help to anchor the Rayleigh–Jeans tail of the FIR component, the fits would otherwise predict artificially increased dust temperatures (Leipski et al. 2013).

The rest frame UV/optical and infrared SEDs of these 10 objects can be fitted well with a combination of these 4 components. The best fitting model combinations are shown in Figure 3 and Table 6 summarizes some basic properties determined from the fitting. Using these fits we also determine the relative contributions of the different components to the total infrared SED. For this we combine the dust component in the NIR and the torus model, both of which are likely to be powered by the AGN. We compare this AGN related emission to the additional FIR component and show their relative contributions to the total infrared emission as a function of wavelength in Figure 4. We see that in the presence of luminous FIR emission (L_FIR ∼ 10¹³ L_☉), this component dominates the total infrared SED at rest frame wavelengths above ∼50 μm for all 10 objects. This means that in such cases of strong FIR/submillimeter emission, rest frame wavelengths ≳50 μm isolate the additional FIR component without the need for full SED fits (at least in our modeling approach). The possible heating source for the additional FIR component (AGN versus star formation) is further discussed in Section 4.4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also extend a similar SED fitting approach to objects with fewer Herschel detections. In cases where two PACS detections are available (nine sources), these data provide sufficient constraints for the torus model, while the upper limits in the SPIRE bands (and in the millimeter where available; see Table 4) limit the contribution of the additional FIR component (fixed to a temperature of 47 K). These fits are presented in Figure 5 and some basic properties derived from the fitted components are presented in Table 6. From this table we use the UV/optical luminosity and the AGN-dominated dust luminosity to show that the ratio of the AGN-dominated dust-to-accretion disk emission decreases with increasing UV/optical luminosity (Figure 6). This behavior may reflect the increase of the dust sublimation radius for more luminous UV/optical continuum emitters (e.g., Barvainis 1987) which, under the assumption of a constant scale height, is often explained in terms of a decreasing dust covering factor with increasing luminosity in the context of the so-called receding torus model (Lawrence 1991).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The measured FIR fluxes for our 10 FIR-detected objects fall only moderately above the 3σ confusion noise limit (Table 5). Thus, the photometric upper limits for the nine FIR non-detections (i.e., only detected in PACS) yield upper limits on L_FIR that do not differ significantly from the detection on an individual basis (Table 6). Further constraints on the average FIR properties of the PACS-only sources are provided by a stacking analysis as presented in Section 4.4.

For two-thirds of the sample, the upper limits in the Herschel observations do not provide strong constraints to MIR or FIR components to allow full SED fitting. We therefore chose to limit the fitting to rest frame wavelengths corresponding to the MIPS 24 μm band (∼3–4 μm rest frame) and shorter where the majority of the sources is well detected. For these data we fit a combination of a power-law in the UV/optical and a hot blackbody in the NIR. To minimize the influence from emission lines (e.g., Lyα, Hα) and the small blue bump on the fitted power-law slope, we limit the data points to Spitzer bands at λ_obs ⩾ 5.8 μm and only using the y-band photometry in the rest frame UV. In those cases where no y-band photometry is available (five objects), we use the z-band instead. For selected sources the UV part of the rest frame SEDs was excluded from the fitting to avoid broad absorption line features (e.g., J0203+0012, J1427+3312). The resulting fits are shown in Figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We derive UV/optical luminosities for the quasars by integrating the fitted power-law between 0.1 μm and 1 μm. Similarly, for the hot dust component, the NIR dust luminosity is provided by integrating the fitted blackbody between 1 μm and 3 μm. The fitted values for α_UV/opt (F_ν∝ν^α) and T_NIR, as well as the integrated luminosities for the two components, are given in Table 7 and their distributions are shown in Figure 8.⁸

Zoom In Zoom Out Reset image size

In the distributions in Figure 8 we also indicate the Herschel FIR-detected objects (blue) and the partly detected objects (red; see Section 4.4 for the definition of these samples). While no specific trends can be identified for α_UV/opt or T_NIR, Figure 8 reveals that Herschel detections are preferentially found at the high end of the UV/optical luminosity distribution (see also, Netzer et al. 2013). This is even more pronounced for the NIR luminosity L_{NIR, dust} (Figure 8, bottom right).

We also find a group of objects which have very low temperatures of the hot dust component in our fitting approach (Figure 8, bottom left). We caution that the actual temperature values provided by our fitting are not well defined in these cases because the data only poorly constrain T_NIR. Nevertheless, the SEDs clearly demonstrate that these objects have a dearth of very hot dust compared to their UV/optical luminosity, and compared to the remainder of the sample. The reduced contributions from hot dust to the SEDs is also reflected in their lower values for L_{NIR, dust}.⁹ In the individual SED plots (Figure 7) these objects can be identified from their shallow rise in flux between the observed bands at 8 μm and 24 μm (Figure 9). Often, the observed IRAC 8 μm data point is still dominated by the power law and not by the onset of the hot dust emission as traced by the MIPS 24 μm photometry. In such cases the 24 μm photometry itself only very moderately exceeds the predictions from the power law.

Zoom In Zoom Out Reset image size

Altogether we find that ∼12%–16% of the sample have NIR to UV/optical properties that are quite different from the rest of the sample.¹⁰ Such sources have been found in similar proportions in other samples (e.g., Jiang et al. 2010; Hao et al. 2011; Mor & Netzer 2012). Jun & Im (2013) show that the fraction of dust-poor quasars increases with optical luminosity and redshift, and our numbers are consistent with their trends. These authors suggest that the dust-poor phase is a transient phenomenon during the evolution of the quasar (e.g., Jiang et al. 2010), rather than a distinct population of quasars with low covering factors (e.g., Hao et al. 2011).

Zoom In Zoom Out Reset image size

The rest of our objects has L_{NIR, dust}/L_UV/opt between ∼0.08 and 0.5 and we see no trends in this ratio with redshift or L_UV/opt (Figure 10). The latter implies that the luminosities in the UV/optical and NIR are well correlated for most objects (see also Mor & Netzer 2012). This is not surprising considering that the accretion disk emission, here traced by the UV/optical luminosity, is expected to be the primary heating source of the hot dust. Neither L_UV/opt nor L_NIR show any trend with α_UV/opt (Figure 10). Similarly, T_NIR shows no obvious trends with α_UV/opt or L_UV/opt, while α_UV/opt is uncorrelated with redshift.

By including the PACS 100 μm band, we extend the analysis to slightly longer infrared wavelengths and determine the flux at a rest frame wavelength of 6.7 μm through interpolation of the observed photometry at 24 μm and 100 μm using a power law in νF_ν. In combination with the monochromatic luminosity at 5100 Å (rest frame) as provided by our UV/optical power-law fits, we can study the monochromatic ratio of MIR-to-optical emission as a function of (monochromatic) optical luminosity (Figure 11). Under the assumption that the infrared-to-optical luminosity ratio is a proxy of the dust covering factor in (type 1) AGN, a plot as in Figure 11 has been used by Maiolino et al. (2007) to identify a trend where the dust covering factor decreases with increasing optical luminosity. Such a general behavior of the dust covering factor (or obscured fraction of AGN) has been detected for many different samples and using various techniques (e.g., Treister et al. 2008; Hasinger 2008; Lusso et al. 2013, and references therein) and is also seen in our sample for the Herschel-detected objects (Figure 6). However, the question whether the covering factor also changes with redshift remains controversial, with claims for (Treister & Urry 2006; Hasinger 2008) and against (Ueda et al. 2003; Lusso et al. 2013) significant redshift evolution.

Zoom In Zoom Out Reset image size

In this context, our high-redshift QSOs show a systematic, albeit very moderate offset in the MIR-to-optical luminosity ratio with respect to 2.0 ≲ z ≲ 3.5 QSOs (Figure 11). Much of the observed offset is currently driven by the Herschel detections where the 6.7 μm flux is determined as the interpolation between two significantly detected data points at λ_obs of 24 μm and 100 μm. However, about 60% of the z > 5 objects have only upper limits on the MIR-to-optical ratio, mostly due to non-detections in the 100 μm band. While in Figure 11 these objects currently populate the same area as the Herschel 100 μm detected objects (colored symbols), their effect on the observed trends remain unclear. This is in particular emphasized when considering the wide range of intrinsic SED shapes that may be present among the Herschel non-detected sources (see Section 4.4), which could potentially result in a wider range of luminosity ratios than seen currently for the Herschel-detected objects. For example, if the z > 5 objects intrinsically showed the same spread in the MIR-to-optical ratio as the 2.0 ≲ z ≲ 3.5 sample (almost 1 dex in Figure 11) then the resulting distribution would be roughly consistent with the observed trends at lower redshift.

On the other hand Figure 11 reveals that the data points and upper limits with the lowest MIR-to-optical ratio among our sample in Figure 11 almost exclusively belong to the group of objects where the SEDs indicate a dearth of hot dust. These objects may be of different nature or reside in a different evolutionary state (e.g., Jiang et al. 2010; Hao et al. 2011; Mor & Trakhtenbrot 2011; Jun & Im 2013) and possibly cannot be directly compared with the other QSOs from either sample.

Hα is one of the most prominent emission lines in the UV/optical spectra of common quasars (e.g., Vanden Berk et al. 2001). For our high-redshift objects (z > 5), this emission line is redshifted into the observed mid infrared, largely precluding direct spectroscopic observations with current facilities. However, we can use our high signal-to-noise Spitzer photometry to estimate Hα fluxes. For redshifts greater than z ∼ 5.2, the influence of this line can be seen in the individual SEDs (Figure 7) where Hα emission boosts the flux in the 4.5 μm IRAC band compared to a power-law continuum (e.g., J0840+5624). At lower redshift (z < 5.2), Hα falls onto the flanks of the filter transmission or largely into the small gap between the 3.6 μm and 4.5 μm IRAC filter. This makes it difficult to extract reliable emission-line flux estimates from the observed photometry. At z ≳ 5.2 and up to our maximum redshift (z = 6.42), the Hα emission line is fully covered by the filter transmission window and peaks within the plateau region of the 4.5 μm filter.¹¹

For the purpose of estimating Hα line fluxes we fit the SEDs slightly differently as compared to Section 4.2 or as shown in Figure 7. Instead of considering the full UV/optical continuum for a power-law fit we now limit the fit to the neighboring photometric points in order to isolate the local continuum. This means that for an Hα line falling into the 4.5 μm band we fit the power law to the 3.6 μm and 5.8 μm bands only. From the offset of the measured flux in the 4.5 μm filter compared to the local estimate of the power-law continuum we then calculate the Hα emission-line flux and equivalent width (EW).¹² We show the distribution of the estimated Hα EWs in Figure 12 and the derived values are also provided in Table 7. When comparing our high-z results to spectroscopic Hα EWs from low-redshift (z ≲ 0.4) SDSS quasars (Shen et al. 2011), we see that the two distributions are quite similar in width and shape. This similarity between low- and high-redshift quasars indicates a lack of redshift evolution in the Hα EWs, which agrees with similar results for rest-frame UV emission lines (Iwamuro et al. 2004; Juarez et al. 2009; De Rosa et al. 2011).

Zoom In Zoom Out Reset image size

From Figure 12 (left) we can also see that the Herschel-detected objects, and in particular the FIR-detected objects, have preferentially low Hα EWs. This trend is also seen in the Lyα EWs (Figure 12, top right) as taken from the literature (Table 7). Such a prevalence of FIR bright objects among sources with low Lyα EW has previously been indicated for millimeter-detected high-redshift quasars (e.g., Omont et al. 1996; Bertoldi et al. 2003; Wang et al. 2008a). Wang et al. (2008a) speculated that in these objects a special dust geometry that only affects the broad emission line clouds (and not the continuum) could in principle lead to such an effect. Because the impact of dust obscuration on Hα would be much smaller than for Lyα, the persisting trend of Herschel FIR-detected objects to be found at low Hα EW values questions such a scenario. However, while the effects of obscuration are indeed reduced for Hα compared to Lyα, they can still be non-negligible. A more definite answer requires higher precision direct spectroscopic measurements, preferably of (rest frame) NIR emission lines, to further reduce the effect of possible dust obscuration.

Due to the large number of Herschel non-detections in our sample, we have used a stacking approach to study the average infrared properties of the high-redshift quasars. For this purpose we divided the full sample into three subsamples.¹³

• 1.  
  Ten FIR-detected objects with detections in at least three Herschel bands (160 μm, 250 μm, and 350 μm).¹⁴
• 2.  
  Fourteen partly (Herschel) detected objects with significant PACS 100 μm and/or 160 μm flux. We refer to this subsample also as the PACS-only objects.¹⁵
• 3.  
  Thirty-three (Herschel) non-detections.

The remaining objects have been excluded from the stacking analysis on various grounds: 10 objects suffer from confusion with nearby FIR bright sources which would influence the stacked fluxes. The science targets on these images are not detected with Herschel individually. Two additional sources, J1148+5253 and J2245+0024, have been excluded because they are significantly fainter in the optical/UV than the rest of the sample. Both are also Herschel non-detections.

The stacking was performed in flux in the observed frame and the resulting mean SEDs were shifted into rest frame using the median redshift in the respective subsamples. In the y-band and in the Spitzer bands (where virtually all objects in all three subsamples are individually detected) we used the observed photometry as input. For Herschel the final images where stacked pixel by pixel centered on the position of the quasar. Photometry on the resulting stacks was performed as described for the individual images (Sections 2.6.1 and 2.6.2). The mean stacks in the Herschel bands for the three subsamples are presented in Figure 13. To estimate the variation present within these subsamples we followed a bootstrapping approach. For a given subsample we randomly selected as many objects as there are members in that subsample, allowing for replacements, created a new stack and performed photometry. This was done for 1000 random combinations of objects in each subsample. The centroid of the distribution of these 1000 individual stacked photometry values was then taken as the final average flux of the subsample. We use the standard deviation of this distribution, which can be considered a measure for the variety of intrinsic SED shapes present in the subsample, as the uncertainty on the average flux. The overall significance of the final stacked mean value in the Herschel bands was determined as follows: we stacked the images at random positions on the background, following a similar procedure as for the quasar positions. If the mean value of the source stack distribution is larger than three times the mean value of the background stack distribution we consider the stacked quasar signal to be significant.

Zoom In Zoom Out Reset image size

In Figure 14 (left) we compare the average SEDs of the FIR-detected objects (blue SED) with that of the partly detected objects (red SED). The SEDs are very similar in absolute scaling and spectral shape up to and including the observed 100 μm band (∼15 μm rest frame). At longer wavelengths, however, the SEDs are very different. In the νF_ν representation of Figure 14 (left), the PACS-only objects show a steep drop above ∼20 μm rest frame while the FIR-detected objects display an additional component toward the FIR. This behavior is emphasized in Figure 14 (right) where we show the average SEDs normalized by the shape of the mean SED of the partly detected objects.

Zoom In Zoom Out Reset image size

The partly Herschel-detected sources (red SED in Figure 14) are optically luminous AGN with powerful NIR and MIR emission, but without exceptional FIR brightness, at least on average. The shape of the SED is very similar to the average SDSS quasar SED and beyond ∼20 μm broadly resembles the shape of typical torus models (e.g., Schartmann et al. 2008; Nenkova et al. 2008; Hönig & Kishimoto 2010; Stalevski et al. 2012). In these cases the AGN is likely contributing significantly or even dominantly to the FIR emission (Netzer et al. 2007; Lutz et al. 2008; Wang et al. 2008a). However, the upper limits in the SPIRE bands are not very stringent and would still be consistent with an FIR component of ∼10¹² L_☉ (assuming a modified black body of T = 47 K and β = 1.6). Therefore, it cannot be ruled out that star formation contributes FIR emission on levels of a few tens to a few hundred solar masses per year as found for other high-redshift QSOs (Wang et al. 2008a; Venemans et al. 2012; Willott et al. 2013; Netzer et al. 2013). We note that for some combinations of objects the bootstrapping indeed reveals significant detections in the SPIRE bands, indicating that some sources in this subsample were just below the individual detection limit. In the global mean, however, the partly Herschel-detected subsample only reaches ∼2σ significance in the stacked values at 250 μm and 350 μm.

The comparison of the average SEDs of the two Herschel-detected subsamples emphasizes that AGN with strong emission in the UV/optical (from the accretion disk) throughout the NIR and MIR (from AGN-powered hot and warm dust) do not necessarily show considerable FIR emission as well (both panels in Figure 14; see also Dai et al. 2012). The fact that some optical and MIR luminous QSOs show strong FIR emission and others do not may indicate that star formation is the dominant driver for the additional FIR component observed in the Herschel FIR-detected QSOs. This is consistent with the results of Lutz et al. (2008) who show that for a sample of millimeter bright QSOs at z ∼ 2 the PAH and FIR luminosities correlate, which also supports star formation as the source of the FIR emission in powerful FIR bright AGN (these objects have optical luminosities comparable to our sources). For some high-redshift millimeter bright QSOs it has been shown that the potentially star-formation dominated FIR continuum and line emission (e.g., of [C ii]) is concentrated in the innermost kiloparsecs (e.g., Walter et al. 2009; Wang et al. 2013). At the highest luminosities, the close proximity to the AGN may thus still lead to significant contributions from the AGN to the FIR emission (e.g., Dai et al. 2012, see also Valiante et al. 2011). For such sources, future observation at high resolution may provide additional clues on the relative AGN-to-SF contributions to the dust heating from the spatial distribution of the cold dust emission.

The average SED of the Herschel non-detections (green SEDs in Figure 14) differs from the SEDs of the Herschel-detected sources in several aspects. Even at 100 μm (observed) we only find a barely significant detection in the maximum flux case (as provided by the bootstrapping), and non-detections for most realizations as well as for the mean of this subsample. No detection was obtained for any combination of objects at longer wavelengths. A significant difference between the average SED of this subsample and those of objects with individual Herschel detections is that the UV/optical and NIR/MIR fluxes are systematically smaller, and accordingly also the luminosities because the median redshifts of the three subsamples are very similar. This has already been indicated by the UV/optical and NIR luminosity distributions of the full sample (Figure 8) where the Herschel-detected sources are found to be clustered at the high luminosity end. This flux difference in the mean increases toward longer wavelengths (factor of 1.4 at ∼0.5 μm rest frame and factor of five at ∼15 μm rest frame; see also Blain et al. 2013). From Figure 14 it appears that on average the shape of the infrared SED is changing for fainter QSOs, even though taking into account the full errors bars shown in that figure could reduce the trend seen for the mean.¹⁶ This behavior is in principle supported by the individual objects that have deeper re-observations. These were drawn from the Herschel non-detected subsample, and by selection correspond to sources with high observed 24 μm fluxes in this sample (above the average for all but one of the six sources). Individually they also show a decline in νF_ν between 24 μm and 100 μm, but shallower than the average SED.

The flux levels of the Herschel/PACS 100 μm detections for these deeper data show that they were not far below the sensitivity limit of our standard observations. Even though a number of objects apparently only barely avoided detection in our standard data, the final average SED of the Herschel non-detections (which includes these objects just below the detection limit) shows steeper 24 μm to 100 μm slopes and remains Herschel PACS non-detected in most cases. This supports the idea that the Herschel non-detected subsample includes objects with a wide range of intrinsic SED shapes and that many of the sources in this sub sample are much fainter in the 100 μm band than our detection limit.

We note that four objects in this Herschel non-detected subsample have individual millimeter detections (Bertoldi et al. 2003; Wang et al. 2007, 2011) and were discussed in Leipski et al. (2013). Considering the shape of the SED in the (AGN dominated) MIR, we can speculate that the cold dust emission responsible for the millimeter flux (typically ∼2–3 mJy at 1.2 mm observed) is probably powered by star formation in these cases (see also Wang et al. 2008a; Leipski et al. 2013).