Constraining quenching timescales in galaxy clusters by forward-modelling stellar ages and quiescent fractions in projected phase space

We forward-model mass-weighted stellar ages (MWAs) and quiescent fractions in projected phase space (PPS), using data from the Sloan Digital Sky Survey, to jointly constrain an infall quenching model for galaxies in $\log(M_{\mathrm{vir}}/\mathrm{M}_{\odot})>14$ galaxy clusters at $z\sim 0$. We find the average deviation in MWA from the MWA-$M_\star$ relation depends on position in PPS, with a maximum difference between the inner cluster and infalling interloper galaxies of $\sim 1$ Gyr. Our model employs infall information from N-body simulations and stochastic star-formation histories from the UniverseMachine model. We find total quenching times of $t_\mathrm{Q}=3.7\pm 0.4$ Gyr and $t_\mathrm{Q}=4.0\pm 0.2$ Gyr after first pericentre, for $9<\log(M_{\star}/\mathrm{M}_{\odot})<10$ and $10<\log(M_{\star}/\mathrm{M}_{\odot})<10.5$ galaxies, respectively. By using MWAs, we break the degeneracy in time of quenching onset and timescale of star formation rate (SFR) decline. We find that time of quenching onset relative to pericentre is $t_{\mathrm{delay}}=3.5^{+0.6}_{-0.9}$ Gyr and $t_{\mathrm{delay}}=-0.3^{+0.8}_{-1.0}$ Gyr for our lower and higher stellar mass bins, respectively, and exponential SFR suppression timescales are $\tau_{\mathrm{env}}\leq 1.0$ Gyr and $\tau_{\mathrm{env}}\sim 2.3$ Gyr for our lower and higher stellar mass bins, respectively. Stochastic star formation histories remove the need for rapid infall quenching to maintain the bimodality in the SFR of cluster galaxies; the depth of the green valley prefers quenching onsets close to first pericentre and a longer quenching envelope, in slight tension with the MWA-driven results. Taken together these results suggest that quenching begins close to, or just after pericentre, but the timescale for quenching to be fully complete is much longer and therefore ram-pressure stripping is not complete on first pericentric passage.


INTRODUCTION
Galaxies have long been known to exhibit a bimodality in their star formation rates (SFRs) such that the population of bluer galaxies with higher SFRs make up the 'star-forming main sequence' and the redder galaxies with lower SFRs forming the quiescent 'red sequence' (e.g.Bell et al. 2004;Brinchmann et al. 2004;Brammer et al. 2011;Muzzin et al. 2012).The quiescent fraction increases strongly with increasing galaxy stellar mass, likely due to some internal ('mass') quenching process, as well as increasing with the density of its surrounding environment (e.g.Cooper et al. 2007;Peng et al. 2010), due to some environmental quenching process(es).Despite significant exploration of quenching trends and their statistical relation to internal and environmental properties, it remains a matter of debate which quenching mechanisms are responsible for the bulk of quenching in the densest environments across cosmic time, particularly galaxy clusters (Naab & Ostriker 2017;Wechsler & Tinker 2018).
★ E-mail: andrew.reeves@uwaterloo.caGalaxies falling into groups or clusters may have their star formation affected by a number of physical processes.Initially, upon infall into a cluster, a galaxy may experience enhanced star formation, due to ram pressure compressing and triggering collapse of a galaxies' cold gas reservoir (see e.g.Jaffé et al. 2018;Vulcani et al. 2018;Roberts et al. 2021).If strong enough, ram pressure can strip the cold gas in star-forming galaxies, triggering relatively abrupt quenching (normally < 1 Gyr is assumed to be indicative of this scenario, see e.g.Boselli et al. 2016;Fossati et al. 2018).Alternatively, numerous studies support stripping of just the hot halo gas could be leading to a 'strangulation' scenario, where a star-forming galaxy continues forming stars using its cold gas reservoir, but since its cold gas supply is not replenished, the galaxy then quenches once it runs out of cold gas (such as argued by Larson et al. 1980;Balogh et al. 2000;Bekki et al. 2002;Boselli & Gavazzi 2006;McGee et al. 2009;Taranu et al. 2014;Paccagnella et al. 2016).Such a scenario is often simply parameterized by a 'delayed-then-rapid' quenching prescription (Wetzel et al. 2013), with a range of time until quenching (after 'infall' -definitions for which vary) of 2 − 7 Gyr preferred by observations at  = 0 (see also e.g.De Lucia et al. 2012;Moster et al. 2018;Rhee et al. 2020;Oman et al. 2021, among numerous others).Other effects, like harassment from galaxies passing near each other, as well as mergers, can affect the star formation in a galaxy (Boselli & Gavazzi 2006).Complicating this picture is the 'pre-processing' of galaxies in host haloes of a lower mass groups prior to their infall into a cluster, which can additionally enhance the quiescent fraction in clusters compared to accreting isolated galaxies from the field.Although widely recognized as playing an important role in explaining the environmental dependence of star formation, work remains in characterizing the significance of pre-processing for quenching, where exactly it occurs, and how it evolves with redshift (Fujita 2004;McGee et al. 2009;De Lucia et al. 2012;Webb et al. 2020).It is at least clear that quenching is significantly enhanced relative to the field at least up to  ∼ 1.5 in rich galaxy groups (Gerke et al. 2007;Wetzel et al. 2012;Mok et al. 2013;Reeves et al. 2021) and clusters (Van Burg et al. 2018;Pintos-Castro et al. 2019;Van Burg et al. 2020).
Constraining timescales using forward modelling of observable quantities in projected phase space (PPS) in clusters can provide significant information for discerning when and where different physical quenching mechanisms may be dominant.In observational works, galaxy properties have shown a strong link to their position in PPS (for a variety of distinct examples see e.g.Mahajan et al. 2011;Smith et al. 2012;Muzzin et al. 2014;Jaffé et al. 2015;Noble et al. 2016;Gavazzi et al. 2018;Jaffé et al. 2018;Kelkar et al. 2019;Kim et al. 2020;Rhee et al. 2020;Roberts et al. 2021).Little detailed forward modelling work has been done on constraining quenching timescales in PPS, aside from that developed in the series of papers by Oman, Hudson and collaborators (Oman et al. 2013;Oman & Hudson 2016;Oman et al. 2021).The second and third of these papers use the quiescent fraction of cluster galaxies in PPS coordinates as well as corresponding distribution of infall times inferred from Nbody simulations.The PPS coordinates used in this context, namely the radial distance from cluster center and line-of-sight velocity offset from its host system, only carry magnitude (i.e.no sign) information in practice, as it's generally not possible to distinguish between a satellite receding towards a background host or a background satellite receding away from a foreground host.This method provides a nearby interloper population on the cluster outskirts (outside some 3D radial cut) against which to compare cluster satellites.This has an advantage over using a generic field population: using these interlopers in the cluster outskirts provides an already pre-processed infalling galaxy sample, allowing better isolation of cluster physics from the physics involved in pre-processing.
Whereas many previous works that modeled environmental quenching have relied on quiescent fractions, observational indicators that are sensitive to the star formation history over longer timescales should provide more information.Luminosity-or massweighted ages (MWAs) in particular not only provide complementary quantitative information about the star formation history of galaxies, but are also insensitive to whether a given galaxy's star formation history is assumed to be smooth or stochastic (see discussion in Section 4.2.2).
Studies using spectroscopic indicators of stellar age have long established that more massive galaxies, or, more accurately, galaxies with deeper potential wells, have older ages (Nelan et al. 2005;Graves et al. 2007).Dependence on environment is weaker: early studies focused on the age difference between field and cluster ellipticals after controlling for mass or velocity dispersion (Thomas et al. 2005;Bernardi et al. 2006).Radial trends of the stellar ages of galaxies within clusters (at fixed velocity dispersion) were found (Smith et al. 2006) with a stronger environmental dependence on the ages of dwarf galaxies (Smith et al. 2012).
Studies of stellar age in PPS are less common.Pasquali et al. (2019) considered luminosity-weighted stellar ages in zones of PPS around low redshift clusters, finding ages that are older by ∼ 2 Gyr for the innermost lower-velocity regions of PPS compared to the outermost regions.Recently, Khullar et al. (2022) studied D4000-derived ages in PPS in 0.3 <  < 1.1 clusters, again finding older ages for PPS regions associated with the core of the cluster (see also Kim et al. 2022).
There has been even less effort modelling the effect of quenching on the stellar ages.Taranu et al. (2014) used infalling subhalo orbits from -body simulations to forward-model the age-sensitive Balmer line Lick indices as a function of cluster-centric radius.Their models preferred long quenching timescales occurring for galaxies that have passed within ∼ 0.5 200 (i.e.preferring slow 'strangulation' over rapid ram-pressure stripping).
As little work has been done forward modelling spectroscopicallyderived ages in clusters and none has been carried out in full PPS with a large sample (but see Upadhyay et al. 2021, for a first attempt with a sample of 11 galaxies in the Coma cluster), the goal of this work is to extend the modelling of Oman & Hudson (2016) and Oman et al. (2021), who studied the quiescent fraction,  Q , in PPS.Oman et al. (2021) found that using  Q alone is not sufficient to constrain both the time of quenching onset and the timescale over which quenching occurs on average and so they marginalized over this latter parameter to find the average time at which 50 per cent of infalling star-forming galaxies have quenched.Extending this and similar work could make use of, for example, the entire sSFR distribution, but as we will describe, the shape of the quiescent bump is sensitive to observational systematics, while the depth of the green valley is sensitive to choice of star formation history (smooth versus stochastic).In this paper, we will focus on extending Oman et al. (2021) by forward modelling the observed MWAs of galaxies in PPS, in addition to using quiescent fractions.This additional observable will enable joint constraints on both the time of onset of infall quenching and the timescale for the duration of quenching preferred by both observables, for the population of galaxies as a whole.
Our paper is structured as follows.In Section 2, we describe how we parameterise PPS and the SDSS-based datasets, from which we get our  Q and mean MWAs.We then, in Section 3, introduce our forward modelling approach (including sources of star formation histories and infall histories) and present our quiescent fraction and mean MWA modelling results in PPS as well as the constraints on quenching timescales they provide.In Section 4, we discuss the robustness of our results, how they compare to literature, and what they imply for the dominant infall quenching processes.We then conclude in Section 5.
Unless otherwise specified, the following assumptions and conventions are used.Uncertainties are estimated from the 16 th -84 th percentile interval (equivalent to 1- for a Gaussian-distributed variable).Logarithms with base 10 (log 10 ) are written simply as 'log' throughout this work.A flat ΛCDM cosmology consistent with the Planck 2015 cosmological parameters (Planck Collaboration et al. 2016) is assumed, namely  0 = 68 km s −1 Mpc −1 , Ω  = 0.31, and Ω Λ = 0.69.We define a virial overdensity at  = 0 in this work as ∼ 360 times the background density, Ω    , corresponding to the density of a recently virialized spherical top-hat solution (see e.g. Bryan & Norman 1998)1 .A Chabrier (2003) initial mass function (IMF) is assumed throughout.

