THE PHASE SPACE AND STELLAR POPULATIONS OF CLUSTER GALAXIES AT z ∼ 1: SIMULTANEOUS CONSTRAINTS ON THE LOCATION AND TIMESCALE OF SATELLITE QUENCHING^*

It is well known that galaxies in high-density environments such as galaxy groups and clusters (i.e., satellite galaxies) exhibit a higher fraction of quiescent galaxies at a fixed stellar mass than more isolated "field" galaxies (i.e., central galaxies). This is true both in the local universe (e.g., Baldry et al. 2006; van den Bosch et al. 2008; Gavazzi et al. 2010; Peng et al. 2010; Wetzel et al. 2012; Rasmussen et al. 2012; Haines et al. 2013) and at z ∼ 1 (e.g., Patel et al. 2009; Cooper et al. 2010; Sobral et al. 2011; Muzzin et al. 2012; Raichoor & Andreon 2012; van der Burg et al. 2013; Mok et al. 2013; Woo et al. 2013; Nantais et al. 2013; Kovač et al. 2014; Tal et al. 2014). While the correlation between galaxy quiescence and environment is well established, and heuristic models that can explain the quenching rates and timescales exist (e.g., Peng et al. 2010; Wetzel et al. 2012, 2013), at present there is little direct observational evidence linking the satellite quenching mechanism to a specific physical process that occurs in clusters/groups such as ram-pressure stripping of cold gas (e.g., Gunn & Gott 1972) or hot gas (e.g., "strangulation"; Larson et al. 1980; Balogh et al. 1999), mergers, or harassment (e.g., Moore et al. 1996). One way to make further progress in identifying the dominant satellite quenching mechanism will be to better constrain both the timescale over which quenching occurs and its location within the cluster/group (e.g., Treu et al. 2003; Moran et al. 2007).

Semi-analytic models are physically motivated and recent works have argued that in order to properly reproduce the fraction of quiescent galaxies as a function of clustercentric radius, long quenching timescales, of the order of 3–7 Gyr, are necessary (e.g., Weinmann et al. 2010; McGee et al. 2011; De Lucia et al. 2012). Taken at face value, these long timescales are difficult to reconcile with observations, where a weak dependence of the specific star formation rates (SSFRs) of star-forming galaxies on environment is found (e.g., Kauffmann et al. 2004; Peng et al. 2010; Vulcani et al. 2010; Muzzin et al. 2012; Wetzel et al. 2013). Indeed, the observations suggest a rapid quenching timescale (e.g., Muzzin et al. 2012; Wetzel et al. 2013; Mok et al. 2013, although see Taranu et al. 2014 for evidence of longer timescales), and this hypothesis is supported by an abundance of poststarburst galaxies found in clusters at higher redshifts (e.g., Poggianti et al. 2004, 2009; Tran et al. 2007; Balogh et al. 2011; Muzzin et al. 2012; Mok et al. 2013; Wu et al. 2014). Interestingly, hydrodynamical simulations predict much faster quenching timescales than semi-analytic models (e.g., McCarthy et al. 2008; Bahé et al. 2013; Cen 2014), and similar to the observations, they do not see the strong dependence of SSFR on environment. It has been argued by Wetzel et al. (2013) that one way to reconcile the long quenching timescales required by semi-analytic models with the short quenching timescales required by observations is to have a "delayed-then-rapid" quenching, where galaxies experience a ∼2–4 Gyr delay after infall into a cluster where they behave as normal star-forming galaxies before quenching in <1 Gyr.

While some progress toward defining the timescale of satellite quenching is clearly being made, there are still few observational constraints on where within the cluster/group the process begins, and this is key information for identifying the physical processes involved (e.g., Treu et al. 2003). Our previous work on z ∼ 1 clusters showed that poststarbursts are more common in the cluster core than in the outskirts (Muzzin et al. 2012); however, most previous studies simply compared the poststarburst fraction between cluster and field. A better way of identifying populations within the cluster is to use the velocity versus position phase space of clusters (e.g., Biviano et al. 2002; Mahajan et al. 2011; Haines et al. 2013; Oman et al. 2013). Recently, Noble et al. (2013) performed such an analysis on the most massive cluster in the GCLASS sample, and showed that this combined space is a more effective way of identifying sub-populations within the cluster than simply using clustercentric radius. In particular, they showed that while properties such as SSFR and D_n(4000) show little dependence on clustercentric radius, once galaxies are separated in phase space there is a dependence of those properties as a function of environment.

In this paper, we continue the phase space analysis approach of Noble et al. (2013) using the full GCLASS sample and examine the location of cluster galaxies at z ∼ 1 in phase space. In particular, we focus on the poststarburst population in order to try and identify the satellite quenching timescale and location at z ∼ 1. Throughout this paper, we assume a cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_m = 0.3, Ω_Λ = 0.7.

