# Hot Extended Galaxy Halos Around Local L* Galaxies From Sunyaev-Zeldovich
Measurements

Joel N. Bregman Department of Astronomy
University of Michigan
Ann Arbor, MI 48109 USA  Edmund Hodges-Kluck  NASA/GSFC
Mail Code: 662
Greenbelt , MD 20771 USA  Zhijie Qu  Department of Astronomy
University of Michigan
Ann Arbor, MI 48109 USA  Cameron Pratt  Department of Astronomy
University of Michigan
Ann Arbor, MI 48109 USA  Jiang-Tao Li  Department of Astronomy
University of Michigan
Ann Arbor, MI 48109 USA  Yansong Yun  Department of Astronomy
University of Michigan
Ann Arbor, MI 48109 USA

###### Abstract

Most of the baryons in $L^{*}$ galaxies are unaccounted for and are predicted
to lie in hot gaseous halos (T $\sim 10^{6.5}$ K) that may extend beyond
$R_{200}$. A hot gaseous halo will produce a thermal Sunyaev-Zeldovich signal
that is proportional to the product of the gas mass and the mass-weighted
temperature. To best detect this signal, we used a Needlet Independent Linear
Combination all-sky Planck map that we produced from the most recent Planck
data release, also incorporating WMAP data. The sample is 12 L* spiral
galaxies with distances of $3-10$ Mpc, which are spatially resolved so that
contamination from the optical galaxy can be excluded. One galaxy, NGC 891,
has a particularly strong SZ signal, and when excluding it, the stack of 11
galaxies is detected at about $4\sigma$ (declining with radius) and is
extended to at least 250 kpc ($\approx R_{200}$) at $>\,$99% confidence. The
gas mass within a spherical volume to a radius of 250 kpc is $9.8\pm 2.8\times
10^{10}$ M⊙, for $T_{avg}=3\times 10^{6}$ K. This is about 30% of the
predicted baryon content of the average galaxy ($3.1\times 10^{11}$ M⊙), and
about equal to the mass of stars, disk gas, and warm halo gas. The remaining
missing baryons ($\approx 1.4\times 10^{11}$ M⊙, 40-50% of the total baryon
content) are likely to be hot and extend to the $400-500$ kpc volume, if not
beyond. The result is higher than predictions, but within the uncertainties.

Spiral galaxies, Sunyaev-Zeldovich effect, Circumgalactic medium

##  1 Introduction

One of the striking realizations of the past 30 years is that galaxies are
missing most of their baryons (McGaugh et al., 2010; Dai et al., 2010). This
comes about because we know the mass of the dark matter halo from the rotation
curves of galaxies. The rotation curves do not extend to the virial radius (we
use R200 henceforth, as it is similar to the virial radius), but a number of
rotation curves extend to radii where dark matter dominates the potential. The
dark matter distribution is understood from models (Salucci, 2019), so we can
determine the halo mass, Mh(R200). The cosmic baryon fraction is 15.7% (Planck
Collaboration et al., 2018), corresponding to a dark matter to baryon ratio of
5.37, so the associated baryonic mass is directly predicted.

An accounting of the mass within the optical galaxy requires one measure the
stellar mass of stars, plus corrections for no longer visible stellar remnants
(e.g., Loewenstein 2013). To this is added the mass of the various ISM phases,
such as the HI, HII and H2 disks of spirals, or the hot atmosphere for early-
type galaxies. The ratio of the sum of these baryonic components to the
predicted baryonic mass is $\sim 30$% for an L* galaxy and can be as little as
a few percent for less massive galaxies (McGaugh et al., 2010). This
constitutes the missing baryon problem for galaxies, which has motivated a
number of observational and theoretical studies.

Theoretical models provide context for this problem and show that it is not
just an accounting curiosity, but a crucial clue to understanding galaxy
formation and evolution. The simplest model for the accretion of baryons were
spherical in nature, which led to an accretion shock at large radius,
producing a hot halo at $\sim$Tvirial that might contain significant mass (Mo
et al., 2010). This picture is modified by radiative losses, leading to
significant cooling in the inner region (50-100 kpc) while the outer part of
the halo has a cooling time greater than the Hubble time (Oppenheimer et al.,
2018; Hopkins et al., 2018; Nelson et al., 2018; Sokołowska et al., 2018).

These numerical simulations show that this model is modified in several major
ways for galaxies near L* (e.g., Sokołowska et al. 2018; Lochhaas et al. 2020;
Pandya et al. 2021). The accretion of gas is generally highly asymmetric,
flowing inward along low entropy streams that terminate within 0.1R200. The
ensuing star formation leads to supernovae that heat their surroundings,
driving winds from the galaxies, so some of the gas is rendered unbound and
some remains bound as a hot Circumgalactic Medium (CGM; for spiral galaxies
AGN heating is likely less important). Models are in agreement that the CGM is
turbulent and multi-temperature due to this feedback process. However, models
differ significantly as to the mass of the CGM, its temperature distribution,
the fraction of baryons that lie beyond R200, and the rate at which gas is
cooling onto the disk of the galaxy (e.g., Pandya et al. 2020; Davies et al.
2020). The reason for these broad uncertainties is that the observational
constraints on the hot gas properties of the CGM are limited.

Gas at Tvirial for a L* spiral has a temperature typically of $2-3\times
10^{6}$ K, such as in the Milky Way (Henley & Shelton, 2012, 2013) or NGC 891
(Hodges-Kluck et al., 2018). Some less massive galaxies, such as NGC 4631 have
a lower temperature halo (Wang et al., 2001; Yamasaki et al., 2009), as
expected. At these temperatures, a significant fraction of the thermal
emission falls below the Galactic absorption threshold (0.2-0.3 keV), so these
halos can be X-ray faint. Also, the emission from external galaxies occurs in
the energy band where the Galactic X-ray background is highest (mostly due to
local hot gas processes), adding to the challenge of observing such objects.

Despite these challenges, diffuse X-ray emission is commonly detected around
nearby edge-on L* spirals (Tüllmann et al., 2006; Li & Wang, 2013). This X-ray
gas is found to extend perpendicularly to the disk, and can be traced to about
5-10 kpc above the plane (e.g., Hodges-Kluck et al. 2018). Around the Milky
Way, such hot gas is detected in both emission and absorption, where it is
dominated by material within 50 kpc, about $0.2R_{200}$ (Miller & Bregman,
2015).

Another approach to detecting the extended hot gas halo is through the thermal
Sunyaev-Zel’dovich (SZ) effect(Sunyaev & Zeldovich, 1980), whereby the CMB
photons are scattered by the hot electrons, thus distorting the spectrum.
Previous works have used ensembles of galaxies to identify a signal, using
Planck data. The effort by the Planck Collaboration et al. (2013) and Greco et
al. (2015) stacked $\sim 10^{4}$ galaxies per bin and detected a signal down
to about logM∗ = 11.3. These are massive galaxies (generally early-type) and
the stacking procedure sometimes includes a signal from the larger scale
environment in which they are found (e.g., galaxy clusters). Another stack of
slightly less massive galaxies ($10^{12.6}\lesssim M_{500}$ M⊙ $\lesssim
10^{13}$ M⊙) by Singh et al. (2018) at typical redshifts of 0.1-0.14 combine
X-ray and SZ data and infer that $\sim 20-30$% of baryons in these galaxies
are hot and produce the observed signal. Their inferred $Y$ value is less than
$2\sigma$ significance in most M∗ bins. The study of Ma et al. (2015) cross-
correlated a galaxy sample with weak gravitational lensing (Heymans et al.,
2012) with the Planck data and obtain a signal in broad mass bins (i.e.,
$10^{12}-10^{14}$ M⊙) at about $3\sigma$, with a signal extending beyond the
virial radius. These studies do not provide information about the more
ordinary L* spirals, nor do they address the emission region within R200,
issues that are addressed in this work.

The total SZ signal for a single object, Y, is roughly proportional to
$M_{gas}T_{vir}$, and when $M_{gas}\propto M_{h}$, simple scaling
relationships show that the total SZ signal scales as $Y\propto M_{h}^{5/3}$,
which is why rich galaxy clusters are most easily detected. However, as Taylor
et al. (2003) point out, the signal that one detects depends on the angular
size distance (DA) in the usual D${}_{A}^{-2}$ manner and that nearby lower
mass objects can also have a detectable signal. The signal is measured against
a background that is significant and reduces the S/N of the halos from nearby
galaxies, which are extended at Planck resolution ($10\arcmin$; discussed
below). Despite this obstacle, we find that the SZ signal from nearby L*
galaxies (D${}_{A}<$ 10 Mpc) are at near-detectable levels, and when stacked
together, reveal a clear signal to distances of about $R_{200}$ (250 kpc; we
note that for our sample, DA equals the comoving radial distance, to high
accuracy). This nearby galaxy sample has a major advantage in that the SZ
signal can be resolved at the resolution of Planck, allowing one to avoid
contamination from the optical galaxy and providing spatial information about
the extent of the hot medium.

##  2 Sample Selection

We based our target selection on the predicted strength of the SZ signal, as
given by the $M_{h}^{5/3}D_{A}^{-2}$ scaling relationship. We chose galaxies
close enough that, at Planck resolution, the optical galaxy can easily be
distinguished from the extended SZ halo. The warm dust in the optical galaxy
can cause a signal that contaminates the SZ signal and must be separated. One
can determine $M_{h}$ with a sufficient amount of dynamical data, but that
does not exist uniformly for each of the target galaxies or can be too
uncertain. Instead, we adopted the stellar mass in the K band as a proxy for
$M_{h}$ (Ktot from 2MASS; Skrutskie et al. 2006). The range in the K band
luminosity, or equivalently, $M_{*}$, is modest, with a median of
log$M_{*}=10.83$ and where the quartile values are at 10.76 and 10.90 (0.07
difference with the median; Table 1). The difference between the median and
the maximum or minimum values is 0.19 in log$M_{*}$. For the K band
luminosity, the difference between the highest and lowest quartile boundaries
is 36%.

The distance is an important parameter, but all galaxies have multiple
distance measurements, so we adopted consensus values from the NASA
Extragalactic Database, as given in Table 1. The target selection was the 12
L* spiral galaxies with the largest values of $M_{h}^{5/3}D_{A}^{-2}$, for
distances less than 10 Mpc but beyond the Local Group, and without significant
dust contamination from the midplane (e.g., NGC 6946 is excluded, with a
Galactic latitude of 11.7°). The angular resolution issue led us to choose a
maximum distance of 10 Mpc, at which the resolution of the maps, 10$\arcmin$,
is 29 kpc. This allows us to exclude the galaxy, which can have contaminating
dust emission, while still having enough area to search for a signal from an
extended halo.

Table 1: Sample Galaxies: Basic Properties

