Quenching low-mass satellite galaxies: evidence for a threshold ICM density

TL;DR

This study investigates how the dense intra-cluster medium influences the quenching of low-mass satellite galaxies, revealing a threshold ICM density where rapid quenching occurs due to ram pressure stripping.

Contribution

It provides the first direct observational evidence of a threshold ICM density triggering rapid quenching in low-mass satellite galaxies, supported by a simple analytic model.

Findings

Quenched fraction increases sharply beyond a certain ICM density threshold.

A broken power-law describes the relation between quenched fraction and ICM density.

Ram pressure stripping explains the rapid quenching at high ICM densities.

Abstract

We compile a sample of SDSS galaxy clusters with high-quality Chandra X-ray data to directly study the influence of the dense intra-cluster medium (ICM) on the quenching of satellite galaxies. We study the quenched fractions of satellite galaxies as a function of ICM density for low- ($10^9 \lesssim M_\star \lesssim 10^{10}\,\mathrm{M_\odot}$), intermediate- ($10^{10} \lesssim M_\star \lesssim 10^{10.5}\,\mathrm{M_\odot}$), and high-mass ($M_\star \gtrsim 10^{10.5}\,\mathrm{M_\odot}$) satellite galaxies with $>\!3000$ satellite galaxies across 24 low-redshift ($z < 0.1$) clusters. For low-mass galaxies we find evidence for a broken powerlaw trend between satellite quenched fraction and local ICM density. The quenched fraction increases modestly at ICM densities below a threshold before increasing sharply beyond this threshold toward the cluster center. We show that this increase in quenched fraction at high ICM density is well matched by a simple, analytic model of ram pressure stripping. These results are consistent with a picture where low-mass cluster galaxies experience an initial, slow-quenching mode driven by steady gas depletion, followed by rapid quenching associated with ram pressure of cold-gas stripping near (one quarter of the virial radius, on average) the cluster center.

Tables

Table 1. Table 1: Galaxy cluster sample

Name	Yang IDa	$z_{\mathrm{cluster}}$b	$M_{\mathrm{halo}}$c	$R_{500}$	$L_{X,\mathrm{bol}}$d	$N_{\mathrm{gal}}/N_{\mathrm{Yang}}$e	Exp. time	Chandra ObsID
			$\left({10}^{14}{M}_{\odot }\right)$	(kpc)	$\left({10}^{44}\mathrm{erg}{s}^{-1}\right)$		(ks)
Coma	1	0.024	10.8	1150	7.1∗	517/652	479	13993,13994,13995
								13996,14406,14410
								14411,14415
Abell 2147	2	0.036	9.6	1078	0.7	329/383	18	3211
Abell 1367	3	0.022	6.2	955	1.2∗	223/352	402	514,17199,17200
								17201,17589,17590
								17591,17592
								18704,18705,18755
Abell 2199	5	0.030	5.8	920	0.8	221/302	158	497,498,10748
								10803,10804,10805
Abell 85	11	0.056	8.4	1030	5.3	152/178	198	904,15173,15174
								16263,16264
Abell 2063	18	0.035	3.8	799	4.3	132/156	50	4187,5795,6262
								6263
Abell 2670	19	0.076	12.6	1115	2.3∗	107/154	40	4959
Abell 2029	20	0.077	9.4	1041	15.6	114/154	128	891,4977,6101
Abell 2065	21	0.072	8.6	1013	3.0	104/154	55	3182,7689
Abell 2142	22	0.090	17.6	1259	75.9	149/153	205	5005,15186,16564
								16565
MKW 8	23	0.027	1.2	697	0.9	132/150	104	4942,18266,18850
Abell 2107	24	0.041	3.8	787	3.0	119/149	36	4960
Abell 2052	25	0.035	4.0	815	5.1	103/141	654	890,5807,10477
								10478,10479,10480
								10879,10914,10915
								10916,10917
Abell 2255	27	0.082	13.6	1167	3.7	107/130	44	894,7690
Abell 2061	29	0.078	5.1	1072	2.5∗	115/129	55	4965
Abell 1795	32	0.063	5.8	895	7.8	103/125	105	493,494,10432,17228
ZwCl 1215	45	0.077	7.3	954	3.5	93/103	12	4184
Abell 1991	52	0.058	5.2	862	0.5	96/99	38	3193
MKW 3S	57	0.045	3.1	743	1.1	77/95	57	900
MKW 4	62	0.021	1.7	625	0.03	58/88	30	3234
Abell 1775	71	0.075	5.7	880	2.8∗	70/82	99	12891,13510
Abell 1650	86	0.084	7.3	947	4.4	56/72	251	4178,5822,5823
								6356,6357,6358
								7242,7691
AWM 4	145	0.032	0.7	564	0.4∗	44/55	74	9423
Abell 2244	191	0.098	4.7	805	5.5	44/47	65	4179,7693
								13192,13193
NGC 4325	611	0.026	0.3	336	0.2∗	18/24	30	3232

Equations

$$R_{500}=R_{200m}/2.7$$

$$R_{200m}=1.61\mathrm{Mpc}{\left(\frac{M_{\mathrm{halo}}}{10^{14}{\mathrm{M}}_{\odot }}\right)}^{1/3}(1+z_{\mathrm{cluster}}{)}^{-1}$$

$$n_H=\sqrt{\frac{4\pi D_A^2(1+z{)}^2\eta \cdot 10^{14}}{1.2V}}$$

$$n_H=n_{H,0}{\left[1+{\left(\frac{R}{R_c}\right)}^2\right]}^{-\frac{3}{2}\beta }$$

$$\begin{array}{rl}{n}_H& =n_{H,01}{\left[1+{\left(\frac{R}{R_{c1}}\right)}^2\right]}^{-\frac{3}{2}{\beta }_1}\\ & +n_{H,02}{\left[1+{\left(\frac{R}{R_{c2}}\right)}^2\right]}^{-\frac{3}{2}{\beta }_2}\end{array}$$

$$\begin{array}{rl}{\rho }_{\mathrm{ICM}}(R)v^2>[g_{\mathrm{DM}}(r)+& g_{d,\star }(r)+g_b(r)+\\ g_{\mathrm{HI}}(r)+g_{{\mathrm{H}}_2}(r)]{\Sigma }_{\mathrm{gas}}(r)& \end{array}$$

$$v=\sqrt{3}\times \frac{|z-z_{\mathrm{cluster}}|}{1+z_{\mathrm{cluster}}}\times c$$

$$g_{d,\star }=2G{\Sigma }_{d,\star }(r)=2G{\Sigma }_0{e}^{-r/R_d}$$

$$M_d=2\pi {\int }_0^{\infty }{\Sigma }_{d,\star }(r)r\mathrm{d}r=2\pi {\Sigma }_0{R}_d^2.$$

$${\varphi }_b(r,z)=-\frac{G{M}_b}{(r^2+{\mathcal{Z}}^2{)}^{1/2}+a}$$

$$\begin{array}{rl}{g}_b(r)& =\underset{\mathcal{Z}}{\max }\frac{\partial {\varphi }_b(r,\mathcal{Z})}{\partial \mathcal{Z}}\\ & =\underset{\mathcal{Z}}{\max }\frac{G{M}_b}{[(r^2+{\mathcal{Z}}^2{)}^{1/2}+a{]}^2}\frac{\mathcal{Z}}{(r^2+{\mathcal{Z}}^2{)}^{1/2}}.\end{array}$$

$$M_{\mathrm{strip}}=M_{\mathrm{gas}}{e}^{-r_{\mathrm{strip}}/R_{\mathrm{gas}}}\left(\frac{r_{\mathrm{strip}}}{R_{\mathrm{gas}}}+1\right)$$

$$f_{\mathrm{strip}}=\frac{M_{\mathrm{strip},\mathrm{HI}}+M_{\mathrm{strip},{\mathrm{H}}_2}}{M_{\mathrm{HI}}+M_{{\mathrm{H}}_2}}.$$

$$t_{\mathrm{slow}}=\left(\frac{R_{200m}-R_{\mathrm{thresh}}}{v_{\mathrm{slow}}}\right)$$

$$f_Q=1-(1-f_{Q,0})e^{-t/{\tau }_Q}.$$

$${\tau }_{Q,\mathrm{rapid}}=t_{\mathrm{rapid}}\times {\left[\ln \left(\frac{1}{1-f_{Q,\mathrm{rescaled}}({\rho }^{\prime })}\right)\right]}^{-1}$$

$$t_{\mathrm{rapid}}=\left(\frac{R_{\mathrm{thresh}}-R({\rho }^{\prime })}{{\bar{v}}_{\mathrm{rapid}}}\right).$$

$$\begin{array}{r}{R}_{\mathrm{proj},\mathrm{sim}}=[0.0,0.1,0.2,0.3,0.4,0.5,0.6,\\ 0.7,0.8,0.9,1.0]R_{200m},\end{array}$$