Our analysis is based on a spectroscopic sample of 424 cluster galaxies in 9 clusters at z ∼ 1 from the GCLASS survey (see Muzzin et al. 2012). The GCLASS clusters were selected from the 42 deg² SpARCS survey (Muzzin et al. 2009; Wilson et al. 2009) using an optical/IR adaptation of the red-sequence method (Gladders & Yee 2000) that is discussed in Muzzin et al. (2008). The clusters have halo masses between M₂₀₀ ∼ (1.0–20.0) × 10¹⁴ M_☉ (van der Burg et al. 2014; G. Wilson et al., in preparation) which have been inferred from the cluster line-of-sight velocity dispersion (σ_v).

We classify each galaxy in the spectroscopic sample as either star-forming, quiescent, or poststarburst⁷ using the [O ii] emission line and D_n(4000) as diagnostics. Star-forming galaxies are classified as those galaxies with detected [O ii] emission, where the detection limit is ∼1–3 Å equivalent width (EW), depending on signal to noise. Quiescent galaxies are defined as those without detected [O ii] emission. Similar to Muzzin et al. (2012), poststarburst galaxies are defined as the subset of quiescent galaxies that have D_n(4000) < 1.45 (i.e., those that are quiescent but have young stellar populations). This is somewhat different than the more typical EW(H_δ) >3–5 Å definition used for K+A galaxies in many studies (e.g., Dressler et al. 1999; Balogh et al. 1999; Poggianti et al. 2009); however, as discussed in Muzzin et al. (2012), H_δ is a weak line and is difficult to measure consistently in all spectra at z ∼ 1, and a weak D_n(4000) serves as a good proxy for strong H_δ in galaxies without [O ii] emission. The average stacked spectrum of our poststarburst definition does have EW(H_δ) >5 Å (Muzzin et al. 2012, see also Section 5), and strong Balmer lines, so on average our selection of poststarbursts is comparable to the K+A selection criteria.

We note that the poststarbust classification may not be 100% complete for all galaxies with strongly truncated star formation. One population that would be missed are old galaxies that experienced a recent rejuvenation of star formation and then a subsequent truncation of that star formation. If the total stellar mass formed in that event was modest, they will have D_n(4000) > 1.45 and remained classified as quiescent. It is unclear if such a population exists in clusters at z ∼ 1; however, if so, then it would be absent from the current sample.

In the left panel of Figure 1, we plot the velocity versus projected clustercentric radius phase space (hereafter referred to simply as "phase space") for all nine clusters. The clusters are combined by normalizing each relative to its σ_v and R₂₀₀ (see van der Burg et al. 2014, Table 1).

Zoom In Zoom Out Reset image size

Figure 1 shows that there is a segregation in phase space between quiescent galaxies (red triangles) and star-forming galaxies (blue triangles). Quiescent galaxies are typically found at smaller clustercentric radii and lower clustercentric velocities; meanwhile, the star-forming population is more extended in both position and velocity. Some star-forming galaxies are located in the central region of phase space; however, it should be noted that phase space is not "projection free", and hence many of these are likely star-forming galaxies with low velocities at large radii (in three-dimensional-position space) but are in projection of the cluster core in two-dimensional-position space. Even with some projection effects, the segregation between the quiescent and star-forming populations is clear and a similar segregation between these types in phase space has also been seen in lower-redshift cluster samples (e.g., Carlberg et al. 1997; Biviano et al. 2002). What is surprising about Figure 1 is the phase space location of the poststarburst galaxies (green stars with circles). Similar to the quiescent galaxies, these tend to lie at smaller clustercentric radii; however, they typically have higher velocities. Most strikingly, these galaxies seem to avoid the "core" region in phase space, where the majority of the quiescent galaxies are located. Indeed, they appear to form roughly a coherent "ring" structure around the core in phase space, although some of the population does extend as far out as R ∼ R₂₀₀.

The phase space segregation is illustrated in the top right panel of Figure 1, where we plot the number of quiescent and poststarburst galaxies relative to the number of star-forming galaxies in three radial bins (dotted lines in the left panel), and the bottom right panel, where we plot these ratios in three phase space bins (shaded ring-shaped regions in the left panel). The shaded regions in Figure 1 have been arbitrarily defined to enhance the contrast of the poststarbursts in phase space; however, they are similar to the variable-slope "chevrons" that have been shown to correlate with the infall histories of galaxies in N-body simulations (see e.g., Mahajan et al. 2011; Taranu et al. 2014; Oman et al. 2013). Within the main region of phase space, they are also similar to the "trumpet"-shaped curves of constant r/R₂₀₀ × v/σ_v that have been shown to correlate with infall times by Noble et al. (2013) based on simulated accretion histories from Haines et al. (2012). Therefore, while the precise definition of the phase space regions is arbitrary (i.e., rings, chevrons, or curves), all three are quite similar and correlate with time since infall into the cluster, making them physically motivated demarcations.

