An Empirical Determination of the Dependence of the Circumgalactic Mass Cooling Rate and Feedback Mass Loading Factor on Galactic Stellar Mass

Huanian Zhang; Dennis Zaritsky; Karen Pardos Olsen; Peter Behroozi; Jessica Werk; Robert Kennicutt; Lizhi Xie; Xiaohu Yang; Taotao Fang; Gabriella De Lucia; Michaela Hirschmann; Fabio Fontanot

doi:10.3847/1538-4357/ac0433

1. Introduction

The interstellar medium (ISM) alone cannot sustain the star formation activity of a galaxy for its entire lifespan, thus necessitating the inflow of cool gas from the circumgalactic medium (CGM) to provide additional fuel (Spitzer 1956). Indeed, the COS-GASS survey found a correlation between interstellar and circumgalactic H i, implying a physical connection between the H i disk and the CGM such that the H i disk is nourished by accretion of gas from the CGM (Borthakur et al. 2015). Through observations and simulations, we have now arrived at the conclusion that the baryon cycling between the CGM and central galaxy is the dominant factor driving the observable properties of galaxies (see Tumlinson et al. 2017, for a review). As such, measuring the net gas mass inflow rate is key to our understanding of galaxy evolution.

Hydrodynamic simulations have shed light on the delicate process of precipitating cold gas out of the hot medium in the CGM. Small-box simulations of a region of a galactic wind flow or the CGM can reach the resolution needed to track cooling/precipitation and have been used to study the possible origin of phenomena such as cold molecular clumps found in quasar outflows (e.g., Ferrara & Scannapieco 2016; Gronke & Oh 2018; Liang & Remming 2020). Larger cosmological simulations provide the information required to understand more global phenomena such as the observed O vi content of galactic halos (e.g., Suresh et al. 2017; Oppenheimer et al. 2018; Peeples et al. 2019; van de Voort et al. 2019; Hummels et al. 2019; Nelson et al. 2020). Among others, one of the ultimate goals of these hydrodynamic simulations, combined with more phenomenological treatments, is to determine the inflow rate of cold gas into the central portion of galaxies and the fraction of that gas that must be removed or reheated, so that together they fuel and throttle the process of star formation on galactic scales (Muratov et al. 2015; Hayward & Hopkins 2017; Xie et al. 2017, 2020).

Most observational studies attempt to trace the CGM and establish its basic characteristics (e.g., total mass, radial density profile, metallicity, kinematics). There are relatively few empirical constraints on the cooling rate itself, which is more closely related to the star formation process in galaxies (for one exception to this statement see Voit et al. 2019, which models the edges of cold clouds to estimate the precipitation rate). The reason for the absence of constraints on the cooling rate is that studies of the CGM have come primarily from the study of absorption lines in the spectra of bright background objects (e.g., Chen et al. 2010; Steidel et al. 2010; Bordoloi et al. 2011; Ménard et al. 2011; Zhu & Ménard 2013a, 2013b; Werk et al. 2013, 2014, 2016; Johnson et al. 2013, 2014, 2015, 2017; Croft et al. 2016, 2018; Prochaska et al. 2017; Cai et al. 2017; Chen et al. 2010, 2019; Chen 2017a, 2017b; Lan & Mo 2018; Joshi et al. 2018; Zahedy et al. 2019). However, a growing set of studies is focusing on emission lines. For example, there are now measurements of UV or soft X-ray emission lines in either low-redshift starburst galaxies (Hayes et al. 2016) or high-redshift galaxies (Vikhlinin 2006; Cantalupo et al. 2014; Prescott et al. 2015; Cai et al. 2019).

In emission lines, we are directly viewing the escaping energy from the CGM. As such, the measurement of emission lines originating from the CGM offers the potential opportunity to explore questions regarding the cooling rates. Of course, the problem is more complex than any simple scenario. In the local universe, galactic coronae can be studied via their X-ray emission, but X-ray observations of the Milky Way and nearby galaxies have shown that charge exchange X-ray emission, generated at the interface between the cool and hot gas in the halo, may also make a non-negligible contribution to the cooling (e.g., Henley et al. 2015; Jiang et al. 2019). As such, the reader should bear in mind that our discussion of the mass cooling process in the following is likely to be oversimplified.

Optical emission lines provide an opportunity to explore the cooling rate of the cool CGM, but are extremely challenging to measure. In the local universe, Hα is detected in individual galaxies only when the systems are extreme and difficult to model (such as in the starburst/merger NGC 6240; Yoshida et al. 2016). However, Zhang et al. (2016) presented the first detection of the emission line flux of Hα and [N ii] λ6583, from low-redshift, normal galaxies extending out to ∼100 kpc projected radius by stacking a sample of over 7 million lines of sight from the Sloan Digital Sky Survey (SDSS DR12; Alam et al. 2015). Building on that first result, we have applied the technique to characterize the nature of the emission of the CGM within 50 kpc in low-redshift galaxies (Zhang et al. 2016, 2018a, 2018b, 2019, 2020b, 2020a, hereafter Papers I, II, III, IV, V, VI, respectively). We restrict ourselves to the inner 50 kpc because the measured emission beyond that is strongly affected by emission from nearby galaxies (Paper II). In Papers IV and V, we studied and modeled the environmental dependence and the stellar mass dependence, respectively.

Theoretical studies in the literature (e.g., Mathews & Bregman 1978; Schure et al. 2009; Lykins et al. 2013; Vasiliev 2013) have examined the cooling of the plasma within the temperature range 10⁴–10⁸ K, which can be assumed to be in collisional ionization equilibrium, using the hydrogen and metal line emission fluxes of the elements H, He, C, N, O, Ne, Mg, Si, and Fe for temperatures of 10⁴–10⁸ K, among which the hydrogen emission flux can be dominant for temperatures of 10⁴–10⁵ K. These existing, extensive studies of cooling rates and the computational infrastructure developed to study the ISM and CGM, such as Cloudy (Ferland et al. 1998, 2013, 2017), are key to our study.

We aim to convert the observed Hα emission flux radiating from the CGM of galaxies into an estimate of the mass cooling rate in units of M_⊙ yr⁻¹. We will then use this information to infer the required rate at which infalling gas must be reheated or ejected back into the CGM, under the hypothesis of steady-state equilibrium, by comparing the mass cooling rate to the star formation rate. This will be compared to theoretical estimates of these quantities derived from requiring a match between other observable properties of the galaxies (stellar mass, star formation rates) in their models. Multiple approaches in estimating these quantities are critical given the likely presence of systematic errors in any individual approach. Before doing all of this, we revisit previous empirical results and models of the CGM dependence on stellar mass and environment using both an expanded data set and a slightly more complicated model to ensure that we have the latest estimates of CGM properties with which to estimate the mass cooling rate and feedback mass loading factor.

The paper is organized as follows. In Section 2 we present the data analysis and results, with a focus on how the emission flux depends on galactic stellar mass and local environment. In Section 3.1 we discuss our joint modeling of the CGM cool gas mass fraction on both galactic stellar mass and local environment. We describe our calculation of the cooling rate in Section 3.2, including the application of Cloudy in Section 3.2.1. We present our estimates of the mass cooling rates and mass loading factors in Section 3.3. Lastly, we summarize in Section 4. Throughout this paper, we adopt a ΛCDM cosmology with parameters Ω_m = 0.3, Ω_Λ = 0.7, Ω_k = 0, and the dimensionless Hubble constant h = 0.7 (see Riess et al. 2018; Planck Collaboration et al. 2020).

2. Data Analysis

We follow the approach developed in Papers I through VI (Zhang et al. 2016, 2018a, 2018b, 2019, 2020b, 2020a, respectively) by selecting galaxies that meet our criteria in redshift (0.02 < z < 0.2), half-light radius (1.5 < R₅₀/kpc < 10), and luminosity (10^9.5 < L/L_⊙ < 10¹¹) as primary galaxies. For these galaxies, we extract measurements of the Sérsic index (n) and r-band absolute magnitude (M_r) from Simard et al. (2011), stellar mass (M_*) from Kauffmann et al. (2003a, 2003b) and Gallazzi et al. (2005), and star formation rates (SFRs) from the MPA-JHU catalog (Brinchmann et al. 2004). The SFR estimates are aperture-corrected to account for the light outside the SDSS fiber aperture, which only collects ∼1/3 of the total light for a typical galaxy at the median redshift of the survey (for details see Brinchmann et al. 2004). Our requirements for these ancillary measurements limit the primary sample to galaxies from the seventh major data release of SDSS (DR7), but we utilize spectra from DR16 to probe the CGM of these galaxies.

We collect the spectra from the SDSS DR16 (Ahumada et al. 2020) for all lines of sight projected within a specified range of scaled projected radii from any of our primary galaxies. As described in Paper V, we use scale-projected radii to account for the range in primary galaxy sizes. We define the scaled projected radius, r_s, as the ratio between the physical projected separation and the virial radius of the primary galaxy (r_s ≡ r_p/r_vir). To estimate the virial radius of the primary galaxy, we use the scaling relation between the luminosity and the virial radius obtained by fitting a high-order polynomial to the results drawn from the UniverseMachine (Behroozi et al. 2019). As discussed in more detail in Paper III, we set a physical lower limit on the projected radius (10 kpc) we explore to mitigate possible contamination of the line-of-sight spectra by the central galaxy. We confirm, but discuss no further below, that increasing this lower limit on the projected separation to 20 kpc produces statistically consistent results except for those involving the innermost radial bin of the smallest primary galaxies. Those measurements are affected because this change removes almost all lines of sight from this bin.