PPS selection and conventions
We adopt the convention of Oman et al. (2013) for selecting objects in PPS to allow direct comparison to observations: there are two sets of cluster-centred coordinates: (, ) for the 2 × 3D physical phase space coordinates and (, ) for the normalized PPS coordinates.Specifically, the PPS coordinates consist of the distance to cluster centre perpendicular to the line-of-sight and the velocity component along the line of sight, which we choose to be the Cartesian -axis of a given simulation dataset.The projected radius between the cluster centre and a galaxy is then and the projected velocity relative to the motion of the cluster centre of mass is Using an absolute value for the projected velocity definition captures an assumption that the distances to clusters and satellites are not measured with enough accuracy to determine the sign of their relative velocity.The additional  term corrects for the Hubble flow.To allow easy comparison between clusters, for both simulations and observations, we have normalized the PPS coordinates.For ease of comparison to Oman & Hudson (2016) and Oman et al. (2021), radial coordinates are normalized by dividing by the 3D virial radius of the cluster,  vir,3D , and velocity coordinates by the 3D velocity dispersion,  3D , of a given cluster.As we are averaging over many clusters in both the orbit library and observed sample, we can assume that clusters are approximately spherically symmetric.This in turn allows relation of the 3D velocity dispersion to the observable 1D velocity dispersion by  3D ≈  1D / √ 3. We note that in order to ensure consistency between the normalization factors used for  and , they should not be independent, i.e. the velocity dispersion should be derived from the virial radius or vice-versa.See Section 2.4 for details of how these are calculated for the SDSS data.

Infall histories: N-body simulation orbit library
We use the PPS orbital libraries from Oman et al. (2021); they base their methodology on Oman & Hudson (2016) to produce a catalogue of orbit parameter probability distributions from an N-body simulation.As described in Oman et al. (2021), the 'level 0' voidsin-void-in-voids simulation (VVV;Wang et al. 2020) is used, which has a box size of 500 ℎ −1 Mpc on a side, mass resolution elements of 10 9 ℎ −1 M ⊙ , and a force softening scale of 4.6 ℎ −1 kpc.They run the simulation to scale factor  = 2, corresponding to redshift  = −0.5 or ∼ 10 Gyr into the future (i.e.negative lookback times), allowing them to find the time of first pericentre for > 99.9 per cent of resolved satellite galaxies.Host halo masses and satellite halo masses  The fraction of interlopers in (R,V) projected phase space bins determined using the orbital library from Oman et al. (2021). is distance from cluster centre, in units of  vir,3D and  is the velocity offset from cluster centre, with units of  3D .The red diagonal line,  = − (4/3)  +2, guides the eye as to where ∼ 50 per cent of galaxies are interlopers, aside from a portion around (,  ) ∼ (1.2, 0.4), which has many galaxies at or close to their first apocentre after entering the cluster.The dotted line delineates a region in the upper right that is heavily dominated by the interloper population, which we use to compare between SDSS and UniverseMachine interlopers.
are estimated using the rockstar halo finder; infall times and times of first pericentre are derived from halo merger trees generated by the Consistent Trees utility.These merger trees include all haloes with > 30 particles, corresponding to  vir ≳ 4 × 10 10 M ⊙ , allowing resolution of  vir ∼ 2.5 × 10 11 M ⊙ haloes, which host the lowest stellar mass galaxies in our observed sample - ★ ∼ 10 9.5 M ⊙ (at least until they are stripped of ≳ 85 per cent of their halo mass).
Satellites are identified as haloes within a 3D aperture of 2.5  vir at  = 0 for log( vir /M ⊙ ) > 12 host systems (although we only use hosts with log( vir /M ⊙ ) > 14 in this work).Satellite primary progenitors/descendants are traced backward/forward through time and orbits relative to the  = 0 primary progenitor/descendent of their host are recorded.The time of first pericentre is not interpolated (it is very challenging to do this properly and is unnecessary for our purposes), so it is worth noting here that this results in non-uniform output times, with a median timestep of 220 Myr and a maximum timestep of 380 Myr, sufficient to resolve the characteristic quenching and stripping timescales fit in Oman et al. (2021).Worth noting is that the satellite halo mass,  sat , is its maximum mass at  ≥ 0, since maximum halo mass is better correlated with stellar mass in 'moderately stripped' satellites; see e.g.Conroy et al. (2006) and Appendix A of Wetzel et al. (2013) for further elaboration on this.
For the orbit library itself, since satellites are all galaxies within 2.5 vir , interlopers are then naturally defined as all galaxies that fall within 2.5 vir in projection, but outside 2.5 vir in 3D.Only a vanishingly small fraction, ≪ 1 per cent, of interlopers would have been classified as satellites at an earlier time.All satellites have a recorded time of first infall into the final cluster (2.5 vir in 3D) and time of first pericentre ( fp ) in the final cluster; we only use  fp in our modelling and analysis.We make use of their relative abundance compared to selected satellites to define the probability that an 'observed' galaxy is an interloper.For the sake of illustration and intuition for the reader, we plot the statistical fraction of galaxies that are interlopers at a given position in PPS, rather than cluster satellites, in Fig. 1.In the plot it's clear that at both high  and  the vast majority of galaxies are interlopers and vice-versa for lower  and .We note that around (, ) = (1.5 − 2, 0 − 0.5) there are many galaxies on their first infall into a cluster, whereas galaxies at around (, ) = (0 − 0.5, 1 − 2) are primarily galaxies at their first pericentre.Over several orbits galaxies settle to low (, ).As well, we plot a straight line indicating where approximately 50 per cent of galaxies are interlopers -this should provide a helpful rule-of-thumb for the reader in the plots of various quantities in PPS that will follow.

Star formation histories: the UniverseMachine
For individual galaxies' SFR histories we use the UniverseMachine semi-analytic model (Behroozi et al. 2019) which parametrises galaxy SFRs as a function of halo potential well depth, redshift, and assembly history.Their halo potentials and assembly history are derived from the Bolshoi-Planck dark matter simulation (Klypin et al. 2016;Rodríguez-Puebla et al. 2016), which features a periodic comoving volume 250ℎ −1 Mpc on a side with 2048 3 (∼ 8 × 10 9 ) particles.The simulation employs a flat ΛCDM cosmology compatible with Planck15 results (Planck Collaboration et al. 2016), with stored snapshots equally spaced in log() (180 intervals).The rockstar halo finder (Behroozi et al. 2013a) and Consistent Trees (Behroozi et al. 2013b) are used to construct merger trees, the same codes used in the construction of the merger trees for the orbit library employed in this work (Section 2.2).
For our purposes, it was important that their model reliably reproduce the observed SFRs and quiescent fractions over time, especially in overdense regions in SDSS (although not necessarily cluster cores, as we are only using UniverseMachine interloper galaxies, described in the following paragraph).A key aspect of the Uni-verseMachine model is that SFRs are stochastic -they fluctuate over time -and are linked to the merger history of the halo.Behroozi et al. (2019) found that without such scatter the quiescent fractions of satellites would be too high (relative to centrals).Specifically, higher SFRs are assigned to halos with higher levels of halo growth, parameterised by the relative logarithmic growth in  max , over the past dynamical time (where  max is the maximum circular velocity of the halo).The stochasticity has two components: a short timescale component (∼ 10 − 100 Myr) representing internal processes affecting local galactic cold gas, and the second component linked to the dynamical time (∼ 1 Gyr).This stochasticity will be a key difference in our comparison to other infall models (which use smooth star formation histories).A galaxy can transition between quiescent and star-forming, as long as the overall quiescent fraction and other properties are matched to the various observations they used, in a self-consistent manner.
We construct a PPS catalogue using the UniverseMachine simulation results at  = 0, assuming arbitrarily that the -axis is the observational line-of-sight.For each of the 144 UniverseMachine simulation sub-volumes, galaxies were selected from around central haloes ('UPID'==-1) with halo masses log( vir /M ⊙ ) > 14.All galaxies (including centrals) are selected so that their projected radius is  < 2.5 and projected velocity is  < 2, where the PPS coordinates (, ) are defined as described previously in Section 2.1.True interlopers and satellites, as opposed to interlopers and satellites as defined in PPS, are identified with 3D cuts of  > 2.5 and  < 2.5, respectively.

SDSS dataset
For the observed sample of galaxies, we build on the sample used by Oman et al. (2021) for their quenching analysis.As in Oman et al. (2021), we use the SDSS Data Release 7 catalogue (Abazajian et al. 2009), supplemented with SFRs (Brinchmann et al. 2004;Salim et al. 2007) and stellar masses (Mendel et al. 2014), which were used in Oman et al. (2021) for estimating galaxy subhalo masses using subhalo abundance matching.We also use the spectroscopic sample of galaxies, which we note suffers from fiber collisions, although Oman et al. (2021) determine through a detailed exploration that the impact of this should be minimal.
We use the mass-weighted ages, as well as the accompanying stellar mass estimates, of Comparat et al. (2017), who performed full spectral fitting of galaxy properties using FIREFLY (Wilkinson et al. 2017).We compare their stellar masses with those of Mendel et al. (2014) and confirm that offsets in the stellar masses do not affect our results.In particular, we use the fits done using the stellar population models of Maraston & Strömbäck (2011), a Chabrier (2003) IMF, and the MILES stellar library (Sánchez-Blázquez et al. 2006;Falcón-Barroso et al. 2011;Beifiori et al. 2011) for the mass-weighted ages ('CHABRIER_MILES_age_massW'), and stellar masses ('CHABRIER_MILES_stellar_mass').
For simplicity, we focus on high-mass cluster haloes, log( vir /M ⊙ ) > 14, where environmental effects should be most extreme.As in Oman et al. (2021), we select candidate satellite galaxies around groups and clusters from the Lim et al. (2017) andVon Der Linden et al. (2007) catalogues, with cluster virial masses peaking at ∼ 3 × 10 13 M ⊙ and ∼ 3 × 10 14 M ⊙ , respectively.We exclude groups with redshifts  < 0.01, as their members are generally too bright to be covered by SDSS spectroscopy.We note that the two group catalogues are constructed with different algorithms.Von Der Linden et al. ( 2007) look for overdensities of galaxies with similar colours, while Lim et al. (2017) use a friends-of-friends group finding algorithm.This does not significantly affect our analysis as we only use the cluster masses, centres, and line of sight velocities from the catalogues, rather than membership information.Satellite velocity offsets are normalized using cluster velocity dispersions, calculated from the virial masses of the clusters in the catalogues.The velocity dispersions are similar to the dark matter particle velocity dispersion of the systems (within ≲ 10 per cent; see Munari et al. 2013, table 1); this is important since the velocity dispersion of dark matter particles in the host halo is used to normalize the velocity offsets of satellite haloes in the simulations that the Oman et al. (2021) orbit library is derived from (see Section 2.2).Host halo masses are calculated following eq. 1 of Von Der Linden et al. (2007) and eq. 4 of Lim et al. (2017).These are converted to virial masses assuming an NFW density profile (Navarro et al. 1997) and mean mass-concentration relation from Ludlow et al. (2014), with differences between assumed cosmologies accounted for.The mean enclosed density is then used to define the virial radii, where Δ vir () is the virial overdensity in terms of the mean matter density Ω  ()  crit (), with critical density  crit = 3 2 /(8).
Group velocity dispersions are defined following Biviano et al. (2006), but with redshift dependence from Bryan & Norman (1998); explicitly this is We use two stellar mass bins, 9 < log( ★ /M ⊙ ) < 10 and 10 < log( ★ /M ⊙ ) < 10.5.We note that for log( ★ /M ⊙ ) > 10.5 there are too few star forming galaxies to yield robust results, so we do not consider this mass range in this paper.
Constraining satellite quenching with stellar ages 5

