CONSTRAINING AGN FEEDBACK IN MASSIVE ELLIPTICALS WITH SOUTH POLE TELESCOPE MEASUREMENTS OF THE THERMAL SUNYAEV–ZEL'DOVICH EFFECT

The tSZ effect describes the process by which CMB photons gain energy when passing through ionized gas (Sunyaev & Zeldovich 1970, 1972). The photons are shifted to higher energies by thermally energetic electrons through inverse Compton scattering, and the resulting CMB anisotropy has a distinctive frequency dependence which causes a deficit of photons at frequencies below ${\nu }_{{\rm{null}}}=217.6\;{\rm{GHz}}$ and an excess of photons above ${\nu }_{{\rm{null}}}$, with no change at ${\nu }_{{\rm{null}}}$. For the nonrelativistic plasma we will be interested in here, the change in CMB temperature as a function of frequency due to the tSZ effect is

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}T}{{T}_{{\rm{CMB}}}}=y\left(x\displaystyle \frac{{e}^{x}+1}{{e}^{x}-1}-4\right),\end{eqnarray} \tag{ 1 }$$

where the dimensionless Compton-y parameter is defined as

$$\begin{eqnarray}&&y\equiv \int {dl}\;{\sigma }_{{\rm{T}}}\displaystyle \frac{{n}_{e}k({T}_{e}-{T}_{{\rm{CMB}}})}{{m}_{e}{c}^{2}},\end{eqnarray} \tag{ 2 }$$

where ${\sigma }_{{\rm{T}}}$ is the Thomson cross section, k is the Boltzmann constant, m_e is the electron mass, c is the speed of light, n_e is the electron number density, T_e is the electron temperature, ${T}_{{\rm{CMB}}}$ is the CMB temperature (we use ${T}_{{\rm{CMB}}}=2.725$ K), and the integral is performed over the line of sight distance l. Finally, the dimensionless frequency x is given by

$$\begin{eqnarray}&&x\equiv \displaystyle \frac{h\nu }{{{kT}}_{{\rm{CMB}}}}=\displaystyle \frac{\nu }{56.81\;{\rm{GHz}}},\end{eqnarray} \tag{ 3 }$$

where h is the Planck constant.

We can calculate the total excess thermal energy associated with a source by integrating Equation (2) over a region of sky around the source as (e.g., Scannapieco et al. 2008; Ruan et al. 2015),

$$\begin{eqnarray}\displaystyle \int d{\boldsymbol{\theta }}\;y({\boldsymbol{\theta }}) & = & \displaystyle \int d{\boldsymbol{\theta }}\displaystyle \int {dl}\;{\sigma }_{{\rm{T}}}\displaystyle \frac{{n}_{e}{{kT}}_{e}}{{m}_{e}{c}^{2}}\\ & = & \displaystyle \frac{{\sigma }_{{\rm{T}}}}{{m}_{e}{c}^{2}}{l}_{{\rm{ang}}}^{-2}\displaystyle \int {dV}\;{n}_{e}{{kT}}_{e},\end{eqnarray} \tag{ 4 }$$

where ${\boldsymbol{\theta }}$ is a two-dimensional vector in the plane of the sky in units of radians, ${l}_{{\rm{ang}}}$ is the angular diameter distance to the source, V is the volume of interest around the source, and we have restricted our attention to hot gas with ${T}_{e}\gg {T}_{{\rm{CMB}}}.$ In Equation (4), the Compton-y integral has become a volume integral of the electron pressure (i.e., ${P}_{e}={n}_{e}{{kT}}_{e}$), which is related to the associated thermal energy as

$$\begin{eqnarray}&&\int {dV}\;{n}_{e}{{kT}}_{e}=\left(\displaystyle \frac{2}{3}\right)\left(\displaystyle \frac{1+A}{2+A}\right){E}_{{\rm{therm}}},\end{eqnarray} \tag{ 5 }$$

where $A=0.08$ is the cosmological number abundance of helium, and ${E}_{{\rm{therm}}}$ is the total thermal energy associated with the source: that gained from the initial collapse of the baryons, plus the contribution from the AGN, minus the losses due to cooling and the ${PdV}$ work done during expansion. We can combine Equations (4) and (5) and solve for ${E}_{{\rm{therm}}}$ to get

$$\begin{eqnarray}{E}_{{\rm{therm}}} & = & 2.9\displaystyle \frac{{m}_{e}{c}^{2}}{{\sigma }_{{\rm{T}}}}{l}_{{\rm{ang}}}^{2}\displaystyle \int d{\boldsymbol{\theta }}y({\boldsymbol{\theta }})\\ & = & 2.9\times {10}^{60}\;{\rm{erg}}\;{\left(\displaystyle \frac{{l}_{{\rm{ang}}}}{{\rm{Gpc}}}\right)}^{2}\displaystyle \frac{\displaystyle \int d{\boldsymbol{\theta }}y({\boldsymbol{\theta }})}{{10}^{-6}\;{{\rm{arcmin}}}^{2}}.\end{eqnarray} \tag{ 6 }$$

Finally, we can combine Equations (1) and (6) to get ${E}_{{\rm{therm}}}$ in terms of ${\rm{\Delta }}T$ at a given dimensionless frequency x,

$$\begin{eqnarray}&&{E}_{{\rm{therm}}}=\displaystyle \frac{1.1\times {10}^{60}\;{\rm{erg}}}{x\frac{{e}^{x}+1}{{e}^{x}-1}-4}\;{\left(\displaystyle \frac{{l}_{{\rm{ang}}}}{{\rm{Gpc}}}\right)}^{2}{\displaystyle \frac{\int {\rm{\Delta }}T({\boldsymbol{\theta }})d{\boldsymbol{\theta }}}{\mu {\rm{K}}\text{arcmin}}}^{2}.\end{eqnarray} \tag{ 7 }$$

To compare the energies and angular sizes above with the expectations from models of feedback, we can construct a simple model of gas heating with and without AGN feedback. To do this we first compute ${R}_{{\rm{vir}}},$ the virial radius of a (spherical) dark matter halo defined as the physical radius within which the density is 200 times the mean cosmic value. As a function of redshift z and mass $M,$ this is

$$\begin{eqnarray}{R}_{{\rm{vir}}} & = & {\left[\displaystyle \frac{M}{(4\pi /3)200{{\rm{\Omega }}}_{0}{\rho }_{{\rm{crit}}}{(1+z)}^{3}}\right]}^{1/3}\\ & = & 0.67\;{\rm{Mpc}}\;{M}_{13}^{1/3}{(1+z)}^{-1},\end{eqnarray} \tag{ 8 }$$

where ${\rho }_{{\rm{crit}}}$ is the critical density at z = 0, and M₁₃ is the mass of the halo in units of ${10}^{13}{M}_{\odot }.$ This can be compared to the angular scales above, using the fact that at an angular diameter distance of 1 Gpc, 1 arcmin corresponds to 0.29 Mpc.

If the gas collapses and virializes along with the dark matter, it will be shock-heated during gravitational infall to the virial temperature,

$$\begin{eqnarray}&&{T}_{{\rm{vir}}}=\displaystyle \frac{{GM}}{{R}_{{\rm{vir}}}}\displaystyle \frac{\mu {m}_{p}}{2k}=2.4\times {10}^{6}\;{\rm{K}}\;{M}_{13}^{2/3}(1+z),\end{eqnarray} \tag{ 9 }$$

where G is the gravitational constant, m_p is the proton mass, and $\mu =0.62$ is the average particle mass in units of ${m}_{p}.$ If we approximate the gas distribution as isothermal at this temperature, its total thermal energy can then be estimated as

$$\begin{eqnarray}{E}_{{\rm{therm,gravity}}} & = & \displaystyle \frac{3{{kT}}_{{\rm{vir}}}}{2}\displaystyle \frac{{{\rm{\Omega }}}_{b}}{{{\rm{\Omega }}}_{0}}\displaystyle \frac{M}{\mu {m}_{p}}\\ & = & 1.5\times {10}^{60}\;{\rm{erg}}\;{M}_{13}^{5/3}(1+z).\end{eqnarray} \tag{ 10 }$$

To relate the stellar masses of the galaxies we will be stacking to the dark matter halo masses, we can take advantage of the observed relation between black hole mass and halo circular velocity, ${v}_{c},$ from Ferrarese (2002; see also Merritt & Ferrarese 2001; Tremaine et al. 2002), and convert the black hole mass to its corresponding bulge dynamical mass using a factor of 400 (Marconi & Hunt 2003). This gives