The right panels of Figure 1 show that the trend for the fraction of quiescent galaxies to increase toward the inner bin is roughly the same in both the phase space and radial bins. In radial bins, the fraction of poststarburst galaxies has a peak in the middle bin; however, in phase space bins there are no poststarbursts in the core region, and the majority are confined to the middle phase space bin, with many being at low radii and high velocities.

In order to test if the three populations have distributions in phase space that are different, and whether the difference is statistically significant, we perform a two-dimensional Kolmogorov–Smirnov test of the distribution (hereafter the "2D KS" test). Comparing the two-dimensional distribution of poststarbusts to the quiescent and star-forming galaxies, we find P values of 0.038 and 0.009, respectively. We therefore reject the null hypothesis that they are likely to be drawn from the same distribution at ∼2σ and 3σ, respectively. This demonstrates that the poststarbursts have a distribution that is distinctive compared to the other galaxy types in phase space.

The coherent distribution of the poststarburst galaxies in phase space suggests that it may be possible to use their distribution to constrain the location and timescale of satellite quenching. In this section, we use a set of dark-matter-only zoom simulations of clusters to test if the phase-space distribution of the poststarbursts can be described using a simple quenching model. We note that this modeling is a simplified interpretation of "satellite quenching" as a whole, and also a simplified analysis of the satellite quenching location and timescale. This analysis is performed with the caveat that there may also be a population of satellites that quench without displaying the poststarburst signature. While such a population can exist, many authors have found that the SSFRs of galaxies do not depend on environment (e.g., Peng et al. 2010; Muzzin et al. 2012; Wetzel et al. 2013), and that poststarbursts are substantially more common in high-density regions compared to low-density regions (e.g., Dressler et al. 1999; Balogh et al. 2011; Muzzin et al. 2012). Our assumption in this modeling is that if there are satellites that quench slowly without displaying the poststarburst signature, they are probably satellites that self-quenched (as in, e.g., Peng et al. 2010), whereas the majority that quench from environmental processes do display the poststarburst signature. As has been noted by other studies, this rapid quenching is necessary to preserve the lack of dependence of SSFRs on environment as well as the observed excess of poststarbursts in high-density environments. Therefore, throughout the rest of the paper we note that when we refer to "satellite quenching," we are specifically referring to the quenching of satellites that occurs as a result of environmental processes. With the current analysis we cannot constrain the timescale for satellites that quench through similar processes as central galaxies, e.g., self-quenching.

Our approach is to assume that satellite quenching may begin at a particular clustercentric radius. We then follow the evolution of cluster subhalos in phase space after they make their first crossing of that radius to test if at some time step (hereafter T) later they resemble the phase space of the observed poststarburst galaxies. If so, this timestep would be indicative of the satellite quenching timescale (hereafter τ_Q). For this simple experiment we use three possible clustercentric radii where quenching could begin: the first time a subhalo passes R = 0.25, 0.5, and 1.0 R₂₀₀. These radii roughly correspond to first passage of the cluster core, first passage of the dense intra-cluster medium (ICM), and first passage of the virial radius. We follow the galaxies in time steps of ΔT = 0.2 Gyr up to 1.1 Gyr after crossing the quenching location. This time is approximately the longest time that we can expect to detect the poststarburst signature in galaxies (e.g., Balogh et al. 1999; Poggianti et al. 1999).

We stress that while these models are based on cosmological simulations and therefore the orbits and infall rates should be correct, they are extremely simplified toy models for the satellite quenching process. As discussed above, the satellite quenching process is likely complex, possibly with multiple timescales, locations, and even a dependence on galaxy properties such as stellar mass; all of which is neglected in the current analysis. Furthermore, some of the poststarbursts found at larger radii are almost certainly falling into the cluster already quenched in the field, and we have not attempted to account for "pre-quenched" galaxies. The goal of this modeling is not to make a comprehensive descriptive model of phase space and quenching, but simply to test if the coherent distribution of the poststarburst galaxies in phase space can be reproduced at all, and if so, what it may imply for the location and timescale for satellite quenching.

For the models we use a set of N-body zoom simulations of four clusters that were first presented in Taranu et al. (2014). We refer the reader to that paper for a detailed description of the simulations. In brief, the clusters were selected for re-simulation from a larger low-resolution simulation covering 512h⁻¹ Mpc³. Each re-simulation has 12.5 million particles with a mass of 6.16 × 10⁸M_☉, and resolves Milky-Way-like halos with ∼1000 particles and Magellanic-Cloud-like halos with ∼30 particles. The simulated clusters have virial masses of M_vir ∼ 0.9–1.8 × 10¹⁵ M_☉ at z = 0, and therefore are comparable to massive clusters such as Coma. Their masses are a factor of a few smaller at z ∼ 1, making them an excellent match in halo mass to the GCLASS clusters.