No. | Galaxy | Alt Name | RA | DEC | Gal b | Maj Diam | Dist | Ktot | log$M_{*}$ | Scaling
---|---|---|---|---|---|---|---|---|---|---
|  |  | (degrees) | (degrees) | (degrees) | (arcmin) | (Mpc) |  | (M⊙) |
1 | NGC 253 |  | 11.8929 | -25.2922 | -87.97 | 37.20 | 3.22 | -23.75 | 10.81 | 14.85
2 | NGC 891 |  | 35.6367 | 42.3467 | -17.42 | 13.18 | 9.52 | -23.92 | 10.88 | 14.02
3 | NGC 1291 |  | 49.3275 | -41.1081 | -57.04 | 14.45 | 8.60 | -24.02 | 10.92 | 14.18
4 | NGC 3031 | M81 | 148.8882 | 69.0653 | 40.90 | 28.20 | 3.67 | -23.97 | 10.90 | 14.88
5 | NGC 3627 | M66 | 170.0624 | 12.9915 | 64.42 | 9.10 | 9.59 | -23.98 | 10.90 | 14.06
6 | NGC 4258 | M106 | 184.7396 | 47.3040 | 68.84 | 18.60 | 7.28 | -23.83 | 10.84 | 14.19
7 | NGC 4736 | M94 | 192.7233 | 41.1194 | 76.01 | 15.14 | 5.11 | -23.41 | 10.67 | 14.22
8 | NGC 4826 | M64 | 194.1818 | 21.6830 | 84.42 | 13.80 | 5.40 | -23.29 | 10.62 | 14.09
9 | NGC 5055 | M63 | 198.9554 | 42.0292 | 74.29 | 16.22 | 7.72 | -23.78 | 10.82 | 14.11
10 | NGC 5194 | M51 | 202.4696 | 47.2344 | 68.52 | 15.85 | 7.19 | -23.64 | 10.76 | 14.08
11 | NGC 5236 | M83 | 204.2540 | -29.8654 | 31.97 | 18.62 | 6.42 | -24.22 | 11.00 | 14.56
12 | NGC 5457 | M101 | 210.8023 | 54.3490 | 59.77 | 30.20 | 6.86 | -23.63 | 10.76 | 14.11

Note. — In the last column (11), Scaling, is, to a constant,
log($M_{*}^{5/3}D^{-2}$), which was used in the sample selection process,
along with the distance being closer than 10 Mpc.

##  3 Map Selection

An all-sky map for the thermal SZ effect can be constructed from the nine
Planck channel maps, in concert with extracting the CMB, which has a known
spectrum (Planck Collaboration et al., 2016a). The thermal SZ effect produces
a distortion in the CMB with a known spectral shape but where the amplitude
and angular size of the distortion is object-dependent. This distortion causes
a negative signal for $\nu<$ 217 GHz and a positive signal for $\nu>$ 217 GHz,
and it is for this reason that there is a largely SZ-free band centered at 217
GHz. The noise increases at low frequencies due to Galactic free-free
emission, synchrotron emission, and spinning dust emission. At high
frequencies, contamination by dust and detector properties lead to larger
noise. Therefore, much of the weight of the SZ signal is given by the 100-353
GHz bands.

From the Planck data sets, thermal SZ maps, covering most of the sky were
produced and made public. The Needlet Independent Linear Combination y map
(NILC; Delabrouille et al. 2009) uses the Internal Linear Combination (ILC) of
maps to fit for the CMB plus other features, which can include the thermal SZ
and contamination by dust or radio sources. To this, 10 needlet function
(window functions) are added, which are filters in multipole space to capture
(and exclude) most of the contamination on various angular scales. The maps
are smoothed to the same $10\arcmin$ resolution and they use the full set of
high frequency maps (100 - 857 GHz) and the lower frequency data (30-70 GHz)
but only for the large scale multipole space.

These maps were constructed so that the SZ signal is net positive relative to
the local background and that one can treat the NILC product as an image for
the purpose of measuring surface brightness. For our initial studies, we used
this public map rather than the public map produced by the Modified Internal
Linear Combination Algorithm (MILCA; Hurier et al. 2013), which, in our tests,
had significantly higher excess kurtosis (i.e., longer winged tails in the
pixel distribution), leading to more outliers and possible signal
contamination.

These two public SZ maps were constructed from the Planck Data Release 2 (PR2)
(Planck Collaboration et al., 2016b). However, there were significant
improvements in the processing of the Planck data in the two subsequent data
releases, Data Release 3 (PR3) (Planck Collaboration et al., 2020a) and Data
Release 4 (PR4) (Planck Collaboration et al., 2020b). For PR3, the
improvements include better calibration and pipeline processing, along with
changes in the map-making algorithms. However, PR3 contain striping residuals,
which arise in an effort to remove the thermal fluctuations and far sidelobes
as the satellite scans the sky. Striping residuals were largely eliminated in
PR4, which had less systematic noise on all angular scales when compared to
previous releases.

We used the PR4 data, along with the WMAP 9 year data (Bennett et al., 2013)
to construct SZ maps along the lines of the NILC approach, which is more fully
discussed in Pratt et al. (2021). There are a number of choices that need to
be made for the needlet window functions and Gaussian weights, so there is not
a unique set of choices for all purposes. For example, in the public NILC PR2
map, features on the scale of a few degrees was suppressed below its true
value, due to the ILC bias and striping. One cannot account for the ILC bias
in the public maps due to limitations in the information provided regarding
the window functions. In the maps described by Pratt et al. (2021), the SZ
signals on the scales relevant to nearby galaxies are not suppressed, which
was tested by comparing input and extracted signals (below).

##  4 Data Extraction and Contamination

The process of constructing an all-sky SZ map leads to a map in which most
regions away from the Galactic plane are close to a net zero value. However,
on an angular scale in which we extract an image (a few degrees), the local
value of the background is statistically different than zero. The mean
background value can be either positive or negative, so we use a local
background for signal extraction. We exclude strong positive or negative
sources that appear in the field and when we identify these sources, they are
either clusters of galaxies (often in the Planck catalog of SZ sources; Planck
Collaboration et al. 2016a), strongly IR-bright galaxies (easily identified in
the 857 GHz maps), or radio-loud AGNs (also evident in the 100 GHz Planck map,
as well as in other catalogs). When we identify the sources by eye, they are
typically point sources with S/N $>$ 5\. Automated source identification
algorithms (e.g., SExtractor; Bertin & Arnouts 1996) yielded nearly identical
source lists. This typically led to $\sim 10$ excluded regions per 50 square
degrees, as seen in Figure 2. Removing subsequently weaker sources did not
change the outcome of the SZ signal for a target galaxy.

The size of the exclusion regions around contaminating point sources were
circles of 20′ or 30′, depending on source brightness. Beyond these exclusion
regions, the contamination by a point source to the surface brightness of the
target galaxy was well below the signal being sought or the secular
fluctuations in the background. Also, within $R_{200}$ of a sample galaxy, we
excluded other galaxies above 0.2L* (e.g., NGC 3628, near NGC 3627).

Extended sources are rare but sometimes occur in our fields; these are nearby
rich galaxy clusters that are in the Planck SZ Catalog (Planck Collaboration
et al., 2016a), the Abell Cluster catalog (Abell et al., 1989), or in a a
meta-catalogue of X-ray detected clusters of galaxies (MCXC, Piffaretti et al.
2011). In such cases, we would fit annular profiles to the SZ data to
determine the radii at which the surface brightness contribution fell below
that of the background fluctuations, which was typically at 1-2 Mpc. There
were many cases where the exclusion radius could not be determined from the SZ
data, due to signal weakness. In this case, if $R_{200}$ was reliably
determined from the X-ray data, we excluded the region within $R_{200}$.
However, there were a number of clusters identified in the MCXC catalog with
weak SZ signals and for which $R_{200}$ was not reliably estimated from the
X-ray data, due to few X-ray source counts. In those cases, we excluded a
region of radius 1 Mpc.

Radio emission, especially from an AGN, can also contaminate the SZ map and
almost always produce a negative signal. These are not usually important
within the galaxies in our sample, which are generally late-type galaxies with
weak radio emission. We excluded all bright sources identified in the 143 GHz
band and the 353 GHz band (Planck Collaboration et al., 2017), many of which
were similar. An example of a galaxy with strong radio contamination is the
early-type system NGC 5128 (Cen A), which was not a candidate because it is
not a spiral galaxy. In addition, we excluded point like objects, positive or
negative, that were not otherwise identified. These outliers were more than
about $3\sigma$ from the local field mean and are likely due to artifacts of
the NILC method from minimizing the local covariance. Typical examples of
target regions, the annuli, and the exclusion regions are shown in Figure 1
and Figure 2.

The mean SZ surface brightness values, $y$, were obtained in annuli defined in
both angular and physical scales. When using angular scales, the width of a
region was 10$\arcmin$, extending to at least 500 kpc ($4.1^{\circ}$ at D = 7
Mpc, the median of the sample), about $2R_{vir}$ for an L* galaxy.

Figure 1: Left panel: The SZ map for NGC 4258 with annular regions and
excluded regions, where the annular regions extend to a radius of 500 kpc
(236$\arcmin$ ), with significant negative and positive regions excluded.
Among the excluded objects are Abell 1550, NGC 4220 (labeled) and some X-ray
discovered clusters (MCXC designation). Unlabeled positive brightness targets
are possible cluster candidates and negative point sources are often AGNs. The
elliptical HI region around NGC 4258 is also excluded ($13\arcmin\times
27\arcmin$ in semi-major axis).  Figure 2: Right panel: The SZ map for NGC
5194 with annular regions and excluded regions as described in Fig. 1; the
annular regions extend to a radius of 500 kpc (239$\arcmin$). Known QSOs and
galaxy clusters are labeled, while other bright excluded regions may be
cluster candidates. The exclusion region around the optical galaxy
(20$\arcmin$ diameter), also contains NGC 5195.

###  4.1 Excluding Contamination from the Optical Galaxy and Gaseous Disk

We studied the contaminating effect on the SZ maps by dust emission from the
galaxy and central point source emission. Low mass galaxies were chosen for
this study as the SZ signal from a hot halo would be well below our detection
threshold. One of the best cases is M33, with a dynamical mass about five
times lower than our target galaxies (Corbelli et al., 2014), the anticipated
net $Y$ signal would be about 14 times lower than for the L* galaxies in our
sample, assuming that the SZ signal scales as $M_{h}^{5/3}$.

