A Gradual Decline of Star Formation since Cluster Infall: New Kinematic Insights into Environmental Quenching at 0.3 < z < 1.1

Our sample of clusters is a subsample of the clusters identified in the 2500 deg² SPT-SZ (Reichardt et al. 2013; Bleem et al. 2015) and ACT-SZ surveys (Marriage et al. 2011; Hasselfield et al. 2013). We refer the reader to the aforementioned publications for full descriptions of these surveys. In short, these cluster candidates were first detected through their SZ signal (i.e., the SZ signal-to-noise threshold ξ ≳ 4.5) and subsequently followed up with optical/near-IR imaging to confirm the clusters associated with the SZ signal. About 500 clusters (415 and 68 from the SPT and ACT, respectively, with some overlapping clusters between the surveys) have been discovered in the surveys at redshifts z = 0.1–1.5. The cluster masses—M₅₀₀, the mass measured within the radius r₅₀₀ at which the mean density of the cluster is 500 times the critical density at the cluster redshift—are ≳3 × 10¹⁴ M_⊙, nearly independent of redshift due to the selection function based on the SZ effect.

From the surveys, we adopt information about the cluster center, redshift (z_cl), and mass (M₅₀₀) of our 105 sample clusters at 0.26 < z_cl < 1.13 (see Section 3.1 for details about sample selection). Specifically, 99 clusters from SPT and nine from ACT satisfy these criteria, with three clusters appearing in both catalogs. The result is a total sample of 105 SZ-selected galaxy clusters. For the three clusters in both catalogs, we use the SZ information from SPT alone (there are no published scaling relations for a combination of SPT and ACT SZ data, but see Hilton et al. 2018 for a discussion of ACT and SPT mass consistency). In particular, the cluster mass (M₅₀₀) is derived from the SZ signal–cluster mass scaling relation (for 99 SPT clusters measured from Reichardt et al. 2013 and Bleem et al. 2015 and six ACT clusters measured from Hasselfield et al. 2013). For the SPT clusters, the following scaling relation is used:

$$\begin{eqnarray}&&\zeta ={A}_{\mathrm{SZ}}{\left(\displaystyle \frac{{M}_{500}}{3\times {10}^{14}{M}_{\odot }{h}^{-1}}\right)}^{{B}_{\mathrm{SZ}}}{\left(\displaystyle \frac{H(z)}{H(0.6)}\right)}^{{C}_{\mathrm{SZ}}},\end{eqnarray} \tag{ 1 }$$

where A_SZ is the normalization factor, B_SZ is the slope, and C_SZ is the redshift evolution term associated with the Hubble parameter H(z). Here ζ is the unbiased SZ significance ²⁰ associated with the SZ signal-to-noise threshold ξ as follows:

$$\begin{eqnarray}&&\zeta =\sqrt{\langle \xi {\rangle }^{2}-3}.\end{eqnarray} \tag{ 2 }$$

The values of A_SZ, B_SZ, and C_SZ are 4.14, 1.44, and 0.59, respectively, as determined in Reichardt et al. (2013). For the six ACT clusters, we adopt the SZ signal based–cluster mass that is calibrated with the cluster physics model of Bode et al. (2012) in Hasselfield et al. (2013; see Section 3.4 and Table 10 of Hasselfield et al. 2013 for more details).

Further details on the cluster surveys, such as cluster identification and mass estimates, can be found in Bleem et al. (2015) and Hasselfield et al. (2013) for the SPT and ACT clusters, respectively.

We obtain the i-band luminosity (L _i ) of our sample galaxies from the optical follow-up observations of the SPT cluster surveys (Bleem et al. 2015, 2020, and references therein) that are conducted to confirm the clusters associated with the SZ signal. Several different telescopes are used for the follow-up observations (see Table 2 of Bleem et al. 2015; for more recent observations with the PISCO instrument, Stalder et al. 2014, see Section 4 of Bleem et al. 2020), including the Blanco/MOSAIC-II; Magellan/Baade IMACS f-2; Magellan/Clay LDSS3, Megacam, and PISCO; Swope/SITe3; MPG/ESO WFI; and New Technology Telescope/EFOSC2.

While the aperture size of the telescopes is not the same (1–6.5 m), we note that the follow-up campaign is optimally designed to observe clusters with uniformly sufficient depth at both low (z ≲ 0.3) and high (z ≳ 0.75) redshifts. Typically, the observations are required to detect 0.4L* galaxies with 5σ depth, where L* is the characteristic luminosity of galaxies at a given redshift. The majority of the photometry uses the Sloan Digital Sky Survey (SDSS) i-band filter. For clusters observed with the older Johnson–Cousins photometric BVRI system, the filter transformation of the I band into the SDSS i band is applied in Bayliss et al. (2016; see their Section 5.2). As a result, all of our sample galaxies' L _i is measured with respect to the SDSS i-band filter.

We use the i-band luminosity (L _i ) scaled by L*, L _i /L*, as our sample galaxies' photometric brightness. We adopt the i-band L* values computed in the SPT-SZ surveys (High et al. 2010; Bleem et al. 2015). Specifically, L* is derived using the stellar population synthesis model of Bruzual & Charlot (2003), and the stellar population modeling assumes a k-corrected, passively evolving, instantaneous-burst stellar population with a formation redshift of z = 3. The Salpeter initial mass function (Salpeter 1955) and the Padova 1994 stellar evolutionary tracks (Fagotto et al. 1994) are used. Metallicities are selected based on analytic fits to the Red-sequence Cluster Survey 2 data (Gilbank et al. 2011). Interpolation using cubic splines is applied to generate L* values at redshifts where the stellar synthesis model did not directly compute L*. The models are further calibrated to the actual SPT spectroscopic subsample. A comparison shows that the computed L* for our sample cluster galaxies is in good agreement with that of the maxBCG cluster sample of Rykoff et al. (2012). Specifically, the L* values between the two studies are consistently better than ≲5% over the overlapping redshift range at z = 0.05–0.35. As we will describe in Section 3.1, we use galaxies brighter than 0.35 L _i /L* as our sample galaxies across all redshifts of interest (0.26 < z < 1.13).