$$\begin{array}{r}\frac{R_{3D,\mathrm{sim}}}{R_{\mathrm{proj},\mathrm{sim}}}=[2.49,1.55,1.32,1.21,1.16,1.11,\\ 1.08,1.07,1.04,1.01].\end{array}$$

$$R_{\mathrm{deproj},\mathrm{obs}}=S(R_{\mathrm{proj},\mathrm{obs}}/R_{200m})\times R_{\mathrm{proj},\mathrm{obs}}$$

Full text

Galaxies: clusters: general, Galaxies: clusters: intracluster medium, Galaxies: evolution, Galaxies: star formation, X-rays: galaxies: clusters

††software: This work was made possible due to a large number of open-source software packages, including: AstroPy (Astropy Collaboration et al., 2013), Colossus (Diemer, 2017), Matplotlib (Hunter, 2007), NumPy (van der Walt et al., 2011), Pandas (McKinney, 2010), Photutils (Bradley et al., 2016), SciPy (Jones et al., 2001), and Topcat (Taylor, 2005)

1 Introduction

It is now firmly established that local environment plays a pivotal role in dictating the properties of galaxy populations. Galaxies located in the underdense field tend to be blue in colour, with young stellar populations, disc (late-type) morphologies, and high star formation rates (SFRs). In contrast, dense environments such as galaxy clusters, embedded within massive ($\gtrsim 10^{14}\,\mathrm{M_{\odot}}$) dark matter halos, host galaxy populations which are on average red, with bulge-dominated (early-type) morphologies, old stellar populations, and little ongoing star formation. The first clear evidence for such a paradigm was presented in early seminal works (e.g. Oemler, 1974; Butcher & Oemler, 1978; Davis & Geller, 1976; Dressler, 1980; Postman & Geller, 1984), and has since been cemented by more recent studies using large, detailed spectroscopic surveys of galaxies across a variety of environments (e.g. Blanton & Moustakas, 2009; Kimm et al., 2009; Peng et al., 2010; Wetzel et al., 2012; Wilman & Erwin, 2012). Even within individual clusters, galaxy properties are a strong function of environment. Galaxies that inhabit the dense cluster interior are preferentially red, early-type, and quiescent relative to galaxies at large cluster-centric radius (e.g. Postman et al., 2005; Blanton & Roweis, 2007; Prescott et al., 2011; Rasmussen et al., 2012; Fasano et al., 2015; Haines et al., 2015). Ultimately, if we are to understand galaxy evolution we must understand which physical mechanisms drive these observed trends by quenching star formation and transforming morphology as a function of galaxy environment.

In addition to environmental quenching occuring in a galaxy’s present day cluster, there is increasing evidence that a large fraction of galaxies may have their star formation quenched in smaller groups prior to cluster infall. It has been estimated that nearly half of present day cluster galaxies may have infallen as a part of a smaller groups (e.g. McGee et al., 2009), and quenched fractions at the cluster virial radius tend to be enhanced relative to the field – indicating that star formation is influenced environmentally prior to infall (von der Linden et al., 2010; Haines et al., 2015; Roberts & Parker, 2017). This “pre-processing” of galaxy properties is a natural consequence of a hierachical growth of dense structures in the Universe.

Many physical mechanisms have been proposed which, in principle, are capable of quenching star formation in dense environments. These mechanisms can be broadly divided into two classes: hydrodynamic interactions between galaxies and the intracluster medium (ICM); and dynamical interactions between member galaxies, or between galaxies and the cluster halo potential. Examples of hydrodynamic mechanisms include: ‘starvation’ (e.g. Larson et al., 1980; Balogh et al., 2000; Peng et al., 2015), where the high virial temperature of the cluster ($\gtrsim 10^{7}\,\mathrm{K}$) prevents cold-flow accretion of gas onto the disc of satellite galaxies; and ‘ram pressure stripping’ (e.g. Gunn & Gott, 1972; Quilis et al., 2000), where a galaxy passing through the dense ICM will feel a ram pressure ‘wind’ which is strong enough to directly strip cold-gas from the galactic disc. Dynamical interactions thought to be relevant include: galaxy mergers (e.g. Mihos & Hernquist, 1994a, b), where the final end products in dense environments tend to be quiescent galaxies with early-type morphologies; ‘harassment’ (e.g. Moore et al., 1996), where repeated impulsive interactions between galaxies can induce strong starbursts, thereby quickly exhausting cold-gas reserves; and gravitational tidal forces (e.g. Mayer et al., 2006; Chung et al., 2007), which can directly strip gas from a galaxy, or in the less extreme case, transport gas to less bound orbits where it will be more susceptible to hydrodynamic effects such as ram pressure. While it is understood that all of the aforementioned processes should be affecting galaxies in dense environments, it is the balance between these mechanisms, and the dependence of this balance on environment that fuels substantial debate.

Recently, starvation or ram pressure stripping (or a combination of the two) have been favoured as the primary quenching mechanism in galaxy groups and clusters (Muzzin et al., 2014; Peng et al., 2015; Fillingham et al., 2015; Wetzel et al., 2015; Brown et al., 2017; Foltz et al., 2018). Substantial effort has been devoted to determine how to distinguish between these two quenching pathways observationally. A common technique is to constrain the timescale over which quenching occurs, which is expected to be relatively long ($\gtrsim 3-4\,\mathrm{Gyr}$) for starvation but short ($\lesssim 1\,\mathrm{Gyr}$) for ram pressure stripping. It is important to note, however, that ram pressure stripping of cold-gas will not be immediately efficient upon cluster infall and a delay time is likely necessary for the galaxy to reach the dense interior ICM before quenching begins. Therefore the total quenching timescale for ram pressure (delay + quenching) should be on the order of the cluster dynamical time. Starvation should begin to act immediately after infall as the galaxy encounters the hot, virialized halo, however quenching by starvation will still produce an excess of quenched galaxies in the cluster interior (relative to the outskirts) as the time-since-infall for these central galaxies will be relatively long.

The effects of ram pressure stripping can, in some cases, be studied directly by observing cluster galaxies with extended Hi distributions (Kenney et al., 2004; Chung et al., 2007, 2009; Kenney et al., 2015), Hi deficient discs (post-stripping, Kenney & Young 1989; Boselli & Gavazzi 2006; Jaffé et al. 2016), and by observing “jellyfish galaxies” with extended, stripped “tentacles” of gas and stars (Poggianti et al., 2017; Jaffé et al., 2018). Ram pressure also lends itself well to analytic modelling through the simple balance between the restoring potential of a galactic disc and the strength of ram pressure given by $\rho_{\mathrm{ICM}}\,v_{\mathrm{galaxy}}^{2}$ (Gunn & Gott, 1972). Such an approach has been used to constrain the regions in cluster phase space where stripping should be efficient, finding that the “stripping” regions tend to be populated by galaxies which are Hi-deficient and show morphological signs of ongoing stripping (Jaffé et al., 2015, 2016, 2018).

In this paper we compile a statistical sample ($>3000$) of satellite galaxies in 24 low redshift ($z<0.1$) clusters observed with the Chandra X-ray observatory to directly study the connection between galaxy quenching and ICM density. We can then estimate the physical ICM density (based on density profiles for each cluster) around each satellite galaxy and constrain the impact of ram pressure stripping using a simple analytic model. This represents the first systematic study of environmental quenching directly as a function ICM density for a large, statistical sample of galaxies across many clusters.

The paper is organized as follows: in section 2 we describe the sample of cluster satellite galaxies, as well as the optical and X-ray data for the host clusters; in section 3 we present the dependence of satellite quenched fraction on ICM density; in section 4 we describe an analytic ram pressure stripping model we construct and make comparisons to the observed trends; in section 5 we constrain quenching timescales assuming a “slow-then-rapid” framework for satellite quenching; and finally, in sections 6 and 7 we discuss and summarize our results.

This paper assumes a flat $\mathrm{\Lambda}$ cold dark matter cosmology with $\Omega_{\mathrm{M}}=0.3$, $\Omega_{\mathrm{\Lambda}}=0.7$, and $H_{0}=70\,\mathrm{km}\,\mathrm{s^{-1}}\,\mathrm{Mpc^{-1}}$. Throughout this paper we will use lowercase $r$ to represent galactocentric radii (in cylindrical coordinates), and uppercase $R$ to represent cluster-centric radii.

2 DATA

2.1 Cluster sample