Quiescent fraction trends in PPS around SDSS clusters
Our first source of information for deriving quenching timescales is the quiescent fraction in SDSS clusters.As we are using the same observational data sample as Oman & Hudson (2016) and Oman et al. (2021), we define define star-forming and passive galaxies using the same cut, namely log(sSFR/yr −1 ) = −0.4log( ★ /M ⊙ ) − 6.6, ( with star-forming (quiescent) galaxies lying above (below) the specific star formation rate (sSFR) relation.
We illustrate the quiescent fraction trends in PPS for 9 < log( ★ /M ⊙ ) < 10 galaxies in Fig. 2. In the figure we see wellknown observed trends of decreasing  Q with increasing distance from cluster centre as well as higher  Q for galaxies with lower velocity offsets.We draw a line indicating roughly where ∼ 50 per cent of galaxies are interlopers, with galaxies above this line increasingly dominated by interlopers, moving from lower to higher R and V coordinates.Enhanced quenching largely falls within the region bounded by this 50 per cent interloper line in PPS.Contrasting with the field results of Oman et al. (2021, see upper right panel of e.g.fig. 10 to see star-forming fraction versus stellar mass for the SDSS field), with field defined as the average overall SDSS population, we see that the interloper quiescent fraction is higher than the field by 0.05 and 0.09 for our lower and higher stellar mass bins, respectively.
Having information on how quenching depends on both the radial position and velocity relative to host cluster centre gives us very useful information on when and where quenching occurs.It's clear in Fig. 2 that most quenching of infalling star forming galaxies occurs well within the virial radius.Modestly enhanced quenching up to high velocities indicates quenching is at least beginning to occur within the first infall or just after first pericentre.Although we have some useful information from  Q in PPS, this information alone results in significant degeneracies between time of quenching onset and duration of quenching, as we will clearly show later in this work (see Section 3.2).To break this timescale information degeneracy we must turn to another independent observable as a function of its position in PPS.

Mass-weighted age -stellar mass relation
For the stellar ages of SDSS galaxies we use a simple proxy to capture when the bulk of the stellar mass in a galaxy was formed, namely the mass-weighted age (MWA).We include only quiescent galaxies for our MWA-based measurements, since the emission lines of starforming galaxies fill in the age-sensitive Balmer absorption lines, making them rather unreliable as an age indicator.In this subsection, we note a trend in MWA with stellar mass in both the observationallyderived SDSS values and the UniverseMachine simulated galaxies.As we are only interested in the differential change in MWA relative to an infalling population, we first describe how we control for the MWA- ★ effect in this work.After this, we will briefly illustrate observed MWA trends in PPS in Section 2.4.3.
We show the trend in MWA with stellar mass for both the observationally-derived SDSS values and the UniverseMachine (Behroozi et al. 2019) simulated galaxies side-by-side in Fig. 3.The mean and median trends that we see for SDSS and Uni-verseMachine are qualitatively similar, but the UniverseMachine galaxies are offset in overall MWA such that they are younger than SDSS MWAs by ∼ 3 Gyr on average.We note that mass-weighted ages derived from spectra depend on the assumed star formation history.For example, Allanson et al. (2009) showed that a range of star formation histories can fit the observed stellar absorption lines equally well in quiescent cluster galaxies, yielding mass-weighted ages that range from 5 Gyr (single burst) to 8 Gyr (exponentially declining star formation history) to 12.8 Gyr ('frosting' model) for the lowest-mass cluster galaxies.Since we use ages only in a differential sense in our analysis, i.e. comparing the difference in MWA in the cluster core to the infalling regions in PPS, any systematic offset in the ages due to different star-formation history assumptions will not affect our analysis.For SDSS, the MWA- ★ relation is approximately flat, with MWA ∼ 9 Gyr for 9 < log( ★ /M ⊙ ) ≲ 10.5 galaxies, which then increases for log( ★ /M ⊙ ) ≳ 10.5 to ∼ 12-13 Gyr for the highest stellar mass galaxies, log( ★ /M ⊙ ) > 11.The increase in MWA for UniverseMachine galaxies from lowest to highest stellar masses shown in Fig. 3 is similar in magnitude to that for SDSS (∼ 4 Gyr), albeit with a more steadily increasing monotonic trend, i.e. a non-flat trend in MWA below log( ★ /M ⊙ ) < 10.5.The scatter is also quite noticeably different -the SDSS scatter can be thought of as the natural distribution of galaxies (shown modeled here by the UniverseMachine ) convolved with the scatter induced in estimated MWAs from uncertainties in SED-fitting the galaxy spectra.
To remove the systematic trend in MWA with stellar mass and the ∼ 3 Gyr offset between the SDSS observationally-derived MWAs and the simulated UniverseMachine MWAs so that we can focus solely on the effect of environment on galaxy MWAs, we look at the mean MWA residual from the mean MWA- ★ relation.This allows us to measure the differential effect between inner cluster regions and cluster outskirts (i.e.controlled for stellar mass).In particular, we use interloper galaxies as our reference infalling (already pre-processed) population.
For a simple clean proxy interloper sample that can be easily compared between UniverseMachine and SDSS, we choose galaxies in the region of PPS given by as illustrated in Fig. 1 (see also Section 3.1.1).Based on statistics from the orbit library of Oman et al. (2021), this region of PPS provides a sample that is made up of ≳ 97 per cent true interlopers.Correcting for the ∼ 3 per cent impurity is unnecessary for our purposes.
For SDSS, we parameterise the mean MWA- ★ relation (fit and residual shown as black curve in the left panels of Fig. 3) for (predominantly interloper) galaxies with  < 2.5 and  < 2 in the stellar mass range 10 9 <  ★ /M ⊙ < 10 12 using the function where  = 0.1 Gyr is the normalization,  = 9.16 Gyr is the vertical offset,  crit = 10 10.35 M ⊙ is the location of the break in slopes, and  = 0 and  = 40 specify the lower and higher stellar mass slopes, respectively.The fit is chosen to minimize the residual between the running mean MWA- ★ relation and the curve, and an offset is added so that the interloper population has mean(MWA − MWA SDSS ) = 0 by construction.This was done rather than fit the interloper MWA- ★ relation directly since the interloper population had relatively few galaxies to give a good fit across the whole stellar mass range.
For the UniverseMachine simulated sample, a polynomial is fit to the ∼pure interloper population (Eq.6) to the running mean MWA (fit and residual shown as the black curve in the right panels of Fig. 3).The polynomial fit has the following form: where  ≡ log( ★ /M ⊙ ) and the coefficients have units of Gyr.The quiescent fraction is significantly higher in cluster centres than in their outskirts (≳ 30 per cent).For the lower stellar mass bin, regions within the ∼ 50 per cent interloper (red diagonal) line, show at least modestly elevated quiescent fractions, with little evidence of enhanced quenching immediately outside of (above) this line.For the higher stellar mass bin, some higher velocity galaxies within the virial radius also show somewhat elevated quiescent fractions.

Trends in the mean deviation of mass-weighted age in PPS
We present a binned map of the deviation in mass-weighted age from the mean MWA- ★ relation of the SDSS interloper proxy sample, mean(MWA − MWA SDSS ) (see previous section), as a function of the PPS coordinates  and  in Fig. 4. Errors on the mean values are calculated by bootstrapping over clusters.The figure shows that mean(MWA − MWA SDSS ) is enhanced for galaxies at lower radii and velocities.This is particularly clear for the higher stellar mass bin.Within the 50 per cent interloper line shown in the figure, we calculate mean(MWA − MWA SDSS ) = 0.76 +0.12 −0.14 and mean(MWA − MWA SDSS ) = 0.66 +0.08 −0.09 for the lower and higher stellar mass bins, respectively.Our higher stellar mass bin is consistent with Kim et al. (2022), who find that populations in PPS around clusters at 0.3 <  < 1.4 identified as 'early infall' were older on average by 0.71±0.4Gyr than comparable 'recent infall' galaxies.
In the following section we will show that the distribution of modeled mean(MWA − MWA UM ), ΔMWA, in PPS differs from that of  Q (shown previously in Fig. 2, Section 2.4.1).These differences in observed   and ΔMWA distributions will provide us with a distinct set of preferred values for time of quenching onset and the duration of quenching timescales.

MODELLING AND RESULTS
Our goal is to model the effects of infall quenching on two SDSS observables,  Q and mean ΔMWA, in PPS around  ∼ 0 galaxy clusters, and then use these results to constrain the associated quenching timescales.In this section, we first describe our infall quenching model in detail in Section 3.1.We show our model results for  Q in Section 3.2 and for ΔMWA in Section 3.3.Finally, we combine our modelling results for both  Q and ΔMWA to find a joint constraint on our model's two quenching timescales in Section 3.4.

The model
To model the effects of infall quenching on both  Q and mean ΔMWA, we use the UniverseMachine semi-analytic model, as described in Section 2.3, as a representative model of star formation histories for galaxies in the infall region around a given cluster.A simple, smooth parametric model could be used, but such models don't have any intrinsic scatter in their star formation histories.We will see that this scatter has an impact on the bimodality in the star formation rates of the galaxy population (see Section 4.2).The UniverseMachine simulation also provides us with a pre-processed infall sample.
Using a galaxy cluster PPS orbit library from Oman et al. (2021), as described in Section 2.2, we sample the infall time distribution of galaxies in PPS in galaxy clusters, including the interloper fraction across PPS.With these tools in hand, we describe in Section 3.1.1how we construct an infall quenching model to predict  Q and ΔMWAs as a function of position in PPS.For the sake of intuition, we visually illustrate some basic trends in  Q and ΔMWA with our quenching parameters in Section 3.1.2.

Description of our satellite quenching model
We take the previously defined interloper sample as the infalling population, which will naturally include any pre-processing of galaxies in infalling structures (e.g.groups or filaments), as captured by the UniverseMachine model.Our model assumes that star formation histories of  = 0 cluster satellite galaxies are the same as the  = 0 infalling population of galaxies, but with the addition of a quenching mechanism on infall, as described below.Galaxies that fall into the cluster begin suppression of their SFR when a time interval of  delay has elapsed after their first pericentre passage.Defining the model this way means that we are modeling the mean quenching timescales for the population of galaxies as a whole, over all possible infall orbits (e.g.circular versus eccentric/plunging), giving us robust mean observed properties, insensitive to the choice between stochastic or smooth star formation histories (see Section 4.2.2 for discussion of this point).
For UniverseMachine galaxies, we use the same sSFR cut to separate quiescent/star forming galaxies as Behroozi et al. (2019), namely log(sSFR/yr −1 ) = −11, where star-forming (quiescent) galaxies are defined as lying above (below) the cut.The UniverseMachine interloper population and SDSS interloper population differ slightly in their quiescent fractions (see e.g.Fig. A1), but we are only interested in the differential quenching effect (  Q in the cluster versus some infall population).
To find the timescales associated with infall quenching, we will fit for the infalling  Q .In Appendix A we will confirm that our results are robust to adjusting the UniverseMachine infalling star formation histories such that the infalling galaxies closely match the SDSS In the literature, some variation of a 'delayed-then-rapid' quenching scenario is common, like that proposed and explored in Wetzel et al. (2013).Exponentially declining SFR (e.g.Wetzel et al. 2013;Roberts et al. 2019) or linearly declining quiescent fraction (such as in Oman et al. 2021) are common assumptions.To facilitate com-parison with previous work, we choose a model with an exponential suppression of the SFR on a timescale,  env , that starts after some time delay,  delay , relative to the time of first pericentre  fp (we note that this allows for negative values of  delay , i.e. quenching onset prior to  fp ).
Galaxies with a given stellar mass and time of first pericentre are Monte Carlo sampled from satellites in the orbital library of Oman et al. (2021).This gives, for each part of PPS, some mix of interloper galaxies, as well as in-cluster galaxies, which have a distribution of infall times depending on position in PPS.The star formation history for a random interloper with the closest stellar mass to that just drawn is then also drawn from UniverseMachine , which gives the star formation history of this Monte Carlo-simulated galaxy.Each of these simulated satellite galaxies have quenching applied by starting the decline in the SFR at the time  fp +  delay , i.e.SFR sat = SFR  (),   where () is the multiplicative quenching envelope, where  refers to a given Monte Carlo-sampled interloper (infalling galaxy) star formation history.In Fig. 5, we illustrate the effect of the multiplicative envelope on the stochastic star formation histories of the UniverseMachine galaxies.
We then define our model of the quiescent fraction in a given PPS bin as where Δ   is the increase in quenching from our model and  Q,infall is left as a free parameter of the infalling population's quiescent fraction to account for systematic differences between UniverseMachine and the SDSS data.
To determine the stellar masses of any galaxy at  = 0, the Uni-verseMachine SFRs are integrated numerically according to eq.21 and eq.22 in Behroozi et al. (2019).Explicitly, this gives a final stellar mass of where  loss () = 0.05 ln 1 + /1.4Myr and  now is the age of the Universe at  = 0. Integrated true (simulated) stellar masses are offset to give 'observed' stellar masses, according to the prescription in Behroozi et al. (2019, eq. 25 and eq. 27), namely by offsetting them by  = SM obs −SM true =  0 +  (1−) [dex] (where  = 1/(1+) is the cosmological scale factor) and additionally adding Gaussian scatter of  SM,obs = min( SM,z (1 + ), 0.3) [dex].Best-fitting values associated with UniverseMachine catalogues available online 2 were used for these stellar mass adjustments, namely  0 = 5.6 × 10 −3 and 2 UniverseMachine data release 1 catalogues using the Bolshoi-Planck   = −0.03used for the offset and  SM,z = 0.069 for the Gaussian scatter parameter, respectively.
Mass-weighted ages are integrated numerically from the stellar mass calculations, using the midpoints between the simulation snapshot timesteps and d/d between two snapshots' timesteps, where  true is the total simulated stellar mass at the present time without any corrections for observational effects.Note that MWA is expressed as a lookback time.
The MWA predicted by our quenching model is then Similar to   in Eq. 11 above, we define the ΔMWA observable of our model in a given PPS bin for each galaxy as where MWA  ( delay ,  env ) is the MWA from our quenching model, MWA infall is a floating offset to allow for the possibility that the infalling population's MWA may be different than that of the interloper phase space bins in equation 6, and MWA UM is the MWA expected from the mean MWA- ★ relation as defined in Eq. 8 for UniverseMachine galaxies.We use the mean ΔMWA in a given PPS bin in the rest of this work.

Simple illustration of our quenching model
Now that we have described our infall quenching model, we illustrate our predictions for trends in  Q and ΔMWA in PPS.In particular, we start by exploring the effect of varying the quenching time delay relative to a satellite galaxy's first pericentre passage,  delay , assuming instantaneous quenching (i.e. env = 0).We first illustrate this visually in Fig. 6 for  delay ranging from  delay = −2 Gyr (i.e.prior to first pericentre) to  delay = 4 Gyr.We include contaminating interlopers.
For  delay = −2 Gyr, quiescent fractions are very high, ≳ 0.75 throughout most of the region dominated by satellites (i.e.below the ∼ 50 per cent interloper line, indicated by the same diagonal line as in Figs. 2 and 4).ΔMWA reaches ∼ 1.1 Gyr for  < 0.5 in the cluster core ( < 0.5). Q and ΔMWA values increase with decreasing radius and also generally increase with decreasing velocity offset (albeit less strongly than with radius).
As  delay increases, there is a decrease in both  Q and ΔMWA, as intuitively expected. Q is still clearly enhanced relative to the interloper population within the ∼ 50 per cent interloper line, even for a  delay as long as 4 Gyr.ΔMWA reaches ≲ 0.25 Gyr at all (, ) bins for  delay = 4 Gyr.For an even longer  delay ≳ 6 Gyr (not shown), continues to drop and ΔMWA drops to ∼ 0 Gyr.With this physical intuition in mind, we now carry out a quantitative analysis of our model for  Q and ΔMWA in turn, as well as the joint constraint that observed SDSS   and ΔMWA values place on our model.We will see in the subsections that follow, differences in the trends for these two observables will enable us to break degeneracy in the  delay and  env model parameters.

Quiescent fraction: comparison of models and data
We plot  Q model predictions as a function of time of quenching onset relative to the first pericentre,  delay , for a range of  env values (different coloured curves) in Fig. 7 for 9 < log( ★ /M ⊙ ) < 10 galaxies.For the equivalent plot for 10 < log( ★ /M ⊙ ) < 10.5 galaxies at  ∼ 0, see Fig. B1 in Appendix B. For reference, we show the mean  Q (grey horizontal line) and the bootstrapped error on the mean (shaded gray region) for the SDSS sample in each PPS bin.When finding the best-fitting  delay values (see Section 3.4 for this), we also fit for the infalling population's  Q as a nuisance parameter, to account for possible systematic differences between the quenching from UniverseMachine and the actual infall popu-lation, that we will then marginalize over.For the infalling galaxies, our fits prefer  Q,infall = 0.16 ± 0.01 (  Q,infall = 0.48 ± 0.01) for 9 < log( ★ /M ⊙ ) < 10 (10 < log( ★ /M ⊙ ) < 10.5) galaxies, shown in Table 1, which are offset −0.06 ± 0.01 (0.09 ± 0.01) relative to the infall/interloper population's MWAs predicted by the UniverseMachine model.
We now turn our attention to the trends in  Q with  delay and  env , by focussing on the region of PPS where the infall quenching effect in our model is most pronounced, namely the innermost parts of clusters.In e.g. the  < 0.5,  < 0.5 bin we see no effect on  Q for very long quenching times ( delay ∼ 10 Gyr), as expected. Q then steadily increases with decreasing time delay up until a point at which it begins decreasing again, at e.g. delay ∼ −3 Gyr.This turnover feature comes about because the star formation histories of infalling galaxies are being truncated so aggressively that many quenched galaxies are dropping out of our stellar mass range, with increasing numbers of quenched galaxies having log( ★ /M ⊙ ) < 9. Since there are fewer high stellar mass galaxies, truncated higher stellar mass galaxies dropping into the bin are not enough to compensate for those dropping below log( ★ /M ⊙ ) = 9.The effect is most pronounced for instantaneous quenching,  env = 0, and less pronounced for longer star formation suppression timescales,  env > 0;  env acts to smooth out dependencies of  Q on  delay .This delay in quenching due to longer  env values also results in an offset between the maximum possible  Q as a function of  delay .This offset effect is such that for higher values of  env the maximum possible  Q (the turnover point) requires much earlier quenching, significantly before the time of first pericentre.Features on shorter timescales (i.e.bumps and wiggles) shown in our model curves are due to the discrete snapshot times of the N-body simulation used for the orbit library in Oman et al. (2021).
We contrast this with the methodology of Oman et al. (2021), where stellar masses were not truncated by quenching, since for simplicity they did not make use of star formation histories.In that case, there is no turnover effect in  Q as a function of  delay since galaxies are not dropping out of the stellar mass range.That said, the  Q data do not prefer excessively negative values of  delay and so the shift in preferred delay time compared to our choice of truncating stellar masses is rather minimal: it simply causes an offset of ∼ 0.5 Gyr (earlier) relative to our model.In PPS, the models reproduce the general trends in  Q visible in the SDSS data, but there is clearly significant degeneracy in the preferred model (i.e. which  delay ,  env combination is preferred).In order to break this degeneracy, we need another observable beyond the quiescent fraction to constrain  delay and  env .As mentioned earlier, we will use observationally-derived stellar ages of galaxies for this purpose.We explore model predictions for stellar ages of galaxies in PPS in the next section.

Stellar ages: comparison of models and data
We now compare our models of stellar age to the SDSS data.As in the previous section, in Fig. 8, we plot the mean ΔMWAs as a function of  delay in bins of PPS coordinates  and  for 9 < log( ★ /M ⊙ ) < 10.For the equivalent plot for 10 < log( ★ /M ⊙ ) < 10.5 galaxies, see Fig. B1 in Appendix B. We note that for the infalling galaxies, we find that our fits prefer offsets of MWA infall = 0.59 ± 0.12 Gyr and MWA infall = 0.47 ± 0.07 Gyr for 9 < log( ★ /M ⊙ ) < 10 and 10 < log( ★ /M ⊙ ) < 10.5 galaxies (shown in Table 1), which correspond to infalling MWAs of 6.70 +0.11  −0.13 Gyr and 7.63 ± 0.07 Gyr.As in the previous section, looking at e.g. the innermost, lowvelocity bin ( < 0.5,  < 0.5), we see no effect on ΔMWA for very long quenching delay times, with the effect on ΔMWA usually increasing with decreasing  delay (earlier onset of quenching).This general trend, for the most part, is similar for different values of  env , but with shallower slopes for higher values of  env .Or, put another way, a longer timescale for the suppression of star formation (a higher value of  env ) would require an earlier onset of quenching to have the same degree of impact on the stellar age.Interestingly, the trend in higher ΔMWA for earlier onset of quenching relative to time of first pericentre is not monotonic for all values of  env : for example, for  env ∼ 0 (instantaneous or almost-instantaneous suppression of star formation) around  fp that there is a small ∼ 0.1 − 0.3 Gyr dip.This dip is most pronounced for the higher velocity bins, as this is where galaxies reaching their first pericentre will preferentially be located in PPS -close to the cluster core and moving at a high velocity.Since the figure is only showing the mean ΔMWA of quiescent galaxies, and this region of PPS is dominated by galaxies on their first infall, the influx of recently quenched galaxies lowers the mean ΔMWA relative to values of  delay slightly earlier/later than  delay ∼ −1 Gyr (the deepest part of the dip).Longer suppression timescales blur out this effect.Since the lower velocity regions in the cluster core ( < 0.5) are dominated by galaxies that entered the cluster a long time ago, this effect is washed out.Although the uncertainties on ΔMWA are larger than for  Q , our model still predicts clear trends for ΔMWA across PPS.Because the ΔMWA trends with the timescale parameters differ from those obtained from  Q , we explore the constraints that arise from combining  Q and ΔMWA in the next section.

Fitting and joint constraints of time delay and exponential quenching timescale
By running our model across a range of  delay and  env values we can jointly constrain these timescale parameters by using the information contained in the trends of  Q and ΔMWA as a function of the PPS coordinates.For each  env we perform  2 fitting to determine the obs, and   are the mean observed SDSS data value and its uncertainty, respectively, for either   or MWA, in a given PPS bin . model, is the modeled value, as defined in Eq. 11 and Eq. 16, respectively, and is a function of  delay ,  env , and the infall population.
The sum is over all PPS bins.We find that the uncertainties on the data may be underestimated, since for the 20 − 3 = 17 degrees of freedom (20 PPS bins, three free parameters) we find the following reduced  2 values:  2 red (  Q ) ≈ 2.75 and  2 red (ΔMWA) ≈ 2.47 for 9 < log( ★ /M ⊙ ) < 10.Similarly, for 10 < log( ★ /M ⊙ ) < 10.5 we find  2 red (  Q ) ≈ 1.38 and  2 red (ΔMWA) ≈ 1.25.A high  2 leads to parameter constraints that may be too tight, so to be conservative we inflate the uncertainties on the observed SDSS mean  Q and mean ΔMWA in each PPS bin by a factor of √︃  2 red .We note that uncertainties on the models are negligible as they can be made arbitrarily small simply by sampling more UniverseMachine galaxies.
We present our best-fitting values for  delay as a function of  env for our observables  Q and ΔMWA in Fig. 9, for both 9 < log( ★ /M ⊙ ) < 10 (left panel) and 10 < log( ★ /M ⊙ ) < 10.5 (right panel) stellar mass bins.Consistent with the PPS results presented in Section 3.3, ΔMWA provides a constraint on  delay that depends only weakly on  env .For  Q there is a steeper trend between the fit  delay for a given  env .These relations make sense intuitively: as  Q depends on the distributions of SFRs at a given moment in time, it can be greatly impacted by any change in the timescale of SFR suppression,  env .ΔMWA, on the other hand, will not change significantly with  env since the bulk of stellar mass growth occurred in the past and will mainly depend on the time at which a galaxy has quenched.
Also shown in Fig. 9 are contours showing the confidence intervals for  delay and  env jointly when combining  Q and ΔMWA.We list our best-fitting joint values with uncertainties (16 th and 84 th percentiles) defined by the  delay and  env 's marginal probability distributions in Table 1.
Looking first at the 9 < log( ★ /M ⊙ ) < 10 bin in Fig. 9 (left panel), quenching time delay relative to time of first pericentre is preferred to be  delay = 3.5 +0.6 −0.9 Gyr.We find a relatively rapid suppression of star formation once quenching has started, with a best-fitting  env ≤ 1.0 Gyr (95% confidence level used as an upper bound).We choose to use a one-sided upper limit for  env for the lower stellar mass bin, as the probability distribution for  env is not Gaussian, peaking at  env ∼ 0 and the 50 th -percentile occurring at  env = 0.3 Gyr.
Turning to the higher stellar mass bin, 10 < log( ★ /M ⊙ ) < 10.5 (bottom panel), a different fit is preferred: an earlier onset of quenching,  delay = −0.3+0.8  −1.0 Gyr, corresponding to the onset of quenching beginning close to the time of first pericentre.The corresponding best-fitting  env is 2.3 +0.5 −0.4 Gyr, indicating a significantly slower suppression of star formation.

DISCUSSION
We now discuss the meaning of our modelling results and how they contrast with previous literature.For a discussion of the robustness of our results to changes in the overall UniverseMachine star formation histories, see Appendix A. Mean ΔMWA predictions for quiescent 9 < log(  ★ /M ⊙ ) < 10 galaxies for a range of models where galaxies quench after some delay time, relative to time of first pericentre.We show models run for a range of exponential suppression timescales.The SDSS mean ΔMWA values are shown as grey lines, with the shaded regions showing the bootstrapped (over clusters) uncertainty in the mean.Note that our models include an additional offset to allow for the possibility that the infalling population may be different than that of the interloper phase space bins in Eq. 6.
We first compare with the work that we have directly built upon, namely that of Oman et al. (2021) (see also : Oman et al. 2013;Oman & Hudson 2016).We then compare our results to the alternative framework proposed in the seminal work of Wetzel et al. (2013), followed by other observational studies: Taranu et al. (2014) which made use of stellar age-related spectral indices and Rhee et al. (2020) which used PPS information.Finally, we compare to a study examining the predictions of hydrodynamical simulations (Wright et al. 2022;Lotz et al. 2019).
For the purposes of our discussion, we define a quantity,  Q , as the average time for a star-forming galaxy with median sSFR (of a star-forming galaxy) to exponentially decline in SFR until it crosses the UniverseMachine quenching threshold of 10 −11 yr −1 .For our lower (higher) stellar mass bin this difference in sSFR is ∼ 1.0 dex (∼ 0.8 dex).Assuming a constant stellar mass for simplicity, the quenching timescale is  Q =  delay + 2.33 env ( Q =  delay + 1.82 env ).These timescales are summarized in Table 2 and their contours are also plotted in Fig. 10.

Oman et al. (2021)
The satellite quenching model of Oman et al. (2021), based on modelling quiescent fractions in PPS, employs a maximum likelihood model constraining four parameters: the quiescent fraction of an infalling population (  before ), the final quiescent fraction of galaxies after satellite quenching has taken place (  after , set to zero for their core analysis), the time at which half of the drop in quiescent fraction is complete ( mid , relative to time of first pericentre), and the timescale to go from the initial to final quiescent fraction (Δ).For the purpose of consistency across the following discussion sections, we will refer to the quenching time as  Q rather than  mid .
We've previously discussed how the preferred time of quenching in Oman et al. (2021) is offset to be ≲ 0.5 Gyr later than ours, since for simplicity they assumed galaxy star formation histories are not truncated by premature quenching.
They find a nearly flat  Q relation with stellar mass (excluding log( ★ /M ⊙ ) ≳ 10.5), with  Q = 3.7±0.2Gyr and  Q = 3.4±1.0Gyr for our lower and higher stellar mass bins, respectively.We find consistent overall quenching times ( Q ( delay = 0) = 3.7 ± 0.4 Gyr and  Q ( delay = 0) = 3.9 ± 0.2 Gyr, respectively).Relaxing this assumption, our model is still consistent with theirs, although the lower stellar mass bin  Q increases slightly, by ∼ 0.5 Gyr. 9 < log( ★ /M ⊙ ) < 10 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 env [Gyr] Marginal best-fitting parameter values on our model's quenching timescales for the two observables,  Q (blue), mean ΔMWA for quiescent galaxies (orange), as well as joint (black lines/grey shading) for 9 < log(  ★ /M ⊙ ) < 10 galaxies (left) and 10 < log(  ★ /M ⊙ ) < 10.5 (right).Uncertainties on the best-fitting parameters are shown with corresponding shaded regions (darker and lighter shaded regions corresponding to the 68% and 95% confidence regions, respectively).The parameters fit here are the quenching delay time relative to time of first pericentre,  delay , and the exponential suppression of star formation timescale,  env .Dashed lines indicate the best-fitting  delay at a given .The joint 68% and 95% confidence region from both observables for  delay and  env are shown overlaid with black contours.The joint best fitting model, i.e. the location with the peak joint probability, is indicated with a cross symbol.For each plot, the other stellar mass bin's joint constraint is overlaid (black cross and dotted contours).Marginal best-fitting values and uncertainties from the joint constraints are tabulated in Table 1. 9 < log( ★ /M ⊙ ) < 10 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 env [Gyr]  Quenching time, t Q [Gyr] this work Wetzel+2013 Taranu+2014 Rhee+2020 Oman+2021 Wright+2022 Figure 10.Marginal and joint best-fitting parameter values on overall quenching time,  Q , and  env for the two observables,  Q (blue) and quiescent ΔMWA (orange) for 9 < log(  ★ /M ⊙ ) < 10 galaxies (left) and 10 < log(  ★ /M ⊙ ) < 10.5 (right), analogous to Fig. 9.As well, we have overlaid the  Q literature values from Table 2 with the various points and lines shown in the legend.We omit error bars where it was not straightforward to determine them.Oman et al. (2021) were unable to provide strong constraints on their Δ parameter; they find that it tends to be underestimated, e.g. with the median value (based on the probability density from the model) underestimated by up to ∼ 1 Gyr at lower stellar masses.They note that this effect is due to resolution issues in the N-body simulation resulting in low-mass satellite haloes being disrupted too early.Values of Δ ∼ 1 − 1.5 Gyr were preferred by their modelling, but with very large uncertainties, which they marginalized over for their primary analysis to get their tight constraints on  Q and the quiescent fraction parameters.They note that the resolution effect could mask a decreasing trend in  Q , which could likewise occur in our modelling.Wetzel et al. (2013) found that satellites falling into a cluster experience 'delayed-then-rapid' quenching.One difference between their model and ours is that although Wetzel et al. (2013) fits a delayedthen-rapid model with an exponentially declining sSFR (specifically, their parameters are: time delay,  Q,delay , and 'fading timescale', Table 2. Comparison of timescales from various literature results to our work; in each table section they are listed in the order in which they appear in Section 4. 'Lower' and 'higher' stellar mass bins refer to 9 < log(  ★ /M ⊙ ) < 10 and 10 < log(  ★ /M ⊙ ) < 10.5, respectively.Upper section: average time  Q until quenching, relative to time of first pericentre.Lower section: exponential timescales  compared to our exponential suppression best-fitting timescale as well as with the conditional assumption ( * ) that quenching begins at time of first pericentre ( delay ≡ 0), using the  Q contour only.† indicates use of the upper 95% confidence level as an upper bound.Note that definitions of  in the literature vary; the respective discussion sections should be consulted when making comparisons between  values.

Study
Quantity Lower  ★ bin Higher  ★ bin This work (best)  2013) finds a unimodal distribution for long suppression times (contrast the red solid line for their model with the blue solid line for our model).This arises because of our use of a stochastic star formation history, whereas they use smooth star formation histories.In particular, we are contrasting with fig. 9 in Wetzel et al. (2013), for our 10 < log(  ★ /M ⊙ ) < 10.5 stellar mass bin and their 10.1 < log(  ★ /M ⊙ ) < 10.5 stellar mass bin.A correction of −1 Gyr is added to be comparable to their average infall time, which has the time zero-point set at  vir .Additionally, for ease of comparison with Wetzel et al. (2013), who claimed need for a short exponential suppression timescale to reproduce the bimodality, we include only galaxies within  200c ≈ 0.73 360m and add log-normal scatter with mean log(sSFR/yr −1 ) = −12 and 0.25 dex variance to all galaxies with log(sSFR/yr −1 ) < −12.
fade ), our model relies on the star formation histories of Uni-verseMachine which are stochastic rather than a simple smooth analytical function.
Another important difference in the models is the treatment of preprocessing of galaxies prior to entering the cluster.In our model, a simulated infalling population is used to account for pre-processing, whereas in Wetzel et al. (2013) the infall quenching 'clock' for a given galaxy starts on first infall into any larger halo, i.e. when a central becomes a satellite.Because of this, one would expect their (average) quenching time delay to be longer, as their quenching time delay includes the time spent in smaller groups prior to infalling into the main cluster progenitor.The magnitude of the median difference between time of first infall into any halo and infall into the final cluster is 3.2 Gyr (their fig.2; assuming log  200c/M ⊙ ( = 0) = 14.2 -our median host halo mass).