We use the 4000 Å break strength (D_n4000) as a proxy for the luminosity-weighted mean stellar age of galaxies. The D_n4000 is a spectral feature that arises from the accumulation of the metal absorption lines of stars making a "break" between the blue- and red-side continua around 4000 Å in rest-frame wavelength (Bruzual 1983; Hamilton 1985; Balogh et al. 1999). Because the metal absorption line strength is sensitive to stellar type (i.e., age given surface gravity and temperature), D_n4000 is sensitive to the luminosity-weighted mean stellar age of galaxies. Thus, D_n4000 has been widely adopted as a useful stellar age indicator of galaxies (Kauffmann et al. 2003, 2004; Hernán-Caballero et al. 2013; Geller et al. 2014; Haines et al. 2017; Kim et al. 2018, and reference therein).

The D_n4000 is defined as the average flux density (〈F〉) ratio between the blue-side (〈F^B 〉) and the red-side (〈F^R 〉) continua centered at 4000 Å, such that

$$\begin{eqnarray}&&\langle F\rangle =\displaystyle \frac{{\int }_{{\lambda }_{1}}^{{\lambda }_{2}}{F}_{\nu }\ d\lambda }{{\int }_{{\lambda }_{1}}^{{\lambda }_{2}}d\lambda };\end{eqnarray} \tag{ 3 }$$

thus, D_n4000 is expressed as

$$\begin{eqnarray}&&{{\rm{D}}}_{{\rm{n}}}4000=\displaystyle \frac{\langle {F}^{R}\rangle }{\langle {F}^{B}\rangle }=\displaystyle \frac{({\lambda }_{2}^{B}-{\lambda }_{1}^{B}){\int }_{{\lambda }_{1}^{R}}^{{\lambda }_{2}^{R}}{F}_{\nu }\ d\lambda }{({\lambda }_{2}^{R}-{\lambda }_{1}^{R}){\int }_{{\lambda }_{1}^{B}}^{{\lambda }_{2}^{B}}{F}_{\nu }\ d\lambda },\end{eqnarray} \tag{ 4 }$$

where (${\lambda }_{1}^{B}$, ${\lambda }_{2}^{B}$, ${\lambda }_{1}^{R}$, ${\lambda }_{2}^{R}$,) = (3850, 3950, 4000, 4100) Å. The D_n4000 increases with galaxy age such that galaxies with young stellar populations (≲1 Gyr) show D_n4000 ≲ 1.5, while galaxies with old stars (≳1 Gyr) show D_n4000 ≳ 1.5 based on simple stellar population modelings (see, e.g., Figure 2 of Kauffmann et al. 2003).

We obtain the D_n4000 of our sample galaxies measured from the optical follow-up spectroscopic observations of SPT-SZ (Ruel et al. 2014; Bayliss et al. 2016, 2017, and references therein) and ACT-SZ (Sifón et al. 2013) clusters. Most clusters are observed with the Gemini Multi-Object Spectrograph on Gemini South, IMACS on Magellan/Baade, or FORS2 on the Very Large Telescope. The observations are primarily designed to measure spectroscopic redshifts of galaxies for accurate cluster memberships through the cluster velocity dispersion (σ_cl).

The specific observing strategy (e.g., target selection, multislit mask design, and the choice of grating and filter) is set up to uniformly observe a large number (≳100) of clusters. The spectroscopic integration times for individual masks were computed to ensure that the signal-to-noise ratio is ≳3 spectral pixel^–1 in the continuum around D_n4000 based on a model passive galaxy spectrum with a brightness equal to 0.4 L _i /L* at the cluster redshift. Note that we quite consistently sample this luminosity range (see the right panel of Figure 1).

Zoom In Zoom Out Reset image size

The target selection is based on the color–magnitude diagram of galaxies observed in the field of view of a given cluster. The highest priority is given to cluster galaxies identified by the red sequence, regardless of luminosity, down to ≃0.4 L _i /L*. We first identified the red sequence as an overdensity in color and then fit a tilted red sequence to the data in color–magnitude. Likely red-sequence galaxies were those galaxies within ±0.15 mag (in color) of the best-fit red sequence. The choice of ±0.15 mag corresponds to ±2.5–3σ_RS, where σ_RS is the intrinsic width of the observed red sequence, which has been measured to be ≃0.05 mag for massive galaxy clusters (Hennig et al. 2017). The next highest priority is given to candidate "blue cloud" galaxies, which are identified as those galaxies that are bluer than the red sequence (i.e., star-forming). We also only include galaxies that are fainter than the brightest red-sequence member to filter out foreground galaxies.

Obtaining spectroscopic follow-up of a complete sample of cluster galaxies in hundreds of distant galaxy clusters is not practical given the observational resource cost, so it is not a realistic possibility to measure D_n4000 for every cluster member galaxy in our cluster sample. This limitation prevents us from interpreting our results as an absolute measurement of the D_n4000 of the complete sample of SPT cluster member galaxies. That said, we note that the slit placement strategy that produced our spectroscopic member sample was uniformly applied to all spectroscopic observations, with no systematic evolution in the radial density of the spectroscopic slits or the rest-frame magnitude limits. Because the slit placement strategy was applied uniformly, any systematic D_n4000 difference within our sample that scales with other parameters (e.g., galaxy luminosity, redshift, and environment) should be robust.