Differences in the temperature and size distributions of the dust lead to
different spectral energy distributions. Depending on these properties, dust
emission can cause either a negative or positive signal in the SZ map. This is
because the SZ map is a linear combination of maps while the shape of the dust
emission varies as a function of temperature, particle size distribution, and
composition. We found that when N(HI) $<1\times 10^{20}$ cm-2, dust emission
is quite weak, such as in M33 (Figure 3). This figure shows the dust
contamination of the SZ signal has become unimportant at $R_{optical}$ (8.8
kpc for the major axis), which is 15% smaller than the radius at which N(HI)
$<1\times 10^{20}$ cm-2 (10.6 kpc). This is consistent with sharp declines in
dust emission around spiral galaxies, such as NGC 891 (Bocchio et al., 2016),
or the exponential decline in dust emission in the combined signal of the
spiral galaxies in the Herschel Reference Survey (Smith et al., 2016a); this
exponential decline, is shown in Figure 3. We note that when N(HI) $<1\times
10^{20}$ cm-2 in the outer parts of galaxies, star formation becomes
particularly weak (Elmegreen & Hunter, 2017), so thermal emission by dust
becomes faint.

For the spirals in our sample, the radius at which N(HI) $=1\times 10^{20}$
cm-2 occurs is typically 20 kpc, although we obtained individual values for
all sample galaxies. For seven galaxies, the HI exclusion region lay within
the central 10′ region, which is always excluded in our SZ analysis; these
include NGC 891, NGC 1291, NGC 3627, NGC 4736, NGC 4826, NGC 5055, and NGC
5194 (M51a and M51b are separated only by 4.4′, so they fit within the central
exclusion region). Of the remaining five galaxies, the HI exclusion regions
were a 22.3′ circle for NGC 253, a $10\arcmin\times 18\arcmin$ elliptical
region (semi-major axis) for NGC 4258, a $10\arcmin\times 15\arcmin$
elliptical region for NGC 5236 (M83), and a $18\arcmin\times 20\arcmin$
elliptical region for NGC 5457 (M101). A more complicated set of exclusion
regions were applied to NGC 3031 (M81), as HI extends irregularly around it.
All HI zero moment maps are available from public databases and articles
(Boomsma et al., 2005; Walter et al., 2008; Haan et al., 2008; de Blok et al.,
2018).

Figure 3:  The SZ surface brightness measurement from the NILC map vs
projected radius shows a dust signal interior to $30\arcmin$ (7.6 kpc) of the
center of M33. The change in the SZ surface brightness $y$ from negative
(first point) to positive (second, third points) is likely due to a decline in
the dust temperature plus data processing procedures in map making. The
optical radius marker is along the major axis and is the same value for the
POSS1 103a-O (Nilson, 1973) and the RC3 D0 (blue; de Vaucouleurs et al. 1995);
$R_{vir}\approx 120$ kpc. The HI column density distribution (scale on
rightKam et al. (2017)) is more extended than the CO surface brightness
(scaled to fit the figure (Druard et al., 2014)), which is typical of spirals.
Based on the NILC $y$ values, there is no measurable contamination to the null
SZ surface brightness beyond N(HI) $=1\times 10^{20}$ cm-2, so we exclude
those interior regions in our galaxy sample. A model for dust emission in
spiral galaxies (red line; Smith et al. 2016a), fitted to the NILC $y$ values
$R\geq 15^{\prime}$, also shows a sharp radial decline. The SZ signal from a
gaseous halo around M33 is undetectable and is predicted to be $1.5\times
10^{-8}$ at R $=55^{\prime}$.

###  4.2 Notes on Individual Galaxy Fields

This sample comprises well-known systems with optical luminosities near L*.
Notes on individual galaxies, their environment, and the characteristics of
the SZ maps are given below. If no mention is made of the SZ maps, they are
typical and presented no special issues.

1\. NGC 253: This is a Sbc/Sc starburst galaxy that is inclined and located in
the center of the Sculptor Group; it is the nearest galaxy in the sample. It
is the largest field ($2R_{200}=8.9^{\circ}$) and is normal, with no unusual
residuals.

2\. NGC 891: An edge-on spiral (Sb) that is similar to the Milky Way and lies
in a modest group, whose main member is suggested to be NGC 1023, although NGC
891 has a slightly brighter K band luminosity; they are separated by more than
1 Mpc. NGC 1023, at D = 10.88 Mpc, lies beyond the 10 Mpc limit and is not
included in this survey.

3\. NGC 1291: An early-type spiral (SBa) in a region of typical noise
properties.

4\. NGC 3031 (M 81): This is the second closest spiral, nearly at the same
distance as NGC 253. It is the most massive spiral (Sab) in the M 81 group,
which contains the starburst galaxy M 82. M 82 and other galaxies were
excluded in the SZ extraction, as were the extended HI regions connecting
several group galaxies. Within the 500 kpc radius region (7.7∘), there are
noticeable residuals, presumably from Galactic dust, which were excluded.

5\. NGC 3627 (M 66): This SBb system lies in the eastern cluster of the Leo I
Group of galaxies, a cluster known as the Leo Triplet, which includes NGC
3628, NGC 3623 (M 65), and some less luminous galaxies; M6 6 is the most
luminous. The other galaxies were excluded in the flux extraction, as were two
galaxy clusters in the outer parts of the extraction region.

6\. NGC 4258 (M 106): This Sbc galaxy has a LINER nucleus and elevated star
formation, although there is no unusual central residual in the SZ map. It is
the brightest galaxy in a group with about 20 members (Thilker et al., 2007).

7\. NGC 4736 (M 94): This Sab galaxy has a LINER nucleus and is the brightest
galaxy in the Canes Venatici I Group, which has 13 members, whose velocity
structure has a strong component of Hubble expansion, rather than being a
relaxed group (Karachentsev et al., 2003).

8\. NGC 4826 (M 64): This is a Sb galaxy that has a LINER designation and with
prominent dust features. It is the brightest member of a poor group (Tully et
al., 2008).

9\. NGC 5055 (M 63): This Sbc galaxy is listed as a H II galaxy and with a
LINER designation. It has only one significant nearby galaxy, UGC 8313 and may
be associated with the M 51 group, lying in the outer parts (SW), with a
separation of $\sim 1$ Mpc.

10\. NGC 5194 (M 51): This grand design Sc spiral (with a LINER designation)
lies in a small group with NGC 5195 and UGC 8331 and are both to the SW of the
M 101 group with a projected separation of $\sim 1$ Mpc. They may be part of
the M 101 group.

11\. NGC 5236 (M 83): This Sbc galaxy is located in a loose group of at least
seven significant galaxies and about 1.5 Mpc (in projection) to the north of
NGC 5128. It is often assigned as part of the NGC 5128 group. In the SZ maps,
there are a few background clusters and foreground residuals that were
excluded.

12\. NGC 5457 (M 101): This Sc galaxy lies at the center of a significant
group (Giuricin et al., 2000) with more than ten members (not including the
small M 51 group to the SW).

None of our targets lie in known X-ray emitting galaxy groups.

Figure 4:  The projected SZ surface brightness distribution y (dimensionless)
for the 12 target galaxies with errors. The zero point of the background
varies between target fields and here is defined as the mean flux in the
400-500 kpc region, adjusted to zero (dashed line). The 1$\sigma$ errors shown
here scale approximately as the inverse square root of the solid angle in an
annulus, leading to statistical variations between galaxies. The vertical
scales are the same ($5\times 10^{-7}$) except for NGC 891 and NGC 3627, which
is slightly larger ($6\times 10^{-7}$).

###  4.3 Uncertainties, Bias, and Recovery Tests

These topics were discussed at length in section §4 of Pratt et al. (2021), so
we provide a summary here, which parallels and makes use of that approach.
There are a few ways of estimating the uncertainty in the signal as a function
of annular radius, but using the intrinsic instrumental uncertainties yields
to errors that are unrealistically small. This is because most of the noise
comes from the uncertainties in separating the two stronger signals (CMB and
dust) from the SZ signal, with most of the noise associated with the “dust”
map, which is actually the sum of all signals that are not the CMB or the SZ.

One way of obtaining the uncertainty in annuli is from the dispersion in
annuli between blank fields, which leads to a dependence of
$\sigma_{i}\propto\Omega_{i}^{-0.4}$, where $\Omega$ is the solid angle of
annulus $i$. One might have expected a $\Omega_{i}^{-0.5}$ dependence, but the
modest secular variations of the backgrounds in the SZ maps around the sky
flattens the dependence on solid angle. The other approach is to perform a
bootstrap analysis of multiple objects, such as the objects in our sample.
This leads to a similar but slightly larger set of uncertainties, so to be
conservative, we used these errors.

Another test involves the accuracy by which an injected signal can be
recovered, with the central part excised, which contains the optical galaxy.
Pratt et al. (2021) used simulated galaxy groups, which are about the same
angular size as our target galaxies, but with a stronger signal. The input
signal was put into the Planck (Data Release 4) and WMAP (Year 9) frequency
maps, from which the NILC map is produced. The recovered surface brightness y
as a function of radius is very close to the input signal (Fig. 4 in Pratt et
al. 2021), with the difference being a modest bias term ($<5\%$). When
recovering a stacked signal, there is little difference between the weighted
average and the bootstrap procedure.

##  5 Results

The goal of this program was to determine the SZ surface brightness quantity
$y$, which is defined in the usual way as

$y=\frac{\sigma_{T}k}{m_{e}c^{2}}\int_{-R_{max}}^{R_{max}}\,n_{e}(x)\,T(x)\,dx$
(1)

where $x$ is the line of sight coordinate and $R_{max}$ is the limit on the
path length. With an ideal detector, $R_{max}$ is the entire cosmic path
length, including the monopole SZ signal. However, we subtract a mean local
background, removing a monopole term, and expect that any signal is local to
the galaxy, where $R_{max}$ is several multiples of $R_{200}$, which should be
constrained from radial surface brightness fitting. The total projected signal
for a target, $Y$, is obtained by integrating $y$ radially. We begin by
showing the stacked quantities for the entire sample and then identify that
one galaxy is significantly brighter than the rest, so the stack is also
constituted for the 11 galaxies without the outlier (NGC 891). The radial
shape of the SZ signal informs the radial dependence of the pressure and the
extent of the signal. This can be converted to a gas mass within the
physically relevant radius $R_{200}$, upon adopting a mass-weighted
temperature or temperature profile.

Extraction of the signal from each galaxy was determined to the same outer
radius, 500 kpc, which is approximately $2R_{200}$ (Figure 4). We stacked the
galaxies in several different ways, both with and without weighting and
obtained similar results. The results shown here (Fig. 5) are weighted by the
scaling from Table 1, which is dominated by the inverse square of the distance
($D^{-2}$), as the range in masses is much smaller. Nearly identical results
are obtained for weightings entirely by $D^{-2}$ or by the statistical mean
uncertainty in the signal. The results of the weighted averages and the
bootstrap procedure (using the same weights) also yielded nearly identical
results and we present the weighted bootstrap results for $Y(r)$, as shown in
Figure 5. We expected the 400-500 kpc region to have a null signal, as this
region was used for the zero-point shift of the individual galaxy fields.
However, the signal appear to flatten by about 250-300 kpc, justifying the use
of the 400-500 kpc region to define the background level. An improved approach
to determining the background is given below, where we adopt a parametric
model that is guided by Figure 5.