Quenching delay time
To compare their model with ours, we must account for the two populations of galaxies falling into clusters which will have different quenching times relative to time of first pericentre in the  ∼ 0 cluster: those that are already satellites in some smaller halo which have had their quenching clock start at some earlier time (A) and those that are falling into a larger halo for the first time (B).For case (A), according to the fits of Wetzel et al. (2013) for a 9 < log( ★ /M ⊙ ) < 10 galaxy, the full time until quenching for a galaxy that became a satellite in a log( vir /M ⊙ ) > 14 cluster is ∼ 5.7 Gyr (interpolated).Taking the median time between first infall and infall into the cluster to be ∼ 3.2 Gyr (see fig. 2 of Wetzel et al. 2013) and time from  vir to first pericentre in the final log( vir /M ⊙ ) > 14 cluster to be ∼ 1 Gyr (Oman et al. 2013), we would expect such galaxies to quench t ∼ 1.5 Gyr after first pericentre on average.We note that these pre-processed galaxies are quenching earlier than they would in a scenario where quenching occurs only due to infall into the final cluster host.For case (B), according to their fits for galaxies falling into log( vir /M ⊙ ) > 14 clusters, the full time until quenching is  Q ∼ 4.8 Gyr for 9 < log( ★ /M ⊙ ) < 10 galaxies, or about ∼ 3.8 Gyr after first pericentre.The relative proportions of the (A) and (B) scenarios are ∼ 0.6 and ∼ 0.4, respectively (see Wetzel et al. 2013, fig. 1).
These approximations give an average quenching time, relative to time of first pericentre, of  Q ∼ 0.6(1.5 Gyr) + 0.4(3.8Gyr) ∼ 2.4 Gyr.The same calculation for our higher stellar mass bin gives  Q ∼ 1.5 Gyr.These are shorter, by ∼ 1.9 Gyr and ∼ 2.4 Gyr, than our quenching times of  Q = 4.3 ± 0.4 Gyr and  Q = 3.9 ± 0.2 Gyr, for our lower and higher stellar mass bins, respectively.Wetzel et al. (2013) claim that in order for the bimodality in the sSFR distribution to persist across environments, including in the most massive clusters, infall-induced quenching requires a short exponential timescale,  Q,fade < 1 Gyr (they find  Q,fade = 0.8±0.2Gyr and  Q,fade = 0.5 ± 0.2 Gyr for our lower and higher stellar mass bins, respectively).For longer SFR fading timescales (∼ 2 times their best-fitting  fade ), their model predicts no bimodality in the sSFR distribution, but rather a unimodal distribution peaking near the centre of the 'green valley' (log(sSFR) ∼ −11.2 for log( ★ /M ⊙ ) < 10.5 galaxies), which we illustrate in Fig. 11.This occurs because their model assumes a smooth exponentially-declining star formation history for the infalling galaxies, and a gentle fading leads to galaxies  .5 (right).The left panel shows three scenarios, which are (in the order that they are listed in the legend): (i) long  equal to that preferred for our higher stellar mass bin, but  delay taken from along the  Q error contour, (ii) a scenario where  delay = 0, with  env chosen to minimize  2 (see Fig. 9), (iii) 50 th percentile of our timescale parameters.The right panel is similar, showing: (i) rightward end of the 95% confidence level ellipse, (ii) best-fitting timescale parameters, (iii) leftward end of the 95% confidence level ellipse.It is worth noting that there is an offset in  Q not factored into the above model histograms, equal to that described in Section 3.2, and that the quenching sSFR cut for SDSS has some stellar mass dependence, rather than a simple sSFR = 10 −11 yr −1 cut.

