© ESO 2021

1. Introduction

Galaxies occupy different positions in the cosmic web, and their observable properties – such as luminosity, color, morphology, and gas content – at the current epoch are direct results of their evolutionary histories in different environments. The distribution of galaxies in the local Universe shows a bimodality between “blue-cloud” and “red-sequence” galaxies (e.g., Strateva et al. 2001; Baldry et al. 2004, 2006). Blue-cloud galaxies are typically gas rich and star forming, with late-type morphologies, while red-sequence galaxies are typically gas poor with little or no star formation and exhibit early-type morphologies. The fraction of red galaxies increases with both galaxy mass and environmental density (e.g., Dressler 1980; Kauffmann et al. 2004; Baldry et al. 2004, 2006). Galaxies in the high-density environment of clusters are therefore more likely to be quiescent as compared with their similar-mass counterparts in the field (Hogg et al. 2004), an observation further supported by several studies that have found cluster galaxies to be relatively atomic gas deficient (e.g., Gavazzi 1987; Fumagalli et al. 2009; Boselli et al. 2014). The density dependence is significant both within a cluster and out to several times the nominal radius of such systems because the fraction of quiescent galaxies increases with the increase in density toward the cluster core (e.g., Balogh et al. 2000; Lewis et al. 2002; Gómez et al. 2003).

Within a cluster, there are various mechanisms at play that lead to the quenching¹ of galaxies. In addition to “mass quenching”, which affects all galaxies irrespective of their local environment (e.g., Larson 1974; Binney 1977; Dekel & Silk 1986; Birnboim & Dekel 2003; Croton et al. 2006; Martig et al. 2009; Fabian 2012; Cicone et al. 2014), cluster galaxies are also subject to “environmental quenching”. Peng et al. (2010) have shown that environmental and mass quenching efficiency are (nearly) independent. Lower mass galaxies are typically star forming, unless they are quenched by their environment, while massive galaxies are more prone to being quenched regardless of their environment. We now discuss various quenching mechanisms.

There are two major types of environmental quenching, depending on the type of interaction: one is due to the interaction between gas in the interstellar medium (ISM) of a galaxy and the hot gas of the intra-cluster (or intra-group) medium (ICM), and the other is due to the gravitational interaction of galaxies with other cluster members, including the (typically central) brightest cluster galaxy. Effective galaxy-galaxy interactions are rare in clusters due to the high relative speeds of satellites (see Binney & Tremaine 2008, Ch. 8, especially Sect. 8.2g), but repeated high-speed encounters between satellites, termed “galaxy harassment”, can dynamically heat them, either directly removing some gas or rendering it more susceptible to other stripping processes via gravitational heating (Moore et al. 1996, 1998). Tidal effects due to the cluster itself are strongest near its core (Merritt 1983; Mayer et al. 2006), where satellites can be tidally stripped of gas, stars, and dark matter.

Some satellite galaxies enter their host cluster more susceptible to quenching because the environment of a previous, lower mass host has caused a decline in their star formation. These satellites are said to be preprocessed (Balogh et al. 1999; Poggianti et al. 1999; De Lucia et al. 2012; Wetzel et al. 2013; Bahé et al. 2013; Taranu et al. 2014; Roberts & Parker 2017; Haines et al. 2018; Smith et al. 2019; Pallero et al. 2019; Pasquali et al. 2019; Rhee et al. 2020). Several studies have found a high fraction (∼50−65%) of quiescent galaxies in the immediate surroundings of clusters, which may be attributable to preprocessing (Wetzel et al. 2013; Bahé et al. 2013; Taranu et al. 2014; Haines et al. 2018; Pallero et al. 2019). The hierarchical nature of the fiducial ΛCDM cosmological model means that the environments around massive clusters preferentially contain smaller groups of galaxies (see Fig. 1 of Wetzel et al. 2013), so preprocessing plays a significant role in the observed quiescent fraction in and around clusters. The fraction of preprocessed galaxies rises with the final cluster mass, and therefore group preprocessing is more prominent in the satellites of massive (log(M_⋆/M_⊙) > 14.5) clusters, such as the Coma cluster (Hou et al. 2014; Roberts & Parker 2017; Smith et al. 2019; Pallero et al. 2019).

Apart from the group preprocessed population in the outskirts of the cluster, the rest of the satellite population starts experiencing environmental effects even before entering the cluster virial² radius. Accretion of fresh gas becomes inefficient, and the satellites’ hot coronal gas – which acts as a fuel source for star formation by replenishing the cold gas disk – is slowly stripped away, causing a “strangulation” (McGee et al. 2014) of the gas supply required for star formation.

The deep gravitational potential well of a cluster results in high orbital velocities of its satellites and heats up the gas in the ICM. As the cluster galaxies pass through the hot ICM at high speed, their neutral and atomic gas collides with the ICM and may be removed via ram-pressure stripping (RPS; Gunn et al. 1972; Abadi et al. 1999; Jáchym et al. 2007). Ram-pressure stripping likely plays a key role in the quenching of satellite galaxies (Jaffé et al. 2015). Since it is proportional to the density of the ICM, and to the square of the orbital speed, ram pressure peaks strongly around a satellite’s pericenter, though in some cases it can effectively strip their ISM out to at least the virial radius (Tonnesen et al. 2007). Roberts et al. (2019) found that an increase in the quenched fraction near the centers of clusters (r < 0.25r_vir) can be attributed to RPS. Several other studies (e.g., Jaffé et al. 2018; Lotz et al. 2019; Maier et al. 2019; Roberts & Parker 2020; Rhee et al. 2020; Oman et al. 2021) also suggest that RPS plays an important role in satellite quenching.

Several quenching processes act on cluster satellites, and their strengths are dependent on the cluster-centric distance of the satellite and the time spent inside the cluster. One possible approach to constrain the dominant processes is to establish the dependence of their star formation histories (SFHs) on their orbital histories (e.g., Wetzel et al. 2013; Oman & Hudson 2016; Rhee et al. 2020; Oman et al. 2021).

The orbital trajectories and phases of satellites are not directly observable but can be constrained using their projected cluster-centric radius and velocity offset along the line of sight: concisely, their “projected phase space” (PPS) coordinates. Using simulations, galaxy orbits around clusters can be traced backward/forward in time from their current position in PPS (e.g., Mamon et al. 2004; Mahajan et al. 2011; Oman et al. 2013, 2021; Wetzel et al. 2013; Oman & Hudson 2016; Pasquali et al. 2019; Rhee et al. 2020). The distribution of galaxies in and around clusters can be loosely decomposed into three populations: infall, backsplash, and virialized. Infalling satellites are on their initial approach toward their host cluster; virialized galaxies have settled into its central regions; and the backsplash population consists of satellites on eccentric orbits and near their first apocenter, after an initial orbit through the host. The mapping between PPS and orbital phase is not one-to-one: for instance, the PPS location of the backsplash population partially overlaps that of the infalling population. This overlap can make it difficult to distinguish between backsplash galaxies quenched by the cluster, and preprocessed infalling galaxies.

Previous studies have used different parameters to describe the orbital histories of satellites. Examples include the time since infall (TSI; e.g., Wetzel et al. 2013; Oman & Hudson 2016; Rhee et al. 2020), the time since (or until) first pericenter (Oman et al. 2021), and dividing the PPS into zones respectively containing primarily recent and ancient infallers (Mahajan et al. 2011; Hou et al. 2014; Noble et al. 2016; Pasquali et al. 2019; Smith et al. 2019).