We construct a sample of low redshift galaxy clusters with high quality, archival X-ray observations, starting with all clusters in the Yang et al. (2005, 2007) SDSS DR7 catalogue at $z<0.1$ with ten or more member galaxies. We then query the Chandra data archive at the positions of the luminosity-weighted centres of the Yang clusters in this initial sample. We require that any matches have at least 25000 X-ray counts above the background (after combining all Chandra archival observations for a given cluster), which results in 24 clusters that both are in the Yang et al. catalogue and are observed with Chandra to sufficient depth. Table 1 lists the clusters in this sample, Fig 1a shows the $L_{X}-M_{\mathrm{halo}}$ relation for these clusters, and Fig. 1b shows the projected offset between the position of the X-ray peak and the luminosity-weighted cluster centre, showing that typical offsets are only a small fraction of the virial radius (see equation 2). X-ray centres are calculated as the position of the brightest pixel in the X-ray image after smoothing using a Gaussian kernel with a bandwidth of $40\,\mathrm{kpc}$ (Nurgaliev et al., 2013; Roberts et al., 2018) and luminosity-weighted centres are computed using the positions of member galaxies from the Yang et al. catalogue. Bolometric X-ray luminosities are archival and are taken from the ACCEPT cluster sample (Cavagnolo et al., 2009) when available, and otherwise from the Rosat All Sky Survey (Wang et al., 2014b). Cluster halo masses are taken from the Yang et al. catalogue which are computed using abundance matching based on the ranking of the characteristic group stellar mass, given by equation 13 in Yang et al. (2007). We normalize all cluster-centric radii by $R_{500}$ (the radius at which the interior density is 500 times the critical density of the Universe) or by $R_{200m}$, which are computed as

[TABLE]

where,

[TABLE]

is the radius enclosing an average density equal to 200 time the critical mass density of the Universe (Yang et al., 2007; Tinker et al., 2008; Wang et al., 2014b).

2.1.1 ICM density profiles

The Chandra observations were downloaded, reprocessed, cleaned, and calibrated using ciao (ciao version 4.9, caldb version 4.7.7). We apply charge transfer inefficiency and time-dependent gain corrections and filter for background flares using the ciao script lc_clean with a $3\sigma$ threshold. Point sources are identified and masked using the wavdetect script. Backgrounds were estimated using the blanksky event file output from the ciao script blanksky which is normalized to the ratio of observed-to-blanksky background counts in the $9-12\,\mathrm{keV}$ band. To determine the ICM density as a function of radius we extract X-ray spectra in radial annuli from the source and the background data sets using the ciao script specextract in the $0.5-7\,\mathrm{keV}$ energy band. Annuli are centered on the X-ray peak and extend out to the edge of the chip coverage. Weighted response files and redistribution matrices were generated with a count-weighted map across the extent of the extraction regions. We then fit the X-ray spectra in radial annuli for each cluster in the sample. We set a minimum of 5000 counts (after background subtraction) per annulus, which is based on the merged dataset for clusters with multiple observations (see table 1), and then group spectra to have 25 counts per energy channel (ie. $\mathrm{S/N=5}$). We only consider clusters with at least 5 radial annuli (ie. 25000 counts), and set an upper limit of 20 annuli per cluster – this means that clusters with deep Chandra observations will have far more than 5000 counts per annulus. To gain insight into the physical densities in each annulus, we deproject the spectra using the code dsdeproj (Sanders & Fabian, 2007; Russell et al., 2008), which is a model-independent deprojection method assuming only spherical symmetry. We then fit the deprojected spectrum in each annulus with an absorbed single-temperature apec model111http://www.atomdb.org/, with temperature, abundance, and normalization as free parameters. For clusters with multiple Chandra observations we extract a spectrum for each dataset and then simultaneously fit all spectra in sherpa222http://cxc.harvard.edu/sherpa/. The redshift is fixed at the cluster redshift from the Yang catalogue. The Galactic hydrogen column density is estimated from the spectrum extracted over the entire cluster region, which gives sufficient counts to obtain a good constraint. We then assume that the Galactic hydrogen column density ($N_{H}$) is constant across the cluster and fix $N_{H}$ at this fitted value for the spectral fits within each radial annulus. The fitted hydrogen column densities agree well with the observed values from Kalberla et al. (2005). The spectral normalization is then converted into a hydrogen number density using the following relation,

[TABLE]

where $D_{A}$ is the angular diameter distance to the cluster at redshift, $z$, $\eta$ is the spectral normalization, $V$ is the volume corresponding to a given annulus, and we have assumed $n_{e}=1.2\,n_{H}$. The deprojected density profile for each cluster is then fit with both a single (equation 4) and a double (equation 5) beta model,

[TABLE]

[TABLE]

with the central density, $n_{H,0}$, core radii, $R_{c}$, and beta indices, $\beta$, as free parameters. Given the single and double beta fit for each cluster, we use the fit which gives the lowest Akaike information criterion333The Akaike information criterion is a model selection tool used to quantify the information lost when describing a set of data with a particular model. When comparing two models, the model with the lowest Akaike information criterion value is preferred. The Akaike information criterion also includes a penalty for increasing the number of fits parameters, thereby accounting for obtaining “better” fits by increasing model complexity. (AIC; Akaike 1974) value. Allowing clusters to be parameterized by a double-beta model when necessary allows for the more accurate modelling of systems with strong cool-cores, however we note that using a single-beta model for all clusters in the sample does not change the conclusions of this paper. The data prefer a double-beta model for 8 of the clusters in the sample, with the remaining 16 preferring a single-beta fit. The ICM density profile fits are done using the MCMC code emcee444http://dfm.io/emcee/current/ (Foreman-Mackey et al., 2013) assuming a flat prior; best-fit parameters are taken to be the median values from each chain after burn-in. Fig. 2 shows the best-fit profiles for each cluster as well as 1$\sigma$ error bars corresponding to the measured densities. These fits allow us to determine the local, azimuthally averaged, ICM density for each satellite galaxy at a given cluster-centric radius (extrapolating the fits to larger cluster-centric radii when necessary).

2.2 Galaxy sample

We compile a sample of cluster satellite galaxies beginning with the member galaxies for the clusters in Table 1 from the Yang et al. DR7 catalogue (we remove the most-massive galaxy from each cluster). Galaxies are then matched to star formation rates (SFRs) from the Max Planck Institut für Astrophysik and Johns Hopkins University (MPA-JHU) collaboration555https://www.mpa-garching.mpg.de/SDSS/DR7 (Brinchmann et al., 2004), with the updated prescriptions from Salim et al. (2007). To determine specific star formation rates ($\mathrm{sSFR}=\mathrm{SFR}/M_{\star}$) we use stellar masses from Mendel et al. (2014) derived via fits to galaxy broadband spectral energy distributions (SEDs). Mendel et al. also derive bulge and disc stellar masses using the bulge+disc (B+D) decompositions from Simard et al. (2011) (assuming an exponential disc and a De Vaucouleurs bulge), which we use in our ram pressure stripping model (see section 4). After matching the Yang et al. member galaxies to the SFR, stellar mass, and structural catalogues we are left with a total of 3250 galaxies in 24 clusters.

We also make use of the isolated field sample from Roberts & Parker (2017) for comparison to the cluster sample. The field sample is composed of all galaxies in single-member groups from the Yang et al. catalogue with a minimum separation of $1\,\mathrm{Mpc}$ and $1000\,\mathrm{km}\,\mathrm{s^{-1}}$ from their nearest ‘bright’ neighbour. Bright neighbours correspond to galaxies which are brighter than the $r$-band limiting absolute magnitude of the survey at $z=0.1$, which ensures that the strictness of the isolation criteria is redshift independent. Finally, any galaxies within $1\,\mathrm{Mpc}$ of the survey edge, or $1000\,\mathrm{km}\,\mathrm{s^{-1}}$ of the maximum redshift are removed. Stellar masses and SFRs are obtained for the field sample from the same sources discussed above, resulting in an isolated, field sample of 164193 galaxies.

3 SATELLITE QUENCHING VS. ICM DENSITY

It is well established by previous works that the fraction of quenched galaxies is a strong function of cluster-centric radius, from the local Universe out to at least $z\sim 1$ (Bamford et al., 2009; Muzzin et al., 2012; Presotto et al., 2012; Wetzel et al., 2012; Haines et al., 2015). Given the fact that there is significant scatter in the cluster-centric radius vs. ICM density relation across different clusters (see Fig. 2), it is interesting to explore quenched fraction trends as a function of ICM density directly. These trends can be used to explore quenching mechanisms such as ram pressure stripping (see section 4.1), the strength of which depends explicitly on ICM density. Fig. 3 shows quenched fraction ($\log\,\mathrm{sSFR}<-11\,\mathrm{yr^{-1}}$, Wetzel et al. 2013) versus ICM density for low- (purple), intermediate- (gray), and high-mass (red) galaxies, computed in equally spaced bins of ICM density. Galaxy stellar mass subsamples are defined such that there are an equal number of galaxies in each mass bin. We also mark the quenched fraction for the isolated, field sample (described in section 2.1) in each mass bin as the dashed horizontal lines in Fig. 3.

