© ESO, 2016

1. Introduction

Quasars (or quasi stellar objects, QSOs) occupy a special place in large-scale structure and galaxy evolution (Kormendy & Richstone 1995). They are among the most luminous extra-galactic sources, and, as such, have become the focus of many cosmological surveys such as the 2dF Quasar Redshift Survey (2QZ; Croom et al. 2001) or the successive iterations of the Sloan Digital Sky Survey (SDSS I-IV; York et al. 2000; Eisenstein et al. 2011), where they provide a unique insight into the formation of structure on large scales and at high redshift (White et al. 2012). Quasars are likely powered (Salpeter 1964; Lynden-Bell 1969) by accretion of nearby matter onto supermassive black holes (SMBH, Rees 1984), whose evolution appears closely linked to the general cosmic star formation rate (Madau & Dickinson 2014), in particular in galaxies that contain a massive bulge (and therefore a massive central black hole) and a gas reservoir (Nandra et al. 2007; Silverman et al. 2008). The link is a clue to galaxy formation, its significance emphasized by the observed relation between the SMBH mass and galaxy properties (Ferrarese & Merritt 2000; Gebhardt et al. 2000). In addition, galaxy formation models must evoke strong feedback from active galactic nucleus (AGN) to explain the observed properties of massive galaxies and to avoid the overcooling catastrophe (Croton et al. 2006; Bower et al. 2006; McNamara & Nulsen 2007; Somerville et al. 2008; Fabian 2012; Blanchard et al. 1992)

Quasar environments give us a look at these powerful engines, and millimeter/submillimeter observations can offer a particularly revealing view of their effects: at these frequencies it is possible to study the cooler dust emission associated with star formation in the host, look for synchrotron emission from energetic particles, and, perhaps most pertinently, to measure energy feedback through observation of the thermal Sunyaev-Zeldovich (tSZ, Sunyaev & Zeldovich 1970, 1972) effect, a direct probe of the thermal energy contained in surrounding gas. Samples observed in this waveband with large ground-based facilities typically consist of, at most, several tens of objects (Omont et al. 2001, 2003, 2013; Isaak et al. 2002; Priddey et al. 2003; Beelen et al. 2006; Wang et al. 2007). Operating at these same frequencies, cosmic microwave background (CMB) surveys cover large sky areas and present the opportunity of studying much larger samples, albeit with much less detail of individual objects. We can, however, determine the average properties of large representative populations, a valuable compliment to the more detailed examinations.

The Wilkinson Microwave Anisotropy Probe (WMAP, Bennett et al. 2003) and Planck missions (Tauber et al. 2010; Planck Collaboration I 2011) are ideal for this purpose because of their all-sky coverage, including the entire SDSS (York et al. 2000) area. Cross-correlating WMAP data with photometric QSOs from SDSS data release 3 (DR3), Chatterjee et al. (2010) finds evidence at 2.5σ for the tSZ effect from QSO environments. Ruan et al. (2015) improved the significance using a publicly available tSZ map (Hill & Spergel 2014) constructed from the Planck mission dataset and 26 686 spectroscopic QSOs from SDSS DR7. They concluded that the total thermal energy feedback into surrounding gas was significantly larger than expected according to galaxy formation models, which is consistent with the findings of Chatterjee et al. (2010). Planck’s wide frequency coverage, from 30 GHz to 857 GHz, is a distinct advantage for disentangling the different sources of emission from the QSO environment, and was implicitly exploited by Ruan et al. (2015) when using the tSZ map published by Hill & Spergel (2014).

Although not all-sky, the Atacama Cosmology Telescope (ACT) surveyed several hundred square degrees near the celestial equator, overlapping AGN and QSO samples in the North. Gralla et al. (2014) stacked millimeter maps from ACT over this area of radio-loud AGN seen in the FIRST and NVSS radio surveys to find a tSZ signal at 5σ significance. More recently, and in parallel to our present study, Crichton et al. (2015) combined ACT and Herschel data of radio-quiet quasars in the SDSS DR7 and DR10 samples, detecting the tSZ signal at 3−4σ significance.

In this paper, we use the full-mission Planck frequency maps to examine the emission properties of the SDSS-III Baryon Acoustic Oscillation (BOSS, Dawson et al. 2013) DR12 QSO sample. We extract not only the tSZ signal from the population, but also dust and synchrotron emission by stacking measurements made with a set of matched filters directly applied to the seven highest frequency Planck channels. The joint extraction enables us to study not only the gas thermal environment, but star formation conditions and the production of energetic particles. Moreover, it improves analysis robustness by giving information on correlations between the observed signals, compared to the use of a tSZ map. Ruan et al. (2015) were in fact faced with the difficulty of correcting the SZ map for contamination from dust emission.

We extend our multi-matched filter (MMF, Melin et al. 2006; see also Herranz et al. 2002) formalism to simultaneously extract two or three emission components, following development started in Planck Collaboration Int. XI (2013). Stacking the various signals to study their evolution with redshift, we reconstruct for the first time a picture of the evolution of the dust signal, the synchrotron emission, and the tSZ effect between z ~ 0 and z ~ 5 for the general QSO population.

In Sect. 2, we describe the BOSS and Planck datasets. We detail the analysis methods and tools in Sect. 3. In Sect. 4, we apply them on simulations. Results on Planck data are given in Sect. 5. We discuss our results in Sect. 6 and conclude in Sect. 7. Throughout, we use the spatially flat base ΛCDM cosmology from Planck Collaboration XIII (2015): H₀ = 67.27 = h100 km s^-1/Mpc, Ω_m = 0.3156, Ω_bh² = 0.02225 and σ₈ = 0.831.

2. Data

2.1. BOSS quasars

One of the major goal of the SDSS-III Baryon Oscillation Spectroscopic Survey (BOSS) is the production of a QSO catalogue to detect the BAO scale in the Lyman-α forests at redshift z ~ 2.5. A first detection was made on the DR9 QSO catalogue (Pâris et al. 2012) by Busca et al. (2013), confirmed latter by Delubac et al. (2015). In this study, we use the recently published DR12 QSO catalogue¹ (Alam et al. 2015).

The sources are detected and selected with the CCD imaging Camera installed in the Sloan Foundation 2.5 m Telescope at Apache Point observatory, New Mexico. The spectra of each source is measured with the BOSS spectrograph which covers wavelength between 3600 Å and 10 000 Å. Spectra measurement leads to the computation of the spectroscopic redshift and the checking of the nature of the source.

The DR12 QSO catalogue contains 297 301 objects. We remove 5880 QSO falling outside the Planck 65% mask, in a region strongly contaminated by the Milky Way dust. Finally, we also reject 256 QSO because they are at low redshift z< 0.1 so may be partially resolved by Planck , because they belong to the high redshift (z> 5) population or because they have a bad estimate of the magnitude in g band (g< 0) . We thus use a catalogue of 291 165 QSO in the redshift range 0.1 <z< 5. The distribution of the QSO in term of redshift is not flat and reach two notable maxima at z = 0.8 and z = 2.3 (Fig. 1).

The target selection was developed to select quasars with an observable Ly-α forest (i.e. quasars with z> 2.2). However, degeneracies in the color-redshift relation of quasars led to the selection of low-z quasars in BOSS (Fig. 1). The quasars at z ~ 0.8 have MgII λ2800 line at the same wavelength as Ly-α at redshift z ~ 3.1, giving these objects similar broadband colors, while the large number of objects at z ~ 1.6 is due to the confusion between λ1549 C-IV line and Ly-α at z ≈ 2.3 (Ross et al. 2012). In contrast, the selection based on the intrinsic variability of quasars gives a more uniform distribution in redshift (Palanque-Delabrouille et al. 2011).

Fig. 1 Redshift distribution of the BOSS DR12 QSO sample. The distribution presents two peaks, at z ~ 0.8 and z ~ 2.3, due to the QSO target selection process (see text). The regular binning used throughout this paper (Δz = 0.5) is shown as the vertical dotted lines.

2.2. Planck maps

The Planck satellite was launched in May 2009 (Tauber et al. 2010). After traveling to the Earth-Sun Lagrange point L2, it scanned the entire sky continuously from August 2009 to October 2013 in nine frequency bands ranging from 30 to 857 GHz. The primary goal of the mission was to study the primary CMB anisotropies, but its large frequency coverage also enables unique Galactic and extra-galactic astrophysical studies. In particular, the High Frequency Instrument (HFI, Lamarre et al. 2010; Planck HFI Core Team IV. 2011), covering bands from 100 to 857 GHz, is ideally suited for Sunyaev-Zeldovich science. Moreover, the three highest frequencies (353 to 857 GHz) of the instrument pick up Galactic and extra-galactic dust emission. The Low Frequency Instrument (LFI, Bersanelli et al. 2010; Mennella et al. 2011), with channels at 30, 44, 70 GHz, is sensitive to synchrotron, free-free, and spinning dust emission.

We use the seven highest frequency (70, 100, 143, 217, 353, 545, and 857 GHz) full-sky maps from the 2015 data release (Planck Collaboration I 2015; Planck Collaboration VIII 2015; Planck Collaboration VI 2015). These bands are the best suited to study the dust, gas, and synchrotron emission from high redshift QSO environments. The substantially larger beams of the remaining 30 and 40 GHz channels make them less well adapted for studying what are essentially point sources. We divide each Planck HEALPix² all-sky map into 504 overlapping flat tangential patches of 10 × 10 degrees with a 1.72 arcmin pixel scale in order to apply our extraction algorithms described in Sect. 3.2.

3. Analysis

Determining the physical properties of the QSO environment from the Planck maps faces two major challenges. The first is the faint QSO flux in the Planck maps, far below the noise level; individual detection is not possible. We thus have to use a statistical approach similar to the one developed in Melin et al. (2011) and subsequently in Planck Collaboration X (2011), Planck Collaboration XII (2011), Planck Collaboration Int. XI (2013). Details of the approach are given in Sect. 3.4.1.

The second challenge arises from the superposition of different emission sources from the QSO in the Planck beams: emission from dust located inside the host galaxy or possibly at larger scale in the host halo of the QSO; synchrotron emission from the host galaxy or from relativistic outflows of the central AGN; and the tSZ effect due to the hot gas surrounding the QSO and within the host halo. Planck does not have the spatial resolution to separate these components, but its good spectral coverage enables us to disentangle these sources of emission.

We separate the signals using multi-component matched multi-filters (MMF), an extension of the approach used in Planck Collaboration Int. XI (2013). The detailed description is given in Sect. 3.2.2. We will see that the QSO signal in Planck is dominated by dust emission, but that we also detect both synchrotron and tSZ signals. In practice, we proceed as follows:

• 1.

  assume that the QSO signal is a mixture of one, two or three components;

• 2.

  apply the adapted multi-component MMF at each QSO position to obtain the amplitude of each component;

• 3.

  bin average (as a function of redshift or magnitude) these amplitudes over the QSO sample;

• 4.

  evaluate the ability of our model to describe the data using a χ²-test (see Sect. 3.4.2).