To obtain a better understanding of the physics of quenching it is also essential to examine how the quenching timescale (t_q), gas depletion timescale (t_depl), and dynamical timescale (t_dyn) evolve with redshift. The dynamical timescale – or crossing time – is likely a key driver of the redshift evolution of the quenching timescale (e.g., McGee et al. 2014; Foltz et al. 2018; Rhee et al. 2020). It is proportional to (1 + z)^−1.5, or nearly linearly proportional to lookback time. The gas depletion timescale describes the time it would take to use up all the gas to form stars at current star-formation rate (SFR) if the gas is not removed earlier by some other process. It plays a critical role in regulating the gas reservoir for star formation. Using the xCOLD GASS (Saintonge et al. 2017) and xGASS (Catinella et al. 2018) surveys that provide a census of molecular and neutral gas in the galaxies of the local Universe, Feldmann (2020) found that the scaling relation of SFR and total gas mass (molecular+neutral) is directly linked to the gas depletion timescale (see also Tacconi et al. 2018, for the redshift dependence of the gas depletion timescale). Cosmological simulations show that satellite halos stop accreting substantial amounts of mass when they are inside the virialized region of their host halo (e.g., Behroozi et al. 2013b), so when satellites are accreted onto a cluster, the supply of fresh gas is cut off around the time of infall. After this, satellites should continue to form stars from their limited cold gas supply, and hot gas reservoir provided it can cool, unless these are removed from one or many of the gas removal processes discussed in Sect. 1.

Several studies have found longer t_q in clusters at lower redshifts. For example, Wetzel et al. (2013), Taranu et al. (2014), and Balogh et al. (2016) found t_q ∼ 4.4 ± 0.4 Gyr, 4 ± 2 Gyr, and 5 ± 0.5 Gyr, respectively, at z ∼ 0. Haines et al. (2015) found t_q ∼ 3.7 Gyr at z ∼ 0.2. At z ∼ 1, Muzzin et al. (2014), Balogh et al. (2016), Foltz et al. (2018) found t_q ∼ 1 ± 0.25 Gyr, 1.5 ± 0.5 Gyr, and 1.3 ± 0.5 Gyr, respectively. Foltz et al. (2018) studied quenching timescales at an even higher redshift of z ∼ 1.5, finding t_q  ∼  1.1 ± 0.5 Gyr. The longer t_q found by low-redshift studies is suggestive of quenching well after the first pericentric passage (see also Oman et al. 2021), while the shorter timescales at higher redshifts are more consistent with quenching around the first pericenter.

In this study we explore the effects of environment on cluster galaxies. We attempt to constrain the physics driving quenching by comparing the SFHs of satellites with their probable orbital histories. We extract stellar properties from high-S/N spectra of 11 Coma cluster galaxies, and constrain their orbits by comparing with a library compiled from the VVV-L0 N-body simulation extended into the future, up to z = −0.5 (≈10 Gyr from now). In contrast with other studies, we parametrize the orbit by both the expected infall time (when galaxies fall within 2.5r_vir of the cluster) and the time of first pericenter. By determining the fraction of their final stellar mass assembled by the expected infall and first pericenter times, we reconstruct a simple description of the linked orbital and SFHs of our sample.

2. Data and methods

We present the SFHs of our sample of early-type Coma satellites in Sect. 2.1, and the libraries of satellite orbits drawn from an N-body simulation in Sect. 2.2, and how we use these to estimate the orbital parameters of observed satellites in Sect. 2.3.

2.1. STECKMAP star formation histories

We use spectra from Trager et al. (2008), who analyzed the stellar population histories of 12 early-type galaxies (ETGs) in the Coma cluster. The spectra were taken on 7 April 1997 and consist of three consecutive 30-minute exposures taken with the Low Resolution Imaging Spectrograph (LRIS; Oke et al. 1995) on the Keck-II telescope. The spectra have wavelength coverage between 3500–6000 Å with a spectral resolution of 4.4 Å (FWHM) and signal-to-noise (S/N) values ranging between 100–300, except for GMP 3291, GMP 3534, and GMP 3565, which have slightly lower S/N values of 59, 82, and 35, respectively. General properties of our galaxy sample are listed in Table 1. We ignore GMP 3329 in this work as it is a central and not a satellite galaxy in the Coma cluster.

We use STECKMAP (STEllar Content and Kinematics via Maximum A Posteriori, Ocvirk et al. 2006) to extract the stellar and kinematic properties of our galaxy sample. STECKMAP takes the 1D integrated light spectrum of a galaxy as input and computes its stellar population properties and kinematics. The stellar content is inferred by finding the linear combination of model spectra representing single stellar populations (SSPs) that best describes the observed spectrum. The weights of these model spectra give the age distribution of the stellar populations in the galaxy or, in other words, its SFH.

The most important input parameters required for STECKMAP are the wavelength and age range, number of age bins, and the stellar population spectral library. For all galaxies we took the age range from 0.5 to 13.6 Gyr, spanned by 30 bins, and a wavelength coverage of 4050–5500 Å. We chose to use the MILES (Medium resolution INT³ Library of Empirical Spectra) stellar population models, which are based on an extensive empirical stellar library with flux-calibrated spectra of 985 stars over a large range of stellar parameters (Sánchez-Blázquez et al. 2006). It covers a large spectral range (3525–7500 Å) with a spectral resolution of 2.3 Å (FWHM), but more importantly the stellar library is larger and more homogeneous in stellar parameter coverage than the libraries underlying the other models, and for these reasons MILES has become the community standard for nearby galaxy analysis.

The stellar content outputs of STECKMAP are given as a function of age. The basic quantity from which all other stellar content outputs are determined is the stellar age distribution (SAD), which represents the normalized contribution of flux in each component to the observed spectrum. The stellar mass in each time bin is computed from the SAD and M/L(age_i, Z_i) ratio as a function of age and metallicity Z in the given bin:

$$\begin{array}{c}\hfill {\mathrm{mass}}_i=\frac{{\mathrm{SAD}}_i}{M/L({\mathrm{age}}_i,Z_i)}·\end{array}$$(1)

This stellar mass is the initial mass a given component has at the time of its birth because STECKMAP does not consider any mass loss; it only accounts for the dimming of the population and not for the recycling of or decrease in the stellar mass. The normalization of the stellar masses is arbitrary as the SAD is normalized to Σ_iSAD_i = 1. Hence, the SFRs obtained from STECKMAP are not absolute but rather relative SFRs. They are computed from the fraction of the total mass of stars formed by a given time, not the fraction of the final stellar mass (which is less than the total mass of stars formed due to stellar mass loss) formed by that time.

We obtain the (relative) SFR from the stellar mass (mass_i) in each bin by dividing by the duration of age (Δt_i) in each bin (the extent of age in each bin is computed between midpoints of two adjacent bins: Δt_i = t_i + 1/2 − t_{i − 1/2}):

$$\begin{array}{c}\hfill {\mathrm{SFR}}_i=\frac{{\mathrm{SAD}}_i}{\mathrm{\Delta }{t}_i}·\end{array}$$(2)

The relative SFHs in units of Gyr⁻¹ obtained from STECKMAP for our sample of galaxies are shown in Fig. 1. They are color-coded by their stellar mass. As a visual guide, we also show two exponentially decaying rates SFR ∝ exp(τt), with decay timescales representative of early-type (red line, τ = 0.3 Gyr⁻¹) and late-type (blue line, τ = 0.1 Gyr⁻¹) galaxies in low-density environments (Tojeiro et al. 2013). At older ages the SFRs of massive galaxies are higher than those of their lower mass counterparts; this trend reverses at younger ages. We interpret this trend further below.