We find that the quenched fraction at the lowest ICM densities, which corresponds to the vicinity of the virial radius ($R_{200m}$; equation 2), is significantly enhanced relative to the value for the field, at the same stellar mass. This is particularly true for the two lower mass bins, whereas the trend for the highest mass galaxies approaches the field value. Numerous previous studies (e.g. Lu et al., 2012; Bahé et al., 2013; Haines et al., 2015; Roberts & Parker, 2017) have identified that the quenched fraction in the cluster outskirts can be significantly enhanced relative to the field, a fact which is often attributed to the pre-processing of star formation in less dense environments (ie. small groups) prior to infall onto galaxy clusters. For the low-mass galaxies specifically, the data interpretted this way would require a pre-processed fraction of $\sim 30$ per cent.

For intermediate- and high-mass galaxies the quenched fraction increases smoothly toward high ICM density. The environmental effect on the higher-mass galaxies is weaker than for the lowest-mass galaxies, consistent with previous studies that argue environment most strongly influences low-mass galaxies ($M_{\star}\lesssim\mathrm{few}\times 10^{10}\,\mathrm{M_{\odot}}$, Haines et al. 2006; Bamford et al. 2009; Peng et al. 2010). To ensure that the observed trends are not being driven by our particular binning scheme, we rebin the data from 10 bins up to 20 bins and fit the resulting trend for intermediate- and high-mass galaxies with a single powerlaw. The fits to the 10 rebinnings are shown as the faded lines in Fig. 3. While increasing the number of bins adds noise and increases the statistical uncertainties, the fits remain nearly identical and the underlying trend is robust. Given that the trend with ICM density is seemingly different for low-mass galaxies, these low-mass objects will be the focus of the remainder of the paper; however we will present a discussion of differences between high- and low-mass galaxies in Section 6.

The quenched fraction trend for low-mass galaxies shows signs of a broken powerlaw, with a moderate increase at the lowest ICM densities and a steepening in the densest cluster regions. To quantify this trend, we fit the quenched fraction trend for low-mass galaxies with both a single (SPL) and a broken powerlaw666http://docs.astropy.org/en/stable/modeling/ (BPL) and compare the best-fits. Using the AIC as our model comparison tool, we find a lower value for the AIC for the broken powerlaw fit compared to the single powerlaw fit; this is true for each of the rebinnings (faded lines Fig. 3). This suggests that the AIC prefers a broken powerlaw, despite the penalty for increasing the number of fit parameters. We note that this is not the case for the higher-mass bins, where a single powerlaw is always preferred. For the broken powerlaw the powerlaw slope is $\alpha_{1}=0.04_{-0.03}^{+0.01}$ at low ICM densities, and $\alpha_{2}=0.30_{-0.10}^{+0.09}$ at high ICM densities. For the single powerlaw fit, the reduced chi-squared values range between 1.6 and 2.9 for the various rebinnings with a median of 2.2. In comparison, the reduced chi-squared for the broken powerlaw ranges between 0.5 and 1.5 with a median of 0.8.

This apparent BPL trend for low-mass galaxies is an intriguing result, and if robust has important implications for the quenching of cluster satellites. Further work is required to test the validity of this result. We devote more time to the discussion of single versus double powerlaw along with other tests of robustness in section 6. However, given the preference for the broken powerlaw according to the AIC and the reduced chi-squared, in sections 4 and 5 we will take the broken powerlaw fit at face value in order to explore a potential origin for the shape of this trend and the implications for satellite quenching.

4 ram pressure stripping model

The broken powerlaw behaviour (for low-mass galaxies) in Fig. 3 matches the qualitative expectation for quenching via ram pressure stripping – where a galaxy remains star forming until reaching a threshold density beyond which the ram pressure force becomes strong and quenching proceeds efficiently (see section 4.1, Fig. 5, Fig. 6). The break point in the broken power-law fit provides an estimate for this threshold density, and for our sample of low-mass galaxies we find a break point of $\log\rho_{\mathrm{thresh}}=-28.3_{-0.7}^{+0.2}\,\mathrm{g}\,\mathrm{cm^{-3}}$. We now can test whether a simple, analytic ram pressure stripping model is able to reproduce the trend for low-mass galaxies in Fig. 3.

To directly constrain the fraction of galaxies susceptible to ram pressure stripping we take a simple analytic approach, similar to models used previously in literature (e.g. Rasmussen et al., 2008; Jaffé et al., 2015, 2018). The basis of the model is the balance between ram pressure and the gravitational restoring force felt by the gas disc in a galaxy (Gunn & Gott, 1972). Specifically, gas will be susceptible to stripping when:

[TABLE]

where $\rho_{ICM}(R)$ is the density of the ICM as a function of cluster-centric radius, $v$ is the galaxy speed relative to the cluster centre, $\Sigma_{\mathrm{gas}}$ is the surface mass density of the atomic+molecular gas component, and $g(r)$ is the maximum restoring gravitational acceleration for the dark matter halo ($\mathrm{DM}$), the stellar disc ($d,\star$), bulge ($b$), atomic gas ($\mathrm{HI}$), and molecular gas ($\mathrm{H_{2}}$) in the direction perpendicular to the disc. Given a model for the ICM density and the galaxy restoring force, the galactocentric radius at which stripping is efficient can be constrained.

As described in section 2.1.1, we estimate the local ICM density for each galaxy using the beta profile fits to the deprojected Chandra density profiles. To complete the left-hand-side of equation 6, we estimate the galaxy speed relative to the cluster centre as,

[TABLE]

where $z$ is the galaxy redshift, $z_{\mathrm{cluster}}$ is the cluster redshift from the Yang catalogue, and $c$ is the speed of light. The factor of $\sqrt{3}$ is included to convert from line-of-sight to three dimensional speed, on average. Following equation 6, the restoring gravitational accelerations for both the gas and stellar components have to be modeled for each satellite galaxy as a function of galactocentric radius. For the stellar distribution the bulge and disc components are modelled with GIM2D bulge + disc decompositions (exponential disc, De Vaucouleurs bulge; Simard et al. 2011), giving disc scale lengths and bulge effective radii ($R_{d},\;R_{e}$) in the $r$-band. We also make use of bulge and disc stellar masses ($M_{b},\;M_{d}$) from Mendel et al. (2014). For the stellar disc, we assume an exponential profile and calculate the restoring gravitational acceleration as

[TABLE]

where $\Sigma_{d,\star}(r)$ is the surface density of the stellar disk as a function of radius. Given the disc mass and scale length, the normalization, $\Sigma_{0}$, can be determined by integrating the surface density,

[TABLE]

For the bulge component, we assume a Hernquist profile, as it has a convenient analytic form which is a good approximation to the De Vaucouleurs profile assumed in the bulge+disc decomposition (Hernquist, 1990). The bulge potential, $\phi_{b}$, is then given by

[TABLE]

where $\mathcal{Z}$ corresponds to the distance in direction perpendicular to the disc component and $a$ is related to the bulge effective radius as $a=R_{e}/1.815$ (Hernquist, 1990). We then calculate the maximum (over the $\mathcal{Z}$ direction) restoring gravitational acceleration as

[TABLE]

For the dark matter halo, we also assume a Hernquist profile. The galaxy halo mass is determined using the observed stellar mass and the stellar-to-halo mass relation from Hudson et al. (2015). The scale radius, $a=r_{\mathrm{vir}}/c$, is estimated from the concentration and the galaxy $r_{200}$, assuming a concentration-mass relation from Diemer & Kravtsov (2015). Both the bulge and dark matter halo profiles are implemented using the colossus222https://bdiemer.bitbucket.io/colossus/ package (Diemer, 2017).