Understanding our spectroscopic galaxy sample selection is essential to interpreting our results, so we also perform a direct comparison of our spectroscopic cluster member sample against a published photometric analysis that uses complete cluster galaxy populations. Specifically, we measure the fractions of passive (red sequence) and star-forming (blue cloud) cluster member galaxies in our final spectroscopic catalog in our two redshift bins. We measure passive (star-forming) fractions of ∼70% (30%) and ∼66% (34%) for low-redshift (〈z〉 of 0.41) and high-redshift (〈z〉 of 0.66) clusters, respectively. These fractions are consistent with the observed red/passive galaxy fraction measured for all SPT clusters using photometric data (Hennig et al. 2017; see their Figure 15). The consistency between the passive and star-forming fractions in our spectroscopic catalog and published imaging analyses indicates that the sparse spectroscopic sampling does, on average, recover a representative population of cluster galaxies.

Typically, N ≲ 40 galaxies per cluster are observed. The observed spectra cover a galaxy rest-frame wavelength of ≃3500–5150 Å across redshifts z ≃ 0.25–1.1. The observations have a similar spectral resolution, dR ≃ 5–10 Å, corresponding to R (=λ/d λ) ≃ 500–1200. The observed spectroscopy is used to measure D_n4000 by Bayliss et al. (2016) using the same definition given in Equations (3) and (4).

We refer the reader to the spectroscopic follow-ups of SPT-SZ clusters (Ruel et al. 2014; Bayliss et al. 2016, 2017) for more details about the observing strategy and the D_n4000 measurements.

We select our sample clusters by cross-matching the SZ cluster catalog and the optical photometric and spectroscopic follow-up catalogs described in Section 2. From the cross-match, we obtain the cluster mass (M₅₀₀), galaxy luminosity (L _i /L*), D_n4000, and cluster velocity dispersion (σ_cl). There are 105 clusters and 4089 galaxies in our initial sample with a redshift ranging from 0.26 to 1.13. Of the 105 clusters, 99 of them are SPT clusters, and the remaining six are ACT clusters. The cluster mass range is 3 × 10¹⁴ M_⊙ ≲ M₅₀₀ ≲ 3 × 10¹⁵ M_⊙ with a median mass of 4.97 × 10¹⁴ M_⊙ and an associated 1σ population spread of ±2.24 × 10¹⁴ M_⊙. The mass distribution of our sample clusters is shown in the left panel of Figure 1.

Sampling galaxies with reliable D_n4000 measurements is crucial in our analysis. We apply additional sample selection criteria based on physically reasonable D_n4000 values from reliable measurements, so these cuts remove galaxies whose spectra suffer from data reduction artifacts (e.g., bright sky line residuals), as well as high noise. Specifically, we select galaxies whose D_n4000 satisfies 0.7 < D_n4000 < 2.3, σ(D_n4000)/D_n4000 (relative error) < 0.3, leaving 2709 galaxies. The D_n4000 range selected reasonably covers the observed D_n4000 strengths across all types of galaxies, from young (∼10 Myr) to old (∼15 Gyr) stellar ages based on simple stellar population modeling (e.g., Kauffmann et al. 2003; Hernán-Caballero et al. 2013).

We also apply a luminosity cut of L _i /L* > 0.35 to the remaining sample galaxies to ensure a uniform lower limit of galaxy luminosity at all redshifts. Specifically, we derive the luminosity cut by considering the 10th percentile (i.e., above which 90% of the population lies) of the L _i /L* distribution of our high-redshift (0.53 < z < 1.13) subsamples. This corresponds to 0.35 of L _i /L*, which is very similar to the typical depth (5σ for 0.4L*) of the photometric follow-ups of SPT-SZ clusters (Bleem et al. 2015). We thus remove galaxies with L _i /L* < 0.35, leaving 2125 galaxies. The L _i /L* distribution of our sample galaxies with the luminosity cut (L _i /L* > 0.35) applied is shown in the right panel of Figure 1. The L _i /L* range of the sample galaxies roughly corresponds to $10.3\lesssim \mathrm{log}({M}_{\mathrm{star}}/{M}_{\odot })\lesssim 11.6$, given the typical characteristic stellar mass of $10.8\simeq \mathrm{log}({M}_{\mathrm{star}}/{M}_{\odot })$ (e.g., Adams et al. 2021) at similar redshifts and an I-band mass-to-light ratio of 0.8–2.5 (e.g., McGaugh & Schombert 2014).

Lastly, of 2125 galaxies, 499 galaxies are further removed based on their projected clustercentric radius (r_proj/r₅₀₀) and the normalized line-of-sight peculiar velocity (Δv/σ_cl; Equation (5)). This criterion is applied to select cluster member galaxies with reasonable orbital stages (i.e., ∣Δv∣/σ_cl < 3.5 and r_proj/r₅₀₀ < 3), as we will describe in detail in Section 3.2.

In total, our final sample throughout the paper is 105 clusters and 1626 galaxies at 0.26 < z < 1.13. On average, 15 galaxies cluster^–1 are sampled. The typical uncertainty in cluster mass in our cluster sample is 0.9 × 10¹⁴ M_⊙, and the typical uncertainty in the i-band photometry of cluster member galaxies is 0.06 mag. The sample selection procedure is summarized in Table 1.

We estimate the infall time of galaxies based on their location in the phase space, similar to other studies (Noble et al. 2013, 2016; Pasquali et al. 2019). The phase space is a diagram where the infall stage of cluster galaxies can be kinematically estimated (Figure 2). It uses the clustercentric distance and peculiar velocity normalized by cluster velocity dispersion on its axes. For the projected phase space, the projected clustercentric distance (r_proj/r₅₀₀) and the line-of-sight peculiar velocity normalized by the cluster velocity dispersion (Δv/σ_cl) are employed.