Our notation in the subsequent sections is as follows. We use the letter p as index to denote individual QSOs in the catalogue and the indices i and j to denote observation frequencies. The Greek indices λ and μ specify the nature of the emission components, e.g., λ and μ have as possible values dust, synch, and tSZ for dust, synchrotron, and tSZ signals, respectively.

3.1. Model of the QSO emission

We model the QSO emission as a sum of dust, synchrotron, and tSZ signals. The seven-frequency column vector map, m(x), at position x on the sky is written as ${m}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}\sum _{\mathit{\lambda }}{\mathit{a}}_{\mathit{\lambda }}\hspace{0.17em}{{t}}_{\mathit{\lambda }}\mathrm{\left(}\mathit{z,}{x}\mathrm{\right)}\mathrm{+}{n}\mathrm{\left(}{x}\mathrm{\right)}\mathit{,}$(1)where a_λ is the amplitude of the λ component, t_λ(z,x) is the associated, normalized emission vector, z the redshift of the QSO, and n(x) is the astrophysical and instrumental noise. We center the QSO in the map to simplify the expressions. The ith component of each emission vector is the signal profile, τ_λ (normalized to unity at the origin), convolved by the Planck beam, b_i (normalized to unity at the origin), at frequency ν_i, and then scaled by the expected frequency dependance ${{}^{\mathrm{(}}\mathit{S}_{\mathit{\lambda }}{}^{\mathrm{)}}}_{\mathit{i}}\mathrm{\left(}\mathit{z}\mathrm{\right)}$: $\begin{array}{ccc}{\mathrm{(}{\mathit{t}}_{\mathit{\lambda }}\mathrm{)}}_{\mathit{i}}\mathrm{\left(}\mathit{z,}{x}\mathrm{\right)}\mathrm{=}{\mathrm{(}{\mathit{S}}_{\mathit{\lambda }}\mathrm{)}}_{\mathit{i}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\hspace{0.17em}\mathrm{[}{\mathit{b}}_{\mathit{i}}\mathrm{\ast }{\mathit{\tau }}_{\mathit{\lambda }}\mathrm{]}\mathrm{\left(}{x}\mathrm{\right)}\mathit{.}& & \end{array}$(2)We describe the dust emission frequency dependence by a modified blackbody, the synchrotron emission by a power law, and we use the non-relativistic calculation for the tSZ (e.g., Birkinshaw 1999; Carlstrom et al. 2002): $\begin{array}{ccc}& & {\mathrm{(}{\mathit{S}}_{\mathrm{d}}\mathrm{)}}_{\mathit{i}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\left(\frac{{\mathit{\nu }}_{\mathit{i}}}{\mathrm{857}\mathrm{GHz}}\right)}^{{\mathit{\beta }}_{\mathrm{d}}}\frac{{\mathit{B}}_{{\mathit{\nu }}_{\mathit{i}}}\mathrm{[}{\mathit{T}}_{\mathrm{d}}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{\left)}\mathrm{\right]}}{{\mathit{B}}_{\mathrm{857}\mathrm{GHz}}\hspace{0.17em}\mathrm{[}{\mathit{T}}_{\mathrm{d}}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{\left)}\mathrm{\right]}}\mathit{,}\\ & & {\mathrm{(}{\mathit{S}}_{\mathrm{s}}\mathrm{)}}_{\mathit{i}}\mathrm{=}{\left(\frac{\mathrm{100}\mathrm{GHz}}{{\mathit{\nu }}_{\mathit{i}}}\right)}^{{\mathit{\alpha }}_{\mathrm{s}}}\mathit{,}\\ & & {\mathrm{(}{\mathit{S}}_{\mathrm{sz}}\mathrm{)}}_{\mathit{i}}\mathrm{=}{\mathit{B}}_{{\mathit{\nu }}_{\mathit{i}}}\mathrm{\left(}{\mathit{T}}_{\mathrm{cmb}}\mathrm{\right)}\frac{{\mathit{x}}_{\mathit{i}}{\mathit{e}}^{{\mathit{x}}_{\mathit{i}}}}{{\mathit{e}}^{{\mathit{x}}_{\mathit{i}}}\mathrm{-}\mathrm{1}}\left[\frac{{\mathit{x}}_{\mathit{i}}}{\tanh \mathrm{(}\mathit{x}\mathit{/}\mathrm{2}\mathrm{)}}\mathrm{-}\mathrm{4}\right]\mathit{,}\end{array}$where x_i = hν_i/kT_cmb with T_cmb the temperature of the CMB today, B_ν(T) represents the Planck spectrum, T_d and β_d are the dust temperature and emissivity index, and α_s = 0.7 for the synchrotron index. It is the typical value for radio galaxies (Condon 1992; Peterson 1997). We note that the first two expressions are unitless, so that a_dust and a_synch carry units of brightness; on the other hand, S_tSZ carries brightness units while a_tSZ is unitless and corresponds to the central Compton-y parameter.

We adopt the same spatial profile for all signal templates, τ_λ = τ(x/θ_f), where τ is taken from Arnaud et al. (2010) based on the observed cluster halos gas pressure profile. We fix θ_f = 0.27 arcmin, the value corresponding to a halo of mass M₅₀₀ = 10¹³M_⊙ at z = 2. This value is smaller than the smallest Planck beam (fwhm = 4.2 arcmin at 857 GHz): the bulk of the QSO remain unresolved by Planck. We study the sensitivity of our results to this fixed size in Sect. 5.5. The profiles are cut at θ> 5θ₅₀₀.

In the following we therefore denote this universal QSO profile by τ_qso. We are then able to define the QSO brightness, f, as the sum over the different components as ${\mathit{f}}_{\mathrm{i}}\mathrm{=}\sum _{\mathit{\lambda }}{\mathit{a}}_{\mathit{\lambda }}\hspace{0.17em}{\mathrm{\left(}{\mathit{S}}_{\mathit{\lambda }}\mathrm{\right)}}_{\mathit{i}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{,}$(6)and rewrite the map at frequency i as ${\mathit{m}}_{\mathit{i}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathit{i}}\hspace{0.17em}{\mathrm{\left(}{\mathit{t}}_{\mathrm{qso}}\mathrm{\right)}}_{\mathit{i}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{+}{\mathit{n}}_{\mathit{i}}\mathrm{\left(}{x}\mathrm{\right)}\mathit{,}$(7)where ${{}^{\mathrm{(}}\mathit{t}_{\mathrm{qso}}{}^{\mathrm{)}}}_{\mathit{i}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{\equiv }\mathrm{[}{\mathit{b}}_{\mathit{i}}\mathrm{\ast }{\mathit{\tau }}_{\mathrm{qso}}\mathrm{]}\mathrm{\left(}{x}\mathrm{\right)}$.

3.2. Matched filters

3.2.1. Single frequency matched filter

We use single frequency matched filters to extract individual QSO signals (at S/N ≪ 1) at each Planck frequency. The individual signals are necessary to compute the χ² statistic described in Sect. 3.4.2. We also average them to obtain the mean QSO spectrum in Sect. 5.

The single frequency matched filter is a linear filter designed to return an unbiased signal estimate with minimal error (Haehnelt & Tegmark 1996; Sanz et al. 2001). No assumption is made concerning the frequency dependance of the signal. The estimated value, , of the true signal, f_i, is $\mathit{f̂}\mathrm{i}\mathrm{=}\mathrm{\int }{\mathrm{d}}^{\mathrm{2}}\mathit{x}{\mathrm{\Psi }}_{\mathrm{i}}^{\mathit{t}}\mathrm{\left(}{x}\mathrm{\right)}\hspace{0.17em}{\mathit{m}}_{\mathrm{i}}\mathrm{\left(}{x}\mathrm{\right)}\mathit{,}$(8)where Ψ_i(x) is the matched filter for frequency ν_i. The filter Ψ_i is expressed in Fourier space as ${\mathrm{\Psi }}_{\mathrm{i}}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{=}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left[}\mathit{f̂}\mathrm{i}\mathrm{\right]}\frac{{\mathrm{(}{{\mathit{t}}_{\mathrm{qso}}^{\mathrm{\ast }}}^{\mathrm{)}}}_{\mathit{i}}\mathrm{\left(}{k}\mathrm{\right)}}{{\mathit{P}}_{\mathit{ii}}\mathrm{\left(}{k}\mathrm{\right)}}\mathit{,}$(9)with the error on the estimated signal $\begin{array}{ccc}\mathit{\sigma }\mathrm{\left[}\mathit{f̂}\mathrm{i}\mathrm{\right]}& & \mathrm{\equiv }\sqrt{\mathrm{⟨}\mathit{f̂}\begin{array}{c}\mathrm{2}\\ \mathrm{i}\end{array}\mathrm{⟩}\mathrm{-}\mathrm{⟨}{\mathit{f̂}\mathrm{i}\mathrm{⟩}}^{\mathrm{2}}}\\ & \mathrm{=}& {\left[\mathrm{\int }{\mathrm{d}}^{\mathrm{2}}\mathit{k}\frac{{\left|{\mathrm{(}{\mathit{t}}_{\mathrm{qso}}\mathrm{)}}_{\mathit{i}}\mathrm{\left(}{k}\mathrm{\right)}\right|}^{\mathrm{2}}}{{\mathit{P}}_{\mathit{ii}}\mathrm{\left(}{k}\mathrm{\right)}}\right]}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}\end{array}$(10)where P_ii(k) is the noise power spectrum of the ith Planck map, m_i(k). Because the signal-to-noise of an individual QSO signal is smaller than unity, we estimate P_ii(k) directly as the power spectrum of the map.

3.2.2. Matched multi-filters

The MMF were first introduced by Herranz et al. (2002) for SZ detection. They were further developed by Melin et al. (2006) and extensively used to construct the successive Planck SZ cluster catalogs (Planck Collaboration VIII 2011; Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2015). The MMF target the extraction of a single component across a set of maps, assuming a known frequency dependance and a spatial distribution. The filtered map returns a unbiased estimate, â_λ, of a_λ with minimal variance: $\mathit{â}\mathit{\lambda }\mathrm{=}\mathrm{\int }{\mathrm{d}}^{\mathrm{2}}\mathit{x}{{\Psi }}_{\mathit{\lambda }}^{\mathit{t}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{·}{m}\mathrm{\left(}{x}\mathrm{\right)}\mathit{.}$(11)Similarly to the single frequency case, the MMF for the λ emission, Ψ_λ(x), is expressed in Fourier space as ${{\Psi }}_{\mathit{\lambda }}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{=}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left[}\mathit{â}\mathit{\lambda }\mathrm{\right]}\hspace{0.17em}{{P}}^{-1}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{·}{{t}}_{{\lambda }}^{{\ast }}\mathrm{\left(}\mathit{z,}{k}\mathrm{\right)}\mathit{.}$(12)The error on the amplitude is now given as $\begin{array}{ccc}\mathit{\sigma }\mathrm{\left[}\mathit{â}\mathit{\lambda }\mathrm{\right]}& & \mathrm{\equiv }\sqrt{\mathrm{⟨}\mathit{â}\begin{array}{c}\mathrm{2}\\ \mathit{\lambda }\end{array}\mathrm{⟩}\mathrm{-}\mathrm{⟨}{\mathit{â}\mathit{\lambda }\mathrm{⟩}}^{\mathrm{2}}}\\ & \mathrm{=}& {\left[\mathrm{\int }{\mathrm{d}}^{\mathrm{2}}\mathit{k}{{{t}}_{{\lambda }}^{{\ast }}}^{\mathit{t}}\mathrm{\left(}\mathit{z,}{k}\mathrm{\right)}\mathrm{·}{{P}}^{-1}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{·}{{t}}_{{\lambda }}\mathrm{\left(}\mathit{z,}{k}\mathrm{\right)}\right]}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}\end{array}$(13)with P(k) being the inter-band cross-power spectrum matrix with contributions from (non-λ) sky signal and instrumental noise. It is the effective noise matrix for the MMF and, as for the single frequency power spectrum, can be estimated directly on the data, since the λ signal (dust, synchrotron or tSZ) is small compared to other astrophysical signals over the sky patch.

