Constraining AGN feedback with the SZ effect

Constraints on AGN feedback from its Sunyaev-Zel’dovich imprint on the cosmic background radiation

Bjoern Soergel, Tommaso Giannantonio, George Efstathiou, Ewald Puchwein, Debora Sijacki
Institute of Astronomy & Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
Centre for Theoretical Cosmology, DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK
Universitäts-Sternwarte, Ludwig-Maximilians-Universität München, Scheinerstr. 1, 81679 München, Germany

We derive constraints on feedback by active galactic nuclei (AGN) by setting limits on their thermal Sunyaev-Zel’dovich (SZ) imprint on the cosmic microwave background (CMB). The amplitude of any SZ signature is small and degenerate with the poorly known sub-mm spectral energy distribution of the AGN host galaxy and other unresolved dusty sources along the line of sight. Here we break this degeneracy by combining microwave and sub-mm data from Planck with all-sky far-infrared maps from the AKARI satellite. We first test our measurement pipeline using the Sloan Digital Sky Survey (SDSS) redMaPPer catalogue of galaxy clusters, finding a highly significant detection () of the SZ effect together with correlated dust emission. We then constrain the SZ signal associated with spectroscopically confirmed quasi-stellar objects (QSOs) from SDSS data release 7 (DR7) and the Baryon Oscillation Spectroscopic Survey (BOSS) DR12. We obtain a low-significance () hint of an SZ signal, pointing towards a mean thermal energy of  erg, lower than reported in some previous studies. A comparison of our results with high-resolution hydrodynamical simulations including AGN feedback suggests QSO host masses of , but with a large uncertainty. Our analysis provides no conclusive evidence for an SZ signal specifically associated with AGN feedback.

Quasars: general – galaxies: active – galaxies: clusters: general – cosmic background radiation
pubyear: 2016pagerange: Constraints on AGN feedback from its Sunyaev-Zel’dovich imprint on the cosmic background radiationA.3

1 Introduction

Active galactic nuclei (AGN), powered by accretion of material onto a supermassive black hole (Lynden-Bell, 1969; Rees, 1984), are amongst the most luminous objects in the Universe.

AGN can deposit enormous quantities of energy into their surroundings, driving large-scale outflows (e.g. Blandford, 1990; Silk & Rees, 1998; Fabian, 1999; King, 2003). In fact, it is now widely accepted that AGN feedback plays a crucial role in our understanding of galaxy formation and evolution (Kauffmann & Haehnelt, 2000; Di Matteo et al., 2005; Croton et al., 2006; Bower et al., 2006; Sijacki et al., 2007; Schaye et al., 2010, 2015; Sijacki et al., 2015; Dubois et al., 2016). AGN feedback can also affect structure formation on scales larger than the immediate vicinity of the host galaxy. In clusters and groups of galaxies, AGN are believed to provide a heat source to balance radiative cooling of the hot gas in the central regions (e.g. Binney & Tabor, 1995; Churazov et al., 2001, 2002; Quilis et al., 2001; Sijacki et al., 2007; Puchwein et al., 2008; Dubois et al., 2010; McCarthy et al., 2011).

From a cosmological perspective, AGN feedback can alter the shape of the matter power spectrum on small scales, affecting cosmological parameters determined from cosmic shear measurements (e.g. van Daalen et al. 2011; Semboloni et al. 2011; Eifler et al. 2015). Accurate modelling of the AGN feedback energetics is therefore essential when simulating the formation of structure in our Universe, in both high-resolution simulations of individual haloes (e.g. Bhattacharya et al. 2008) and in simulations of cosmological volumes (e.g. Sijacki et al. 2007; Puchwein et al. 2008; Battaglia et al. 2010; Sijacki et al. 2015).

In this paper, we use quasi-stellar objects (QSOs) as tracers of AGN activity. Because of their high luminosities, it is now possible to construct large samples of QSOs from wide-field spectroscopic surveys (e.g. Pâris et al. 2016), making this subclass of AGN well suited to statistical studies of AGN feedback energetics.

Traditionally most of the evidence for AGN feedback has been obtained from radio and/or X-ray observations (see e.g. Burns, 1990; Boehringer et al., 1993; Carilli et al., 1994; Fabian et al., 2006; Forman et al., 2007; and Miley & De Breuck, 2008; Fabian, 2012 for reviews). It has now become possible to study gaseous outflows from individual AGN host galaxies at UV, optical and submillimetre wavelengths (e.g. Crenshaw et al., 2003; Rupke & Veilleux, 2011; Sturm et al., 2011; Maiolino et al., 2012; Tombesi et al., 2015).

Another technique for constraining AGN feedback, still in its nascent phase, is to exploit the thermal Sunyaev-Zel’dovich (SZ) effect (Sunyaev & Zeldovich, 1970, 1972). When passing through the hot ionized gaseous halo of the AGN host galaxy, a small fraction of the cosmic microwave background (CMB) photons is scattered off electrons with high thermal velocities. This leaves a characteristic, frequency-dependent signature in the CMB. If AGN feedback provides a contribution to the thermal energy of the halo gas, this will affect the pressure profile which can, in principle, be detected via the SZ signature (Chatterjee & Kosowsky, 2007; Bhattacharya et al., 2008; Chatterjee et al., 2008; Scannapieco et al., 2008).

In rich clusters of galaxies, the thermal energy from purely gravitational heating is large enough that the SZ signature of individual objects can be detected at high signal-to-noise by modern CMB experiments (e.g. Hasselfield et al. 2013; de Haan et al. 2016; Planck Collaboration 2016e). However, as we will show explicitly in Section 4, the SZ signal from AGN hosts is far too low to be detected in individual systems. It is therefore necessary to use statistical techniques to constrain the average AGN feedback energetics by, for example, cross-correlating CMB maps with optically selected QSOs. In this way, a tentative indication of a QSO feedback signal was obtained by Chatterjee et al. (2010) by cross-correlating CMB maps from the Wilkinson Microwave Anisotropy Probe (WMAP) with a photometric QSO catalogue from the Sloan Digital Sky Survey (SDSS).

More recently, Ruan et al. (2015) stacked data from a Compton- map (constructed by Hill & Spergel 2014 via internal linear combination of the Planck 2013 maps) at the positions of spectroscopically confirmed SDSS QSOs. They reported a significant detection corresponding to a thermal energy in the halo gas of  erg, significantly larger than suggested by the AGN feedback models adopted in hydrodynamical simulations (e.g. Sijacki et al. 2007; Battaglia et al. 2010). However, the interpretation of this signal as a detection of AGN feedback has been disputed by Crichton et al. (2016) and Verdier et al. (2016), who argued that it was primarily caused by dust from extragalactic sources, including the AGN host galaxy. An alternative explanation has been proposed by Cen & Safarzadeh (2015), who claimed that the Ruan et al. (2015) detection was caused by the cumulative, purely thermal, SZ signature from haloes correlated with the AGN.

The analysis by Verdier et al. (2016) used Planck 2015 maps and a multi-frequency multi-component matched filter to extract the SZ signal at the QSO positions. In addition to fitting to SZ, they included models for synchrotron and dust emission. They found that most of the correlated signal could be accounted for by dust, with no significant detection of an SZ signature for QSOs at . For QSOs with , however, they found evidence for a non-zero SZ signal, with an amplitude corresponding to a total thermal energy in the ionized gas of  erg (see Section 5.2 below), much lower than reported by Ruan et al. (2015). Verdier et al. (2016) concluded that it was not possible to unambiguously attribute the SZ signal seen in the high-redshift sample to AGN feedback because of the large uncertainties associated with the halo gas scaling relations.

Crichton et al. (2016) performed a stacking analysis, cross-correlating a radio-quiet subsample from the SDSS spectroscopic QSO catalogue with CMB maps from the Atacama Cosmology Telescope (ACT) and sub-mm maps from Herschel-SPIRE. This analysis also found evidence for a non-zero SZ signal associated with the QSOs, corresponding to a thermal energy of erg, or equivalently to a feedback efficiency of for a typical QSO activity timescale of  yr.

However, the above-mentioned feedback constraints rely on marginalisation over the amplitude and shape of the dust spectral energy distribution (SED). As significant degeneracies exist between dust and SZ parameters, especially when the fit is performed over a limited wavelength range, the SZ results can be biased by an inaccurate determination of the dust properties. Notably, the dust parameters deduced by Crichton et al. (2016) and Verdier et al. (2016) do not agree. Both studies fitted a modified blackbody dust spectrum, characterised by an amplitude, dust temperature , and spectral index : Verdier et al. (2016) found and for their full QSO sample () and their ‘dust only’ filter, with some evidence for a shift towards higher temperatures and lower values of for QSOs at redshifts . On the other hand, Crichton et al. (2016) found and for their QSO sample covering the redshift range .111For reference, the Planck analysis of Galactic dust emission (Planck Collaboration, 2014a) gives , over the high Galactic latitude sky (), parameters that are quite typical of nearby normal galaxies (Clemens et al., 2013), and close to those of the cosmic infrared background (Mak et al., 2016). As there is partial overlap in the QSO samples used in the Verdier et al. (2016) and Crichton et al. (2016) analyses, it seems unlikely that the disparity in the dust parameters is reflecting a difference in the QSO selection or their average physical properties. However, since dust emission dominates the cross-correlation signal in both analyses, any claim of a non-zero SZ signal is contingent on an accurate dust emission model.

In this paper, we use the Planck 2015 microwave and sub-mm maps in combination with all-sky far-infrared (FIR) data at from the AKARI satellite. With this additional spectral coverage on the falling, high-frequency part of the dust SED, it is possible to constrain the parameters of the dust emission model more robustly. This allows us to break the degeneracy between the dust parameters and the inferred SZ amplitude, providing a more reliable constraint on the average AGN feedback energetics.

Figure 1: CMB and FIR maps (in MJy/sr): we show here the seven Planck maps (LFI 70 GHz and HFI 100-857 GHz), as well as the AKARI 140 and 90 m maps. The 140 m channel was not used in the analysis due to larger calibration uncertainties, and is shown for illustration purposes only. All three maps in each row have the same colour scale, but the scale changes between the rows. The extent of the Galaxy in the various maps gives a good visual impression of the SED of Galactic dust; the SED of the cosmic infrared background (CIB) is relatively similar (e.g. Mak et al. 2016), but redshifted.
Band centre (GHz) 70 100 143 217 353 545 857 1,870 2,140 3,330 4,310 3,000
Band centre (m) 4,280 3,000 2,100 1,380 849 550 350 160 140 90 65 100
Beam FWHM () 13.31 9.68 7.30 5.02 4.94 4.83 4.64 1.47 1.30 1.05 4.3
Calibration uncertainty 1% 1% 1% 1% 1% 6% 6% 13.5%
Table 1: Summary of the map properties relevant for our analysis. Due to their large and partially unknown calibration uncertainties at low flux densities, we discard the AKARI bands displayed in grey, and only use the AKARI 90 m channel for our main analysis. (The resolution of the AKARI 160 m channel was not measured during the calibration process, but is expected to be close to the 140 m band; see Takita et al. 2015. The calibration uncertainty is assumed to be fully correlated between the Planck 70 to 353 GHz, 545 and 857 GHz, and IRIS 100 and AKARI 90 m bands, respectively.)
Figure 2: Instrumental bandpasses: we show the detector responses (normalized to a peak amplitude of unity) for the microwave and sub-mm data used in our analysis. Solid lines denote the Planck bandpasses, whereas dashed (dotted) lines show the AKARI (IRAS) transmission curves.

