SEARCHING FOR FOSSIL EVIDENCE OF AGN FEEDBACK IN WISE-SELECTED STRIPE-82 GALAXIES BY MEASURING THE THERMAL SUNYAEV–ZEL'DOVICH EFFECT WITH THE ATACAMA COSMOLOGY TELESCOPE

Galaxy formation was long expected to proceed hierarchically, with larger galaxies forming at later times when larger dark matter halos coalesce and gas has longer to cool and condense (e.g., Rees & Ostriker 1977; White & Frenk 1991; Richstone et al. 1998; Cattaneo et al. 1999; Kauffmann & Haehnelt 2000; Menci 2006). However, an increasing amount of observational evidence suggests recent anti-hierarchical trends in both galaxy and active galactic nucleus (AGN) evolution (Kang et al. 2016; Rosas-Guevara et al. 2016; Siudek et al. 2016). For example, beyond $z\approx 2$ the typical mass of star-forming galaxies has dropped by more than a factor of ≈3 (Cowie et al. 1996; Brinchmann et al. 2004; Kodama et al. 2004; Bauer et al. 2005; Bundy et al. 2005; Feulner et al. 2005; Treu et al. 2005; Papovich et al. 2006; Noeske et al. 2007; Cowie & Barger 2008; Drory & Alvarez 2008; Vergani et al. 2008), while the characteristic AGN luminosity has dropped by more than a factor of ≈10 (Pei 1995; Ueda et al. 2003; Barger et al. 2005; Buchner et al. 2015). To explain these trends, it is widely suggested that galaxy evolution models require additional heating of the circumgalactic medium by energetic AGN feedback (Merloni 2004; Scannapieco & Oh 2004; Scannapieco et al. 2005; Bower et al. 2006; Neistein et al. 2006; Thacker et al. 2006; Sijacki et al. 2007; Merloni & Heinz 2008; Chen et al. 2009; Hirschmann et al. 2012, 2014; Mocz et al. 2013; Lapi et al. 2014; Schaye et al. 2015; Keller et al. 2016). This typically involves an energetic AGN outflow caused by radiation pressure or jets that blow material out of the galaxy, heating the nearby intergalactic medium enough to suppress further generations of stars and AGNs.

In fact, there is abundant observational evidence of AGN feedback in action in galaxy clusters (Schawinski et al. 2007; Rafferty et al. 2008; Fabian 2012; Farrah et al. 2012; Page et al. 2012; Teimoorinia et al. 2016). Most notably, the centers of clusters are more likely to contain galaxies that host large radio-loud jets of AGN-driven material (Burns 1990; Best et al. 2005; McNamara et al. 2005), whose energies are comparable to those needed to stop the gas from cooling (e.g., Simionescu et al. 2009). Furthermore, AGN feedback from the central cD galaxies in clusters increases in proportion to the cooling luminosity, as expected in an operational feedback loop (e.g., Bîrzan et al. 2004; Rafferty et al. 2006).

Direct measurements of AGN feedback in less dense environments are much less common, primarily because of the relatively high redshifts and faint signals involved, although evidence of AGN feedback has also been seen in relatively nearby galaxies (Tombesi et al. 2015; Lanz et al. 2016; Schlegel et al. 2016). For example, broad absorption-line outflows are observed as blueshifted troughs in the rest-frame spectra of ≈20% of all quasars (Hewett & Foltz 2003; Ganguly & Brotherton 2008; Knigge et al. 2008), but quantifying AGN feedback requires estimating mass-flow and the energy released by these outflows (e.g., Wampler et al. 1995; de Kool et al. 2001; Hamann et al. 2001; Feruglio et al. 2010; Sturm et al. 2011; Veilleux et al. 2013). These quantities can only be computed in cases for which it is possible to estimate the distance to the outflowing material from the central source, which is often highly uncertain. While these measurements have been carried out for a select set of objects (e.g., Chartas et al. 2007; Moe et al. 2009; Dunn et al. 2010; Greene et al. 2012; Borguet et al. 2013; Chamberlain et al. 2015), it is still unclear how these results generalize to AGNs as a whole. Furthermore, it is still an open question as to whether AGN outflows triggered by galaxy interactions in massive, high-redshift galaxies actually quench star formation (e.g., Fontanot et al. 2009; Pipino et al. 2009; Debuhr et al. 2010; Ostriker et al. 2010; Faucher-Giguère & Quataert 2012; Newton & Kay 2013; Feldmann & Mayer 2015; Bongiorno et al. 2016).

One way to get around the high redshifts and faint signals involved in AGN feedback measurements is to stack measurements of the cosmic microwave background (CMB) radiation. At angular scales smaller than ≈5 arcmin, Silk damping washes out the primary CMB anisotropies (Silk 1968; Planck Collaboration et al. 2015d), leaving room for secondary anisotropies such as the Sunyaev–Zel'dovich (SZ) effect, where CMB photons interact with hot, ionized gas (Sunyaev & Zeldovich 1970, 1972). When the gas has a bulk velocity, CMB photons interacting with electrons in the gas experience a Doppler boost, resulting in frequency-independent fluctuations in the CMB temperature, known as the kinematic Sunyaev–Zel'dovich (kSZ) effect. Although the kSZ effect does not measure the energy of the gas, it can be used to detect the ionized gas profile within dark matter halos, thereby providing information on where hot gas is located around galaxies. This can be useful for understanding how AGN feedback heats up gas and moves it around (Battaglia et al. 2010). Several recent studies have made significant measurements of the kSZ effect in galaxy clusters by stacking CMB observations (e.g., Planck Collaboration et al. 2016; Schaan et al. 2016; Soergel et al. 2016).

If the gas is sufficiently heated, inverse Compton scattering coupled with the thermal motions of electrons will shift the CMB photons to higher energies. This thermal Sunyaev–Zel'dovich (tSZ) effect directly depends on the temperature of the free electrons that the CMB radiation passes through, and it has a unique redshift-independent spectral signature that makes it well suited for measuring the gas heated through AGN feedback (Voit 1994; Birkinshaw 1999; Natarajan & Sigurdsson 1999; Platania et al. 2002; Lapi et al. 2003; Chatterjee & Kosowsky 2007; Chatterjee et al. 2008; Scannapieco et al. 2008; Battaglia et al. 2010). Individual tSZ distortions per source are very small, however, so a stacking analysis must be performed on many measurements in order to derive a significant signal.

This method has been used previously by a handful of studies in relation to AGNs and galaxies. Chatterjee et al. (2010) found a tentative detection of quasar feedback using the Sloan Digital Sky Survey (SDSS) and Wilkinson Microwave Anisotropy Probe, though the significance of AGN feedback in their measurements is disputed (Ruan et al. 2015). Hand et al. (2011) stacked >2300 SDSS-selected "luminous red galaxies" in data from the Atacama Cosmology Telescope (ACT) and found a $2.1\sigma \mbox{--}3.8\sigma$ tSZ detection after selecting radio-quiet galaxies and binning them by luminosity. Planck Collaboration et al. (2013) investigated the relationship between tSZ signal and stellar mass using ≈260,000 "locally brightest galaxies" with significant results, especially with stellar masses $\gtrsim {10}^{11}\,{M}_{\odot }$. Gralla et al. (2014) stacked data from the ACT at the positions of a large sample of radio AGNs selected at 1.4 GHz to make a 5σ detection of the tSZ effect associated with the haloes that host active AGNs. Greco et al. (2015) used Planck full-mission temperature maps to examine the stacked tSZ signal of 188,042 "locally brightest galaxies" selected from the SDSS Data Release 7, finding a significant measurement of the stacked tSZ signal from galaxies with stellar masses above $\approx 2\times {10}^{11}\,{M}_{\odot }$. Ruan et al. (2015) stacked Planck tSZ Compton-y maps centered on the locations of 26,686 spectroscopic quasars identified from SDSS to estimate the mean thermal energies in gas surrounding such $z\approx 1.5$ quasars to be $\approx {10}^{62}$ erg, although the significance of AGN feedback in their measurements has also been disputed (Cen & Safarzadeh 2015). Crichton et al. (2016) stacked >17,000 radio-quiet quasars from the SDSS in ACT data and found $3\sigma$ evidence for the presence of associated thermalized gas and $4\sigma$ evidence for the thermal coupling of quasars to their surrounding medium. Spacek et al. (2016) stacked 937 massive elliptical galaxies using the South Pole Telescope (SPT) and made a $3.6\sigma$ detection of the tSZ effect at a magnitude suggesting an excess of non-gravitational thermal energy, possibly due to AGN feedback. These tSZ AGN feedback measurements continue to be promising, and in this paper we will especially focus on the methods and results from Spacek et al. (2016).