Fading time and sSFR bimodality
passing through the green valley too slowly to maintain the bimodality.We find our lower stellar mass bin exponential timescale ( env ≤ 1.0 Gyr) is consistent with theirs, but our higher stellar mass bin value of  env = 2.3 +0.5  −0.4 Gyr is significantly longer than their  Q,fade ∼ 0.5 ± 0.2 Gyr.However, these fading times are different quantities:  Q,fade in Wetzel et al. (2013) represents the fading time of an individual galaxy, whereas our  env is a population-wide fading 'envelope' which multiplies the stochastic star formation histories of individual UniverseMachine galaxies.The stochasticity in star formation histories is due to significant variations in SFR of Uni-verseMachine galaxies, which occur by design: in this model, SFR tracks the accretion rate of baryonic matter, as well having a random variability on short timescales (∼ 10 − 100 Myr).This stochasticity allows the sSFR distribution to be bimodal even in the case of a long  env .Because of the rapid changes in SFR for UniverseMachine galaxies, they do not spend significant amounts of time in the green valley, whether or not their SFRs are forced to decline by our exponential suppression envelope (visually illustrated in Fig. 5).The SFRs of infalling galaxies undergoing satellite quenching drop below the sSFR quenched threshold much more abruptly than a suppressed smooth star formation history, so longer fading/suppression envelopes are allowed.
From a qualitative comparison of the scenarios in Fig. 12, a longer  env ≳ 2.3 Gyr is preferred by the depth of the green valley for 9 < log( ★ /M ⊙ ) < 10 galaxies, rather than a short  env ≤ 1.0 Gyr preferred by the  Q and ΔMWA constraints.We note that the location of the star forming peak of the sSFR distribution is better fit by e.g. delay = 0 Gyr and  env ∼ 1.5 Gyr.For 10 < log( ★ /M ⊙ ) < 10.5 galaxies, we find that the relative depth of the green valley implies a preference for longer  env , e.g. env ≳ 2.3 Gyr, which is our bestfitting result.
A robust quantitative comparison of infall quenching models using additional sSFR distribution features, aside from   , clearly requires controlling for type of star formation history.The depth of the green valley depends on whether the assumed star formation histories are smooth or stochastic and the specific choice of stochasticity made in UniverseMachine is not unique.Additionally, measured star formation rates of quiescent galaxies have high systematic uncertainties regardless of how they are obtained (Brinchmann et al. 2004;Salim et al. 2007;Wetzel et al. 2012;et al. 2014), meaning that much of the shape of the quiescent bump is determined by observational systematics, rather than quenching physics.As such, we do not utilize additional features of the sSFR distribution beyond   , as they do not appear informative or robust in constraining our infall quenching model and comparing it to e.g.Wetzel et al. (2013).