Fig. 1. Star formation histories of galaxies in our sample determined using STECKMAP, using the MILES SSP model. The SFRs returned by STECKMAP are relative, not absolute, so we show them in units of Gyr−1. The curves are color-coded based on the satellite stellar masses. STECKMAP does not cover lookback times < 0.5 Gyr. The red and blue lines show the exponentially decaying form (SFR ∝ exp(τt)) of SFHs of red (early-type) and blue (late-type) galaxies respectively (Tojeiro et al. 2013).

All of the 11 galaxies in our sample were quenched ≳2 Gyr ago and show no signs of recent star formation (Trager et al. 2008). The slight upturn in SFR (see Fig. 1) for some galaxies at z ∼ 0 is likely to be due to the influence of a subtle effect in STECKMAP where young, metal-rich populations are substituted for old, metal-poor populations due to the age–metallicity degeneracy and missing hot-star populations in the underlying stellar population models (see also Trager et al. 2005).

We computed stellar masses for our sample galaxies by multiplying their B-band luminosities by the color-dependent M/L ratios of Bell et al. (2003), using the B − R colors measured by Eisenhardt et al. (2007). We assumed a distance modulus of 34.94 to Coma and a solar absolute B-band magnitude of 5.51.

2.2. Orbit libraries from cosmological simulations

We inferred orbital parameters for our sample of Coma galaxies by extracting orbital parameter distributions drawn from an N-body simulation, following the general approach of Oman et al. (2021). We used the same initial conditions as the Level–0 simulation described in Wang et al. (2020). The simulation assumes a ΛCDM cosmology with Planck 2014 (Ade et al. 2014) parameters: mean matter density Ω_m = 0.307, cosmological constant Ω_Λ = 0.693, and Hubble parameter at the current epoch H₀ = 66.7 km s⁻¹ Mpc⁻¹. We evolved the initial conditions up to a scale factor a = 2 (z = −0.5, ≈10 Gyr into the future) to predict the time of the first pericentric passage in cases where it has not yet occurred. The simulation has a cube length of 738 Mpc, a particle mass of 1.55 × 10⁹ M_⊙, and a maximum time between consecutive outputs of 390 Myr. The simulation was run using the GADGET4 code (Springel et al. 2021). The simulation was processed by first finding the halos using the ROCKSTAR (Behroozi et al. 2013a) halo finder and then linking the halos in merger trees using the CONSISTENT–TREES (Behroozi et al. 2013b) code.

We selected a sample of “Coma-like” host systems on the basis of their halo mass (details of this selection are given in the next subsection). We label as satellites the halos within a 3D aperture of 2.5 times the virial radius r_vir at z = 0 around the Coma-like hosts. The primary progenitors of the hosts and satellites were traced both forward and backward in time using the merger trees in order to obtain the orbital histories of the satellites. We tabulated the “observable coordinates” of the satellites at z = 0, namely: the projected offset from cluster center normalized by the cluster virial radius (R/r_vir) and the line-of-sight velocity offset from cluster center normalized by the velocity dispersion of dark matter particles in the host halo (V/σ_3D). We arbitrarily chose the third spatial axis (z–axis) of the simulation as that parallel to the line of sight. Here, R and V are given in terms of simulated system coordinates by

$$\begin{array}{c}\hfill \left(\frac{R}{r_{\mathrm{vir}}}\right)=\frac{\sqrt{{(r_{\mathrm{host},\mathrm{x}}-r_{\mathrm{sat},\mathrm{x}})}^2+{(r_{\mathrm{host},\mathrm{y}}-r_{\mathrm{sat},\mathrm{y}})}^2}}{r_{\mathrm{vir}}}\end{array}$$(3)

and

$$\begin{array}{c}\hfill \left(\frac{V}{{\sigma }_{3\mathrm{D}}}\right)=\frac{|v_{\mathrm{host},\mathrm{z}}-v_{\mathrm{sat},\mathrm{z}}|+H(r_{\mathrm{host},\mathrm{z}}-r_{\mathrm{sat},\mathrm{z}})}{{\sigma }_{3\mathrm{D}}},\end{array}$$(4)

where r and v are the coordinate and velocity of the satellites in the simulated system coordinates, the x, y and z subscripts denote their components along orthogonal axes, H is the Hubble parameter, r_vir is the virial radius used to normalize R and σ_3D is the velocity dispersion used to normalize V.

For each satellite, we also tabulated the maximum virial mass at any past time (M_max), the time of first infall through 2.5r_vir (t_inf), and the time of first pericenter (t_peri). We also selected “interlopers” in the foreground and background of each Coma-like halo. These fall within 2.5r_vir in projection at z = 0, but actually are outside 2.5r_vir in 3D. We treated these similarly to the selected satellites but did not trace their orbits, instead using them only to determine an “interloper probability” for each of our observed galaxies.

2.3. Orbital parameter estimation for observed satellites

We then proceeded to select possible orbits for each observed galaxy from the distribution of simulated orbits. We selected these by performing cuts on the library of simulation orbits based on the observed properties of the Coma satellites.

We begin by determining the virial mass and velocity dispersion of the Coma cluster, needed to normalize the PPS coordinates of the satellites. We assume a virial radius of $r_{\mathrm{vir},\mathrm{Coma}}=2.9h_{70}^{-1}\mathrm{Mpc}$ (Łokas & Mamon 2003), which assumes the same definition of the virial radius as that used in our orbit libraries. The corresponding virial mass is M_host, Coma = 1.53 × 10¹⁵ M_⊙, computed from the virial radius as $M_{\mathrm{vir}}=\frac{4\pi }{3}{\mathrm{\Delta }}_{\mathrm{vir}}(z){\mathrm{\Omega }}_{\mathrm{m}}(z){\rho }_{\mathrm{crit}}(z)r_{\mathrm{vir}}^3$, where Δ_vir is the virial overdensity in units of the mean matter density Ω_mρ_crit and ρ_crit = 3H²/(8πG). We use the redshift of GMP 3329 (z = 0.023), which is taken to be at the location of the X-ray center of the Coma cluster⁴, as the cluster’s systemic redshift. We then compute the velocity dispersion of the Coma cluster following (Biviano et al. 2006, see also Bryan & Norman 1998 for the redshift dependence), noting that this velocity dispersion corresponds closely to the dark matter particle velocity dispersion used in the orbit libraries (within about 10%, Munari et al. 2013):

$$\begin{array}{c}\hfill \frac{{\sigma }_{1\mathrm{D}}}{\mathrm{km}{\mathrm{s}}^{-1}}=\frac{0.0165}{\sqrt{3}}{\left(\frac{M_{\mathrm{vir}}}{M_{\odot }}\right)}^{\frac{1}{3}}{\left(\frac{{\mathrm{\Delta }}_{\mathrm{vir}}(z)}{{\mathrm{\Delta }}_{\mathrm{vir}}(0)}\right)}^{\frac{1}{6}}{(1+z)}^{\frac{1}{2}}.\end{array}$$(5)

Finally we compute σ_3D from σ_1D assuming an isotropic velocity dispersion, ${\sigma }_{3\mathrm{D}}=\sqrt{3}{\sigma }_{1\mathrm{D}}$, using the computed value of σ_1D = 1008 km s⁻¹, which gives σ_3D = 1745 km s⁻¹.

We turn our attention next to estimating the virial masses of the satellite galaxies, using the stellar-to-halo mass relation (SHMR) of Behroozi et al. (2010, Eq. (21)):