Zoom In Zoom Out Reset image size

Indeed, cluster simulations where the infall stages of galaxies can be traced have demonstrated that galaxies distinctively occupy different locations of phase space depending on their infall time (Mamon et al. 2004; Gill et al. 2005; Mahajan et al. 2011; Haines et al. 2012; Oman et al. 2013; Rhee et al. 2017). Typically, galaxies that were accreted early show a wide range of peculiar velocities (Δv/σ_cl) at small clustercentric distances (r/r₅₀₀) but only populate small peculiar velocities over larger clustercentric distance. In contrast, galaxies that were accreted recently are distributed over ranges of clustercentric distance and peculiar velocity but preferentially distributed following trumpet-like profiles that can be described by lines of constant (r/r₅₀₀) × (∣Δv∣/σ_cl) (see, e.g., Figure 3 of Haines et al. 2012).

Motivated by these phase space trends with infall time in simulations, Noble et al. (2013) observationally utilized the caustic lines—the trumpet-shaped lines satisfying (r_proj/r₅₀₀) × (∣Δv∣/σ_cl) = constant (see also Figure 3)—of the projected phase space to divide galaxies into different infall time stages for a z ∼ 0.9 cluster. Their analysis and later studies (Noble et al. 2016; Pasquali et al. 2019) demonstrated the utility of the caustic lines as a useful infall time proxy for cluster galaxies. In particular, a comparison with cluster simulations in Pasquali et al. (2019) shows that galaxies separated by the caustic lines in the projected phase space indeed show systematically different mean infall times, as described below.

Zoom In Zoom Out Reset image size

We adopt the same definition of infall time zones based on the location of phase space as in Noble et al. (2013), which is as follows: ²¹

• 1.  
  early infall, (r_proj/r₅₀₀) × (∣Δv∣/σ_cl) < 0.1;
• 2.  
  intermediate infall, 0.1 < (r_proj/r₅₀₀) × (∣Δv∣/σ_cl) < 0.4; and
• 3.  
  recent infall, (r_proj/r₅₀₀) × (∣Δv∣/σ_cl) >0.4.

As the names of each infall zone indicate, the early infall zone primarily consists of virialized galaxies that were accreted early. The intermediate infall zone contains a mix of galaxy populations with intermediate infallers, as well as a backsplash population that has passed the first cluster pericenter and is currently outbound. The recent infall zone is mainly populated by galaxies that were recently accreted. The mean infall time of the galaxies is typically ${2.7}_{-1.7}^{+1.4}$, ${4.5}_{-2.2}^{+1.9}$, and ${5.1}_{-2.1}^{+1.8}$ Gyr for the recent, intermediate, and early infall zones, respectively. This is based on the comparison with the cluster simulation at z = 0 in Pasquali et al. (2019). (However, we also note that there is scatter in the mean infall time derived from the projected phase space due to projection effects as discussed in Section 4.4.)

For the projected phase space of the sample clusters, Δv and σ_cl are based on the spectroscopic observations of SPT clusters (Section 2.3). In particular, the line-of-sight peculiar velocity (Δv) is expressed as follows:

$$\begin{eqnarray}&&{\rm{\Delta }}v=c\times \left(\displaystyle \frac{{z}_{\mathrm{gl}}-{z}_{\mathrm{cl}}}{1+{z}_{\mathrm{cl}}}\right),\end{eqnarray} \tag{ 5 }$$

where z_gl, z_cl, and c are the redshifts of a galaxy and a cluster and the speed of light, respectively. The r_proj is derived by measuring the angular separation between the cluster center and a galaxy's projected location in units of R.A. and decl. To obtain the clustercentric radius (r_proj/r₅₀₀), we derive the r₅₀₀ of our sample clusters using the cluster mass (Section 2.1) and assuming a spherical mass density profile,

$$\begin{eqnarray}&&{M}_{500}=\displaystyle \frac{4\pi }{3}{r}_{500}^{3}500{\rho }_{\mathrm{crit}},\end{eqnarray} \tag{ 6 }$$

where ρ_crit is the critical density of the Universe at the cluster redshift.

We select galaxies with ∣Δv∣/σ_cl < 3.5 and r_proj/r₅₀₀ < 3 as infalling cluster galaxies as stated in Section 3.1. Figure 3 shows the projected phase space of the sample galaxies. Galaxies mostly populate regions within ∣Δv∣/σ_cl ≲ 2 and r_proj/r₅₀₀ ≲ 1.5. These regions are closely matched with the inner region of the escape velocity radial profile (dashed lines) of the median mass (4.97 × 10¹⁴ M_⊙) of our sample clusters. This suggests that the bulk of infalling galaxies are gravitationally bound by cluster potential wells. To derive the escape velocity profile, we assume a Navarro–Frenk–White (NFW) dark matter halo density profile (Navarro et al. 1996) and a cluster concentration parameter of 4.

The early and intermediate infall zones are all located within the inner region of the escape velocity profile, while the recent infall zone spans the escape velocity profile. The comparison of infall zones with the escape velocity profile makes intuitive sense in that galaxies that have infallen at earlier times are expected to have virialized earlier than those that have infallen more recently. Our low-z (0.26 < z < 0.53) and high-z (0.53 < z < 1.13) subsamples (i.e., gray and white points, respectively, in Figure 3) show similar phase space distributions.