In order to make a fair comparison between the phase space of the models and the observations, a correction for the incompleteness of the spectroscopic sample must be applied. The GCLASS sample has a high spectroscopic completeness overall; however, there is still a spectroscopic targeting bias, with galaxies near the core being targeted more frequently, and more massive galaxies also being targeted more frequently (see Figure 4 of Muzzin et al. 2012, for the spectroscopic completeness corrections). Rather than correct the observations for completeness, we have made the simulations "incomplete" in the same way as the observations, so as to match the observed phase space. For simplicity, the stellar mass in the simulations has been assigned by assuming a uniform stellar-mass-to-halo-mass ratio of 0.1 and the halo mass is the mass of the subhalo before it is accreted.

In Figure 2, we plot output snapshots of subhalos in the simulation in projected phase space at different timesteps (rows). All subhalos in the simulation are plotted as black points, but those subhalos that first crossed the three quenching radii (the three columns of Figure 2) at T = 0 are plotted as red stars. The different columns then allow us to visually follow the evolution of the subhalos in projected phase space through various time steps after "quenching." We work in projected phase space in order to match the observational data. Generally speaking, the marked galaxies in Figure 2 follow similar orbital histories. Most are falling directly into the cluster, likely because they are accreted through filaments in this early stage of cluster formation. Once they pass the chosen radius they continue to increase their velocity and make a close passage of the cluster core (i.e., first pericenter). As would be expected, the passage of the cluster core happens earlier for subhalos that we mark at smaller clustercentric radii (∼0.1 Gyr after crossing 0.25 R₂₀₀ and ∼0.5 Gyr after crossing R₂₀₀). During their first core passage, they tend to be found at high projected velocities and small clustercentric radii. Because of this they are not found in the central "core" region of the cluster phase space at low velocities and small clustercentric radii where many of the massive quiescent galaxies reside in the observations (e.g., Figure 1). Once they make a high velocity crossing of the core they tend to backsplash out to larger radii, as far out as R ∼ R₂₀₀. After another ∼0.5 Gyr they begin to fall back in and appear to be much more mixed in phase space with the rest of the population. In particular, some galaxies manage to penetrate into the low-velocity, small-radius core region in phase space.

Zoom In Zoom Out Reset image size

We make a quantitative comparison between the red stars in the 18 panels of Figure 2 with the poststarbursts in Figure 1 using the 2D KS test. The P values from that test are listed in Table 1. All but four of the snapshots have P < 0.05, and hence we can reject the null hypothesis that they are drawn from the same distribution at ∼2σ. The four radii/timescales that have P > 0.05 and are therefore consistent with being drawn from the same distribution as the poststarbursts are for R = 0.5 R₂₀₀ at T = 0.0 and 0.4, and for R = R₂₀₀ at T = 0.2 and 0.4. In both cases a short timescale is favored, and quenching outside the cluster core is favored. Interestingly, there is no long timestep from any radius where the distribution of subhalos resembles the poststarburst distribution. This appears to rule out the possibility of long quenching times (τ_Q > 0.5 Gyr) no matter what the quenching radius. This may have been expected, as the poststarburst distribution appears coherent, and coherent structures will become mixed and eventually washed out in phase space on the order of a dynamical time.

Formally, the quenching radii and timescales listed above are the most representative of the overall distribution of poststarbursts, which extends out to larger radii; however, from examination of Figure 2 it appears that the quenching at R = R₂₀₀ panels do not reproduce the inner ring structure of the poststarbursts particularly well. They are likely statistically acceptable distributions because they get the overall radial distribution reasonably correct, not the ring structure. As discussed earlier, some of the poststarburst population at larger radius is likely to have been accreted from the field "pre-quenched." The poststarburst fraction in the field at this redshift is ∼1–3%, (Yan et al. 2009; Muzzin et al. 2012), which for ∼400 cluster members implies a maximum of 4–12 galaxies may be infalling already quenched. This number is of the order of the total number of poststarbursts seen at R > 0.5 R₂₀₀ (see Figure 1). If the poststarbursts at larger clustercentric radii are pre-quenched, then it may be that the inner ring of poststarbursts is the primarily signature of satellite quenching, and is the distribution that should be matched in the simulation. Given this, we used the 2D KS test to determine in which timesteps the distribution of subhalos is consistent with the distribution of the poststarburst population using only subhalos and galaxies at R < 0.5 R₂₀₀, where the ring structure is strongest. These results are quite different from the overall 2D KS test. The two timesteps for R = R₂₀₀ that were favored by the full 2D KS test have much smaller P values when only galaxies at R < 0.5 R₂₀₀ are considered (P = 0.069 and 0.058), and are rejected at the ∼2σ level. The only timesteps that are consistent with the observations are for quenching at R = 0.5 R₂₀₀ at T = 0.0, 0.2, and 0.4, and for R = 0.25 R₂₀₀ at T = 0.0. The strong ring structure from these panels can also be seen quite clearly by eye in Figure 2. The reason some of these are not formally the good descriptions of the data in the full 2D KS test is simply because they fail to reproduce the few poststarburst galaxies at larger radii that are seen in the observations. Therefore, this may actually make these the most likely candidates for the location and timescale of the dominant satellite quenching process, even though several of them fail to reproduce the full phase space distribution of phase space.