Modeling the gas component is more challenging as, unlike the stellar component, these galaxies are generally not observed in atomic or molecular gas. We assume that the gas distribution consists of an atomic and molecular component, each following an exponential disc with different scale lengths – with the restoring acceleration taking the same exponential form as equation 8. We assume that the atomic gas scale length is twice the optical disc scale length (ie. $R_{\mathrm{HI}}=2\,R_{d}$), which is consistent with the observed value for local non-HI deficient galaxy discs (Cayatte et al., 1994; Cortese et al., 2010; Boselli et al., 2014). For the molecular component we assume a scale length equal to the optical disc scale length (ie. $R_{\mathrm{H_{2}}}=R_{d}$), again corresponding to molecular gas-to-stellar size ratio for non-deficient galaxy discs (Boselli et al., 2014). Because the galaxies in this sample are not observed in atomic or molecular gas, we assign gas masses statistically using the atomic/molecular gas fraction - stellar mass ($f_{\mathrm{gas}}$ vs. $M_{\star}$, $f_{\mathrm{gas}}=M_{\mathrm{gas}}/M_{\star}$) scaling relations for representative samples of galaxies. In particular, we assign HI gas fractions for a given stellar mass assuming the $f_{\mathrm{HI}}-M_{\star}$ relation for isolated galaxies. We are therefore approximating the gas fraction that the cluster galaxies in this sample would have had in the field, prior to infall; allowing us to estimate the amount of pre-infall gas that can be stripped by the cluster environment. To do this we construct the $f_{\mathrm{HI}}-M_{\star}$ relation using the spectral stacking technique (Fabello et al., 2011) applied to a sample of 14128 isolated, low-redshift ($0.02<z<0.05$) SDSS galaxies with processed ALFALFA data (Brown et al., 2017). We direct the reader to Brown et al. (2015) for a complete description of the parent sample and stacking technique, however most importantly, spectral stacking allows us to exploit both HI detections and non-detections to obtain the $f_{\mathrm{HI}}-M_{\star}$ relation for an unbiased sample of galaxies from low ($\sim 10^{9}\,\mathrm{M_{\odot}}$) to high ($\sim 10^{11}\,\mathrm{M_{\odot}}$) stellar mass. Fig. 4 shows the $f_{\mathrm{HI}}-M_{\star}$ relation for the sample of isolated SDSS galaxies for the total sample (black) as well as for subsamples of star-forming (blue, 10984 galaxies with $\log\,\mathrm{sSFR}>-11\,\mathrm{yr^{-1}}$) and passive (red, 3144 galaxies $\log\,\mathrm{sSFR}<-11\,\mathrm{yr^{-1}}$) galaxies. A consequence of the spectral stacking technique is that the intrinsic scatter in the $f_{\mathrm{HI}}-M_{\star}$ relation is unknown, therefore we use the difference between the relation for star-forming and passive galaxies as a rough estimate of the scatter. We assign ‘pre-infall’ HI masses to each cluster galaxy by stochastically sampling from the ‘scatter’ in Fig. 4 (the shaded region) at the stellar mass of the galaxy. The stochastic selection is weighted by the inverse distance from the relation for the total sample (black line), therefore the sampling reflects the fact the scatter is not symmetric about the trend for the total sample. This process is iterated 1000 times and final values calculated using the HI mass are taken to be the median of these Monte Carlo re-samplings.

Pre-infall molecular gas masses ($M_{\mathrm{H_{2}}}$) are assigned using the $f_{\mathrm{H_{2}}}-M_{\star}$ relation from xCOLD GASS (Saintonge et al., 2017) obtained from spectral stacking using a representative sample of low redshift galaxies ($0.01<z<0.02$, $M_{\star}>10^{9}\,\mathrm{M_{\odot}}$). Saintonge et al. (2017) find that, for the stacked sample, the molecular gas fraction is approximately constant at $f_{\mathrm{H_{2}}}\sim 0.1$ for masses $\lesssim 10^{10.5}\,\mathrm{M_{\odot}}$ (see table 5 in Saintonge et al. 2017). For the ram-pressure model we focus on low-mass galaxies ($M_{\star}<10^{10}\,\mathrm{M_{\odot}}$), therefore we choose to assign a constant molecular gas fraction of $f_{\mathrm{H_{2}}}=0.1$ to our cluster galaxies. We note that at these masses the gas content of galaxies remains largely dominated by the atomic component (e.g. Saintonge et al., 2017), which is also the component most suceptible to environmental interactions. Therefore our specific assumptions regarding the molecular gas component do not strongly affect the results.

Using equation 6, combined with the galaxy and ICM models described above we can determine the stripping radius, $r_{\mathrm{strip}}$, the galactocentric radius at which $\rho_{\mathrm{ICM}}(R)\,v^{2}>[g_{\mathrm{DM}}(r)+g_{d,\star}(r)+g_{b}(r)+g_{\mathrm{HI}}(r)+g_{\mathrm{H_{2}}}(r)]\Sigma_{\mathrm{gas}}(r)$. Given the stripping radius, the stripped mass (ie. the gas mass outside of the stripping radius) can be calculated as,

[TABLE]

(Binney & Tremaine, 2008). We compute the stripped mass for both the atomic and molecular gas components and then calculate the total stripped gas fraction as,

[TABLE]

Uncertainties on the disk scale lengths and bulge effective radii from Simard et al. (2011) introduce uncertainties on the stripped masses of roughly 10 per cent, but do not bias the results in any way.

We stress that the ram pressure stripping model described above makes many simplifying assumptions. This model is not intended to be a detailed treatment of ram pressure stripping, but instead to provide a rough estimate of where ram pressure is expected to significantly influence satellite galaxies. We highlight and discuss the primary assumptions of this model in appendix A.

4.1 Comparison to observed quenched fractions

Using this analytic model of ram pressure stripping, we can now make direct comparisons to the observed quenched fractions for the low-mass galaxies in Fig. 3. Specifically, we aim to constrain whether the apparent quenched fraction upturn at high ICM density and low stellar mass can be reproduced by a simple ram pressure model. Therefore, we isolate the high density upturn by subtracting off the powerlaw fit to the low ICM density trend ($\rho_{\mathrm{ICM}}\lesssim 10^{-28}\,\mathrm{g}\,\mathrm{cm^{-3}}$). This allows us to make a direct comparison between the observed upturn and the output from our ram pressure model. We calculate the fraction of galaxies in the sample which have $n$-per cent of their cold gas reserves susceptible to stripping (where $n$ is a free parameter in the model).

In Fig. 5 we overlay the re-scaled quenched fraction, with the low density powerlaw subtracted off, on top of output tracks from the ram pressure model. Each track corresponds to the fraction of galaxies in the sample which have at least the given percentage of their cold-gas mass located beyond the stripping radius (and therefore susceptible to stripping). We show tracks ranging between $>$20 and $>$90 per cent and the observed quenched fraction trend is well matched by a model where $\sim$half of a galaxy’s cold-gas mass is available to be stripped. The stripped fractions in Fig. 3 correspond to the removal of a parcel of gas entirely from the galaxy, as they are computed using the maximum restoring accelerations (in the direction perpendicular to the disc, equation 6). Considering instead the removal of gas from a typical gaseous disc (scale height $\sim$ few hundred parsecs), the data are then well fit by a stripped fraction of $\gtrsim 70$ per cent. Regardless of the precise definition of gas stripping, we emphasize that the primary takeaway is the fact that the observed upturn in quenched fraction coincides closely to the onset of a significant ram pressure force relative to the galaxy restoring potential. The results from this simplified ram pressure model are consistent with low-mass cluster galaxies experiencing enhanced quenching due to relatively efficient ram pressure. Below this threshold density, ram pressure seems to be not strong enough, on average, to quench even low-mass galaxies.

5 Slow-then-rapid quenching

A commonly invoked model for environmental quenching is so-called delayed-then-rapid quenching (Wetzel et al., 2013), where satellite quenching does not occur immediately upon infall but proceeds rapidly only after a characteristic delay time. Qualitatively, the low-mass quenched fraction versus ICM density trend found in this work lends itself naturally to a similar interpretation (see Fig. 3), namely a “slow-then-rapid” quenching framework. Again taking the apparent BPL trend at face value (see section 6 for a detailed discussion of the robustness of this trend), we interpret the modest powerlaw slope between infall and a galaxy reaching the threshold ICM density as the slow-quenching portion, beyond which point the galaxy rapidly quenches as it moves to higher ICM density. This interpretation is illustrated schematically in Fig. 6 both as an annotated quenched fraction plot, and as a diagram showing a toy infall track for a star-forming galaxy onto a cluster.

To be more quantitative, and to permit comparisons to previous studies, we derive rough estimates for the time that a satellite galaxy spends on the slow-quenching track after infall ($t_{\mathrm{slow}}$) as well as the quenching timescale associated with the rapid-quenching component ($\tau_{Q,\mathrm{rapid}}$) using a simple exponential model. The primary simplifying assumption made in this model is that galaxies are quenched exclusively on their first infall on a radial orbit. We note that these assumptions are consistent with recent results from hydrodynamic simulations (Lotz et al., 2018; Arthur et al., 2019).

We calculate the time spent in the slow-quenching mode as,

[TABLE]