We now investigate how the D_n4000 (mean stellar age) of galaxies varies with time since infall to study the environmental impacts on the star formation of galaxies. Figure 4 shows D_n4000 versus (r_proj/r₅₀₀) × (Δv/σ_cl) (a proxy for infall time; see Section 3.2 for details). The top panel shows the distribution of galaxies at all redshifts (0.26 < z < 1.13). The colored diamonds indicate the mean D_n4000 of galaxies in the corresponding colored infall time zone, which is derived from 1000 bootstrap realizations accounting for the D_n4000 uncertainties.

Zoom In Zoom Out Reset image size

Notably, there is a continuous increase in the mean D_n4000 of galaxies with the infall time proxy moving from recent to early infall populations. The increase in the mean D_n4000 is from 1.51 ± 0.01 to 1.62 ± 0.01. The net increase (ΔD_n4000) corresponds to an age increase of ∼0.71 ± 0.4 Gyr based on a simple stellar population modeling assuming an instantaneous burst of star formation with solar metallicity from Kauffmann et al. (2003). ²² The quantitative age increase estimated above only describes the relative difference in mean stellar age between recent and early infall populations. This average age difference is estimated from galaxies that span a wide range of redshifts (0.26 < z < 1.13); therefore, the reported age difference does not capture any redshift-dependent age differences between infall zones. We examine in detail the redshift dependence of the age difference between infall zones in Section 4.3.

The continuous increase in the mean D_n4000 of galaxies with infall proxy is also shown by the gradual shift of the D_n4000 histograms of each infall zone in the right panels of Figure 4. That is, early infall galaxies (red solid line) tend to have a larger fraction of high-D_n4000 (D_n4000 ≳ 1.5) galaxies compared to the recent infall counterparts (blue dashed line), and the intermediate infall galaxies show the intermediate distribution between recent and early infall zones. This steady shift of the distributions toward larger D_n4000 strength from recent (blue) to early (red) infall populations suggests that the mean age of cluster galaxies increases with time since infall, although there is a broad distribution of D_n4000 values in each infall zone, which is reflected in the overlapping D_n4000 histograms of different infall populations.

The different D_n4000 distributions between different infall time zones are further supported by the two-sample Kolmogorov–Smirnov (K-S) test (Table 2). Table 2 presents the results of the K-S tests for the null (i.e., false-positive) probability that the two selected distributions are statistically different. The left column of the table shows that the D_n4000 distributions between early infallers and the other two categories of infallers are statistically different, with null probabilities of <0.05.

We also perform a linear fit to the D_n4000 versus infall time proxy relation for galaxies with the following form:

$$\begin{eqnarray}&&{{\rm{D}}}_{n}4000=\alpha \ \mathrm{log}[({r}_{\mathrm{proj}}/{r}_{500})\,\times \,({\rm{\Delta }}v/{\sigma }_{\mathrm{cl}})]+\beta ,\end{eqnarray} \tag{ 7 }$$

where α and β are the slope and intercept of the relation, respectively. We perform 1000 iterations for the fits by accounting for the D_n4000 uncertainties. We find a negative slope (α) of −0.096 ± 0.016 of the relation that is constrained to be negative at high statistical significance. This suggests that the D_n4000 of galaxies increases with lower (r_proj/r₅₀₀) × (Δv/σ_cl), meaning older age with longer time spent in clusters since infall. We note, however, that the quantitative slope that we measure could be subject to small selection biases due to our observational selection effects, specifically that we preferentially target red-sequence (passive) galaxies (see Section 3.1 for details).

The middle panel of Figure 4 shows the distribution of low-redshift (0.26 < z < 0.53) subsamples. The low-z galaxies show similar trends as those of the entire sample (top panel). That is, they also show a constant increase of D_n4000 with infall proxy, with a fitting slope (α) of −0.093 ± 0.020 (Table 3) and statistically different D_n4000 distributions between early infallers and other infallers by the K-S test (middle panel of Table 2). A noticeable difference from the entire sample is that the D_n4000 of low-z galaxies is, on average, larger by ∼0.05, regardless of infall time proxy. This is further shown by the larger intercept value of the fitted D_n4000 versus infall proxy relation in low-z subsamples compared to the entire sample in Table 3 (i.e., 1.54 ± 0.02 versus 1.48 ± 0.01 for low-z subsamples and the entire sample, respectively). The ΔD_n4000 of the mean D_n4000 from recent to early infall is a 0.11 increase from ∼1.56, which corresponds to an ∼0.84 ± 0.6 Gyr age increase.

The bottom panel of Figure 4 shows the distribution of high-redshift (0.53 < z < 1.13) subsamples. Like the other two redshift bins (top and middle panels), the high-z subsamples also show the increasing D_n4000 trends with infall proxy, as shown by the increasing mean D_n4000 and the gradual shift of the D_n4000 histogram toward larger D_n4000 strength when moving from recent to early infallers. We also found a negative slope (α of −0.078 ± 0.025) of the relation for high-z galaxies, which is slightly shallower but still consistent with the other redshift bins within 1σ uncertainties (Table 3).

The K-S test for the D_n4000 distributions of different infall zones of the high-z subsamples does not show statistical significance, as the related null probability that the tested distributions are drawn from the same distribution is larger than 0.05 (right column of Table 2). Nonetheless, like the other redshift bins, the null probabilities systematically vary with infall time zones, such that we see the same qualitative trends in D_n4000 with infall time proxy. These systematic K-S test results with infall zones at least suggest that the D_n4000 distribution of our high-z subsamples is related to the galaxies' infall time proxy. The ΔD_n4000 of the mean D_n4000 from recent to early infall for high-z galaxies is an ∼0.08 increase from ∼1.47, which corresponds to an ∼0.63 ± 0.4 Gyr age increase.

Contrary to the low-z galaxies, the high-z subsamples show, on average, smaller D_n4000 values by ∼0.05 compared to the entire sample at all infall proxies. The smaller D_n4000 value of the high-z subsamples is also seen by their small intercept of the D_n4000 versus infall proxy relation compared to the entire sample and low-z subsamples in Table 3.