Ignoring this possibility for the moment, and attempting to match both the inner ring and the overall distribution of the poststarbursts with a single model, we combined a few of the timescales to produce the best-possible description of the observations. This is shown in Figure 3 where we plot the distribution of galaxies quenched at R = 0.5 R₂₀₀ on a timescale of 0.1 <T < 0.5. A 2D KS test applied to that data set provides values of P = 0.358 and 0.185 for all galaxies and galaxies at R < 0.5 R₂₀₀, respectively. This model, with a combined range of timescales has acceptable P values for both the overall distribution and the inner ring structure. We also considered a range of combined timescales for quenching at R = 0.25 R₂₀₀ and R = 1.0 R₂₀₀; however, no combination provided acceptable P values for both the overall distribution and the ring structure simultaneously. Therefore we adopt the model with quenching at R = 0.5 R₂₀₀ on a timescale of 0.1 <T < 0.5 as our best overall description of the data.

Zoom In Zoom Out Reset image size

The phase space modeling suggests that satellite quenching occurs on a timescale of roughly 0.1 Gyr <τ_Q < 0.5 Gyr after galaxies make their first passage of R ∼0.5 R₂₀₀. The simulations imply that the median time it takes for an accreted galaxy to travel from 2.0 R₂₀₀ (1.0 R₂₀₀) to 0.5 R₂₀₀ is ∼0.7 ± 0.2 Gyr (∼0.5 ± 0.7 Gyr), where the uncertainty is the 1σ rms dispersion in infall times. This scenario is qualitatively similar to the "delayed-then-rapid" quenching model proposed by Wetzel et al. (2013), where the quenching time appears to be roughly similar, but the delay time may be a factor of ∼2–4 faster at z ∼ 1. We note that this comparison of the delay time is qualitative, because our delay time is inferred since first passage of 2.0 R₂₀₀, whereas the Wetzel et al. (2013) delay time is measured as the time since a galaxy is identified as a subhalo of a parent halo in an N-body simulation using friends-of-friends linking. It is therefore not immediately clear how these two definitions are related. Although qualitative, it is unlikely that the longer end of the delay time proposed by Wetzel et al. (2013; e.g., 3–4 Gyr) could be supported by the data at z ∼ 1 for two reasons. First, because the simulations show that almost all subhalos would have already made a passage of R = 0.5 R₂₀₀, and hence should be quenched based on our model, and second, because this long delay time is approaching the age of the universe at z = 1 (∼6 Gyr) which would require clusters to form at very high redshift (z > 3).

If the timescales in our phase space model are correct, they place constraints on the ages of the stellar populations of both the poststarburst galaxies and their star-forming progenitors. Poststarburst galaxies should have star formation histories (SFHs) that are consistent with this rapid quenching timescale, and cluster star-forming galaxies should be at least old enough that they have had time to migrate from the cluster outskirts to R = 0.5 R₂₀₀ while maintaining active star formation.

In order to test the consistency of the phase space model and the stellar populations of the galaxies, we fit the spectra of both the star-forming galaxies and the poststarbursts, with the assumption that the cluster star-forming galaxies are the progenitor population of the poststarburst population. Our methodology is as follows: we fit the age of the star-forming galaxies to define the initial age of the poststarburst galaxies once quenching began, defined as t₁. We then fit the poststarburst spectrum starting from this age and fit for two parameters, τ_Q, the timescale over which quenching occurred, and t₂, the age of the galaxy since it was fully quenched. This fiducial SFH is illustrated in the bottom left panel of Figure 4, and is, again, schematically similar to the "delayed-then-rapid" SFH proposed by (Wetzel et al. 2013, see their Figure 12)

Zoom In Zoom Out Reset image size