$$\begin{array}{c}\hfill log(M_{\mathrm{h}}(M_{\star }))=log(M_1)+\beta log\left(\frac{M_{\star }}{M_{\star ,0}}\right)+\frac{{\left(\frac{M_{\star }}{M_{\star ,0}}\right)}^{\delta }}{1+{\left(\frac{M_{\star }}{M_{\star ,0}}\right)}^{-\gamma }}-\frac{1}{2},\end{array}$$(6)

where M_h(M_⋆) is the halo mass for which the average stellar mass is M_⋆, M₁ is the characteristic halo mass (log(M₁) = 12.35 at z = 0), M_⋆, 0 is the characteristic stellar mass (log(M_⋆, 0) = 10.72 at z = 0), and the constants are power-law coefficients for the relation (at z = 0: β = 0.44, δ = 0.57 and γ = 1.56).

The virial masses estimated in this way do not account for tidal stripping of the satellites’ dark matter halos: we interpret these as the virial masses that the satellites had before entering the cluster, and will match them against the maximum mass of the satellite halos in the simulations. This approximation is imperfect for two reasons. First, we use the z = 0 SHMR. We should in principle use the SHMR at the redshift corresponding to the infall time of the satellites, but to estimate this infall time requires us to know their virial masses. Second, we neglect any stellar mass growth or loss between infall and the present. Although we do have SFH information, a self-consistent accounting for the stellar mass growth again requires knowledge of the orbital history, which we are in the process of attempting to estimate. For stellar mass loss, we assume that any tidal stripping of stars has been moderate (see Smith et al. 2016): once satellites are heavily stripped of stars, they are very quickly destroyed (Bahé et al. 2019). We will explore the extent to which the systematic offsets introduced by this approximation may affect our results in Sect. 4.2.

The PPS coordinates of each satellite, R/r_vir and V/σ_3D, are computed from their measured sky coordinates and redshifts as

$$\begin{array}{c}\hfill \frac{R}{r_{\mathrm{vir}}}=\frac{d_{\mathrm{A}}\mathrm{\Delta }\theta }{r_{\mathrm{vir}}}\end{array}$$(7)

and

$$\begin{array}{c}\hfill \frac{V}{{\sigma }_{3\mathrm{D}}}=\frac{c|z_{\mathrm{g}}-z_{\mathrm{c}}|}{(1+z_{\mathrm{c}}){\sigma }_{3\mathrm{D}}},\end{array}$$(8)

where d_A = 99 Mpc is the angular diameter distance of the Coma cluster and Δθ is the angular separation of the satellites from the Coma center (taken to be the angular distance from GMP 3329, as above), and where c is the speed of light, z_g is the redshift of the satellite, z_c is the redshift of the Coma cluster and σ_3D is the 3D velocity dispersion of the Coma cluster (Eq. (5)). The measured and computed properties for satellites of Coma are listed in Table 2 for reference. Figure 2 shows the Coma galaxies in projected phase space. We note that the spread in the radii is very small. The orbits in our N-body libraries are not significantly different between satellites so close together in radius, so differences between the inferred orbits are driven primarily by the satellite velocity offsets (V/σ_3D), and their halo masses.

Fig. 2. Galaxies used in this study shown in projected phase space. The galaxies are colored by log(Msat/M⊙). The points are labeled with their GMP identifiers.

We select a distribution of likely orbits for each satellite by selecting orbits from the simulation that have the same PPS coordinates and host and satellite virial masses, within the following tolerances. For M_host and M_max, the upper and lower limits of the selections are set 0.5 dex above and below the virial mass of the Coma cluster and the satellite, respectively. The limits for the selections around the PPS coordinates R/r_vir and V/σ_3D are set to ±0.05 of their measured values. The cuts resulted in many orbits for satellites that are similar (in terms of mass and current projected coordinates) to the observed galaxies. We then make the reasonable assumption that the selected orbits represent an approximation of the probability distribution of possible orbits for a given satellite. The parameters used in this study are listed in Table 2.

The Coma cluster has the interesting property of having two similarly bright BCGs, GMP 3329 (NGC 4874) and NGC 4889. We find a reasonable number of Coma like clusters (458) in the simulation that could be split into two roughly even groups based on whether they have a very massive satellite near the center. Of these clusters, 126 have a satellite with mass at infall > 1 × 10¹⁴ M_⊙, which might be reasonable matches for the halo mass of NGC 4889 at infall, though this mass is difficult to constrain observationally. However, as the number of matches for each satellite in our sample is rather small (≲100), splitting the cluster sample would strain our ability to construct orbital parameter PDFs for the sample of satellites. In light of these consideration, we do not explore the topic of Coma’s substructure further in this study; analysis of a larger-volume simulation where a large number of double-BCG systems could be analyzed could be pursued in future work. Another possible complication is that the halo mass function is rather steep, and we have used a symmetrical selection criterion of ±0.5 dex from the target values for cluster and satellite halos; however, we shown in Sect. 4.3 below that this steepness makes no appreciable impact on our results.

We make analogous selections from the library of interloper halos. By comparing the number of interlopers selected to the number of satellite orbits, we obtain an estimate of the probability that the satellite is in fact an interloper. These interloper probabilities are included self-consistently throughout the rest of the analysis described below, but they have little influence because they are very small (see Table 2). We therefore do not discuss this point further.

3. Results

3.1. PDFs of orbital time and SFR comparison

We illustrate our approach to compare the SFHs and orbital histories of our sample of Coma satellites in Fig. 3, using GMP 3254 as a representative example. In the left panel, we show the differential probability distribution for the infall time of this satellite (gray histogram), as determined from the N-body orbit library, and the differential SFH (red line).

Fig. 3. Connecting orbital history with SFH. Left panel: probability distribution of infall time of GMP 3254 is shown in the gray histogram (left axis labels) and its SFH (right axis labels) is shown as the red points. The median value of its infall time probability distribution is shown as a dotted black vertical line. The SFR has dropped significantly by ⟨tinf⟩. Right panel: cumulative SFH – the fractional stellar mass accumulation over time – is shown (red line) along with the cumulative probability distribution of the infall time (black line) for GMP 3254. The median infall time (black filled circle) is marked with its 68% confidence interval (horizontal black line) at the 50th and 16th–84th percentiles, respectively. The corresponding accumulated stellar mass at the median infall time (red filled circle) and the its associated 68% confidence interval (vertical red line) are also shown. The arrows (black and red) indicate the translation of 16th and 84th percentile points on the cumulative distribution curve to boundary of the confidence interval lines. GMP 3254 (log10M⋆/M⊙ = 9.92) has formed around 90% of its stellar mass at the median infall time. Figure 4 (left panel) shows the accumulated stellar mass at the (median) infall and pericenter times for all galaxies in this study.

The infall time probability distribution is broad and has two peaks: the younger peak corresponds to the possibility that GMP 3254 is now near its first pericentric passage in Coma, while the older peak corresponds to the possibility that it is near its second pericentric passage (the tail to even older ages corresponds to multiple orbits having been completed). The separation of the two peaks reflects (twice) the crossing time of the Coma cluster t_cross = r_vir/σ_1D ∼ 2.5 Gyr, keeping in mind that the crossing time would be somewhat shorter for earlier infall times. The SFH is monotonically declining toward the present time, and has no clear or sharp features. This makes it difficult to immediately see any obvious link between the SFH and orbital history.