Introduced in Planck Collaboration Int. XI (2013) for a mixture of dust and tSZ signals, multi-component filtering deals with different MMF to separate dust and tSZ signals. Multi-component filtering has also been successfully used to jointly extract the tSZ and kinetic SZ (kSZ) effects in Planck Collaboration Int. XIII. (2014). For these filters, the recovered amplitude of each component is unbiased with minimal variance if the assumption on the number and type of the components is correct.

We consider three multi-component filtering schemes: 1) dust + tSZ; 2) dust + synchrotron; and 3) dust + tSZ + synchrotron filters. The amplitudes a_λ are estimated using a linear combination of the partial MMFs, $\mathit{â}\mathit{\lambda }\mathrm{=}\sum _{\mathit{\mu }}{\mathrm{\left[}{{D}}^{-1}\mathrm{\right]}}_{\mathit{\lambda ,\mu }}\hspace{0.17em}\mathrm{\int }{\mathrm{d}}^{\mathrm{2}}\mathit{x}{{\Phi }}_{\mathit{\mu }}^{\mathit{t}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{·}{m}\mathrm{\left(}{x}\mathrm{\right)}\mathit{,}$(14)with the partial MMF defined in Fourier space by ${{\Phi }}_{\mathit{\mu }}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{=}{{P}}^{-1}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{·}{{t}}_{{\mu }}\mathrm{\left(}\mathit{z,}{k}\mathrm{\right)}\mathit{,}$(15)and the matrix D computed as ${\mathit{D}}_{\mathit{\lambda ,\mu }}\mathrm{=}\left[\mathrm{\int }{\mathrm{d}}^{\mathrm{2}}\mathit{k}{{{t}}_{{\lambda }}^{{\ast }}}^{\mathit{t}}\mathrm{\left(}\mathit{z,}{k}\mathrm{\right)}\mathrm{·}{{P}}^{-1}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{·}{{t}}_{{\mu }}\mathrm{\left(}\mathit{z,}{k}\mathrm{\right)}\right]\mathit{.}$(16)The estimated amplitudes are combined together in a N-component vector (N being the number of signal components considered), $\hat{{a}}$. They are correlated with covariance matrix $\begin{array}{ccc}\mathit{C}\mathrm{[}\hat{{a}}{\mathrm{]}}_{\mathit{\lambda ,\mu }}& & \mathrm{\equiv }\mathrm{⟨}\mathit{â}\mathit{\lambda }\mathit{â}\mathit{\mu }\mathrm{⟩}\mathrm{-}\mathrm{⟨}\mathit{â}\mathit{\lambda }\mathrm{⟩}\mathrm{⟨}\mathit{â}\mathit{\mu }\mathrm{⟩}\\ & \mathrm{=}& {\mathrm{[}{{D}}^{-1}\mathrm{]}}_{\mathit{\lambda ,\mu }}\mathit{.}\end{array}$(17)In summary, we use the single frequency matched filter plus the following five filters in our study:

• single component MMF for dust;

• single component MMF for tSZ;

• two-component MMF for dust + tSZ;

• two-component MMF for dust + synchrotron;

• three-component MMF for dust + tSZ + synchrotron.

3.3. From observed signals to physical QSO properties

The amplitudes â_λ (determined with profiles cut at θ< 5θ₅₀₀) are converted into spherically integrated quantities within R₅₀₀ by $\begin{array}{ccc}\begin{array}{c}\mathrm{􏽢}\\ {\mathit{A}}_{\mathit{\lambda }}\end{array}\mathrm{=}\mathit{â}\mathit{\lambda }\hspace{0.17em}{\mathrm{\int }}_{\mathit{\theta }\hspace{0.17em}\mathit{<}\hspace{0.17em}\mathrm{5}{\mathit{\theta }}_{\mathrm{500}}}\mathrm{d\Omega }\hspace{0.17em}{\mathit{\tau }}_{\mathit{\lambda }}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{\times }{\mathit{C}}_{\mathit{\lambda }}\mathrm{(}\mathrm{5}{\mathit{R}}_{\mathrm{500}}\mathrm{\to }{\mathit{R}}_{\mathrm{500}}\mathrm{)}\mathit{,}& & \end{array}$(18)where the factor C_λ(5R₅₀₀ → R₅₀₀) is the conversion from the volume within 5R₅₀₀ to R₅₀₀, given the adopted profile. Similarly, we obtain the source flux density vector, $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {F}\end{array}$, as $\begin{array}{c}\mathrm{􏽢}\\ {F}\end{array}\mathrm{=}\hat{{f}}\hspace{0.17em}\hspace{0.17em}{\mathrm{\int }}_{\mathit{\theta }\hspace{0.17em}\mathit{<}\hspace{0.17em}\mathrm{5}{\mathit{\theta }}_{\mathrm{500}}}\mathrm{d\Omega }\hspace{0.17em}{\mathit{\tau }}_{\mathrm{qso}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{\times }{\mathit{C}}_{\mathrm{qso}}\mathrm{(}\mathrm{5}{\mathit{R}}_{\mathrm{500}}\mathrm{\to }{\mathit{R}}_{\mathrm{500}}\mathrm{)}\mathit{.}$(19)We express the quantities $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{A}}_{\mathrm{dust}}\end{array}$, $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {F}\end{array}$, and $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{A}}_{\mathrm{synch}}\end{array}$ in mJy and $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{A}}_{\mathrm{tSZ}}\end{array}$ in arcmin².

Following Beelen et al. (2006) we estimate the total dust mass of the QSO from $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{A}}_{\mathrm{dust}}\end{array}$ as $\begin{array}{c}\mathrm{􏽢}\\ {\mathit{M}}_{\mathrm{dust}}\end{array}\mathrm{=}\begin{array}{c}\mathrm{􏽢}\\ {\mathit{A}}_{\mathrm{dust}}\end{array}\frac{{\mathit{D}}_{\mathrm{L}}\mathrm{(}\mathit{z}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}{\mathit{\kappa }}^{-1}\left[\mathrm{857}\hspace{0.17em}\mathrm{GHz}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\right]{\mathit{B}}_{\mathrm{857}\mathrm{GHz}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}^{-1}\mathrm{\left(}{\mathit{T}}_{\mathit{d}}\mathrm{\right)}\mathit{,}$(20)with $\mathit{\kappa }\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}{\mathit{\kappa }}_{\mathrm{0}}{\left(\frac{\mathit{\nu }}{\mathrm{249.8}\mathrm{GHz}}\right)}^{{\mathit{\beta }}_{\mathrm{d}}}$, κ₀ = 0.4 g^-1 cm² and D_L(z) the luminosity distance. From the synchrotron flux density at 100 GHz, we compute the monochromatic synchrotron luminosity according to $\begin{array}{c}\mathrm{􏽢}\\ {\mathit{L}}_{\mathrm{synch}}\end{array}\mathrm{=}\begin{array}{c}\mathrm{􏽢}\\ {\mathit{A}}_{\mathrm{synch}}\end{array}\hspace{0.17em}\frac{{\mathit{D}}_{\mathrm{L}}\mathrm{(}\mathit{z}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{{\mathit{\alpha }}_{\mathrm{s}}}}\hspace{0.17em}$(21)which is referenced to the rest-frame frequency of 100 GHz assuming our power law with α_s.

We define the intrinsic Compton parameter, $\begin{array}{c}{\mathrm{􏽢}}^{}\\ {\mathit{Y}}_{\mathrm{500}}\end{array}$, as the quantity $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{Y}}_{\mathrm{500}}\end{array}$ from Planck Collaboration Int. XI (2013), $\begin{array}{c}\mathrm{􏽢}\\ {\mathit{Y}}_{\mathrm{500}}\end{array}\mathrm{=}\begin{array}{c}\mathrm{􏽢}\\ {\mathit{A}}_{\mathrm{tSZ}}\end{array}\mathit{E}\mathrm{(}\mathit{z}{\mathrm{)}}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{3}}{\left(\frac{{\mathit{D}}_{\mathrm{A}}\mathrm{\left(}\mathit{z}\mathrm{\right)}}{\mathrm{500}\hspace{0.17em}\mathrm{Mpc}}\right)}^{\mathrm{2}}\mathit{,}$(22)D_A(z) being the angular distance to the QSO in Mpc. Applying Eq. (22) of Arnaud et al. (2010), we estimate the total mass of the halo as $\begin{array}{c}\mathrm{􏽢}\\ {\mathit{M}}_{\mathrm{500}}\end{array}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{b}\mathrm{)}\begin{array}{c}\mathrm{􏽢}\\ {\mathit{M}}_{\mathrm{500}\mathit{,}\mathrm{true}}\end{array}\mathrm{=}\mathrm{3}\mathrm{\times }{\mathrm{10}}^{\mathrm{14}}{\left(\frac{\mathit{h}}{\mathrm{0.7}}\right)}^{-1}{\left[\frac{\begin{array}{c}\mathrm{􏽢}\\ {\mathrm{A}}_{\mathrm{tSZ}}\end{array}\mathit{E}\mathrm{(}\mathit{z}{\mathrm{)}}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{3}}}{{\mathit{A}}_{\mathrm{x}}\mathrm{\left(}\mathit{z}\mathrm{\right)}}\right]}^{\mathrm{1}\mathit{/}\mathit{\alpha }}{\mathit{M}}_{\mathrm{\odot }}\mathit{,}$(23)with α = 1.78 and ${\mathit{A}}_{\mathrm{x}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{2.925}\mathrm{\times }{\mathrm{10}}^{-5}\mathrm{\times }\mathrm{0.6145}\mathrm{\times }{{}^{\mathrm{\right(}}\frac{\mathit{h}}{\mathrm{0.7}}^{\mathrm{\left)}}}^{-1}\mathrm{\times }{\mathit{D}}_{\mathrm{A}}^{-2}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{\times }{{}^{\mathrm{(}}\frac{\mathit{\pi }}{\mathrm{180}}{\frac{\mathrm{1}}{\mathrm{60}}}^{\mathrm{)}}}^{-2}{\mathrm{arcmin}}^{\mathrm{2}}$. The so-called mass bias parameter, b, accounts for any bias between the estimated mass and the true halo mass, M_500,true, for example from violation of hydrostatic equilibrium and/or X-ray instrument calibration. Planck SZ cluster counts require high values of the mass bias, (1 − b) ~ 0.6−0.7, for consistency with the base ΛCDM model favored by Planck’s measurements of the primary CMB anisotropies (see Planck Collaboration XX 2014; Planck Collaboration XXIV 2015, and references therein).

The $\begin{array}{c}{}_{\mathrm{􏽢}}\\ {\mathit{M}}_{\mathrm{500}}\end{array}\mathrm{-}\begin{array}{c}{}^{\mathrm{􏽢}}\\ {\mathit{Y}}_{\mathrm{500}}\end{array}$ relation can be derived from the last two equations: $\frac{\begin{array}{c}\mathrm{􏽢}\\ {\mathit{M}}_{\mathrm{500}}\end{array}}{{\mathrm{10}}^{\mathrm{13}}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}}\mathrm{=}{\left(\frac{\mathit{h}}{\mathrm{0.7}}\right)}^{\mathrm{1}\mathit{/}\mathit{\alpha }\mathrm{-}\mathrm{1}}{\left[\frac{\begin{array}{c}\mathrm{􏽢}\\ {\mathit{Y}}_{\mathrm{500}}\end{array}}{\mathrm{2.00}\mathrm{\times }{\mathrm{10}}^{-6}\hspace{0.17em}{\mathrm{arcmin}}^{\mathrm{2}}}\right]}^{\mathrm{1}\mathit{/}\mathit{\alpha }}\mathrm{·}$(24)