Figure 5:  The integrated projected SZ signal ($Y)$ of the 12 galaxy sample,
produced with the weighted bootstrap analysis, where the errors are the
statistical uncertainties. Using a $\beta$ model, a value of $\beta$ = 0.6 is
preferred, while the two models only differ in the cutoff at large radius. The
red line has a Gaussian cutoff with $r_{cut}=300$ kpc, while the blue line has
a much sharper cutoff (see text), both being equally good fits.

We fitted the data with a $\beta$ model with a smooth cutoff, defined as

$y(r)=S_{1}(1+(r/r_{c})^{2})^{-3\beta/2+1/2}\times
e^{(-(r/r_{cut})^{d})}+S_{0}.$ (2)

One cutoff has $d=2$, a Gaussian, and with $r_{cut}=300$ kpc, while the other
has a sharper cutoff, with $d=10$ and $r_{cut}=330$ kpc. Despite the
difference between these two cutoff terms, either provide an acceptable fit
(shown in integrated form in Fig. 5), indicating that there is too little
information to determine the shape of the outer cutoff. Henceforth, we adopt
$d=2$ and $r_{cut}=300$ kpc for fitting surface brightness profiles. Our fit
to $y(r)$ leads to a best-fit value of $\beta=0.61\pm 0.06$, so we chose to
fix $\beta=0.6$ in much of our analysis, although we also examine the impact
of other $\beta$ values in the gas mass determination. This value of
$\beta\approx 0.6$ can be due to the sum of the decline in the density and the
temperature, an issue considered below.

Figure 6: The amplitude of the $y$ signal, based on the $\beta$ model as a
function of NGC number. The galaxy NGC 891 is an outlier, but the remaining 11
galaxies are consistent with a single population with a modest internal
dispersion. For stacking, we use the 11 galaxies and examine NGC 891
separately.

The results from a stack are valid when the objects are drawn from the same
population, without a large intrinsic dispersion about the mean. We examine
whether this is satisfied by fitting the projected $\beta$ models to $y$ for
each galaxy, where we fix $\beta=0.6$, $r_{core}=10$ kpc, and adopt the
Gaussian cutoff with $d=2$, $r_{cut}=300$ kpc. The result does not depend on
the value of $r_{core}$, provided it is less than the inner excluded region
(20-30 kpc), and for galaxies, the core radius is usually $<10$ kpc (e.g.,
Miller & Bregman 2015; Goulding et al. 2016). There is a clear outlier in the
SZ signal, where NGC 891 lies $4.6\sigma$ from the median line defined by the
weighted mean of the other 11 galaxies (Figure 6).

Table 2: SZ Fitted Amplitude and S/N  No. | Galaxy | Alt Name | S1 | S/N
---|---|---|---|---
|  |  | (10-7) |
1 | NGC 253 |  | 2.94 | 1.56
2 | NGC 891 |  | 15.0 | 3.56
3 | NGC 1291 |  | 0.46 | 0.13
4 | NGC 3031 | M81 | 1.74 | 1.16
5 | NGC 3627 | M66 | 0 | 0
6 | NGC 4258 | M106 | 0 | 0
7 | NGC 4736 | M94 | 0.70 | 0.24
8 | NGC 4826 | M64 | 6.99 | 2.16
9 | NGC 5055 | M63 | 6.22 | 2.12
10 | NGC 5194 | M51 | 0.17 | 0.06
11 | NGC 5236 | M83 | 4.45 | 1.43
12 | NGC 5457 | M101 | 2.00 | 0.48

The weighted mean of the 11 galaxies is $S_{1}=2.40\pm 0.64\times 10^{-7}$,
where individual fits resulting in $S_{1}<0$ were set to zero. This fit has a
$\chi_{\nu}^{2}=2.08$, which corresponds to a chance occurrence of 2.5%,
suggesting that there is a modest intrinsic dispersion that has not been
accounted for. The $\chi_{\nu}^{2}$ can be reduced to lie within the 95%
confidence bounds by adding to each $S_{1}$ value an intrinsic dispersion of
$0.6\times 10^{-7}$, or 25% of the weighted mean. When including this
intrinsic dispersion, the error weighted value becomes $S_{1}=2.43\pm
0.67\times 10^{-7}$. This addition of an intrinsic dispersion is not large
compared to the mean, so we conclude that, after excluding NGC 891, the
remaining 11 galaxies can fairly be used for stacking.

###  5.1 Individual Galaxy Detection Significance

We determined the significance of each galaxy in the sample in a few different
ways, which included fitting a parametric function (above) as well as using
the Pearson, Spearman, and Kendall tests to assess the presence of a
correlation between $y$ and $r$. Among the correlation tests, we choose to
discuss the Kendall test, which is non-parametric, more robust against
outliers, and usually yields a lower statistical significance, so it is the
more conservative choice. The six objects with the most significant
correlations are NGC 253 (chance probability or p-value of 5.9%), NGC 3031 (p
= 9.2%), NGC 5236 (p = 1.3%), and three objects with p $<$ 0.2%: NGC 891, NGC
4826, and NGC 5055. This gives 50% of the sample (six systems) with p $<$ 10%,
whereas 1.2 were expected. There are three systems with p $<$ 0.2%, where
0.024 were expected. This indicates that the sample is inconsistent with a
null signal and that a few galaxies were individually detected, although with
significant uncertainties.

The same galaxies with the smallest Kendall p-values have the highest
significance in the $\beta$ model fits (Table 2), with some minor differences.
The confidence for the fit is highest for NGC 891 (chance probability is
0.02%, or $3.56\sigma$), followed by NGC 4826 (1.5% or $2.16\sigma$), which
agrees with visual inspection. These are followed by NGC 5055 (1.7% or
$2.12\sigma$) and two galaxies with S/N $\approx 1.5$: NGC 253 (5.9%); and NGC
5236 (7.6%). This parametric fitting leads to lower significances than the
Kendall tests, due to fitting for two free parameters. From these fits, 5/12
galaxies are significant above the 90% confidence, whereas 1.2 would have been
expected. Three galaxies lie above the 98% confidence, compared to 0.024
expected.

###  5.2 Signal from the 11 Galaxy Stack

We calculated the weighted bootstrap results for both $y$ and $Y$, which are
given in Figure 7 and Figure 8. The significance of $Y$ requires no parametric
model, while $y$ is fitted with the beta model, with the significance only
weakly dependent on most model parameters, such as the nature of the cutoff at
large radius.

When inspecting the significance of $Y(r)$, we note that there are two
uncertainties, the statistical uncertainty, and the uncertainty in the zero
point offset of the stack, which is determined from the mean uncertainty of
the summed individual extractions. This uncertainty grows in importance with
radius as $y$ becomes comparable to the offset uncertainty. It may seem that
the signal could become small at large radius, but it is actually telling us
that the signal at large radius is poorly defined. Physically, $Y$ cannot
decrease with radius, as that would require a negative pressure, which is a
physical impossibility. At the minimum, one would expect $Y$ to become flat at
large radii, if hot gas becomes greatly reduced beyond some radius, such as
$R_{200}$. There is likely an upper bound on $Y$, because if it continued to
rise at a rapid rate, the inferred baryon fraction would exceed the cosmic
baryon fraction.

Figure 7:  The stacked signal from 11 galaxies (NGC 891 is excluded), obtained
by the weighted bootstrap method and fitted with a $\beta$ model where
$\beta=0.6$ and a Gaussian cutoff with $r_{cut}$ = 300 kpc (red curve). The
zero point of the signal is given by the black dashed line and the $1\sigma$
error given by the blue band ($S_{0}=0.18\pm 2.08\times 10^{-9}$).

For $Y(r)$, we show both the statistical uncertainty and the error including
the uncertainty in the zero point offset. For just the statistical
uncertainty, the signal is more than 99% significant through 500 kpc. When
including the offset uncertainty, the signal is significant above the 99%
confidence out to 350 kpc, above the 98% confidence out to 430 kpc, and above
95% confidence to 500 kpc. At a projected radius of 250 kpc, the signal is
$Y(r)=4.39\pm 0.93(stat)\pm 0.79(sys)\times 10^{-3}$ kpc2; the sum of the two
errors is $1.23\times 10^{-3}$ kpc2. The presence of a signal is known better
than the absolute value due to the offset in the background. The signal at 250
kpc is 4.66$\sigma$ (statistical) and 3.57 $\sigma$ (total).

Determining the integrated value at 500 kpc is uncertain due to the errors.
Relative to $Y(250\,\text{kpc})$, the value at 500 kpc could be as much as a
factor of two larger. Using the last five data points (400-500 kpc),
$Y(500\,\text{kpc})\approx 1.5\,Y(250\,\text{kpc})$, where the likely bounds
are $Y(250\,\text{kpc})<Y(500\,\text{kpc})<2Y(250\,\text{kpc})$.

The $y(r)$ distribution has a S/N = 4.1, based on the same parametric fit
applied to the individual galaxies. The fit places constraints on the data and
also determines the zero point offset independently. The normalization $S_{1}$
of the stack, $2.40\times 10^{-7}$, is nearly identical to the weighted mean
of the $S_{1}$ values from the individual 11 galaxies (discussed above),
$2.42\times 10^{-7}$. The extent of the signal can be judged by eliminating
inner bins until the confidence falls below some critical value. We find that
fits to $y$, using data in the 150-500 kpc range, still leads to a significant
signal at the 99% confidence level. The significance falls below 90% when only
using data beyond 190 kpc. There are still positive contributions beyond 190
kpc, as there is only one point that lies below the zero point (at 350 kpc)
until reaching 450 kpc.

Figure 8:  The integrated projected SZ signal, $Y(r)$ (left panel), based on
the 11 galaxy sample (not including NGC 891) and produced with a weighted
bootstrap analysis. The black error bars show the statistical error while the
total error (magenta) includes the zero point uncertainty errors. The S/N of
$Y(r)$ (right panel) is for statistical errors only (black circles) and the
total error that includes the background uncertainty (blue triangles).

We also consider the signals from the outlier system, NGC 891, where the
normalization $S_{1}$ is six times larger than the stack of 11 galaxies. The
signal $Y(r)$ rises up steadily until 200 kpc, at which point it flattens to a
value of $Y(r)\approx 2.6\times 10^{-2}$ kpc2, about 5 times larger than for
the stack of 11 galaxies (Figure 9). At 250 kpc, the signal is $Y=2.51\pm
0.31(stat)\pm 0.27(sys)\times 10^{-2}$ kpc2; the sum of the two errors is
$0.41\times 10^{-2}$ kpc2. The rise is steeper than for the stack of 11
galaxies and is more consistent with $\beta=0.5$ (for $d=2$ and $r_{cut}=200$
kpc, $\beta=0.53\pm 0.10$).