We also computed the differential probability distribution for pericenter time and compared it with the relative SFR, following the same method illustrated in the left panel of Fig. 3 for the infall time. However, it is difficult to see any clear trends because the probability distributions are broad and the inferred SFHs are smooth.

To bring out subtler trends, we next quantify the probable fraction of GMP 3254’s stellar mass assembled by the time it fell into Coma. This is illustrated in the right panel of Fig. 3. The expected infall time is the time where its cumulative distribution (black curve) crosses 0.5, and the 68% confidence interval on this expectation (thick black bar) is bounded by the locations where the cumulative distribution crosses 0.16 and 0.84. The cumulative fraction of the stellar mass formed, according to the STECKMAP SFH, is shown with the red curve. At the expected infall time (∼6 Gyr), GMP 3254 had formed ∼90% of its final stellar mass, while at the boundaries of the 68% confidence interval on the infall time, it had formed ∼69% and ∼98% of its final stellar mass – we interpret this range as an approximate 68% confidence interval on the stellar mass assembled at the time of infall. We repeat the same measurement replacing the infall time probability distribution with the probability distribution for the time of the first pericentric passage to estimate the fraction of stellar mass assembled by that time, and its associated confidence interval, and repeat the same measurement for each galaxy in our sample. The resulting values are listed in Table 3.

The left panel of Fig. 4 shows the fraction of cumulative stellar mass formed at the expected infall (circles) and pericenter (stars) times for all galaxies in this study. The 68% confidence interval for the expected orbital times and the corresponding fractions of cumulative stellar mass are also shown. The points are colored according to the present-day stellar masses of the satellites. The satellites in our sample likely fell into the Coma cluster at lookback times of around 6–9 Gyr and had their first pericentric passages at lookback times of 1–5 Gyr, with the time between infall and pericenter ranging between 3–5 Gyr.

Fig. 4. The stellar mass formed in the Coma cluster environment as a function of orbital phase. Left panel: fraction of cumulative stellar mass for all the galaxies in this study at their expected infall (circles) and pericenter (stars) times. The galaxies are color-coded based on their present-day stellar masses. The error bars represent 68% confidence intervals on the expectation values. Higher stellar mass galaxies have formed a higher percentage of stellar mass at infall in comparison to the lower stellar mass galaxies. All galaxies have likely formed more than 95% of their stellar mass by the time of their first pericentric passage in Coma. Right panel: fractional increases in stellar mass between infall and pericenter time (black filled diamonds with 68% confidence intervals) are shown for all galaxies. The medians and confidence intervals are determined from the distribution of time differences and corresponding assembled mass differences computed for each individual matching N-body orbit (i.e., we do not plot the difference between the marginalized distributions shown in the left panel). The linear regression through the data (magenta line) hints that lower stellar mass galaxies may have formed relatively more stellar mass as compared to their higher stellar mass counterparts.

We find tentative evidence for a trend in the fractional mass formed after infall as a function of (present-day) galaxy mass, as illustrated in the right panel of Fig. 4. More massive galaxies seem to have formed a larger fraction of their final stellar mass relative to the lower mass galaxies at their expected infall times. The galaxies with higher stellar mass (log(M_⋆/M_⊙) ≳ 10) have typically formed more than ∼90% of their stellar mass by their expected infall time, while the lower mass galaxies have mostly formed around 80–90% of their stellar mass. However, GMP 3414 (log(M_⋆/M_⊙) = 10.73) and GMP 3664 (log(M_⋆/M_⊙) = 10.82), which are near the high-mass end of our sample, are outliers to this trend; they likely fell in ∼2 Gyr earlier than the other satellites. One plausible explanation for this could be that these high mass galaxies formed most of their stars early, so when these two fell in (early), they were at their peak of star formation and therefore somewhat more resilient to the cluster environment. The cluster environment was also perhaps a bit gentler earlier on, but this is probably not the main driver of this trend. We note that Wetzel et al. (2013) have also highlighted that star formation continues after infall.

All the sample galaxies likely formed more than ∼98% of their stellar mass by the time of their first pericentric passage and thus had nearly completely quenched by that time. Some galaxies have a nonzero probability of currently being on their first infall, which is translated into error bars for the pericenter time dropping below zero in some cases (see left panel of Fig. 4). Despite the wide confidence intervals, our measurements hint that these galaxies likely formed a significant fraction (≳5%) of their stellar mass while in the high-density environment of the Coma cluster.

4. Discussion

4.1. Quenching in the Coma cluster

Our results suggest that the galaxies we have studied in the Coma cluster are quenched near their first pericenter passage. Furthermore, the galaxies with higher stellar masses likely entered the cluster with a slightly larger fraction of stellar mass already formed; lower stellar mass galaxies formed roughly twice as much of their total stellar mass after infall. Using our main result shown in Fig. 4 (left panel) and comparing them to previous studies, we now discuss possible quenching processes along the orbits of the satellites.

To guide our interpretation, we first ask, based on our results (left panel of Fig. 4, Table 3), how long the galaxies in our sample could have remained star-forming. We quantify this by dividing the stellar mass yet to be formed at the first pericentric passage (based on the median estimates in Table 3, last column) by either the minimum “active” SFR corresponding to a specific star formation rate (SSFR) of 10⁻¹¹ yr⁻¹, or a representative SSFR for an actively star forming galaxy at z ∼ 0.25 (the typical pericenter time for our sample) of 3 × 10⁻¹⁰ yr⁻¹ (Madau & Dickinson 2014). The resulting values (Table 4) suggest that star formation typically ceases well under a gigayear after first pericenter in most cases, even in the contrived scenario where it occurs steadily at the threshold SFR for active star formation. We note, however, that the uncertainty in the fraction of stellar mass formed by pericenter for our galaxies is substantial, which prevents our ruling out much longer quenching timescales with confidence.

This serves to highlight that our measurements lead to a qualitatively different picture from most other low-redshift determinations of the quenching timescale. For instance, Wetzel et al. (2013), Rhee et al. (2020), Oman et al. (2021) suggest that star formation in satellites of comparable mass to those in our sample continues for several gigayears after the first pericentric approach, and that quenching is usually not complete until about the first apocentric passage.

We have not explicitly measured the quenching timescale in this study, so we have to be careful while comparing the quenching scenario obtained in our results with others. For example, Oman et al. (2021) is sensitive to quenching timescale for galaxies undergoing quenching now (i.e., those that had their first pericentric passage ∼t_q ago – they measure t_q with reference to this pericentric time). Rhee et al. (2020) give a z = 0 value assuming t_q ∼ (1 + z_inf)^−1.5. The method of Wetzel et al. (2013) models the full SSFR distribution, and so likely mixes galaxies being quenched at different times with a nontrivial weighting. However, we have estimated infall redshifts of z ∼ 0.6 for the galaxies in our sample, so we expect to be sensitive to t_q for satellites that fell in around this time. We therefore expect to find shorter quenching times (t_q is shorter at higher redshift in comparison to lower redshift – see Sect. 1) than determinations at z ∼ 0, but it is not clear by how much. For instance, the multiplicative scaling of t_q(z) assumed by Rhee et al. (2020) is ambiguous, as it is sensitive to the definition of the time zero-point. Their reference time is that of first infall across r = 2R_vir. Adopting instead the first pericenter reference time of Oman et al. (2021), for instance, would lead to a qualitatively similar, but quantitatively different prediction.