3.4. Statistics

3.4.1. Average values for the full QSO sample

As described before, the individual QSO signals are expected to be faint, so we need to average them over the population. Here, we detail the computation of this average. We adopt the generic notation ŵ_λ for $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{A}}_{\mathit{\lambda }}\end{array}$, $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{M}}_{\mathrm{dust}}\end{array}$, $\begin{array}{c}{\mathrm{􏽢}}^{}\\ {\mathit{L}}_{\mathrm{synch}}\end{array}$, $\begin{array}{c}{\mathrm{􏽢}}^{}\\ {\mathit{Y}}_{\mathrm{500}}\end{array}$ or $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{F}}_{\mathit{i}}\end{array}$, arranging them into a vector $\hat{{w}}$, and assume that they are Gaussian distributed.

If these different components of $\hat{{w}}$ are not correlated, we employ the usual inverse-variance weighted average: $\begin{array}{ccc}& & \mathrm{⟨}\mathit{ŵ}\mathit{\lambda }\mathrm{⟩}\mathrm{=}{\left[\sum _{\mathit{p}}\frac{\mathrm{1}}{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left[}{\mathrm{(}\mathit{ŵ}{\mathit{\lambda }}^{\mathrm{)}}}_{\mathit{p}}\mathrm{\right]}}\right]}^{-1}\left[\sum _{\mathit{p}}\frac{{\mathrm{(}\mathit{ŵ}{\mathit{\lambda }}^{\mathrm{)}}}_{\mathit{p}}}{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left[}{\mathrm{(}\mathit{ŵ}{\mathit{\lambda }}^{\mathrm{)}}}_{\mathit{p}}\mathrm{\right]}}\right]\mathit{,}\\ & & \mathit{\sigma }\mathrm{[}\mathrm{⟨}\mathit{ŵ}\mathit{\lambda }\mathrm{⟩}\mathrm{]}\mathrm{=}\sqrt{{\left[\sum _{\mathit{p}}\frac{\mathrm{1}}{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left[}{\mathrm{(}\mathit{ŵ}{\mathit{\lambda }}^{\mathrm{)}}}_{\mathit{p}}\mathrm{\right]}}\right]}^{-1}}\mathit{,}\end{array}$with p the index over the QSO catalogue.

In the correlated case, the average values are obtained according to $\begin{array}{ccc}& & \mathrm{⟨}{ŵ}\mathrm{⟩}\mathrm{=}{\left[\sum _{\mathit{p}}{{C}}^{-1}\mathrm{\left[}{\mathrm{(}\hat{{w}}\mathrm{)}}_{\mathit{p}}\mathrm{\right]}\right]}^{-1}\mathrm{·}\left[\sum _{\mathit{p}}{{C}}^{-1}\mathrm{\left[}{\mathrm{(}\hat{{w}}\mathrm{)}}_{\mathit{p}}\mathrm{\right]}\mathrm{·}{\mathrm{(}\hat{{w}}\mathrm{)}}_{\mathit{p}}\right]\mathit{,}\\ & & {C}\mathrm{[}\mathrm{⟨}{ŵ}\mathrm{⟩}\mathrm{]}\mathrm{=}{\left[\sum _{\mathit{p}}{{C}}^{-1}\mathrm{\left[}{\mathrm{(}\hat{{w}}\mathrm{)}}_{\mathit{p}}\mathrm{\right]}\right]}^{-1}\mathit{.}\end{array}$

3.4.2. χ²-test for the hypothesis

If the assumptions on the nature and number of components are correct, the measured flux densities, $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{F}}_{\mathit{i}}\end{array}$, and the reconstructed model flux densities, $\sum \begin{array}{c}{}_{\mathrm{􏽢}}\\ {\mathit{A}}_{\mathit{\lambda }}\end{array}\hspace{0.17em}{{}^{\mathrm{(}}\mathit{S}_{\mathit{\lambda }}{}^{\mathrm{)}}}_{\mathit{i}}\mathrm{\left(}\mathit{z}\mathrm{\right)}$, must be consistent within the uncertainties. We introduce the residual flux density for a given QSO as ${\mathit{R}}_{\mathit{i}}\mathrm{=}\begin{array}{c}\mathrm{􏽢}\\ {\mathit{F}}_{\mathit{i}}\end{array}\mathrm{-}\sum _{\mathit{\lambda }}\begin{array}{c}\mathrm{􏽢}\\ {\mathit{A}}_{\mathit{\lambda }}\end{array}\hspace{0.17em}{\mathrm{\left(}{\mathit{S}}_{\mathit{\lambda }}\mathrm{\right)}}_{\mathit{i}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{.}$(29)Correlation between the measured flux densities, $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {F}\end{array}$, and the model amplitudes, $\begin{array}{c}{{􏽢}}_{}\\ {A}\end{array}$, leads to the following expression for the variance of the residuals: $\begin{array}{ccc}\mathit{C}\mathrm{[}{R}{\mathrm{]}}_{\mathit{ij}}& & \mathrm{\equiv }\mathrm{⟨}{\mathit{R}}_{\mathit{i}}{\mathit{R}}_{\mathit{j}}\mathrm{⟩}\mathrm{-}\mathrm{⟨}{\mathit{R}}_{\mathit{i}}\mathrm{⟩}\mathrm{⟨}{\mathit{R}}_{\mathit{j}}\mathrm{⟩}\\ & \mathrm{=}& \mathit{C}\mathrm{[}\begin{array}{c}\mathrm{􏽢}\\ {F}\end{array}{\mathrm{]}}_{\mathit{ij}}\mathrm{-}\sum _{\mathit{\lambda ,\mu }}{\mathrm{(}{\mathit{S}}_{\mathit{\lambda }}\mathrm{)}}_{\mathit{i}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{C}\mathrm{\left[}\begin{array}{c}{􏽢}\\ {A}\end{array}{\mathrm{]}}_{\mathit{\lambda ,\mu }}{\mathrm{(}{\mathit{S}}_{\mathit{\mu }}\mathrm{)}}_{\mathit{j}}\mathrm{\right(}\mathit{z}\mathrm{)}\mathit{,}\end{array}$(30)where $\begin{array}{ccc}\mathit{C}\mathrm{[}\begin{array}{c}\mathrm{􏽢}\\ {F}\end{array}{\mathrm{]}}_{\mathit{ij}}& & \mathrm{\equiv }\mathrm{⟨}\begin{array}{c}\mathrm{􏽢}\\ {\mathit{F}}_{\mathit{i}}\end{array}\begin{array}{c}\mathrm{􏽢}\\ {\mathit{F}}_{\mathit{j}}\end{array}\mathrm{⟩}\mathrm{-}\mathrm{⟨}\begin{array}{c}\mathrm{􏽢}\\ {\mathit{F}}_{\mathit{i}}\end{array}\mathrm{⟩}\mathrm{⟨}\begin{array}{c}\mathrm{􏽢}\\ {\mathit{F}}_{\mathit{j}}\end{array}\mathrm{⟩}\\ & \mathrm{=}& {\mathit{N}}_{{\mathit{R}}_{\mathrm{500}}}{\left[\mathrm{\int }{\mathrm{d}}^{\mathrm{2}}\mathit{k}\frac{{\mathrm{(}{{\mathit{t}}_{\mathrm{qso}}^{\mathrm{\ast }}}^{\mathrm{)}}}_{\mathit{i}}\mathrm{\left(}{k}\mathrm{\right)}\hspace{0.17em}{\mathrm{(}{{\mathit{t}}_{\mathrm{qso}}}^{\mathrm{)}}}_{\mathit{j}}\mathrm{\left(}{k}\mathrm{\right)}\hspace{0.17em}{\mathit{P}}_{\mathit{ij}}\mathrm{\left(}{k}\mathrm{\right)}}{{\mathit{P}}_{\mathit{ii}}\mathrm{\left(}{k}\mathrm{\right)}\hspace{0.17em}{\mathit{P}}_{\mathit{jj}}\mathrm{\left(}{k}\mathrm{\right)}}\right]}^{-1}\mathit{,}\end{array}$(31)with N_R500 = ^∫_r<R500dr4πr²τ_qso(r). The residuals R_i are expected to be Gaussian distributed with zero mean. We employ Eqs. (27) and (28) with R, instead of ŵ, to compute the average value, and define the χ² for the average residual as ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\mathrm{⟨}{{R}\mathrm{⟩}}^{\mathit{t}}\mathrm{·}{{C}}^{-1}\mathrm{[}\mathrm{⟨}{R}\mathrm{⟩}\mathrm{]}\mathrm{·}\mathrm{⟨}{R}\mathrm{⟩}\mathit{.}$(32)The average residual vector ⟨ R ⟩ has seven components, one per frequency. We expect the χ² to follow Student’s law with seven degrees of freedom, illustrated with simulations in Sect. 4.2. The value of χ² depends on the number of components and their frequency dependance. In particular, the dust emission is parametrized by T_d and β_d, so the χ² also depends on their values. If we leave these two parameters free, the χ² would be expected to follow Student’s law with 7−2 = 5 degrees-of-freedom.

We study the average residual vector, ⟨ R ⟩, because it is the most useful quantity for identifying deviations between observed signals and model predictions (dust, dust + tSZ, dust + synchrotron or dust + tSZ + synchrotron). Individual ${\mathit{\chi }}_{\mathit{i}}^{\mathrm{2}}$ obtained from the residuals R_i can also be computed, but given the large measurement uncertainties, study of the individual ${\mathit{\chi }}_{\mathit{i}}^{\mathrm{2}}$ distributions do not efficiently discriminate models.

3.5. Marginalization over dust properties

The previous section considered computation of ⟨ ŵ_λ ⟩ and χ² for fixed dust properties, described by the parameters T_d and β_d. Although $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{A}}_{\mathrm{dust}}\end{array}$ is only weakly dependent on these parameters, the derived dust mass, $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{M}}_{\mathrm{dust}}\end{array}$ (Eq. (20)), is more sensitive to them. We use Powell’s algorithm to minimize the χ² relative to T_d, β_d to find their best-fit values. We also marginalize over T_d, β_d when determining the physical properties of the QSO environment, $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{M}}_{\mathrm{dust}}\end{array}$, $\begin{array}{c}{\mathrm{􏽢}}^{}\\ {\mathit{L}}_{\mathrm{synch}}\end{array}$ and $\begin{array}{c}{\mathrm{􏽢}}^{}\\ {\mathit{Y}}_{\mathrm{500}}\end{array}$.