Figure 9:  The integrated SZ signal, $Y(r)$ for NGC 891, where the errors and
lines have the same meaning as the above figures. The signal rises to 200 kpc
and is then flattens at about 0.026 kpc2, about six times larger than from the
stack of the other 11 galaxies.

##  6 Determination of Deprojected SZ Signal and Gas Mass

Beginning with the projected SZ signal, we seek to determine the total SZ
signal within a sphere of 250 kpc and also to place bounds on the total $Y$
signal within a sphere of 500 kpc. This deprojection requires knowledge of the
three dimensional distribution of the pressure, which we will express by a
power-law model with some sort of cutoff, which can be very sharp or a more
gradual steepening. Most of the correction from a projected to spherical SZ
signal occurs at radii greater than 250 kpc, where the SZ signal is most
poorly known, so we calculate the correction for a variety of conditions.

One limiting case, which leads to the largest correction, is where we assume a
power-law density distribution for the pressure, $nT$, to 500 kpc (a $\beta$
model), with $\beta=0.5,0.6,0.7$. For these three cases, the ratio of the SZ
signal within a 250 kpc sphere to that projected on the sky is 0.66, 0.72, and
0.77. Reducing the cutoff from 500 kpc to 400 kpc makes the ratios 0.68, 0.74
and 0.78. Another model is where the signal cuts off more slowly, with a
Gaussian cutoff in three-dimensional radius $R$, with a multiplicative term
added to the power-law of the form $exp(-R^{2}/400^{2})$, where $R$ is in kpc.
We note that this is different than the form used above, where the projected
radius $r$ was used. That cutoff, applied to the $\beta=0.5,0.6,0.7$
distributions, leads to ratios of 0.73, 0.78, and 0.82. At fixed $\beta=0.6$,
decreasing this Gaussian cutoff radius from 400 kpc to 300 kpc raises the
ratio to 0.89, while increasing the radius to 500 kpc, decreases the ratio to
0.74. The data seem to be most consistent with $\beta=0.6$ and a Gaussian
cutoff, although other models cannot be ruled out. Here we adopt an
intermediate ratio of $0.76\pm 0.10$, where the $1\sigma$ range encompasses
nearly all likely values. The observed projected signal within 250 kpc for the
stack of 11 galaxies is $Y(r)=4.39\pm 1.23\times 10^{-3}$ kpc2, so the signal
within a sphere of 250 kpc is $3.34\pm 0.94\times 10^{-3}$ kpc2 (this includes
the uncertainty in the conversion ratio).

The relationship between $Y$ and gas mass is given by

$\displaystyle YD_{A}^{2}=\frac{\sigma_{T}k}{m_{e}c^{2}}\int_{0}^{R_{200}}4\pi
R^{2}\,n_{e}(R)\,T(R)\,dR\ $ (3) $\displaystyle=\
\frac{\sigma_{T}k}{m_{e}c^{2}}M_{e}(R_{200})\,T_{avg}$ (4)

where $M_{e}$ is the mass of electrons and $T_{avg}$ is the mass-weighted
average. When expressed in the gas mass $M_{gas}$, it is

$Y=2.36\times 10^{-3}\,\frac{T_{avg}}{3\times
10^{6}\,\textrm{K}}\,\frac{M_{gas}}{10^{11}\textrm{M}_{\odot}}\,\frac{(7\,\textrm{Mpc})^{2}}{D_{A}^{2}}\,\text{kpc}^{2}.$
(5)

The median distance for the galaxies in the sample is 7 Mpc, but we need to
use the weighted distances, which leads to $D_{A}=5.81$ Mpc, and for the above
value of $Y$, one obtains a gas mass of

$M_{gas}=0.98\pm 0.28\times 10^{11}\,\frac{T_{avg}}{3\times
10^{6}\,\textrm{K}}\ \textrm{M}_{\odot}\ .$ (6)

The temperature of the hot halo for galaxies of this mass is known only close
to the galaxy (e.g., within 10 kpc; e.g., Hodges-Kluck et al. 2018), and is
unknown beyond about 50 kpc even in the Milky Way (Bregman et al., 2018), so
an accurate mass-weighted temperature is not available directly from X-ray
data. We can estimate the temperatures from the virial mass, from X-ray
observations, or from simulations. The total galaxy mass (Mh) within the
virial radius is measured from dynamics or abundance matching, and the
expected baryon mass is calculated by multiplying by the cosmic baryon
fraction of 0.158 (Planck Collaboration et al., 2018). For our sample, we used
the abundance matching approach(Moster et al., 2010), which gives
$M_{h}\approx 2\times 10^{12}$ M⊙, and a total baryon mass of $\approx 3\times
10^{11}$ M⊙. This $M_{h}$ is about 50-100% larger than the value inferred for
the Milky Way (Watkins et al., 2019; Karukes et al., 2020), and as the Milky
Way hot halo gas is about $2\times 10^{6}$ K (Bregman et al., 2018), a value
of about $3\times 10^{6}$ K might be expected for our galaxy stack.
Temperature profiles are determined in simulations, which yield a value near
the virial temperature (Oppenheimer et al., 2018; Hopkins et al., 2018; Nelson
et al., 2018; Sokołowska et al., 2018). However, the mass, temperature, and
extent of simulated halos strongly depend on the uncertain feedback
prescriptions (Davies et al., 2020). Here we adopt $T_{avg}=3\times
10^{6}\,\textrm{K}$ and hope that future work will lead to a more accurate
treatment.

###  6.1 Halo Gas Mass of NGC 891

The $Y$ signal around NGC 891 rises up sharply and turns over near 200 kpc
(Figure 9), so foreground-background contamination from beyond 250 kpc is not
expected to be large. With either a sharp cutoff at 200 kpc, or a Gaussian
decline of $exp(-R^{2}/200^{2})$, the signal within a sphere of 250 kpc is
within a few percent of the projected signal (for $\beta=0.5$, which is a good
fit to the data). The signal of $Y(r)=2.57\pm 0.37\times 10^{-2}$ kpc2 and a
distance of 9.52 Mpc, this corresponds to a mass of $M_{gas}=2.10\pm
0.30\times 10^{12}\,(T/3\times 10^{6}\text{K})$ M⊙.

NGC 891 was observed with long XMM-Newton and Chandra observations (Hodges-
Kluck et al., 2018), where the X-ray emission is detected within 15 kpc of the
plane. Near the star forming regions of the disk, the gas temperatures can be
about 0.71 keV ($8.2\times 10^{6}$ K), but at distances from the plane of 2-10
kpc, the temperature is closer to 0.2 keV ($2.3\times 10^{6}$ K), with a
significantly subsolar metallicity of $Z/Z_{\odot}\approx 0.14$. The halo does
not have an abnormally high gas temperature, so the above mass estimate should
be correct. However, this implies that $M_{gas}\approx M_{h}$, which means
that the baryon mass is about six times larger than the cosmic value if the
galaxy contained all the baryons associated with the dark matter. This excess
baryon mass is most likely associated with a galaxy group, although NGC 891 is
not at the center of a rich group (see object notes above).

##  7 Contribution to the Baryon Budget

We consider the major contributions to the baryon budget for the mean galaxy
in the stack of 11 systems, being the stellar mass, the disk gas, the warm
halo gas, and the hot ambient gas producing the SZ signal. The total baryon
mass is $\approx 3\times 10^{11}$ M⊙, which is about a factor of three higher
than the mass of the hot gas creating the SZ signal within 250 kpc ($0.98\pm
0.28\times 10^{11}\,(T_{avg}/3\times 10^{6}\,\textrm{K})\
\textrm{M}_{\odot}$). The combination of stellar and disk gaseous mass is
typical for an L* galaxy, is $\approx 6\times 10^{10}$ M⊙, which is 20% of the
cosmic baryon mass. To this, we need to consider the mass of warm-hot gas in
the halo, which is detected from UV absorption line measurements.

A census of the warm ionized gas (T $<10^{5}$ K) in the halo of the Milky Way
(Zheng et al., 2019; Qu & Bregman, 2019) and Andromeda (Lehner et al., 2020)
(similar in mass to our sample) leads to masses of $\sim 3\times 10^{9}$ M⊙,
insufficient to close the baryon deficit (Bregman et al., 2018). These masses
do not make a significant contribution to the missing baryon issue.

At higher temperatures, near $10^{5.5}$ K, gas may be traced by O VI. In the
Milky Way, this gas has an estimated halo mass component of $3\times 10^{9}$
M⊙ (Qu & Bregman, 2018a). However, it is about an order of magnitude higher in
M31 (Lehner et al., 2020). We wish to model these O VI columns as a function
of radius, so we fit a model to the declining median projected N(O VI) with
increasing radius $r$ of the form N(OVI) = $5.0\times
10^{14}\times(1+(r/100\,\text{kpc})^{2})^{-0.25}$. We note that this projected
distribution, at large radius, would correspond to a three dimensional density
that decreases as $R^{-3/2}$ (if fit with a $\beta=0.5$ model). Our form for
the O VI density is $n_{OVI}(R)=4.7\times
10^{-12}\,(1+(R/100\,\text{kpc})^{2})^{3/4}$ cm-3 . We note that the integral
of the O VI mass diverges at large radii (500 kpc), increasing as $R^{3/2}$,
so there must be an outer limiting radius. For calculating the gas mass
contribution as a function of radius, we assume an ionization fraction of 0.20
for O VI and a metallicity of 0.3 Solar. We assume that this mass associated
with O VI is typical of our stack of galaxies and show the sum of the
different mass components in Figure 10.

Figure 10:  The baryon masses of different components in a L* galaxy,
including the stars (star symbol), cool disk gas, the hot gas that produces
the SZ signal (black line), and the gas at T $\sim 10^{5.3}$ responsible for
the O VI absorption line gas (green; after Lehner et al. 2020); the sum is
given by the blue line. The uncertainties in the SZ gas are given at 250 kpc
and 500 kpc by red vertical lines. The sum of the stars, disk gas, hot gas
from the SZ study, and the OVI bearing gas is shown as the dark blue line,
which approaches the cosmic baryon mass at about 400–500 kpc. The brown box
(lower left) shows the constraint that the hot gas from the SZ signal exceeds
the X-ray limit.

The contributions to the total expected baryon mass, given in Figure 10, shows
that the $M_{gas,SZ}$ equals the sum of the stars plus cold-warm gas at about
170 kpc, and becomes the dominant component at larger radii. At a radius of
250 kpc, the sum of the mass components constitutes about 50-60% of the total
expected baryon mass, which implies that a significant amount of mass lies
beyond the virial radius. The baryon deficit is nearly gone by 400-500 kpc,
given the uncertainties in the measurements. It is possible that the O VI
contribution is greater, if the metallicity or the ionization fraction is
lower than assumed. Without additional information, it is challenging to put
uncertainties on the O VI contribution, so we do not include any.