This redshift dependence of the D_n4000 strengths across infall proxies is likely attributed to the redshift evolution of the star-forming main sequence, such that the average D_n4000 strength of galaxies at fixed stellar mass decreases with increasing redshift (e.g., Whitaker et al. 2012; Haines et al. 2017; Pandya et al. 2017). We will further discuss the redshift dependence of the D_n4000 distributions of our sample galaxies associated with infall time in Section 4.3.

The continuous galaxy age (D_n4000) increase with infall time proxy (Figure 4) is also consistent with the increasing fraction of quiescent galaxies (f_Q) with infall proxy shown in Figure 5. In the figure, f_Q is calculated as the fraction of galaxies whose D_n4000 is >1.55 in bins of infall proxy. Indeed, the f_Q increases from 0.4 ± 0.05 to 0.7 ± 0.07 while moving from recent to early infall zones across all redshifts (z = 0.26–1.13).

Zoom In Zoom Out Reset image size

These steady increases in D_n4000 and f_Q with infall time proxy would suggest that galaxies that have infallen into clusters earlier are quenched earlier due to longer exposure to environmental effects. Our findings are qualitatively consistent with previous studies for several high-z (z ∼ 1; Noble et al. 2013, 2016; Werner et al. 2022) and local (z < 0.2) clusters (Pasquali et al. 2019; Smith et al. 2019; Upadhyay et al. 2021), as well as recent cluster simulations showing systematic suppression of star formation of galaxies since infall (Wetzel et al. 2013; Oman & Hudson 2016; Rhee et al. 2020; Coenda et al. 2022; Oman et al. 2021). Our time-averaged kinematic analysis spans a wide range of redshift, z = 0.26–1.13, with a sample of clusters that are uniformly selected above an approximately constant mass threshold, ²³ and reveals a remarkably consistent picture of environmental quenching extending out to z ∼ 1.

For instance, environmental mechanisms such as ram pressure stripping (Gunn & Gott 1972; Abadi et al. 1999; Ebeling et al. 2014), starvation (Larson et al. 1980), harassment (Moore et al. 1996), and tidal interactions (Byrd & Valtonen 1990) have been suggested as the ones enabling suppression of star formation of galaxies in clusters. Indeed, a systematic gas stripping process has been observed in the Virgo and A963 clusters based on the projected phase space analysis (Vollmer et al. 2001; Boselli et al. 2014b; Jaffé et al. 2015; Yoon et al. 2017; Morokuma-Matsui et al. 2021), where the level of H i gas stripping is closely related to the infalling stages identified in the phase space (see, e.g., Figure 5 of Yoon et al. 2017). Furthermore, the stripping of other components of galaxies, such as warm dust (Noble et al. 2016), molecular gas (CO; Fumagalli et al. 2009; Boselli et al. 2014a; Lee & Chung 2018), and dust (Cortese et al. 2016; Longobardi et al. 2020), is also reported in cluster environments. These results clearly show that gas stripping is actively occurring in cluster galaxies (likely) through the ram pressure by the hot ICM in the deep gravitational potential well of clusters. This is further supported by the recent cosmological cluster simulation that shows the systematic gas depletion of galaxies by ram pressure stripping since infall (Jung et al. 2018).

It is not unreasonable to expect that such ram pressure stripping has occurred in our sample cluster galaxies as well, considering the similar mass range (and thus gravitational potential well) of our sample clusters with those reported for gas stripping (i.e., 1–8 × 10¹⁴ M_⊙). Indeed, the significant X-ray detection (X-ray luminosity > 10⁴⁴ erg s⁻¹) and the associated hot ICM temperature (>2 × 10⁷ K) are found in our cluster subsamples across all redshifts (McDonald et al. 2014; Bulbul et al. 2019). However, ram pressure stripping should be most effective on low-mass systems (e.g., Janz et al. 2021), and the luminosity cut of our analysis (i.e., L _i /L* > 0.35; Section 3) restricts our sample to moderately massive galaxies, suggesting that we should not be highly sensitive to the effects of the ram pressure mechanism.

A slow quenching process, such as the gradual shutoff of gas supply ("starvation"; see, for example, Larson et al. 1980), seems to be a plausible quenching process driving the gradual age (D_n4000) trend with infall proxy in our sample galaxies (Figure 4). Physically, this process is one that acts via the truncation of the galaxy's circumgalactic medium, rather than a faster mechanism that achieves the rapid removal of interstellar medium gas. A slow starvation quenching process would take effect over the course of several dynamical crossing times after a galaxy has fallen into the galaxy cluster environment, which is in good agreement with our measurement of a small average net age increase from recent to early infall populations (i.e., 0.71 ± 0.4 Gyr) compared to the typical dynamical timescales of 1–3 Gyr of cluster galaxies (e.g., Rhee et al. 2020).

A slow quenching trend is also consistent with the approximately constant slope (α) that we measure for the mean stellar population age versus infall proxy relation (i.e., the D_n4000 versus α log [(r_proj/r₅₀₀) × (Δv/σ_cl)] + β relation; Equation (7)). Specifically, the slope α does not significantly depend on the galaxies' L _i /L* (stellar mass), as both bright (L _i /L* > 1) and and faint (L _i /L* < 1) subsamples have statistically similar α (∼−0.09) at all redshifts (Table 3). The similar α we measure across galaxy luminosity suggests that the environmental quenching mechanism acting on these galaxies is not strongly dependent on either galaxy stellar mass or orbital velocity, which matches the expectation for the starvation effect in cluster environments. Qualitatively, our results are in broad agreement with the known "delayed-then-rapid" quenching processes (e.g., Wetzel et al. 2013; Haines et al. 2015; Gallazzi et al. 2021).