Taranu et al. (2014)
Using an approach broadly similar to our work, Taranu et al. (2014) combined a star formation and quenching model with information from a library of subhalo orbits from N-body simulations of 4 rich clusters with log( 200 /M ⊙ ) > 15, with median stellar masses in their sample of log( ★ /M ⊙ ) = 10.4.They compared model predictions to observed bulge and disc colours and are one of the only previous studies (see also Upadhyay et al. 2021, for an attempt with a spectroscopic sample of 11 galaxies) to make use of age-sensitive stellar absorption line-strength indices for constraining cluster infall quenching times using radial information.Since galaxy disc colours are particularly affected by environment, they can be used to constrain the quenching timescale and location of quenching onset in the cluster using their median trends in cluster-centric radius.They used the age-sensitive Balmer lines of luminous passive cluster galaxies, which are more sensitive to older stellar populations than disc colours.
From their set of models, they found that quenching starting shortly before the time of first pericentre (within a quenching radius of 0.5 200 , which is ∼ 0.5 Gyr prior to the time of first pericentre for a galaxy on first infall) produced the best fits to disc colours as a function of radius.They found that delayed-then-rapid models lead to an excessively large slope in galaxy disc colour versus cluster-centric radius, and also overpredict the strength of the Balmer lines.Instead, their preferred model is one with a gentle exponential quenching timescale of  post ∼ 3 Gyr, consistent with the  env that we find for the higher stellar mass bin.This corresponds to a total quenching time of  Q ∼ 5 Gyr, longer than ours by ∼ 0.7 Gyr (see Table 2).

Rhee et al. (2020)
Rhee et al. ( 2020) fit a quenching model to disc galaxies in  = 0.08 clusters.Infall is defined as when a galaxy first crosses 1.5 vir .They follow an approach similar to ours, parametrizing quenching with a time delay followed by an exponential suppression of star formation, but instead of modelling the quiescent fraction they model the full SFR distribution of galaxies as a function of position in PPS.
They find that the time delay from infall until the onset of quenching is 2 Gyr for all stellar masses log( ★ /M ⊙ ) > 9.5.The time for a galaxy to fall from their infall definition of 1.5 vir to first pericentre is ∼ 1.3 Gyr.We can then subtract this from their quenching time since infall (as presented as  Q in their table 1 for similar stellar mass bins), giving 5.45 Gyr − 1.3 Gyr ∼ 4.2 Gyr and ∼ 3.8 Gyr(interpolated) −1.3 Gyr = 2.5 Gyr for our lower and higher stellar mass bins, respectively.This is consistent with our lower stellar mass bin (we find  Q ∼ 4.3 ± 0.4 Gyr) but earlier for our higher stellar mass bin by ∼ 1.4 Gyr (our  Q ∼ 3.9 ± 0.2 Gyr).
This difference in trend with stellar mass between our results and those of Rhee et al. (2020) is not due to a difference in delay times, but rather is due to their SFR suppression timescale,  cluster ( = 0), which declines with increasing stellar mass over our stellar mass range, whereas we find the opposite trend.We note that they add a redshift-dependent factor to their quenching timescale  cluster ( inf ) =  cluster,0 (1 + ) −  and to the quenching time delay timescale (due to the redshift evolution of the dynamical time in clusters),  d ( inf ) =  d,0 (1 +  inf ) −1.5 .This only results in a small correction (a decrease by ∼ 10 per cent if we use e.g. = 0.5, the time at which an average galaxy has fallen into the cluster).For  = 0, they find  cluster,0 = 1.7 +0.2 −0.3 Gyr and  cluster,0 ∼ 1.1 +0.3 −0.2 Gyr for our lower and higher stellar mass bins, respectively.This is consistent with our lower stellar mass bin's  env ( delay = 0) = 1.6 ± 0.2 Gyr, but more rapid than our higher stellar mass bin's  env ( delay = 0) = 2.2 ± 0.1 Gyr, respectively.