Our procedure for the processing of the line-of-sight probes follows from our previous papers. For each such spectrum, we fit and subtract a tenth-order polynomial to a 300 Å wide section surrounding the wavelength of Hα at the related primary galaxy redshift to remove the continuum. We only present and discuss Hα fluxes in this work because the theoretical modeling of [N ii] is more complicated and beyond the question we focus on here. We then measure the residual Hα flux within a velocity window corresponding to the recessional velocity of the associated primary galaxy to capture the majority of the emission flux from the gas surrounding that galaxy. We vary the velocity window as a function of stellar mass as described in Paper V.

We apply this procedure to all of the line-of-sight spectra. Specifically, we have a set of spectral selection criteria. First, the continuum level must be <3 × 10⁻¹⁷ erg cm⁻² s⁻¹ Å⁻¹ to limit the noise introduced by the actual SDSS spectral target. Second, we require that the measured emission line flux be within 3σ of the mean of the whole sample so as to remove spectra of interloping strong emitters such as satellite galaxies or outer disk H i i regions. We have confirmed that using the mean or median of the resulting set of flux measurements in our subsequent analysis produces consistent results. We present the results of mean flux in the following study.

We estimate the uncertainties in the mean flux values by randomly selecting half of the individual spectra in the relevant subsample, calculating the mean emission line flux, and repeating the process 1000 times to establish the distribution of measurements from which we quote the values corresponding to the 16.5, 50.0 and 83.5 percentiles. We compensate for using only half the sample in each measurement by dividing the 1σ estimated uncertainties by a factor of $\sqrt{2}$ .

2.1. Dependence of Hα Emission on the Galaxy's Stellar Mass

We begin our analysis by re-evaluating the dependence of the Hα flux on the primary galaxy stellar mass. We follow the procedure outlined in Paper V and measure the Hα flux for galaxies in different stellar mass bins. The set of primary galaxies used here, taken from SDSS DR7, is somewhat different from that used in Paper V, which was taken from the group catalog from Yang et al. (2007). As in the previous study, we select lines of sight with 0.05 < r_s < 0.25 around each primary galaxy and split the data into two bins with equal ${\rm{\Delta }}\mathrm{log}{r}_{s}$ . We present the Hα flux results, as well as the average SFR, which will be used in the following discussion, in Table 1, and in Figure 1. N in Tables 1 and 2 refers to the number of galaxies in the corresponding subsample. We use the average value of logM_* and r_s to reference each bin.

**Figure 1.** Measured average Hα flux, f, in two r_s bins (0.05 < r_s ≤ 0.11 and 0.11 < r_s ≤ 0.25) as a function of the stellar mass of the primary galaxy. For easier visualization, we apply slight horizontal offsets for the data points of the inner and outer r_s bins. The blue dashed horizontal line represents no net emission or absorption at Hα. The solid black lines represent the modeled flux, and the associated shaded regions correspond to 1σ and 2σ uncertainties of the model described in Section 3.1.
Download figure:
Standard image High-resolution image

Table 1. Stacked Hα Emission Flux vs. Stellar Mass and Radius

${\rm{log}}({M}_{\ast }/{M}_{\odot })$	r_s ^a	N	f	SFR
			[10⁻¹⁷ erg cm⁻² s⁻¹ Å⁻¹]	[M_⊙ yr⁻¹]
9.44	0.09	53	0.038 ± 0.015	0.70 ± 0.06
	0.19	611	0.020 ± 0.004
10.05	0.09	290	0.033 ± 0.006	1.34 ± 0.14
	0.19	2059	0.006 ± 0.002
10.64	0.09	1265	0.013 ± 0.002	1.86 ± 0.04
	0.19	6486	−0.001 ± 0.001
11.07	0.09	966	0.000 ± 0.002	1.59 ± 0.06
	0.19	4236	0.004 ± 0.001

Note.

^a r_s is the ratio between the physical projected separation and the virial radius of the primary galaxy, r_s ≡ r_p/r_vir.

Download table as: ASCII Typeset image

Table 2. Hα Emission Flux vs. Scaled Distance from Nearest Massive Neighbor^a

R_p/R_vir	r_s	N	f	SFR
			[10⁻¹⁷ erg cm⁻² s⁻¹ Å⁻¹]	[M_⊙ yr⁻¹]
0.15	0.09	8	0.044 ± 0.065	0.54 ± 0.07
	0.19	86	−0.015 ± 0.009
0.33	0.09	39	0.070 ± 0.014	0.84 ± 0.06
	0.19	219	−0.012 ± 0.005
0.68	0.09	46	0.017 ± 0.011	0.93 ± 0.05
	0.19	350	0.005 ± 0.046
1.41	0.09	67	0.031 ± 0.010	1.14 ± 0.06
	0.19	625	0.007 ± 0.003
2.92	0.09	131	0.027 ± 0.008	1.56 ± 0.20
	0.19	1022	0.010 ± 0.003
≳10.0	0.09	1137	0.028 ± 0.003	1.17 ± 0.02
	0.19	6572	0.011 ± 0.001

Note.

^aThe scaled distance refers to the projected separation between the primary and the massive neighbor divided by the estimated virial radius of the massive neighbor. A massive neighbor is defined to be at least five times as massive as the primary galaxy and to have M_* ≥ 10¹⁰ M_⊙.

Download table as: ASCII Typeset image

In broad terms, we find that the measured Hα flux decreases as the galaxy stellar mass increases in both r_s bins, and that the flux is consistently larger in the inner r_s bin. These fluxes are quantitatively consistent with those presented in Paper V except for the smallest stellar mass bin. There are two substantive differences between this study and the previous one. We have already mentioned the first, which is that the primary galaxy samples come from slightly different sources. The second is that the line-of-sight spectra came from DR12 for the previous study and come from DR16 for the current study. We confirm that using DR12 with the current primary galaxy sample leads to insignificant differences, so we attribute the discrepancy to differences in the primary galaxy samples. This difference may be a statistical fluctuation or related to the nature of the group catalog. For example, low-mass galaxies in the group catalog are exclusively satellites, while those in the more general SDSS population are likely to be a more heterogeneous population.

2.2. Dependence of Hα Emission on the Galaxy's Local Environment

Quantifying environment is always a challenge, and part of that challenge is whether to focus on a local or global measurement of environment. In Paper IV, one of the environmental measurements we adopted was the projected proximity to a massive (M_* > 10¹¹ M_⊙) galaxy. We justified this choice in the context of the scenario in which the CGM of a primary galaxy is removed once it enters the sphere of influence of a massive neighbor. In that work, we also explored other scenarios, such as one where the CGM of a primary galaxy is removed once it enters a group or cluster environment. Although we also found some evidence in favor of this latter scenario, the number of galaxies in rich environments is small, and exploring this scenario further with the current sample proved difficult. Therefore, we proceed by further exploring only the scenario of a single massive neighbor and caution that environmental effects beyond that explored here may be at play.

We adjust our massive neighbor scenario slightly by lowering the required minimum massive neighbor mass to 10^10.5 M_⊙ and 10¹⁰ M_⊙, in two separate trials, both to increase the sample size and to probe whether the physical effect continues for lower-mass "massive" neighbors. We counter the absolute decrease in neighbor mass by requiring that every affected primary galaxy be substantially less massive (a factor of 5) than its massive neighbor. This additional criterion is set to ensure that we are still considering systems where the primary galaxy is dominated by its neighbor. We chose a factor of 5 because in our previous study (Paper VI) we found that galaxy pairs with a mass ratio less than 5 and low average stellar mass (<10^10.4 M_⊙) exhibit enhanced rather than diminished Hα emission, perhaps due to contributions from tidal material.

As before, we quantify the effect using the projected separation, R_p, between the affected primary and its associated massive neighbor, for potential massive neighbors that satisfy the mass criteria described above and are within recessional velocities offset from the primary by 1000 km s⁻¹. As we did before for the projected separation of lines of sight to the primary, we now scale the projected radius to the massive neighbor using the virial radius of the massive neighbor. Also as before, we estimate the virial radius of the massive galaxy using our derived scaling relation between the virial radius and galactic stellar mass. We confirm that lowering the maximum velocity difference to 500 km s⁻¹ leads to consistent results, although it decreases the sample size by approximately 30%.