##  8 Discussion and Conclusions

This work showed that, with NILC SZ maps made from the most recent Planck Data
Release 4, we are able to see a clear SZ signal in a stack of 11 nearby
($D<10$ Mpc) spiral galaxies. One galaxy was excluded from the original stack
of 12 galaxies because it was about 6 times brighter than the stack, but the
signals from the 11 galaxies do not suggest a large internal dispersion. This
implies that the L* galaxies possess a substantial circumgalactic medium at an
inferred temperature near the virial value. When converted to a gas mass at a
constant mean $T_{avg}=3\times 10^{6}$ K, it equals the stellar mass plus cold
disk gas at 170 kpc, it has $M_{gas}(250kpc)=0.98\pm 0.28\times 10^{11}$ M⊙
and probably rises toward larger radii. At a radius of 250 kpc, the sum of all
baryons is about half of the cosmic value, so a significant amount of “missing
baryons” lie at larger radii. In the following, we examine the relationship
between this SZ signal and a variety of predictions and expectations.

###  8.1 Relationship to X-Ray Emission

The combination of these SZ measurements with X-ray measurements further shows
that some of the baryons must lie beyond the virial radius for spiral galaxies
of similar mass (Fig. 10). The X-ray luminosity is proportional to the volume
integrated emission measure, $\int_{0}^{R_{max}}\\!4\pi
r^{2}n_{e}^{2}\mathrm{d}r$, so at fixed gas mass the emission measure rises
for declining volumes. Diffuse X-ray emission around spiral galaxies is rarely
detected beyond 20 kpc, which implies a limiting line of sight emission
measure of about $\sim 1-3\times 10^{-5}$ cm-6 kpc, for typical Chandra and
XMM-Newton observations (Anderson & Bregman, 2010; Li et al., 2018; Bogdán et
al., 2013, 2017). Our fiducial SZ model is close to the X-ray detection limit
at 20 kpc, and a somewhat lower metallicity (0.1 Solar) or a hotter
temperature ($T>3\times 10^{6}$ K) would reduce the emission further. These
surface brightness limits place a maximum value for the extended hot halo
within 250 kpc of $\approx 1\times 10^{11}$ M⊙, for T $=3\times 10^{6}$ K and
0.3 Solar metallicity. This excludes a hot halo within 250 kpc that contains
all missing baryons.

###  8.2 Comparison to Other SZ Studies

The extended nature of the hot halo is consistent with a study of more distant
systems, where the Canada-France-Hawaii Telescope Lensing Survey (Heymans et
al., 2012) was cross-correlated with the Planck SZ data (Ma et al., 2015),
returning a signal at about $3\sigma$. They divided the signal into two mass
bins, the lower one being M${}_{h}=10^{12}-10^{14}$ M⊙, and two radii bins of
$0.01-1\,R_{vir}$ and the other being $1-100\,R_{vir}$. Using both bins, they
return a SZ signal at about $3\sigma$. Their data suggest that about 40% of
the SZ signal in this bin is extended beyond the virial radius, about a
$2\sigma$ result. This is consistent with our finding, although we use a
narrow mass bin at M${}_{h}=10^{12.3}$ M⊙. Also, there is no discussion as to
whether this extended emission is due to the low mass galaxy clusters that
would lie in the mass range of their bin.

Previous investigations with Planck data reported extended hot halos around
stacks of significantly more massive galaxies (Planck Collaboration et al.,
2013; Greco et al., 2015; Singh et al., 2018; Pratt et al., 2021), mostly
elliptical systems about $\sim 10^{2}$ times more distant than our average
galaxy. At masses of logM∗ $\approx 11.3$, Planck Collaboration et al. (2013)
found that $y$ is not point-like and appears to be larger than 2 Mpc in
radius. This is much larger than $R_{200}$ for their galaxies and we suspect
it is related to the contamination by the surrounding cluster. The luminous
early-type galaxies in their study are rare and often lie in rich
environments, such as the center of a galaxy cluster. In contrast, the L*
spiral galaxies used here are much more common, occur in more ordinary
environments, and define the critical “knee” of the galaxy luminosity
function.

###  8.3 SZ Scaling Relationships Relative to Galaxy Clusters and Massive
Galaxies

At large masses, such as galaxy clusters, the gas is nearly virialized value
and nearly all of the baryons are contained within $R_{200}$ (Laganá et al.,
2013) due to the depth of the potential well. The thermal contribution due to
feedback is small compared to the thermal energy from gravitational collapse,
so under these conditions, the mass weighted temperature $T\propto
M_{h}^{2/3}$. This leads to a total signal $Y\propto M_{gas}\,T_{avg}\propto
M_{h}^{5/3}$, so if an object has an SZ signal that lies on this relationship,
one might conclude that the system contains all of its baryons.

There is a problem when extending this conclusion to galaxies, as they are not
self-similar to galaxy clusters. For example, the slopes of the X-ray surface
brightness distributions are shallower in galaxies, with $\beta\approx 0.5$
(Goulding et al., 2016), than in rich clusters, where $\beta\approx 0.7-0.8$
(Vikhlinin et al., 2006). This difference in $\beta$ is likely due to feedback
having a greater effect on the temperature, which flattens the hydrostatic
density distribution. The gas temperature in spiral and elliptical galaxies is
often $\sim$50% higher than the virial temperature (e.g., Bregman et al.
2018), whereas the rich cluster temperatures are close to the virial
temperature and scale as $M^{2/3}$ (e.g., Mantz et al. 2016 and references
therein). Another difference between clusters and galaxies is that the
fractional stellar mass, $M_{*}/M_{h}$, is typically 0.2, such as for this
sample, while in rich clusters, a value of 0.05 is more typical (Laganá et
al., 2013).

A comparison that can be made is from the scaling of the SZ signal from entire
rich galaxy clusters, for which we can assume that 90% of the cosmic value of
baryons lie within a spherical volume bounded by $R_{200}$ (Laganá et al.,
2013) and at the virial temperature. From the SZ signal from our galaxy stack,
about 30% of the baryons are in a hot halo within $R_{200}$, a factor of three
below the value for rich clusters. However, the hot gas temperature in galaxy
halos is often above the virial temperature, so the net result is that one
would expect $Y$ from our galaxy stack to lie about 50% below the scaling
relationship from rich clusters. That is approximately the result we find,
although it is difficult to determine this with precision, given the scatter
and normalization of the scaling relationship and the uncertainties in our
stack $Y$ value. However, one can make more accurate and physically meaningful
$Y-\text{M}_{*}$ relationships from detailed semi-analytic modeling and
numerical simulation.

Figure 11:  The $Y(R_{500})$ signal vs the log of the stellar mass, where the
red symbol is for the stack of nearby spirals in this work (placed at the
common distance $D_{A}=500$ Mpc). The purple circles are the SZ signal
inferred from X-ray observations, mostly of massive spiral galaxies (Li et
al., 2018), while the black squares are the stacked galaxies from Planck
Collaboration et al. (2013), mostly massive early type galaxies in the centers
of galaxy groups and clusters. The blue line is our scaling relationship for
$Y$ for the hot halo gas mass of Qu & Bregman (2018a). This relationship falls
below our stacked value by a factor of three, it passes within a factor of two
of the massive spirals, and lies an order of magnitude below the stacks from
Planck Collaboration et al. (2013).

###  8.4 Comparison to Semi-Analytic Models

There have been a few semianalytical models for gaseous halos applied to the
Milky Way and to galaxies of a range of masses (Faerman et al., 2017; Qu &
Bregman, 2018a, b; Voit et al., 2018; Faerman et al., 2020). These
investigations include calculations of halo gas masses, SZ signals and entropy
distributions. This type of modeling requires that investigators make
assumption for the temperature properties of the halo, as they are not yet
sufficiently constrained by observations. For example, Faerman et al. (2020)
model the Milky Way halo, where they assume an isentropic (constant entropy)
halo, but with a ratio of specific heats $\gamma=4/3$ to account for non-
thermal pressure support (e.g., cosmic rays, magnetic fields). The mass of
their hot gaseous halo is $5.5\times 10^{10}$ M⊙, compared to the value from
the stack, $9.8\pm 2.8\times 10^{10}$ M⊙. This is about a $1.5\sigma$
difference, so not in significant disagreement, especially considering that
the Milky Way galaxy is less massive than our stack of galaxies.

Another work on hot halos is by Qu & Bregman (2018a), who adopt a maximum
temperature, such as 1-2Tvir, which does not change with radius. This is not a
true isothermal models because the gas radiatively cools, so there is a
significant temperature distribution where the cooling time is less than a
Hubble time; the mass-weighted temperature is not independent of radius. The
model was developed for comparison with the columns of ionic species, so the
cooling function is corrected for photoionization and here we use the results
of this TPIE case. At R200, the model gas mass is about half of the value for
the stack, but the better comparison may be the SZ $Y$ value, which we have
calculated here, based on the Qu & Bregman (2018a) work. This work provides
M∗, Mgas and mass-weighted temperatures as a function of Mh, and the results
are shown in Figure 11. This model line lies about a factor of three below the
measured value from our stack, although it does pass within a factor of two of
the massive spiral $Y$ values determined from X-ray data in Figure 11.

It is not surprising that $Y$ from the stacks of Planck Collaboration et al.
(2013) lie above our galaxy scaling, as the environments are extremely
different. Most galaxies in their study lie in galaxy clusters and galaxy
groups, which have hot gas with a mass that is often large compared to the
central galaxy if it were relatively isolated. This pressure of the
cluster/group medium would likely raise the temperature and pressure of the
gas within $R_{500}$ of the central galaxy, thereby biasing $Y_{500}$ to
larger values relative to our relationship for non-cluster galaxies. The
magnitude of the offset would require extensive simulations and is beyond the
scope of this paper.

Voit et al. (2018) infer the CGM pressure profile around galaxies by assuming
that the UV absorption line clouds (Werk et al., 2014) are in pressure
equilibrium with a hot ambient medium. They find that this ambient pressure
can be fit with a $R^{-1.8}$ profile, which is consistent with the pressure
profile of our stack of galaxies.

###  8.5 Comparison to Predicted Signals from Large Scale Structure
Simulations

Hot gaseous halos are a generic feature of simulated galaxies (Oppenheimer et
al., 2018; Hopkins et al., 2018; Nelson et al., 2018; Sokołowska et al.,
2018). These simulations show that a fraction of the baryons form a stable
rotating disk of stars and cool or warm gas in the central $0.1R_{vir}$. Much
of the remaining gas forms a hot halo near the virial temperature that extends
to near (or beyond) the virial radius, but the mass and extent of simulated
halos strongly depend on the uncertain feedback prescriptions (Davies et al.,
2020).