We therefore turn to higher-redshift determinations of the quenching timescale for guidance. We compile the quenching timescales (t_q) at different redshifts from several studies and show them in Fig. 5 as solid points. We note that the quenching timescales measured at higher redshift (1.0 ≲ z ≲ 1.5; Muzzin et al. 2014; Balogh et al. 2016; Foltz et al. 2018) are much shorter, by ∼3 Gyr, than lower-redshift (0 ≲ z ≲ 0.2) determinations (Wetzel et al. 2013; Taranu et al. 2014; Haines et al. 2015). Given that the galaxies in our sample fell in around z ∼ 0.6 (Fig. 5), we expect our results to be consistent with a timescale bounded between the higher and lower redshift measurements.

Fig. 5. Redshift dependence of the quenching timescale tq. The solid black points show the quenching timescales in clusters at different redshifts as measured by Wetzel et al. (2013, W13), Taranu et al. (2014, T14), Muzzin et al. (2014, M14), Haines et al. (2015, H15), Balogh et al. (2016, B16), Foltz et al. (2018, F18). The dashed and dotted lines show the free fall timescales tff from r = 2.5rvir and r = rvir, respectively, for the Coma cluster. The dash-dotted line shows the molecular gas depletion timescale (tdepl) with redshift (Tacconi et al. 2018), normalized to ∼1 Gyr at z = 0. The red and black ticks show the expected infall and pericenter redshifts, respectively, for the galaxies in this study.

We also show in Fig. 5 the redshift evolution of the free-fall timescale $t_{\mathrm{ff}}=0.5/\sqrt{G\overline{\rho }(r,z)}$ of the Coma cluster, where $\overline{\rho }(r,z)$ is the mean enclosed density within radius r at redshift z, from two different starting radii, r = r_vir (dotted line) and r = 2.5r_vir (dashed line⁵). Finally, we show with a dash-dotted line the evolution of the molecular gas depletion timescale with redshift,

$$\begin{array}{c}\hfill t_{\mathrm{depl}}\sim {(1+z)}^{-0.6}{({\delta }_{\mathrm{MS}})}^{-0.44},\end{array}$$

taken from Tacconi et al. (2018). They estimate the molecular gas mass by measuring CO line fluxes of a large sample of star-forming galaxies in the PHIBSS survey between 0 < z < 4. The depletion timescale depends on the offset in SFR from the star forming main sequence δ_MS = SFR/SFR_MS(M_⋆); we plot the curve corresponding to δ_MS = 0. This gives an indicative trend, and in any case the uncertain absolute calibration of our STECKMAP SFHs means that we can do little more than assume that our sample of galaxies had typical SFRs around the time they fell into Coma.

We note that t_depl and t_ff (r = r_vir) are similar at higher redshift, including at z ∼ 0.6, which corresponds to the expected infall times for our sample of satellites. The accretion of fresh gas for star formation in cluster satellites is cut off somewhere between r ∼ 2.5r_vir and r_vir, likely toward the latter, and from here the gas reservoir depletes without replenishment. For gas depletion to lead to quenching (SSFR < 10⁻¹¹ yr⁻¹) likely takes a few times t_depl as this timescale is that of an approximately exponential decline in the gas content⁶. At z ∼ 1, the quenching times t_q ∼ t_depl therefore seem somewhat too short to be consistent with a starvation scenario, suggesting that some gas removal by other mechanisms is likely required. This hypothesis is further reinforced by the fact that t_q ∼ t_ff, such that these galaxies are quenching around the time when gas removal by ram pressure or tidal stripping is most likely.

In contrast, at low redshift (z ∼ 0–0.2), the quenching timescale is significantly longer than both the free-fall time (from ∼r_vir), and the depletion timescale. This suggests that recently infalling satellites survive their first pericentric passage with a substantial gas reservoir, enough to sustain SF for several gigayears more, and that they quench after a few gas depletion timescales. This leads to a qualitatively different picture than at z ∼ 1: z ∼ 0 satellites seem to undergo starvation-driven quenching, with most of their gas supply being consumed by star formation rather than being stripped away.

As our sample of galaxies correspond to a time intermediate between the observations constraining the two scenarios discussed above, it is interesting to consider which one(s) they may be consistent with. Our finding that sustained star formation after first pericenter is unlikely for our sample of satellites is reminiscent of the z ∼ 1 “stripping-assisted” quenching outlined above. We therefore suggest that the transition to z ∼ 0 starvation-driven quenching occurs at z < 0.6. A key caveat to this argument is that if satellites in our sample had their supply of fresh gas cut off much earlier than around when they entered the Coma virial radius (e.g., if they were preprocessed in infalling groups), they may then plausibly quench by starvation by the time of their first pericentric passage in Coma. Bahé & McCarthy (2015) reported a similar prediction from simulations that predict a transition from stripping to consumption-driven quenching progressing from higher to lower redshift. Disentangling these possibilities will benefit from larger samples of satellites with resolved SFHs to help strengthen (or rule out) the tentative link we find between the quenching time and the first pericentric passage.

In the right panel of Fig. 4, we highlight that higher mass galaxies (log(M_⋆/M_⊙) > 10) have likely formed a lower fraction of their stellar mass (∼5%) between infall and pericenter than lower mass galaxies, which form almost double (∼10%) the fractional amount. Furthermore, all of the galaxies in our sample have nearly completed their star formation around the pericenter (see left panel of Fig. 4). This could simply be a reflection of the somewhat higher SSFRs of lower mass galaxies (see Fig. 5 of Ilbert et al. 2015, the SSFR evolution with redshift up to z ∼ 1.4 for galaxies in mass range of 9.5 ≤ log(M_⋆/M_⊙) ≤ 11.5). Alternately, this could be interpreted in terms of the quenching timescale, suggesting that massive satellites quench earlier than their lower mass counterparts. We note that this interpretation is at least qualitatively consistent⁷ with those of Wetzel et al. (2013), Rhee et al. (2020), who also suggest a shorter quenching timescale for massive satellites (see also Contini et al. 2020, who found similar results for cluster satellites at higher redshift).

The infall time probability distribution for GMP 3254, shown in the left panel of Fig. 3, is bimodal. This begs the question of whether our simple treatment of the distribution, characterizing it by its width, biases our conclusions. We have examined the case of GMP 3254 in more detail by decomposing the distribution into an earlier (older) and a later (younger) peak, and we have determined the percentage of stellar mass formed at the expected infall time of each peak. We find that the earlier peak is consistent with our overall qualitative conclusions drawn from the full distribution, while the later peak, considered alone, may hint that the galaxy was quenched – or nearly quenched – before it entered the cluster. However, extending this analysis to our entire sample is prone to a two caveats: first, we do not observe a consistent bimodal trend in PDFs of all the galaxies; second, the conclusions are not highly sensitive to our choice of ignoring the bimodality, and we have no clear reason to prefer one peak over the other. In addition, the cumulative distribution shown in the right panel of Fig. 3 explicitly accounts for the correlation across the two probability distributions.

Finally, we note that many Coma cluster galaxies in our study might have spent time in overdense regions prior to infall: at least some of them probably spent some time as satellites of infalling groups. This may be an interesting topic for future work as the approach we have adopted applied to larger samples may be able to find evidence for (or against) two distinct populations of satellites, corresponding to preprocessed and direct infall objects.

4.2. The stellar-to-halo mass relation

The Behroozi et al. (2010) SHMR (see Sect. 2.3) uses subhalo abundance matching (SHAM) to assign stellar masses to subhalos. The SHAM technique is very successful in matching observed galaxy statistics in spite of the differing evolutionary histories of satellite and central galaxies. There are several ways in which SHAM can be implemented; one particular way is based on the subhalo mass at the current epoch. However, since satellite subhalos are in general stripped of a fraction of their halo mass, before the stellar mass begins to be stripped, this introduces a bias for satellites. Wetzel et al. (2013, see their Appendix A) argue that M_max is better correlated with the stellar mass of subhalos since it always occurs before infall: it is not altered during satellite orbits. In our study, we have implemented SHAM based on M_max rather than the z = 0 M_sat.