We begin by stacking the spectra of all star-forming galaxies (details of the stacking process are discussed in Muzzin et al. 2012), with an additional cut requiring that D_n(4000) < 1.45. There are an additional 26 star-forming galaxies with D_n(4000) > 1.45; however, given that this cut is a requirement for the poststarburst selection, it is impossible that star forming galaxies with D_n(4000) > 1.45 can be the progenitors of the poststarburst population. The majority of the galaxies with D_n(4000) > 1.45 that have been excluded have 1.45 < D_n(4000) < 1.6, and so are still young star-forming galaxies, just slightly older than can be allowed to be the progenitors of the poststarburst population. A few have 1.6 < D_n(4000) < 2.0 and clear calcium H+K lines. Most likely these are quiescent galaxies with older population but with an active galactic nucleus which is the source of the [O ii] emission. The top panel of Figure 4 shows the mean stacked spectrum of the star-forming galaxies with prominent absorption features labeled.

To fit the spectrum, we employ the high-resolution models from Bruzual & Charlot (2003) assuming solar metallicity, a Calzetti et al. (2000) dust law (A_v = 0–4), and a Chabrier (2003) IMF. We degrade the spectral resolution of the models to match that of the data which is 17 Å in the observed frame, and corresponds to ∼9 Å rest-frame. We assume a declining SFH based on the decline of the global star-formation rate (SFR) from z = 3 to z = 0 compiled by Bouwens et al. (2012). We also tried an SFH with a constant SFR; however, this continued high level of star formation cannot reproduce the strength of the calcium K absorption line in the star-forming galaxies even at old ages and therefore some form of a declining SFH is required. The Bruzual & Charlot (2003) models do not contain emission lines, so the [O ii] emission line is not fit by the model.

The black solid line in Figure 4 shows the best-fit model to the star-forming population, which has a time since star-formation began of t₁ = 2.3$^{+0.5}_{-0.1}$ Gyr. This implies a mass-weighted age of the stellar population of 1.5$^{+0.2}_{-0.1}$ Gyr. The strongest constraint on the age comes from the relatively strong calcium K line, which cannot be reproduced by a very young stellar population. This intermediate-age population is old enough that it would be continuously forming stars over the infall time from R = 2.0 R₂₀₀ to R = 0.5 R₂₀₀, and therefore is consistent with the delay time implied by the simulation.

To fit the poststarburst population, we create a set of new model grids with the same range of ages and dust extinctions as the declining SFH of the star-forming galaxies. We then add a linear quenching of star formation at t₁ = 2.3 Gyr, and create 11 new grids with τ_Q ranging from 0 to 1.0 Gyr in steps of 0.1 Gyr, where τ_Q is the time between when star formation starts to decline and SFR = 0 (see Figure 4). We then fit the poststarbust spectrum to each of these grids and fit for t₂ (the time since SFR = 0) and A_v. We then find the minimum χ² across the 11 grids which gives us a best fit τ_Q, t₂, and A_v.

The best-fit model for the SFH of the poststarbursts (t₁ = 2.3$^{+0.5}_{-0.1}$ Gyr, τ_Q = 0.4$^{+0.3}_{-0.4}$ Gyr, t₂ = 0.1$^{+0.4}_{-0.1}$ Gyr) is plotted in Figure 4 as the black model on top of the green spectrum. The best fit reproduces the line strengths of the calcium K line, the G-band, and H_δ impressively well. This is remarkable in the sense that calcium K and the G-band are typical of older stellar populations, whereas strong H_δ is from recently quenched galaxies. It would be difficult to reproduce such a spectrum with any form of an SFH other than that in Figure 4.

Remarkably, the implied SFH, τ_Q, and best-fit ages of the stellar populations inferred from the spectral fitting are consistent with the timescales implied by the phase space modeling. It is possible that the SFH and ages could be fit by a wider range of models which we have not explored in exhaustive detail; however, the consistency between the two with basic modeling is encouraging. It is particularly encouraging considering that the stellar population measurements of the satellite quenching timescale are completely independent of the velocities and positions of the galaxies within the cluster.

The overall picture of satellite quenching that arises from the modeling is that star-forming galaxies at z ∼ 1 are not immediately quenched once they are accreted by a cluster/group. Instead, they evolve as normal star-forming galaxies as they fall into the central regions of the cluster. Based on the N-body simulation, it takes ∼0.7 Gyr for a galaxy to first cross R ∼ 0.5 R₂₀₀, and this is where satellite quenching begins. Once the satellite quenching process starts, it proceeds on a short timescale 0.1 <τ_Q < 0.5 Gyr. Moreover, not only is the quenching timescale short, but most of the observed poststarbursts are consistent with having their star formation fully ended only recently. This observation of most poststarbursts having been recently quenched is not a selection effect, as our selection criteria for the poststarbursts of D_n(4000) < 1.45 and [O ii] < 3 Å could select galaxies as old as ∼1.0–1.5 Gyr, depending on the τ_Q.