One simulation study, from the Feedback in Realistic Environments program
(FIRE; van de Voort et al. 2016) calculated the various SZ signals and baryon
fractions at $z=0-0.5$, within $R_{500}$, as a function of $M_{500}$. We note
that FIRE includes stellar but not AGN feedback in this work. Their hot baryon
fraction within $R_{500}$ is similar to our value, $0.21\pm 0.06$. However,
the more direct comparison is to their SZ parameter, where our signal lies
toward the upper $1\sigma$ envelope of their simulation results.

A study of halo masses with $EAGLE$ (Davies et al., 2019) has a similar factor
of two spread in the prediction of $Y$ and the gaseous halo baryon fraction at
fixed $M_{h}$. Our hot halo mass (from the SZ signal) at $R_{200}$ is $0.31\pm
0.09$ of the cosmic baryon fraction, indistinguishable from the $0.27$ value
of the model. For $Y$, our value is about a factor of two above their median
value at logM200 = 12.3 but within the dispersion. This could be explained if
there is a temperature difference of a factor of two between the simulations
and our data.

In a subsequent paper, Davies et al. (2020) considered the difference between
EAGLE and IllustrisTNG, simulations that include feedback from both AGNs and
stars. Their baryon fractions for the CGM differ by a factor of 2-3 for
$11.5<\text{logM}_{200}<12$ (IllustrisTNG is higher), although at
logM${}_{200}=12.3$, they happen to be very similar, with baryon fractions for
the CGM $f_{CGM}/(\Omega_{b}/\Omega_{0})\approx 0.25$. These values are
consistent with our stack value of $0.31\pm 0.09$ at R200. Overall, our halo
gas properties and $Y$ value are in agreement with these three simulation
codes.

###  8.6 Future Prospects

The primary limitation in extracting a SZ signal for nearby galaxies (resolved
by a 10$\arcmin$ beam) is not the sensitivity of the instrumentation – it is
the uncertainty introduced by the separation of the Galactic signal. Pratt et
al. (2021) improved upon the production of SZ maps by using the most recent
Planck data release and window functions that do not suppress the SZ signal on
the scale of relevance to this study. We continue to develop better algorithms
to address the separation of the dust contamination, using existing data.
Further improvement should be possible by introducing more channels, such as
the probe concept PICO, which is designed to have 21 channels (Hanany et al.,
2019). Another improvement might result from higher angular resolution, as
would be possible with such as CMB-S4 (Abazajian et al., 2019). Higher
resolution could help to more accurately remove contaminating features. It
might be possible to improve results by combining the Planck and WMAP data
with ground based programs, such as SPT (Carlstrom et al., 2011) and ACT
(Thornton et al., 2016), and we hope to explore that in the future.

Until new instruments are possible, we can expand the size of the stack of
resolved galaxies by going to greater distances. The disadvantage of using
more distant galaxies is that a larger central exclusion region is required to
avoid contamination by the optical galaxy signal. The goal would be to improve
the signal in the $0.5-2R_{200}$ region so that we can understand the
properties of the SZ signal in this outer region. The practical distance limit
is about 50 Mpc, where the 10′ beam has a linear size of 145 kpc.

For the conversion of the SZ signal into a mass, the primary uncertainty is
the temperature profile. The temperature profile of gas at $\sim 10^{6.5}$ K
can only be obtained from X-ray observations. Such halos should be detectable
in absorption lines from the resonance lines of O VII, O VIII, C VI, and Ne
IX, to name a few (Bregman et al., 2015). Future missions, such as Athena
(Barcons et al., 2017), Arcus (Smith et al., 2016b), or Lynx (Gaskin et al.,
2018) will have the sensitivity to make major contributions in such absorption
line studies. Absorption lines from such hot halo gas should be detectable to
at least $R_{200}$, and likely beyond.

Emission line studies with future instruments (calorimeters) will probe the
hot halos as well, although to smaller radii, such as $\sim 0.5R_{200}$ for L*
galaxies. The primary missions that could survey galaxies to significant radii
are Athena, Lynx, the Hot Universe Baryon Surveyor (HUBS) Cui et al. (2020),
or Super DIOS (Sato et al., 2020). These emission studies will give us
temperature and metallicity profiles for individual galaxies. When combined
with absorption lines, permit greater insight into the metallicity and the
filling factors. The X-ray observations, combined with SZ studies, will
provide the definitive description of the hot halo gas.

We thank the many individuals who offered guidance and insight, including
Mathieu Remazeilles, Renato Dupke, Chris Miller, Oleg Gnedin, Ben Oppenheimer,
Mike Anderson, and Jeff McMahon. We are grateful for support from NASA through
the Astrophysics Data Analysis Program, awards NNX15AM93G and 80NSSC19K1013.
The work is based on observations obtained with Planck,
http://www.esa.int/Planck, an ESA science mission with instruments and
contributions directly funded by ESA Member States, NASA, and Canada. This
research has made use of data and software provided by the High Energy
Astrophysics Science Archive Research Center (HEASARC), which is a service of
the Astrophysics Science Division at NASA/GSFC and the High Energy
Astrophysics Division of the Smithsonian Astrophysical Observatory. This
research has made use of the NASA/IPAC Extragalactic Database (NED), which is
operated by the Jet Propulsion Laboratory, California Institute of Technology,
under contract with NASA.