As evidenced above, quasars are a popular target for measuring AGN feedback due to their brightness and their active feedback processes, but they have drawbacks due to their relative scarcity and the contaminating emission they contain that obscures the tSZ signatures of AGN feedback. In this paper, we follow Spacek et al. (2016) and focus on measuring co-added tSZ distortions in the CMB around massive ($\geqslant {10}^{11}\,{M}_{\odot }$) quiescent elliptical galaxies at moderate redshifts ($0.5\leqslant z\leqslant 1.5$). To accomplish this, we use data from the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), SDSS (Alam et al. 2015), and ACT (Dünner et al. 2013). These galaxies contain almost no dust and are numerous on the sky, making them well-suited for co-adding in large numbers in order to obtain good constraints on the energy stored in the surrounding gas.

The structure of this paper is as follows: in Section 2, we explain our method of measuring the thermal energy around galaxies using the tSZ effect and how that might relate to non-gravitational heating by AGN feedback. In Section 3, we discuss the data that we use to both select galaxies and make our tSZ measurements. In Section 4, we outline our galaxy selection process and how we estimate the physical parameters of the galaxies. In Section 5, we explain how we filter the ACT images. In Section 6, we discuss our stacking procedure and results. In Section 7, we use a ${\chi }^{2}$ analysis to model and remove contaminant signal, and in Section 8 we do the same but with data from Planck included. In Section 9, we summarize our results, discuss the implications for AGN feedback, and provide conclusions.

Throughout this work, we adopt a ΛCDM cosmological model with parameters (from Planck Collaboration et al. 2015d), h = 0.68, ${{\rm{\Omega }}}_{0}=0.31$, ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.69$, and ${{\rm{\Omega }}}_{b}=0.049$, where h is the Hubble constant in units of 100 km s⁻¹ Mpc⁻¹, and ${{\rm{\Omega }}}_{0}$, ${{\rm{\Omega }}}_{{\rm{\Lambda }}}$, and ${{\rm{\Omega }}}_{b}$ are the total matter, vacuum, and baryonic densities, respectively, in units of the critical density. All of our magnitudes are quoted in the AB system (i.e., Oke & Gunn 1983).