$$\begin{eqnarray}{M}_{{\rm{stellar}}} & = & {6.6}_{-3.2}^{+5.5}\times {10}^{10}\;{M}_{\odot }\;{\left(\displaystyle \frac{{v}_{{\rm{c}}}}{300\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{5}\\ & = & {2.8}_{-1.4}^{+2.4}\times {10}^{10}{M}_{\odot }\;{M}_{13}^{5/3}{(1+z)}^{5/2},\end{eqnarray} \tag{ 11 }$$

where we have used the fact that ${v}_{c}={({GM}/{R}_{{\rm{vir}}})}^{1/2}=254\;\mathrm{km}\;{{\rm{s}}}^{-1}\;{M}_{13}^{1/3}{(1+z)}^{1/2},$ and taken ${M}_{{\rm{stellar}}}\propto {v}_{c}^{{\alpha }_{c}}$ with the power law index ${\alpha }_{c}=5,$ which is near the center of the allowed range of $5.4\pm 1.1,$ and we take our uncertainties from Ferrarese (2002). Substituting Equation (11) into (10) gives

$$\begin{eqnarray}&&{E}_{{\rm{therm,gravity}}}={5.4}_{-2.9}^{+5.4}\times {10}^{60}\;{\rm{erg}}\;\displaystyle \frac{{M}_{{\rm{stellar}}}}{{10}^{11}{M}_{\odot }}{(1+z)}^{-3/2}.\end{eqnarray} \tag{ 12 }$$

This is the total thermal energy expected around a galaxy of stellar mass ${M}_{{\rm{stellar}}}$ due purely to gravitational heating, and ignoring both radiative cooling, which will decrease ${E}_{{\rm{therm}}},$ and AGN feedback, which will increase it.

While there are many models of AGN feedback, each of which will lead to somewhat different signatures in our data, we can estimate the overall magnitude of this effect by making use of the simple model described in Scannapieco & Oh (2004) (see also Thacker et al. 2006; Scannapieco et al. 2008). In this case, AGN feedback is described as tapping into a small fraction, ${\epsilon }_{k},$ of the total bolometric luminosity of the AGN to heat the surrounding gas. In particular, black holes are assumed to shine at the Eddington luminosity ($1.2\times {10}^{38}$ erg ${{\rm{s}}}^{-1}$ ${M}_{\odot }^{-1}$ ) for a time $0.035\;{t}_{{\rm{dynamical}}}$, where

$$\begin{eqnarray}&&{t}_{{\rm{dynamical}}}\equiv {R}_{{\rm{vir}}}/{v}_{c}=2.6\;{\rm{Gyr}}\;{(1+z)}^{-3/2}.\end{eqnarray} \tag{ 13 }$$

This choice of timescale gives a good match to the observed evolution of the quasar luminosity function (Wyithe & Loeb 2002, 2003; Scannapieco & Oh 2004). This gives

$$\begin{eqnarray}&&{E}_{{\rm{therm,feedback}}}=4.1\times {10}^{60}\;{\rm{erg}}\;{\epsilon }_{k,0.05}\;\displaystyle \frac{{M}_{{\rm{stellar}}}}{{10}^{11}{M}_{\odot }}\;{(1+z)}^{-3/2}.\end{eqnarray} \tag{ 14 }$$

Here ${\epsilon }_{k,0.05}\equiv {\epsilon }_{k}/0.05$, such that the kinetic energy input is normalized to a typical value needed to achieve antiheirarchical galaxy evolution through effective feedback (e.g., Scannapieco & Oh 2004; Thacker et al. 2006; Costa et al. 2014). Note that the uncertainty in this equation is completely dominated by the ${\epsilon }_{k,0.05}$ term, which is uncertain to within an order of magnitude.

This energy input is equal in magnitude to the errors in ${E}_{{\rm{therm,gravity}}},$ meaning that the differences between models with and without AGN feedback will not be dramatic. Thus only detailed simulations will be able to make precise predictions at the level needed to rule out or lend support to a particular model of AGN feedback. Although carrying out such simulations is beyond the scope of this paper, Equations (12) and (14) are roughly consistent with such sophisticated models (e.g., Thacker et al. 2006; Chatterjee et al. 2008), meaning that they can be used as an approximate guide to interpreting our results. Thus, we will use them to provide a general context for thinking about our observational results in terms of AGN feedback.

Finally, we note that the sound speed c_s of the gas in the gravitationally heated case is similar to the circular velocity (i.e., ${c}_{s}={[\gamma {kT}/(\mu {m}_{p})]}^{1/2}={(\gamma /2)}^{1/2}{v}_{c}$, where γ is the adiabatic index), and the expected energy input from the AGN is similar to the energy input from gravitational heating. This means that the energy input from the AGN will take a timescale $\approx {t}_{{\rm{dynamical}}}$ to impact gas on the scale of the halo, and it is unlikely to affect scales much larger than $\approx 2{R}_{{\rm{vir}}}.$ These sizes and timescales mean that at the moderate redshifts we will be exploring, the majority of the gas heating we are interested in will occur on scales $\lesssim 2$ arcmin.

The BCS (Desai et al. 2012) was a National Optical Astronomy Observatory (NOAO) Large Survey project that observed ≈80 deg² of the southern sky over 60 nights between 2005 November and 2008 November on the 4 m Víctor M. Blanco telescope at the Cerro Tololo Inter-American Observatory in Chile using the Mosaic II imager with g, r, i, and z bands. The filter centers, effective widths, and magnitude limits are given in Table 1.

The BCS data is split up into many smaller 36 × 36 arcmin (8192 × 8192 pixel) images called tiles, with $\approx 1$ arcmin overlap between neighboring tiles. Each pixel subtends 0.27 arcsec on the sky. As described in Desai et al. (2012), each raw data tile is put through a detrending pipeline, which consists of crosstalk corrections, overscan, flatfield, bias and illumination correction, and astrometric calibration. The average FWHM of the seeing disk in the single epoch images ranges between 0.7 and 1.6 arcsec. Each tile is then put through a co-addition pipeline that combines data taken over the same locations on the sky to build deeper single images. This results in a co-added tile image and an inverse-variance weightmap (confidence image) for each tile region of the survey. We use the area of the BCS that overlaps with the VHS and SPT-SZ, known as the 5 hr field (referring to its right ascension). This dataset consists of 135 tiles and their associated weightmaps for each band, covering ≈45 deg².

The VHS (McMahon 2012) is a large-scale near-infrared survey whose goal is to survey the entire southern celestial hemisphere (≈20,000 deg²). The survey component in which we are interested is called the VHS DES (DES because it overlaps with the Dark Energy Survey), a 5000 deg² region that is imaged with 120 s exposure times in the J, H, and K_s bands (see Table 1). The data was obtained from 2009 to 2011 on the 4.1 m Visible and Infrared Survey Telescope for Astronomy (VISTA) at the Paranal Observatory in Chile.

The VHS data is split up into ≈2 × 1.5 deg ($\approx$ 12,770 × 15,660 pixel) tiles. Each pixel subtends 0.33 arcsec on the sky. As described on the VISTA data processing web page¹ , the raw VHS data go through a pipeline that involves reset correction, dark correction, linearity correction, flat field correction, sky background correction, destripe, jitter stacking, astrometric and photometric calibration, and tile generation. Tiles are generated from six smaller, stacked pawprints, each containing 16 even smaller detector-level images² , and the median image seeing as measured from stellar FWHM on VHS pawprints ranges from 0.89 arcsec in K_s to 0.99 arcsec in J. The stacked paw prints then result in a science-ready tile image and inverse-variance weightmap for each tile region of the survey. We are interested in the area of the VHS that overlaps with the BCS (see Figure 1). This results in 20 tiles and their associated weightmaps, covering ≈55 deg².

The SPT-SZ survey (Schaffer et al. 2011) used the 10 m SPT at the National Science Foundation's (NSF) Amundsen–Scott South Pole Station to survey a large area of the sky at millimeter and sub-millimeter wavelengths with arcminute angular resolution and low noise (Ruhl et al. 2004; Padin et al. 2008; Carlstrom et al. 2011). The survey observed 2500 deg² of the southern sky during the austral winter seasons of 2008 through 2011. Data from the 2011 release that we are using covers ≈95 deg² to a depth of 17 and 41 μK arcmin at 150 GHz and 220 GHz, respectively, centered at (R.A., decl.) = (82.7, −55) degrees (see Table 1).

The SPT-SZ data is contained in a single image per band, ≈20° × 10° ($\approx 3000$ × 3000 pixels) projected as either a Sanson–Flamsteed projection or an oblique Lambert equal-area azimuthal projection. The Sanson–Flamsteed projection is most useful for cluster-finding and contains masked point-sources, while the Lambert projection is most useful for point-source analysis. Since we are interested in individual sources that are undetected and within the noise level, we use the Sanson–Flamsteed projection with point-sources masked. Each pixel subtends 15 arcsec on the sky. As described in Schaffer et al. (2011), the raw data goes through a pre-processing stage where the data from a single observation of the field is calibrated, data selection cuts are applied, and initial filtering and instrument characterization are performed. A map-making stage with additional filtering is performed on the pre-processed data and the data are binned into single-observation maps used for final co-adds. The final data products include a co-added image, two-dimensional beam functions, filter transfer functions, and noise power spectral densities for each band.

It is worth noting that 220 GHz is very close to the frequency at which there is no change in the CMB due to the tSZ effect (${\nu }_{{\rm{null}}}$ = 217.6 GHz), while 150 GHz, which is close to the peak of the undistorted CMB spectrum (160 GHz), will see a decrement in radiation. Equations (1) and (3) can be rewritten for these bands, though the equations must now involve integration over the SPT filter curves. Once this is done, we can write the Compton-y parameter as

$$\begin{eqnarray}&&y=-0.41\;\displaystyle \frac{{\rm{\Delta }}{T}_{150}}{1\;{\rm{K}}}\quad {\rm{and}}\quad y=9.9\;\displaystyle \frac{{\rm{\Delta }}{T}_{220}}{1\;{\rm{K}}},\end{eqnarray} \tag{ 15 }$$

where ${\rm{\Delta }}{T}_{150}$ and ${\rm{\Delta }}{T}_{220}$ are the temperature anisotropies at 150 and 220 GHz, respectively. Here we can explicitly see that, for the same y, the increase in ${\rm{\Delta }}{T}_{220}$ is about 24 times less than the decrease in ${\rm{\Delta }}{T}_{150}$. A measurement of the tSZ effect is therefore expected to give us a clear decrement at 150 GHz and no detectable change at 220 GHz.

We can also use Equations (3) and (7), integrated over the SPT filter curve as mentioned above, to compare the tSZ decrement at 150 GHz to the total thermal energy of an object as

$$\begin{eqnarray}&&{E}_{{\rm{therm}}}=-1.2\times {10}^{60}\;{\rm{erg}}{\left(\displaystyle \frac{{l}_{{\rm{ang}}}}{{\rm{Gpc}}}\right)}^{2}\displaystyle \frac{\int {\rm{\Delta }}{T}_{150}({\boldsymbol{\theta }})d{\boldsymbol{\theta }}}{\mu {\rm{K}}{\text{arcmin}}^{2}}.\end{eqnarray} \tag{ 16 }$$

Given the arcminute angular resolution and 17 μK arcmin sensitivity of the SPT 150 GHz data, this means that for stacks of several thousand sources we can hope to derive constraints on the order of ${\rm{\Delta }}{E}_{{\rm{therm}}}\approx {10}^{60}$ erg. This is sufficient to derive constraints that are interesting for discriminating between models of AGN feedback, as discussed in Section 2.2.

As seen in Figure 2, the sources in which we are interested in are brightest in the K_s band, and thus we use it to make all of our detections. Because the BCS and VHS tiles are different sizes and in different locations (see Figure 1), we consider every possible overlap between images when aligning the other data to the K_s tiles. We then match pixel sizes and locations, and to insure that fixed aperture flux measurements are consistent between bands, we also match the seeing between the K_s tiles and the other bands.

If the K_s tile has worse seeing than the other band, we simply degrade the other image with a Gaussian filter until it matches the FWHM of the K_s image. On the other hand, if the K_s image has better seeing, we degrade it to match the other band, compute the ratio of 3 arcsec diameter aperture fluxes between the two bands described below, and finally multiply the ratio by the 3 arcsec diameter flux measured from the unconvolved ${K}_{s}$ image. That is, we compute and apply an aperture correction as Flux_grizJH = Flux${}_{{grizJH},0}$ × (Flux${}_{{K}_{s},0}$/Flux${}_{{K}_{s},{\rm{degraded}}}$), where 0 denotes the non-degraded measurement. This is done to preserve both accurate colors and the best possible K_s flux in every case, since that is the most important band for our purposes. Note that, since the seeing is nearly the same for all tiles ($\approx 1$ arcsec, see Table 1), this correction is minor, with a mean ratio (Flux${}_{{K}_{s},0}$/Flux${}_{{K}_{s},{\rm{degraded}}}$) of 1.026 across all tiles.

To detect and measure every object in our field, we use the SExtractor software package, version 2.8.6 (Bertin & Arnouts 1996).^{3 ,4} The code detects and measures sources in an image through the following five-step process: (i) it creates a background map that estimates the noise at every pixel in the image; (ii) it detects sources using a thresholding technique; (iii) it uses a multiple isophotal analysis technique to deblend objects; (iv) it throws out spurious detections made in the wings of larger objects; and (v) it estimates the flux of each remaining object. Each of these steps can be adjusted by the user through configuration parameters, and we list our choice of these parameters for both BCS and VHS tiles in Table 2.

In all cases, we use SExtractor's dual-image mode, which allows us to make flux measurements in all bands from the same source detections in the K_s band. This results in a catalog of K_s-detected sources with MAG_AUTO and 3 arcsec diameter aperture flux measurements in every band. We use corrected MAG_AUTO for our final catalogs and the 3 arcsec diameter aperture MAG_APER to compute aperture corrections as described in Section 5.1. Finally, the overlap between tiles within both the BCS and VHS images leads to some sources being detected in multiple tiles. To correct for this, we match our catalog with itself and remove multiple occurrences of sources within 1 arcsec of each other. At this point in the analysis, our full catalog contains 565,561 sources, 168,944 (30%) of which are identified as duplicates. This leaves 396,617 total sources.

To confirm the reliability of our measurements, we compare our J, H, and K_s magnitudes with the source catalog released with the public VHS data. In particular, we select stars from our catalog using Equation (17), with < instead of $\geqslant ,$ and SExtractor FLAGS = 0 and match them with a random 10,000 source subset of the pre-made VHS catalog, where we define a match as two sources within 0.5 arcsec of each other. Note that our magnitudes are measured within 3 arcsec diameter apertures while the pre-made catalog uses 2.83 arcsec diameter apertures. A plot of the difference in magnitudes between our catalog and the pre-made VHS catalog can be seen in Figure 3. To remove extreme outliers, magnitudes from both catalogs are cut at the depths given in Table 1. The mean offsets from 0 in the magnitude differences are −0.11, −0.05, and −0.03 mag for J, H, and K_s, respectively. The mean photometric uncertainty in those offsets across all magnitudes are ±0.09, ±0.12, and ±0.18 mag, respectively. The solid lines in Figure 3 represent the mean uncertainty within magnitude bins of width 1 mag, and they are plotted as positive offsets from 0 on the y-axis. The uncertainty includes any differences in the measurement process between this paper and McMahon (2012), as well as the inherent uncertainty in the SExtractor measurements. As can be seen, the difference in magnitudes is reasonably within the uncertainty.

Zoom In Zoom Out Reset image size

We are not able to carry out a similar comparison for the BCS bands because the data do not come with reliable zeropoints, which are required to convert the measured image-level fluxes to actual fluxes. Instead, we compute the BCS band zeropoints ourselves using the stellar locus regression (SLR) code Big MACS (Kelly et al. 2014).⁵ This code calibrates the photometric zeropoints by creating a model stellar locus for every input filter and fitting them simultaneously to a selection of input stars. To input the best possible selection of stars for this purpose, we use a combination of several criteria that are fine-tuned for each tile to balance between quality and quantity. These include selecting stars using Equation (17), SExtractor FLAGS = 0, CLASS_STAR $\;\geqslant$ 0.9, A_IMAGE/B_IMAGE $\;\geqslant$ 0.8, FWHM_IMAGE within a certain range from the point-source limit, and selecting bright, but unsaturated fluxes.

Our star selection results in a mean of 525 stars used per tile. We run the code using the BCS bands (g, r, i, z) plus J and H. Since the VHS bands (J, H, K_s) already have accurate zeropoints, we use the code to compute the zeropoints of the other 5 bands relative to H. This allows us to do an independent check on the code by comparing the code's value for the J zeropoint with the actual J zeropoint. We find that the mean difference between the two is 0.0078 mag, which is close to the uncertainty of the code. The mean uncertainties in the derived zeropoint calibrations are 0.043, 0.037, 0.018, 0.012, and 0.0052 mag for g, r, i, z, and J, respectively.

Having obtained a calibrated catalog of sources, we then apply an initial set of cuts to remove cases that are too uncertain to be suitable for stacking. Our goal here is not to select a statistically complete set of large, old, passive, $z\geqslant 0.5$ galaxies in the survey area, but rather to select a subset of such galaxies that can be cleanly identified. To count any source as reliable, we first require that it triggers no SExtractor output flags (FLAGS = 0). This choice excludes: (i) sources that have neighbors bright enough to bias the photometry; (ii) sources that were originally blended with another source; (iii) sources with at least one saturated pixel; (iv) sources with incomplete or corrupted data; and (v) sources for which a memory overflow occurred when measuring their flux. Furthermore, we remove all sources with a measured FLUXERR_APER $\;\leqslant \quad 0$ in any band, and any source within 3 × FWHM_IMAGE from the edge of a tile, since the data become unreliable near these boundaries due to dithering.

Next, we separate stars from galaxies by making use of the ${{gzK}}_{s}$ method given by Equation (17). As in the plot of model galaxies (lower-right panel of Figure 2), our data (Figure 4) shows a clear division between the galaxy locus and the stellar locus along this limit. Note however that Arcila-Osejo & Sawicki (2013) proposed a star cut of $(z-{K}_{s})\lt 0.45(g-z)-0.57$, which differs from ours slightly. Furthermore, we apply Equation (18) to separate out young, lower-redshift galaxies from the $z\geqslant 0.5$ old, passive galaxies we are interested in. After applying these criteria, we are left with a catalog of 332,037 sources consisting of 123,567 stars (37%) and 208,470 galaxies (63%), 195,426 (59%) of which satisfy Equation (18). We then correct for Galactic dust extinction using the Schlegel et al. (1998) dust map and the extinction curve of Fitzpatrick & Massa (1999). Source-count histograms of the K_s magnitudes for stars and the corrected K_s magnitudes for galaxies are shown by the solid black and dashed blue lines, respectively, in Figure 6.

Zoom In Zoom Out Reset image size

With our catalog of galaxies we use the EAZY software package (Brammer et al. 2008) to estimate photometric redshifts and the FAST software package (Kriek et al. 2009) to estimate various characteristics such as redshift, age, mass, and star formation rate (SFR). First, EAZY steps through a grid of redshifts, fits linear combinations of template spectra to our photometric data, and ultimately finds the best estimate for redshift, including optional flux- and redshift-based priors. We allow for fits to make use of linear combinations of up to two of the default template spectra, and also apply the default K-band flux- and redshift-based prior derived from the GOODS-Chandra Deep Field-South (Wuyts et al. 2008).

The resulting redshifts are then fed into the FAST code, along with our seven-band photometric data, to fit for six additional parameters: age, mass, star formation timescale τ, SFR, dust content, and metallicity. FAST allows for a range of parameters when generating model fluxes, and in this analysis we choose: (i) a stellar population synthesis model as in Conroy & Gunn (2010); (ii) a Chabrier (2003) initial stellar mass function; (iii) an exponentially declining star formation history $\propto \mathrm{exp}(-t/\tau );$ and (iv) a dust extinction law as given by Kriek & Conroy (2013). To determine the best-fit parameters, the code simply determines the ${\chi }^{2}$ of every point of the model cube and finds the minimum. While the code allows for confidence intervals calibrated using Monte-Carlo simulations, here we simply make use of the best-fit values for each galaxy, recording its ${\chi }^{2}$ for use in our final galaxy selections, described in Section 6.

As discussed in Section 2.2, the signal we are looking for occurs on arcminute scales, comparable to the resolution of the SPT-SZ data we are working with. On the other hand, the overall anisotropy of the CMB is dominated by the primary signal, which is strongest on degree scales. For this reason it is essential for us to filter our maps before obtaining our measurements. Since we are making measurements on the smallest scales (approaching the beam size), we apply a filter to the SPT-SZ data in order to optimize point-source measurements. This optimal filter in Fourier-space, ψ, is (Schaffer et al. 2011)

$$\begin{eqnarray}&&\psi =\displaystyle \frac{\tau }{P}{\left[\int {d}^{2}k\displaystyle \frac{{\tau }^{2}}{P}\right]}^{-1},\end{eqnarray} \tag{ 19 }$$

where τ is the Fourier-space source profile and P is the Fourier-space noise power spectrum, which is the sum of the (squared) instrument-plus-atmosphere power spectral density and the primary CMB power spectrum. For a point source $\tau =B\times F$, where B is the Fourier-space beam function and F is the Fourier-space filter transfer function. We then scale ψ in order to preserve the total flux within a 1 arcmin radius circle in each map. Thus we expect our primary signal, which we measure in a 1 arcmin radius aperture around our stacked galaxies, to be minimally affected by our filtering. The resulting optimal point source filters for the 150 and 220 GHz bands are shown in Figure 8.

Zoom In Zoom Out Reset image size

We carried out our final co-add measurements by averaging the SPT-SZ maps around the galaxies in both our final low- and high-redshift galaxy samples. Before we are able to measure a signal from these averages, however, we first need to correct for a bias introduced by our removal of all sources within 4 arcmin of contaminating sources. Because the SPT-SZ maps themselves are normalized to a mean of 0, and all of the contaminating sources introduce positive signal into the maps, the average value in the uncontaminated regions of the maps is slightly biased to negative values. We therefore calculate a bias for the "contaminant-free" images by choosing 140,000 random points in our field (chosen so that there are not more random points than possible beams on the sky) and subjecting the points to the same contaminating-source cuts as our galaxies. We then take the resulting 107,561 random points and compute the mean sums within a 1 arcmin radius around each point. With these values we calculate an offset value needed to re-normalize the mean to 0. These offset values are 0.24 ± 0.09 and 0.58 ± 0.13 μK arcmin² at 150 and 220 GHz,respectively.

We then sum and average the total signal within 0.5, 1, 1.5, and 2 arcmin radius apertures around our sources and add the offset, scaling them appropriately for the different aperture sizes. The 0.5 arcmin radius aperture represents roughly the size of the 150 GHz and 220 GHz beam FWHMs, which are 1.15 arcmin and 1.05 arcmin, respectively. Additionally, we calculate the standard deviation for each of these measurements by finding the standard deviation of the same size co-added region around an equal number of random points in our field, subjected to the same contaminating source cuts. The offset uncertainties are also included but are negligible. The final co-add values for each aperture size and redshift range are given in Table 4. The final galaxy co-add images for both redshift ranges are shown in Figure 9.

Zoom In Zoom Out Reset image size

The upper left panel of this figure shows a clear $\approx 1$ arcmin size $\approx 2\sigma$ negative feature centered directly on our stack of low-z galaxies, with a magnitude consistent with a significant tSZ signal. Moreover, the low-z 220 GHz measurements show a strong positive signal centered on our co-added sources. Because the tSZ effect has a negligible impact at this frequency, this indicates that despite our cuts on detected contaminating sources, there still remains a significant positive contaminating signal at 220 GHz, made up of the sum of fainter sources. Looking at the emission by the typical range of dust temperatures at z = 1 (light blue and dark blue curves in Figure 10), the CMB spectrum (green curve in Figure 10), and the frequencies of our SPT measurements (rightmost red hatched regions in Figure 10), it appears likely that this contaminating signal at 220 GHz extends into the 150 GHz band, meaning that the true tSZ signal for the lower-redshift galaxies we have selected is even more negative than the values in Table 4.

Zoom In Zoom Out Reset image size

Moving to the higher redshift stacks shown in the bottom panel of Figure 9, we find that the stacked emission of our galaxies in the 150 GHz map is consistent with zero signal. However, as this band measures the sum of the (negative) tSZ and the (positive) contaminating sources, it is difficult to interpret these results without also considering the high-z measurements at 220 GHz. As in the lower redshift case, this band shows a clear excess, but now at a magnitude that is roughly twice that seen in the low-z stack. This suggests that, because it is more difficult to identify contaminating sources at higher redshift, the high-z 150 GHz measurement is more contaminated than the lower-redshift measurement, covering up the negative signal in which we are interested. In both redshift ranges, however, it is clear that obtaining the best possible constraints on AGN feedback requires making the best possible separation between the tSZ signal and the contaminating signal, a topic we address in detail below. Finally, we can convert our co-added ${\rm{\Delta }}T$ signal into gas thermal energy using Equation (16). These values (using a 1 arcmin radius aperture) are shown in Table 6, under "Data only."

As mentioned in Section 6, we performed the same co-adding method while also removing galaxies near bright radio sources. After these additional cuts, the number of final galaxies becomes 2219 for our low-z subset and 614 for our high-z subset. The resulting co-add values for a 1 arcmin radius aperture are: −1.7 ± 0.9 μK arcmin² for 150 GHz low-z, 2.4 ± 1.7 μK arcmin² for 220 GHz low-z, 0.9 ± 1.7 μK arcmin² for 150 GHz high-z, and 10.4 ± 3.3 μK arcmin² for 220 GHz high-z. We find that the additional radio source cuts do not significantly change our results except for an increased positive signal at 220 GHz in the high-z subset, though they do increase the noise in our measurements because we are co-adding a smaller number of galaxies. As a result, we do not use the radio source cuts in our modeling and analysis below, though the higher 220 GHz signal in the high-z subset may suggest a more significant tSZ detection in our high-z results.

In order to constrain the impact that undetected contaminating sources have on our tSZ measurements, we built a detailed model of contaminants based on an extrapolation of the SPT source counts measured in Mocanu et al. (2013). Our approach is to extend these counts to fainter values by modeling a random population of undetected sources that follow the trend of the detected sources into the unresolved region, which we then relate to the contaminating signal in both our 150 and 220 GHz measurements.

Following Mocanu et al. (2013) we separate contaminants into synchrotron sources, which emit most at lower frequencies, and dusty sources, which emit most at higher frequencies. For each source population we model the number counts as a power law,

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dS}}=\displaystyle \frac{{N}_{0}}{{S}_{{\rm{max}}}}{\left(\displaystyle \frac{S}{{S}_{{\rm{max}}}}\right)}^{\alpha },\end{eqnarray} \tag{ 20 }$$

where dN/dS is the number of sources between flux S and S + dS, N₀ is an overall amplitude, α is the power-law slope, and ${S}_{{\rm{max}}}$ is the flux at which all brighter sources have a 100% completeness level in the source count catalog. We then compute a range of allowed source count slopes from the Mocanu et al. (2013) data, by carrying out a ${\chi }^{2}$ fit in log-space. Our best-fit slopes at 150 GHz were ${\alpha }_{s}=-2.05\pm 0.04$ for the synchrotron sources and ${\alpha }_{d}=-2.70\pm 0.19$ for the dusty sources.

Note that our calculated values for ${\alpha }_{d}$ are much steeper than ${\alpha }_{s}$, meaning that while the number density of detected sources is dominated by synchrotron sources, the number density of undetected sources is likely to be dominated by dusty emitters. Note also that ${\alpha }_{d}$ and ${\alpha }_{s}$ are sufficiently steep that the number of sources diverges as S goes to 0, meaning that the source count distribution must fall off below some as-yet undetected flux. For simplicity, we model this fall-off as a minimum flux ${S}_{{\rm{min}}}$ below which there are no contaminating sources associated with the galaxies we are stacking.

For any choice of ${\alpha }_{d},$ ${\alpha }_{s},$ and ${S}_{{\rm{min}}}$ (which we will call a "source-count model"), we are then able to construct a model population of contaminating source fluxes through a four-step procedure as follows: (i) for each model source, we randomly decide whether it is a synchrotron source or a dusty source, such that the overall fraction of detectable dusty sources to synchrotron sources matches the observed source counts; (ii) we then assign the source a random 150 GHz flux, ${S}_{{\rm{150,rand}}},$ by inverting

$$\begin{eqnarray}&&{\displaystyle \int }_{{S}_{{\rm{150,min}}}}^{{S}_{{\rm{150,rand}}}}{dS}\;\displaystyle \frac{{dN}}{{dS}}=R{\displaystyle \int }_{{S}_{{\rm{150,min}}}}^{{S}_{{\rm{150,max}}}}{dS}\;\displaystyle \frac{{dN}}{{dS}},\end{eqnarray} \tag{ 21 }$$

where $R\in [0,1]$ is a random number, such that their overall population matches the source count slopes. This gives

$$\begin{eqnarray}&&{S}_{{\rm{150,rand}}}={[(1-R){S}_{{\rm{150,min}}}^{\alpha +1}+R{S}_{{\rm{150,max}}}^{\alpha +1}]}^{\displaystyle \frac{1}{\alpha +1}};\end{eqnarray} \tag{ 22 }$$

(iii) to obtain a corresponding flux for the source at 220 GHz we use the ${\alpha }_{220}^{150}$ spectral index distributions from Mocanu et al. (2013), which we assume to have normalized Gaussian shapes with the properties (center, σ) = (−0.55, 0.55) for synchrotron sources and (3.2, 0.89) for dusty sources. We then randomly choose ${\alpha }_{220}^{150}$ values that fit these distributions and calculate the 220 GHz flux (following Mocanu et al. 2013) as

$$\begin{eqnarray}&&{S}_{220,{\rm{rand}}}={S}_{150,{\rm{rand}}}\times 0.82\times {1.43}^{{\alpha }_{220}^{150}},\end{eqnarray} \tag{ 23 }$$

where we use units of μK arcmin² for all S; and (iv) finally, if the source has a detectable 150 or 220 GHz flux, we randomly discard it with a probability chosen to match the completeness percentages in the source count catalog.

For any single source-count model, we repeat the process 100,000 times, resulting in a large catalog of contaminating fluxes in both bands. From these, we compute the mean flux per contaminating source in each band, $\langle {S}_{150,{\rm{cont}}}\rangle$ and $\langle {S}_{220,{\rm{cont}}}\rangle$, which represents the contamination we are measuring in our stacks. To account for variations in the input parameters, we compute model contamination signals for a wide range of source-count models, with ${S}_{150,{\rm{max}}}=260\;\mu {\rm{K}}$ arcmin². We vary both ${\alpha }_{d}$ and ${\alpha }_{s}$ from $-2\sigma$ to $+2\sigma$ in steps of σ, and we let ${\mathrm{log}}_{10}$(${S}_{{\rm{min}}}$) vary from ${\mathrm{log}}_{10}$(0.01 μK arcmin²) to ${\mathrm{log}}_{10}$(30 μK arcmin²) in steps of 0.2 in log-space.

For each source-count model, we compute best-fit tSZ values by varying our two free parameters, tSZ signal (${S}_{{\rm{SZ}}}$) and the fraction of our measured sources that are contaminated (${f}_{{\rm{cont}}}$). We vary the tSZ signal from −50 to 50 μK arcmin² in steps of 0.1 μK arcmin², and we vary the fraction contaminated from −3 to 9 in steps of 0.01. For every combination of these parameters we compute

$$\begin{eqnarray}&&{\chi }^{2}({f}_{{\rm{cont}}},{S}_{{\rm{SZ}}})={ \mathcal B }\times {{ \mathcal A }}^{-1}\times {{ \mathcal B }}^{T},\end{eqnarray} \tag{ 24 }$$

where ${ \mathcal B }$ is the signal array,

$$\begin{eqnarray}&&{ \mathcal B }=\left(\begin{array}{c}{f}_{{\rm{cont}}}\times \langle {S}_{150,{\rm{cont}}}\rangle +{S}_{{\rm{SZ}}}-{S}_{150}\\ {f}_{{\rm{cont}}}\times \langle {S}_{220,{\rm{cont}}}\rangle -{S}_{220}\end{array}\right),\end{eqnarray} \tag{ 25 }$$

and ${ \mathcal A }$ is the noise matrix containing the noise for each band plus the covariance terms between each band,

$$\begin{eqnarray}{ \mathcal A }=\left(\begin{array}{cc}{\sigma }_{150}^{2} & {\sigma }_{150}{\sigma }_{220}\\ {\sigma }_{150}{\sigma }_{220} & {\sigma }_{220}^{2}\end{array}\right).\end{eqnarray} \tag{ 26 }$$

Here, S₁₅₀, S₂₂₀, ${\sigma }_{150}$, and ${\sigma }_{220}$ are our measured 1 arcmin radius values from Table 4. As discussed in Section 7.2, the σ values are computed using random point measurements. To be explicit here,

$$\begin{eqnarray}&&{\sigma }_{i}{\sigma }_{j}=\displaystyle \frac{{\displaystyle \sum }_{a=0}^{{N}_{{\rm{rand}}}}({S}_{i,{\rm{rand}}}^{a}-\langle {S}_{i,{\rm{rand}}}\rangle )\times ({S}_{j,{\rm{rand}}}^{a}-\langle {S}_{j,{\rm{rand}}}\rangle )}{{N}_{{\rm{rand}}}{N}_{{\rm{source}}}},\end{eqnarray} \tag{ 27 }$$

where i and j represent the bands, ${S}_{a,{\rm{rand}}}$ and ${S}_{b,{\rm{rand}}}$ represent the 1 arcmin radius aperture values for the random points, ${N}_{{\rm{rand}}}$ = 107,561 is the number of random points used, and ${N}_{{\rm{source}}}$ is the number of galaxies used (3394 for low-z and 924 for high-z). We then convert the ${\chi }^{2}$ values to Gaussian probabilities P by taking

$$\begin{eqnarray}&&P({S}_{{\rm{SZ}}})=\displaystyle \frac{{\displaystyle \sum }_{{f}_{{\rm{cont}}}\in [0,1]}\mathrm{exp}[-{\chi }^{2}({f}_{{\rm{cont}}},{S}_{{\rm{SZ}}})/2]}{{\displaystyle \sum }_{{f}_{{\rm{cont}}}\in [-3,9]}{\displaystyle \sum }_{{S}_{{\rm{SZ}}}\in [-50,50]\mu {\rm{K}}{\mathrm{arcmin}}^{2}}\mathrm{exp}[-{\chi }^{2}({f}_{{\rm{cont}}},{S}_{{\rm{SZ}}})/2]}.\end{eqnarray} \tag{ 28 }$$

Our approach is thus to marginalize over values of ${f}_{{\rm{cont}}}$ in the full physical range from 0 to $1,$ but normalize the overall probability by the sum of ${f}_{{\rm{cont}}}$ over a much larger range, including unphysical values. This excludes models in which a good fit to the data can only be achieved by moving ${f}_{{\rm{cont}}}$ outside the range of physically possible values.

Equation (28) then gives us a function $P({S}_{{\rm{SZ}}})$ for each combination of ${\alpha }_{d},$ ${\alpha }_{s},$ and ${S}_{{\rm{min}}}$. We can convert the corresponding ${S}_{{\rm{SZ}}}$ value to the gas thermal energy, ${E}_{{\rm{therm}}},$ using Equation (16) and the average ${l}_{{\rm{ang}}}^{2}$ from Table 3. Note that a positive detection of the tSZ effect is seen as a negative ${\rm{\Delta }}T$ signal at 150 GHz, and it represents a positive injection of thermal energy into the gas around the galaxy. Additionally, we compute a corresponding range for ${E}_{{\rm{grav}}}$ using Equation (12) and values from Table 3. The peak of each $P({S}_{{\rm{SZ}}})$ curve is shown as the colored points in Figure 11, where ${\alpha }_{s}$ (represented by point size) is increasing (i.e., becoming more positive) downward, and ${\alpha }_{d}$ (represented by point color) is increasing upward. The 1σ and 2σ contours are computed for each ${S}_{{\rm{min}}}$ by averaging $P({S}_{{\rm{SZ}}})$ across ${\alpha }_{d}$ and ${\alpha }_{s}.$ The resulting probability distribution depends only on ${E}_{{\rm{therm}}}$ and ${S}_{{\rm{min}}}$, and 1σ and 2σ are represented by the values $P({S}_{{\rm{SZ}}})=0.61$ and 0.13, respectively (i.e., $\mathrm{exp}[-{\sigma }^{2}/2]$). These contours are shown in Figure 11, along with the $\pm 1\sigma$ range for ${E}_{{\rm{grav}}}$. From this figure we see that there is a $\gt 2\sigma$ tSZ detection for every source model at low-z, and a $\approx 1\sigma$ detection of a signal exceeding the range that can be explained without feedback. At high-z, where the contaminants are harder to constrain, there is a $\approx 1\sigma$ tSZ detection for every source model.

Zoom In Zoom Out Reset image size

Finally, we average the probability distribution across ${S}_{{\rm{min}}}$ to get a final distribution as a function of only ${E}_{{\rm{therm}}}$. The significance values of this curve are shown in the "SPT only" part of Table 6. We see a $\gt 3\sigma$ total tSZ detection at low-z, and a nearly $1\sigma$ tSZ detection at high-z. Furthermore, in both redshift bins, the best fit values are higher than expected from models that do not include AGN energy input.

In order to better constrain the impact of dusty contaminating sources on our measurements, we make use of the 2015 public data release from the Planck mission, focusing on the high-frequency bands at 857, 545, 353, and 217 GHz (see the rightmost black hatched regions in Figure 10). While the FWHM beam size for this data is about 5 arcmin (Planck Collaboration et al. 2015a), and thus too low-resolution to detect the tSZ signal in which we are interested, the data provides information at the higher frequencies at which the dusty sources should be much brighter (i.e., the light blue and dark blue curves in Figure 10). Therefore these measurements have the potential to discriminate between contaminant models, allowing us to better remove this contribution from the tSZ signal.

Our goal is to use this data to add terms to our ${\chi }^{2}$ fit that quantify, for each source model, how consistent a given choice of ${f}_{{\rm{cont}}}$ is with the Planck measurements. To compute these extra ${\chi }^{2}$ terms, we again construct stacks of the data over each of our galaxies, but now, because of the lower resolution of the Planck data, we extend our contaminant source cuts to within 10 arcmin of known potentially contaminating sources. This results in a decreased number of final galaxies, now 937 at low-z and 240 at high-z. In order to filter out the primary CMB signal, we convolve each map with a 7 arcmin FWHM Gaussian and subtract the resulting map from the original. We then stack the central pixels of every source to get co-added values for our galaxies in each of the Planck bands. In addition, we degrade the SPT 150 and 220 GHz maps to match the beam size of Planck, apply the same 7 arcmin FWHM filtering, and stack the galaxies on those images as well.

As was the case in Section 7.2, in all of these stacks there is an offset we need to correct for since we are purposely avoiding positive contaminations in the maps. To do this we also make measurements at 3000 random points on the sky that were restricted to the same contaminating-source cuts as our galaxies. These measurements allow us to compute offset values needed to re-normalize each band to a mean of 0, which we applied to our final measurements.

Finally, we compute our measurement errors by using the random point measurements (since the proper noise covariance matrix is not provided, i.e., Planck Collaboration et al. 2015b), corrected in two ways. First, because we account for the residual CMB primary signal in the ${\chi }^{2}$ calculations as discussed below, we remove the error due to the CMB primary itself. To estimate this contribution, we take 95% of the minimum covariance between the SPT 150 and 220 GHz bands (filtered to match the Planck bands) and the Planck 217 GHz band, since these are mostly dominated by the CMB primary signal, which will therefore be correlated between them. The minimum covariance is between the two SPT bands, and it is 7.85 μK. Second, there is an error introduced due to our offset corrections because they are made from a large, but finite number of points. We then get the corrected error from

$$\begin{eqnarray}&&{({\sigma }_{i}{\sigma }_{j})}_{{\rm{corr}}}=\sqrt{\displaystyle \frac{{\sigma }_{i}{\sigma }_{j}-{\sigma }_{{\rm{cov}}}^{2}}{{N}_{{\rm{source}}}}+\displaystyle \frac{{\sigma }_{i}{\sigma }_{j}}{{N}_{{\rm{random}}}}},\end{eqnarray} \tag{ 29 }$$

where ${\sigma }_{i}{\sigma }_{j}$ is given by Equation (27), with i and j representing the various bands used, ${\sigma }_{{\rm{cov}}}=7.85\;\mu {\rm{K}}$ is the minimum CMB covariance discussed above, ${N}_{{\rm{source}}}$ is the number of sources used for the measurements (937 for low-z and 240 for high-z), and ${N}_{{\rm{random}}}=3000$ is the number of random points used.

Note that this represents both the error due to detector noise in each band, as well as the error due to contributions from foregrounds on the sky. In fact, the majority of the variance at the highest frequencies is correlated between the bands and likely due to contributions from Galactic dust emission. However, unlike the primary CMB signal, the spectral shape of this foreground is similar to that of the dusty sources we are trying to constrain, and it cannot be removed by fitting it separately.

In the same manner as the previous section, we model SPT-SZ 150 and 220 GHz contaminant source fluxes using a range of different source-count models (i.e., ${\alpha }_{d},$ ${\alpha }_{s},$ and ${S}_{{\rm{min}}}$), resulting again in 100,000 modeled contaminating source fluxes in each SPT band, ${S}_{150,{\rm{cont}}}$ and ${S}_{220,{\rm{cont}}}$. We also model what the contaminating signal would be in the Planck bands and the SPT bands filtered to match Planck. For each modeled contaminating source, if it is chosen to be a synchrotron source we simply extrapolate the Planck-based fluxes as

$$\begin{eqnarray}&&{S}_{\nu ,{\rm{sync}}}={S}_{150,{\rm{cont}}}\times {\left(\displaystyle \frac{\nu }{150}\right)}^{{\alpha }_{220}^{150}}\times {C}_{\nu }\times F,\end{eqnarray} \tag{ 30 }$$

where ${\alpha }_{220}^{150}$ is the same used in the previous section, C_ν is a frequency-dependent factor involved in the conversion from Jy/sr to μK, and $F=0.021$ is the factor required to preserve the signal within a 1 arcmin radius aperture after applying the Planck filtering we used.

In order to accurately describe thermal dust emission across the Planck frequencies, we adopt a modified blackbody with a free emissivity index, β, and dust temperature, ${T}_{{\rm{dust}}},$ often referred to as a graybody (Planck Collaboration et al. 2015c). This requires us to add another free parameter, the temperature of the contaminant dust, ${T}_{{\rm{dust}}}$. This slope of each dusty source as a function of frequency is then

$$\begin{eqnarray}&&{\left.\displaystyle \frac{d\mathrm{ln}{S}_{\nu }}{d\mathrm{ln}\nu }\right|}_{\nu =185{\rm{GHz}}}=3+\beta -{x}_{185}{[1-\mathrm{exp}(-{x}_{185})]}^{-1},\end{eqnarray} \tag{ 31 }$$

where ${x}_{185}\equiv (185\;\mathrm{GHz})\times h/({kT})=(185/416)(1+z)/{T}_{20}$ and T₂₀ is the dust temperature in units of 20 K, and we use the slope of the blackbody function at $\nu =185\;{\rm{GHz}}$ because it is halfway between our two SPT bands (150 and 220 GHz). This can be related, in turn, to the power law index ${\alpha }_{220}^{150},$ from Section 7.3, as

$$\begin{eqnarray}&&\beta +3={\alpha }_{220}^{150}+{x}_{185}{[1-\mathrm{exp}(-{x}_{185})]}^{-1}.\end{eqnarray} \tag{ 32 }$$

This then gives

$$\begin{eqnarray}{S}_{\nu ,{\rm{dust}}} & = & \;{S}_{150,{\rm{cont}}}\times {\left(\displaystyle \frac{\nu }{150}\right)}^{{\alpha }_{220}^{150}+{x}_{185}{[1-\mathrm{exp}(-{x}_{185})]}^{-1}}\\ & & \times \displaystyle \frac{\mathrm{exp}[(150/416)(1+z)/{T}_{20}]-1}{\mathrm{exp}[(\nu /416)(1+z)/{T}_{20}]-1}\times {C}_{\nu }\times F,\end{eqnarray} \tag{ 33 }$$

where we vary ${T}_{{\rm{dust}}}$ from 20 to 50 K in steps of 3 K.

With these expressions, we are able to compute ${\chi }^{2}$ values for each source-count model accounting for the Planck measurements. This time, in addition to varying ${f}_{{\rm{cont}}}$ and ${S}_{{\rm{SZ}}}$, we also vary ${T}_{{\rm{dust}}}$ (as discussed above) and a parameter ${\rm{\Delta }},$ which represents the offset due to the CMB primary signal, which we vary from −3 to 3 μK in steps of 0.1 μK. Computing ${\chi }^{2}$ now involves the original SPT terms plus the new Planck terms, and it follows the same process as in the previous section (e.g., Equation (24)),

$$\begin{eqnarray}&&{\chi }^{2}({f}_{{\rm{cont}}},{S}_{{\rm{SZ}}},{T}_{{\rm{dust}}},{\rm{\Delta }})={ \mathcal B }\times {{ \mathcal A }}^{-1}\times {{ \mathcal B }}^{T},\end{eqnarray} \tag{ 34 }$$

where ${ \mathcal B }$ is the signal array and ${ \mathcal A }$ is the noise matrix containing the noise for each band plus the covariance terms between each band. We will denote each element of the signal array ${{ \mathcal B }}_{i}$, where i runs over the two SPT bands (i.e., 150 and 220 GHz) and then every Planck-filtered band (i.e., the Planck bands at 857, 545, 353, and 217 GHz, plus the SPT bands at 220 and 150 GHz filtered to match the Planck images). We then have ${{ \mathcal B }}_{1}={f}_{{\rm{cont}}}\times \langle {S}_{150,{\rm{cont}}}\rangle +{S}_{{\rm{SZ}}}-{S}_{150}$, ${{ \mathcal B }}_{2}={f}_{{\rm{cont}}}\times \langle {S}_{220,{\rm{cont}}}\rangle -{S}_{220}$, and ${{ \mathcal B }}_{3-8}={f}_{{\rm{cont}}}\ \times \langle {S}_{3-8,{\rm{cont}}}\rangle +{\rm{\Delta }}-{S}_{3-8}$. As before, S_i represents the final values of our galaxy stacks for each band. We similarly define the elements of the noise matrix as ${{ \mathcal A }}_{{ij}}={\sigma }_{i}{\sigma }_{j}$, where i and j run over all of the bands and ${\sigma }_{i}{\sigma }_{j}$ is given by Equation (29).

As in the previous section, we then convert the ${\chi }^{2}$ values to Gaussian probabilities by taking

$$\begin{eqnarray}P({S}_{{\rm{SZ}}}) & = & \displaystyle \sum _{{f}_{{\rm{cont}}}\in [0,1],{T}_{{\rm{dust}}},{\rm{\Delta }}}\displaystyle \frac{\mathrm{exp}[-{\chi }^{2}({f}_{{\rm{cont}}},{S}_{{\rm{SZ}}},{T}_{{\rm{dust}}},{\rm{\Delta }})/2]}{{\displaystyle \sum }_{{\rm{all}}}\mathrm{exp}[-{\chi }^{2}/2]},\end{eqnarray} \tag{ 35 }$$

where the whole function is normalized to a total of 1, and each final SZ value contains the sum over the corresponding ${T}_{{\rm{dust}}}$, Δ, and fractions from 0 to 1. Since in this case there are eight terms contributing to ${\chi }^{2}$ and four fit parameters, this leaves us with 4 degrees of freedom. Thus the minimum ${\chi }^{2}$ will not be 0 in every case as they were previously with just two parameters and two fit parameters, and so for each model we scale the final probabilities by $\mathrm{exp}(-{\chi }_{{\rm{min}}}^{2}/2)$, where ${\chi }_{{\rm{min}}}^{2}$ is the minimum ${\chi }^{2}$ value for that model.

This again gives us a function $P({S}_{{\rm{SZ}}})$ for each combination of ${\alpha }_{d},$ ${\alpha }_{s},$ and ${S}_{{\rm{min}}}$, which we can convert to an energy ${E}_{{\rm{therm}}}$. The peaks of each $P({S}_{{\rm{SZ}}})$ curve are shown as the colored points in Figure 12, where the points are colored by the minimum ${\chi }^{2}$ value for each model. Note that the the best fit ${\chi }^{2}$ values are smaller than expected for 4 degrees of freedom due to the correlations between the errors of the various frequency bands due primarily to foreground contamination by Galactic dust. The 1σ and 2σ contours are created for each ${S}_{{\rm{min}}}$ by averaging $P({S}_{{\rm{SZ}}})$ across ${\alpha }_{d}$ and ${\alpha }_{s}$ and then dividing the final result by the single maximum value. 1σ and 2σ are again represented by the values 0.61 and 0.13, respectively. These contours are shown in Figure 12, along with the $\pm 1\sigma$ range for ${E}_{{\rm{grav}}}$. From this figure we can see that the σ values have slightly decreased at both low-z and high-z, especially at higher ${S}_{{\rm{min}}}$ values for high-z. This is where the contaminants are hardest to constrain with just SPT, and where Planck data helps us the most.

Zoom In Zoom Out Reset image size

Finally, we average the probability distribution across ${S}_{{\rm{min}}}$, divided by the maximum value again, and get a final distribution as a function of only ${E}_{{\rm{therm}}}$. The significance (σ) values of this curve are shown in the "With Planck" part of Table 6. Planck has helped to constrain the tSZ signal, especially at high-z, although it is clear that the gain in sensitivity has been limited by the decrease in the number of galaxies in each redshift bin due to the much larger beam size of Planck compared to SPT.

Alternatively, we can also characterize the total tSZ signal for our coadds with the angularly integrated Compton-y parameter, Y (e.g., Ruan et al. 2015). While we cannot directly compare peak Compton-y values with past measurements, as these are beam-dependent quantities, we can compare the angularly integrated Y values between our results and past experiments (see Table 5). Using Equation (15) at 150 GHz, this is

$$\begin{eqnarray}Y & \equiv & {l}_{{\rm{ang}}}^{2}\displaystyle \int y({\boldsymbol{\theta }})d{\boldsymbol{\theta }}\\ & = & -3.5\times {10}^{-8}\;{{\rm{Mpc}}}^{2}{\left(\displaystyle \frac{{l}_{{\rm{ang}}}}{{\rm{Gpc}}}\right)}^{2}\displaystyle \frac{\displaystyle \int {\rm{\Delta }}{T}_{150}({\boldsymbol{\theta }})d{\boldsymbol{\theta }}}{\mu {\rm{K}}{\text{arcmin}}^{2}},\end{eqnarray} \tag{ 36 }$$

such that $Y=2.9\times {10}^{-8}\;{{\rm{Mpc}}}^{2}{E}_{{\rm{60}}}$, where ${E}_{{\rm{60}}}$ is ${E}_{{\rm{therm}}}$ in units of 10⁶⁰ erg. In these units, the mean Y values computed directly from the 150 GHz maps from our co-added galaxies are $1.2(\pm 0.6)\times {10}^{-7}$ Mpc² for low-z and $-1.7(\pm 1.5)\times {10}^{-7}$ Mpc² for high-$z.$ When these measurements are corrected for contamination using the 220 GHz SPT data, the mean Y values become ${2.3}_{-0.7}^{+0.9}\times {10}^{-7}$ Mpc² for low-z and ${1.9}_{-2.0}^{+2.4}\times {10}^{-7}$ Mpc² for high-$z,$ and once the Planck data is incorporated, the mean Y values become ${2.2}_{-0.7}^{+0.9}\times {10}^{-7}$ Mpc² for low-z and ${1.7}_{-1.8}^{+2.2}\times {10}^{-7}$ Mpc² for high-$z.$ These values are also given in Table 6.

Table 5.  Previous tSZ Measurements

Study	N	Type	z	Mass (M⊙)	tSZ Measurement	Type
Chatterjee et al. (2009)	500,000	SDSS quasars	0.08–2.82	⋯	(7.0 ± 3.4) $\times {10}^{-7}$	y
Chatterjee et al. (2009)	1,000,000	SDSS LRGs	0.4–0.6	⋯	(5.3 ± 2.5) $\times {10}^{-7}$	y
Hand et al. (2011)	1732	SDSS radio-quiet LRGs	0.30 (mean)	8.0 × 1013	(7.9 ± 6.2) × 10−7 Mpc2	${Y}_{200\bar{\rho }}$
Gralla et al. (2014)	667	SDSS radio-loud AGN	0.3 (median)	2 × 1013	(1.5 ± 0.5) × 10−7 Mpc2	Y200
Gralla et al. (2014)	4,352	FIRST AGN	1.06 (median)	⋯	(5.7 ± 1.3) × 10−8 Mpc2	Y200
Greco et al. (2014)	188,042	SDSS LBGs	0.05–0.3	1.4 × 1011 ⋆	$({0.6}_{-0.6}^{+5.4})$ × 10−6 arcmin2	${\tilde{Y}}_{c}^{\mathrm{cyl}}$
Ruan et al. (2015)	≈14,000	SDSS quasars	1.96 (median)	5.0 × 1012	(4.8 ± 0.8) × 10−6 Mpc2	Y
Ruan et al. (2015)	≈14,000	SDSS quasars	0.96 (median)	5.0 × 1012	(2.2 ± 0.9) × 10−6 Mpc2	Y
Ruan et al. (2015)	81,766	SDSS LBGs	0.54 (median)	3.2 × 1011 ⋆	(1.4 ± 0.4) × 10−6 Mpc2	Y
Crichton et al. (2015)	17,468	SDSS radio-quiet quasars	0.5–3.5	⋯	(6.2 ± 1.7) × 1060 erg	Eth

Note. LRGs—luminous red galaxies; LBGs—locally brightest galaxies. Masses refer to halo masses, except for those of Greco et al. (2014) and Ruan et al. (2015) LBGs which refer to stellar masses ($\star$). We select Hand et al. (2011) and Greco et al. (2014) values that have the most similar masses to our galaxies.

Download table as:  ASCIITypeset image

Table 6.  Our Final tSZ Measurements Using Various Methods for Removing Contamination

Model	N	z	$\int {\rm{\Delta }}{T}_{150}({\boldsymbol{\theta }})d{\boldsymbol{\theta }}$	Y	Etherm (±1σ)	Etherm (±2σ)	S/N
			(μK arcmin2)	(10−7 Mpc2)	(1060 erg)	(1060 erg)	(S/σ)
Data only	3394	0.5–1.0	−1.5 ± 0.7	1.2 ± 0.6	4.1 ± 1.9	4.1 ± 3.8	2.2
	924	1.0–1.5	1.6 ± 1.4	−1.7 ± 1.5	−5.8 ± 5.1	−5.8 ± 10.2	−1.1
χ2 (SPT only)	3394	0.5–1.0	$-{2.9}_{-1.1}^{+0.9}$	${2.3}_{-0.7}^{+0.9}$	${8.1}_{-2.5}^{+3.0}$	${8.1}_{-4.8}^{+6.8}$	3.5
	924	1.0–1.5	$-{1.8}_{-2.3}^{+1.9}$	${1.9}_{-2.0}^{+2.4}$	${6.7}_{-7.0}^{+8.3}$	${6.7}_{-13.3}^{+18.6}$	0.9
χ2 (With Planck)	937	0.5–1.0	$-{2.8}_{-1.1}^{+0.8}$	${2.2}_{-0.7}^{+0.9}$	${7.6}_{-2.3}^{+3.0}$	${7.6}_{-4.3}^{+7.1}$	3.6
	240	1.0–1.5	$-{1.7}_{-2.1}^{+1.7}$	${1.7}_{-1.8}^{+2.2}$	${6.0}_{-6.3}^{+7.7}$	${6.0}_{-12.3}^{+18.0}$	0.9

Note. The last three columns represent the best fit E_therm values with ±1σ values and ±2σ values and the E_therm signal-to-noise ratio (S/σ), respectively.

Download table as:  ASCIITypeset image

Our $0.5\leqslant z\leqslant 1.0$ value of ${2.2}_{-0.7}^{+0.9}\times {10}^{-7}$ Mpc² is more than three times smaller than the Hand et al. (2011) $z\approx 0.3$ SDSS radio-quiet LRG result. If we estimate the stellar mass for the Hand et al. (2011) results using Table 5 and Equation (11) we get $1.7\times {10}^{12}{M}_{\odot }$, which is about an order of magnitude greater than the average stellar mass of our galaxies. Our smaller values could be indicative of the relation that tSZ signal increases with halo (and stellar) mass (i.e., Gralla et al. 2014). Our low-z result is within about 1σ of the Gralla et al. (2014) $z\approx 0.3$ SDSS radio-loud AGN result, though our $1.0\lt z\leqslant 1.5$ result of ${1.7}_{-1.8}^{+2.2}\times {10}^{-7}$ Mpc² is $\gt 3\sigma$ higher than their $z\approx 1.1$ FIRST AGN result. This discrepancy is not too significant because our high-z detection is only at a 0.9σ confidence. Our results are also within $\approx 1\sigma$ of Greco et al. (2014) when comparing their results for galaxies with masses similar to ours. At smaller masses our results are consistent with theirs, while at larger masses they find an even greater tSZ signal, following the relation that tSZ signal increases with stellar mass. Ruan et al. (2015) also obtain values from stacks of SDSS quasars about an order of magnitude ($\gt 2\sigma$) higher than ours, although, according to Cen & Safarzadeh (2015b), the maximum AGN feedback signal from Ruan et al. (2015) can only be 25% of their quoted values. Furthermore, the $\approx {10}^{11}{M}_{\odot }$ galaxy results from Ruan et al. (2015) are consistent with zero signal, while their $\approx 3\times {10}^{11}{M}_{\odot }$ galaxy results are $\gt 2\sigma$ larger than ours. Their high mass sample represents almost three times the mean mass of our galaxies, though, so the larger values may be indicative of the stellar mass−tSZ signal relation as well as the potential overestimation of the tSZ signal. We can compare our results to Cen & Safarzadeh (2015b) by multiplying their average Compton-y values over 1 arcmin by $\pi \times {l}_{{\rm{ang}}}^{2}$ (where ${l}_{{\rm{ang}}}^{2}=3.18\;{{\rm{Gpc}}}^{2}$ for $z=1.5$) to get Y values of $\approx 6.8\times {10}^{-7}$ Mpc² for their halo occupation distribution (HOD) model, and $\approx 4.2\times {10}^{-7}$ Mpc² for their Cen & Safarzadeh (2015a, CS) model, which places quasars in lower mass dark matter halos. This CS model value is within ≈2σ of our results, and our measurements would favor the lower estimates of their CS model over their HOD model. Finally, both our low-z and high-z ${E}_{{\rm{therm}}}$ results are well within $1\sigma$ of the Crichton et al. (2015) SDSS radio-quiet quasar results.

With Equations (12) and (14) and the redshifts and masses from Table 3, we can also investigate theoretical thermal energies of the gas around elliptical galaxies due to both gravity and AGN feedback. We estimate the gravitational heating energy to be ${E}_{{\rm{therm,grav}}}={3.6}_{-1.9}^{+3.6}\times {10}^{60}$ erg for our low-z sample and ${E}_{{\rm{therm,grav}}}={3.0}_{-1.6}^{+3.0}\times {10}^{60}$ erg for our high-z sample. We therefore measure excess non-gravitational energies (for our results using Planck) of ${E}_{{\rm{therm,feed}}}={4.0}_{-4.3}^{+3.6}\times {10}^{60}$ erg for low-z and ${E}_{{\rm{therm,feed}}}={3.0}_{-7.0}^{+7.9}\times {10}^{60}$ erg for high-z. Plugging these into the theoretical AGN feedback energy equation and solving for ${\epsilon }_{k}$, we get feedback efficiencies of ${1.5}_{-1.6}^{+1.3}$% for low-z and ${1.3}_{-3.1}^{+3.5}$% for high-z. These values are very uncertain, though they do hint at feedback efficiencies slightly lower than the suggested 5% (i.e., Scannapieco et al. 2005; Ruan et al. 2015).

The measurements described above depend on co-adding data from a large number of sources. In principle, however, there is additional information in the distribution of measured values that is lost through this process. For example, imagine a set of 1001 150 GHz measurements, 1000 of which contributed a negative signal of −2 μK arcmin² and one of which contributed a positive signal of 2000 μK arcmin². While the average value of these measurements would be zero, looking at the distribution of values would indicate strong evidence of a negative tSZ signal, offset by contamination from a single, overpowering positive source.

To quantify the additional information available by the full distribution of SPT data, we apply the same contaminant source modeling described above (i.e., Section 7.3) and use a goodness-of-fit test, the A–D test (Anderson & Darling 1954), to find models that poorly fit the data. In this case we restrict our attention purely to the 1 arcmin resolution data used to construct Figure 9, without folding in the lower resolution Planck data as described in Section 7.4. To perform the test, we run through every pair of galaxy measurements (i.e., the 1 arcmin radius aperture sums at both 150 and 220 GHz) and find the fraction of galaxy measurements in each of the four quadrants around the pair of measurements (i.e., ($x\lt {x}_{i},y\lt {y}_{i}$); ($x\lt {x}_{i},y\gt {y}_{i}$); ($x\gt {x}_{i},y\lt {y}_{i}$); ($x\gt {x}_{i},y\gt {y}_{i}$); where x and y represent the co-add sum in each band, and i runs through all the galaxies). We will call these fractions ${f}_{i,j}$, where i specifies the galaxy ($i=[1,2,...,n]$, with n being the number of galaxies) and j specifies the quadrant ($j=[1,2,3,4]$). In the same four quadrants we also find the fraction of model measurements, ${F}_{i,j}$. We therefore define our A–D statistic as

$$\begin{eqnarray}&&{S}_{\mathrm{AD}}=n\;\sum \displaystyle \frac{{({f}_{i,j}-{F}_{i,j})}^{2}}{{F}_{i,j}(1-{F}_{i,j})},\end{eqnarray} \tag{ 37 }$$

where a smaller S_AD indicates a better fit between the model and the data.

Each model, as a function of ${S}_{{\rm{min}}}$, ${\alpha }_{{\rm{dust}}}$, and ${\alpha }_{{\rm{sync}}}$, will have a corresponding S_AD. In order to interpret our results, we follow the same process, but instead of using data from our selected galaxies we use a random subset of the modeled sources, with the subset containing the same number of elements as the number of galaxies we are using for both redshift ranges. We do this 200 times each for 4 different combinations of ${S}_{{\rm{min}}}$, ${\alpha }_{{\rm{dust}}}$, and ${\alpha }_{{\rm{sync}}}$. We then define our confidence of the model fits, $P(\lt {S}_{\mathrm{AD}})$, as the fraction of these subset calculations that are less than the corresponding S_AD. When defined in this way, P indicates our confidence that the model is not a good fit with the data. The fractions correspond to σ in the standard way (i.e., 1σ = 0.68, 2σ = 0.95, etc.).

As mentioned above, we only do these subset calculations for four sets of model parameters, which is due to their time-intensive nature. To get confidence values for every other set of model parameters, we simply do a linear interpolation between the four sets of model parameters that we did use. The results are shown in Figure 13. They reveal that the models with both the highest ${S}_{{\rm{min}}}$ and the highest ${\alpha }_{{\rm{dust}}}$ are disfavored up to $\approx 1.5\sigma$ (87%) confidence, with the trend much more pronounced in the high-z range. On the other hand, unlike the analysis using the Planck data, the A–D test primarily serves to constrain ${S}_{{\rm{min}}}$ rather than ${E}_{{\rm{thermal}}},$ meaning that it does not allow us to obtain significantly better constraints on feedback itself. Finally, we note that we also carried out a two-dimensional Kolmogorov–Smirnov test (e.g., Peacock 1983) but it was far less constraining than the A–D test, and so we do not present it here.

Zoom In Zoom Out Reset image size