A series of steps is required in order to assemble our environmentally affected sample and control samples. When a number, i, of massive galaxies lie within R_p,i/R_vir,i < 4 of the primary galaxy, we consider only the one with the smallest R_p,i/R_vir,i as the massive neighbor. To construct one control sample, we identify a sample of "isolated" primary galaxies that have no neighbors that satisfy our massive neighbor criteria projected within 5 Mpc and lie within 1000 km s⁻¹ in velocity offset. For reference, the average calculated virial radius of the massive neighbors in our sample is ∼500 kpc, which implies that the typical galaxy in the "isolated" galaxy sample is free of massive neighbors to R_p/R_vir ∼ 10. Next, we require R_p/R_vir > 0.05 to remove pairs in which emission from the massive neighbor is likely to contaminate our measurement of the primary galaxy. We find that adopting twice as large a limit makes an insignificant difference in the results. Finally, because the primary galaxy mass distributions are different between the data sample and the control sample, we force the mass distribution of the control sample to be the same as that of the data sample at R_p/R_vir > 1. We do this by randomly drawing from the larger control sample until we match the mass distribution of the data. The Hα fluxes of the mass-matched control sample are (0.028 ± 0.003) × 10⁻¹⁷ erg cm⁻² s⁻¹ Å⁻¹ and (0.011 ± 0.001) × 10⁻¹⁷ erg cm⁻² s⁻¹ Å⁻¹ for the radial bins 0.05 < r_s ≤ 0.11 and 0.11 < r_s ≤ 0.25, respectively.

We find that the results are insensitive to our choice of the lower mass limit for the neighbors, within the range we explore, so we adopt our lowest mass limit, 10¹⁰ M_⊙, to maximize the sample size. In Table 2 and Figure 2 we present the stacked Hα fluxes as a function of R_p/R_vir for both the inner and outer r_s bins around our primary galaxies and isolated galaxy control sample. We refer to this control sample as having R_p/R_vir ≳ 10. We also present the mean SFRs for the primary galaxies in each subsample. Finally, we construct a second control sample in which we insert two artificial primaries at the same R_p/R_vir as each primary, but at different angles around the massive neighbor. We refer to this second control sample as the "azimuthal" control sample with ∼6000 lines of sight, in contrast to the "isolated" control sample with ∼8000 lines of sight.

We draw several conclusions from these results (Table 2 and Figure 2). First, and most significantly, for the outer r_s bin around the primary galaxies, the fluxes for primaries at the two smallest values of R_p/R_vir are significantly lower than for either those at larger values of R_p/R_vir or those that are isolated. To confirm this impression, we calculate the nonparametric rank-order Spearman correlation and find that the null hypothesis of no correlation between flux and R_p/R_vir can be rejected at >99% confidence in the outer r_s bin. Second, less significantly, there is the appearance that the effect exists even beyond R_p/R_vir > 1. Third, the fluxes in the inner r_s bin appear to be unaffected by the presence of the massive neighbor. If anything, there may be a slight flux increase rather than decrease at small R_p/R_vir, but the results for the inner r_s bin have larger uncertainties due to the smaller number of lines of sight. Fourth, for both r_s bins, the flux values for the outermost R_p/R_vir bin are consistent with those of the isolated control sample, indicating that we have measured out to the asymptotic value. Lastly, as expected if there is no contamination from the massive neighbor, we find no mean Hα flux and no systemic dependence of the flux on R_p/R_vir for the azimuthal control sample for either radius bin (Figure 2).

The one remaining puzzle is the slight systematic tendency to produce negative fluxes at small R_p/R_vir in the larger r_s bin. Because we do not find the same trend in the azimuthal control sample, we exclude the possibilities that the effect is due to a systematic overestimate of the continuum in these environments or to the presence of a contribution to the absorption line spectrum from the massive neighbor. At this point, given its somewhat marginal statistical significance, we attribute the result to noise but note that it could also reflect stellar absorption from tidal material around the primary.

3. Models and Interpretation

We now present two different modeling efforts. First, we expand our previous work of defining broad characteristics of the CGM by setting basic rules within the construct of the UniverseMachine models and optimizing the parameter values for those rules. Our effort here is distinguished from our previous work in that we are jointly fitting the mass and environment dependences to determine the interdependence of those, and fitting separately to the measurements in both r_s bins. The latter is important given that it is only the flux in the outer r_s bin that appears affected by the nearby massive neighbor. Second, we present a new analysis where we use the observed Hα fluxes and our derived CGM characterization to constrain the gas mass cooling rate in the halos. In turn, when combined with the independently measured star formation rates, we will use these results to determine how much gas mass needs to be reheated or ejected in order that the central galaxy not systematically gain or lose cool gas. We then compare our measurement to independent theoretical predictions as a demonstration of where future work might lead.

3.1. UniverseMachine and Global Trends

The UniverseMachine is a dark matter-only simulation with comoving box length of 250 Mpc h⁻¹, mass resolution of 1.8 × 10⁸ M_⊙ (2048³ particles), and force resolution of 1 kpc h⁻¹. Halos were found with the Rockstar phase-space halo finder (Behroozi et al. 2013a), and merger trees were generated with the Consistent Trees code (Behroozi et al. 2013b). The code empirically determines the dependence of galaxy star formation rates on the host dark matter halo mass, the halo accretion rate, and redshift, and produced observables are matched to SDSS galaxies in stellar mass and star formation rate. In particular, galaxy colors are matched based on both stellar mass and specific SFR.

In previous studies, Papers IV and V, we studied the dependence of the CGM on the environment and stellar mass separately. However, due to the common multiparameter correlations present among galaxies, it is possible that effects attributed to either environment or mass in those studies should rightly be ascribed to the other. Here we perform a joint modeling exercise to avoid "double-counting" any effect.

The central tenet of this model is the following description of the cool gas fraction that includes a simple dependence on the stellar mass of the primary galaxy, M_*, and environment taken in large part from our previous work:

$\begin{eqnarray}\begin{array}{c}{{ \mathcal F }}_{{\rm{c}}{\rm{o}}{\rm{o}}{\rm{l}}}(R,{M}_{\ast })=\left\{\begin{array}{ll}a{\left({M}_{\ast ,10}\right)}^{b}\left({\textstyle \tfrac{R}{x{R}_{{\rm{v}}{\rm{i}}{\rm{r}}}}}\right) & \,{\rm{i}}{\rm{f}}\,R\lt x{R}_{{\rm{v}}{\rm{i}}{\rm{r}}}\\ a{\left({M}_{\ast ,10}\right)}^{b} & \,{\rm{i}}{\rm{f}}\,R\geqslant x{R}_{{\rm{v}}{\rm{i}}{\rm{r}}}\end{array}\right.\end{array}\end{eqnarray} \tag{ 1 }$

where (x, a, b) are free parameters, R is the three-dimensional distance between the primary galaxy and the massive neighboring galaxy, R_vir is the virial radius of the massive galaxy, which is defined to be the radius of a spherical volume within which the mean density is Δ_c times the critical density at that redshift (where Δ_c = 18π² + 82x − 39x² and x = Ω_m(z) − 1), ${{ \mathcal F }}_{\mathrm{cool}}$ describes the fraction of the cosmological baryon fraction in the halo that is in the cool (T ∼ 10⁴⁻⁵ K) gas phase, and M_*,10 is the stellar mass in units of 10¹⁰ M_⊙. Our assumption that the baryon fraction of galaxy halos is the universal one has not yet been empirically verified. While it is observed to hold roughly for higher-mass halos (Gonzalez et al. 2013; Lim et al. 2020), it does not appear to hold universally for the lower-mass halos of galaxy groups (e.g., Vikhlinin et al. 2009; Gonzalez et al. 2013). On galaxy scales we have few constraints, but the rough concurrence in the Milky Way's dynamical mass and that inferred assuming the cosmological baryon fraction indicates that the baryon fraction for large galaxies, at least for our own galaxy, does not deviate grossly from the universal value (Zaritsky & Courtois 2017). To the degree that the baryon fraction is lower than the universal value in our sample of galaxies, our results underestimate ${{ \mathcal F }}_{\mathrm{cool}}$ . The value of the parameter x affects both the extent of the environmental effect and its amplitude of the environmental effect at a given physical radius. The environmental effect is only applied to the gas in the outer r_s bin, because there is no evidence for stripping of the CGM in the inner r_s bin. For gas in the inner r_s bin we adopt ${{ \mathcal F }}_{\mathrm{cool}}$ as given by the condition R ≥ x R_vir in Equation (1). Additionally, only ∼30% of galaxies have neighbors that are at least five times as massive as the primary galaxy within 4R_vir of the massive neighbor. For those galaxies without such a massive neighbor, ${{ \mathcal F }}_{\mathrm{cool}}$ in both r_s bins is defined using the condition R ≥ x R_vir in Equation (1). Our functional form for the cool gas fraction thus accounts for both stellar mass and environmental dependences, which is notionally consistent with what numerical simulations find when investigating star formation/quenching (e.g., Wang et al. 2018).

We follow our earlier approach in Papers IV and V using the UniverseMachine (Behroozi et al. 2019) catalog, and adopt a gas mass density profile for the primary galaxy consistent with that describing gas in a Navarro–Frenk–White (NFW) potential (Navarro et al. 1996, 1997). As we discussed in Paper II, the Hα emission flux we detect comes from two contributions: the emission from the central galaxy itself and that from the spatially associated halos around the central galaxy. By constraining our measuring and modeling to within projected radii of 50 kpc, or r_s ∼ 0.25, from the primary galaxy, we naturally focus on the first of these components. In Paper IV, we roughly approximated the effect of the second component of Hα emission in those calculations by assigning half of the measured emission to the second component. Now, we account for both components in our models because we include the flux contributions from nearby neighbors. Given the simplicity of this approach, we neglect the higher-order effect represented by the interaction of massive galaxies with the CGM of the other neighbors. Moreover, the agreement in measured fluxes between our outermost R_p/R_vir radial bin and the isolated sample suggests that this effect is not measurable at the precision available in the current data. Furthermore, by modeling the radial profile of the emission, in addition to the total amplitude, we are able to attempt to reproduce the measurements in our two radial bins.

We perform a Bayesian analysis to derive confidence intervals on each of the three parameters in Equation (1). Based on the results from our previous studies we adopt uniform priors over the given ranges as follows:

$\begin{eqnarray}&&0.25\lt x\lt 4.0,\quad 0\lt a\lt 0.5,\quad -1.0\lt b\lt 0,\end{eqnarray} \tag{ 2 }$

and utilize the measured Hα fluxes presented in Tables 1 and 2. There is an implicit physical restriction on the model parameters that the cool gas fraction cannot exceed one. Therefore, we add that constraint by setting any ${{ \mathcal F }}_{\mathrm{cool}}\gt 1$ to 1. As standard, the posterior distribution, p(Θ∣data), is described as follows:

$\begin{eqnarray}&&p({\rm{\Theta }}| {\rm{data}})=\displaystyle \frac{p({\rm{\Theta }})\times {\rm{\Pi }}\,p({\rm{data}}| {\rm{\Theta }})}{p({\rm{data}})}\end{eqnarray} \tag{ 3 }$

where Θ is the parameter space Θ = (x, a, b), p(Θ) is the prior distribution described in Equation (2), p(data∣Θ) is the likelihood of the data, and p(data) is the marginal probability for the data. We define the likelihood of obtaining the data given a specific model using the difference between the actual and model data, $\exp [-0.5\times {({\rm{actual}}-{\rm{model}})}^{2}/{\sigma }^{2}]$ , where σ is the observational uncertainty. We use "emcee" (Foreman-Mackey et al. 2013) to implement a Markov Chain Monte Carlo sampling of the likelihoods across the parameter space to calculate the posterior distribution. In Figure 3 we display the posterior distributions of the three parameters, the median and 1σ confidence levels (based on the 16th and 84th percentiles) for each, and the correlations between the probability distributions.

**Figure 3.** Posterior distribution for the cool gas fraction in the CGM of the primary galaxy for the two radial bins of 0.05 < r_s ≤ 0.11 and 0.11 < r_s ≤ 0.25, as a function of distance, R, from the massive neighbor, and the stellar mass of the primary galaxies. The red vertical line indicates the median value of the model parameters and the shaded region highlights the 1σ uncertainties.
Download figure:
Standard image High-resolution image

Regarding the modeling of the environmental dependence, the preferred value for x is ${1.92}_{-0.88}^{+1.05}$ , demonstrating that the effects of a massive neighbor on the CGM of the primary galaxy are evident even when the primary lies beyond the neighbor's virial radius. Unfortunately, the precision on x is low, presumably because the sample is still relatively small (we only have 1447 data points for R_p /R_vir < 2) and because we only have projected distances between primaries and neighbors. Even so, the finding that x > 1 suggests either that the physical effect begins at these large radii or that our sample includes galaxies that have 'splashed' back, having already passed within the virial radius, and currently lie outside. Various modeling and sample differences complicate the comparison to our previous result in Paper IV, but there we estimated the effect to extend to a projected separation of ∼724 kpc for a sample of neighbors with a mean virial radius of ∼550 kpc, or alternatively an x value of ∼1.32, which is consistent with our current estimate.

Regarding the modeling of the mass dependence, we find a = 0.23 ± 0.01, which is only slightly smaller than, but consistent with, our previous estimate of ${0.28}_{-0.06}^{+0.07}$ in Paper V, and b = −0.55 ± 0.03, which is significantly lower than our previous estimate of −0.33 ± 0.06 in Paper V. To explore the cause of this marked change, we model the stellar mass dependence and the environmental effect separately, as done in our previous studies. Our new best fit for b is −0.55 ± 0.03, demonstrating that the difference from our previous result does not lie with the new joint analysis of mass and environment.

We ascribe the significant difference in b to two factors. First, there is a statistically significant difference in the measured flux for the lowest-mass galaxies in the innermost r_s bin between the two studies. The flux was lower in the previous study and this might be related to the selection of group galaxies in that previous work, most of which would be satellites at these low masses and might differ from a more general sample of low-mass galaxies. Second, the basic power-law model we use is manifestly overpredicting the measured flux for these low-mass galaxies at the smallest separations (Figure 1), suggesting that this model may not be an accurate representation and that the dependence of ${{ \mathcal F }}_{\mathrm{cool}}$ on mass may not be quite as steep as indicated by our new estimate of b.

We show in Figures 1 and 2 that the most likely model is also one that reasonably reproduces the observational data. With regard to the mass dependence of ${{ \mathcal F }}_{\mathrm{cool}}$ , we that find a strong decline in the cool gas fraction with increasing stellar, and presumably total, mass is needed.

With regard to the environmental effects, we interpret the model as indicative of the removal of the outer portions of the CGM of the primary galaxy by the massive neighbor (Putman et al. 2021). The effect is significant even when the primary galaxy falls in the potential well within 2 × R_vir of the massive neighbor, and most severe when the primary galaxy lies within the virial radius of the massive neighbor. The model suggests that satellites rapidly lose some of their gaseous halos, or that their hot gaseous halos quickly become unable to cool after reaching the virial radius of a larger halo. The apparent onset of the effect at radii larger than the virial radius could reflect the nature of the splashback radius (for discussion see More et al. 2015), where galaxies that have fallen inside the virial radius return to larger radii as they execute their orbits, or the affected nature of mass accretion onto halos at these radii (Behroozi et al. 2014). The significant removal of the gas but only small decrease in the SFR might indicate the delayed quenching of those satellites in rich environments (Wetzel et al. 2013).

3.1.1. An Aside on Hα Emission and Galaxy's Color

As is always the case for galaxy studies, a variety of properties are related to each other and this complicates the analysis and interpretation. We just measured that Hα flux from the CGM is inversely related to the stellar mass of the galaxy, but stellar mass is related to a number of other galaxy properties, such as color and morphology. Is it possible that a relationship between Hα flux and galaxy color, which is a parameter that is tabulated for our sample, provides an easy-to-access additional constraint on our modeling?

At a given stellar mass, one might expect blue, star-forming galaxies to have a richer gas reservoir. Indeed, we do find that the bluer galaxies have larger Hα emission, as shown in Figure 4, and that color is highly correlated to the galaxy stellar mass (for g − r bins centered on 0.34, 0.52, 0.73, and 0.83, the corresponding mean stellar masses are 10^9.8 M_⊙, 10^10.34 M_⊙, 10^10.76 M_⊙, and 10^10.91 M_⊙, respectively).

**Figure 4.** Measured average Hα flux for 0.05 < r_s < 0.25 as a function of g − r. The solid black line represents the model flux and the corresponding shaded regions represent to 1σ and 2σ uncertainties of the model described in Section 3.1. Again, zero net flux is presented as a blue dashed line to guide the eye.
Download figure:
Standard image High-resolution image

We now use the UniverseMachine models to see if there is potential additional information in the relationship between Hα flux and galaxy color. Galaxy color is included in the output of our previous modeling, so we now compare the model predictions to the data in Figure 4. We find good agreement, which is at a minimum a reassuring consistency check and suggests that galaxy color is highly degenerate with galaxy stellar mass. The result does not inform us as to whether stellar mass or color is the primary driver of the observed trend in flux, only that ascribing it to one or the other is likely to produce similar results. Given the expected rise in the virial temperature of galaxy halos with halo mass, which is presumably tracked by stellar mass, we prefer to base our simple modeling of ${{ \mathcal F }}_{\mathrm{cool}}$ on stellar mass rather than on galaxy color.

3.2. Modeling the Mass Cooling Rate and Inferring the Mass Loading Factor

We do not have a quantitative, complete, empirically validated understanding of how star-forming galaxies sustain their star formation. On one hand, ongoing star formation in a large fraction of galaxies demonstrates that star formation lasts at least ∼10 billion years in those systems, but on the other hand, the depletion time of the molecular gas for low-redshift normal galaxies is in the range 1–2 Gyr (Saintonge et al. 2011; Tacconi et al. 2018). The natural solution is that the galaxy disk is being refueled and that the source of that fuel is, by consensus, the CGM (Borthakur et al. 2015; Tumlinson et al. 2017).

Our measurement of the Hα flux coming from the CGM provides an opportunity to measure the rate at which the CGM is radiatively losing energy, and hence cooling. Because the radiative energy loss must lead to the equivalent loss of thermal energy in the gas, we will estimate the rate at which CGM gas is cooling and, presumably, sinking into the central disk. A complication that we address below is that Hα is only one of the emission lines through which the gas is cooling. Assuming that the disk maintains a certain gas mass on average over some current modest time interval, the excess of the mass inflow rate over the star formation rate is the amount of gas that is reheated or ejected. This is the mass loading factor of whatever feedback mechanism is at work regulating star formation.

Once we have calculated the total radiative losses across the entire wavelength range (Section 3.2.1), we will set that equal to the thermal energy loss in the gas. Assuming an ideal gas, the energy radiated per particle as it changes in temperature by ΔT is

$\begin{eqnarray}&&{\rm{\Delta }}E=\displaystyle \frac{3}{2}k{\rm{\Delta }}T.\end{eqnarray} \tag{ 4 }$

In our calculation, we divide the entire temperature range into a few bins and ΔT is the change in temperature of the gas parcel that is considered to contribute to the measured Hα flux. For a mean particle mass of μ m_p, where m_p is the proton mass, the rate at which mass is cooling (the mass cooling rate), $\dot{m}$ , is

$\begin{eqnarray}&&\dot{m}=(\dot{E}/{\rm{\Delta }}E)\times \mu {m}_{p},\end{eqnarray} \tag{ 5 }$

where $\dot{E}$ is the radiated energy rate for the whole halo and is calculated as follows:

$\begin{eqnarray}&&\dot{E}={R}_{\mathrm{fiber}}\times {f}_{{\rm{H}}\alpha }\times {C}_{\mathrm{Cl}}\end{eqnarray} \tag{ 6 }$

where R_fiber is the ratio of the flux contained within the fiber aperture to that contributed by the full halo, f_Hα is the Hα flux measurements described above in units of erg cm⁻² s⁻¹, and C_Cl is the unitless correction factor derived from Cloudy modeling (Section 3.2.1), which converts the Hα radiative energy loss to that across all emission lines. R_fiber is defined as $\pi {({D}_{{\rm{A}}}{r}_{\mathrm{fiber}})}^{2}\times {C}_{{\rm{G}}}$ , where D_A is the angular diameter distance at the redshift of our primary galaxy sample, which on average is approximately 0.1, r_fiber is the eBOSS angular fiber size, and C_G is a unitless geometric correction factor that we describe in Section 3.2.2.

3.2.1. Determining C_Cl

We model the dependence of the emission spectrum from the CGM on the ionization state (logU), the metallicity ([X/H]), the gas temperature (T), the total hydrogen column density ( $\mathrm{log}{N}_{{\rm{H}}}$ ), and the hydrogen density (ρ_H) using Cloudy (Ferland et al. 1998, 2013, 2017) version 17.02. For the ionization spectrum we adopt either the background radiation field from quasars and galaxies of Haardt & Madau (2001), or the O-star spectrum of Kurucz (1979) corresponding to an effective temperature of the star of 30,000 K. We assume that the gas is in ionization equilibrium.

We find that the results are relatively insensitive to the choice of ionization source spectrum (resulting in <10% changes to the mass cooling rate estimated below) and so we only present results from calculations using the O-star spectrum as the ionization source.

One complication that arises when converting the Hα flux into an energy loss rate is that the halo gas emitting Hα is not at a single temperature, but rather consists of gas at a large range of temperatures, approximately 10^4–5 K (e.g., Schure et al. 2009; Lykins et al. 2013), all of which can contribute to the observed flux. We present results from a calculation to illustrate how the flux and inferred cooling rate depend on the gas temperature T. For this exercise we set the other parameters to those inferred by Werk et al. (2014) and Prochaska et al. (2017): ionization parameter logU = −3.0, metallicity [X/H] = −0.75, total hydrogen column density $\mathrm{log}({N}_{{\rm{H}}}$ /cm⁻²) = 19.5, and hydrogen density ρ_H = 10^−3.0 cm⁻³. We do not include dust in the Cloudy modeling. The interplay of total hydrogen column density and the hydrogen density determines the geometry of the clouds; in other words, it sets the total length of the gas cloud along the line of sight. We vary the gas temperature T from 10^3.5 K to 10⁶ K. The selected temperature range roughly covers all the different phases for the halo gas in the CGM for a Milky Way–like galaxy that will emit Hα.

The fraction of all the cooling radiation that comes from Hα emission (which we refer to as the Hα fraction), the Hα radiation intensity, and the total cooling/heating rate are shown as functions of temperature in Figure 5. At higher ionization ( $\mathrm{log}U\gt -3.5$ ) and temperatures <10⁴ K the cooling rate drops by almost two orders of magnitude and the heating rate exceeds the cooling rate, suggesting that this simple model is missing cooling channels, such as molecular lines. Therefore, we set the lower boundary of the temperature range we will consider at 10⁴ K. At the high-temperature end (T > 10⁵ K), the Hα intensity and the Hα fraction drop rapidly, so we also consider the contribution of this gas to Hα to be negligible. At lower ionization ( $\mathrm{log}U\lt -3.5$ ), the heating rate is orders of magnitude lower than the cooling rate at all temperature. Although we have limited the temperature range of the gas that is relevant to our calculation, the range remains large, 10^4–5 K, and Hα fraction and intensity vary significantly across this range, suggesting that the estimate of the distribution of gas across this range of temperatures is fundamental to understanding the origin of Hα properties.

**Figure 5.** The fraction of the total cooling rate attributed to the Hα emission line (top), the Hα line intensity (middle), and the total cooling/heating rate (bottom) as functions of the gas temperature as calculated by Cloudy for logU = −3.0, $[{\rm{X}}/{\rm{H}}]=-0.75,\mathrm{log}$ (N_H/cm⁻²) = 19.5, and ρ_H = 10^−3.0 cm⁻³.
Download figure:
Standard image High-resolution image

**Figure 5.** The fraction of the total cooling rate attributed to the Hα emission line (top), the Hα line intensity (middle), and the total cooling/heating rate (bottom) as functions of the gas temperature as calculated by Cloudy for logU = −3.0, $[{\rm{X}}/{\rm{H}}]=-0.75,\mathrm{log}$ (N_H/cm⁻²) = 19.5, and ρ_H = 10^−3.0 cm⁻³.
Download figure:
Standard image High-resolution image

Unfortunately, this information is currently unavailable to the required accuracy and, in order to make progress, we need to introduce additional assumptions. Because the cooling rates vary with temperature, there are naturally temperature regimes that the gas transitions through faster and slower as it cools. To avoid gas piling up at certain temperatures, we enforce a continuous cooling mass flow through the temperature regimes by scaling the gas mass at the different temperatures inversely by its corresponding mass cooling rate. As expected, this implies that most of the gas ends up at temperatures near 10⁴ K, where the cooling rate falls steeply. This approach also means that if one calculates the mass cooling rate at any single temperature, that value is the mass cooling rate at all temperatures.

In addition, we also require the integrated Hα flux in the temperature range 10⁴–10⁵ K to match the observational values shown in Tables 1 and 2. After scaling, the integrated Hα intensity over the temperature range 10⁴–10⁵ K is 0.81 × 10⁻⁹ erg cm⁻² s⁻¹. The aperture size of eBOSS is 2 arcsec, which gives a conversion factor from intensity to flux of 9.41 × 10⁻¹¹; then the integrated Hα flux over the SDSS fiber size in the temperature range 10⁴–10⁵ K is 0.06 × 10⁻¹⁷ erg cm⁻² s⁻¹. The mass cooling rate corresponding to this Hα flux is estimated to be 0.15 M_⊙ yr⁻¹ before geometric correction.

This approach allows us to estimate the mass cooling rate from the observed Hα fluxes, building upon the derived scaling relations between the Cloudy Hα flux and the Cloudy mass cooling rate. Nonetheless, this result may be sensitive to our assumptions in the Cloudy modeling, as we discuss below.

3.2.2. Determining C_G

To determine the geometric correction factor, C_G, we follow a two-step process. First, we correct for the fraction of the projected area covered by the eBOSS fiber within the radial range probed by our study (0.05 < r_s ≤ 0.25). Because we average results from thousands of fibers scattered throughout this area with no correlation to the primary galaxy, the mean flux in a fiber aperture accurately represents a sampling of the flux within the full area, modulo the ratio of the fiber aperture to the full aperture. As such, we simply correct by the ratio of the two areas $({A}_{{r}_{s}\lt 0.25}/{A}_{\mathrm{fiber}})$ . The correction factor varies for different galaxies because the physical radii are different even for the same range of scale radius of 0.05 < r_s < 0.25. Then the mean correction factor is 14.3, 36.5, 128.9, and 682.5 for the four stellar mass bins described above. Second, we consider the correction for the ratio of the full halo to that which the 0.05 < r_s ≤ 0.25 annulus represents, which we define as ${{ \mathcal F }}_{\mathrm{corr}}$ . We find in our model that although there is nearly as much baryonic mass in the gaseous halo outside r_s = 0.25 as inside, there is a negligible contribution to the Hα flux (see Figure 7 in Zhang et al. 2018). Therefore, ${{ \mathcal F }}_{\mathrm{corr}}$ is approximately equal to 1. Our calculation for C_G is

$\begin{eqnarray}&&{C}_{G}=\displaystyle \frac{{A}_{{r}_{s}\lt 0.25}}{{A}_{\mathrm{fiber}}}\times {{ \mathcal F }}_{\mathrm{corr}},\end{eqnarray} \tag{ 7 }$

where we set ${{ \mathcal F }}_{\mathrm{corr}}=1$ .

3.2.3. Uncertainties in the Derived Mass Cooling Rate

To develop intuition regarding the sensitivity of our results to various parameter choices, we investigate the dependence of the total mass cooling rate on the adopted ionization parameter (U), the metallicity ([X/H]), and the hydrogen density (ρ_H), while keeping the column density fixed. Prochaska et al. (2017) presented a radial distribution of hydrogen column density measurements in L_* galaxies, roughly ranging from $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})=20$ to $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})=19$ for projected radius between 10 and 50 kpc, respectively, so we adopt $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})=19.5$ . For the other parameters we allow:

1.
logU to vary from −6 to −1 (our preferred parameter range is from −4 to −2),
2.
[X/H] to vary from −2.5 to 0.5 (our preferred parameter range is from −1.5 to 0), and
3.
ρ_H to vary from 10^−2.0 cm⁻³ to 10^−4.0 cm⁻³ (our preferred parameter range is from 10^−2.5 to 10^−3.5 cm⁻³).

Our preferred parameter ranges are based on the constraints from absorption line measurements (Werk et al. 2014; Prochaska et al. 2017) using QSO spectra from the COS-Halos Survey (Werk et al. 2013). Given the fixed column density, the range of hydrogen densities we probe corresponds to a column length of the simulated halo between 10 kpc and 300 kpc.

Dependence on the ionization parameter— The ionization parameter, which measures the number of ionizing photons per atom, affects the calculated heating and cooling rates, and therefore the resulting Hα flux. Here, we fix the other parameters to be [X/H] = −0.75 and ρ_H = 10^−3.0 cm⁻³, and vary only logU and the temperature. As described above, for each ionization parameter, we enforce continuity in the gas mass flow rate at different temperatures (Section 3.2.1) and then estimate the mass cooling rate.

The result is shown in Figure 6. We find that changes of about a factor of three in the estimated mass cooling rate are possible, although the change is rather sudden at about logU ∼ −3.5. At our preferred value of logU = −3 we expect to be calculating the comparably lowest likely values of the mass cooling rate, and plausible changes in logU should only produce higher estimates of the mass cooling rate.