where $R_{200m}-R_{\mathrm{thresh}}$ is the radial distance traveled between infall and reaching the threshold ICM density, and $v_{\mathrm{slow}}$ is the mean galaxy velocity over the ‘slow-quenching’ portion. For each cluster, $R_{200m}$ is given by equation 2, and we measure $R_{\mathrm{thresh}}$ given the observed density profile (Fig. 2). To determine galaxy velocities we make use of the high-resolution dark matter only simulations from Joshi et al. (2016). Specifically, we consider galaxy cluster subhalos on their first infall, located between $R_{200m}$ and the median $R_{\mathrm{thresh}}$ for the clusters in our sample, $0.25\times R_{200m}$. We only consider subhalos in $>10^{14}\,\mathrm{M_{\odot}}$ clusters with peak dark matter masses between $10^{11.1}<M_{\mathrm{peak}}<10^{11.75}\,\mathrm{M_{\odot}}$ which corresponds to our low-mass galaxy stellar mass range of $10^{9}\lesssim M_{\star}\lesssim 10^{10}\,\mathrm{M_{\odot}}$ assuming a stellar-to-halo mass relation from Hudson et al. (2015). We normalize the subhalo velocities by the one-dimensional velocity dispersion of the host clusters in the $z=0$ snapshot, which is directly comparable to the measured line-of-sight velocity dispersions for the clusters in our sample. The cluster velocity dispersions measured from the dark matter simulations and the velocity dispersions measured for the observed clusters are both calculated using the biweight estimator (Beers et al., 1990). For subhalos on first infall, located between $0.25<R_{\mathrm{3D}}<1\,R_{200m}$, we find the median velocity is $v=1.8\times\sigma_{\mathrm{1D}}$. To account for the spread in infall velocities, we calculate $t_{\mathrm{slow}}$ in a Monte Carlo sense by sampling velocities, $v_{\mathrm{slow}}$, from the full distribution extracted from the simulations. We measure $t_{\mathrm{slow}}$ as the median of 1000 random samplings. Using equation 14, we obtain an estimate of $t_{\mathrm{slow}}=1.8_{-0.3}^{+0.6}\,\mathrm{Gyr}$. Therefore, infalling galaxies spend $\sim 1.5-2.5\,\mathrm{Gyr}$ before reaching the threshold ICM density, beyond which quenching proceeds rapidly. The value for $t_{\mathrm{slow}}$ is, unsurprisingly, close to the dynamical time for cluster-mass halos.

For the rapid portion we model the quenched fraction as increasing exponentially over a characteristic e-folding time, $\tau_{Q}$. The quenched fraction is then given by,

[TABLE]

Specifically, the e-folding timescale for rapid-quenching is estimated as

[TABLE]

with,

[TABLE]

We take $\rho^{\prime}=10^{-27}\,\mathrm{g}\,\mathrm{cm^{-3}}$ and $f_{Q,\mathrm{rescaled}}(\rho^{\prime})=0.35$, where $f_{Q,\mathrm{rescaled}}$ is the rescaled quenched fraction where the contribution from the low-ICM density powerlaw has been subtracted in order to isolate the high-density upturn (see Fig. 5). The characteristic infall velocity is again determined from dark matter simulations, now for subhalos at $R_{\mathrm{3D}}<0.25\times R_{200m}$. In this inner cluster region we find the median velocity $v=2.8\times\sigma_{\mathrm{1D}}$, and again estimate the timescale by randomly sampling the simulated velocity distributions. Using equation 16, this gives a median quenching e-folding time of $\tau_{Q,\mathrm{rapid}}=0.6_{-0.1}^{+0.1}\,\mathrm{Gyr}$.

Our estimates for $t_{\mathrm{slow}}$ and $\tau_{Q,\mathrm{rapid}}$ suggest a total quenching time of $\sim 2-3\,\mathrm{Gyr}$ for low-mass galaxies in clusters. Studying cluster dwarf galaxies in the Illustris simulation, Mistani et al. (2016) measure the time ellapsed between infall and the first time a galaxy’s sSFR falls below $10^{-11}\,\mathrm{yr^{-1}}$, finding timescales ranging between $\sim 3-5.5\,\mathrm{Gyr}$ for stellar masses $10^{9}\lesssim M_{\star}\lesssim 10^{10}\,\mathrm{M_{\odot}}$. Haines et al. (2015) employ a model where galaxies in LoCuSS clusters are quenched instantaneously after a delay time $\Delta t$ since infall, and find that the surface density of star-forming galaxies based on infrared (UV) observations is best-fit by a delay time of $2.1_{-0.7}^{+0.8}$ ($3.2\pm 0.4$) Gyr. Haines et al. (2015) also derive quenching times of $\sim 1.5-2\,\mathrm{Gyr}$ assuming that star formation declines exponentially upon cluster infall. Wetzel et al. (2013) obtain total quenching times of $4-5\,\mathrm{Gyr}$ for low-mass satellites of cluster-mass halos, which are somewhat larger than the estimates derived in this work. Using semi-analytic models applied to the Millenium Simulation (Springel et al., 2005) along with observations of galaxies in groups and clusters from the Yang et al. (2007) catalogue, De Lucia et al. (2012) argue that on average galaxies spend $5-7\,\mathrm{Gyr}$ in halos $>10^{13}\,\mathrm{M_{\odot}}$ before quenching. These timescales are longer than the total quenching timescale that we derive, however given that our sample is dominated by large clusters (median halo mass, $4\times 10^{14}\,\mathrm{M_{\odot}}$) it is difficult to make a direct comparison. In fact we do see evidence for pre-processing, meaning that many galaxies in our cluster sample were likely members of smaller groups ($\sim 10^{13}\,\mathrm{M_{\odot}}$) prior to cluster infall. Including potential additional time spent as satellites of smaller groups could bring our quenching timescales closer to the De Lucia et al. (2012) estimates. In general, the quenching times that we derive for low-mass cluster galaxies are roughly consistent, if somewhat smaller, than previous estimates from the literature.

6 DISCUSSION

6.1 Pre-processing

In Fig. 3 we show that the quenched fraction in the cluster outskirts is significantly enhanced relative to the field, especially for low-mass galaxies. A natural interpretation for this result is that a fraction of infalling galaxies have been pre-processed prior to infall onto the current cluster. The data for low-mass galaxies require a pre-processed fraction of $\sim 0.3$ to account for the difference from the field value. This pre-processed fraction estimate, however, is quite crude and does not account for any contamination from backsplashing galaxies which will have the effect of artificially increasing the apparent level of pre-processing. Therefore, it is more precise to treat this value as an upper limit. With that in mind, this fraction is roughly consistent with estimates for the pre-processed fraction from previous studies. For example, estimates of the fraction of galaxies which infall onto clusters as members of smaller groups (where the galaxies would be susceptible to pre-processing) range from $\sim\!25$ to $\sim\!60$ per cent using both simulations (McGee et al., 2009; De Lucia et al., 2012; Bahé et al., 2013) and observations (Hou et al., 2014). These constraints provide an upper limit to the pre-processed fraction, as not all galaxies infalling as a member of a group will necessarily have been quenched. More direct constraints on the fraction of pre-processed galaxies range from $\sim\!10-30$ per cent for cluster galaxies (Haines et al., 2015; Roberts & Parker, 2017; van der Burg et al., 2018). The fraction derived in this work falls on the upper of this range, however, it is again important to note that a portion of the satellite population near the virial radius are actually backsplash galaxies that have already made a pericentric passage (e.g. Mahajan et al., 2011; Bahé et al., 2013; Oman et al., 2013; Hirschmann et al., 2014). Therefore it is likely that quenched backsplash galaxies at the virial radius are contaminating the ‘infalling’ population, and artificially increasing the apparent pre-processed fraction.

6.2 Is a broken powerlaw required for the low-mass data?

A key question to address in this study is not only whether the BPL gives a better fit than the SPL to the low-mass data, but also whether or not the SPL provides an acceptable fit (regardless of the quality of the BPL fit). We have already made use of the AIC as a model discriminator, and as is discussed in section 3, the AIC favours the BPL fit over the SPL fit for the low-mass data and all of the re-binnings of these data. This is certainly clear evidence that the BPL fit is statistically preferred over the simpler SPL, however the AIC says nothing about the quality of individual fits, only whether one fit is better than another.

We also test the significance of the single and broken powerlaw fits individually by applying a simple chi-squared test. In Fig. 7 we plot the p-value from the chi-squared test for both the single (crosses) and broken (circles) powerlaw fits for each of the re-binnings, ranging from 10 to 20 bins on the y-axis. For all of the re-binnings we see evidence (at the 90 per cent level) that the SPL does not provide a sufficient fit to describe the low-mass galaxy trend. This suggests that the SPL fit may not be sufficient based on the chi-squared test. One aspect of the fit that a simple chi-squared test does not take into account is any structure in the fit residuals. In Figs 7b and 7c we show the fit residuals (normalized by the uncertainty on each data point) for the single and double powerlaw fits, respectively. We show the residuals for the fits to all of the re-binnings and mark the number of bins in each fit with the colourbar in Figs 7b and 7c. Two trends are clear by comparing Figs 7b and 7c: 1. The amplitude of the residuals are smaller for the BPL fit relative to the SPL. 2. For the BPL fit the residuals seem to be randomly scattered around the zero-line, whereas there is apparent structure in residuals for the SPL fit with a consistent excess of positive residuals between $10^{-31}$ and $10^{-30}\,\mathrm{g}\,\mathrm{cm^{-3}}$ and an excess of negative residuals between $10^{-29}$ and $10^{-28}\,\mathrm{g}\,\mathrm{cm^{-3}}$. This structure in the residuals for the SPL fit is evidence that the SPL model does not fully describe the low-mass quenched fraction data.