We also note that while the average net age increase (0.7 ± 0.4 Gyr) seen in our sample galaxies appears to be broadly consistent with the quenching age (0.1–1.3 Gyr) by the ram pressure gas stripping in local galaxy clusters (Crowl & Kenney 2008; Boselli et al. 2016), one should be cautious in directly comparing the age increase we measure with the quenching ages from the literature. This is due to the difference in how the quenching-related ages are measured between the studies such that our measurement is based on the average net age increase of stellar populations since infall, whereas the SFR quenching age often measured in the literature is considered to be the time taken for a galaxy to transform from gas-rich and star-forming to a totally quenched passive system.

We also note that there are other environmental mechanisms that are likely at play in our sample galaxies. Of them, preprocessing—i.e., pre-exposure of group-scale environmental effects before infalling to the main cluster (Balogh et al. 2000; Fujita 2004; Han et al. 2018; Lee et al. 2022)—is known to be able to quench the star formation of infalling galaxies even in the cluster outskirts (≳R_vir; Balogh et al. 2000), the region where cluster environmental effects are expected to be insignificant. Indeed, recent observations for high-redshift clusters (z ∼ 1) suggest that the majority of massive galaxies are quenched during infall in the cluster outskirts (1 < R/R₂₀₀ < 3), which might imply preprocessing effects prior to infall (Werner et al. 2022).

While our analysis alone cannot pin down the most likely environmental mechanisms responsible for the quenching of cluster galaxies, we emphasize that the steady D_n4000 increase with infall proxy in Figure 4 strongly indicates the environmental quenching at play over a wide span of redshift (0.26 < z < 1.13). Also, our results are qualitatively consistent with recent work at similar redshifts (0.3 < z < 1.5; Webb et al. 2020; Khullar et al. 2021), where the authors used a stellar population fitting analysis to measure the star formation histories and ages of quiescent cluster galaxies.

In the previous section, we have shown that there is a continuous increase of D_n4000 with infall time proxy, which suggests that galaxies that have infallen earlier are quenched earlier due to environmental effects. In this section, we investigate a galaxy luminosity dependence of the environmental quenching trend. We divide our sample galaxies into faint (sub-L*; L _i /L* < 1) and bright (super-L*; L _i /L* > 1) subsamples based on their i-band luminosity (Section 2.2).

Figure 6 shows the same D_n4000 versus infall time proxy diagram as the top panel of Figure 4 but focusing on the galaxy luminosity dependence by splitting it into faint and bright subsamples. The black diamonds and associated error bars are the bootstrapped mean D_n4000 strength and the 1σ spread of the full sample, same as in the top panel of Figure 4. Orange and blue indicate the same symbols as black but for faint and bright subsamples, respectively.

Zoom In Zoom Out Reset image size

Notably, both the faint and bright subsamples qualitatively show the same environmental quenching trend with infall time, as both subpopulations show the continuous increase of D_n4000 with the infall time proxy. The faint and bright subsamples also have slopes (α) of D_n4000 versus infall proxy relation consistent within 1σ, with α = −0.082 ± 0.022 and −0.108 ± 0.024 for faint and bright galaxies, respectively, as in Table 3. This means that the environmental quenching affects galaxies regardless of galaxy luminosity (stellar mass). ²⁴

However, there is a noticeable difference between the faint and bright subsamples, such that the bright subsamples always have a larger mean D_n4000 strength than their faint counterparts across infall zones. This trend does not arise from any potential dependence of luminosity on infall proxy in the sense that bright galaxies might be preferentially located toward earlier infall zones and thus have a larger mean D_n4000 than the faint counterparts. We show that the luminosity distribution of sample galaxies is similar across all infall proxies in the left panel of Figure 7.

Zoom In Zoom Out Reset image size

Thus, at a fixed infall zone (environmental effect), this quantitative D_n4000 difference between faint and bright subsamples is likely attributed to the mass-dependent "internal" quenching of galaxies (e.g., Gavazzi & Scodeggio 1996; Boselli et al. 2001; Peng et al. 2010; Kim et al. 2016, 2018, and references therein). That is, more massive (bright) galaxies become quenched earlier than their less massive (faint) counterparts. This mass quenching may be attributed to mechanisms such as virial shock heating of infalling gas in massive galaxy halos (e.g., Dekel & Birnboim 2006) and/or secular AGN feedback from the supermassive black holes in the galaxy center (e.g., Choi et al. 2015; Bluck et al. 2022).

The individual effects of environment-related infall time proxy and galaxy luminosity (L _i /L*) on quenching as traced by mean galaxy age are demonstrated in Figure 8. Figure 8 is generated by creating a density map of the spatially averaged D_n4000 value within a grid of infall proxy and luminosity. Each pixelated grid position is 0.025 in L _i /L* by 0.05 in $\mathrm{log}[({r}_{\mathrm{proj}}/{r}_{500})\,\times \,({\rm{\Delta }}v/{\sigma }_{\mathrm{cl}})]$, and within each grid pixel, we measure the mean D_n4000 strength of the nearest 100 sample galaxies. Due to some regions of the grid space being very sparsely sampled, we then apply a 2D Gaussian smoothing kernel with a width of 3 grid pixels to improve the visualization of the resulting density map. Clearly, the mean D_n4000 (age) of the galaxies is a function of both environment and luminosity, as the colored D_n4000 strength continuously increases with infall proxy (y-axis) and luminosity (x-axis). However, when fixing one parameter (i.e., x- or y-axis), the mean D_n4000 strength increases along with the other parameter. This means that a continuous D_n4000 increase along with the y-axis (infall proxy) is due to the environmental quenching effect at any given luminosity. This trend is also consistent with the "knee"-shaped white dashed lines that are the contours of constant D_n4000 strengths tracing the mean D_n4000 distributions of sample galaxies.