The use of this relation relies on the assumptions that (i) the satellites have the same stellar mass now as they did at infall, in other words tidal stripping of stars and stellar mass growth through star formation are negligible, and (ii) the SHMR has not evolved significantly since the time when our sample of galaxies became satellites. Based on Table 3, the satellites have grown by up to about 25% (or 70% at 68% confidence) since infall. This corresponds to ∼0.1 dex (0.2 dex at 68% confidence). While this is several times smaller than the 0.5 dex interval in halo mass that we use to match simulated and observed satellites, the bias introduced is clearly systematic, causing all the halo masses to be overestimated.

The degree to which our sample of satellite galaxies may have been tidally stripped of stars is more difficult to assess; however, we think it unlikely that any satellites in our sample are “heavily” stripped of stars. This is first because the stellar portion of the galaxies are the most tightly bound, so that once the system begins to lose substantial amounts of stellar mass, it is already close to being completely disrupted (Bahé et al. 2019). This makes satellites with significant tidal stellar mass loss short-lived, so finding multiple examples in our sample is unlikely. Second, our satellites have a distribution of shapes typical of central ETGs with similar masses, with ellipticity (e = 1 − b/a) in the range 0.2–0.5 (Hoyos et al. 2011), while significantly stripped satellites are expected to be noticeably rounder, on average (e.g., Barber et al. 2014). We also note that any hypothetical stellar mass loss to tides would partially compensate the stellar mass growth reflected in the SFHs.

The SHMR is redshift dependent. For our sample of satellites, using the z = 1 SHMR instead of that at low redshift (z = 0.1) – the bulk of our sample have likely infall times around z ∼ 0.6 – would cause us to revise their halo masses upward by up to 0.2 dex (Behroozi et al. 2010). We note that this does not linearly combine with the systematic offset due to their stellar mass growth, described above since when using a higher redshift SHMR, the correct stellar mass to use is that “at the same redshift”. However, the biases introduced in our analysis by these two effects act in opposite senses, mitigating their overall severity.

We chose to use the Behroozi et al. (2010) SHMR, but there are many other published relations that could have been used in its stead. The Behroozi et al. (2010) relation is convenient in that it lies approximately in the middle of the scatter between the various proposed relations (e.g., the compilation in Behroozi et al. 2019, their Fig. 34). However, the SHMR does include a ∼0.15 dex scatter, which we have ignored; we find it unlikely that this additional scatter will make a significant change to our results. At the stellar mass scale of interest, the (minimum-to-maximum) scatter between the various proposed mean relations is about 0.3 dex. This represents a possible systematic bias affecting our analysis; unlike the others discussed in this subsection, the sign of this bias is unknown.

Based on the above discussion, we (conservatively) estimate that there may be up to a ∼0.5 dex systematic mismatch between the actual halo masses at infall of our sample of Coma satellites and those of the satellite halos selected from the simulations.

4.3. Systematic errors

We checked whether such a systematic bias could influence our qualitative interpretation of Fig. 4 by increasing all the halo mass estimates for our sample of satellites (Table 2) by 0.5 dex and repeating the rest of the analysis. The resulting shifts of the points in Fig. 4 are illustrated in Fig. 6 with dash-dotted black arrows. We repeated the same exercise but instead reduced all halo mass estimates by 0.5 dex; the resulting shifts are shown with black dashed arrows in Fig. 6. While many individual measurements move significantly, the overall qualitative picture is unchanged: these satellites form up to ∼15% of their stellar mass between infall and their first pericentric passage, and very little thereafter.

Fig. 6. Estimation of systematic errors. Left panel: change in fractional stellar mass formed at infall and pericenter times incurred by introducing a systematic 0.5 dex offset in the satellite halo mass (black arrows) or a 10% change to the Coma virial radius, propagated to all quantities that depend on this measurement (red arrows). Dash-dotted (dashed) arrows indicate the result of a positive (negative) shift in the relevant quantity. These shifts are typically small and hence do not affect the overall qualitative interpretation of our work. We note that the axes of the figure are zoomed-in relative to Fig. 4 to show the very small shifts. See Sect. 4.3 for further details. Right panel: change in fractional stellar mass formed between infall and pericenter due to the introduction of the same systematic offsets as in the left panel.

We also checked the sensitivity of our results for another possible bias, unrelated to those discussed in Sect. 4.2, by following a similar procedure. We adjusted the virial radius of the Coma cluster (R_vir, Coma) by ±10%, consistent with the plausible interval of 2.6−2.9 Mpc given in Łokas & Mamon (2003), and propagated the change to M_vir, Coma and σ_3D, Coma. We then repeated the selection of host systems from the N-body simulation, made new selections of satellite halos, and repeated the determination of the infall and pericenter times and the corresponding stellar mass fractions formed at those times. The relative shifts of the points in Fig. 4 are again shown in Fig. 6, with that due to an increased (decreased) R_vir, Coma shown with a dash-dotted (dashed) red arrow. We find again that, although some points may move substantially, the overall interpretation is robust against moderate uncertainty in the properties of the Coma cluster.

Finally, the fiducial result shown in Fig. 4 is based on STECKMAP SFHs using the MILES SSP models, as discussed in Sect. 2.1. STECKMAP provides four options for stellar population libraries: BC03 (Bruzual & Charlot 2003), MILES (Vazdekis et al. 2010), PHR (Le Borgne et al. 2004) and GD05 (González Delgado et al. 2005). Their general properties like spectral resolution and range of wavelength, age, and metallicity are listed in Table 5. We check the sensitivity of our results to the choice of SSP model by repeating our analysis using the two other SSP models available in STECKMAP: BC03 and PHR. Figure 7 shows the cumulative stellar mass formed at infall and pericenter, and the fraction of stellar mass formed between infall and pericenter, for all three SSP models. The 68% confidence intervals are somewhat wider when BC03 or PHR are used, but the results are otherwise very similar. There is a hint in Fig. 7 that using MILES leads to the weakest gradient in the stellar mass formed between infall and pericenter as a function of stellar mass, which mildly increases our confidence that this may be a real trend.

Fig. 7. Fraction of stellar mass formed at infall and pericenter time (top panel) and increase in fractional stellar mass between infall and pericenter (bottom panel) using different SSP models: MILES (left panel), BC03 (center panel), and PHR (right panel) are shown. MILES provides the most robust results, with the smallest error bars, in comparison to the other SSPs.

5. Conclusions

In this study, we have attempted to understand the environmental quenching of satellites in the high-density environment of galaxy clusters by linking their SFHs with their orbital histories. Our observed sample consists of high-S/N spectra of 11 ETGs (Trager et al. 2008) around the Coma cluster center, from which SFHs were inferred based on the relative SFR values obtained using STECKMAP (see Sect. 2.1). The intra-cluster orbits of these galaxies have been constrained based on their projected phase space coordinates and a library of orbits extracted from an N-body simulation (see Sects. 2.2 and 2.3). We match the cumulative distributions of orbital parameters against the satellite SFHs to infer the fraction of stellar mass formed by the expected infall and first pericenter times. Our main results (see Sect. 3.1, Table 3 and Fig. 4) are summarized as follows.