The tSZ effect is the process by which CMB photons gain energy when passing through ionized gas (Sunyaev & Zeldovich 1970, 1972) through inverse Compton scattering with energetic electrons. The resulting CMB anisotropy has a distinctive frequency dependence, which causes a deficit of photons at frequencies below ${\nu }_{\mathrm{null}}=217.6\,\mathrm{GHz}$ and an excess of photons above ${\nu }_{\mathrm{null}}$. The change in CMB temperature ${\rm{\Delta }}T$ as a function of frequency due to the (nonrelativistic) tSZ effect is

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}T}{{T}_{\mathrm{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 }_{T}\displaystyle \frac{{n}_{e}k({T}_{e}-{T}_{\mathrm{CMB}})}{{m}_{e}{c}^{2}},\end{eqnarray} \tag{ 2 }$$

where ${\sigma }_{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}_{\mathrm{CMB}}$ is the CMB temperature (we use ${T}_{\mathrm{CMB}}=2.725$ K), the integral is performed over the line-of-sight distance l, and the dimensionless frequency x is given by $x\equiv \tfrac{h\nu }{{{kT}}_{\mathrm{CMB}}}=\tfrac{\nu }{56.81\,\mathrm{GHz}}$, where h is the Planck constant. We can calculate the total excess thermal energy ${E}_{\mathrm{therm}}$ associated with a source by integrating Equation (2) over a region of sky around the source, as detailed in Spacek et al. (2016), and combining the result with Equation (1) to get ${E}_{\mathrm{therm}}$ as a function of x and ${\rm{\Delta }}T$. This gives

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

To compare the measured energies above with the expectations from models of feedback, we use the simple models of gas heating with and without AGN feedback worked out in Spacek et al. (2016). For purely gravitational heating, we can assume that the gas collapses and virializes along with an encompassing spherical halo of dark matter. The gas is shock-heated during gravitational infall to a virial temperature ${T}_{\mathrm{vir}}$, and if we approximate the gas as isothermal at this temperature we can estimate its total thermal energy as

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

where m_p is the proton mass, $\mu =0.62$ is the average particle mass in units of m_p, and M₁₃ is the mass of the halo in units of ${10}^{13}\,{M}_{\odot }.$ We can convert from halo mass to the stellar mass of the galaxies we will be measuring if we use the observed relation between black hole mass and halo circular velocity v_c from Ferrarese (2002), convert the black hole mass to its corresponding bulge dynamical mass using a factor of 400 (Marconi & Hunt 2003), use the fact that ${v}_{c}={({GM}/{R}_{\mathrm{vir}})}^{1/2}$ = $254\,\mathrm{km}\ {{\rm{s}}}^{-1}\,{M}_{13}^{1/3}{(1+z)}^{1/2}$, where G is the gravitational constant, and take ${M}_{\mathrm{stellar}}\propto {v}_{c}^{5}$. As shown in Spacek et al. (2016), this gives ${M}_{\mathrm{stellar}}={2.8}_{-1.4}^{+2.4}\,\times {10}^{10}\,{M}_{\odot }\,{M}_{13}^{5/3}{(1+z)}^{5/2}$, and substituting this into Equation (4) yields

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

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

For heating due to AGN feedback it is difficult to be precise because little is known about the dominant mechanism by which AGN feedback operates, and as a result there are many models, each of which leads to somewhat different signatures in our data. We can, however, try to estimate the overall magnitude of AGN feedback heating by making use of the simple model described in Scannapieco & Oh (2004). This is characterized as the heating of gas by a fraction ${\epsilon }_{k}$ of the total bolometric luminosity of the AGN, where the black hole shines at the Eddington luminosity ($1.2\times {10}^{38}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{M}_{\odot }^{-1}$) for a time $0.035\,{t}_{\mathrm{dynamical}}$, with ${t}_{\mathrm{dynamical}}\equiv {R}_{\mathrm{vir}}/{v}_{c}$ $\,=2.6\,\mathrm{Gyr}\,{(1+z)}^{-3/2}$, where ${R}_{\mathrm{vir}}$ is the halo virial radius. This gives

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

where ${\epsilon }_{k,0.05}\equiv {\epsilon }_{k}/0.05$. In this case, 5% is a typical, though still very uncertain, efficiency needed to achieve anti-heirarchical galaxy evolution through effective feedback (e.g., Scannapieco & Oh 2004; Thacker et al. 2006; Costa et al. 2014).

It is evident that our simple model of feedback energy falls within the errors of our model for gravitational energy, indicating that the differences between models with and without AGN feedback are subtle. Detailed simulations beyond the scope of this paper will be needed to make precise predictions regarding particular AGN feedback models. Even so, our models above are roughly consistent with more sophisticated models (e.g., Thacker et al. 2006; Chatterjee et al. 2008), and we will therefore use them to provide a general context for our results.

In order to select a large number of galaxies for our stacking analysis, we wanted to use a large region of the sky that was covered with a wide wavelength range of telescope surveys and included microwave data for our tSZ measurements. We therefore chose the SDSS Stripe-82 region, which is covered by many surveys, including ultraviolet, visible, and infrared data from the SDSS, infrared data from the WISE All-Sky Data Release, and microwave data from the ACT. We also used the extensive pre-existing source catalogs corresponding to the SDSS and WISE data.

Our SDSS data were taken from Data Release 12 (DR12) of the third generation (SDSS-III; Alam et al. 2015). Since 2000 the SDSS has used a 2.5 m wide-field telescope at Apache Point Observatory in New Mexico to image roughly one-third of the sky (31,637 deg²), and the catalog contains information on over 1 billion objects. The SDSS bands we used are u, g, r, i, and z, with average wavelengths of 355.1, 468.6, 616.5, 748.1, and 893.1 nm, respectively, and an average point-spread function (PSF) width of 1.4 arcsec in r-band.¹ Stripe-82 is a 25 wide stripe along the celestial equator that was imaged multiple times, resulting in deeper SDSS data than the main survey. The general Stripe-82 region runs from −125 to 125 decl. and −65° to 60° R.A., with an area of ≈312 deg².²

The WISE All-Sky Data Release contains data from the full WISE mission in 2010 using the 0.4 m space telescope (Wright et al. 2010). The four WISE infrared bands are labelled W1, W2, W3, and W4, and are centered at 3.4, 4.6, 12, and 22 μm, respectively, with average PSF full-width-at-half-maximum (FWHM) values of 6.1, 6.4, 6.5, and 12.0 arcsec, respectively. The corresponding source catalog contains over 500 million sources with a signal-to-noise ratio (S/N) $\gt 5$.³

The ACT is a 6 m telescope on Cerro Toco in Chile which started observing in 2007. It was equipped with the Millimeter Bolometric Array Camera (MBAC), with bands at 148, 220, and 277 GHz (Dünner et al. 2013). The data used in this paper covering the equatorial Stripe-82 region are from ACT seasons 3 and 4 (2009 and 2010) using the 148 and 220 GHz bands.⁴ These have beam FWHM values of $\approx 1.44$ and $\approx 1.08$ arcmin, respectively. We used the data designated as "src_free," where flux from point sources has been removed. Since both seasons cover the same region of sky, we are able to average them together to increase our S/N.

The ACT bands at 148 and 220 GHz are ideal for our tSZ measurements because 148 GHz is close to the peak of the undistorted CMB spectrum (160 GHz) and will see a significant decrement, while 220 GHz is very close to a frequency where the tSZ effect has no effect (${\nu }_{\mathrm{null}}=217.6\,\mathrm{GHz}$). Equation (1) can be rewritten for the ACT bands after integrating over their filter curves and solving for the Compton-y parameter. The filter curves are taken from Swetz et al. (2011). This gives

$$\begin{eqnarray}&&y=-0.38\,\displaystyle \frac{{\rm{\Delta }}{T}_{148}}{1\,{\rm{K}}}\hspace{1cm}\mathrm{and}\hspace{1cm}y=5.4\,\displaystyle \frac{{\rm{\Delta }}{T}_{220}}{1\,{\rm{K}}},\end{eqnarray} \tag{ 7 }$$

where ${\rm{\Delta }}{T}_{148}$ and ${\rm{\Delta }}{T}_{220}$ are the temperature anisotropies at 148 and 220 GHz. We can similarly integrate Equation (3) over the ACT 148 GHz filter curve to get the total thermal energy surrounding a galaxy as a function of the 148 GHz tSZ decrement,

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

In order to select galaxies best suited for our tSZ measurements, we have followed the selection criteria in Spacek et al. (2016). We therefore have restricted our attention to massive elliptical galaxies with redshifts 0.5 ≤ z ≤ 1.5. Galaxies are initially selected from the WISE catalog to have equatorial coordinates that lie within SDSS Stripe-82 ($306^\circ \lt {\alpha }_{2000}\lt 62^\circ$ and $-1\buildrel{\circ}\over{.} 27\lt {\delta }_{2000}\lt 1\buildrel{\circ}\over{.} 27$). A cut was made requiring ${\rm{S}}/{\rm{N}}\gt 5$ in both W1 and W2.

We incorporated the SDSS bands in order to perform a color selection analogous to the passive BzK selection of Daddi et al. (2004) using g − z and $z-W1$ colors. This gzW1-selection is illustrated for several example spectral energy distributions (SEDs) from Bruzual & Charlot (2003) in the left plots of Figure 1. The selection lines (in the AB magnitude system) are

$$\begin{eqnarray}&&(z-W1)\leqslant (g-z)-0.02,\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&(z-W1)\geqslant 2.0,\end{eqnarray} \tag{ 10 }$$

which define the wedge in the upper right of the bottom plots in Figure 1. The color-selection was necessary to reduce the sample to a useable size. Galaxies whose colors lie in this wedge should be old and passively evolving galaxies at $1\lesssim z\lesssim 2$. Following the WISE and gzW1 criteria, and selecting only WISE sources with unique SDSS-DR12 matches, we arrived at a sample of ≈30,000 galaxies which were further pared down using redshift and SED parameters. We emphasize that we are after a highly reliable sample and are willing to sacrifice completeness in the interest of purity.

Zoom In Zoom Out Reset image size

Photometric redshifts were computed with extinction-corrected SDSS ugriz and WISE W1 and W2 photometry using EAZY (Brammer et al. 2008). The W3 and W4 bands were omitted since they are not comparably deep and are dominated by a warm/hot dust component that obscures the stellar component of the target galaxies that we are after. We then applied the constraint that 0.5 ≤ z ≤ 1.5. The seven-band SEDs were fit with Bruzual & Charlot (2003) exponentially declining star formation rate (with timescale τ) models. We only used those objects with reliable SED fits by selecting galaxies with reduced ${\chi }^{2}\leqslant 5$. We applied further selections based on the results of the SED fits by taking only galaxies with ages $\geqslant 1\,\mathrm{Gyr}$ and specific star formation rates ≤0.01 Gyr⁻¹. This ensured that we were choosing older galaxies that were not actively forming stars, especially excluding dusty starbursts which made it through the color-cuts. This gave a sample of ≈10,000 galaxies, which were then pruned of known contaminants to the tSZ signal.

In order to estimate the possible contamination of our sample, we can appeal to morphological measurements where the same SED selection criteria as above are applied. Unfortunately, we cannot use the sample we have selected from Stripe-82 because Hubble Space Telescope (HST) resolution is required to classify galaxies at $0.5\lesssim z\lesssim 1.5$. Instead we use the largest HST survey available, CANDELS (Grogin et al. 2011; Koekemoer et al. 2011). We use the SED fit parameters from Skelton et al. (2014) and Sérsic function fits from van der Wel et al. (2012), who used GALFIT (Peng et al. 2002) to fit the 2D light profiles to the CANDELS HST images in the ${H}_{F160W}$-band. Since galaxies with stellar mass greater than ${10}^{11}\,{M}_{\odot }$ are rare on the sky and even the largest HST survey is significantly smaller than the $\approx 300$ deg² Stripe-82 survey, we lower our mass limit to $M\geqslant {10}^{10}\,{M}_{\odot }$ for this exercise. We choose all galaxies with age greater than 1 Gyr and $\mathrm{SSFR}\lt 0.01$ Gyr⁻¹. Additionally, we require a reliable Sérsic fit from van der Wel et al. (2012) and choose Sérsic index $n\gt 2.5$ to be representative of "elliptical" galaxies. For samples of 364 and 346 galaxies, we find that 82% and 78% are $n\gt 2.5$ galaxies at $0.5\leqslant z\leqslant 1.0$ and $1.0\leqslant z\leqslant 1.5$, respectively. If we remove the age and SSFR requirements we get only 45% and 30% $n\gt 2.5$ ellipticals at low and high redshift, respectively. Therefore, assuming morphology and SED parameters are correlated in this way, we estimate that our sample of old galaxies that are not actively forming stars is $\simeq 80$% pure. Removing this 20% contamination is part of the focus of Sections 7 and 8.

The most significant contaminants that must be removed are known AGN and galaxy clusters. We therefore removed sources from the ROSAT Bright and Faint Source catalogs (BSC and FSC; Voges et al. 1999). We additionally removed known clusters from ROSAT (Piffaretti et al. 2011). Clusters selected via the SZ effect would completely counter our measurements, so we removed both known Planck (Planck Collaboration et al. 2015e) and ACT (Marriage et al. 2011; Hasselfield et al. 2013) clusters. X-ray sources from XMM-Newton and Chandra (LaMassa et al. 2013), sources from the AKARI/FIS Bright Source Catalog (Yamamura et al. 2010), and sources from the AKARI/IRC Point Source Catalog (Ishihara et al. 2010) were also removed. We also removed Galactic molecular clouds by cross-matching with the Planck Catalogue of Galactic Cold Clumps (Planck Collaboration et al. 2015f), compact sources from the nine-band Planck Catalog of Compact Sources (Planck Collaboration et al. 2014), and all sources from the IRAS Point Source Catalog (Helou & Walker 1988, pp. 1–265) and radio sources from Best & Heckman (2012). We also verified that none of our sources satisfied the "$W1W2$-dropout" criteria for extremely luminous infrared galaxies of Eisenhardt et al. (2012). In all cases, sources with a possible contaminant within 4.0 arcmin, approximately double our region of interest around each source, were flagged and those sources were removed from further consideration. This left ≈7200 massive, quiescent, $0.5\leqslant z\leqslant 1.5$ galaxies that are away from known potential contaminants.

Finally, to make sure we were selecting galaxies with the most reliable parameters, we limited the AB magnitude errors in the SDSS bands (ugriz), with $\mathrm{mag}\_\mathrm{error}$ <1.5 mag, required ${\mathrm{log}}_{10}(\mathrm{SSFR})$ to be finite, and limited the galaxy stellar mass as $M\lt {10}^{13}\,{M}_{\odot }$. This resulted in a final selection of 4453 galaxies to include in our tSZ stacks. To narrow down our measurements in redshift space we split our galaxy sample into two redshift bins: a "low-z" bin with 1179 $0.5\leqslant z\lt 1.0$ galaxies and a "high-z" bin with 3274 $1.0\leqslant z\leqslant 1.5$ galaxies. We show how our final galaxy selection fits in with our original SED color selection in the right plots of Figure 1. The locations of the final selection of galaxies is shown in Figure 2. The mass, redshift, and age distribution of the final sample is shown in Figure 3. Mean and mass-averaged values for redshift, angular diameter distance, mass, age, and W1 luminosity in both redshift bins are given in Table 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Mean and Mass-averaged Values for Several Relevant Galaxy Parameters in the Two Final Redshift Ranges

Cut	z	N	$\langle z\rangle$	$\langle {l}_{\mathrm{ang}}^{2}\rangle$	$\langle {\rm{M}}\rangle$	$\langle \mathrm{Age}\rangle$	$\langle {{\rm{L}}}_{W1}\rangle$	${\langle z\rangle }_{{\rm{M}}}$	${\langle {l}_{\mathrm{ang}}^{2}\rangle }_{{\rm{M}}}$
				(Gpc2)	(${{\rm{M}}}_{\odot }$)	(Gyr)	(erg s−1 Hz−1)		(Gpc2)
All	0.5–1.0	1179	0.83	2.56	7.81 × 1011	3.80	7.83 × 1030	0.86	2.61
All	1.0–1.5	3274	1.20	3.04	10.1 × 1011	3.56	12.8 × 1030	1.21	3.05
Planck	0.5–1.0	227	0.83	2.55	6.93 × 1011	3.63	7.04 × 1030	0.86	2.60
Planck	1.0–1.5	529	1.21	3.05	9.68 × 1011	3.44	12.4 × 1030	1.21	3.05

Note. "All" represents our complete, final galaxy sample, and "Planck" represents our final galaxy sample with further cuts applied, as discussed in Section 8.

Download table as:  ASCIITypeset image

Before stacking the ACT data around our selected galaxies, we needed to filter the ACT maps to remove the primary CMB anisotropy and maximize the signal-to-noise at the spatial scales we are measuring. An ideal Fourier-space point source filter is given by $\psi =\tfrac{\tau }{P}{\left[\int {d}^{2}k\tfrac{{\tau }^{2}}{P}\right]}^{-1}$, where τ is the Fourier-space source profile and P is the Fourier-space noise covariance matrix (e.g., Haehnelt & Tegmark 1996). For the source profile we assumed a slightly extended source such that $\tau =B\times G$, where B is the Fourier-space beam function and G is a Fourier-space Gaussian function. We have also approximated the noise P as the ACT noise power spectra given in Das et al. (2014) for seasons 3 and 4 plus the CMB power spectrum. We therefore have approximated the filter as

$$\begin{eqnarray}&&\psi \approx \displaystyle \frac{B\times G}{N}{\left[\int {d}^{2}k\displaystyle \frac{{(B\times G)}^{2}}{N}\right]}^{-1},\end{eqnarray} \tag{ 11 }$$

where N is the Fourier-space CMB+noise power spectrum. We did this for each band in each season and then averaged the two seasons together, resulting in an averaged filter for each band. For G we chose a Gaussian with a FWHM of 1.5 arcmin to represent a slightly extended source. This is because our signal of interest is from hot gas within and surrounding the galaxies that likely represents the cumulative heating due to multiple cycles of AGN activity. We have no way of knowing the true shape of this gas but we expect the tSZ signal to be greatest near the galaxy and decrease away from it, meaning it is simplest to assume a slightly extended Gaussian profile. We note that this Gaussian we use is slightly larger than the ACT 148 GHz beam, which has a FWHM of 1.44 arcmin.

We then scaled the filters so that the flux within a 1 arcmin radius aperture is preserved in our maps after filtering, representing the regions of interest we measure around the selected galaxies. This choice of aperture follows from Spacek et al. (2016), where it is noted that the energy input from AGN feedback is unlikely to affect scales much larger than twice the dark matter halo virial radius, which corresponds to about 2 arcmin at the redshifts we are investigating. Although the gas surrounding the galaxies will have both intrinsically different angular sizes and different angular sizes due to their differing redshifts, we have no way of knowing these sizes and the best we can do is take an aperture that is not unnecessarily large and that we expect to contain most of the signal in all cases. Additionally, we want these measurements to be comparable with cosmological AGN feedback simulations, and that is easiest to do with a constant measurement aperture. The factors used to scale the filters are 0.0167 and 0.0162 for 148 and 220 GHz, respectively. The final averaged, scaled filters are applied to the corresponding ACT maps for both seasons. Pictures of the final scaled Fourier-space filters are shown in Figure 4. The lack of smoothness is due primarily to the CMB power spectrum.

Zoom In Zoom Out Reset image size

To stack the CMB data we first made a 8.4 × 8.4 arcmin (17 × 17 pixel) stamp around each galaxy at 148 and 220 GHz in the filtered ACT data for seasons 3 and 4. Then we averaged the individual stamps together to make two stacked stamps for each band in each season, split into low-z and high-z galaxies. Finally, we averaged the seasons together. The resulting stamps are shown in Figure 5, with scales centered around 0 in units of μK. Any pixels >4 arcmin away from the center were set to 0 since that was the distance of our potential contaminant cuts. We get our final measurements integrated over the sky by summing the stacked signal within a 1 arcmin radius aperture (corresponding to a 2 pixel radius), shown as black circles in Figure 5.

Zoom In Zoom Out Reset image size

The upper left panel of this figure shows a signal close to zero, while the lower left panel shows a clear positive signal in the center. These are the low-z and high-z 148 GHz stamps, respectively, and there is no tSZ detection in our initial stacks, which would be a negative signal at 148 GHz. In fact, at least at high-z, there is a significant contaminant signal. Looking at the right panels, at 220 GHz, we see that the stamps are even more dominated by positive contaminant signals. Since the tSZ effect has a negligible impact at this frequency, indicated by Equation (7), this indicates that our galaxy selection process was imperfect, and there still remains a positive contaminating signal composed of faint sources that we were unable to account for. Looking at a typical range of emission by dust at z = 1 (light and dark blue curves in Figure 6), the CMB spectrum (green curve in Figure 6), and the ACT bands (rightmost red hatched region in Figure 6), it seems likely that this contaminating signal at 220 GHz also extends into the 148 GHz band. It is therefore likely that we are in fact seeing a significant tSZ signal that is obscured by contaminant emission.

Zoom In Zoom Out Reset image size

In order to estimate the uncertainty in our final measurements we generated 429,571 random points in our field on the sky, a number chosen by dividing the area of our field by the area of the 148 GHz beam which we approximated as $2\pi {\sigma }^{2}$, where σ is the Gaussian beam standard deviation. We then applied the same 4 arcmin contaminant source cuts as we applied to our galaxy selection, leaving us with 294,176 random points. We stacked these random points on the sky for each band in the same way as we stacked our galaxies, and we computed the corresponding 1 arcmin radius aperture sums. First we computed an overall offset from the random points by getting the mean value of these sums. Since we remove galaxies anywhere near potential contaminants that might have positive or negative signal, we inherently bias the zeropoint of the ACT maps. We therefore corrected our aperture sums by subtracting off these mean offset values. These season-averaged offsets are 0.014 and −0.35 μK arcmin² for 148 and 220 GHz, respectively. Next we used the random point sums to compute the variance for an individual measurement in each case. If we first define the normal variance of N random values x as

$$\begin{eqnarray}&&\mathrm{var}=\displaystyle \frac{1}{N-1}\displaystyle \sum _{i=1}^{N}{({x}_{i}-\langle x\rangle )}^{2},\end{eqnarray} \tag{ 12 }$$

where N = 294,176, we can then define our overall uncertainty σ as a combination of measurement error and offset error, given by

$$\begin{eqnarray}&&\sigma ={\left(\displaystyle \frac{\mathrm{var}}{n}+\displaystyle \frac{\mathrm{var}}{N}\right)}^{1/2},\end{eqnarray} \tag{ 13 }$$

where n is the number of galaxy measurements (1179 for low-z and 3274 for high-z). The result is final co-added sums and uncertainties for each individual season-averaged band and redshift bin, and these are given in Table 2.

From this table, we can directly see that there is a $\approx 1\mbox{--}2\sigma$ contaminant signal at 148 GHz, and a $\approx 3\mbox{--}6\sigma$ contaminant signal at 220 GHz. It is clear that obtaining the best possible constraints on non-gravitational heating and AGN feedback requires making the best possible separation between the tSZ signal and the contaminating signal, which is addressed in the following sections. Finally, we convert our co-added ${\rm{\Delta }}T$ signal into gas thermal energy using Equation (8). These values are shown in Table 3, under "Data only."

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

Model	N	z	$\int {\rm{\Delta }}{T}_{150}({\boldsymbol{\theta }})d{\boldsymbol{\theta }}$	Y	${E}_{\mathrm{therm}}(\pm 1\sigma )$	${E}_{\mathrm{therm}}(\pm 2\sigma )$	S/N
			(μK arcmin2)	(10−7 Mpc2)	(1060 erg)	(1060 erg)	((Etherm/1σ)
Data only	1179	0.5–1.0	1.0 ± 1.4	−0.8 ± 1.1	−2.8 ± 3.9	−2.8 ± 7.9	−0.72
	3274	1.0–1.5	2.1 ± 0.9	−1.9 ± 0.8	−7.0 ± 3.0	−7.0 ± 6.0	−2.33
${\chi }^{2}$ (ACT only)	1179	0.5–1.0	$-{1.6}_{-1.9}^{+1.8}$	${1.2}_{-1.4}^{+1.4}$	${4.5}_{-5.1}^{+5.4}$	${4.5}_{-10.7}^{+11.0}$	0.85
	3274	1.0–1.5	$-{2.1}_{-1.2}^{+1.1}$	${1.9}_{-1.0}^{+1.1}$	${7.0}_{-3.7}^{+4.0}$	${7.0}_{-7.7}^{+8.4}$	1.78
${\chi }^{2}$ (With Planck)	227	0.5–1.0	$-{2.0}_{-2.1}^{+2.0}$	${1.5}_{-1.5}^{+1.6}$	${5.6}_{-5.6}^{+5.9}$	${5.6}_{-11.6}^{+12.4}$	0.97
	529	1.0–1.5	$-{2.1}_{-1.4}^{+1.3}$	${1.9}_{-1.2}^{+1.3}$	${7.0}_{-4.4}^{+4.7}$	${7.0}_{-9.4}^{+9.7}$	1.50

Note. The last three columns represent the best-fit ${E}_{\mathrm{therm}}$ values with $\pm 1\sigma$ values and $\pm 2\sigma$ values and the ${E}_{\mathrm{therm}}$ signal-to-noise ratio (E_therm/1σ), respectively.

Download table as:  ASCIITypeset image

As evidenced by Table 2, there appears to still be a significant contamination signal indicated by the large positive 220 GHz values, which is likely contributing to the 148 GHz values that we are interested in for our tSZ measurements. This is illustrated by the blue lines in Figure 6, where it is clear that dust at reasonable temperatures around z = 1 will have significant emission in the ACT bands we are using. In order to constrain and subtract out this undetected contamination, we have followed the process described in Spacek et al. (2016) and built a detailed model of contaminants based on extrapolations of the source counts measured for SPT data in Mocanu et al. (2013). We extended these source counts to fainter values by modeling a random population of undetected sources that follow the trend of the detected sources into the unresolved region. We then related these models to the contaminating signal in our 148 and 220 GHz measurements.

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

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

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}_{\max }$ is the flux at which we expect all brighter sources to have a 100% completeness level in the source count catalog. We then computed 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 220 GHz were ${\alpha }_{s}=-2.08\pm 0.09$ for the synchrotron sources and ${\alpha }_{d}=-2.91\pm 0.17$ 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 modeled this fall-off as a minimum flux ${S}_{\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}_{\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 decided 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. We used a maximum flux cutoff of 305.7 μK arcmin², corresponding to the faintest bin of detected sources found by Mocanu et al. (2013). (ii) We then assigned the source a random 220 GHz flux, ${S}_{220,\mathrm{rand}},$ by inverting

$$\begin{eqnarray}&&{\int }_{{S}_{220,\min }}^{{S}_{220,\mathrm{rand}}}{dS}\,\displaystyle \frac{{dN}}{{dS}}=R{\int }_{{S}_{220,\min }}^{{S}_{220,\max }}{dS}\,\displaystyle \frac{{dN}}{{dS}},\end{eqnarray} \tag{ 15 }$$

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

$$\begin{eqnarray}&&{S}_{220,\mathrm{rand}}={[(1-R){S}_{220,\min }^{\alpha +1}+R{S}_{220,\max }^{\alpha +1}]}^{\tfrac{1}{\alpha +1}};\end{eqnarray} \tag{ 16 }$$

(iii) To obtain a corresponding flux for the source at 148 GHz we used 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 chose ${\alpha }_{220}^{150}$ values that fit these distributions and calculated the 148 GHz flux (following Mocanu et al. 2013) as

$$\begin{eqnarray}&&{S}_{148,\mathrm{rand}}=\displaystyle \frac{{S}_{220,\mathrm{rand}}}{{C}_{1}\times {C}_{2}},\end{eqnarray} \tag{ 17 }$$

where C₁ is the conversion factor between Jy and μK arcmin² integrated over the band filter curves, C₂ is the relating factor ${({\nu }_{220}/{\nu }_{148})}^{{\alpha }_{220}^{150}}$ integrated over the band filter curves, and we used units of μK arcmin² for all S. (iv) Finally, we estimated the completeness of our 220 GHz measurements and randomly discarded modeled sources to match the estimated fraction of 220 GHz sources detected per flux. To do the estimation, we assumed a cutoff S_cut of 250 μK arcmin², representing the $3\sigma$ limit of our 220 GHz source measurements. We discarded galaxies with signals greater than S_cut or signals less than $-{S}_{\mathrm{cut}}$. We then determined the completeness fractions for the modeled sources, and accounted for our measurement uncertainty, by adding the 220 GHz random point distribution onto S_cut in a cumulative manner using the fraction of random point measurements below a given flux. This means S_cut has a 50% completeness, while fainter fluxes are increasingly less complete and brighter fluxes are increasingly more complete.

For each source-count model, we repeated the process 100,000 times, resulting in a large catalog of contaminating fluxes in both bands. From these, we computed the mean flux per contaminating source in each band, $\langle {S}_{150,\mathrm{cont}}\rangle$ and $\langle {S}_{220,\mathrm{cont}}\rangle$, which represents the contamination we are measuring in our stacks. To account for variations in the input parameters, we computed model contamination signals for a wide range of source-count models, with ${S}_{220,\max }=305.7\,\mu {\rm{K}}$ arcmin². We varied ${\alpha }_{s}$ from −2.26 to −1.90 in steps of 0.09 and we varied ${\alpha }_{d}$ from −3.25 to −2.57 in steps of 0.17, representing ranges of $\pm 2\sigma$ in steps of σ. We let log₁₀(${S}_{\min }$) vary from log₁₀(0.01 μK arcmin²) to log₁₀(30 μK arcmin²) in steps of 0.2 in log-space.

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

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

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

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

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{ 20 }$$

Here, S₁₅₀, S₂₂₀, ${\sigma }_{150}$, and ${\sigma }_{220}$ are our measured 1 arcmin radius values from Table 2. The σ values are computed using random point measurements, given by

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

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

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

Our approach was thus to marginalize over values of ${f}_{\mathrm{cont}}$ in the full physical range from 0 to $1,$ but normalize the overall probability by the sum of ${f}_{\mathrm{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}_{\mathrm{cont}}$ outside the range of physically possible values.

Equation (22) then gives us a function $P({S}_{\mathrm{SZ}})$ for each combination of ${\alpha }_{d},$ ${\alpha }_{s},$ and ${S}_{\min }$. We can convert the corresponding ${S}_{\mathrm{SZ}}$ value to the gas thermal energy, ${E}_{\mathrm{therm}},$ using Equation (8) and the average ${l}_{\mathrm{ang}}^{2}$ from Table 1. Note that a positive detection of the tSZ effect is seen as a negative ${\rm{\Delta }}T$ signal at 148 GHz, and it represents a positive injection of thermal energy into the gas around the galaxy. Additionally, we compute a corresponding range for ${E}_{\mathrm{grav}}$ using Equation (5) and values from Table 1. The peak of each $P({S}_{\mathrm{SZ}})$ curve is shown as the colored points in Figure 7, where ${\alpha }_{s}$ (represented by point size) is increasing (i.e., becoming more positive) downwards, and ${\alpha }_{d}$ (represented by point color) is increasing upwards. The 1σ and 2σ contours are computed for each ${S}_{\min }$ by averaging $P({S}_{\mathrm{SZ}})$ across ${\alpha }_{d}$ and ${\alpha }_{s}.$ The resulting probability distribution depends only on ${E}_{\mathrm{therm}}$ and ${S}_{\min },$ and 1σ and 2σ are represented by the values $P({S}_{\mathrm{SZ}})=0.61$ and 0.13, respectively (i.e., $\exp [-{\sigma }^{2}/2]$). These contours are shown in Figure 7, along with the $-1\sigma$ range for ${E}_{\mathrm{grav}}$. Peak probability values are chosen using $P({S}_{\mathrm{SZ}})\gt 0.99$, shown as the black regions in Figure 7. From this figure we see that there is a $\approx 1\sigma$ tSZ detection for ${S}_{\min }\gtrsim 3\,\mu {\rm{K}}$ arcmin² at low-z. At high-z we see a $\approx 2\sigma$ tSZ detection for ${S}_{\min }\gtrsim 5\,\mu {\rm{K}}$ arcmin².

Zoom In Zoom Out Reset image size

Finally, we average the probability distribution across ${S}_{\min }$ to get a final distribution as a function of only ${E}_{\mathrm{therm}}$. The significance values of this curve are shown in Table 3 under "ACT only." We see a $0.9\sigma$ tSZ detection at low-z and a $1.8\sigma$ detection at high-z.

As indicated by Figure 6, we should be able to better constrain the contamination due to undetected dusty sources by incorporating data at higher frequencies than the ACT bands we are using to make our measurements. We therefore made use of the 2015 public data release from the Planck mission, and focused on the high-frequency bands at 217, 353, 545, and 857 GHz. These data, with a ≈5 arcmin FWHM beam (Planck Collaboration et al. 2015a), is too low-resolution to be useful in our direct tSZ measurements, but can still provide useful information about our galaxies at higher frequencies where the contaminant emission should be much brighter (see the light blue and dark blue curves in Figure 6). By incorporating these Planck measurements we should be able to better discriminate between contaminant models, allowing us to better identify the true tSZ signal.

To utilize the Planck data, we followed the same process as in the previous section by computing ${\chi }^{2}$ values for a number of modeled contaminants, but now we had several extra terms in each computed ${\chi }^{2}$ relating to the Planck measurements. In order to stack our galaxies in the Planck data, we first extended our contaminant source cut distance from 4 to 10 arcmin due to the much lower resolution. This resulted in a significant decrease in our number of galaxies, with 227 at low-z and 529 at high-z. In order to filter out the primary CMB signal, we convolved each Planck map with a 7 arcmin FWHM Gaussian and subtracted the resulting image from the original. We then stacked the central pixels of each galaxy to get co-added values in each of the Planck bands. In addition, we degraded the ACT 148 and 220 GHz maps to match the Planck beam, applied the same 7 arcmin FWHM filtering, and stacked the central pixels of galaxies in those images as well.

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

Finally, we computed 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, discussed below, we removed the corresponding uncertainty term, taken to have a covariance of 7.85 μK as used in Spacek et al. (2016). 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})}_{\mathrm{corr}}=\sqrt{\displaystyle \frac{{\sigma }_{i}{\sigma }_{j}-{\sigma }_{\mathrm{cov}}^{2}}{{N}_{\mathrm{source}}}+\displaystyle \frac{{\sigma }_{i}{\sigma }_{j}}{{N}_{\mathrm{random}}}},\end{eqnarray} \tag{ 23 }$$

where ${\sigma }_{i}{\sigma }_{j}$ is given by Equation (21) with i and j representing the various bands used, ${\sigma }_{\mathrm{cov}}=7.85\,\mu {\rm{K}}$ is the minimum CMB covariance discussed above, ${N}_{\mathrm{source}}$ is the number of sources used for the measurements (227 for low-z and 529 for high-z), and ${N}_{\mathrm{random}}={\rm{54,962}}$ is the number of random points used. 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. 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 modeled ACT 148 and 220 GHz contaminant source fluxes using a range of different source-count models (i.e., ${\alpha }_{d},$ ${\alpha }_{s},$ and ${S}_{\min }$), resulting again in 100,000 modeled contaminating source fluxes in each ACT band, ${S}_{148,\mathrm{cont}}$ and ${S}_{220,\mathrm{cont}}$. We also modeled what the contaminating signal would be in the Planck bands and the ACT bands filtered to match Planck. For each modeled contaminating source, if it was chosen to be a synchrotron source we simply extrapolated the Planck-based fluxes as

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

integrated over the relevant band filter curves, where ${\alpha }_{220}^{150}$ is the same used in the previous section, C_ν is a frequency-dependent factor involved in the conversion from Jy to μK arcmin², 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 adopted a modified blackbody with a free emissivity index, β, and dust temperature, ${T}_{\mathrm{dust}},$ often referred to as a gray-body (Planck Collaboration et al. 2015c). This requires us to add another free parameter, the temperature of the contaminant dust, ${T}_{\mathrm{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\mathrm{GHz}}=3+\beta -{x}_{185}{[1-\exp (-{x}_{185})]}^{-1},\end{eqnarray} \tag{ 25 }$$

where ${x}_{185}\equiv (185$ 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\,\mathrm{GHz}$ because it is halfway between our two ACT bands (148 and 220 GHz). This can be related, in turn, to the power law index ${\alpha }_{220}^{150}$ as

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

This then gives

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

integrated over the relevant band filter curves, where we vary ${T}_{\mathrm{dust}}$ from 20 to 50 K in steps of 3 K.

With these expressions, we were able to compute ${\chi }^{2}$ values for each source-count model accounting for the Planck measurements. This time, in addition to varying ${f}_{\mathrm{cont}}$ and ${S}_{\mathrm{SZ}}$, we also varied ${T}_{\mathrm{dust}}$ (as discussed above) and a parameter Δ, 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 involved the original ACT terms plus the new Planck terms, and it followed the same process as in Equation (18),

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

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 denote each element of the signal array ${{ \mathcal B }}_{i}$, where i runs over the two ACT bands (i.e., 148 and 220 GHz) and then every Planck-filtered band (i.e., the Planck bands at 857, 545, 353, and 217 GHz, plus the ACT bands at 220 and 148 GHz filtered to match the Planck images), such that ${{ \mathcal B }}_{1}={f}_{\mathrm{cont}}\times \langle {S}_{148,\mathrm{cont}}\rangle$ + ${S}_{\mathrm{SZ}}-{S}_{148}$, ${{ \mathcal B }}_{2}={f}_{\mathrm{cont}}$ × $\langle {S}_{220,\mathrm{cont}}\rangle -{S}_{220}$, and ${{ \mathcal B }}_{3-8}={f}_{\mathrm{cont}}\times \langle {S}_{3-8,\mathrm{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 (23).

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

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

where the whole function is normalized to a total of 1, and each final SZ value contains the sum over the corresponding ${T}_{\mathrm{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 four degrees of freedom. Thus the minimum ${\chi }^{2}$ was not 0 in every case as it was above with just 2 measurements and 2 fit parameters, and so for each model we scale the final probabilities by $\exp (-{\chi }_{\min }^{2}/2)$, where ${\chi }_{\min }^{2}$ is the minimum ${\chi }^{2}$ value for that model.

This again gave us a function $P({S}_{\mathrm{SZ}})$ for each combination of ${\alpha }_{d},$ ${\alpha }_{s},$ and ${S}_{\min },$ which we can convert to an energy ${E}_{\mathrm{therm}}$. The peak of each $P({S}_{\mathrm{SZ}})$ curve is shown as the colored points in Figure 8, where the points are colored by the minimum ${\chi }^{2}$ value for each model. The 1σ and 2σ contours are created for each ${S}_{\min }$ by averaging $P({S}_{\mathrm{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, with peak probability values represented by $P({S}_{\mathrm{SZ}})\gt 0.99$. These contours are shown in Figure 8, along with the $-1\sigma$ range for ${E}_{\mathrm{grav}}$. From this figure we can see that, for low-z, including Planck has slightly increased the estimated tSZ effect, though it now favors higher ${S}_{\min }$ values. The high-z result has not changed much besides an increased uncertainty due to fewer galaxies. For both redshift bins, the σ values have increased due to the large decrease in the number of stacked galaxies because of Planck's much larger beam.

Zoom In Zoom Out Reset image size

Finally, we average the probability distribution across ${S}_{\min }$, divided by the maximum value again, and get a final distribution as a function of only ${E}_{\mathrm{therm}}$. The significance (σ) values of this curve are shown in the "With Planck" part of Table 3. At low-z, the significance of our tSZ detection has increased from $0.9\sigma$ to $1.0\sigma$, while at high-z the tSZ detection significance decreases from $1.8\sigma$ to $1.6\sigma$. It is clear that the gain in sensitivity with Planck 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 the ACT. To clearly show this, we followed the methods of Section 7 (i.e., only ACT data were used), using this limited galaxy sample. The result is a $0.1\sigma$ tSZ detection at low-z and a $0.5\sigma$ tSZ detection at high-z. It is apparent, then, that adding the Planck data helps immensely with constraining the tSZ signal when the same number of galaxies are used, but since we have to limit our galaxy sample size so much to avoid contaminants in the Planck beam, the loss of accuracy due to fewer measurements just about offsets the gain in accuracy given by the added Planck data.

Alternatively, we can also characterize the total tSZ signal for our co-adds with the angularly integrated Compton-y parameter, Y. 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. Using Equation (7) at 148 GHz, this is

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

such that $Y=2.7\times {10}^{-8}\,{\mathrm{Mpc}}^{2}\,{E}_{60}$, where ${E}_{60}$ is ${E}_{\mathrm{therm}}$ in units of 10⁶⁰ erg. We can use this to compare the results in this paper with the similar work in Spacek et al. (2016), with a detailed comparison shown in Table 4. Comparing the two results, we see a decrease in Y at low-z in this work compared to Spacek et al. (2016), and similar Y values at high-z, but we also see a significant increase in the average galaxy mass. This is contrary to the expected trend of higher Y with higher mass seen in previous work (e.g., Planck Collaboration et al. 2013; Greco et al. 2015; Ruan et al. 2015).

Table 4.  A Comparison between Spacek et al. (2016) and the Current Work

Study	N	Type	z (mean)	Mass (${M}_{\odot }$)	tSZ $Y({10}^{-7}$ Mpc2)
Spacek et al. (2016)	3394	SPT	0.5–1.0 (0.72)	$1.51\times {10}^{11}$	${2.3}_{-0.7}^{+0.9}$
Spacek et al. (2016)	924	SPT	1.0–1.5 (1.17)	$1.78\times {10}^{11}$	${1.9}_{-2.0}^{+2.4}$
Spacek et al. (2016)	937	SPT+Planck	0.5–1.0 (0.72)	$1.51\times {10}^{11}$	${2.2}_{-0.7}^{+0.9}$
Spacek et al. (2016)	240	SPT+Planck	1.0–1.5 (1.17)	$1.78\times {10}^{11}$	${1.7}_{-1.8}^{+2.2}$
Current	1179	ACT	0.5–1.0 (0.83)	$7.81\times {10}^{11}$	${1.2}_{-1.4}^{+1.4}$
Current	3274	ACT	1.0–1.5 (1.20)	$10.1\times {10}^{11}$	${1.9}_{-1.0}^{+1.1}$
Current	227	ACT+Planck	0.5–1.0 (0.83)	$6.93\times {10}^{11}$	${1.5}_{-1.5}^{+1.6}$
Current	529	ACT+Planck	1.0–1.5 (1.21)	$9.68\times {10}^{11}$	${1.9}_{-1.2}^{+1.3}$

Note.  Y is the angularly integrated Compton-y parameter given by Equation (30). Mass refers to stellar mass.

Download table as:  ASCIITypeset image

We can compare with this previous work by defining $\widetilde{Y}\equiv \tfrac{Y}{{l}_{\mathrm{ang}}^{2}}\times {E}^{-2/3}\times {\left(\tfrac{{l}_{\mathrm{ang}}}{500\mathrm{Mpc}}\right)}^{2}$, where E(z) is the Hubble parameter, and this is shown in Figure 9. The circles represent Spacek et al. (2016) (black) and this work (red), both using measurements of massive quiescent elliptical galaxies with average redshifts greater than 0.7 for the lower-mass values of either color and 1.1 for the higher-mass values of either color. The blue squares represent Planck Collaboration et al. (2013) measurements of locally brightest galaxies, with redshifts less than ≈0.3. The orange diamonds represent Greco et al. (2015) measurements of locally brightest galaxies, with redshifts less than ≈0.3. The light blue triangles represent Ruan et al. (2015) measurements of locally brightest galaxies, with median redshifts of $\approx 0.5$ and redshifts less than ≈0.8. Looking at this figure, we see that the results of Spacek et al. (2016) are roughly consistent with the previous tSZ measurements, while the results of this work are over $\approx 2$ orders of magnitude smaller than previous tSZ measurements at the same mass. We note that the measurements of locally brightest galaxies from Planck Collaboration et al. (2013), Ruan et al. (2015), and Greco et al. (2015) are of significantly lower redshifts than our galaxies, with the highest overall redshift being no more than $\approx 0.8$ from Ruan et al. (2015) while our average redshifts range from 0.7 to 1.2. In addition, our selection criteria involve choosing quiescent elliptical galaxies and removing any galaxies in or around all detectable clusters, while the low-redshift locally brightest galaxies of the previous studies are more likely to be found in the centers of massive galaxy groups and clusters. Redshift alone cannot account for the $\approx 2$ orders of magnitude difference in measurements, as is made clear by the black and red lines representing our models at z = 0.8 (without and with AGN feedback, respectively) and the blue and orange lines representing our models at z = 1.2 (without and with AGN feedback, respectively). However, the differences in redshift combined with the different galaxy selection methods suggest that the galaxies used in this paper could be from fundamentally different populations and environments.

Zoom In Zoom Out Reset image size

With Equations (5) and (6) and the redshifts and masses from Table 1, we can also investigate theoretical thermal energies of the gas around elliptical galaxies due to both gravity and AGN feedback. Without Planck, we estimate the gravitational heating energy to be ${E}_{\mathrm{therm},\mathrm{grav}}={17.0}_{-9.1}^{+17.0}\times {10}^{60}$ erg for our low-z sample and ${E}_{\mathrm{therm},\mathrm{grav}}={16.7}_{-9.0}^{+16.7}\times {10}^{60}$ erg for our high-z sample. We therefore measure excess non-gravitational energies of ${E}_{\mathrm{therm},\mathrm{feed},\mathrm{dat}}=-{12.5}_{-10.6}^{+17.7}\times {10}^{60}$ erg for low-z and ${E}_{\mathrm{therm},\mathrm{feed},\mathrm{dat}}=-{9.7}_{-9.8}^{+17.1}\times {10}^{60}$ erg for high-z, both of these values consistent with zero detection. For completeness, we can plug these into the theoretical AGN feedback energy equation and solving for ${\epsilon }_{k}$, we get feedback efficiencies of $-{4.8}_{-4.1}^{+6.8}$% for low-z and $-{3.8}_{-3.9}^{+6.7}$% for high-z.

With Planck, we estimate the gravitational heating energy to be ${E}_{\mathrm{therm},\mathrm{grav}}={15.1}_{-8.1}^{+15.1}\times {10}^{60}$ erg for our low-z sample and ${E}_{\mathrm{therm},\mathrm{grav}}={15.9}_{-8.5}^{+15.9}\times {10}^{60}$ erg for our high-z sample. We therefore measure excess non-gravitational energies of ${E}_{\mathrm{therm},\mathrm{feed},\mathrm{dat}}=-{9.5}_{-10.0}^{+16.1}\times {10}^{60}$ erg for low-z, and ${E}_{\mathrm{therm},\mathrm{feed},\mathrm{dat}}=-{8.9}_{-9.7}^{+16.5}\times {10}^{60}$ erg for high-z, both consistent with 0. Plugging these into the theoretical AGN feedback energy equation and solving for ${\epsilon }_{k}$, we get feedback efficiencies of $-{4.1}_{-4.4}^{+7.0}$% for low-z and $-{3.7}_{-4.0}^{+6.8}$% for high-z. These values are very uncertain and consistent with a detected AGN feedback signal of 0. They also do not rule out the suggested and measured ≈5% (e.g., Scannapieco et al. 2005; Ruan et al. 2015; Spacek et al. 2016). We also note that the feedback efficiencies stated in Spacek et al. (2016) are mistakenly off by a factor of 5, and they should be ${7.3}_{-7.8}^{+6.6}$% for low-z and ${6.6}_{-15.3}^{+17.3}$% for high-z.

In this paper we have performed a stacking analysis of the tSZ signal around 4453 massive elliptical galaxies that are promising candidates for containing relic heating due to past episodes of AGN feedback. We split our selected galaxies into two redshifts bins, with 1179 galaxies in our "low-z" bin ($0.5\leqslant z\leqslant 1.0$) and 3274 galaxies in our "high-z" bin ($1.0\leqslant z\leqslant 1.5$). Our initial stacks were dominated by considerable contaminating emission which was much stronger at 220 GHz than at 148 GHz. Since dusty contaminant emission emits at both 148 and 220 GHz, as suggested by Figure 6, the large signals at 220 GHz, where the tSZ effect is expected to be negligible, indicate a corresponding large contaminant signal at 148 GHz, where the tSZ effect causes a decrement. We therefore performed an analysis of the contaminating signal by modeling potential undetected sources and running a ${\chi }^{2}$ probability test on the models. This revealed the underlying tSZ signal, with a 0.9σ significance at low-z and a 1.8σ significance at high-z. Finally, in order to better constrain the stacked contaminating signal, we incorporated high-frequency Planck measurements of a subset of 227 low-z galaxies and 529 high-z galaxies. These results indicated tSZ detections with a 1.0σ significance at low-z and a 1.6σ significance at high-z. The values for each of these analyses are given in Table 3.

The work done here is complementary to the work done in Spacek et al. (2016), which stacked 4318 galaxies in a southern patch of sky using SPT data, while we stacked 4453 galaxies in the Stripe-82 equatorial band using ACT data. Both analyses used similar galaxy selection criteria, though that of Spacek et al. (2016) favored lower-mass, $0.5\leqslant z\leqslant 1.0$ "low-z" galaxies while our selection favored $1.0\leqslant z\leqslant 1.5$ "high-z" galaxies with higher stellar masses. Their most significant low-z and high-z tSZ detections were at $3.6\sigma$ and $0.9\sigma$ levels, respectively, while ours were at $1.0\sigma$ and $1.8\sigma$, respectively. A detailed comparison between the two studies can be seen in Table 4, where we see similar tSZ Y measurements in this work, although we use galaxies with higher masses. A plot comparing these results, as well as results from other previous galaxy tSZ measurements from Planck Collaboration et al. (2013), Greco et al. (2015), and Ruan et al. (2015), is shown in Figure 9. These last three results appear to be significantly higher than the results of this paper. This may be due to several factors, including inherent differences in the measurements due to different galaxy populations. The previous studies focus on lower-redshift, locally brightest galaxies, while this work looks at quiescent elliptical galaxies at significantly higher redshifts. We also perform extensive cuts to avoid clusters and dusty galaxies. There is the additional possibility that we are not completely accounting for and removing the contamination signal in this work despite our best efforts, though this seems unlikely to be the main reason for the discrepancy. Also shown in Figure 9 are lines representing our simple gravitational and AGN feedback heating models given by Equations (5) and (6). These simple models indicate that the types of galaxies and redshifts that we are looking at are expected to produce significantly lower tSZ measurements. The measurements presented here are also unique, with a review of the literature revealing no other similar measurements of the tSZ signal around such massive, high-redshift, quiescent elliptical galaxies. It therefore may not be completely appropriate to make direct comparisons between these measurements and measurements of less massive, lower-redshift, locally brightest galaxies. While Spacek et al. (2016) estimate AGN feedback efficiencies of around ≈7%, close to the suggested 5% (e.g., Scannapieco et al. 2005; Ruan et al. 2015), this work sees an AGN feedback heating signal consistent with 0, with efficiencies of $-{4.1}_{-4.4}^{+7.0}$% for low-z and $-{3.7}_{-4.0}^{+6.8}$% for high-z. It is important to note, however, that we use simple, general models of gravitation and AGN feedback in this paper to estimate the corresponding energies, and that specific, detailed galaxy simulations are needed to draw more precise conclusions from these measurements.

tSZ measurements of galaxies and AGNs are likely to improve significantly in the near future. More data and an additional band at 277 GHz will be released from ACT observations (Dünner et al. 2013), while an upcoming full survey release of SPT data will include a 2500 deg² field using bands at 95, 150, and 220 GHz (Schaffer et al. 2011). These much larger fields with more bands will allow for a much larger set of galaxies to be co-added at more frequencies, vastly improving the signal-to-noise of the measurements and allowing for further constraints on contaminating signals. Separating out such contaminants will also become more effective with future surveys such as those to be carried out by the upgraded ACT telescope (Advanced ACTPol) and the proposed Cerro Chajnantor Atacama Telescope (CCAT).⁵ Another approach to constraining AGN feedback is through deep measurements of smaller samples of galaxies, identified as the most interesting, using large radio telescopes. The Goddard IRAM Superconducting Two Millimeter Camera (GISMO) and the New IRAM KIDs Array (NIKA) are powerful new instruments mounted on the Institute de Radioastronome Millimetrique (IRAM) 30 m telescope⁶ that may prove useful for this purpose. Also promising is the National Radio Astronomy Observatories (NRAO) Green Bank Telescope (GBT), whose Continuum Backend operates at lower frequencies where the tSZ signal is roughly three times larger. Finally, although tSZ observations reveal the total thermal heating around galaxies, they must be complimented by theoretical models and simulations in order to best distinguish between heating due to gravitation, AGN feedback, and other effects. Observations can therefore be combined with tSZ simulations of the same types of objects with the same average parameters (e.g., mass, redshift, age) to produce weighted stacks that are adapted to be as sensitive as possible to the differences between AGN feedback models. The tSZ effect provides a promising tool for future measurements to improve our understanding of AGN feedback and galaxy evolution.

We would like to thank Arthur Kosowsky for helpful discussions. This publication makes use of data products from the WISE, which is a joint project of UCLA, JPL, and Caltech funded by NASA. It also makes use of data produced from SDSS-III, which was funded by the Alfred P. Sloan Foundation, the National Science Foundation (NSF), the U.S. Department of Energy Office of Science, and the SDSS-III Participating Institutions. Finally, we make use of data from the ACT project, which operates in the Parque Astronómico Atacama in northern Chile under the auspices of the Comisión Nacional de Investigación Científica y Tecnológica de Chile (CONICYT), and is funded by the NSF, Princeton University, U. Penn., and a Canada Foundation for Innovation award to UBC. AS, ES, & SC were supported by the NSF under grant AST14-07835.