6.3 Robustness tests for the broken powerlaw trend

6.3.1 SFR indicator

When calculating quenched fractions we use the MPA-JHU SFRs, which are derived from $\mathrm{H\alpha}$ emission (when detected; Brinchmann et al. 2004). $\mathrm{H\alpha}$ has the advantage of being a tracer of star formation on very short timescales ($\lesssim 10\,\mathrm{Myr}$, Kennicutt & Evans 2012) which is especially important when investigating rapid quenching mechanisms (such as ram pressure). At low SFRs (where emission lines are not detected) the MPA-JHU values are based on the $\mathrm{D_{n}4000}$ break and are therefore less precise. Given the large number of value-added catalogues available for the Sloan Digital Sky Survey, we are able to reproduce Fig. 5 for different star formation rate estimators to test whether the observed shape is driven by our choice of star formation indicator. In Fig. 8 we show the re-scaled quenched fraction (where we have subtracted the low density powerlaw, as in Fig. 5) versus ICM density for SFRs derived from UV+Optical+mid IR SED fitting (Salim et al., 2016, 2018), Optical+IR SED fitting (Chang et al., 2015), rest-frame222k-corrected to $z=0.1$ $u-r$ colours (Blanton & Roweis, 2007) (where we assume that the red fraction corresponds to the quenched fraction), along with the primarily $\mathrm{H\alpha}$ SFRs (Brinchmann et al., 2004) from the main text for reference. For all SFR estimators we define quenched galaxies to have $\mathrm{sSFR}<10^{-11}\,\mathrm{yr^{-1}}$, and for $u-r$ colours we define quenched galaxies to have ${}^{0.1}(u-r)>2.4$ which corresponds to the intersection between the red sequence and blue cloud using a double-gaussian fit. Also plotted for reference are the ram pressure stripping tracks for stripped fractions of $>$40, $>$50, and $>$60 per cent. Fig. 8 shows that the general trend presented in sections 3 & 4.1 persists regardless of SFR indicator, and that in all cases the trend is well matched by a model where quenching becomes efficient when $\sim$half of the galactic cold-gas reservoir is susceptible to ram pressure stripping.

6.3.2 ICM density uncertainty at large radius

Given that the X-ray data used in this work only reach sufficient depth ($>5000$ counts per annulus) within $R_{500}$, it is necessary to extrapolate the density profiles to obtain local ICM density estimates for galaxies in the cluster outskirts. Previous work has shown that ICM density profiles tend to steepen at large radius ($\gtrsim R_{500}$) compared to the inner regions (Morandi et al., 2015). In addition, clumpy gas distributions are common in the cluster outskirts and can introduce biases affecting estimates of the gas density (e.g. Walker et al., 2013; Morandi & Cui, 2014; Ichinohe et al., 2015). This uncertainty in the ICM density profiles at large radius is an important source of uncertainty for this analysis, however it is difficult to quantify on a case-by-case basis without deep X-ray observations out to the virial radius.

Given the observed steepening of density profiles at large radii, by extrapolating profiles we may be overestimating the ICM density in the cluster outskirts. We employ a few different methods to test what effect this could have on the broken power-law trend we observe. As a simple first test we arbitrarily decrease the local ICM density estimate for each galaxy beyond $R_{500}$ by a constant factor. The lines labelled ‘x2’, ‘x5’, and ‘x10’ in Fig. 9 correspond to the quenched fraction trend after assuming a decrease in ICM density by a factor of two, five, and ten beyond $R_{500}$. Note that lines in Fig. 9 have been offset by 0.3 in the vertical direction for readability. For a decrease of a factor of two and five the broken power-law shape is still evident (and is still preferred over a single power-law by the AIC), however for a factor of ten decrease there is no evidence for a broken power-law trend. Therefore unless we are systematically overestimating the local ICM density for each galaxy beyond $R_{500}$ by a factor of $\geq 10$, the broken power-law trend appears robust. The line labelled ‘x1-10’ in Fig. 9 corresponds to decreasing the local ICM density randomly by a factor between one and ten. Specifically, the plotted line shows the median of 1000 random Monte Carlo trials, and the presence of the broken power-law shape is evident and confirmed by the AIC. Finally, Morandi et al. (2015) study a large sample ($>$300) of galaxy clusters observed by Chandra with coverage out to and beyond the virial radius. For their low-redshift sample ($z<0.3$) they find that the slope of the gas density profile ($\beta$) steadily steepens by a factor of 1.75 between $R_{500}$ and $R_{100}$ ($R_{100}$ is similar to the $R_{200m}\simeq R_{\mathrm{vir}}$ that we use). Knowing this, we linearly interpolate between a steepening of a factor of 1 at $R_{500}$ and a factor of 1.75 at $R_{200m}$, and integrate the resulting density profile out to each galaxy’s clustercentric radius to give an updated local ICM density decreased according to Morandi et al. (2015). The result of this is shown with the labelled line in Fig. 9. Again, the broken-powerlaw trend is apparent and is preferred over a single power-law fit to the data by the AIC.

6.3.3 Effect of projected radii

The local ICM density estimates in the paper are determined using observed projected cluster-centric radii and observed deprojected ICM density profiles. Observing galaxy positions in projection is unavoidable, however we can attempt to gauge the effect of observing galaxy positions in projection with the aid of simulations. In particular, observed projected radii are a lower limit to the true cluster-centric radius, and small cluster-centric radii (where ICM densities are highest) is where projection effects are most severe. Given that strong mass segregation is not observed in low-redshift clusters (e.g. Roberts et al., 2015; Kafle et al., 2016), projection effects should not bias the low-mass galaxies differentially relative to the higher-mass bins. Therefore, the fact that we see trends which differ in shape between mass bins suggests that projection effects alone are not driving these differences. However, in the interest of completeness, we explore (in an approximate manner) the influence of projected cluster-centric radii on the BPL trend seen for low-mass galaxies.

To roughly deproject observed, projected cluster-centric radii we make use of dark matter simulations to obtain three-dimensional radii information. In particular, we again make use of high-resolution dark matter only simulations from Joshi et al. (2016) that were used to estimate subhalo infall velocities in section 5. We extract $R_{\mathrm{3D,sim}}/R_{\mathrm{proj,sim}}$ for galaxy-mass subhalos in these simulated galaxy clusters and then measure the median value of $R_{\mathrm{3D,sim}}/R_{\mathrm{proj,sim}}$ in bins of projected cluster-centric radius:

[TABLE]

for these radial bins, the median $R_{\mathrm{3D,sim}}/R_{\mathrm{proj,sim}}$ are:

[TABLE]

When computing projected radii from the simulations we project clusters along a random axis 1000 times and then use median projected radii from these random trials. We perform a linear spline fit to the $R_{\mathrm{3D,sim}}/R_{\mathrm{proj,sim}}$ vs. $R_{\mathrm{proj,sim}}$ relationship and then deproject the observed cluster-centric radius, $R_{\mathrm{proj,obs}}$, as follows:

[TABLE]

where $S(R_{\mathrm{proj,obs}}/R_{200m})$ is the spline fit to the $R_{\mathrm{3D,sim}}/R_{\mathrm{proj,sim}}$ vs. $R_{\mathrm{proj,sim}}$ relationship interpolated to the observed normalized cluster-centric radius, $R_{\mathrm{proj,obs}}$ is the observed physical cluster-centric radius, and $R_{\mathrm{deproj,obs}}$ is the resulting estimate of the deprojected cluster-centric radius.

In Fig. 10 we show the quenched fraction as a function of ICM density, using local ICM densities determined using the deprojected cluster-centric radii outlined above. For low-mass galaxies the BPL shape is still apparent, despite increased scatter at low-densities. In Fig. 10 we also show the best-fit single (gray) and broken (purple) powerlaw fits to the “deprojected” data. The AIC prefers the BPL fit. We emphasize that these are rough deprojections and are only appropriate in an average sense. Unfortunately, obtaining accurate deprojected radii on a galaxy-by-galaxy basis is not possible. A more detailed analysis should include phase space information to incorporate deprojection not only as a function of radial position, but also as a function of velocity offsets. We note that if we deproject only along the z-axis of the simulation, instead of many random halo projections, then the BPL trend in Fig. 10 becomes less apparent and is only marginally preferred over a single power law fit. It is clear that robustly deprojecting observed radii is still an outstanding issue, and we do not rule out that projection effects could contribute to the observed BPL trend.

6.4 A mass dependence of quenching mechanisms?