**Figure 6.** Mass cooling rate as a function of the ionization parameter, logU, for fixed $[{\rm{X}}/{\rm{H}}]=-0.75,{\rm{log}}({N}_{{\rm{H}}}$ /cm⁻²) = 19.5, and ρ_H = 10^−3.0 cm⁻³. What we consider to be the plausible range of logU, based on Werk et al. (2014), is identified by vertical dashed lines and the corresponding range of plausible mass cooling rates by the shaded region.
Download figure:
Standard image High-resolution image

Dependence on metallicity— Metals contribute significantly to the total cooling in the temperature range 10^4–5 K (e.g., Schure et al. 2009; Lykins et al. 2013). Here we investigate how the adopted gas metallicity affects the inferred mass cooling rate. We adopt logU = −3.0 and ρ_H = 10^−3.0 cm⁻³, and vary the gas metallicity and the temperature. Again, we enforce continuity of gas cooling across temperature ranges and estimate the mass cooling rate, as described in Section 3.2.1.

The result is shown in Figure 7. As expected, the mass cooling rate increases as the metallicity increases for metallicity below solar metallicity, resulting in the inferred mass cooling rate increasing significantly as we adopt a higher metallicity of the gas. At our preferred value of [X/H] = −0.75, we are near the lower limit of the inferred mass cooling rates. If the CGM has solar metallicity we would be underestimating the mass cooling rate by about a factor of 2 for our standard assumed [X/H], but it is highly unlikely that the CGM reaches solar metallicity (Prochaska et al. 2017 find that the typical CGM metallicity is ∼−0.5). For more plausible errors in the metallicity, the corresponding errors are significantly smaller.

**Figure 7.** Mass cooling rate as a function of the metallicity, [X/H], for fixed logU = − 3.0, ${\rm{log}}({N}_{{\rm{H}}}$ /cm⁻²) = 19.5, and ρ_H = 10^−3.0 cm⁻³. What we consider to be the plausible range of [X/H], based on Werk et al. (2014) and Prochaska et al. (2017), is identified by vertical dashed lines, and the corresponding range of plausible mass cooling rates by the shaded region.
Download figure:
Standard image High-resolution image

Dependence on the density— Both the cooling rate and the Hα flux are proportional to the gas density squared. Because we have fixed the total hydrogen column density, our inferred cooling rates and Hα fluxes are proportional to the gas density, not the density squared. In the Cloudy calculation, we have adopted a uniform gas density as a simple approximation to the halo gas density profile, which roughly follows a power law according to Werk et al. (2014) in other simulations. Adopting logU = −3.0 and [X/H] = −0.75, we obtain the results shown in Figure 8. As expected, the mass cooling rate is proportional to the adopted gas density.

**Figure 8.** Dependence of the mass cooling rate on the gas density (ρ_H) for fixed logU = − 3.0, [X/H] = −0.75, and ${\rm{log}}({N}_{{\rm{H}}}$ /cm⁻²) = 19.5. What we consider to be the plausible range of ρ_H, based on Werk et al. (2014) and Prochaska et al. (2017), is identified by vertical dashed lines and the corresponding range of plausible mass cooling rates by the shaded region.
Download figure:
Standard image High-resolution image

We close this section by noting that the model we present is a highly oversimplified version of the CGM. Highlighting a few specific shortcomings, we note that (1) the density we model with Cloudy is independent of the galactocentric distance, whereas we know that the true density of the galaxy halo decreases significantly with radius, (2) we do not consider density variations within a given column (i.e., clumpiness), and (3) we only consider the variation of ${{ \mathcal F }}_{\mathrm{cool}}$ with galaxy mass and ignore other potential mass dependences, such as one between CGM metallicity and mass. The model we present is best at establishing broad behavior and testing our conceptual understanding rather than at providing a quantitative understanding. Ultimately, the data presented here will provide the strongest constraints when compared with forward-modeling of physically motivated, high-resolution CGM simulations.

3.3. The Mass Cooling Rate

3.3.1. Empirical Results

We integrate the flux within the radial range 0.05 < r_s ≤ 0.25 over the velocity window described in Section 2 and then apply the appropriate geometric correction described in Equation (7) to obtain the CGM Hα emission within the virial radii of galaxies in the different stellar mass bins. We also calculate the CGM Hα flux using our fitted model for ${{ \mathcal F }}_{\mathrm{cool}}$ in Equation (1) and the best-fit parameter values for the assumed T = 12,000 K. We will present both results.