Comparison with hydrodynamical simulations
A study directly comparable to ours, which examines quenching in the EAGLE hydrodynamical simulations, is that of Wright et al. (2022).Their work provides an orbital analysis of galaxies' gas inflow, stripping, star formation and quenching.Stripping of infalling galaxies' hot gas begins at 2-3 virial radii from the host and takes longer for high-mass satellites (log( ★ /M ⊙ ) > 10).This begins the process of starvation and the removal of the hot gas 'buffer', resulting in galaxies becoming vulnerable to cold gas stripping.They refer to the hot gas halo as having a protective effect, as they observe that the onset of significant cold gas stripping only begins after stripping of the hot gas halo is complete.In their work, they include both H I and molecular hydrogen in the mass of cold gas in a galaxy.
For low-mass satellites (log( ★ /M ⊙ ) < 10), they find that suppression of gas cooling onto the galaxy becomes permanent after hot gas stripping, normally occurring around the time of first pericentre.Some high-mass satellites, on the other hand, retain small hot gas reservoirs, and continue to cool gas for star-formation after first pericentre.Cold gas stripping is shown to be periodic, being strongest for galaxies near pericentre, as that is when density and velocity is at a maximum, hence maximizing the ram-pressure force ( ram ∝  icm  2 ).All of this results in the following: low mass satellites experience very efficient ram-pressure stripping of cold gas, leading to rapid quenching, whereas high mass satellites experience less efficient stripping and a more gradual starvation-like scenario after their first pericentre.
In terms of quenching, we note that Wright et al. ( 2022) remove pre-processed infalling satellites, namely those galaxies which were satellites in a host halo with log( 200 /M ⊙ ) ≥ 12 prior to falling into the current host (≈ 30 per cent of infalling galaxies).This handling of pre-processing is somewhat different from ours, as we fit the preprocessed quiescent fraction for infalling galaxies and focus on the differential quenching (relative to the infalling population) in the final  = 0 cluster.
To compare our quenching timescales with Wright et al. (2022), we compare the time required for a galaxy of average sSFR (see their fig.7) to cross the quenching threshold and assume  delay = 0 to find an approximate average quenching time relative to the time of first pericentre.Wright et al. (2022) finds median quenching times 0.25⟨ orb ⟩ = 0.25(4.5 Gyr) = 1.1 Gyr and 0.65⟨ orb ⟩ = 0.65(3.5 Gyr) = 2.3 Gyr relative to the time of first pericentre, for their 9 < log( ★ /M ⊙ ) < 10 and 10 < log( ★ /M ⊙ ) < 11 bins, respectively.For the higher stellar mass galaxies, the majority (∼ 80 per cent) are quenched by second pericentre.We find longer quenching times of  Q = 3.7 ± 0.4 Gyr for 9 < log( ★ /M ⊙ ) < 10 galaxies and  Q = 4.0±0.2Gyr for 10 < log( ★ /M ⊙ ) < 10.5 galaxies.Using similar reasoning, their SFR suppression is equivalent to an exponential suppression time of  ∼ 0.4 Gyr and ∼ 0.9 Gyr for our lower and higher stellar mass bins, respectively.Taking  delay = 0, our  env is significantly longer, by 1.2 Gyr and 1.3 Gyr, respectively.Relaxing  delay = 0, we find little to no change in preferred  Q , but results in our preferred exponential suppression timescale being consistent with theirs for lower stellar mass galaxies.
Based on this discussion, we conclude that we prefer significantly longer total quenching timescales than Wright et al. (2022).If we assume  delay = 0, as their models predict, then we find that our star formation suppression timescales,  env , are longer than theirs by 1.2 − 1.3Gyr.
A similar previous study is that of Lotz et al. (2019), which instead examined  ∼ 0 quiescent fractions in the Magneticum Pathfinder hydrodynamical simulation.They found that most log( ★ /M ⊙ ) < 10.5 galaxies are quenched within ∼ 1 Gyr of crossing  200 (i.e.around the time of first pericentre), with the relatively small fraction of galaxies with tangential orbits and very high stellar masses able to maintain star formation after first pericentre.This quenching is significantly earlier than we find for both of our stellar mass bins, and all of the results examined in our discussion above.

Towards a consistent model of quenching
There is ample observational evidence that ram pressure stripping of the cold gas starts at or just before first pericentre.For example, Smith et al. (2010) found that ram-pressure stripped tails of Coma cluster galaxies were prevalent within half of the cluster virial radius and that most of these tails pointed away from the cluster, indicating that the stripping was occurring on infall, i.e. just before pericentre.Studies of the distribution of H I abundance as a function of the PPS coordinates (Jaffé et al. 2015;Oman et al. 2021) find that H I depletion begins close to pericentre, i.e within 0.5 200 , in agreement with Smith et al. (2010).However, quenching is not instantaneous and likely proceeds from outside inwards if it is due to ram-pressure stripping (Boselli et al. 2022).Owers et al. (2019) studied cluster galaxies with spatially resolved spectroscopy and uncovered a population of galaxies with strong H absorption -indicative of recent quenchingin their outskirts, but these same galaxies had ongoing star formation in their centres.They modeled the distribution of these galaxies in PPS, finding indications that galaxies with a recent quenching event in their outskirts are within 1 Gyr of entering within 0.5 200 of the cluster centre.
That quenching should start at (or just before) pericentre is supported by our results for our higher stellar mass bin ( delay = −0.3±+0.8  −1.0 Gyr) but is at odds with our results for low stellar mass galaxies.For these we find that the onset of quenching occurs well past pericentre  delay = 3.5± +0.6 −0.9 Gyr and with a fast quenching envelope ( env ≤ 1.0 Gyr).Such a short quenching timescale, however, would predict a deep 'green valley' in the sSFR distribution, inconsistent with the observed shallow depth (see Fig. 12).For the lower stellar mass bin, quenching starting at pericentre ( delay = 0) is permitted by the quiescent fraction but disfavoured by the stellar ages at the ≳ 2 level (see Fig. 9).This preference is driven by there being little-to-no gradient in ΔMWA between galaxies in the cluster core and those infalling, as shown in Fig. 4.This result, however, appears to be driven partly by the very lowest stellar mass galaxies: if we restrict analysis to 9.5 < log( ★ /M ⊙ ) < 10, we find somewhat better compatibility with  delay = 0.If we fix  delay = 0, as suggested by other observational and theoretical evidence, then, for the lower stellar mass bin, we find slower quenching with  env = 1.6 ± 0.2 Gyr.This timescale is in better agreement with the sSFR distribution in Fig. 12.For the higher stellar mass sample, the quenching timescale is longer, with  env = 2.2 ± 0.1 Gyr, and is in reasonable agreement with the sSFR distribution.While the SFR suppression timescales that we find are longer than those found in hydrodynamical simulations (see Section 4.5), they are shorter than the gas depletion timescales of 3.5-4 Gyr for field galaxies of comparable stellar mass (Boselli et al. 2014).
Taking all of these results together suggests a picture in which ram pressure stripping starts close to pericentre and is effective in a satellite galaxy's outskirts (where the restoring force is low compared to the force due to ram pressure) but may not be fully effective in their more tightly bound inner regions.For gas in the inner regions, while there may be no inflow of new cold gas (due to the complete stripping of the hot gas in the halo), the remaining cold gas will then be consumed by star formation.This consumption timescale is shorter than in the field for two reasons: first, because there is less cold gas available due to the stripping in the outskirts; and second, because ram pressure stripped galaxies have modestly-enhanced star formation rates at fixed stellar mass, compared to similar galaxies in the field, of 0.2-0.3dex (Roberts & Parker 2020;Roberts et al. 2022).

CONCLUSIONS
We have combined SDSS photometry and spectroscopy, orbital information from tracking haloes in an N-body simulation, and simulated galaxies from the UniverseMachine empirical model to constrain a simple model that suppresses star-formation histories of galaxies that have fallen into clusters.In the model, an infalling galaxy will have its star formation rate suppressed by an exponential envelope with timescale  env after some delay time  delay , relative to the time of the first pericentre.The parameter fits from this modeling give the mean quenching timescales for the infalling population as a whole.We jointly fit these model parameters using both the quiescent fraction and mean deviation from the mean MWA- ★ relation (ΔMWA) of the infalling population in projected phase space.Doing so allows us to break the degeneracy between time of quenching onset and quenching duration that was present in previous models.The method accounts for interloper galaxies (which appear in projection, but are not physically in the cluster), and the pre-processing of infalling galaxies, allowing us to isolate the quenching effect of the (most recent) infall into a massive cluster.Our main results can be summarized as follows: • The mean mass-weighted stellar age depends on location in projected phase space, with cluster member-dominated regions being older (by ≲ 1 Gyr) relative to an interloper-dominated region.
• Overall quenching times for our two stellar mass bins are driven by the quiescent fraction and are consistent with Oman et al. (2021), whose methodology we build on directly:  Q = 4.3 ± 0.4 Gyr and  Q = 3.9 ± 0.2 Gyr for our lower and higher stellar mass bins, respectively.We find longer overall quenching timescales than other works in the literature where only star-formation rates are modeled, but agree with Taranu et al. (2014), who make additional use of the age-sensitive Balmer lines of quiescent galaxies.
• Using mass-weighted ages allows us to break the degeneracy between  delay and  env .We find that the onset of quenching occurs at  delay = 3.5 +0.6 −0.9 Gyr and  delay = −0.3+0.8 −1.0 Gyr, relative to time of first pericentre, for galaxies in our 9 < log( ★ /M ⊙ ) < 10 and 10 < log( ★ /M ⊙ ) < 10.5 stellar mass bins, respectively.The models prefer a short SFR suppression timescale,  env ≤ 1.0 Gyr (consistent with ram-pressure stripping), for our lower stellar mass bin, and a longer 2.3 +0.5 −0.4 Gyr (consistent with strangulation) for our higher stellar mass bin.
• In contrast to Wetzel et al. (2013), our model is able to reproduce the SFR bimodality even with long exponential suppression timescales, thanks to the stochasticity of the UniverseMachine star formation histories that we employ (as opposed to using smooth analytic star formation histories).We note that, for our lower stellar mass bin, the depth of the green valley prefers values of  env ≳ 1.5 Gyr and  delay ∼ 0, in slight tension with the later quenching onset preferred by ΔMWAs.
Based on these findings and on our detailed discussion of the literature, we argue that satellites infalling into clusters experience ram pressure stripping of cold gas starting close to pericentre, which is only effective in the galaxy's outskirts, at least on the first pericentre passage.This leaves reduced cold gas available for continued star formation, resulting in star formation suppression timescales of  env ∼ 2 Gyr -longer than if galaxies were fully stripped on their first pericentre passage, but shorter than a simple starvation scenario.
Future surveys like the Bright Galaxy Survey of the Dark Energy Spectroscopic Instrument (DESI) will provide a large increase in sample size over that of the one million galaxies in the SDSS DR14 main galaxy sample (Ruiz-Macias et al. 2021).An increase in sample size of galaxies with spectroscopically-derived MWAs should reduce the errors on ΔMWA significantly.Using spectroscopically-derived quantities more sensitive to recent SFR suppression or quenching, such as the time at which 90 per cent of the stellar mass has formed (Webb et al. 2020;Upadhyay et al. 2021), could also provide additional constraining power on infall-related and general quenching models.The models could also be improved by using a physicallymotivated model of ram-pressure stripping, rather than a generic timescale for SFR decline.Such a model could involve radius and velocity at first pericentre, as suggested in a simple model by Owers et al. (2019, see also Roberts et al. 2019 for a quenching model depending on ICM density).With these various improvements, tighter constraints on time of quenching onset and duration via infall quenching as a function of stellar mass should be possible in the future.SDSS sSFR cut.The value of this correction at  = 0 is such that the increase or decrease in SFR forms a broken power law with slope −1 for log( ★ /M ⊙ ) > 10 and slope −0.7 for log( ★ /M ⊙ ) < 10, with the value of the correction set to 1.0 at log( ★ /M ⊙ ) = 10.3.This corresponds to a boost in the star formation histories for log( ★ /M ⊙ ) < 10.3 and suppression for log( ★ /M ⊙ ) > 10.3.The correction is applied to the 'true' rather than 'observed' Uni-verseMachine SFHs as well as the 'true' final SFRs and stellar masses (see Behroozi et al. 2019, for these definitions).A shift plus scatter is then applied to approximate observed versions of the quantities afterwards, according to the prescription used in Behroozi et al. (2019).
We find that this adjustment only results in minor effects on the best-fitting  delay and  env parameters.Specifically, the preferred  delay changes by +0.2 Gyr and  env by −0.2 Gyr for our lower stellar mass bin, or in terms of  Q , a shift of 0.5 Gyr.For the 10 < log( ★ /M ⊙ ) < 10.5 stellar mass bin, the corresponding shifts in  delay and  env are +0.1 Gyr and −0.5 Gyr, respectively.Given our uncertainties of ∼ 0.5 Gyr for  env and ≲ 1.0 Gyr for  delay (see Table 1), we conclude that any impact on our parameter constraints from changes in overall star formation history is negligible and therefore does not influence our discussion and qualitative conclusions.