The two most commonly invoked mechanisms to quench satellite star formation in galaxy clusters are ram pressure stripping and starvation (Gunn & Gott, 1972; Quilis et al., 2000; Wetzel et al., 2013; Muzzin et al., 2014; Fillingham et al., 2015; Jaffé et al., 2015; Peng et al., 2015; Wetzel et al., 2015). It is often argued that the timescale over which these two mechanisms act provides a method of disentanglement. For instance, starvation should quench star formation on relatively long timescales ($\gtrsim 3\,\mathrm{Gyr}$ with the timescale becoming longer for low-mass galaxies) dictated by the gas depletion time, whereas efficient ram pressure stripping will deplete the galaxy of gas, and therefore quench star formation, on much shorter timescales ($\lesssim 1\,\mathrm{Gyr}$ for efficient ram pressure stripping, Quilis et al. 2000; Roediger & Hensler 2005; Steinhauser et al. 2016). However, this argument can be complicated by the fact that ram pressure stripping may require a significant delay before infalling galaxies encounter the densest regions of the ICM, this delay (which should be on the order of the dynamical time of the cluster) can lead to a total quenching time since infall (delay + quenching) which is similar to that of starvation. The quenched fraction trend for low-mass galaxies in this work is consistent with gas depletion (starvation) driving the slow-quenching phase at low ICM density, and ram pressure driving the rapid-quenching phase in the cluster interior. The high-mass galaxies in this sample, though not investigated in detail in this paper, lack the same ram pressure signature. It has been shown that gas depletion times are shortest for high-mass galaxies (Davé et al., 2011; Fillingham et al., 2015; Saintonge et al., 2017), it is therefore plausible that high-mass galaxies consume their gas reserves and quench via starvation prior to reaching the densest cluster interior where ram pressure becomes efficient. We also note that the deeper potential wells of high-mass galaxies will make them more resistent to ram pressure stripping in general. Indeed, Yun et al. (2018) show that the fraction of jellyfish galaxies undergoing strong ram pressure stripping in the Illustris-TNG simulation is strongly dependent on stellar mass, with the jellyfish fraction being highest for low-mass galaxies. Furthermore, the quenching of high-mass galaxies may be largely driven by internal mechanisms, irrespective of environment (Peng et al., 2010). Star formation in low-mass galaxies should persist for much longer after cluster infall, due to the long total gas depletion times, and therefore low-mass galaxies can be still actively forming stars when reaching the densest region of the ICM; where any residual star formation may be quickly quenched due to ram pressure stripping. Substantial ram pressure stripping of atomic gas (the more concentrated molecular component is left largely un-stripped in our models) can disconnect a galaxy from its cold-gas supply, leaving the galaxy to quench via gas depletion (Cen, 2014). Molecular gas depletion timescales for star-forming, low-mass galaxies are $\sim 1\,\mathrm{Gyr}$ (Saintonge et al., 2017), consistent with the rapid quenching timescale we derive. The picture that we present here is consistent with van der Burg et al. (2018) who suggest a quenching scenario where ram pressure is able to “finish the job” when starvation does not quench satellites rapidly enough.

A mass dependent transition between starvation and ram pressure stripping has been previously advocated, where it has been argued that dwarf galaxies ($M_{\star}\lesssim 10^{8}\,\mathrm{M_{\odot}}$) are primarily quenched through ram pressure, whereas the quenching of galaxies with $M_{\star}\gtrsim 10^{8}\,\mathrm{M_{\odot}}$ is dominated by starvation (Fillingham et al., 2015; Wetzel et al., 2015; Rodriguez Wimberly et al., 2018). These conclusions are derived from observations of galaxies primarily in group-mass systems, which are significantly lower mass than the sample of large clusters in this work. Here we find evidence that this transition mass may be larger ($\sim 10^{9}-10^{10}\,\mathrm{M_{\odot}}$) in dense clusters, where both the ICM density and relative velocities are large leading to a strong ram pressure force.

7 SUMMARY

In this work we have used a sample of 24 low redshift SDSS galaxy clusters observed by Chandra to present the first direct study, using a large sample of cluster galaxies, of the relationship between satellite quenching and measured ICM density. The main results of this paper are the following:

Comparing quenched fractions of galaxies at the lowest ICM densities to those for isolated, field galaxies, we find evidence that approximately one third of cluster galaxies may have been pre-proccd essed prior to infall. 2. 2.

The quenched fractions of intermediate- and high-mass cluster galaxies show a modest, continuous increase with ICM density. 3. 3.

The quenched fraction vs. ICM density trend for low-mass galaxies shows evidence of a broken powerlaw trend. The quenched fraction increases modestly at low ICM density, before increasing sharply beyond a threshold ICM density. We show that a broken powerlaw gives a statistically better fit (even after accounting for extra parameters) than a single powerlaw. 4. 4.

The observed broken powerlaw trend is still apparent after observed cluster-centric radii are deprojected using galaxy cluster dark matter simulations, however the strength of the broken powerlaw trend shows some dependence on how galaxy positions are deprojected. We do not rule out that projection effects may contribute to the observed trend. 5. 5.

The quenched fraction upturn at high ICM density, for low-mass galaxies, is well matched by a simple analytic model of ram pressure stripping, where quenching is efficient when more than $\sim$half of a galaxy’s cold gas reservoir becomes susceptible to stripping. 6. 6.

These results are consistent with a slow-then-rapid picture of satellite quenching. We argue that the slow-quenching portion is consistent with quenching via steady gas depletion (starvation) and the rapid-quenching portion is consistent with ram pressure stripping “finishing off” the quenching of low-mass satellites.

IDR, LCP, JHL, and JW are supported by the National Science and Engineering Research Council of Canada. We thank Adam Muzzin for his helpful comments on an early draft, we also thank Katy Rodriguez Wimberly for sharing the colour scheme used in figures throughout the manuscript. The scientific results reported in this article are based on data obtained from the Chandra Data Archive. This research has also made use of software provided by the Chandra X-ray Center (CXC) in the application packages CIAO and Sherpa. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University. .

Appendix A Ram pressure model assumptions & considerations

The ram pressure model described in the manuscript makes a number of simplifying assumptions. Below we highlight and briefly discuss the primary assumptions that went in to the ram pressure model.

• $-$

Both the spectral deprojection software and the fact that we fit to azimuthally averaged density profiles, assume that the clusters in the sample are spherically symmetric. The majority of clusters in this sample do indeed show relaxed, symmetric X-ray morphologies, and if we exclude the 5/24 clusters (containing $\sim 15$ per cent of the low-mass galaxy sample) with clear signs of disturbed morphologies we find that the observed trends are unchanged.

• $-$

We assume that the fits to the cluster density profile are valid out to the virial radius. Previous work (Morandi et al., 2015) has shown that ICM density profiles for galaxy clusters tend to steepen beyond $R_{500}$, therefore we may be overestimating the local ICM density for galaxies at $R>R_{500}$. We address this point in detail in section 6.3.2.

• $-$

We assume that the atomic gas component is distributed in an exponential disc. While this is a common assumption, previous work has suggested that Hi profiles may in fact flatten at small radii (Wang et al., 2014a) and flare at large radii (Kalberla & Kerp, 2009; O’Brien et al., 2010). We employ an exponential profile due to the ease of constructing the analytic surface density given disc masses and scale lengths, however this is at the cost of a more realistic atomic gas distribution.

• $-$

By assigning ‘pre-infall’ gas masses (atomic and molecular) based on $z\approx 0$ observations, we are assuming that the gas fractions for low-mass isolated galaxies at $z\approx 0$ are representative of typical gas fractions at the redshift of infall. Detailed observations of total gas content for large samples of galaxies out to high redshift are not currently feasible, however semi-analytic work has found no strong evolution in the $f_{\mathrm{gas}}-M_{\star}$ relationship out to at least $z\sim 0.5-1$ (Popping et al., 2014) and the lookback time to these redshifts are well in excess of typical crossing times for clusters in our sample.

• $-$

We are assuming that galaxies have not undergone significant stellar stripping since infall, consistent with simulations which show that stellar mass loss after infall is only a minor effect for cluster galaxies (Mistani et al., 2016; Joshi et al., 2018).

• $-$

We are assuming that the present-day structural properties (ie. the outputs of the bulge+disc decompositions) of the galaxies have not evolved since infall. Given that the fraction of bulge-dominated galaxies is enhanced in dense clusters (e.g. Wilman & Erwin, 2012), the galaxies in our sample may have had smaller bulge-to-total ratios at infall. This being the case would lead to our ram-pressure model overestimating the bulge contribution to the restoring potential of the galaxy.

• $-$

We are assuming that galaxies interact with the ICM face-on, leading to the maximal ram pressure efficiency. In reality, inclined interactions will reduce the efficiency of stripping, in particular for near edge-on interactions (Jáchym et al., 2009). This means that our model overestimates the amount of stripping in this respect, however the edge-on interactions where this effect plays a significant role, are relatively rare.

• $-$

We are assuming that satellite galaxies are quenched exclusively on their first infall. While some infalling satellites will certainly survive over multiple orbits, in particular those on tangential orbits, simulations suggest that most cluster satellites are indeed quenched during first infall (e.g. Lotz et al., 2018).