By “data” in the following, we will mean the combination of the Planck maps and the BOSS QSO positions. The ⟨ ŵ_λ ⟩ being Gaussian distributed (Sec. 3.4.1), we write $\mathit{P}\mathrm{\left(}{\mathit{ŵ}}_{\mathit{\lambda }}\mathrm{\right|}{\mathit{T}}_{\mathit{d}}\mathit{,}{\mathit{\beta }}_{\mathit{d}}\mathit{,}\mathrm{data}\mathrm{)}\mathrm{\propto }\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{{\mathit{ŵ}}_{\mathit{\lambda }}\mathrm{-}\mathrm{⟨}\mathit{ŵ}\mathit{\lambda }\mathrm{⟩}}{\mathit{\sigma }\mathrm{[}\mathrm{⟨}\mathit{ŵ}\mathit{\lambda }\mathrm{⟩}\mathrm{]}}\right)}^{\mathrm{2}}\right)\mathrm{·}$(33)With $\mathit{P}\mathrm{\left(}\mathrm{data}\mathrm{\right|}{\mathit{T}}_{\mathrm{d}}\mathit{,}{\mathit{\beta }}_{\mathrm{d}}\mathrm{)}\mathrm{\propto }\exp \left(\mathrm{-}\frac{{\mathit{\chi }}^{\mathrm{2}}}{\mathrm{2}}\right)$, Bayes’ theorem gives $\mathit{P}\mathrm{\left(}{\mathit{T}}_{\mathrm{d}}\mathit{,}{\mathit{\beta }}_{\mathrm{d}}\mathrm{\right|}\mathrm{data}\mathrm{)}\mathrm{\propto }\mathit{P}\mathrm{(}{\mathit{T}}_{\mathrm{d}}\mathit{,}{\mathit{\beta }}_{\mathrm{d}}\mathrm{)}\mathrm{\times }\exp \left(\mathrm{-}\frac{{\mathit{\chi }}^{\mathrm{2}}}{\mathrm{2}}\right)\mathrm{·}$(34)Assuming a flat prior on P(T_d,β_d), we compute the probability of ŵ_λ, T_d, and β_d given the data as, $\mathit{P}\mathrm{\left(}{\mathit{ŵ}}_{\mathit{\lambda }}\mathit{,}{\mathit{T}}_{\mathrm{d}}\mathit{,}{\mathit{\beta }}_{\mathrm{d}}\mathrm{\right|}\mathrm{data}\mathrm{)}\mathrm{\propto }\mathit{P}\mathrm{(}{\mathit{ŵ}}_{\mathit{\lambda }}\mathrm{\left|}{\mathit{T}}_{\mathrm{d}}\mathit{,}{\mathit{\beta }}_{\mathrm{d}}\mathit{,}\mathrm{data}\mathrm{\right)}\mathrm{\times }\mathit{P}\mathrm{\left(}{\mathit{T}}_{\mathrm{d}}\mathit{,}{\mathit{\beta }}_{\mathrm{d}}\mathrm{\right|}\mathrm{data}\mathrm{)}\mathit{.}$(35)Averages are now calculated by marginalizing over T_d and β_d and via an integration over ŵ_λ: ${\overline{\mathit{w}}}_{\mathit{\lambda }}\mathrm{=}\${\mathit{ŵ}}_{\mathit{\lambda }}\mathit{P}\mathrm{\left(}{\mathit{ŵ}}_{\mathit{\lambda }}\mathit{,}{\mathit{T}}_{\mathrm{d}}\mathit{,}{\mathit{\beta }}_{\mathrm{d}}\mathrm{\right|}\mathrm{data}\mathrm{)}\mathrm{d}{\mathit{T}}_{\mathrm{d}}\mathrm{d}{\mathit{\beta }}_{\mathrm{d}}{\mathit{ŵ}}_{\mathit{\lambda }}\mathit{,}$(36)with variance ${\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left[}{\overline{\mathit{w}}}_{\mathit{\lambda }}\mathrm{\right]}\mathrm{=}\$\mathrm{(}{\mathit{ŵ}}_{\mathit{\lambda }}\mathrm{-}{\overline{\mathit{w}}}_{\mathit{\lambda }}{\mathrm{)}}^{\mathrm{2}}\mathit{P}\mathrm{(}{\mathit{ŵ}}_{\mathit{\lambda }}\mathit{,}{\mathit{T}}_{\mathrm{d}}\mathit{,}{\mathit{\beta }}_{\mathrm{d}}\mathrm{\left|}\mathrm{data}\mathrm{\right)}\mathrm{d}{\mathit{T}}_{\mathrm{d}}\mathrm{d}{\mathit{\beta }}_{\mathrm{d}}\mathrm{d}{\mathit{ŵ}}_{\mathit{\lambda }}\mathit{.}$(37)

4. Simulations

We validate our methods by injecting simulated QSOs into the Planck maps and re-extracting their properties. The simulations are described in Sect. 4.1. Section 4.2 shows that our χ² statistic carries the expected number of degrees-of-freedom. We demonstrate that the MMFs properly recover the properties of the mock catalogs in Sect. 4.3, and in Sect. 4.4 we examine the sensitivity our results to the dust model.

4.1. Mock catalogs and injected maps

We build mock catalogs directly from the original BOSS sample, keeping the same QSO redshift distribution (Fig. 1) but drawing the sky positions at random outside the Planck point-source mask. We fix the QSO host halo mass M₅₀₀ = 2.7 × 10¹³M_⊙, close to the value found from the data in Sect. 5.4. We then compute the expected tSZ flux for each host from (z,M₅₀₀) using Eq. (23).

We assign to each host halo a constant dust mass (M_dust)_input = 2.5 × 10⁸M_⊙ emitting with a modified blackbody spectrum at (T_d)_input = 25 K, (β_d)_input = 2.5. When specified, we also implement variations in the dust temperature (T_d)_input using a Gaussian distribution with mean (T_d)_input = 25 K and standard deviation (σ_T)_input = 5K.

With this model we compute the signal of each QSO in the Planck bands (from tSZ and dust) and we inject them directly into the maps using our source profile, τ_qso, convolved with the individual channel beams. We artificially fix θ_s to 0.27 arcmin for the injected profile to avoid possible biases due to the mismatch between the extraction profile and the halo profile at low redshift. The sensitivity of our result to this assumption is studied in Sect. 5.5.

4.2. Degrees of freedom

We test the number of degrees-of-freedom for the χ² statistic defined in Sect. 3.4.2 by injecting into the Planck maps a mock catalogue with only dust emission from the QSOs and then extracting the signal with the dust-only MMF at fixed T_d = (T_d)_input and β_d = (β_d)_input. We bin our catalogue into 148 sub-catalogs of 2000 QSOs each to compute 148 independent χ² values. The cumulative distribution of these values is given in the top panel of Fig. 2 as the solid red line. Student’s cumulative distribution with varying degrees-of-freedom (d.o.f.) are shown as the dashed lines. The red line follows Student’s law with seven degrees-of-freedom, corresponding to the seven Planck maps. When leaving T_d and β_d free and choosing their best-fit values as the point minimizing the χ², we obtain the solid blue line with 7−2 = 5 degrees-of-freedom, as expected.

Fig. 2 Cumulative χ2 distributions when Td and βd are fixed in the filtering (solid red curve) or when they are left free (solid blue curve). Top panel: distribution for a dust-only filtering on a dust-only mock catalogue. Bottom panel: distribution for a dust + tSZ filtering on a dust + tSZ mock catalogue. In both cases (dust-only or dust + tSZ), fixing Td and βd leads to a χ2 distribution with seven degrees-of-freedom (corresponding to our seven maps), while leaving them free reduces the number of degrees-of-freedom to five.

We next inject a mock catalogue with both dust and tSZ emission and re-extract the QSO signals using a dust + tSZ filter to compute the χ². The results are shown in the bottom panel of Fig. 2.

Fixing T_d and β_d to the input values also leads to a χ² with seven d.o.f. Leaving the two dust parameters free lowers the χ² to five d.o.f., as for the dust-only case. Although we are searching for an additional component in the data between the first and the second tests, the χ² conserves 7/5 d.o.f. for (T_d, β_d) fixed/free respectively. This is because the correlations between our residuals are taken into account in the covariance matrix C [ ⟨ R ⟩ ] in Eq. (32). In other words, C [ ⟨ R ⟩ ] changes when considering the dust-only MMF or the dust + tSZ MMF, so the χ² statistic defined in Eq. (32) does not depend on the number of components assumed for the MMF.

4.3. Extraction of simulated QSO

We first consider null tests where the signals are extracted at random positions in the Planck maps. Table 1 shows the output for three filters (dust, dust + tSZ, dust + synch) assuming fixed values of T_d = 25 K and β_d = 2.5 in the extraction. The output dust masses, SZ masses and synchrotron luminosities are compatible with zero, as expected, and the χ² do not show any significant deviation from Student’s law with the expected seven d.o.f.

Results of the extraction on mock catalogs injected into the Planck maps are summarized in Table 2. The recovered dust masses, SZ masses, and synchrotron luminosities are marginalized over T_d and β_d as described in Sect. 3.5.

For the injected dust + tSZ catalogue, the dust + tSZ filter recovers a unbiased estimate of the input dust temperature and spectral index. The input dust mass and the halo mass from the tSZ are also recovered at S/N ~ 8 and S/N ~ 14. The dust mass recovered using the dust-only filter is biased by + 3.4σ due to contamination of the dust emission by the tSZ. The recovered dust spectral index is also biased by −3.7σ. Adopting the dust + synch filter leads to a significant detection of negative synchrotron luminosity at + 4.6σ, which mimics the tSZ effect at low frequencies.

The dust-only extraction provides a reduced χ²/ 5 value of 12.1, larger than the dust + tSZ value of 1.1, indicating that the dust + tSZ model must be preferred over the dust-only model. The dust + synch model provides a better χ²/ 5 than the dust-only case (4.0), but the significantly negative value for the synchrotron luminosity discards the model. The (T_d,β_d) contours for the three filters are displayed in Fig. 3. The dust + tSZ filter shows contours enclosing the input values while the two other filters are biased.