To infer a mass cooling rate (MCR) from those emission fluxes, we adopt our preferred parameter set (logU = −3.0, [X/H] = −0.75, ${\rm{log}}({N}_{{\rm{H}}}$ /cm⁻²) = 19.5, and ρ_H = 10^−3.0 cm⁻³ (from Werk et al. 2014; Prochaska et al. 2017). From Figures 6–8, we know that the mass cooling rate depends significantly on the adopted Cloudy model parameters. To estimate the effect of these uncertainties, we randomly uniformly distribute 10,000 points in the 3D parameter space within the range (−4 < logU < −2, −1.5 < [X/H] < 0, ${10}^{-3.5}\lt {\rho }_{{\rm{H}}}/{{\rm{cm}}}^{-3}\lt {10}^{-2.5})$ . We follow the steps discussed in Section 3.2.1 to calculate the mass cooling rate from the Hα flux derived from our model for ${{ \mathcal F }}_{\mathrm{cool}}$ using Cloudy. We adopt the standard deviation among those 10,000 estimates of the mass cooling rate as the 1σ uncertainty of the modeled mass cooling rate.

We present our results for the mass cooling rates in Figure 9 and compare those rates to the current SFR as a function of galaxy stellar mass. As expected, the mass cooling rates are much higher across the board than the SFR of the central galaxies. We draw several conclusions from this comparison: (1) for our preferred parameters the mass cooling rate is well above what is needed to maintain the current rate of star formation, thereby addressing the gas consumption problem; (2) this estimated mass cooling rate significantly exceeds what is needed and therefore supports the theoretical understanding that the conversion of cooling gas into stars is inefficient and that a large fraction of the gas must be reheated or ejected; and (3) the strength of this "feedback" is the largest for the galaxies with the largest stellar masses, but otherwise appears roughly constant as a function of stellar mass across this limited mass range. The approximate linearity in log–log space indicates that the relation between logMCR and $\mathrm{log}{M}_{* ,10}$ , where M_*,10 is the stellar mass in units of 10¹⁰ M_⊙, can be approximated by a second-order polynomial function as follows:

$\begin{eqnarray}&&\mathrm{logMCR}=0.81+0.48\mathrm{log}{M}_{* ,10}+0.35{\mathrm{log}}^{2}{M}_{* ,10}.\end{eqnarray} \tag{ 8 }$

**Figure 9.** Total mass cooling rates (MCRs) and star formation rates (SFRs) of the primary galaxies as a function of stellar mass. The error bars on the blue circles (the cooling rates calculated from our Hα flux measurements) reflect only the flux measurement uncertainties and do not reflect any uncertainty in the adopted Cloudy parameters. In contrast, the mass cooling rates calculated from our fitted model for ${{ \mathcal F }}_{\mathrm{cool}}$ , the solid green line, include the systematic modeling uncertainties (shaded region). We also plot the mass cooling rates predicted by the GAEA simulation (Xie et al. 2020) as triangles. The SFRs, both empirical (red stars) and as predicted in the GAEA realizations (yellow stars), provide benchmark comparisons. Finally, we also plot our polynomial fit to our model as the red dashed line (Equation (8)).
Download figure:
Standard image High-resolution image

**Figure 9.** Total mass cooling rates (MCRs) and star formation rates (SFRs) of the primary galaxies as a function of stellar mass. The error bars on the blue circles (the cooling rates calculated from our Hα flux measurements) reflect only the flux measurement uncertainties and do not reflect any uncertainty in the adopted Cloudy parameters. In contrast, the mass cooling rates calculated from our fitted model for ${{ \mathcal F }}_{\mathrm{cool}}$ , the solid green line, include the systematic modeling uncertainties (shaded region). We also plot the mass cooling rates predicted by the GAEA simulation (Xie et al. 2020) as triangles. The SFRs, both empirical (red stars) and as predicted in the GAEA realizations (yellow stars), provide benchmark comparisons. Finally, we also plot our polynomial fit to our model as the red dashed line (Equation (8)).
Download figure:
Standard image High-resolution image

When drawing quantitative conclusions from our results, such as inferring the fraction of the gas that needs to be reheated or ejected to maintain the SFR at a steady state, one must bear in mind the large systematic uncertainties that exist. The error bars as plotted on the individual data points reflect only the uncertainties due to the measurement errors on the Hα fluxes. The uncertainties plotted on the model predictions reflect systematic uncertainties in the adopted parameters. In either case, our previous qualitative conclusions drawn from Figure 9 hold, but these uncertainties, particularly the systematic ones, obviously significantly affect a quantitative estimate of the ratio of the cooling mass that must be reheated or ejected to that resulting in new stars.

3.3.2. Comparison to Previous Theoretical Results

Given the importance of gas cooling to galaxy evolution, there is a large set of existing literature examining this topic. Of particular relevance to our results, there is theoretical work estimating the gas cooling flow. For example, Stern et al. (2019) derived the gas cooling flow rate for a Milky Way–like galaxy. Concisely put, their simulations are controlled numerical experiments in which gas that is initially in hydrostatic equilibrium within an external gravitational potential is allowed to cool radiatively. They estimate the ratio between the cooling flow rate and the SFR ( $\dot{M}$ /SFR) to be ∼10 within 100 kpc for a Milky Way–like galaxy of halo mass ∼10¹² M_⊙. This is broadly consistent with what we have shown are our results in Figure 9. Again, broadly, this estimate is consistent with results obtained utilizing the radiatively cooling gas traced by O vi (Mathews & Prochaska 2017; McQuinn & Werk 2018; Stern et al. 2018). In particular, McQuinn & Werk (2018) obtained an estimation of the mass flow of cooling gas of ∼10–30 M_⊙ yr⁻¹ from the O vi absorption for low-redshift L_* galaxies.

To pursue a more quantitative comparison, and in particular to examine the trend in cooling rate with galaxy mass, we explore a more direct comparison to simulations. In particular, we consider the predictions of the Galaxy Evolution and Assembly semianalytic model that has been developed over a series of studies (GAEA; De Lucia & Blaizot 2007; De Lucia et al. 2014; Hirschmann et al. 2016; Xie et al. 2017, 2020; Zoldan et al. 2018; Fontanot et al. 2018, 2020). It is characterized by an accurate treatment of chemical enrichment and energy recycling, a realistic stellar feedback model, a self-consistent partition of multiphase cold gas, and an explicit treatment of the environmental effects on satellite galaxies. The model produces good agreement with the observed stellar mass function up to z ∼ 7, and with the cosmic star formation history up to z ∼ 10 (Fontanot et al. 2017). It is also able to reproduce the H i and H₂ fractions, as well as the quenched fractions in the local universe for central and satellite galaxies respectively (Xie et al. 2020). In particular, the theoretical predictions in the following have been taken from the model realization applied to the Millennium Simulation (Springel et al. 2005) and discussed in Xie et al. (2020). In this semianalytic model, baryons in galaxies are distributed in discrete components: hot gas, cold gas, stellar disk, stellar bulge, and central massive black hole. In detail, Xie et al. (2020) assume that cold gas is partitioned into H i and H₂, and that stars form from the H₂ component following the empirical law proposed by Blitz & Rosolowsky (2006). This approach implies that the SFR is therefore a function of midplane pressure of the gas disk. The gas disk consists only of cold gas that either cools from a hydrostatic hot atmosphere when the cooling time is longer than dynamical time or comes directly from outside the dark matter halo when the cooling time is shorter. The former cooling regime dominates at low redshift. When cooling from a hydrostatic hot gas halo, the cooling rate is calculated using the cooling function (Sutherland & Dopita 1993; De Lucia et al. 2010). In the latter case, all of the matter accreted into the dark matter halo cools onto the gas disk. Cooling is also allowed for satellite galaxies in GAEA. The cold gas amassing in galactic disks is then supposed to fuel star formation. The resulting feedback reheats the cold gas at a rate that depends on the SFR, on the circular velocity, and on redshift (Hirschmann et al. 2016), using functional forms that have been inspired from the results of high-resolution N-body simulations (Muratov et al. 2015). Reheated gas is moved from the cold gas component to the hot gas component, leading to a complex interplay between cooling and heating. Moreover, radio-mode feedback from active galactic nuclei (AGNs) (due to inefficient gas accretion onto the central massive black holes) is assumed to provide an additional heating channel (Croton et al. 2006; Xie et al. 2017), which further prevents the hot gas from cooling in the more massive systems.

The resulting mass cooling rates plotted in Figure 9 are the combination of the simulated mass of the cooling gas retained by the disk and the mass of the reheated gas. In addition, the models present calculated SFRs. The calculated SFRs are naturally a good fit to the measured ones because the models are in part tuned to the reproduce related constraints such as galaxy stellar masses. More impressive, however, is the independent agreement between the model predictions for the mass cooling rates and our measurements. The only significant disagreement appears at the lowest galaxy masses, where we have already noted that our model may have a problem (Section 3.1) and the GAEA model may also have a problem, partially due to limited resolution (Xie et al. 2017, 2020).

3.4. The Feedback Mass Loading Factor

3.4.1. Empirical Results

As discussed above, the mass cooling rate is significantly higher than the SFR, demonstrating that the conversion of cooling gas into stars is inefficient and that a large fraction of the gas must be reheated or ejected. We define the mass loading factor, η, to be the ratio between the net mass cooling rate (the difference between the mass cooling rate and the SFR) and the SFR. Note that this definition may differ from that in other studies that focus on physical mass outflow. Our estimate includes the effects both of outflow and of any gas reheating that prevents that gas from participating in star formation.

The mass loading factors as calculated from our model of CGM Hα fluxes show a rise of about a factor of 10 across the galaxy mass range we explore (Figure 10) and in all cases have values that are significantly larger than 1. The rise is steeper at the larger M_* in our sample, in the log–log plot, indicating that the dependence of η on M_* rises faster than a simple power law. The relation in the limited mass range between $\mathrm{log}\eta$ and $\mathrm{log}{M}_{* ,10}$ , where M_*,10 is again the stellar mass in units of 10¹⁰ M_⊙, can be approximated by a second-order polynomial function as follows:

$\begin{eqnarray}&&\mathrm{log}\eta =0.66+0.27\mathrm{log}{M}_{* ,10}+0.46{\mathrm{log}}^{2}{M}_{* ,10}.\end{eqnarray} \tag{ 9 }$

**Figure 10.** The feedback mass loading factor, η, calculated using our model for ${{ \mathcal F }}_{\mathrm{cool}}$ , with resulting systematic uncertainties plotted as the shaded region, and the results from two different published theoretical studies are presented vs. primary galaxy stellar mass. The dashed red line presents the second-order polynomial fit to our results (Equation (9)). For comparison, we plot the mass loading factors from the GAEA (Xie et al. 2020), the EAGLE (Mitchell et al. 2020), and FIRE simulations (Muratov et al. 2015; Hayward & Hopkins 2017) as triangles, diamonds, and squares, respectively. The apparent discrepancy with the FIRE results is resolved in the text.
Download figure:
Standard image High-resolution image

3.4.2. Comparison to Previous Theoretical Results

In Figure 10 we compare our estimates of the mass loading factors to three theoretical predictions. First, we compare to the same GAEA models we discussed previously. The agreement is excellent, which is as expected given the agreement in both the predicted versus observed SFRs and mass cooling rates (Figure 9).

Next we compare to the results of Mitchell et al. (2020), who present the mass loading factor predicted using the Virgo Consortium's Evolution and Assembly of Galaxies and their Environments (EAGLE; Crain et al. 2015; Schaye et al. 2015) model suite. EAGLE is a suite of hydrodynamical simulations that simulate the formation and evolution of galaxies within the context of the ΛCDM cosmological model, integrated down to z = 0. The models account for important physical processes that are not resolved by the simulation (radiative cooling, star formation, stellar mass loss and metal enrichment, growth of a supermassive black hole, and energy injection from stellar and AGN feedback). They calibrate the feedback to the present-day galaxy stellar mass function and the amplitude of the galaxy–central black hole mass relation. When doing that, they find good agreement with the observed galaxy specific star formation rates, passive fractions, Tully–Fisher relation, total stellar luminosity of galaxy clusters, and column density distributions of intergalactic C iv and O vi. Their predictions are in excellent agreement with our results across the mass range we cover (Figure 10).

Finally, we compare to results presented from the FIRE (Feedback in Realistic Environments) collaboration (Hopkins et al. 2014, 2018). They executed a series of high-resolution cosmological simulations of galaxy formation to z = 0, spanning halo masses ∼10⁸–10¹³ M_⊙, and stellar masses ∼10⁴–10¹¹ M_⊙. These simulations explicitly treat the multiphase ISM with heating and cooling physics from gas at a range of temperatures T ∼ 10–10¹⁰ K, star formation restricted only to self-gravitating, self-shielding, molecular, high-density (n_H ≳ 5–50 cm⁻³) gas, and (most importantly) explicit treatment of stellar feedback including the energy, momentum, mass, and metal fluxes from supernovae of Types Ia and II, stellar mass loss (O star and asymptotic giant branch), radiation pressure (UV and IR), and photoionization and photoelectric heating. These sources of feedback reproduce the observed relation between stellar and halo mass up to M_halo ∼ 10¹² M_⊙, roughly corresponding to M_* ∼ 10^10.5 M_⊙. Star formation rates agree well with the observed Kennicutt relation at all redshifts.

In Muratov et al. (2015) and Hayward & Hopkins (2017), the investigators identified whether a particle is considered as outflowing or inflowing according to its radial velocity. Once outflowing particles are identified, the authors calculate outflow rates in M_⊙ yr⁻¹ by computing the instantaneous mass flux through a thin spherical shell. Combining with the SFR from the simulation, they calculate the mass loading factors for low-mass galaxies with M_* < 10¹⁰ M_⊙ and an upper limit for galaxies with M_* > 10¹⁰ M_⊙ at redshift ∼0.25, which are shown as squares in Figure 10. A more recent analysis with FIRE-2 simulations (Pandya et al. 2021) is able to convert these limits to measurements, and those measurements lie slightly below the previous limits. Note that these studies focus on outflowing gas, and therefore that the estimate of the mass loading factor, as defined here, using their results is a lower bound on what we should be measuring.

At the low-mass end with M_* < 10¹⁰ M_⊙, the FIRE simulation is statistically consistent with our estimation. At M_* > 10¹⁰ M_⊙, however, the results of the FIRE simulation significantly diverge from our estimates. However, this disagreement does not necessarily reflect a fault in either our measurements or the FIRE results. Again recall that their results only represent true outflow and do not include reheating. Furthermore, their simulations are intended to reproduce only galaxies with M_halo < 10¹² M_⊙, which for our stellar mass–halo mass relation corresponds to a stellar mass of 10^10.5 M_⊙. Lastly, the divergence in Figure 10 is likely also reflecting the lack of feedback other than from stars in these models, e.g., AGN feedback, which is included in the GAEA and EAGLE modeling. The agreement among the mass loading factors from our empirical results, the semianalytic results from GAEA, the results from hydrodynamical simulations of EAGLE, and the more detailed, physically motivated simulations from FIRE at low masses shows an encouraging convergence on our understanding of mass loading factors at those masses.

The agreement between simulations with sub-grid physics prescriptions and our own highly simplified analysis of our Hα observations is highly encouraging. The concurrence of the GAEA and EAGLE models, which include AGN feedback, with the empirical results, and the discrepancy of the FIRE results—at the high-mass end—with the same empirical results, is a clear indication of the importance of AGN feedback in understanding the character of the CGM for massive galaxies. We had speculated the same on the basis of the change in emission line ratios measured for the CGM for low- and high-mass galaxies (Paper III), but the comparison here is an even stronger indication of such.

4. Summary

We present measurements of the mass cooling rates for CGM gas and infer feedback mass loading factors from the comparison of those rates to the measured star formation rates as a function of galactic stellar mass.

To calculate these estimates, we first revisited our measurements of the Hα emission line flux as a function of stellar mass, M_*, and as a function of the local environment using both newer data and an improved treatment. We track the behavior of Hα flux within two scaled annular bins around our primary galaxies as a function of both stellar mass (10^9.5 < M_*/M_⊙ < 10¹¹) and scaled projected distance from a massive neighbor (0.05 < R_p/R_vir < 4). The flux in both annular bins drops as M_* increases, indicating a strong inverse relationship between the cool gas fraction of the CGM, ${{ \mathcal F }}_{\mathrm{cool}}$ , and M_*. Here ${{ \mathcal F }}_{\mathrm{cool}}$ describes the fraction of the cosmological baryon fraction in the halo that is in the cool (T ∼ 10⁴⁻⁵ K) gas phase. We also find a strong environmental effect for the outer annular bin at the inner two smallest values of R_p/R_vir. We fit a model that includes both mass and environment dependences within the framework of the UniverseMachine modeling (Behroozi et al. 2019) and find

$\begin{eqnarray}&&{{ \mathcal F }}_{\mathrm{cool},\mathrm{isolated}}({M}_{* })=(0.23\pm 0.01){\left({M}_{* ,10}\right)}^{(-0.55\pm 0.03)}\end{eqnarray} \tag{ 10 }$

for galaxies beyond the influence of a nearby, massive neighbor (as defined in Section 3.1), and

$\begin{eqnarray}&&{{ \mathcal F }}_{\mathrm{cool}}(R,{M}_{* })=\displaystyle \frac{R}{{1.92}_{-0.88}^{+1.05}{R}_{\mathrm{vir}}}\cdot {{ \mathcal F }}_{\mathrm{cool},\mathrm{isolated}}\end{eqnarray} \tag{ 11 }$

otherwise, where R is the physical distance between the primary and the massive neighbor and R_vir is the virial radius of the massive neighbor. Our prediction for the cool gas fraction is as high as 90% for the smallest halos we probe (although we may be overestimating the fraction at the low-mass end, see Section 3.1), 23 ± 1% for a Milky Way–like galaxy, and only a few percent for the most massive galaxies in our sample. These results, except for the lowest-mass galaxies, are consistent with our previous results in Paper IV and V, despite changes to the data and treatment. The only category where there is some tension with the previous work is for the inner annular bin for the lowest-mass galaxies.

Our measurement of the Hα flux originating from the CGM, and our models, provide us the opportunity to measure the rate at which the CGM is radiatively losing energy, and hence cooling. We use the Cloudy modeling suite (Ferland et al. 1998, 2013, 2017) to understand what fraction of the cooling is contributed by the Hα flux we measure. In doing this, we need to account for gas at different temperatures and present arguments for why we consider only gas in the range 10⁴–10⁵ K. By enforcing a continuous cooling mass flow through the different temperatures we are able to calculate the gas mass at the different temperatures within this range and complete the calculation.

We then estimate the rate at which the CGM gas is cooling and presumably sinking into the central disk. We explore the dependence of the CGM emission spectra on the ionization state, the metallicity, and the hydrogen density, and estimate the uncertainties in our recovered mass cooling rates. We present both the gas mass cooling rate and effective mass loading factor required to avoid growing the mass of disk gas. The mass cooling rate is a factor of several higher than the star formation rate of the central galaxies. The simple fitted relations between the mass cooling rate, MCR, and the mass loading factor, η, with galactic stellar mass are as follows:

$\begin{eqnarray}&&\mathrm{logMCR}=0.81+0.48\mathrm{log}{M}_{* ,10}+0.35{\mathrm{log}}^{2}{M}_{* ,10}\end{eqnarray} \tag{ 12 }$

$\begin{eqnarray}&&\mathrm{log}\eta =0.66+0.27\mathrm{log}{M}_{* ,10}+0.46{\mathrm{log}}^{2}{M}_{* ,10}\end{eqnarray} \tag{ 13 }$

where M_*,10 is the galactic stellar mass in units of 10¹⁰ M_⊙.

Although we have used highly idealized models to transform our measured Hα flux values into estimates of the mass cooling rates and feedback mass loading factors, it is heartening to find excellent agreement between our estimates of those quantities and those derived from three theoretical treatments that are themselves quite different in their approach (Xie et al. 2018, 2020; Hopkins et al. 2014, 2018; Mitchell et al. 2020). The general agreement we find suggests that we have a relatively robust quantitative understanding of the cooling and feedback rates for galaxies in this mass range. The empirical confirmation of these rates provided here supports further modeling of the role of the CGM in galaxy evolution. It is perhaps somewhat surprising that our highly idealized model produces results in such good agreement given the likely complexity of the CGM (Tumlinson et al. 2017), but perhaps our measurement, which involves massively averaging over thousands and thousands of galaxy lines of sight, provides a bulwark against the complexity of any individual galaxy's CGM. The field critically needs maps of the CGM in individual galaxies to further refine our understanding of this key galactic component.

D.Z. and H.Z. acknowledge financial support from NSF grant AST-1311326. K.P.O. is funded by NASA under award No 80NSSC19K1651. X.Y. is supported by the national science foundation of China (grant Nos. 11833005, 11890692, 11621303) and 111 project No. B20019. T.F. is supported by the National Key R&D Program of China No. 2017YFA0402600. The authors gratefully acknowledge the SDSS-III team for providing a valuable resource to the community. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III website is http://www.sdss3.org/.

SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.

An Empirical Determination of the Dependence of the Circumgalactic Mass Cooling Rate and Feedback Mass Loading Factor on Galactic Stellar Mass

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Data Analysis

2.1. Dependence of Hα Emission on the Galaxy's Stellar Mass

2.2. Dependence of Hα Emission on the Galaxy's Local Environment