## References

  * Abazajian et al. (2019) Abazajian, K., Addison, G., Adshead, P., et al. 2019, arXiv e-prints, arXiv:1908.01062.  https://arxiv.org/abs/1908.01062
  * Abell et al. (1989) Abell, G. O., Corwin, Harold G., J., & Olowin, R. P. 1989, ApJS, 70, 1, doi: 10.1086/191333
  * Anderson & Bregman (2010) Anderson, M. E., & Bregman, J. N. 2010, ApJ, 714, 320, doi: 10.1088/0004-637X/714/1/320
  * Barcons et al. (2017) Barcons, X., Barret, D., Decourchelle, A., et al. 2017, Astronomische Nachrichten, 338, 153, doi: 10.1002/asna.201713323
  * Bennett et al. (2013) Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20, doi: 10.1088/0067-0049/208/2/20
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393, doi: 10.1051/aas:1996164
  * Bocchio et al. (2016) Bocchio, M., Bianchi, S., Hunt, L. K., & Schneider, R. 2016, A&A, 586, A8, doi: 10.1051/0004-6361/201526950
  * Bogdán et al. (2017) Bogdán, Á., Bourdin, H., Forman, W. R., et al. 2017, ApJ, 850, 98, doi: 10.3847/1538-4357/aa9523
  * Bogdán et al. (2013) Bogdán, Á., Forman, W. R., Kraft, R. P., & Jones, C. 2013, ApJ, 772, 98, doi: 10.1088/0004-637X/772/2/98
  * Boomsma et al. (2005) Boomsma, R., Oosterloo, T. A., Fraternali, F., van der Hulst, J. M., & Sancisi, R. 2005, A&A, 431, 65, doi: 10.1051/0004-6361:20041715
  * Bregman et al. (2015) Bregman, J. N., Alves, G. C., Miller, M. J., & Hodges-Kluck, E. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 045003, doi: 10.1117/1.JATIS.1.4.045003
  * Bregman et al. (2018) Bregman, J. N., Anderson, M. E., Miller, M. J., et al. 2018, ApJ, 862, 3, doi: 10.3847/1538-4357/aacafe
  * Carlstrom et al. (2011) Carlstrom, J. E., Ade, P. A. R., Aird, K. A., et al. 2011, PASP, 123, 568, doi: 10.1086/659879
  * Corbelli et al. (2014) Corbelli, E., Thilker, D., Zibetti, S., Giovanardi, C., & Salucci, P. 2014, A&A, 572, A23, doi: 10.1051/0004-6361/201424033
  * Cui et al. (2020) Cui, W., Bregman, J. N., Bruijn, M. P., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11444, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 114442S, doi: 10.1117/12.2560871
  * Dai et al. (2010) Dai, X., Bregman, J. N., Kochanek, C. S., & Rasia, E. 2010, ApJ, 719, 119, doi: 10.1088/0004-637X/719/1/119
  * Davies et al. (2019) Davies, J. J., Crain, R. A., McCarthy, I. G., et al. 2019, MNRAS, 485, 3783, doi: 10.1093/mnras/stz635
  * Davies et al. (2020) Davies, J. J., Crain, R. A., Oppenheimer, B. D., & Schaye, J. 2020, MNRAS, 491, 4462, doi: 10.1093/mnras/stz3201
  * de Blok et al. (2018) de Blok, W. J. G., Walter, F., Ferguson, A. M. N., et al. 2018, ApJ, 865, 26, doi: 10.3847/1538-4357/aad557
  * de Vaucouleurs et al. (1995) de Vaucouleurs, G., de Vaucouleurs, A., Corwin, H. G., et al. 1995, VizieR Online Data Catalog, VII/155
  * Delabrouille et al. (2009) Delabrouille, J., Cardoso, J.-F., Le Jeune, M., et al. 2009, A&A, 493, 835, doi: 10.1051/0004-6361:200810514
  * Druard et al. (2014) Druard, C., Braine, J., Schuster, K. F., et al. 2014, A&A, 567, A118, doi: 10.1051/0004-6361/201423682
  * Elmegreen & Hunter (2017) Elmegreen, B. G., & Hunter, D. A. 2017, in Astrophysics and Space Science Library, Vol. 434, Outskirts of Galaxies, ed. J. H. Knapen, J. C. Lee, & A. Gil de Paz, 115, doi: 10.1007/978-3-319-56570-5_4
  * Faerman et al. (2017) Faerman, Y., Sternberg, A., & McKee, C. F. 2017, ApJ, 835, 52, doi: 10.3847/1538-4357/835/1/52
  * Faerman et al. (2020) —. 2020, ApJ, 893, 82, doi: 10.3847/1538-4357/ab7ffc
  * Gaskin et al. (2018) Gaskin, J. A., Dominguez, A., Gelmis, K., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10699, Space Telescopes and Instrumentation 2018: Ultraviolet to Gamma Ray, 106990N, doi: 10.1117/12.2314149
  * Giuricin et al. (2000) Giuricin, G., Marinoni, C., Ceriani, L., & Pisani, A. 2000, ApJ, 543, 178, doi: 10.1086/317070
  * Goulding et al. (2016) Goulding, A. D., Greene, J. E., Ma, C.-P., et al. 2016, ApJ, 826, 167, doi: 10.3847/0004-637X/826/2/167
  * Greco et al. (2015) Greco, J. P., Hill, J. C., Spergel, D. N., & Battaglia, N. 2015, ApJ, 808, 151, doi: 10.1088/0004-637X/808/2/151
  * Haan et al. (2008) Haan, S., Schinnerer, E., Mundell, C. G., García-Burillo, S., & Combes, F. 2008, AJ, 135, 232, doi: 10.1088/0004-6256/135/1/232
  * Hanany et al. (2019) Hanany, S., Alvarez, M., Artis, E., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 194.  https://arxiv.org/abs/1908.07495
  * Henley & Shelton (2012) Henley, D. B., & Shelton, R. L. 2012, ApJS, 202, 14, doi: 10.1088/0067-0049/202/2/14
  * Henley & Shelton (2013) —. 2013, ApJ, 773, 92, doi: 10.1088/0004-637X/773/2/92
  * Heymans et al. (2012) Heymans, C., Van Waerbeke, L., Miller, L., et al. 2012, MNRAS, 427, 146, doi: 10.1111/j.1365-2966.2012.21952.x
  * Hodges-Kluck et al. (2018) Hodges-Kluck, E. J., Bregman, J. N., & Li, J.-t. 2018, ApJ, 866, 126, doi: 10.3847/1538-4357/aae38a
  * Hopkins et al. (2018) Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, 480, 800, doi: 10.1093/mnras/sty1690
  * Hurier et al. (2013) Hurier, G., Macías-Pérez, J. F., & Hildebrandt, S. 2013, A&A, 558, A118, doi: 10.1051/0004-6361/201321891
  * Kam et al. (2017) Kam, S. Z., Carignan, C., Chemin, L., et al. 2017, AJ, 154, 41, doi: 10.3847/1538-3881/aa79f3
  * Karachentsev et al. (2003) Karachentsev, I. D., Sharina, M. E., Dolphin, A. E., et al. 2003, A&A, 398, 467, doi: 10.1051/0004-6361:20021598
  * Karukes et al. (2020) Karukes, E. V., Benito, M., Iocco, F., Trotta, R., & Geringer-Sameth, A. 2020, J. Cosmology Astropart. Phys, 2020, 033, doi: 10.1088/1475-7516/2020/05/033
  * Laganá et al. (2013) Laganá, T. F., Martinet, N., Durret, F., et al. 2013, A&A, 555, A66, doi: 10.1051/0004-6361/201220423
  * Lehner et al. (2020) Lehner, N., Berek, S. C., Howk, J. C., et al. 2020, ApJ, 900, 9, doi: 10.3847/1538-4357/aba49c
  * Li et al. (2018) Li, J.-T., Bregman, J. N., Wang, Q. D., Crain, R. A., & Anderson, M. E. 2018, ApJ, 855, L24, doi: 10.3847/2041-8213/aab2af
  * Li & Wang (2013) Li, J.-T., & Wang, Q. D. 2013, MNRAS, 428, 2085, doi: 10.1093/mnras/sts183
  * Lochhaas et al. (2020) Lochhaas, C., Bryan, G. L., Li, Y., Li, M., & Fielding, D. 2020, MNRAS, 493, 1461, doi: 10.1093/mnras/staa358
  * Loewenstein (2013) Loewenstein, M. 2013, ApJ, 773, 52, doi: 10.1088/0004-637X/773/1/52
  * Ma et al. (2015) Ma, Y.-Z., Van Waerbeke, L., Hinshaw, G., et al. 2015, J. Cosmology Astropart. Phys, 2015, 046, doi: 10.1088/1475-7516/2015/09/046
  * Mantz et al. (2016) Mantz, A. B., Allen, S. W., Morris, R. G., et al. 2016, MNRAS, 463, 3582, doi: 10.1093/mnras/stw2250
  * McGaugh et al. (2010) McGaugh, S. S., Schombert, J. M., de Blok, W. J. G., & Zagursky, M. J. 2010, ApJ, 708, L14, doi: 10.1088/2041-8205/708/1/L14
  * Miller & Bregman (2015) Miller, M. J., & Bregman, J. N. 2015, ApJ, 800, 14, doi: 10.1088/0004-637X/800/1/14
  * Mo et al. (2010) Mo, H., van den Bosch, F. C., & White, S. 2010, Galaxy Formation and Evolution (Cambridge University Press, Cambridge, UK)
  * Moster et al. (2010) Moster, B. P., Somerville, R. S., Maulbetsch, C., et al. 2010, ApJ, 710, 903, doi: 10.1088/0004-637X/710/2/903
  * Nelson et al. (2018) Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624, doi: 10.1093/mnras/stx3040
  * Nilson (1973) Nilson, P. 1973, Nova Acta Regiae Soc. Sci. Upsaliensis Ser. V, 0
  * Oppenheimer et al. (2018) Oppenheimer, B. D., Schaye, J., Crain, R. A., Werk, J. K., & Richings, A. J. 2018, MNRAS, 481, 835, doi: 10.1093/mnras/sty2281
  * Pandya et al. (2020) Pandya, V., Somerville, R. S., Anglés-Alcázar, D., et al. 2020, ApJ, 905, 4, doi: 10.3847/1538-4357/abc3c1
  * Pandya et al. (2021) Pandya, V., Fielding, D., Anglés-Alcázar, D., et al. 2021, arXiv e-prints, arXiv:2103.06891.  https://arxiv.org/abs/2103.06891
  * Piffaretti et al. (2011) Piffaretti, R., Arnaud, M., Pratt, G. W., Pointecouteau, E., & Melin, J. B. 2011, A&A, 534, A109, doi: 10.1051/0004-6361/201015377
  * Planck Collaboration et al. (2013) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2013, A&A, 557, A52, doi: 10.1051/0004-6361/201220941
  * Planck Collaboration et al. (2016a) —. 2016a, A&A, 594, A27, doi: 10.1051/0004-6361/201525823
  * Planck Collaboration et al. (2016b) Planck Collaboration, Adam, R., Ade, P. A. R., et al. 2016b, A&A, 594, A1, doi: 10.1051/0004-6361/201527101
  * Planck Collaboration et al. (2017) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2017, VizieR Online Data Catalog, J/A+A/594/A26
  * Planck Collaboration et al. (2018) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2018, arXiv e-prints, arXiv:1807.06209.  https://arxiv.org/abs/1807.06209
  * Planck Collaboration et al. (2020a) —. 2020a, A&A, 641, A1, doi: 10.1051/0004-6361/201833880
  * Planck Collaboration et al. (2020b) Planck Collaboration, Akrami, Y., Andersen, K. J., et al. 2020b, arXiv e-prints, arXiv:2007.04997.  https://arxiv.org/abs/2007.04997
  * Pratt et al. (2021) Pratt, C. T., Qu, Z., & Bregman, J. N. 2021, arXiv e-prints, arXiv:2105.01123.  https://arxiv.org/abs/2105.01123
  * Qu & Bregman (2018a) Qu, Z., & Bregman, J. N. 2018a, ApJ, 856, 5, doi: 10.3847/1538-4357/aaafd4
  * Qu & Bregman (2018b) —. 2018b, ApJ, 862, 23, doi: 10.3847/1538-4357/aaccec
  * Qu & Bregman (2019) —. 2019, ApJ, 880, 89, doi: 10.3847/1538-4357/ab2a0b
  * Salucci (2019) Salucci, P. 2019, A&A Rev., 27, 2, doi: 10.1007/s00159-018-0113-1
  * Sato et al. (2020) Sato, K., Ohashi, T., Ishisaki, Y., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11444, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 114445O, doi: 10.1117/12.2561681
  * Singh et al. (2018) Singh, P., Majumdar, S., Nath, B. B., & Silk, J. 2018, MNRAS, 478, 2909, doi: 10.1093/mnras/sty1276
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163, doi: 10.1086/498708
  * Smith et al. (2016a) Smith, M. W. L., Eales, S. A., De Looze, I., et al. 2016a, MNRAS, 462, 331, doi: 10.1093/mnras/stw1611
  * Smith et al. (2016b) Smith, R. K., Abraham, M. H., Allured, R., et al. 2016b, in Proc. SPIE, Vol. 9905, Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray, 99054M, doi: 10.1117/12.2231778
  * Sokołowska et al. (2018) Sokołowska, A., Babul, A., Mayer, L., Shen, S., & Madau, P. 2018, ApJ, 867, 73, doi: 10.3847/1538-4357/aae43a
  * Sunyaev & Zeldovich (1980) Sunyaev, R. A., & Zeldovich, I. B. 1980, ARA&A, 18, 537, doi: 10.1146/annurev.aa.18.090180.002541
  * Taylor et al. (2003) Taylor, J. E., Moodley, K., & Diego, J. M. 2003, MNRAS, 345, 1127, doi: 10.1046/j.1365-2966.2003.07025.x
  * Thilker et al. (2007) Thilker, D. A., Bianchi, L., Meurer, G., et al. 2007, ApJS, 173, 538, doi: 10.1086/523853
  * Thornton et al. (2016) Thornton, R. J., Ade, P. A. R., Aiola, S., et al. 2016, ApJS, 227, 21, doi: 10.3847/1538-4365/227/2/21
  * Tüllmann et al. (2006) Tüllmann, R., Pietsch, W., Rossa, J., Breitschwerdt, D., & Dettmar, R. J. 2006, A&A, 448, 43, doi: 10.1051/0004-6361:20052936
  * Tully et al. (2008) Tully, R. B., Shaya, E. J., Karachentsev, I. D., et al. 2008, ApJ, 676, 184, doi: 10.1086/527428
  * van de Voort et al. (2016) van de Voort, F., Quataert, E., Hopkins, P. F., et al. 2016, MNRAS, 463, 4533, doi: 10.1093/mnras/stw2322
  * Vikhlinin et al. (2006) Vikhlinin, A., Kravtsov, A., Forman, W., et al. 2006, ApJ, 640, 691, doi: 10.1086/500288
  * Voit et al. (2018) Voit, G. M., Ma, C. P., Greene, J., et al. 2018, ApJ, 853, 78, doi: 10.3847/1538-4357/aaa084
  * Walter et al. (2008) Walter, F., Brinks, E., de Blok, W. J. G., et al. 2008, AJ, 136, 2563, doi: 10.1088/0004-6256/136/6/2563
  * Wang et al. (2001) Wang, Q. D., Immler, S., Walterbos, R., Lauroesch, J. T., & Breitschwerdt, D. 2001, ApJ, 555, L99, doi: 10.1086/323179
  * Watkins et al. (2019) Watkins, L. L., van der Marel, R. P., Sohn, S. T., & Evans, N. W. 2019, ApJ, 873, 118, doi: 10.3847/1538-4357/ab089f
  * Werk et al. (2014) Werk, J. K., Prochaska, J. X., Tumlinson, J., et al. 2014, ApJ, 792, 8, doi: 10.1088/0004-637X/792/1/8
  * Yamasaki et al. (2009) Yamasaki, N. Y., Sato, K., Mitsuishi, I., & Ohashi, T. 2009, PASJ, 61, S291, doi: 10.1093/pasj/61.sp1.S291
  * Zheng et al. (2019) Zheng, Y., Peek, J. E. G., Putman, M. E., & Werk, J. K. 2019, ApJ, 871, 35, doi: 10.3847/1538-4357/aaf6eb