This scenario and the inferred timescales are remarkably similar to that presented by Smith et al. (2010) for the Coma cluster. They found that galaxies with asymmetric UV morphologies were concentrated within the central 500 kpc (i.e., R ∼ 0.25 R₂₀₀) of the Coma cluster. With numerical modeling they showed that this population could be explained if they were quenched at R ∼ 1 Mpc (i.e., R ∼ 0.5 R₂₀₀) and they were viewable for ∼0.5 Gyr after quenching began. This quenching location and timescale are nearly identical to that derived from the phase space analysis at z ∼ 1.

Our inferred timescale is also consistent with the quenching timescale at z ∼ 1 proposed by Mok et al. (2013) based on their analysis of the fraction of red, blue, and "green" galaxies in groups. They showed that in order to properly reproduce the fractions of these galaxy types, delay times before quenching begins of <2 Gyr were required along with a τ_Q < 1 Gyr. Our quenching timescale is also similar to the quenching timescale proposed by Wetzel et al. (2013) at z = 0 (τ_Q = 0.2–0.8 Gyr); however, it appears that the delay time before quenching may be a factor of 2–4 longer at z = 0.

One potential issue with the proposed scenario for satellite quenching is that we have pre-selected only poststarburst galaxies as our tracer of quenched/quenching population. If there are galaxies that quench on long timescales, they will have been omitted from this analysis. While we cannot formally rule this out, we note that a significant number of studies have shown that the SSFR of star-forming galaxies has little dependence on clustercentric radius at both high and low redshift (e.g., Patel et al. 2009; Vulcani et al. 2010; Muzzin et al. 2012; Wetzel et al. 2013). Although this does not rule out some population of slowly quenching satellite galaxies, it strongly suggests that environmentally quenched satellite galaxies probably quench rapidly, whereas satellites that quench through the same processes as central galaxies could possibly quench slowly.