We also implemented the triple dust + tSZ + synch filter and applied it to the dust + tSZ simulation. The triple filter transforms part of the tSZ emission into a negative (i.e., unphysical) synchrotron emission, as shown in Table 2. It recovers biased T_d, β_d, and M_dust. The reduced χ²/ 5 increases with respect to the dust + tSZ filter, showing that the dust + tSZ fit must be preferred over the dust + tSZ + synch fit. We show in Sect. 5 that the triple filter exhibits the same behavior on the data. We will thus present our main results using the dust, the dust + tSZ, and the dust + synch filters, depending on the value of the reduced χ²/ 5.

Fig. 3 Contours of χ2 for the mock QSO dust + tSZ catalogue using the dust, the dust + tSZ, and the dust + synch filters. The solid black line indicates the 1σ deviation from the minimal χ2 and the dotted black line the 2σ deviation. The horizontal and vertical dashed lines mark the position of the input Td and βd.

Fig. 4 Contours of χ2 for the mock QSO dust (σT = 5 K) + tSZ catalogue, obtained using the dust(σT = 2 K) + tSZ, the dust (σT = 5 K) + tSZ, and the dust (σT = 10 K) + tSZ filters. The solid black line indicates the 1σ deviation from the minimal χ2, and the dotted black line the 2σ deviation. The horizontal and vertical dashed lines mark the position of the input Td and βd.

4.4. Possible systematics due to the dust characteristics

We test the sensitivity of the extraction method to improper modeling of the dust by injecting mock QSOs with a dust temperature following a Gaussian distribution, described in Sect. 4.1, and then recovering the parameters with the different MMF. Results are shown in the bottom part of Table 2. The dust (σ_T = xK) + tSZ filters are a dust + tSZ filter for which a Gaussian scatter of x Kelvin in the dust temperature is added. As expected, the dust (σ_T = 5 K) + tSZ filter provides an unbiased estimate of the parameters. The derived dust and SZ masses are not very sensitive to the dispersion in the dust temperature: the dust + tSZ filter returns a satisfactory estimate of the two quantities. However, for this latter filter the recovered dust parameters T_d and β_d are significantly biased. This sensitivity is illustrated in Fig. 4 displaying the (T_d, β_d) contours for the three last filters of Table 2.

5. Results

We apply the method described in Sect. 3 to the real Planck data. In Sect. 5.1, we show the average total flux of the QSO sample across the Planck channels. We then estimate the contribution of the various components in Sects. 5.2 (dust), 5.3 (synchrotron), and 5.4 (hot gas).

5.1. The QSO signal at Planck frequencies

We first estimate the mean flux density, $\mathrm{⟨}\begin{array}{c}{}_{\mathrm{􏽢}}\\ {F}\end{array}\mathrm{⟩}$, of the 291 165 QSOs selected in Sect. 2.1 by averaging the individual flux densities in each band as extracted with a single frequency matched filter centered on the QSO positions. The average total flux is shown as the black diamonds in Fig. 5 (main panel and inset). It is strictly positive in all Planck bands above 100 GHz, and positive but compatible with zero at 70 GHz. The signal continuously increases from the lowest to the highest frequencies: emission in the direction of the QSO is dominated by dust. Since none of the frequencies below 217 GHz presents a negative flux, there is no obvious indication of strong tSZ emission.

Fig. 5 Average flux density for the QSO population across Planck channels (black diamonds), and dust signal estimated with the dust-only filter (red triangles) for Td and βd at the minimum χ2. The bottom panel shows the average residuals.

Using the dust-only MMF with T_d and β_d adjusted to minimize the χ² defined in Eq. (32) (T_d = 18.7 K, β_d = 2.79), we compute the mean dust flux density $\begin{array}{c}{\mathrm{􏽢}}_{}\\ {\mathit{A}}_{\mathrm{dust}}\end{array}$ and plot the resulting mean dust spectrum as the red triangles in Fig. 5. Dust accounts for essentially all the observed signal. The residual (difference between the black diamonds and the red triangles) is shown in the bottom panel of the figure, and is consistent with zero except at 353 GHz. This deviation explains why the reduced χ²/ 5 highlighted in the first gray line of Table 3 is not good (χ²/ 5 = 5.2).