1. We find that infall across 2.5r_vir for galaxies in our sample occurred around z ∼ 0.6, and they reached their first pericenter ∼4 Gyr later.

2. Massive galaxies (log(M_⋆/M_⊙) ≳ 10) in our sample have formed a higher fraction of stellar mass (≳90%) before infall than less massive galaxies (∼80–90%), although the uncertainties are large.

3. Massive galaxies in the Coma cluster quenched on a shorter timescale than their lower mass counterparts.

4. Conversely, our sample galaxies formed a significant fraction of stellar mass (≳5%) within the high-density Coma cluster environment.

5. The galaxies in our sample are likely to have formed nearly all of their stellar mass (≳98%) by the time of their first pericentric passage. It is therefore likely that these galaxies quenched around (or within ≲1 Gyr after) first pericenter.

In contrast with studies sensitive to the current quenching timescale at z ∼ 0 that find that star formation in satellites continues well after their first pericentric passage (e.g., their Fig. 34 of Wetzel et al. 2013; Rhee et al. 2020; Oman et al. 2021), our measurements are consistent with truncation of star formation in present-day quiescent cluster galaxies around first pericenter. This is reminiscent of higher-redshift (z  ∼  1) determinations of the quenching timescale for satellites (e.g., Muzzin et al. 2014; Balogh et al. 2016; Foltz et al. 2018). Putting this together with our results, the typical infall redshift for our sample (z  ∼  0.6), the relevant free-fall and gas consumption timescales and their scalings with redshift, our interpretation is that galaxies in our sample likely lost substantial amounts of gas to ram pressure and/or tidal stripping. Put another way, their quenching was “not” primarily starvation-driven. We caution that this interpretation would change if most of the satellites in our sample were satellites of smaller groups before falling into the Coma cluster (preprocessed).

This exploratory analysis demonstrates the information-rich nature of SFHs when combined with orbital information from simulations. Repeating our study with the SFHs replaced with only a present-day star formation indicator (e.g., broadband color) would lead to only trivial conclusions. Our approach paves the way for future work leveraging large, deep spectroscopic surveys of clusters. WEAVE (Dalton et al. 2014), with its 2° field of view and nearly 1000 individual fibers, currently being installed on the William Herschel Telescope (WHT), will be an ideal instrument for such a study, and the Infall Regions subsurvey of the WEAVE Clusters survey will determine SFHs of tens of thousands of individual galaxies out to several virial radii in more than a dozen clusters.

Acknowledgments

We thank S. Bose and A. Jenkins for providing us the code and initial conditions for the N-body simulation and R. Peletier for a careful reading of an early version of the text. We also thank the anonymous referee for providing valuable (and rapid) feedback on the manuscript. KAO acknowledges support by the Netherlands Foundation for Scientific Research (NWO) through VICI grant 016.130.338 to M.A.W. Verheijen, and support by the European Research Council (ERC) through Advanced Investigator grant to C.S. Frenk, DMIDAS (GA 786910). This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.

References

All Tables

All Figures

Fig. 1. Star formation histories of galaxies in our sample determined using STECKMAP, using the MILES SSP model. The SFRs returned by STECKMAP are relative, not absolute, so we show them in units of Gyr−1. The curves are color-coded based on the satellite stellar masses. STECKMAP does not cover lookback times < 0.5 Gyr. The red and blue lines show the exponentially decaying form (SFR ∝ exp(τt)) of SFHs of red (early-type) and blue (late-type) galaxies respectively (Tojeiro et al. 2013).

In the text

	Fig. 2. Galaxies used in this study shown in projected phase space. The galaxies are colored by log(Msat/M⊙). The points are labeled with their GMP identifiers.
In the text

Fig. 3. Connecting orbital history with SFH. Left panel: probability distribution of infall time of GMP 3254 is shown in the gray histogram (left axis labels) and its SFH (right axis labels) is shown as the red points. The median value of its infall time probability distribution is shown as a dotted black vertical line. The SFR has dropped significantly by ⟨tinf⟩. Right panel: cumulative SFH – the fractional stellar mass accumulation over time – is shown (red line) along with the cumulative probability distribution of the infall time (black line) for GMP 3254. The median infall time (black filled circle) is marked with its 68% confidence interval (horizontal black line) at the 50th and 16th–84th percentiles, respectively. The corresponding accumulated stellar mass at the median infall time (red filled circle) and the its associated 68% confidence interval (vertical red line) are also shown. The arrows (black and red) indicate the translation of 16th and 84th percentile points on the cumulative distribution curve to boundary of the confidence interval lines. GMP 3254 (log10M⋆/M⊙ = 9.92) has formed around 90% of its stellar mass at the median infall time. Figure 4 (left panel) shows the accumulated stellar mass at the (median) infall and pericenter times for all galaxies in this study.

In the text

Fig. 4. The stellar mass formed in the Coma cluster environment as a function of orbital phase. Left panel: fraction of cumulative stellar mass for all the galaxies in this study at their expected infall (circles) and pericenter (stars) times. The galaxies are color-coded based on their present-day stellar masses. The error bars represent 68% confidence intervals on the expectation values. Higher stellar mass galaxies have formed a higher percentage of stellar mass at infall in comparison to the lower stellar mass galaxies. All galaxies have likely formed more than 95% of their stellar mass by the time of their first pericentric passage in Coma. Right panel: fractional increases in stellar mass between infall and pericenter time (black filled diamonds with 68% confidence intervals) are shown for all galaxies. The medians and confidence intervals are determined from the distribution of time differences and corresponding assembled mass differences computed for each individual matching N-body orbit (i.e., we do not plot the difference between the marginalized distributions shown in the left panel). The linear regression through the data (magenta line) hints that lower stellar mass galaxies may have formed relatively more stellar mass as compared to their higher stellar mass counterparts.

In the text

Fig. 5. Redshift dependence of the quenching timescale tq. The solid black points show the quenching timescales in clusters at different redshifts as measured by Wetzel et al. (2013, W13), Taranu et al. (2014, T14), Muzzin et al. (2014, M14), Haines et al. (2015, H15), Balogh et al. (2016, B16), Foltz et al. (2018, F18). The dashed and dotted lines show the free fall timescales tff from r = 2.5rvir and r = rvir, respectively, for the Coma cluster. The dash-dotted line shows the molecular gas depletion timescale (tdepl) with redshift (Tacconi et al. 2018), normalized to ∼1 Gyr at z = 0. The red and black ticks show the expected infall and pericenter redshifts, respectively, for the galaxies in this study.

In the text

Fig. 6. Estimation of systematic errors. Left panel: change in fractional stellar mass formed at infall and pericenter times incurred by introducing a systematic 0.5 dex offset in the satellite halo mass (black arrows) or a 10% change to the Coma virial radius, propagated to all quantities that depend on this measurement (red arrows). Dash-dotted (dashed) arrows indicate the result of a positive (negative) shift in the relevant quantity. These shifts are typically small and hence do not affect the overall qualitative interpretation of our work. We note that the axes of the figure are zoomed-in relative to Fig. 4 to show the very small shifts. See Sect. 4.3 for further details. Right panel: change in fractional stellar mass formed between infall and pericenter due to the introduction of the same systematic offsets as in the left panel.

In the text

	Fig. 7. Fraction of stellar mass formed at infall and pericenter time (top panel) and increase in fractional stellar mass between infall and pericenter (bottom panel) using different SSP models: MILES (left panel), BC03 (center panel), and PHR (right panel) are shown. MILES provides the most robust results, with the smallest error bars, in comparison to the other SSPs.
In the text