APPENDIX B: DETAILED MODELLING PREDICTIONS FOR HIGHER STELLAR MASS BIN
Analogous to Fig. 7 and Fig. 8, we show our detailed model predictions for 10 < log( ★ /M ⊙ ) < 10.5 galaxies in bins of PPS for  Q and ΔMWA in Fig. B1.Similar to our lower stellar mass bin, a somewhat different range of  delay and  env values are preferred by  Q and ΔMWA.The joint best-fit for these parameters using both observables is given in Table 1, with the contours shown in Fig. 9.
Figure1.The fraction of interlopers in (R,V) projected phase space bins determined using the orbital library fromOman et al. (2021). is distance from cluster centre, in units of  vir,3D and  is the velocity offset from cluster centre, with units of  3D .The red diagonal line,  = − (4/3)  +2, guides the eye as to where ∼ 50 per cent of galaxies are interlopers, aside from a portion around (,  ) ∼ (1.2, 0.4), which has many galaxies at or close to their first apocentre after entering the cluster.The dotted line delineates a region in the upper right that is heavily dominated by the interloper population, which we use to compare between SDSS and UniverseMachine interlopers.

Figure 2 .
Figure2.Quiescent fraction of SDSS galaxies in the stellar mass ranges 9 < log(  ★ /M ⊙ ) < 10 (left) and 10 < log(  ★ /M ⊙ ) < 10.5 (right) in clusters with host halo mass log(  vir /M ⊙ ) > 14, binned in PPS coordinates.The quiescent fraction is significantly higher in cluster centres than in their outskirts (≳ 30 per cent).For the lower stellar mass bin, regions within the ∼ 50 per cent interloper (red diagonal) line, show at least modestly elevated quiescent fractions, with little evidence of enhanced quenching immediately outside of (above) this line.For the higher stellar mass bin, some higher velocity galaxies within the virial radius also show somewhat elevated quiescent fractions.

Figure 3 .
Figure 3. Top panels: SDSS observationally-derived mass-weighted ages for quiescent galaxies (left) exhibit a high degree of scatter and trend older with increasing galaxy stellar mass.UniverseMachine simulated galaxies (right) show a similar trend in mass-weighted age with stellar mass, albeit with a ∼ 3 Gyr offset and significantly less scatter (see Section 2.4.2 for discussion).SDSS points indicate the best-fitting mass-weighted ages for SDSS galaxies from Wilkinson et al. (2017, stellar masses are from Comparat et al. 2017), with the running mean (median) shown as a dashed (dotted) line.The functional fits of the mean MWA- ★ relation used in the analysis of this work is fit to the running mean mass-weighted age at a given stellar mass and is shown with the solid black curves.Bottom panels: Residual MWA from the respective MWA-M ★ functional fits for SDSS (left) and UniverseMachine (right) is shown along with the running mean (dashed) and median (dotted) of the residual.

Figure 4 .
Figure 4. Mean deviation in MWA from the MWA- ★ relation of an SDSS interloper proxy sample (top right three bins; outlined by dotted black lines), mean(MWA − MWA SDSS ), for quiescent 9 < log(  ★ /M ⊙ ) < 10 (left) and 10 < log(  ★ /M ⊙ ) < 10.5 (right) SDSS galaxies as a function of PPS coordinates.Errors on the mean values were calculated by bootstrapping over clusters, rather than individual galaxies.We note that the zero-point for this figure, defined by the three PPS bins in the top right corner are 97 per cent interloper galaxies according the orbit library ofOman et al. (2021).The dashed grey diagonal line indicates where approximately 50 per cent of galaxies are interlopers.

Figure 6 .
Figure 6.Delayed-then-rapid (instantaneous quenching, i.e.  env = 0) model predictions for a range of delay times,  delay , after time of first pericentre (a negative  delay corresponds to quenching prior to  fp ).For these illustrative plots, we use the UniverseMachine 'observed' interloper proxy sample as defined in Section 2.4.2, for clusters with log(  vir /M ⊙ ) > 14.Top row: predicted quiescent fractions,  Q , in PPS.The diagonal dashed line,  = − (4/3)  + 2, marks the location where ∼ 50 per cent of galaxies are interlopers.Bottom row: quiescent mass-weighted age predictions in PPS.

Figure 7 .
Figure7.Quiescent fraction,  Q , predictions for 9 < log(  ★ /M ⊙ ) < 10 galaxies for a range of models where galaxies quench after some delay time,  delay , relative to the time of first pericentre.We show models for a range of exponential suppression timescales.The SDSS mean  Q values are shown as grey lines, with the shaded regions showing the bootstrapped (over clusters) uncertainty in the mean.
Figure8.Mean ΔMWA predictions for quiescent 9 < log(  ★ /M ⊙ ) < 10 galaxies for a range of models where galaxies quench after some delay time, relative to time of first pericentre.We show models run for a range of exponential suppression timescales.The SDSS mean ΔMWA values are shown as grey lines, with the shaded regions showing the bootstrapped (over clusters) uncertainty in the mean.Note that our models include an additional offset to allow for the possibility that the infalling population may be different than that of the interloper phase space bins in Eq. 6.

Figure 11 .
Figure 11.Our model displays a persistent bimodality in sSFR distribution, regardless of SFR suppression timescale, whereasWetzel et al. (2013) finds a unimodal distribution for long suppression times (contrast the red solid line for their model with the blue solid line for our model).This arises because of our use of a stochastic star formation history, whereas they use smooth star formation histories.In particular, we are contrasting with fig.9inWetzel et al. (2013), for our 10 < log(  ★ /M ⊙ ) < 10.5 stellar mass bin and their 10.1 < log(  ★ /M ⊙ ) < 10.5 stellar mass bin.A correction of −1 Gyr is added to be comparable to their average infall time, which has the time zero-point set at  vir .Additionally, for ease of comparison withWetzel et al. (2013), who claimed need for a short exponential suppression timescale to reproduce the bimodality, we include only galaxies within  200c ≈ 0.73 360m and add log-normal scatter with mean log(sSFR/yr −1 ) = −12 and 0.25 dex variance to all galaxies with log(sSFR/yr −1 ) < −12.

Figure A1 .
Figure A1.Plots illustrating the adjusted UniverseMachine interlopers (red) relative to the SDSS data (black).Above log(  ★ /M ⊙ ) = 10.2, galaxies have had their star formation history suppressed to increase their  Q to match the SDSS data, whereas galaxies below log(  ★ /M ⊙ ) = 10.2 have had their star formation boosted.Left: specific star formation rate plotted against stellar mass, with the adjusted UniverseMachine sSFR multiplied by 2.1 (vertical shift) for the purposes of comparing with SDSS data on this plot.The quiescent/star-forming cut used for both is shown with the solid black line.Right: quiescent fraction versus stellar mass.The adjusted UniverseMachine interloper sSFRs increase their  Q (dashed red) up to that of the SDSS data for higher stellar masses.There is a small offset in the two curves as only the shape of the  Q − M * relation was fit in this stellar mass binning and the absolute  Q values were fit for stellar mass bins with edges at log(  ★ /M ⊙ ) = (9, 10, 10.5, 11, 11.5).
Illustration of the effect of our delayed-then-exponential envelope (thick blue line) on SFR, where at some time delay,  delay , after time of first pericentre upon infall into a cluster,  fp , SFR is quenched exponentially by a multiplicative envelope with timescale,  env .
dark matter simulation, including catalogues for complete star formation histories, are available on the following page of Peter Behroozi's website: https://www.peterbehroozi.com/data.html.

Table 1 .
Best-fitting quenching model parameters, given as the 50 th percentile for each respective marginal posterior probability distribution.Upper and lower uncertainties are the 84 th and 16 th percentiles of the marginal distributions.Top rows: The time delay until the onset of quenching relative to the time of a galaxy's first pericentre ( delay ), and the exponential star formation suppression timescale after the onset of quenching ( env ).These are the best-fitting values indicated by the black diamond markers for  delay and  env on the subplots of Fig.9.Bottom rows: Best-fitting infalling population values, as described for  Q in Eq. 11 and ΔMWA in Eq. 16. † indicates that the value is a 95 th -percentile upper limit.Parameter 9 < log(  ★ /M ⊙ ) < 10 10 < log(  ★ /M ⊙ ) < 10.5

9
< log( ★ /M ⊙ ) < 10 Similar to Fig.11, except we now demonstrate how we can reproduce the bimodality in the sSFR distributions of the SDSS data (solid black histogram) for a range of relevant models with different exponential quenching timescales (coloured lines) for our two stellar mass bins, 9 < log(  ★ /M ⊙ ) < 10 (left) and 10 < log(  ★ /M ⊙ ) < 10