By marginalizing over T_d and β_d, we obtain an average dust mass, $\overline{{\mathit{M}}_{\mathrm{dust}}}\mathrm{=}\mathrm{(}\mathrm{0.84}\mathrm{\pm }\mathrm{0.07}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{\mathrm{8}}{\mathit{M}}_{\mathrm{\odot }}$. This value is ~25 times smaller than the estimated dust mass in optically-selected clusters, around 2 × 10⁹M_⊙ (Gutiérrez & López-Corredoira 2014). It is also five times smaller than the lowest estimated QSO dust mass of the Beelen et al. (2006) sample. We note, however, that our average mass is dominated by low redshift objects (z ~ 0.5), as shown in Fig. 7, while the QSO sample studied in Beelen et al. (2006) resides at higher redshift.

The recovered dust mass also depends strongly on the assumed value for T_d and β_d. If we fix β_d = 1.6, as in Beelen et al. (2006), we find $\overline{{\mathit{M}}_{\mathrm{dust}}}\mathrm{=}\mathrm{(}\mathrm{1.70}\mathrm{\pm }\mathrm{0.08}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{\mathrm{8}}{\mathit{M}}_{\mathrm{\odot }}$, in better agreement with the above-cited study, but the χ²/ 5 increases to 16.3, as shown in Table 3. This result demonstrates that the bulk of the BOSS QSO population contains significantly less dust than QSOs examined in these previous studies.

For the dust temperature, we find a marginalized T_d = 19.1 ± 0.8 K, close to the dust temperature of normal galaxies (e.g., Clemens et al. 2013). This is significantly lower than the temperature found by Beelen et al. (2006) and Dai et al. (2012), ranging between 33 and 55 K, and 18.1 and 79.7 K respectively. We note that Dai et al. (2012) use the two component model from Blain et al. (2003) with β_d = 2, and Beelen et al. (2006) fix β_d = 1.6 for the majority of their QSOs. We find β_d = 2.71 ± 0.13, significantly higher than these values. When fixing β_d = 1.6 in our analysis, our result moves along the degeneracy line in the (T_d, β_d)-plane to reach T_d = 28.2 ± 0.5 K, in better agreement with the values published in the previously mentioned analyses. Restricting this test to high redshift QSOs leaves our conclusions unchanged (line 6 of Table 3).

Fig. 6 Top panel: average Td as a function of redshift z. Middle panel: average βd as a function of redshift z. There is a clear evolution between the low-z (z< 1.5) , high-temperature and low-β and the high-z, low-temperature and high-β QSO populations. Bottom panel: average dust ${\mathit{T}}_{\mathrm{d}}\mathrm{\times }{\mathit{\beta }}_{\mathrm{d}}^{\mathrm{0.6}}$ (combination describing the degeneracy line) as a function of redshift z. There is no evidence for evolution of the dust properties between z = 0 and 1.5. The product increases between z = 1.5 and 4.

Fig. 7 Average dust flux density after marginalization over Td and βd (black diamonds) as a function of redshift z. The corresponding dust mass is shown as blue triangles (Eq. (20)). The extraction is achieved using the dust-only filter.

5.2. Redshift dependence of the dust properties

In order to adapt the MMF to a possible variation of the QSO dust properties across redshift, we divide our QSO sample into nine regular bins of size Δz = 0.5, shown as the vertical dashed lines in Fig. 1. We include in the last bin the QSOs with 4 <z< 5 because there are not enough statistics to fill a tenth bin.

Fig. 8 Top panel: average dust flux density (black diamonds) for QSOs with a radio counterpart in the FIRST catalogue. The signal is extracted with the dust + synch filter and marginalized over Td and βd. The corresponding dust mass is shown as the blue triangles. Bottom panel: average synchrotron flux density at 100 GHz observed frequency (black diamonds) and the corresponding monochromatic luminosity (magenta triangles).

We minimize the χ² for each bin individually and plot the value for the combination ${\mathit{T}}_{\mathrm{d}}\mathrm{\times }{\mathit{\beta }}_{\mathrm{d}}^{\mathrm{0.6}}$ in Fig. 6. This quantity follows the degeneracy line for the dust parameters. There is no evidence for significant evolution across redshift between z = 0 and 2. For z> 2, ${\mathit{T}}_{\mathrm{d}}\mathrm{\times }{\mathit{\beta }}_{\mathrm{d}}^{\mathrm{0.6}}$ increases with redshift. If we suppose that β_d remains constant, this implies that the dust was hotter at high redshift and cooled until z = 2 before stabilizing at the current values.

The dust flux density extracted using the dust-only MMF is shown in Fig. 7 as the black diamonds. It varies with redshift, increasing from z = 0 to z = 2 (~8 mJy within R₅₀₀ at 857 GHz) and decreasing with redshift for z> 2. The corresponding dust mass (Eq. (20)) is shown by the blue diamonds and follows the same trend as the flux density, reaching a maximum of ~5 × 10⁸M_⊙ at z ~ 2. This is of the same order as dust masses determined from high signal-to-noise observations of individual QSOs (e.g., Beelen et al. 2006; Wang et al. 2007).

The IR luminosity, or equivalently the dust mass (Eq. (20)), is a tracer of star formation, and it is remarkable that the dust mass evolution in Fig. 7 follows that of the general star formation rate with a peak around z ~ 2 (see, e.g., Fig. 9 in Madau & Dickinson 2014). The accretion rate onto AGN SMBHs follows the same trend (Hopkins et al. 2007; Shankar et al. 2009; Aird et al. 2010; Delvecchio et al. 2014). This is strong evidence that star formation in QSO environments is typical of the general galaxy population, and also that it is linked to the QSO central engines. Our result clearly shows this trend for a large, representative sample of QSOs, and is consistent with the study by Wang et al. (2015) of the correlation between QSOs and Herschel measurements of the cosmic infrared background.

Table 3 also contains results for the dust + tSZ, dust + synch, and the dust + tSZ + synch filters on the full QSO population (0.1 <z< 6.5). All three of the multiple component filters increase the χ² with respect to the dust-only filter.

5.3. Synchrotron emission

We now select only QSOs with at least one FIRST³ radio source counterpart at 1.4 GHz using the FIRST_MATCH keyword in the DR12 catalogue (FIRST_MATCH = 1). This provides a sub-catalogue of 9983 QSOs. Applying the dust + synch filter to this sub-catalogue using the same redshift binning, we plot the average dust and synchrotron flux densities, and , in Fig. 8. Dust emission from the sub-catalogue is compatible with that of the full catalogue, although with larger uncertainties due to the smaller sample (blue diamonds).

This population of radio-loud QSOs exhibits strong synchrotron emission, as shown in the lower panel. The monochromatic synchrotron luminosity (green diamonds) increases from 0 to reach 0.2 L_⊙ Hz^-1 at the higher redshifts; We note that this is the luminosity at the rest-frame frequency 100 GHz, assuming our power-law in Eq. (21).

Figure 9 compares the synchrotron emission determined from Planck (expressed at 100 GHz) to the FIRST measurements at 1.4 GHz. For this comparison, we bin-average the FIRST signal with the same MMF weights used to extract the synchrotron signal. The emission seen by Planck steadily decreases with redshift relative to the signal measured by FIRST, demonstrating a variation in synchrotron index with redshift.

We show the effective spectral index (for an assumed power law) between FIRST and Planck in Fig. 10. The spectrum is quite flat at low redshift with values for the spectral index that are much smaller than our adopted value of 0.7. However, these are values describing the emission over a very large frequency range, while our adopted value is assumed to describe the emission around 100 GHz.

We also see that the effective spectral index steepens with redshift. This could be due to evolution intrinsic to the sources, or more simply the effect of curvature in the synchrotron spectrum. It is important to note in this light that although our synchrotron luminosity measurements (Eq. (21)) are referenced to 100 GHz rest-frame, assuming our adopted power law, the effective index shown in Fig. 10 is taken between the observed 1.4 and 100 GHz bands. In other words, it is the effective index between rest-frame frequencies of 1.4(1 + z) and 100(1 + z) GHz. The apparent evolution with redshift in the figure could therefore be steepening in the synchrotron emission with frequency that would be expected given the greater energy loss of the higher energy electrons contributing to the signal as we move up in rest-frame frequency.

Fig. 9 Average FIRST flux density at 1.4 GHz observed frequency compared to the average synchrotron flux density at 100 GHz observed frequency as a function of redshift z. The magenta triangles correspond to the black diamonds in the bottom panel of Fig. 8.

Fig. 10 Effective spectral index of power law in frequency between 1.4 GHz (FIRST) and 100 GHz (Planck) observation frequencies.

Fig. 11 Top panel: average dust flux density (black diamonds) of the radio-quite QSO subsample, i.e., without a radio counterpart in the FIRST catalogue. The signal is obtained with the dust + tSZ filter and is marginalized over Td and βd. The corresponding dust mass is shown by the blue triangles. Bottom panel: average Compton parameter (black diamonds) and corresponding integrated Compton parameter (red triangles) as a function of redshift.

5.4. Hot halo gas

The strong synchrotron emission of QSOs with FIRST counterparts complicates extraction of the tSZ signal from the hot gas in the host halo. We therefore construct a sub-sample without FIRST counterparts (FIRST_MATCH = 0) of 252 888 QSOs and apply the dust + tSZ filter. Results are shown in Fig. 11, and summarized in the bottom segment of Table 3 for this and other filters.

Applied to the same redshift range, the other filters yield larger χ² values. In particular, the dust + synch and dust + tSZ + synch filters show the same behavior as on simulations: they convert part of the tSZ signal into negative synchrotron luminosity. The dust(σ_T = x K) + tSZ filters presented in Sect. 4.4 also increase the χ², favoring a low scatter in the dust temperature.

The filter does not detect the tSZ effect at z< 2.5. The tSZ signal at z< 1.5 is slightly negative at ~2σ, pointing towards a possible contamination of the extracted signal at low z by a weak synchrotron emission. The filter does show a clear signal rising with redshift over the range 2.5 <z< 4. The global significance of the signal in the latter redshift range is 7.4σ when expressed in terms of the population mean, $\overline{{\mathit{Y}}_{\mathrm{500}}}$.

The dust + tSZ + synch filter applied on the radio-loud sub-catalogue, i.e., with at least one FIRST source (third segment of Table 3), provides a slightly better fit (χ²/ 5 = 0.3) than the dust + synch filter (χ²/ 5 = 0.4) over the 2.5 <z< 4 range. The recovered SZ signal remains compatible with the value derived from the radio-quiet sub-catalogue, but the uncertainties are large with the signal appearing at only 1.5σ.

When applying the tSZ-only filter to the 2.5 <z< 4 population, we find SZ-based masses that are twice as large as the mass provided by the dust + tSZ filter; the fit, however, is quite bad at χ²/ 7 > 100. (We note that there are seven degrees-of-freedom in this case because the dust temperature and emissivity do not come into play.) This simply reflects the fact that a tSZ-only signal is a poor description of the spectrum of the source. Details of these results are given in Table 4.

We stacked the ROSAT All-Sky Survey maps⁴ of the radio-quiet QSO sub-sample to find a significant signal, corresponding to a population mean luminosity of L₅₀₀ = (3.30 ± 0.30) × 10⁴⁴erg / s in the [0.1–2.4] keV band, the error taken as the standard deviation of the stacked maps. Pratt et al. (2009) established a power-law scaling relation between X-ray luminosity and halo mass based on analysis of a representative sample of X-ray galaxy clusters employing Bivariate Correlated Errors and intrinsic Scatter estimators with a minimization of the residuals in X-ray luminosity (BCES (YlX)). If we adopt their Malmquist-bias corrected relation between L₅₀₀, in this energy band, and mass (line 6 of their Table B.2), we infer a mean halo mass of M₅₀₀ = 5.87 × 10¹³M_⊙. This is 2.3 times higher than our SZ-based mass; alternatively, it implies an expected $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{Y}}_{\mathrm{500}}\end{array}\mathrm{=}\mathrm{(}\mathrm{51.4}\mathrm{\pm }\mathrm{4.6}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-6}{\mathrm{arcmin}}^{\mathrm{2}}$, or 4.7 times larger than our measurement. In other words, the tSZ and X-ray signals do not follow the empirical relation for thermal halo gas emission seen at low redshift (Planck Collaboration XI 2011). This is not a surprise given that we expect a non-negligible contribution to the X-ray emission from the central engine of the QSO itself.

5.5. Sensitivity to the assumed size θ_s

We fix θ_s = 0.27 arcmin throughout our analysis. The value corresponds to a halo mass of M₅₀₀ = 10¹³M_⊙ at z = 2. To test the robustness of our results to this assumption, we fix alternatively the size to θ_s = 0.12 and θ_s = 0.57 arcmin, corresponding to halos of mass M₅₀₀ = 10¹²M_⊙ and M₅₀₀ = 10¹⁴M_⊙, respectively, at z = 2. Results are shown in Table 5. Changing the filtering scale slightly increases the χ², and only weakly impacts the derived parameters (<0.5σ), leaving our results essentially unchanged.

6. Discussion

Submillimeter observations are a valuable source of information on the physical conditions in QSO environments because they probe star formation in and around the hosts and energy injection into the surrounding medium by the SMBHs. Due to atmospheric absorption, this waveband is however difficult to access from the ground, and this has limited studies to relatively small samples of several tens of objects (Omont et al. 2001, 2013; Isaak et al. 2002; Priddey et al. 2003; Beelen et al. 2006; Wang et al. 2007).

The space-based platforms WMAP, Planck and Herschel now open this window for much more extensive studies. Taking advantage of Planck’s all-sky and wide-frequency coverage, we have determined several mean characterstics of the QSO population through binning measurements of dust, synchrotron and tSZ signals from the very large BOSS QSO sample.

6.1. Dust emission

We find that the QSO spectrum is dominated by graybody dust emission with a temperature of T_d = 19.1 ± 0.8 K and an emissivity spectral index of β_d = 2.71 ± 0.13 when averaged over the entire sample spanning the full redshift range, 0.1 <z< 5; these are the one-dimensional marginalized constraints for the dust-only filter that gives the best χ² (see Table 3). While there is no clear evidence that these quantities evolve with redshift at z< 2, there is a rise at higher redshifts in the product ${\mathit{T}}_{\mathrm{d}}{\mathit{\beta }}_{\mathrm{d}}^{\mathrm{0.6}}$, describing the degeneracy direction, at higher redshifts (Fig. 6).

Curiously, the temperature is notably lower than found in ground-based observations of high-redshift quasars, which find temperatures of 40–50 K, while our value is more typical of ordinary galaxies like the Milk Way. Our dust emissivity index is significantly steeper than the values of 1–2 for normal Galactic dust, and should be compared to the value of 1.6 found in the ground-based surveys. It is notably steeper than the limiting case of β = 2 for ideal insulators and conductors far from a resonance. These characteristics remain in the other sub-samples and component filters summarized in Table 3, apart the unusual behaviour of the radio-loud sub-sample with FIRST counterparts. Fixing β_d = 1.6, we obtain a somewhat higher mean temperature, T_dust = 28.2 ± 0.5 K, but at the cost of a much worse spectral fit according to the large increase in χ² (see Table 3).

We calculate dust mass using the simple prescription of Eq. (20), which can also be viewed as a measure of far-IR luminosity. An important result is our finding that the dust mass, a tracer of local star formation, evolves with redshift in the same way as the global cosmic star formation rate, peaking at z ~ 2 (see e.g. Fig. 9 in Madau & Dickinson 2014). Star formation in QSO environments apparently does not distinguish itself from others. It is also known, from measures of AGN luminosity functions, that the SMBH accretion rate follows the same evolutionary track (Hopkins et al. 2007; Shankar et al. 2009; Aird et al. 2010; Delvecchio et al. 2014).

6.2. Synchrotron emission

A subsample of our QSOs are radio loud with counterparts at 1.4 GHz in the FIRST catalogue. We detect synchrotron emission from this subsample with an monochromatic luminosity, referenced to rest frame frequency 100 GHz (adopting our power-law of Eq. (21)), increasing with redshift between z = 0 and z = 3. Comparing the FIRST and Planck flux densities, we find an effective spectral index between the two corresponding rest-frame frequencies that increases with redshift from α_s = 0.3 to α_s = 0.7. This steepening is likely the result of curvature in the spectrum as we observe progressively higher rest-frame frequencies.

6.3. tSZ and feedback

Restricting analysis to the radio-quite sub-sample, we detect the tSZ effect (with the dust + tSZ filter) at S/N ~ 7.4 over the redshift range 2.5 <z< 4, with an integrated Compton-y parameter of $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{Y}}_{\mathrm{500}}\end{array}\mathrm{=}\mathrm{(}\mathrm{10.83}\mathrm{\pm }\mathrm{1.46}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-6}{\mathrm{arcmin}}^{\mathrm{2}}$ (one-dimensional, marginalized over dust properties). We do not, however, find evidence of tSZ signal at the lower redshifts. At z< 1.5, we obtain a mean value of $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{Y}}_{\mathrm{500}}\end{array}\mathrm{(}\mathit{z}\mathit{<}\mathrm{1.5}\mathrm{)}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{4.93}\mathrm{\pm }\mathrm{2.48}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-6}$ arcmin², which we convert simply into an upper limit at two sigma above zero as $\begin{array}{c}\mathrm{˜}\\ {\mathit{Y}}_{\mathrm{500}}\end{array}\mathrm{(}\mathit{z}\mathit{<}\mathrm{1.5}\mathrm{)}\mathit{<}\mathrm{4.97}\mathrm{\times }{\mathrm{10}}^{-6}\hspace{0.17em}{\mathrm{arcmin}}^{\mathrm{2}} \mathrm{\left(}\mathrm{2}\mathit{\sigma }\mathrm{\right)}\mathit{.}$(38)Similarly, we find an upper limit at z< 2.5 of $\begin{array}{c}\mathrm{˜}\\ {\mathit{Y}}_{\mathrm{500}}\end{array}\mathrm{(}\mathit{z}\mathit{<}\mathrm{2.5}\mathrm{)}\mathit{<}\mathrm{2.93}\mathrm{\times }{\mathrm{10}}^{-6}\hspace{0.17em}{\mathrm{arcmin}}^{\mathrm{2}} \mathrm{\left(}\mathrm{2}\mathit{\sigma }\mathrm{\right)}\mathit{,}$(39)even stronger, as expected given the behaviour seen in Fig. 11.