Our paper is organized as follows. Section 2 describes the various data products used in this work. Section 3 provides a detailed description of our analysis methods; in particular, we discuss the filtering of the input maps, the SED measurement, and our multi-component modelling and fitting of the QSO emission. Section 4 presents our main results. As an important consistency check of our methodology, we demonstrate using optically-selected galaxy clusters that our pipeline leads to an unambiguous detection of the SZ signal together with a correlated dust signal. We then report our results for the QSO cross-correlations. Section 5 discusses our results and compares then against high-resolution hydrodynamical simulations of both QSO host haloes and clusters that incorporate AGN feedback. Our conclusions are summarized in Section 6. In the course of this project, we further experimented extensively with IRAS maps over the frequency range and with AKARI maps at , and . Our reasons for excluding these maps from our main analysis are described in Appendix A.

2 Data sets

2.1 Planck CMB maps

The Planck satellite has surveyed the full sky in nine microwave and sub-millimetre bands, with the Low Frequency Instrument (LFI, Planck Collaboration 2016d) contributing three bands between 30 and 70 GHz, and the High Frequency Instrument (HFI, Planck Collaboration 2015a, b) featuring six channels between 100 and 857 GHz. Here we use the LFI 70 GHz map and the six HFI single-frequency maps222These maps were publicly released by the Planck Collaboration as part of their 2015 data release and are available in the Healpix ( format with at the Planck Legacy Archive:; the two other LFI maps (30 GHz, 44 GHz) do not add additional information because of their significantly poorer angular resolution (30 GHz) and higher noise level (44 GHz); see e.g. Planck Collaboration (2016d, 2015a). For a detailed description of the time-ordered data filtering, mapmaking and calibration process we refer to the HFI and LFI papers, and just summarize the properties of the maps relevant for our analysis in Table 1. For completeness, Fig. 1 shows a plot of the maps used in this analysis, and Fig. 2 shows the respective instrumental bandpasses.

The Planck maps up to 353 GHz used the orbital dipole of the CMB for calibration, whereas the 545 and 857 GHz maps were calibrated on planets (Planck Collaboration, 2016c). Here we assume the following conservative and simplified model for the calibration uncertainties: for the bands up to 353 GHz, we assume 1% absolute calibration uncertainty correlated between the individual bands. Similarly, we assume a 6% uncertainty for the 545 and 857 GHz bands, as in Mak et al. (2016). These calibration uncertainties are conservative estimates, but adopting tighter errors would have negligible effects on the results presented in this paper.

The broad spectral range makes the Planck maps well suited to investigating a possible Sunyaev-Zel’dovich signature at the lower frequencies, with the higher frequencies providing constraints on dust emission.

However, for all Planck maps the beam is significantly larger than the angular extent of the typical QSO host galaxy (for reference, a length of proper size 0.1 Mpc at subtends an angle of ); this is especially true for the Planck low-frequency channels that are important for the SZ signature. Nonetheless, although we are not able to resolve individual sources, the Planck maps still allow us to constrain the average emission properties of the QSO sample.

2.2 AKARI far-infrared data

At the median QSO redshifts of and for typical dust parameters, we expect the dust SED to peak at around  THz. Therefore, additional data at frequencies  THz adds valuable information with which to fix the dust emission parameters. In turn, this allows a more reliable extrapolation of the dust SED down to low frequencies at which we might hope to see an SZ signal.

To complement the Planck maps at higher frequencies, we use all-sky FIR data from the AKARI satellite. For additional tests we have also used data from the Infrared Astronomical Satellite (IRAS, Neugebauer et al. 1984), as reprocessed by Miville-Deschênes & Lagache (2005) in the form of the IRIS maps. However, we do not use the IRIS maps in our main analysis for reasons discussed in Appendix A. We summarize the key properties of the FIR data in Table 1 and show the maps and instrumental bandpasses in Figs. 1 and 2, respectively.

The AKARI satellite has surveyed almost the entire far-infrared sky in four bands with central wavelengths between 160 and 65 .333The full-sky AKARI maps at Healpix  resolution are available from the Centre d’Analyse de Données Etendues: In comparison to IRAS/IRIS, the nominal noise levels are broadly similar, but AKARI has a higher angular resolution of compared to the resolution of the 100 IRIS map. The relatively high resolution of AKARI is particularly useful to disentangle ‘actual’ point sources from objects that are merely unresolved in the high-frequency Planck and the IRIS 100 maps. The data processing, mapmaking and calibration of the AKARI maps is described in detail by Doi et al. (2015) and Takita et al. (2015).

The AKARI map at 90 has the lowest noise level and is the only AKARI map for which a reliable calibration is available down to low intensities (see Table 1 and Takita et al. 2015). The other AKARI bands have large and uncertain calibration errors at low intensities, making their robust inclusion into the analysis problematic. If appropriate calibration uncertainties were folded into the analysis, these bands would carry little statistical weight. If, however, they were included with underestimated calibration errors, they could indeed bias our results (as we demonstrate in Appendix A.3). We therefore include only the 90 map in our main analysis.

2.3 SDSS QSO catalogue

In this paper we make use of catalogues of spectroscopically confirmed QSOs created from SDSS-II (York et al., 2000) and SDSS-III (Eisenstein et al., 2011) Baryon Oscillation Spectroscopic Survey (BOSS, Dawson et al. 2013) data. In particular, we use the QSO catalogues from the seventh (DR7, Schneider et al. 2010 and twelfth (DR12, Pâris et al. 2016) SDSS data releases.444The catalogues are publicly available at and, respectively. The QSO target selection process is described in detail by Ross et al. (2012); here we simply summarize the points relevant to this work. Whereas the DR7 sample mostly contains ‘low’ redshift () objects, the DR12 sample specifically targeted QSOs with (Pâris et al., 2016). However, as a result of a colour degeneracy in the target selection from photometric data, lower-redshift QSOs were also observed, leading to a secondary maximum in the redshift distribution around . From these two catalogues, we create a merged sample. It is worth noting that there is a non-zero overlap between the two catalogues, mostly because previously confirmed QSOs were re-observed for DR12 (Ross et al., 2012). When combining the two samples, we remove any duplicates. The sky coverage of the merged sample is shown in Fig. 3.

Figure 3: Sky coverage of the merged SDSS QSO sample in Galactic coordinates: for plotting purposes we have smoothed the QSO density on a scale of .

Both catalogues provide several redshift estimates. These include the SDSS pipeline -estimate, as well as results obtained via principal component analysis, the position of specific emission lines, and visual inspection. Any of these redshift estimates is accurate enough for our purposes. For the main analysis we use the SDSS pipeline redshifts, but have verified that our results do not change if we use any of the other methods. Fig. 4 shows the redshift distribution of the DR7 and DR12 QSO samples and the combined sample. Since we are interested in the average emission properties of the QSOs, we remove the sparsely populated low- and high- tails of the redshift distribution ( and ). These cuts remove less than a percent of the sample. For our main analysis we use the full redshift range of the merged sample, containing 377,136 spectroscopically confirmed QSOs. To test for redshift evolution of the measured SZ and dust parameters, we also split the sample into five redshift bins with approximately equal number of QSOs (see Section 4.4).

Figure 4: Redshift distribution of the QSO samples: the red and green solid lines refer to the individual DR7 and DR12 samples, whereas the blue histogram shows the merged sample used in the main analysis (duplicates have been removed). The dashed vertical lines at the edges denote the redshift cuts used in the main analysis.

2.4 RedMaPPer galaxy clusters

For validation and testing purposes, we also use a catalogue of galaxy clusters identified in the SDSS DR8 data with the red-sequence Matched-filter Probabilistic Percolation (redMaPPer) cluster finder (Rykoff et al., 2014). Most of the optically identified clusters do not have an SZ signal that is strong enough for a blind detection in the Planck data. By stacking them, however, we should obtain a significant SZ detection, which will serve as a test for our measurement pipeline. For similar stacking analyses of clusters and CMB data see e.g. Afshordi et al. (2007); Diego & Partridge (2010); Planck Collaboration (2011); Hajian et al. (2013).

Briefly, the redMaPPer cluster finding algorithm uses a spectroscopic training set to calibrate red-sequence models at various redshifts, which are then used to iteratively assign membership probabilities to galaxies in the vicinity of a given cluster candidate. The resulting number of member galaxies — the optical richness — has been found to provide a low-scatter mass proxy (e.g. Rykoff et al. 2012; Saro et al. 2015); in addition the redMaPPer algorithm provides accurate photometric redshift estimates. The original SDSS DR8 catalogue has been compared to CMB data from Planck by Rozo et al. (2015). Subsequently, the algorithm has been further developed and applied to Dark Energy Survey data (Rykoff et al., 2016); Saro et al. (2015), Soergel et al. (2016), and Saro et al. (2016) have used these catalogues for SZ studies in conjunction with CMB data from the South Pole Telescope (SPT) SZ survey.

In this work, we use an updated version (v6.3) of the SDSS catalogue that includes the improvements to the algorithm described by Rykoff et al. (2016).555The catalogue is publicly available at Prior to the cuts we apply during the analysis, the catalogue contains 26,111 clusters with in the redshift range is , corresponding to a lower mass threshold of (Saro et al., 2015). The median richness and redshift of the catalogue are and , respectively.

3 Analysis methods

3.1 Selection of a clean sample

A crucial part of this analysis is the selection of a sample that is not affected by the various possible contaminants such as point sources, radio emission at low frequencies, and Galactic dust at higher frequencies. Here we describe and motivate these cuts and discuss their impact on the sample size.

Radio-loud cut for QSOs:

Some SDSS QSOs are radio-loud (see e.g. Kellermann et al. 1989), that is, they show synchrotron emission in the radio bands. Typically, the synchrotron contribution to the microwave SED can be modelled as a falling power law in frequency, thus contributing mainly at the lowest frequencies used in the analysis. Radio emission has the potential to hide an SZ signature, as its positive contribution to the low-frequency SED can partially cancel a negative SZ signal. We therefore remove from the sample all QSOs that have a radio counterpart detected in the Faint Images of the Radio Sky at Twenty-Centimeters (FIRST) survey (White et al., 1997) with a nominal source detection threshold of 1 mJy at 1.4 GHz; for details about the matching process see Pâris et al. (2016). The FIRST sky coverage mostly overlaps with the SDSS/BOSS footprint; however small regions at the edges were not covered. To avoid contamination from radio sources in these regions, we additionally remove all QSOs in parts of the sky that have not been observed by FIRST. These two cuts reduce our sample size by approximately 14%.

Even after the application of this cut, it is still possible for radio sources below the FIRST detection threshold to contribute to our flux density measurement at the lowest frequencies. We derive an estimate for the contamination by undetected radio sources in Section 4.3, and explicitly test for it by fitting a radio component to our measured SED.

Radio contamination for clusters:

As cluster galaxies can also host radio sources (see e.g. Gupta et al. 2016), we have tested removing clusters from the redMaPPer sample if they are close to a FIRST-detected radio source However, because of the spatial extent of clusters and the large Planck beams, this leads to a large number of chance associations. For this reason, we do not perform a ‘radio-loud’ cut for the clusters. It is therefore possible that residual synchrotron contamination might affect the lowest frequency, as discussed in Section 4.1 below. Because of the much higher SZ signal, synchrotron contamination is a much smaller problem for the clusters than for the QSOs.


We next remove objects in regions with high contamination by Galactic dust. In our main analysis we only include objects within the 40% of the sky with the lowest Galactic dust emission, as quantified by the masks produced by the Planck Collaboration. AS SDSS/BOSS mostly covers the Northern and Southern Galactic caps, this conservative cut only reduces our sample size by approximately . To ascertain that our results do not depend significantly on the choice of mask, we repeat our analysis with an even more conservative mask (using 30% of the sky), and a slightly less conservative mask (using 50% of the sky). In both of these tests, we obtain results that are consistent with the main analysis, demonstrating that contamination by Galactic dust is not a significant issue.

Point source cut:

Resolved point sources in the Planck sub-mm bands typically have flux densities of Jy, with the brightest sources reaching 10 Jy or more (e.g. Planck Collaboration 2014b); this is two to three orders of magnitude larger than the typical sub-mm flux densities of a QSO host (, see Section 4.2 below). Even a relatively small number of point sources associated with QSOs in our sample can therefore affect our results.

To identify point sources, we use the Planck 2015 Catalogue of Compact Sources (Planck Collaboration, 2015c), the AKARI Far-Infrared Surveyor (FIS) Bright Source Catalogue (Ishihara et al., 2010), and for additional tests also the IRAS Point Source Catalogue (Helou & Walker, 1988).

We also include unconfirmed point source candidates and sources with unreliable flux measurements from the point source catalogues. We then remove all QSOs that are closer than to a point source in the respective catalogue in any of the bands used for our analysis; this conservative selection removes approximately 20% of the remaining QSOs. It is worth noting that the high resolution of the AKARI 90  band enables us to perform a more stringent point source removal than possible with the Planck bands alone.

3.2 Map filtering

Having defined a QSO sample, we turn to estimating the flux densities in the CMB and IR maps at the QSO positions. The first step is to construct a matched spatial filter (e.g. Haehnelt & Tegmark 1996; Schäfer et al. 2006) that is designed to recover the QSO signal in the presence of much larger contaminants (primary CMB, dust emission, instrumental noise). For an object centred at , we write the data as


where is the unknown amplitude of a known (beam-convolved) profile with rotational symmetry, , and are all other signals (which we consider as effective ‘noise’). We perform the filtering in the spherical harmonic domain, which avoids having to project submaps into an approximate flat-sky geometry.

Using the symmetry of the input profile and the convolution theorem on the sphere, the harmonic coefficients of the filtered map, , are related to the unfiltered ones according to (e.g. Schäfer et al. 2006)




Here and are the harmonic coefficients of the filter function and input profile, respectively. All and vanish due to the symmetry; is the ‘noise’ power spectrum.

As the QSO contribution to the full power spectrum is completely negligible compared to the primary CMB (at low frequencies), foregrounds (mostly dust emission at higher frequencies), and instrumental noise, we set and measure directly as the power spectrum of the unfiltered input maps. We use the same mask (apodized with a Gaussian with 2 deg FWHM) as for the QSO sample selection above, and correct for the effect of the mask and the pixel window function on the measured .

We further include an additional smooth high-pass filter, precluding any large-scale features () from biasing our flux density measurements. This also removes any monopole from the maps, removing any additive overall calibration uncertainty. The filtering is furthermore limited to , where for Planck and IRIS maps, and for AKARI. We also ensure that the filter smoothly rolls off to zero at to avoid ringing artefacts in the filtered maps.

Filtering for QSOs:

Except for the lowest- objects (that we have already removed from our sample), the SZ and IR emission by QSO host haloes is completely unresolved in the Planck maps. Verdier et al. (2016) constructed a filter assuming a Navarro-Frenk-White (NFW, Navarro et al. 1996) profile matched to a halo with at , resulting in . This is not necessarily an optimal choice, as indeed most studies point to somewhat lower QSO host masses (see Section 5.1 below); nonetheless it is significantly smaller than the Planck beam at any frequency (). Even at the higher spatial resolution of AKARI (), the QSO host haloes still remain unresolved. Therefore we assume in this work that the input profile for QSOs is simply given by the beam (normalized to unit amplitude in pixel space).

In this case, the matched filter reduces to a simple Wiener filter on the sphere (e.g. Tegmark & de Oliveira-Costa 1998; Chatterjee et al. 2010) with


where , is the instrumental beam in harmonic space and (the extra normalization factor occurs because we have normalized the profile to unit amplitude in pixel space).

This construction of the filter is model-independent since we do not need to assume a specific profile. As an example for the resulting filters, we show in Fig. 5 the filtering kernels for the first five Planck bands from 70–353 GHz; these are the bands with sensitivity to the SZ signal.

Figure 5: Filtering: we show here the filter kernels for the unresolved QSOs in the first five Planck bands, which are those sensitive to a potential SZ signal. From low to high frequencies, the resolution of the Planck bands improves; therefore the peak of the filtering kernels are shifted to higher . For plotting purposes only, we have smoothed the kernels with a running average of .

Filtering for clusters:

The redMaPPer clusters that we use as a test of our analysis pipeline are partially resolved in the higher-frequency Planck maps and mostly resolved in the AKARI maps. To account for this, we modify the way we construct the filter as follows: we use a projected isothermal -profile (Cavaliere & Fusco-Femiano, 1976) with the index (not to be confused with the dust spectral index ) set to unity, i.e.


as an input profile for the filter. This is a widely used choice, both for blind SZ cluster detection (e.g. Bleem et al. 2015b), and matched-filter estimates for the CMB temperature (or intensity) at the positions of optically selected clusters (e.g. Soergel et al. 2016).

The profile is scaled with the projected core radius ; here we use . For a cluster sample with similar richness limits, but with slightly higher median redshift , Soergel et al. (2016) found to be a good match to the actual cluster profile, so that should be reasonably well matched in our case. As we mostly intend the clusters to be a test case for our SZ extraction algorithm, we do not attempt to further optimize this choice or to construct a more sophisticated (e.g. adaptive) matched filter. Nonetheless we have checked that our recovered emission parameters do not depend strongly on this choice by repeating the filtering with and .

Using this particular profile has the advantage that there is an analytic expression for , the spherical harmonic coefficients of the input profile. This can be seen as follows: (1) In the flat sky case, the Hankel transform of zeroth order of the profile is given by , where is the modified Bessel function of the second kind and is the wavenumber. (2) In the case of a function with circular symmetry (such as our input profile), the Hankel transform of order zero reduces to the 2D Fourier transform (modulo a normalization constant of ). (3) Generalizing this to the sphere (see equation 2), we compute via


(4) We finally account for the effect of the instrumental beam with .

This analytic construction of the filter has the advantage that we do not need to compute spherical harmonic transforms of the input profile numerically, making the filter construction both computationally efficient and robust against numerical errors.

3.3 Flux density measurement

The flux density contributed by a QSO or cluster is then


where is the pixel value in the filtered map at its position and is the filter profile for the band . In the case of the unresolved QSOs , so that the integral in equation 7 simply yields the beam area. For the clusters, we have normalized the -profile to unit amplitude in real space before convolving with the instrumental beam for the respective frequency. Therefore the matched filter returns the amplitude of the unconvolved profile, which is thus also used in the integral in equation 7. In this case, we have to evaluate the integral numerically, truncating the integration at a suitably high value . In this analysis, we choose , but we find that our results are insensitive to the precise choice of .

From the individual we then compute the mean flux densities at every frequency and their frequency-frequency covariance matrix as


where is the number of QSOs or clusters in the sample. We have also estimated from jack-knife or bootstrap resampling from the respective catalogue. The results from both of these approaches are in good agreement with the covariance estimated via equation 8.

It is worth noting that we do not use an inverse-variance weighting for estimating the mean flux density in equation 8; this is because the noise properties of our all-sky maps are relatively uniform once the Galactic plane is masked. The only notable exceptions are the ecliptic poles, where the Planck scanning strategy produces caustic-like structure in the number of observations per pixel. These ‘caustic’ regions are prone to systematics such as destriping errors; so upweighting them via strict inverse-variance weighting could potentially cause a bias.

We finally add the contributions to the covariance matrix from the absolute and relative calibration of the individual bands, which we have modelled as described in Section 2. For visualisation purposes, we define the correlation matrix and show it in Fig. 6 for the merged QSO sample. The off-diagonal correlations partially originate from the inclusion of correlated calibration uncertainties.

Figure 6: Correlation matrix for the QSO flux density measurement.

3.4 Emission modelling

We now describe how we model the flux density data with a multi-component emission model.

SZ effect:

In the non-relativistic approximation, the thermal SZ contribution to the total flux density can be written as (e.g. Birkinshaw 1999; Carlstrom et al. 2002)


where ; the SZ has a distinctive frequency dependence of the form


where is a dimensionless frequency. The SZ template itself, , has units of specific intensity; its amplitude, , has units of solid angle, and it corresponds to the average integrated Compton- parameter


In the above equation, and are the electron number density and temperature, and is the Thomson cross section.

Our QSO sample spans a large range in redshift, therefore the projected angular sizes vary in our sample. Even in the absence of any intrinsic evolution of the QSO host properties, this could lead to a redshift evolution of the integrated Compton- parameter of individual objects, . Variations in the intrinsic properties of the host haloes can be disentangled from size variations by binning the QSO data in narrow redshift bins (see Section 4.4 below). This can, of course, increase the statistical noise. To co-add the data over a wide redshift range, we therefore define the rescaled (or intrinsic) parameter (e.g. Planck Collaboration 2013)


where is the angular diameter distance to redshift and is the dimensionless Hubble expansion rate. Unless otherwise noted, we use a flat CDM cosmology with and to compute these quantities. In the case of a purely self-similar scaling (Kaiser, 1986) and fixed halo mass, is independent of redshift. For the actual measurement we absorb this rescaling in the template, such that the amplitude of a recovered SZ signal directly measures .

It is worth noting that measures the SZ signal integrated over a cylinder along the line of sight. Under the assumption of a specific gas pressure profile, it is possible to convert this into a spherically integrated quantity, which is however not directly observable. As gas pressure profiles for the low-mass haloes that typically host QSOs are not well known, we have kept our analysis as model-independent as possible, and therefore we do not perform such a conversion; this should be kept in mind when comparing our results to those of previous studies. We note further that for the comparison to simulations presented in Section 5.1 such a conversion is not necessary either, as it is possible to directly predict the cylindrically integrated signal from the simulations.

Once we have measured , we can relate it to the thermal energy content of the hot halo gas as follows: for an individual halo at redshift , the corresponding thermal energy is given by


where is the average particle weight per electron and is estimated from the mean redshift-independent via equation 12. The mean thermal energy of our sample, accounting for the redshift distribution, is then given by


where is the normalized redshift distribution of the sample.

Dust emission (CIB):

The second contribution to our model is the sub-mm emission from the QSO host galaxy or other galaxies along the same line of sight. This is mostly optical and UV light that has been absorbed by dust grains and re-emitted in the infrared or sub-mm. Essentially all of this emission is unresolved, hence it is also known as the cosmic infrared background (CIB).

We model the CIB with a single-component modified blackbody (e.g. Blain et al. 2002),


where is the Planck spectrum, is the dust spectral index, and is the dust temperature. We normalize the dust emission at a frequency GHz. Its amplitude, , has units of solid angle, analogously to . As we are mainly interested in a potential SZ signature, we treat as a purely phenomenological parameter quantifying the dust emission, and do not relate it to quantities such as the dust mass by adopting a specific model. We note that the parameters governing the dust model can be significantly degenerate in SED fits: the parameters and follow a well-known degeneracy, such that typically only the combination is well constrained by the data (e.g. Verdier et al., 2016).

The standard modified blackbody template of equation 15 is a simple parametrization of the complex dust physics. With the limited number of frequencies used in this paper, it is difficult to test more complex models, such as a two-temperature dust component (see e.g. Appendix C of Crichton et al. 2016).

Dust emission (Galactic):

Galactic dust emission is also well approximated by a modified blackbody (Planck Collaboration, 2014a). This makes it hard to disentangle Galactic dust from the CIB based on their spectral properties only (e.g. Planck Collaboration 2015d, 2016a; Mak et al. 2016). However, emission from Galactic dust is uncorrelated with the QSO sample, adding noise but not bias to our measurement of the CIB amplitude. Furthermore, we focus on the northern and southern Galactic caps, where we do not expect a strong contamination by Galactic dust. We therefore make no attempt to model Galactic dust, but instead test the sensitivity of our results by repeating our analysis with different masks (see Section 3.1).

Synchrotron emission:

We have removed all QSOs with a known radio counterpart, and those outside the footprint of the FIRST radio observations, from our sample. In Section 4.3 we estimate that the contribution from radio sources below the FIRST detection threshold is mJy at 100 GHz; this is significantly smaller than the uncertainties on the lowest-frequency points. We therefore do not expect a significant level of contamination by synchrotron emission. Nevertheless, we test for residual synchrotron contamination by fitting an additional component with


where the spectral index is typically found to be between 0.5 and 1 (e.g. Planck Collaboration 2016b). Here we either fix the spectral index to the most conservative choice , or fit it jointly with the radio amplitude .

Other contributions:

There are several other sources of anisotropies in the microwave and sub-mm sky, in particular the primary CMB at the lower frequencies, the Poisson component of the CIB at the higher frequencies, and instrumental noise at all frequencies. All of these are uncorrelated with the positions of the QSOs. Primary CMB and instrumental noise can be positive or negative, therefore they cancel out when averaging over all QSOs, adding noise but no bias to our flux density measurements. On the other hand, the CIB Poisson component is always positive and could in principle bias our measured SED. In practice, we find that after filtering the maps as described in Section 3.2 above, this contribution is negligible. We explicitly demonstrate this by replacing the QSO positions with randomly drawn points with the same distribution on the sky as the original objects. These are then processed with the same pipeline; we show the results in Appendix A.2. We note that this test also demonstrates that our results for the Compton- parameter are not affected by a stacking bias.

Full model:

The model for the emission of QSOs at redshift and frequency is thus


where the sum over labels components. For the main analysis , while for additional tests we also consider radio emission. Furthermore, contains the parameters for component and is the full parameter vector that we are aiming to constrain; i.e.  for the main analysis. For the radio contamination test we also include and .

In this analysis, we do not measure the emission at a fixed central frequency , but rather over a range of frequencies defined by the instrumental bandpasses. In our case, the latter are relatively wide, especially for the AKARI bands, some of which span almost a factor of two in frequency (see Fig. 2). Additionally, the modified blackbody dust spectrum rises steeply for GHz, and then falls again sharply for THz. For these reasons, we must include the instrumental bandpasses in our model prediction.

As the models for the individual components depend on , one in principle has to evaluate them for every single QSO. To speed up the evaluation of the model, we group the QSOs into redshift bins, evaluate the model at the redshifts of the bin centres, and weight the bins with the redshift distribution . We have tested that with the binned model only differs from the full model at the sub-percent level, which is well below the statistical uncertainties.

Our full model is thus


3.5 Model fitting

Now we have all the ingredients in place to estimate the parameters of the QSO emission model. We write the posterior probability distribution of the parameters in Bayesian fashion as


where is the vector of mean observed flux densities, is the likelihood and is the prior assumed on the parameters. In our case, we assume Gaussian likelihoods, so that


We restrict all parameters to positive values to avoid unphysical results; apart from this restriction, we choose the priors to be flat and uninformative. We then sample from the posterior using the emcee implementation (Foreman-Mackey et al., 2013) of an affine-invariant Markov Chain Monte Carlo (MCMC) ensemble sampler (Goodman & Weare, 2010). To assess the convergence of our chains, we use the Gelman-Rubin criterion (Gelman & Rubin, 1992; Brooks & Gelman, 1998).

We then report the results for the QSO emission parameters and their uncertainties as the 50th, 84th and 16th percentile of the one-dimensional (i.e. marginalized over all other parameters) posterior distribution of the respective parameter. If a parameter is consistent with zero at the level, we also report the upper limit, defined as the 95th percentile of the respective 1D posterior.

4 Results

4.1 Galaxy clusters

Figure 7: Results for the RedMaPPer galaxy clusters. Top: we show here the measured flux densities and the best-fitting model together with its individual components. For clarity, we have included both the (solid) and the (dashed) error bars on the measured flux densities. Furthermore, the red shaded region shows the uncertainty on the amplitude of the SZ signal. Note that the vertical scale is linear for flux densities below , and then switches to a logarithmic scale; the dashed horizontal line marks the transition. Bottom: results of the MCMC parameter estimation. The panels on the diagonal show the marginalized 1D posteriors on the four parameters, whereas the off-diagonal panels show marginalized 2D constraints. Here, the colour scale denotes the density of MCMC samples, while the contours show the and confidence regions. Note that and are displayed in units of and , respectively, to allow for a direct comparison with the QSO results in Fig. 8 below.

We begin with discussing the results from the RedMaPPer clusters, which we provide as a test of our analysis pipeline, demonstrating that we are able to recover a significant SZ signal associated with optically selected clusters, together with physically sensible dust emission parameters. We show in Fig. 7 the flux densities measured from the cluster sample together with our best-fitting model from the MCMC analysis, and the constraints on the emission parameters.

For the integrated, rescaled Compton- parameter of the clusters we obtain


which is a highly significant detection (), as expected. From equation 14, this corresponds to a mean thermal energy in the hot halo gas of


Furthermore, we measure an amplitude of the CIB dust emission and modified blackbody parameters of


These results are also summarized in Table 2. The strong correlated emission in the sub-mm bands indicates that the dust emission from the galaxies within redMaPPer clusters is clearly detected in our data. The levels of sub-mm emission that we find are compatible with those found by Bleem (2013) when stacking maps from the SPT-SZ survey and Herschel-SPIRE at the positions of optically selected clusters from the Blanco Cosmology Survey (Bleem et al., 2015a). Evidence for dust in clusters has also been reported by Chelouche et al. (2007) and Gutiérrez & López-Corredoira (2014).

Our dust emissivity index is in good agreement both with the theoretical expectation of in the FIR (e.g. Blain et al. 2002; Franceschini 2000) and the value used to fit sub-mm SEDs of individual galaxies (e.g. Calzetti et al. 2000). Furthermore, the dust temperature of is consistent with the temperature seen in normal quiescent galaxies (e.g. Smith et al., 2012), as we would expect for member galaxies of the relatively low-redshift SDSS RedMaPPer clusters.

Our simple two-component model nominally does not provide a good fit to all data points: we find when evaluated at the parameters given by the 50-th percentile of the MCMC chain.666Note that this is not the nominal minimum , but should be close to it in the case of Gaussian posterior distributions. However, this unusually high value is driven largely by the 70 GHz data point, which deviates from the best-fitting model by . This is mainly caused by radio emission from cluster members, which is known to be a potential bias for SZ studies (e.g. Bleem 2013; Gupta et al. 2016). The value of decreases significantly if we also fit for the amplitude of a radio component () or perform a radio cut that removes all clusters within of a FIRST-detected radio source (). However, the former ignores the uncertainty on the radio spectral index (see Section 3.4), whereas the latter removes roughly half of our clusters (mostly due to chance associations) and introduces an additional selection effect (see Section 3.1). In these cases, the constraint on (or ) is affected at the level. We therefore include an additional systematic uncertainty when comparing our results to simulations in Section 5.1. Another option would be to completely discard the 70 GHz data point for the cluster analysis. This would also lead to a significantly improved , while not changing our results significantly.

4.2 QSOs

We proceed with the constraints on the QSO emission model, which are the main result of this paper. Similarly to Fig. 7 for the clusters, Fig. 8 shows the flux density measurement, the best-fitting model, and the MCMC results for the merged QSO sample. Here our dust+SZ model provides an acceptable fit to the data ( when evaluated at the 50-th percentile as before).

We measure an SZ amplitude of


The data indicate a mild preference for a non-zero SZ signal, but its amplitude is consistent with zero at the level. This translates into a thermal energy of


These numbers are significantly smaller than those reported by Ruan et al. (2015), and more in line with the feedback energetics suggested by hydrodynamical simulations. Section 5 below presents a more detailed discussion of the implications of these results for QSO feedback energetics and a comparison to previous studies.

Regarding the dust parameters, we measure an amplitude of ; the modified blackbody parameters are


These QSO results are also summarized in Table 2. Here we find a dust emissivity index that is marginally higher than that for the clusters, but it is consistent with at the level. The best-fitting dust temperature determined from the QSOs is significantly higher than that for the clusters. This is perhaps not surprising, as the dust in the QSO host galaxy may be heated by the UV and optical emission from the central QSO. In addition, as the typical redshift of the QSOs is — close to the peak epoch of cosmic star formation (e.g. Madau & Dickinson 2014) — the typical QSO host galaxy should have a larger fraction of O/B-type stars than a quiescent low- cluster galaxy, providing additional UV/optical flux. Both of these effects will lead to an increase of the dust temperature compared to normal galaxies. We note, however, that the dust temperature of equation 26 is lower than the temperature of seen in optically luminous, higher-redshift QSOs. (e.g. Priddey & McMahon, 2001; Omont et al., 2003).

For this main result, we have only used the bands that have a well-determined calibration and satisfy a successful stacking null test (see Sections 2.2 and Appendix A.2, respectively). In Appendix A.3 we repeat the analysis adding in additional AKARI and IRAS bands. In this case, the AKARI bands with uncertain calibration pull the dust solution away from the best-fit values determined above. This leaves more room for an SZ signal, but at the same time provides a poor overall fit to the data. We provide this example in order to demonstrate the sensitivity of the SZ result to an accurate dust solution.

Figure 8: Results for the merged QSO sample. Top: we show here the flux density measurement and best-fitting emission model for the merged QSO sample. The vertical scale is linear for flux densities below , and logarithmic above this value. As before, the red shaded region shows the confidence region for the SZ contribution to the measured SED. Unlike the case for clusters, there is only weak evidence for an SZ contribution to the measured SED. Bottom: this panel shows the MCMC results for the QSO sample. It is worth noting that there is no strong degeneracy between and the dust parameters, showing that our result on is robust against uncertainties in the dust solution.

4.3 QSOs: constraints on residual radio emission

redMaPPer clusters 0.37
merged QSO sample 2.15
   + radio amplitude 2.15
Table 2: We summarize here our best-fitting model parameters from Section 4. Except for the 95% upper limit on , all error bars are , obtained from the 16th and 84th percentile of the marginalized posterior distributions. Note that the conversion between and depends on the redshift distribution of the sample.

We have already removed all QSOs associated with radio sources above the FIRST 1.4 GHz detection threshold of  mJy, which removes of the QSOs. In our main analysis, we have therefore assumed that synchrotron emission does not contribute significantly to our flux density measurements. In this section we check this assumption by estimating the contribution to our flux density measurements from faint radio sources with FIRST flux density .

Naively one could approach this by simply fitting an additional radio component to the flux density measurements. However, only the two to three lowest-frequency points are sensitive to radio emission; therefore the SZ and synchrotron amplitudes are highly degenerate, introducing the risk of overfitting. This was also noted by Verdier et al. (2016), who found that their three-component multi-frequency filter could not distinguish between SZ and a synchrotron component with ‘negative’ amplitude, and vice versa. We therefore first estimate an upper limit for residual radio contamination from the emission properties of QSOs with a radio counterpart detected in FIRST. We find that the expected level of radio contamination is well below the statistical uncertainties. For completeness, we then repeat our MCMC analysis including a radio component, with a prior for the radio amplitude based on our estimate of residual radio contamination.

We begin by estimating the radio luminosity function of the SDSS/BOSS QSOs that have a radio counterpart detected in FIRST, finding that it is well approximated by


where the flux density range up to  mJy contains more than 95% of all detected radio counterparts in our sample; for larger flux densities the luminosity function cuts off exponentially. We assume that the scaling of equation 27 approximately holds true for as well. Then the ratio of QSOs associated with radio sources below the FIRST detection threshold to those with detected counterparts, , and the mean flux density of the undetected sources, , are given by


where is the lower flux density limit that we consider for our estimate of radio contamination. For concreteness we use  mJy, leading to and  mJy. The lower limit was chosen to be around 10% of the uncertainty of our 100 GHz flux density measurement; even fainter radio sources would contribute negligibly to our measurement.

The radio spectrum of synchrotron sources detected in microwave bands can typically be approximated with a broken power law , which is relatively flat () for and steepens () for , with  GHz (e.g. Planck Collaboration 2016b). Using these values, we estimate the residual radio contamination at 100 GHz as


We consider this a conservative estimate, for the two following reasons: (1) We have not accounted for a flattening of at the faint end. (2) By construction our radio counterparts are associated with high-redshift sources, and therefore likely at higher than the sample studied in Planck Collaboration (2016b). If the rest-frame spectra of these objects are nevertheless broadly comparable, we have likely underestimated the ‘downscaling’ in frequency in equation 29. Both of these effects would cause us to overestimate . We thus estimate the radio contamination at 100 GHz to be  mJy, which is at most a third of the purely statistical uncertainties on our flux density measurement. We therefore conclude that residual synchrotron emission does not contribute significantly to our measurement and can safely be neglected.

As an additional test, we repeated our MCMC analysis with an additional radio component. Based on the estimates presented above, we chose a prior . The results of the three-component fits are given in Table 2. The constraint on is weak and dominated by the prior choice ( at 84% C.L.). We find no evidence for a radio component, in agreement with our estimate in equation 29. The addition of a radio component has almost no effect on the parameters of the dust component, but does, however, slightly degrade our constraints on the SZ component. We now measure , or an upper limit of (95% C.L.).

We have also tested allowing a completely free radio amplitude (i.e. without the physically motivated prior from above), and/or to fit the radio spectral index jointly with the amplitude. In none of these cases we have found evidence for radio contamination or a significant detection of an SZ signal, consistently with our physical argument given above.

4.4 QSOs: redshift evolution

Figure 9: Redshift split: we show here two representative examples of the measured SED and best-fitting model. Left: Second redshift bin (). Right: Fourth redshift bin (). The redshifting of the SED is clearly visible when comparing the two panels.
Figure 10: Redshift split: we show here the redshift evolution of the parameters of the emission model (plotted at their respective median ). From top to bottom, the panels refer to the total thermal energy, CIB amplitude, the well-constrained dust parameter combination , and further and separately.

In the analysis presented so far, we have rescaled the Compton- parameter so that it is independent of redshift in the case of purely self-similar scaling (see equation 12). Nonetheless, especially given the large redshift range of the QSO sample, there could still be a redshift evolution caused by deviations from self-similarity. Furthermore, any change in the AGN activity over cosmic time would also affect the measured thermal energies. Therefore we proceed to split the QSO sample into five redshift bins with approximately equal number of objects, leading to splits at . We then repeat the measurement of the stacked SED and of the MCMC sampling; as examples, we show the results of the second and fourth -bin in Fig. 9, and we present in Fig. 10 the full redshift evolution of the measured emission parameters.

Similarly to the main sample, we also do not obtain a strong detection of the SZ signal in any redshift bin; only in the highest -bin the significance is . We further observe some redshift evolution in the dust parameters, with a particularly low and high in the second redshift bin. However, the combination , which is best constrained by the data, stays relatively stable (except for the highest -bin). Qualitatively, a similar trend in the dust parameters was also found by Verdier et al. (2016). Quantitatively, however, they reported a much larger variation between the low- and high- dust solutions, which is in some tension with the physical expectation of in the low-frequency limit of optically thin dust emission (e.g. Blain et al. 2002).

The evolution of the dust parameters could indicate changes in the dust properties, which could potentially be linked to the redshift dependence of the star formation rate and/or AGN accretion rate. However, our data does not allow for any more conclusive statements on this hypothesis.

5 Interpretation and Discussion

In this section, we first compare our results for both clusters and QSOs to hydrodynamical simulations that include AGN feedback. We then discuss the implications of our measurements for QSO host halo masses and feedback energetics. Finally, we compare our results and conclusions with those of previous studies.

5.1 Comparison to simulations

Figure 11: Projected Compton- maps: we show here the maps obtained from the snapshot of a halo with . The left panel shows the case without AGN feedback, whereas the right-hand side includes AGN feedback. While we see evidence of gas being pushed towards the outskirts of the main halo, the overall thermal energy only increases by in this representative example. The white circle denotes the beam of the Planck 70 GHz band, which is broadly representative of the ‘aperture’ probed by our analysis.

We use high-resolution zoom-in simulation of individual haloes spanning a mass range from to at , covering the expected masses for QSO hosts and the lower end of the galaxy cluster mass scale. The simulations have been performed with the smoothed particle hydrodynamics (SPH) code P-GADGET3 (which is based on GADGET-2, Springel 2005) using a zoomed initial conditions technique (Tormen et al., 1997). 777The simulations we use have an input cosmology that is slightly different from our fiducial cosmology (, , ), but given the uncertainties in our estimates of thermal energy and mass this does not affect our comparison significantly.

These simulations have previously been used to demonstrate the influence of AGN feedback on X-ray group and cluster scaling relations (Puchwein et al., 2008) and on the different stellar components of clusters (Puchwein et al., 2010). AGN feedback in the simulations is modelled as described by Sijacki et al. (2007). All of the simulations were run twice, resulting in a control set of simulations that include cooling and star formation but no AGN feedback together with a second set in which AGN feedback was switched on. This makes it relatively easy to disentangle the effects of AGN feedback from purely gravitational heating.

Since these simulations were run, there have been incremental changes to the AGN feedback model (e.g. Sijacki et al. 2015) and significant developments in cosmological numerical hydrodynamics. To check that our results are not affected significantly by these changes, we compute the integrated Compton- parameter and total thermal energy for one massive cluster (the AREPO-IL run in Sembolini et al. 2016) simulated with the moving-mesh code AREPO (Springel, 2010), using the same AGN feedback model as that adopted in the Illustris simulation (Vogelsberger et al., 2014; Genel et al., 2014). The resulting thermal energy agrees with the scaling relation determined from the P-GADGET3 simulations shown in Fig. 12 below. Changes in the hydrodynamic scheme and AGN feedback model are therefore unlikely to affect the conclusions presented in this section, especially given the large uncertainties of the observational results presented in this paper.

Depending on the halo mass, the diameter of the high-resolution zoomed-in region ranges from 7 to 30 comoving Mpc: this is much larger than the QSO host halo, as shown in Fig. 11, and it therefore contains a large fraction of the surrounding correlated large-scale structure. Due to the use of zoom-ins, our SZ estimate from the simulations will miss any contribution from the line of sight outside this region, which would otherwise be present in a full cosmological lightcone. However, we show explicitly in Appendix A.2 that our measurement on the data is not affected by stacking bias from uncorrelated structure along the same line of sight. Given the present level of accuracy, a precise quantification of residual contribution from correlated structure outside the zoom region is beyond the scope of this work.

To compare with our observational results, we choose the simulation snapshot that is closest to the median redshift of the corresponding sample, leading to for the QSOs, and for the clusters. We then create pixel maps of the Compton- parameter from the zoom simulations by projecting the gas particles along the line of sight, which we choose to be parallel to one axis of the simulation box. Perpendicular to the line of sight, we project all particles within a square centred on the respective halo of interest. For the latter, we use a fixed physical size of .

Following Springel et al. (2001), we make use of the fact that for every gas (SPH) particle in the simulation, its contribution to the integrated Compton- parameter is proportional to its pressure times volume, given by . Here , , , and are the adiabatic index, primordial hydrogen fraction, mean molecular weight and electron-to-hydrogen number density ratio (the latter is traced dynamically in the simulations); and are the mass and internal energy per unit mass of the particle. We evaluate the integral in equation 11 by projecting all particles as described above, weighting the contribution of particle to the map with , where is the projected SPH smoothing kernel (e.g. Springel et al. 2001).

Fig. 11 shows the resulting -maps, both with and without AGN feedback, for a halo at with a mass of (defined as the mass that is enclosed within a sphere that has 200 times the critical density of the universe at the given redshift). In this representative example, AGN feedback affects the distribution and temperature of the gas around the central halo, leading to a more extended and more diffuse Compton- signature. However, the total thermal energy is only increased by , an effect too small to be seen in our data since we find only marginal evidence for an SZ signal in the QSO sample.

As we use different filters for the clusters and the QSOs, we need to process the -maps in different ways, as described below.

Figure 12: Thermal energy as estimated from the simulations: the blue points denote the results from the control-run simulations without AGN feedback; red are simulations of the same haloes with AGN feedback. Left: we show here the result at the median of the redMaPPer cluster catalogue (low redshift). The purple horizontal band shows our estimate of from the redMaPPer clusters, where we have conservatively assumed uncertainty to account for systematic effects such as the choice of the filter profile (see Section 3.2) and radio contamination (see Section 4.1). The grey vertical region shows an estimate of the typical of the clusters. Right: this panel displays the progenitors of the same objects, at the median redshift of the merged QSO sample. Again, the purple and grey shaded bands show our estimate of , and the range of QSO host masses from the literature, respectively.

Low redshift (clusters):

For the clusters, we filtered the maps using a -profile matched filter with ; to estimate the total , we have then integrated the profile out to . For the comparison to simulations, we therefore integrate the -maps obtained from the low-redshift snapshot at the median of the redMaPPer clusters within a circular aperture of the same radius. We then use equation 13 to convert the integrated Compton- parameter to the total thermal energy. The left-hand panel of Fig. 12 shows the dependence of on the halo mass.

We find a tight scaling relation between halo mass and thermal energy over a large range in halo masses, with a slope that is approximately consistent with the self-similar relation (Kaiser, 1986). At larger halo masses (), AGN feedback leads to a small increase in the overall thermal energy. At lower masses, however, AGN feedback leaves the total thermal energy almost unchanged, while reducing the halo mass. This arises because strong AGN feedback can remove some of the gas from low-mass haloes. At higher mass, the gas stays bound and receives extra heating from the feedback.

To compare our results from the data to the simulations, we convert richness to mass using the relation from Simet et al. (2016), who have calibrated the masses of SDSS redMaPPer clusters from stacked weak lensing measurements. Melchior et al. (2016) have extended this to using Dark Energy Survey (DES) Science Verification cluster and lensing data, finding good agreement with the Simet et al. (2016) relation.888Saro et al. (2015, 2016) have used clusters detected both in the SPT-SZ survey and DES to calibrate the masses from the SZ detection significance, also yielding a broadly consistent mass calibration. The redMaPPer clusters used for our measurement have a median richness of . Including a factor , this leads to a typical halo mass of . To account for statistical and systematic uncertainties in this estimate, we conservatively assume a uncertainty on this value.

We have measured a thermal energy of for the SDSS redMaPPer clusters, which is fully consistent with the scaling relation inferred from the simulations. This agreement demonstrates that our pipeline is able to robustly estimate the thermal energy of a sample of objects from CMB and sub-mm maps, and that the results are compatible with theoretical expectations.

High redshift (QSOs):

We now turn to the high-redshift snapshots at the median of the QSO sample. Here the QSO emission is completely unresolved by the Planck maps and so we simply integrate the entire simulated -map, before using equation 13 to compute the thermal energy. The results are shown in the right-hand panel of Fig. 12. As for the low-redshift clusters, we see a scaling relation between and halo mass, but with higher scatter. This increased scatter is caused by the different environments of the simulated haloes, which contribute to estimated within the Planck beam (see Fig. 11). We have explicitly verified this by integrating the -map only within the projected of the respective halo, finding a significantly lower scatter in the scaling relation.

AGN feedback slightly increases the thermal energy over the entire mass range by about , with no significant dependence on mass. Unlike for the case of low-mass haloes at low redshift where AGN feedback leads to an appreciable halo mass reduction due to the powerful “radio mode" feedback, here this is not the case as energetically less efficient “quasar mode" feedback is operating (for further details see Sijacki et al. 2007).

As we are not sampling the higher-luminosity part (i.e. ) of the QSO luminosity function very well because of the limited size of our sample, we have performed several tests to establish that our comparison between observed and simulated QSOs is adequate. For every black hole in the simulation accreting in the radiatively efficient way, we have estimated its bolometric luminosity as , where is the current accretion rate and we have set the radiative efficiency to . We have then computed the luminosity function and have checked that it is sufficiently well populated at QSO luminosities of . Repeating the same procedure at , we find that the luminosity function is well populated out to , consistent with the expectation of a higher accretion rate at this redshift. However, we find that the impact of AGN feedback on the total thermal energy is lower at than at and that the relative increase in does not strongly depend on the current accretion rate of the central black hole. This indicates that measured from the SZ effect probes the integrated past accretion history, thus making our comparison to observations meaningful.

Mass estimates for AGN host haloes are usually derived from clustering measurements or halo occupation distribution (HOD) modelling and vary considerably, with values ranging from (White et al., 2012), to (Croom et al., 2005; Richardson et al., 2012), and up to (Porciani et al., 2004). Most studies find no evidence for significant evolution of the host mass with redshift (note however the tentative indication for an upturn in the typical mass at reported by Richardson et al. 2012). Galaxy formation simulations suggest that masses of host haloes inferred from clustering measurements may be overestimated by a factor of at redshifts (DeGraf & Sijacki, 2016). The masses of the host haloes for the QSOs in our sample are, therefore, extremely uncertain. The grey band in Fig. 12 shows the wide range of halo masses quoted in the literature.

The results presented in Section 4.2 suggest thermal energies of , but the uncertainties are so large that we cannot claim a statistically significant detection. Comparing our results to the simulations indicates halo masses of , consistent with the clustering or HOD measurements of QSO host masses discussed above. Our upper limit on restricts the allowed mass range for QSO host haloes to . Given these large uncertainties in both the observational constraint on the SZ amplitude and the QSO host halo masses, it will be difficult to detect the small changes in thermal energy associated with AGN feedback via an SZ signature in the foreseeable future (unless our simulations have grossly underestimated the effects of AGN feedback, which seems unlikely).

5.2 Comparison with previous work

Ruan et al. (2015) found a strong correlation between the Hill & Spergel (2014) Compton- map and the SDSS DR7 QSO sample, inferring a total thermal energy in the hot halo gas of . Such a high value is inconsistent with the scaling relations determined from our simulations (see Fig. 12) for any plausible value of the QSO host halo masses. Our upper limit on the thermal energy associated with the QSO host haloes is approximately an order of magnitude lower than the Ruan et al. (2015) estimate. Subsequent studies by Verdier et al. (2016) and Crichton et al. (2016) showed that the emission around QSOs at frequencies is dominated by dust; here we confirm this finding. The signal detected by Ruan et al. (2015) is therefore most likely related to residual dust contamination in the Hill & Spergel (2014) Compton- map. This dust emission must be subtracted accurately in order to recover a constraint on the SZ amplitude. Our results differ from those of Verdier et al. (2016) and Crichton et al. (2016) in the following important aspects.

For low-redshift QSOs (), Verdier et al. (2016) found a indication for a negative Compton- parameter, whereas we obtain a preference for a positive SZ amplitude at all redshifts. At higher redshifts (), they found a detection of an SZ signal, with an amplitude corresponding to a spherically-integrated Compton- parameter of . This can be converted to the cylindrically-integrated total Compton- parameter we measure by assuming a specific pressure profile for the halo gas. For the Arnaud et al. (2010) profile adopted by Verdier et al. (2016), the conversion factor is , and so their result translates to . This value is almost a factor of two larger than our measurement of for QSOs at .

The negative SZ amplitude at and the relatively high SZ at determined by Verdier et al. (2016) correlate with large changes in the parameters of the dust model: low values of and high dust temperatures at and high values of and low values of at . Although we find qualitatively similar trends for the dust parameters (see Section 4.4 and Fig. 10), the variations of and with redshift found in our analysis are less extreme than those found by Verdier et al. (2016). The differences between our and their SZ results may therefore arise from instabilities in their dust solutions caused by the lack of far-infrared data in their analysis: only the introduction of the AKARI FIR data allows us to firmly anchor the dust SED.

On the other hand, Crichton et al. (2016) reported a evidence for an SZ signal (averaged over their entire sample spanning the redshift range ), corresponding to a measured thermal energy of . This is marginally higher than the value we find, but consistent within the measurement uncertainties. However, their estimate of the SZ amplitude (as quantified by the feedback efficiency), is again strongly degenerate with the dust parameters, , , which differ from those found here: , . A shift in the Crichton et al. (2016) dust parameters towards our values would reduce the amplitude of their recovered SZ signal (see their Fig. 3). The interpretation of their results as a detection of SZ signal therefore depends on the fidelity of their dust solution. Crichton et al. (2016) then used the - scaling relation of Planck Collaboration (2013) to estimate the expected thermal energy from purely gravitational heating. Assuming QSO halo masses in the range of , they concluded that gravitational heating accounts for a small fraction () of their measured signal. However, this conclusion is incompatible with the results from our simulations shown in Fig. 12, where we find that the AGN contribution to is always a small fraction of that from gravitational heating.

Finally, Spacek et al. (2016b, a) investigated constraints on AGN feedback energetics using samples of massive quiescent elliptical galaxies with redshifts in the range , with the aim of finding signatures of past AGN activity. By stacking CMB maps from SPT and ACT, respectively, they found low-significance indications of an SZ signal associated with these galaxies, with corresponding thermal energies of , comparable to the values found for QSOs by Crichton et al. (2016) and in this paper. However, the SZ signal is consistent (with large uncertainties) with their estimates for gravitational heating and therefore does not require additional energy input from AGN feedback.

6 Conclusions

The analysis discussed in this paper was motivated by the results of Ruan et al. (2015) who presented evidence for a high-amplitude SZ signal in the vicinity of QSOs, suggestive of strong AGN feedback. Subsequent work by Verdier et al. (2016) and Crichton et al. (2016) showed that dust emission dominates at microwave and sub-mm frequencies and must be subtracted to high accuracy to recover an SZ signal. Both of these groups found evidence for an SZ signal, though at a much lower amplitude than found by Ruan et al. (2015). However, the corresponding thermal energies and dust emission parameters that they recovered differ substantially.

Our main contribution has been to analyse the cross-correlation of QSO catalogues from SDSS/BOSS with maps from Planck and AKARI. In particular, the inclusion of the AKARI far-infrared data at extends beyond the peak of the dust emission and helps to break the degeneracy between dust parameters and the amplitude of any SZ signal. We have paid careful attention to the selection of a clean QSO sample, removing objects associated with radio sources, extragalactic point sources, and in areas of high contamination by Galactic dust. We then applied a filter to the CMB and sub-mm maps that optimally recovers the flux density of unresolved objects such as the QSO hosts. The emission around the QSO was modelled using two components, consisting of thermal SZ and dust emission (CIB) approximated by a single-temperature modified blackbody spectrum. We also experimented with adding a synchrotron component at low frequencies, which has little effect on our solutions, consistent with our expectations based on source counts.

We find indications for an SZ signal at low significance (). In particular, we do not reproduce the strong detection of an SZ signal for QSOs at found by Verdier et al. (2016). We also do not find the strong trends in the dust parameters with redshift reported by Verdier et al. (2016). The redshift dependence of the dust parameters in our analysis is more gentle, though we do find evidence for a rise in the dust temperature to and a lowering of the spectral index to at . Averaging over the entire redshift range, the best-fitting modified blackbody parameters are and , similar to the dust parameters of normal nearby galaxies. Our dust parameters disagree strongly with those determined by Crichton et al. (2016). Since the dust parameters are highly correlated with the SZ amplitude, our results suggest that the detections of an SZ signal reported by Verdier et al. (2016) and Crichton et al. (2016) should be treated with caution.

We have compared our results with hydrodynamic simulations of haloes run with and without AGN feedback. In these simulations, the effects of AGN feedback lead to small () enhancements to the SZ signal both at the low redshifts of redMaPPer clusters () and at the median redshift of the QSO sample (). For QSO hosts, our upper limits to the SZ signal are consistent with the simulations, provided the typical QSO host halo masses are , in agreement with the halo masses inferred by other techniques.

As an aside, we also analysed the SZ signal of redMaPPer clusters of galaxies using the Planck and AKARI data. The results show a clear () detection of a SZ signal, with an amplitude consistent with theoretical expectations from our numerical simulations. Interestingly, our analysis also shows a strongly correlated dust signal in both the high-frequency Planck and AKARI maps, with dust parameters of , .

In our analysis we have assumed that dust emission is described by a single-temperature modified blackbody. This assumption provides good fits to the data and we find no evidence to support using a more complex model. With only 4 frequency bands above , it is simply not possible to constrain the parameters of a more complicated model for dust emission reliably.

According to our simulations it will be difficult to measure the small () enhancement of the SZ signal associated with AGN feedback. With larger spectroscopic samples of QSOs (or galaxies that have hosted AGN in the past; see Spacek et al. 2016b) and CMB data with higher resolution and sensitivity (e.g. Benson et al., 2014; De Bernardis et al., 2016), the detection of a statistically significant SZ signal (disentangled from dust emission) is a more realistic goal. Together with numerical simulations, it may then be possible to accurately determine QSO host halo masses, which are poorly constrained at present.


The authors thank Nick Battaglia, Lindsey Bleem, Anthony Challinor, Suet-Ying Mak, and Glenn White for helpful discussions during this work, and further Tom Crawford and Kyle Story for helpful conversations about matched filtering during a previous study. BS acknowledges support from an Isaac Newton Studentship at the University of Cambridge and from the Science and Technology Facilities Council (STFC). TG acknowledges support from the Kavli Foundation and STFC grant ST/L000636/1. DS acknowledges support by the STFC and the ERC Starting Grant 638707 “Black holes and their host galaxies: coevolution across cosmic time”. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration, 2013). Some of the results in this paper have been derived using the HEALPix (Górski et al., 2005) package. Furthermore, the use of the following tools is acknowledged: healpy (, CosmoloPy (, and (

This research is based on observations obtained with Planck (, an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. Furthermore, it is in parts based on observations with AKARI, a JAXA project with the participation of ESA.

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


  • Afshordi et al. (2007) Afshordi N., Lin Y.-T., Nagai D., Sanderson A. J. R., 2007, MNRAS, 378, 293
  • Arnaud et al. (2010) Arnaud M., Pratt G. W., Piffaretti R., Böhringer H., Croston J. H., Pointecouteau E., 2010, A&A, 517, A92
  • Astropy Collaboration (2013) Astropy Collaboration 2013, A&A, 558, A33
  • Battaglia et al. (2010) Battaglia N., Bond J. R., Pfrommer C., Sievers J. L., Sijacki D., 2010, ApJ, 725, 91
  • Benson et al. (2014) Benson B. A., et al., 2014, in Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VII. p. 91531P (arXiv:1407.2973), doi:10.1117/12.2057305
  • Bhattacharya et al. (2008) Bhattacharya S., Di Matteo T., Kosowsky A., 2008, Mon. Not. Roy. Astron. Soc., 389, 34
  • Binney & Tabor (1995) Binney J., Tabor G., 1995, MNRAS, 276, 663
  • Birkinshaw (1999) Birkinshaw M., 1999, Phys.Rept., 310, 97
  • Blain et al. (2002) Blain A. W., Smail I., Ivison R. J., Kneib J.-P., Frayer D. T., 2002, Phys. Rep., 369, 111
  • Blandford (1990) Blandford R., 1990, in Courvoisier T.-L., Mayor M., eds, Active Galactic Nuclei. Saas-Fee Advanced Course 20. pp 161–269
  • Bleem (2013) Bleem L. E., 2013, PhD thesis, The University of Chicago
  • Bleem et al. (2015a) Bleem L. E., Stalder B., Brodwin M., Busha M. T., Gladders M. D., High F. W., Rest A., Wechsler R. H., 2015a, ApJS, 216, 20
  • Bleem et al. (2015b) Bleem L. E., et al., 2015b, ApJS, 216, 27
  • Boehringer et al. (1993) Boehringer H., Voges W., Fabian A. C., Edge A. C., Neumann D. M., 1993, MNRAS, 264, L25
  • Bower et al. (2006) Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645
  • Brooks & Gelman (1998) Brooks S. P., Gelman A., 1998, Journal of computational and graphical statistics, 7, 434
  • Burns (1990) Burns J. O., 1990, AJ, 99, 14
  • Calzetti et al. (2000) Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storchi-Bergmann T., 2000, ApJ, 533, 682
  • Carilli et al. (1994) Carilli C. L., Perley R. A., Harris D. E., 1994, MNRAS, 270, 173
  • Carlstrom et al. (2002) Carlstrom J. E., Holder G. P., Reese E. D., 2002, Ann.Rev.A&A, 40, 643
  • Cavaliere & Fusco-Femiano (1976) Cavaliere A., Fusco-Femiano R., 1976, A&A, 49, 137
  • Cen & Safarzadeh (2015) Cen R., Safarzadeh M., 2015, Astrophys. J., 809, L32
  • Chatterjee & Kosowsky (2007) Chatterjee S., Kosowsky A., 2007, Astrophys. J., 661, L113
  • Chatterjee et al. (2008) Chatterjee S., Di Matteo T., Kosowsky A., Pelupessy I., 2008, Mon. Not. Roy. Astron. Soc., 390, 535
  • Chatterjee et al. (2010) Chatterjee S., Ho S., Newman J. A., Kosowsky A., 2010, Astrophys. J., 720, 299
  • Chelouche et al. (2007) Chelouche D., Koester B. P., Bowen D. V., 2007, ApJ, 671, L97
  • Churazov et al. (2001) Churazov E., Brüggen M., Kaiser C. R., Böhringer H., Forman W., 2001, ApJ, 554, 261
  • Churazov et al. (2002) Churazov E., Sunyaev R., Forman W., Böhringer H., 2002, MNRAS, 332, 729
  • Clemens et al. (2013) Clemens M. S., et al., 2013, MNRAS, 433, 695
  • Crenshaw et al. (2003) Crenshaw D. M., Kraemer S. B., George I. M., 2003, ARA&A, 41, 117
  • Crichton et al. (2016) Crichton D., et al., 2016, MNRAS, 458, 1478
  • Croom et al. (2005) Croom S. M., et al., 2005, MNRAS, 356, 415
  • Croton et al. (2006) Croton D. J., et al., 2006, MNRAS, 365, 11
  • Dawson et al. (2013) Dawson K. S., et al., 2013, AJ, 145, 10
  • De Bernardis et al. (2016) De Bernardis F., et al., 2016, preprint, (arXiv:1607.02120)
  • DeGraf & Sijacki (2016) DeGraf C., Sijacki D., 2016, preprint, (arXiv:1609.06727)
  • Di Matteo et al. (2005) Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604
  • Diego & Partridge (2010) Diego J. M., Partridge B., 2010, MNRAS, 402, 1179
  • Doi et al. (2015) Doi Y., et al., 2015, PASJ, 67, 50
  • Dubois et al. (2010) Dubois Y., Devriendt J., Slyz A., Teyssier R., 2010, MNRAS, 409, 985
  • Dubois et al. (2016) Dubois Y., Peirani S., Pichon C., Devriendt J., Gavazzi R., Welker C., Volonteri M., 2016, MNRAS, 463, 3948
  • Eifler et al. (2015) Eifler T., Krause E., Dodelson S., Zentner A. R., Hearin A. P., Gnedin N. Y., 2015, MNRAS, 454, 2451
  • Eisenstein et al. (2011) Eisenstein D. J., et al., 2011, AJ, 142, 72
  • Fabian (1999) Fabian A. C., 1999, MNRAS, 308, L39
  • Fabian (2012) Fabian A. C., 2012, ARA&A, 50, 455
  • Fabian et al. (2006) Fabian A. C., Sanders J. S., Taylor G. B., Allen S. W., Crawford C. S., Johnstone R. M., Iwasawa K., 2006, MNRAS, 366, 417
  • Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  • Forman et al. (2007) Forman W., et al., 2007, ApJ, 665, 1057
  • Franceschini (2000) Franceschini A., 2000, ArXiv Astrophysics e-prints,
  • Gelman & Rubin (1992) Gelman A., Rubin D. B., 1992, Statistical science, pp 457–472
  • Genel et al. (2014) Genel S., et al., 2014, MNRAS, 445, 175
  • Goodman & Weare (2010) Goodman J., Weare J., 2010, Comm. App. Math. Comp. Sci., 5, 65
  • Górski et al. (2005) Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759
  • Gupta et al. (2016) Gupta N., et al., 2016, preprint, (arXiv:1605.05329)
  • Gutiérrez & López-Corredoira (2014) Gutiérrez C. M., López-Corredoira M., 2014, A&A, 571, A66
  • Haehnelt & Tegmark (1996) Haehnelt M. G., Tegmark M., 1996, MNRAS, 279, 545
  • Hajian et al. (2013) Hajian A., Battaglia N., Spergel D. N., Bond J. R., Pfrommer C., Sievers J. L., 2013, J. Cosmology Astropart. Phys., 11, 064
  • Hasselfield et al. (2013) Hasselfield M., et al., 2013, J. Cosmology Astropart. Phys., 7, 008
  • Helou & Walker (1988) Helou G., Walker D. W., eds, 1988, Infrared astronomical satellite (IRAS) catalogs and atlases. Volume 7: The small scale structure catalog   Vol. 7
  • Hill & Spergel (2014) Hill J. C., Spergel D. N., 2014, JCAP, 1402, 030
  • Ishihara et al. (2010) Ishihara D., et al., 2010, A&A, 514, A1
  • Kaiser (1986) Kaiser N., 1986, MNRAS, 222, 323
  • Kauffmann & Haehnelt (2000) Kauffmann G., Haehnelt M., 2000, MNRAS, 311, 576
  • Kellermann et al. (1989) Kellermann K. I., Sramek R., Schmidt M., Shaffer D. B., Green R., 1989, AJ, 98, 1195
  • King (2003) King A., 2003, ApJ, 596, L27
  • Lynden-Bell (1969) Lynden-Bell D., 1969, Nature, 223, 690
  • Madau & Dickinson (2014) Madau P., Dickinson M., 2014, ARA&A, 52, 415
  • Maiolino et al. (2012) Maiolino R., et al., 2012, MNRAS, 425, L66
  • Mak et al. (2016) Mak D. S. Y., Challinor A., Efstathiou G., Lagache G., 2016, preprint, (arXiv:1609.08942)
  • McCarthy et al. (2011) McCarthy I. G., Schaye J., Bower R. G., Ponman T. J., Booth C. M., Dalla Vecchia C., Springel V., 2011, MNRAS, 412, 1965
  • Melchior et al. (2016) Melchior P., et al., 2016, preprint, (arXiv:1610.06890)
  • Miley & De Breuck (2008) Miley G., De Breuck C., 2008, A&ARv, 15, 67
  • Miville-Deschênes & Lagache (2005) Miville-Deschênes M.-A., Lagache G., 2005, ApJS, 157, 302
  • Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  • Neugebauer et al. (1984) Neugebauer G., et al., 1984, ApJ, 278, L1
  • Omont et al. (2003) Omont A., Beelen A., Bertoldi F., Cox P., Carilli C. L., Priddey R. S., McMahon R. G., Isaak K. G., 2003, A&A, 398, 857
  • Pâris et al. (2016) Pâris I., et al., 2016, preprint, (arXiv:1608.06483)
  • Planck Collaboration (2011) Planck Collaboration 2011, A&A, 536, A12
  • Planck Collaboration (2013) Planck Collaboration 2013, A&A, 557, A52
  • Planck Collaboration (2014a) Planck Collaboration 2014a, A&A, 571, A11
  • Planck Collaboration (2014b) Planck Collaboration 2014b, A&A, 571, A28
  • Planck Collaboration (2015a) Planck Collaboration 2015a, preprint, (arXiv:1502.01586)
  • Planck Collaboration (2015b) Planck Collaboration 2015b, preprint, (arXiv:1502.01587)
  • Planck Collaboration (2015d) Planck Collaboration 2015d, preprint, (arXiv:1502.01588)
  • Planck Collaboration (2015c) Planck Collaboration 2015c, preprint, (arXiv:1507.02058)
  • Planck Collaboration (2016b) Planck Collaboration 2016b, preprint, (arXiv:1606.05120)
  • Planck Collaboration (2016a) Planck Collaboration 2016a, preprint, (arXiv:1605.09387)
  • Planck Collaboration (2016c) Planck Collaboration 2016c, A&A, 594, A1
  • Planck Collaboration (2016d) Planck Collaboration 2016d, A&A, 594, A2
  • Planck Collaboration (2016e) Planck Collaboration 2016e, A&A, 594, A27
  • Porciani et al. (2004) Porciani C., Magliocchetti M., Norberg P., 2004, MNRAS, 355, 1010
  • Priddey & McMahon (2001) Priddey R. S., McMahon R. G., 2001, MNRAS, 324, L17
  • Puchwein et al. (2008) Puchwein E., Sijacki D., Springel V., 2008, ApJ, 687, L53
  • Puchwein et al. (2010) Puchwein E., Springel V., Sijacki D., Dolag K., 2010, MNRAS, 406, 936
  • Quilis et al. (2001) Quilis V., Bower R. G., Balogh M. L., 2001, MNRAS, 328, 1091
  • Rees (1984) Rees M. J., 1984, ARA&A, 22, 471
  • Richardson et al. (2012) Richardson J., Zheng Z., Chatterjee S., Nagai D., Shen Y., 2012, ApJ, 755, 30
  • Ross et al. (2012) Ross N. P., et al., 2012, ApJS, 199, 3
  • Rozo et al. (2015) Rozo E., Rykoff E. S., Bartlett J. G., Melin J. B., 2015, MNRAS, 450, 592
  • Ruan et al. (2015) Ruan J. J., McQuinn M., Anderson S. F., 2015, Astrophys. J., 802, 135
  • Rupke & Veilleux (2011) Rupke D. S. N., Veilleux S., 2011, ApJ, 729, L27
  • Rykoff et al. (2012) Rykoff E. S., et al., 2012, ApJ, 746, 178
  • Rykoff et al. (2014) Rykoff E. S., et al., 2014, Astrophys. J., 785, 104
  • Rykoff et al. (2016) Rykoff E. S., et al., 2016, ApJS, 224, 1
  • Saro et al. (2015) Saro A., et al., 2015, MNRAS, 454, 2305
  • Saro et al. (2016) Saro A., et al., 2016, preprint, (arXiv:1605.08770)
  • Scannapieco et al. (2008) Scannapieco E., Thacker R. J., Couchman H. M. P., 2008, ApJ, 678, 674
  • Schäfer et al. (2006) Schäfer B. M., Pfrommer C., Hell R. M., Bartelmann M., 2006, MNRAS, 370, 1713
  • Schaye et al. (2010) Schaye J., et al., 2010, MNRAS, 402, 1536
  • Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  • Schneider et al. (2010) Schneider D. P., et al., 2010, AJ, 139, 2360
  • Sembolini et al. (2016) Sembolini F., et al., 2016, MNRAS, 459, 2973
  • Semboloni et al. (2011) Semboloni E., Hoekstra H., Schaye J., van Daalen M. P., McCarthy I. G., 2011, MNRAS, 417, 2020
  • Sijacki et al. (2007) Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MNRAS, 380, 877
  • Sijacki et al. (2015) Sijacki D., Vogelsberger M., Genel S., Springel V., Torrey P., Snyder G. F., Nelson D., Hernquist L., 2015, MNRAS, 452, 575
  • Silk & Rees (1998) Silk J., Rees M. J., 1998, A&A, 331, L1
  • Simet et al. (2016) Simet M., McClintock T., Mandelbaum R., Rozo E., Rykoff E., Sheldon E., Wechsler R. H., 2016, preprint, (arXiv:1603.06953)
  • Smith et al. (2012) Smith M. W. L., et al., 2012, ApJ, 748, 123
  • Soergel et al. (2016) Soergel B., et al., 2016, MNRAS, 461, 3172
  • Spacek et al. (2016a) Spacek A., Scannapieco E., Cohen S., Joshi B., Mauskopf P., 2016a, preprint, (arXiv:1610.02068)
  • Spacek et al. (2016b) Spacek A., Scannapieco E., Cohen S., Joshi B., Mauskopf P., 2016b, ApJ, 819, 128
  • Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  • Springel (2010) Springel V., 2010, MNRAS, 401, 791
  • Springel et al. (2001) Springel V., White M., Hernquist L., 2001, ApJ, 549, 681
  • Sturm et al. (2011) Sturm E., et al., 2011, ApJ, 733, L16
  • Sunyaev & Zeldovich (1970) Sunyaev R. A., Zeldovich Y. B., 1970, Ap&SS, 7, 3
  • Sunyaev & Zeldovich (1972) Sunyaev R. A., Zeldovich Y. B., 1972, Comments on Astrophysics and Space Physics, 4, 173
  • Takita et al. (2015) Takita S., et al., 2015, PASJ, 67, 51
  • Tegmark & de Oliveira-Costa (1998) Tegmark M., de Oliveira-Costa A., 1998, ApJ, 500, L83
  • Tombesi et al. (2015) Tombesi F., Meléndez M., Veilleux S., Reeves J. N., González-Alfonso E., Reynolds C. S., 2015, Nature, 519, 436
  • Tormen et al. (1997) Tormen G., Bouchet F. R., White S. D. M., 1997, MNRAS, 286, 865
  • Verdier et al. (2016) Verdier L., Melin J.-B., Bartlett J. G., Magneville C., Palanque-Delabrouille N., Yèche C., 2016, A&A, 588, A61
  • Vogelsberger et al. (2014) Vogelsberger M., et al., 2014, Nature, 509, 177
  • White et al. (1997) White R. L., Becker R. H., Helfand D. J., Gregg M. D., 1997, ApJ, 475, 479
  • White et al. (2012) White M., et al., 2012, MNRAS, 424, 933
  • York et al. (2000) York D. G., et al., 2000, AJ, 120, 1579
  • de Haan et al. (2016) de Haan T., et al., 2016, preprint, (arXiv:1603.06522)
  • van Daalen et al. (2011) van Daalen M. P., Schaye J., Booth C. M., Dalla Vecchia C., 2011, MNRAS, 415, 3649

Appendix A Additional data and data quality tests

In this appendix, we describe the additional FIR data that we have used for testing purposes. We further provide details about the tests that we have performed on the microwave and FIR data at the various frequencies before including them into our analysis.

a.1 IRAS/IRIS data

IRAS has surveyed the infrared sky in four bands, with band centres between 12 and 100 , and a resolution of approximately . Here we make use of the Improved Reprocessing of the IRAS Survey (IRIS) maps by Miville-Deschênes & Lagache (2005), which — amongst others — feature improvements in calibration and destriping over the original IRAS maps.999The IRIS maps are available at both as cutouts and in the Healpix format with . The IRAS 100  band is already measuring the falling part of the dust spectrum. The higher-frequency bands of IRIS are more sensitive to dust at higher temperatures and therefore do not add additional information to our measurement. As additionally the 60, 25 and 12  bands show a stronger contamination by residual zodiacal light, we only consider the IRIS 100  map.

a.2 Random points null test

We perform a null test by stacking on random positions in the same footprint. The latter are generated as follows: (1) We generate a low-resolution Healpix map of the QSO density, similar to the one in Fig. 3. (2) Next, we create a large number of ‘candidate points’ which are randomly distributed on the sphere. (3) We finally select final points randomly from the candidates. The probability of any given candidate point to be selected is given by the appropriately normalized QSO surface density at its position in the original catalogue.

We then process the random points through the same measurement pipeline as the original catalogue and show the resulting flux density measurement in Fig. 13. From the results of this tests, we draw the following conclusions: all Planck frequency points pass the null test; the same holds true for the AKARI 90  band. The results from the AKARI 160, 140, and 65  bands are also consistent with no bias, but these bands show large scatter between individual realizations, reflecting their high noise levels and significant calibration uncertainties. This provides additional justification to our decision not to include them in the main analysis.

Most notably, however, the IRIS 100  band fails the null test, revealing a statistically significant offset of around . Amongst the possible reasons for this result are residual striping errors or calibration uncertainties. However, a more in-depth investigation of this, for example by a detailed comparison of AKARI and IRIS at the map level, is beyond the scope of this paper. We therefore do not include the IRIS data in our main analysis. This does not have a significant impact on our ability to constrain the dust parameters, as there is significant overlap in the instrumental bandpasses of the IRIS 100  and AKARI 90  bands (see Fig. 2).

Figure 13: Random points null test: we show here the results of the stacking on random positions in the map. The blue points show one representative realization, with uncertainties estimated as in the original measurement. The shaded dashed lines indicate the results of 60 further random realizations.

It is worth noting that the successful null test at low frequencies also demonstrates that our result for the Compton- parameter is not significantly affected by a stacking bias. The latter can, in principle, arise because the measured receives contributions from all objects along the same line of sight, not only from the QSO (or cluster) in our catalogue. However, this would result in a signal also in sightlines that do not contain such an object. Our null test demonstrates that this bias is well below the statistical uncertainties in our analysis.

a.3 Impact of discarded data

We next demonstrate the impact that the discarded data can have on our measurement. As a reminder, we have excluded the AKARI 160, 140 and 90 bands because of their uncertain calibration at low intensities (see Table 1 and Takita et al. 2015). Furthermore, we have not used the IRIS 100 band because it does not pass the stacking null test. Here we now repeat our analysis including these frequencies.

Adding only the IRIS 100 band only marginally affects our result, as it is consistent with the SED we have determined from the seven Planck bands and AKARI 90 . When additionally adding the three discarded AKARI bands, we assume that their flux calibration uncertainty is twice the lower limit we quote in Table 1. Despite their large uncertainties, we find that including these three frequencies pulls the dust solution away from the best fit in the main analysis: we now obtain and . This changes the shape of the dust spectrum at low frequencies, so that now there is room for an SZ component with , or . At the same time, however, the fit both to the higher-frequency Planck points and the lowest-frequency AKARI points is poor, which would already cause us to reject this solution on goodness-of-fit considerations only. We show this in Fig. 14. Therefore, we stress that this should not be interpreted as a significant detection of an SZ signal when adding additional data. Rather, we consider it a demonstration how sensitive the SZ result is on the dust solution, and how easily the latter is affected by the quality of the used data.

Figure 14: Impact of the discarded data on our analysis: we show here the best-fitting SED if the discarded frequencies are included in the analysis. We denote the latter by the open triangles, from low to high frequencies they are: AKARI 160 and 140 , IRIS 100 , and AKARI 65 . Including these frequencies significantly affects the dust solution, and leaves more room for an SZ component. However, the quality of the fit is poor.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question