Zoom In Zoom Out Reset image size

Therefore, Figures 6 and 8 suggest that environmental quenching impacts galaxies independent of galaxy luminosity, although the exact age (D_n4000) depends on the galaxies' luminosity due to the internal mass quenching mechanism. Our kinematic environmental results through infall time proxy show a consistent picture, where stellar mass and environment are the two main drivers of galaxy quenching (e.g., Peng et al. 2010; Smith et al. 2012; Kawinwanichakij et al. 2017; Sobral et al. 2022).

So far, we have demonstrated that environmental quenching occurs since infall such that the galaxies' D_n4000 continuously increases when moving from recent to early infallers (Figure 4). Also, environmental quenching seems to impact galaxies regardless of their luminosity (Figures 6 and 8). We now control for redshift and luminosity simultaneously to further isolate environmental quenching effects, as both parameters can independently alter the D_n4000 strengths (similarly sSFR) of the galaxies regardless of environment (e.g., Kauffmann et al. 2003; Haines et al. 2017; Kim et al. 2018).

Figure 9 shows the same D_n4000 trends with infall time proxy as in other figures but further divided by specific redshift and luminosity bins.

Zoom In Zoom Out Reset image size

First, in all panels showing different luminosity bins, low-z galaxies (blue points) clearly show larger mean D_n4000 values than their high-z counterparts (red points) across infall zones. At the same time, both redshift bins show a continuous increase of mean D_n4000 with infall time proxy due to the environmental quenching at play at all redshift bins. This systematic redshift dependence of D_n4000 at a fixed infall proxy (environmental effect) and luminosity is likely attributed to the redshift-dependent star-forming main sequence of galaxies (e.g., Whitaker et al. 2012; Sobral et al. 2014; Pandya et al. 2017). As the sSFRs of galaxies for a given stellar mass decrease with decreasing redshift, the D_n4000 strength of galaxies given luminosity is also expected to increase with decreasing redshift. Note that this trend does not arise from any potential dependence of galaxy redshift on infall proxy in the sense that low-z galaxies might be preferentially located toward earlier infall zones and thus have a larger mean D_n4000 than their high-z counterparts. This is shown in the right panel of Figure 7, where the redshift distribution of sample galaxies is mostly similar across all infall proxies, although the median redshift appears to slightly increase with infall proxy within a 1σ population spread (i.e., Δz ∼ 0.1 across infall proxies).

Also, by fixing the redshift bins (i.e., focusing on the same colored points across panels), it is clear that bright (super-L*) galaxies have larger mean D_n4000 values than the faint (sub-L*) galaxies at all infall zones, while both luminosity subsamples show an environmental quenching effect (a continuous increase of D_n4000 with infall time proxy). As noted in Section 4.2, this luminosity dependence of D_n4000 is likely explained by the mass quenching effect in that more massive (bright) galaxies tend to stop star formation earlier than those with less mass (faint; e.g., Peng et al. 2010; Kim et al. 2016, 2018; Haines et al. 2017). By splitting into low- and high-redshift bins, we further show that this luminosity dependence of D_n4000 at a fixed infall proxy exists independent of specific redshift bins.

The trends shown in Figure 9 further suggest that environmental quenching indeed affects galaxies since infall regardless of redshift and luminosity, as clearly shown by the systematic increase of the mean D_n4000 with infall time proxy across all redshift and luminosity bins. This is also supported by the similar gradient (slope α) of D_n4000 versus infall proxy relation for all subsamples (Table 3) that are consistent within the uncertainties. However, it is also evident that the exact D_n4000 strength (age) at a fixed environmental effect shows the systematic dependency on redshift and luminosity due to their individual effects on D_n4000 as discussed above.

Despite the exceptional usefulness of the projected phase space for investigating environmental effects, it is also worth noting that there is a spread in our infall time proxy due to sources of uncertainty, such as projection effects and variations in the orbital parameters (e.g., velocity dispersion anisotropy profiles; see Capasso et al. 2019). These include a difference between the projected 2D clustercentric distance and the actual 3D distance and a mix of infalling and outfalling (e.g., backsplash and mere interlopers) populations in the projected phase space.

Cluster simulations (Gill et al. 2005; Rhee et al. 2017) have demonstrated that the line-of-sight velocity distribution of backsplash galaxies, which have passed the first pericenter and are currently outbound, overlap with the infalling galaxies at a low-velocity regime (≲400 km s⁻¹) in the cluster outskirts (1–2 R₂₀₀; see, i.e., Figure 8 of Gill et al. 2005). These backsplash populations are thus likely mixed with the infalling galaxies in the projected phase space.

As a result, such sources of uncertainty in the projected phase space possibly cause a moderate spread (≳1.5 Gyr) in the mean infall time of each infall zone, based on the comparison with the cluster simulation at z = 0 (see, i.e., Table 1 and Figure A1 of Pasquali et al. 2019; Smith et al. 2019; Rhee et al. 2017, 2020). Nevertheless, such a comparison clearly shows that the mean infall time of each infall zone continuously increases from recent (∼2.7 Gyr) to early (∼5 Gyr) infallers. This suggests that (1) the projected phase space is indeed a powerful tool to study the average properties of cluster galaxies as a function of infall time, and (2) projected radial studies can only be more affected by these projection effects than phase space studies. More importantly, we note that these sources of uncertainty only make the trends in Figures 4, 6, 8, and 9 shallower, which means that the trends would be more significant if there were no such sources of uncertainty.