Ruan et al. (2015) recently reported a detection of the tSZ signal by stacking the Hill & Spergel (2014) tSZ map, constructed from the Planck frequency maps, on 26 686 spectroscopic DR7 QSOs. For their low redshift bin at z< 1.5, they find $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{Y}}_{\mathrm{500}}^{\mathrm{c}}\end{array}\mathrm{=}\mathrm{(}\mathrm{74}\mathrm{\pm }\mathrm{30}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-6}$ arcmin², quoted in terms of the Compton-y parameter integrated along the line-of-sight cylinder. The conversion factor between the cylindrical and spherical integrated Compton-y parameters for our adopted profile is $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{Y}}_{\mathrm{500}}^{\mathrm{c}}\end{array}\mathrm{=}\mathrm{1.52}\mathrm{\times }\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{Y}}_{\mathrm{500}}\end{array}$, which translates the upper limit of Eq. (38) to $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{Y}}_{\mathrm{500}}^{\mathrm{c}}\end{array}\mathit{<}\mathrm{7.55}\mathrm{\times }{\mathrm{10}}^{-6}$ arcmin² (2σ), well below their detection level.

One possible explanation for the difference is dust contamination. Ruan et al. (2015) took particular pains to correct for any contamination by dust in the tSZ map. Our filters directly measure the individual local signals from dust and tSZ, as well as any covariance between them, allowing us to separate and marginalize over the dust influence when quoting our tSZ signal strengths. The importance of separating out the dust signal is easily appreciated from the results of applying the tSZ-only filter given in Table 4. Contamination by dust in the tSZ-only filter can produce tSZ signals ~3 times larger.

The tSZ signal directly measures the thermal energy stored in diffuse ionized gas in the QSO environment. We expect a signal from the intra-halo gas at the virial temperature of the QSO host halo, and there may be additional heating powered by feedback from the SMBH. Chatterjee et al. (2010) and Ruan et al. (2015) both interpreted their tSZ signals as evidence for feedback into the surrounding gas.

To look for this, we estimate the tSZ signal expected from intra-halo gas using the scaling relation between halo mass and tSZ signal in Eq. (23). We first establish a mean mass scale for the sample of 77 655 QSOs in 2.5 <z< 4 over the 14 555 sq. deg. SDSS footprint. Using the Tinker et al. (2008) mass function, we find the lower mass bound above which there are as many predicted halos as observed QSOs in this redshift range and for this sky area. We then calculate the mean mass by integrating over the mass function above this mass limit and over this redshift range. This gives us an estimate of the maximum mean mass characterizing the sample; QSOs could be hosted by lower mass halos if the QSO selection function is not 100% complete, as is most certainly the case.

We find a mean halo mass of M_500,true = 2.18 × 10¹³M_⊙, a value in reasonable agreement with Richardson et al. (2012), who found a QSO halo mass $\mathrm{14.}{\mathrm{1}}_{-6.9}^{\mathrm{+}\mathrm{5.8}}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}}{\mathit{h}}^{-1}{\mathit{M}}_{\mathrm{\odot }}$ at z ~ 3.2 for the SDSS sample. Our characteristic mass corresponds to a predicted tSZ signal of Y₅₀₀ = 8.47(1 − b)^1.78 × 10^-6 arcmin², only one sigma below our measured signal of Y₅₀₀ = (10.02 ± 1.34) × 10^-6 arcmin² if b = 0. This would not leave much room for additional heating by feedback. With the higher mass bias values suggested by Planck’s cluster counts, e.g., (1 − b) = 0.6, we instead have Y₅₀₀ = 3.4 × 10^-6 arcmin², implying that ~60% of the gas’s thermal energy is generated by feedback.

X-ray emission from the hot gas is also expected. Stacking the RASS maps at the QSO position, we find a luminosity L₅₀₀ = (3.30 ± 0.30) × 10⁴⁴ erg / s in the [0.1–2.4] keV band. The tSZ – X-ray luminosity scaling relation established by Planck Collaboration X (2011) says that this should correspond to an SZ signal of $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{Y}}_{\mathrm{500}}\end{array}\mathrm{=}\mathrm{(}\mathrm{51.4}\mathrm{\pm }\mathrm{4.6}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-6}{\mathrm{arcmin}}^{\mathrm{2}}$, much higher than our actual tSZ measurement, a conclusion that is independent of the mass bias.

The cluster halo gas scaling relations appear therefore to be violated in the QSO environment. This could be due to feedback from the QSO SMBH. It is likely that X-ray emission from the QSO itself also affects our X-ray measurements. Moreover, we are extrapolating the relation to redshifts much higher than those used to establish its validity. The redshift evolution of these relations is in fact not well known, and we have adopted self-similar evolution as a reasonable assumption given the lack of empirical constraints. Measurements of these relations on non-QSO systems at the same redshifts would provide a valuable comparison sample to judge if the QSO’s activity is in fact responsible for violating the scaling relations.

Finally, we note that interpretation of the measured tSZ signal is also clouded by possible contribution from gas in structures correlated with the QSO halos and falling with the relatively large effective beam of Planck, as suggested by Cen & Safarzadeh (2015). We leave careful modeling of this effect to future work.

7. Conclusion

Planck’s coverage of the SDSS survey area gives us the opportunity for in-depth study of quasar environments through submillimeter observations of optically selected QSOs. This spectral band provides valuable information on star formation and energy production by the QSO central engines. Filtering Planck’s seven spectral bands from 70–857 GHz with multi-component MMFs, we have measured dust, synchrotron and tSZ signals from the large BOSS QSO sample spanning redshifts out to 0.1 <z< 5. While measurements of individual QSOs are well below Planck’s sensitivity, the MMFs enable us to determine the population’s mean properties as a function of redshift.

Dust emission dominates the average QSO spectrum, and its evolution with redshift suggests that star formation in QSO environments follows the general cosmic star formation rate with a peak at z ~ 2. We did not detect a tSZ signal from QSOs at z< 2.5, but did find a signal at higher redshifts. The signal is clearly seen out to z ~ 4, the highest redshift to which the tSZ signal has been measured to date. A non-negligible fraction of the gas’s thermal energy could be supplied by QSO feedback, but interpretation is difficult because of modeling uncertainties associated with halo gas scaling relations. Finally, we observe synchrotron emission from a radio-loud sub-sample at rest-frame frequencies ν> 100 GHz and whose luminosity rises monotonically with redshift.

These results are based on a large sample of nearly 300 000 BOSS spectroscopic QSOs falling outside of the Planck mask. Our findings present a fresh and representative view of the average properties characterizing the QSO population.

Acknowledgments

Some of the results in this paper have been derived using the HEALPix (Górski et al. 2005) package. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. We thank T. Marriage for helpful discussions.

References

All Tables

All Figures

	Fig. 1 Redshift distribution of the BOSS DR12 QSO sample. The distribution presents two peaks, at z ~ 0.8 and z ~ 2.3, due to the QSO target selection process (see text). The regular binning used throughout this paper (Δz = 0.5) is shown as the vertical dotted lines.
In the text

Fig. 2 Cumulative χ2 distributions when Td and βd are fixed in the filtering (solid red curve) or when they are left free (solid blue curve). Top panel: distribution for a dust-only filtering on a dust-only mock catalogue. Bottom panel: distribution for a dust + tSZ filtering on a dust + tSZ mock catalogue. In both cases (dust-only or dust + tSZ), fixing Td and βd leads to a χ2 distribution with seven degrees-of-freedom (corresponding to our seven maps), while leaving them free reduces the number of degrees-of-freedom to five.

In the text

	Fig. 3 Contours of χ2 for the mock QSO dust + tSZ catalogue using the dust, the dust + tSZ, and the dust + synch filters. The solid black line indicates the 1σ deviation from the minimal χ2 and the dotted black line the 2σ deviation. The horizontal and vertical dashed lines mark the position of the input Td and βd.
In the text

	Fig. 4 Contours of χ2 for the mock QSO dust (σT = 5 K) + tSZ catalogue, obtained using the dust(σT = 2 K) + tSZ, the dust (σT = 5 K) + tSZ, and the dust (σT = 10 K) + tSZ filters. The solid black line indicates the 1σ deviation from the minimal χ2, and the dotted black line the 2σ deviation. The horizontal and vertical dashed lines mark the position of the input Td and βd.
In the text

	Fig. 5 Average flux density for the QSO population across Planck channels (black diamonds), and dust signal estimated with the dust-only filter (red triangles) for Td and βd at the minimum χ2. The bottom panel shows the average residuals.
In the text

Fig. 6 Top panel: average Td as a function of redshift z. Middle panel: average βd as a function of redshift z. There is a clear evolution between the low-z (z< 1.5) , high-temperature and low-β and the high-z, low-temperature and high-β QSO populations. Bottom panel: average dust ${\mathit{T}}_{\mathrm{d}}\mathrm{\times }{\mathit{\beta }}_{\mathrm{d}}^{\mathrm{0.6}}$ (combination describing the degeneracy line) as a function of redshift z. There is no evidence for evolution of the dust properties between z = 0 and 1.5. The product increases between z = 1.5 and 4.

In the text

	Fig. 7 Average dust flux density after marginalization over Td and βd (black diamonds) as a function of redshift z. The corresponding dust mass is shown as blue triangles (Eq. (20)). The extraction is achieved using the dust-only filter.
In the text

Fig. 8 Top panel: average dust flux density (black diamonds) for QSOs with a radio counterpart in the FIRST catalogue. The signal is extracted with the dust + synch filter and marginalized over Td and βd. The corresponding dust mass is shown as the blue triangles. Bottom panel: average synchrotron flux density at 100 GHz observed frequency (black diamonds) and the corresponding monochromatic luminosity (magenta triangles).

In the text

	Fig. 9 Average FIRST flux density at 1.4 GHz observed frequency compared to the average synchrotron flux density at 100 GHz observed frequency as a function of redshift z. The magenta triangles correspond to the black diamonds in the bottom panel of Fig. 8.
In the text

	Fig. 10 Effective spectral index of power law in frequency between 1.4 GHz (FIRST) and 100 GHz (Planck) observation frequencies.
In the text

Fig. 11 Top panel: average dust flux density (black diamonds) of the radio-quite QSO subsample, i.e., without a radio counterpart in the FIRST catalogue. The signal is obtained with the dust + tSZ filter and is marginalized over Td and βd. The corresponding dust mass is shown by the blue triangles. Bottom panel: average Compton parameter (black diamonds) and corresponding integrated Compton parameter (red triangles) as a function of redshift.

In the text