We do note that some studies of clusters at low-to-intermediate redshift (0.0 <z < 0.5) suggest that there is a weak dependence of the SSFR of star forming galaxies on environment, and this would be consistent with longer quenching timescales than those measured at higher redshift (e.g., Vulcani et al. 2010; von der Linden et al. 2010; Haines et al. 2013. It is not entirely clear why a difference of quenching timescales exists. In a recent work, McGee et al. (2014) compiled quenching timescales from the literature and showed that there does appear to be a redshift dependence of the quenching timescale, with longer timescales at lower redshift. McGee et al. (2014) argue that if this is real, then it may be that because of an change in the dominant physical process causing satellite quenching with redshift. Specifically, it may be that at higher redshift, the high gas outflow rates caused by the higher SSFRs of star-forming galaxies cause them to quench very fast from "overconsumption," i.e., that they lose a significant amount of their gas from outflows that cannot be replaced because cosmological gas accretion has been shut off in groups/clusters. However, McGee et al. (2014) also point out that the measurement of quenching timescales at different redshifts comes from very different data sets, different SFR indicators, and different modeling assumptions, and that some of the redshift evolution may be because of these systematic effects. Similar modeling at multiple redshifts in the future will be essential to understand whether this is the case.

Given the concerns about systematic effects between studies that use different methods, it is remarkable that the quenching timescales determined from four independent methods (Smith et al. 2010; Wetzel et al. 2013; Mok et al. 2013; and the current phase space analysis) are consistent. The phase space constraints are particularly useful because they also provide an additional constraint on where within the cluster/group the satellite quenching begins. With both an inferred location and timescale, we can attempt to infer the physical process that may be responsible.

We first consider the location within the cluster where quenching begins. As derived by Treu et al. (2003) and Moran et al. (2007), quenching from mergers or high-speed galaxy interactions ("harassment") occurs preferentially at R > R₂₀₀, which is inconsistent with our derived quenching radius. Likewise, tidal processes such as halo stripping or disruption occur close to the cluster core (R < 0.25 R₂₀₀), and also seem inconsistent with our data, where the implied quenching radius is R ∼0.5 R₂₀₀. Hot halo gas stripping ("strangulation") and cold gas stripping ("ram-pressure") are most effective where the ICM is dense, at roughly R < 2.0 R₂₀₀ and R < 0.5 R₂₀₀, respectively (e.g., Bahé et al. 2013). Based on the location where we expect quenching to occur, these appear to be the best candidates for the physical process.

If the timescale for quenching is considered at face value, it would suggest that complete removal of the galaxy cold gas via ram-pressure stripping will be necessary to quench galaxies so rapidly. However, Carilli & Walter (2013) have compiled the latest measurements of gas fractions (i.e., f_gas = M_gas/(M_gas+M_star)) of galaxies at z ∼ 1 and these are of the order of 0.3. Our measurements of the SSFRs of galaxies in the GCLASS sample (see Muzzin et al. 2012) show that galaxies with stellar masses of a few times 10¹⁰ M_☉ (the typical stellar mass of the poststarburst galaxies) have log(SSFRs) ∼ −8.8 Gyr⁻¹, which implies that if they were cut off from their hot gas supply completely, they would consume their cold gas in ∼0.5 Gyr. This is consistent with the long end of the quenching timescale that we derive; therefore, it means that we cannot formally rule out hot gas stripping as a plausible mechanism for quenching the poststarbursts in the cluster.

It is also worth considering that the cold gas consumption timescales tend to evolve to longer values at lower redshift. This is simply because SSFRs decline faster than f_gas with decreasing redshift (e.g., Carilli & Walter 2013). This evolution of the cold gas consumption timescale is fairly weak (evolving from ∼0.5 Gyr at z ∼ 1, to ∼1 Gyr at z ∼ 0); however, it is potentially measurable. If cold gas stripping (i.e., ram-pressure) is the dominant satellite quenching process, we might expect that the satellite quenching timescale may not evolve with redshift. However, if hot gas stripping dominates, then the evolution of the cold gas consumption timescale would suggest that the satellite quenching timescale may also evolve to longer values at lower redshift. Therefore, measuring the redshift evolution of the satellite quenching timescale may be a useful approach for identifying whether hot or cold gas stripping is the dominant physical process for satellite quenching. Interestingly, lower-redshift studies such as Wetzel et al. (2013) seem to suggest slightly longer quenching timescales (τ_Q = 0.2–0.8 Gyr), although this has significant overlap with the timescale derived at z ∼ 1 (τ_Q = 0.1–0.5 Gyr). If this could be shown to be statistically significant, then it would imply that hot gas stripping is likely the most important satellite quenching process. However, given the large uncertainties at present, this clearly requires more detailed investigation. It would also benefit from a comparison of timescales measured using similar techniques to avoid systematic errors.

In this paper, we have shown that the population of poststarburst galaxies in clusters at z ∼ 1 has a distribution in phase space that is distinctive from both the quiescent and star-forming cluster galaxy populations. Using a set of dark-matter-only zoom simulations of clusters, we showed that this distribution can be recovered if galaxies quench on a rapid timescale (0.1 < τ_Q < 0.5 Gyr) after first passage of R ∼ 0.5 R₂₀₀. The simulations also show that longer quenching timescales (τ_Q > 0.5 Gyr), or quenching at R ∼ R₂₀₀ provide poor descriptions of the poststarburst phase space distribution. Fitting of the stacked spectra of the star-forming and poststarburst galaxies shows that the SFH, τ_Q, and ages of their stellar populations are consistent with the timescales derived from the phase space analysis, and the similarity between these independent indicators provides the strongest consistency check of the overall model.

The derived quenching location and timescale suggest that gas stripping processes are most likely responsible for quenching the satellite population; however, the current constraints are not strong enough to distinguish between hot or cold gas stripping as the dominant quenching mechanism. A measurement of the evolution of the quenching timescale and possibly the location could be extremely valuable for determining the dominant physical process. Cluster samples with high-quality spectroscopic data exist at lower redshift, and so this is a tractable problem for the future.

An additional conclusion of this work is that the approach of using phase space and simulations to constrain the location and timescale of satellite quenching seems to be a promising new way forward on this problem (see also Noble et al. 2013). We note that better constraints at z ∼ 1 using this approach could be made with larger samples of spectroscopic cluster members. The results of most of the 2D KS tests in this work provide constraints on the timescales at the ∼2σ level, with the limitation being the total number of poststarbursts in the sample (only 28 galaxies out of 424 cluster members). Increasing the number of spectra of poststarbursts by a factor of a few would allow constraints at ∼3σ level or better, which would be useful for further refining the quenching model. Also, more detailed modeling of the infall process that includes tracking self-quenching of the infalling field population would be useful for understanding the population of poststarbursts at larger radii, and putting tighter constraints on the quenching timescale. This will be addressed in future papers.

We acknowledge the work of Chris Power and Brad Krane in developing the N-body simulations used in this analysis. A.M. and M.F. acknowledge support from an NWO Spinoza grant. R.F.J.v.d.B. and H.H. acknowledge support from the Netherlands Organisation for Scientific Research grant number 639.042.814. M.L.B. acknowledges support from an NSERC Discovery Grant, and from NOVA and NWO visitor grants that supported his sabbatical at Leiden Observatory, where this work was completed. M.J.H. acknowledges support from an NSERC Discovery Grant. H.H. acknowledges support from ERC FP7 grant number 279396.