© The Authors 2026

1. Introduction

Galaxy clusters have long been recognized as key environments for studying the impact of large-scale structure on galaxy evolution. As the most massive gravitationally bound systems in the Universe, clusters host deep potential wells filled with a hot intracluster medium (ICM), high galaxy densities, and a complex accretion history of groups and field galaxies. These conditions make them natural laboratories to investigate the ensemble of external (environmental) processes that drive galaxy evolution: observations over the past decades have revealed that the galaxy population in clusters differs markedly from that in the field both in terms of morphology and star formation activity. Furthermore, these differences are found to evolve significantly over cosmic time.

One of the most prominent signatures of cluster-driven evolution is the transformation of late-type spirals into lenticular (S0) galaxies. Early studies of the morphology-density relation (Oemler 1974; Dressler 1980) established that local clusters are dominated by early-type systems, with spirals contributing only a minor fraction of the population. In contrast, later studies based on Hubble Space Telescope (HST) observations (e.g. Dressler et al. 1997; Fasano et al. 2000; Postman et al. 2005; Desai et al. 2007) have shown that spirals are proportionally much more common (a factor of ∼2 − 3), while S0s are much rarer in distant (up to z ≃ 0.8) clusters, as compared to those nearby. This strongly supports an evolutionary scenario whereby the cluster environment acts on newly accreted late-type galaxies, driving a morphological transition towards earlier galaxy types on timescales of a few billion years (Gyr; e.g. Poggianti et al. 2009b; D’Onofrio et al. 2015). The main mechanism responsible for this transformation has yet to be identified and figures among the subjects of this study.

The other key effect associated with galaxy clusters is environment-driven quenching, whereby star-forming spirals infalling from the field are deprived of their cold gas reservoirs, first reducing their star formation (Vulcani et al. 2010; Paccagnella et al. 2016; Old et al. 2020) and eventually driving their evolution into passive systems. Amongst the various quenching mechanisms, ram-pressure stripping (RPS; Gunn et al. 1972) has emerged as one of the most effective: as galaxies plunge through the dense ICM at high velocities, their cold gas is displaced and removed, often on timescales shorter than a Gyr (e.g. Boselli et al. 2006; Marasco et al. 2016; Jaffé et al. 2018). The presence of a high fraction of post-starburst galaxies in cluster centres strongly supports RPS as a key evolutionary mechanism (e.g. Poggianti et al. 1999; Paccagnella et al. 2017, 2019). The most spectacular manifestation of ongoing RPS is the population of so-called ‘jellyfish galaxies’, which exhibit one-sided gaseous and stellar tails stretching tens of kiloparsecs (kpc; e.g. Smith et al. 2010b; Poggianti et al. 2017; Boselli et al. 2022; Poggianti et al. 2025), providing unique snapshots of environmental quenching in action. The observational picture across redshifts indicates that environmental quenching is already well underway by z ≃ 1 and beyond, when clusters are known to host a substantial fraction of passive galaxies, while still retaining a higher proportion of spirals compared to the present day (Poggianti et al. 2009a; Webb et al. 2020; McNab et al. 2021; Baxter et al. 2022). This implies that morphological transformation and quenching act on different timescales, with the former process being slower than the latter, or simply becoming more efficient at later times (e.g. Poggianti et al. 1999; Vulcani et al. 2015).

Observational and theoretical arguments indicate that RPS and morphological transformation in clusters are not independent processes. Once a galaxy is deprived of its dynamically coldest component (i.e. the interstellar medium), resonances induced by spiral arms are expected to drive up the level of random motion in the stellar disc, which becomes unable to support the spiral pattern itself (e.g. Sellwood 2014, and references therein). Ultimately, this ‘secular’ evolution alters the stellar structure of the galaxy and pushes it towards the S0 type. However, there is no consensus on the characteristic timescale of this mechanism, with values quoted in the literature ranging from ∼2 to 10 Gyr (e.g. Sellwood & Carlberg 1984; Fujii et al. 2011).

More recently, Marasco et al. (2023, hereafter M23) have shown that the simple ageing of the stellar populations in a spiral galaxy that is abruptly quenched by the cluster environment leads to a spatial and spectral redistribution of its stellar light that might be sufficient to alter its morphological classification. Once the star formation is halted, the galaxy becomes progressively redder, more concentrated, and smoother as the blue light originating from its spiral arms fades away. M23 quantified the timescale of this transformation by building and classifying a collection of ‘artificially aged’ synthetic images of 91 galaxies taken from the GAs Stripping Phenomena in galaxies (GASP) project (Poggianti et al. 2017), an integral-field spectroscopic survey with VLT/MUSE aimed at studying gas removal processes in galaxies. They found that the morphology transformation into an S0 is completed after 1.5 − 3.5 Gyr, proceeding faster in more efficient quenching scenarios. These results have revealed a promising channel for the spiral-to-S0 transformation of cluster galaxies, providing a natural link between RPS, star formation quenching and morphological evolution. However, to which extent this process alone is responsible for the transformation of the cluster galaxy population over the last several Gyr is yet to be determined.

The goal of this study, which complements the work of M23, is to infer the characteristic timescales of cluster-driven quenching and morphological transformation required to match the photometric and morphological properties of cluster galaxies in the redshift range 0 < z ≲ 1. We pursue our goal by building simplified models for the evolution of cluster galaxies in the lambda cold dark matter (ΛCDM) framework, where the main environmental mechanisms leading to galaxy quenching and morphological transformation are included in a simple, parametric form. We compare this model with the data to infer the evolutionary timescales, which ultimately enables us to address the question of whether the spectrophotometric evolution is indeed the main channel for the spiral-to-S0 transformation in clusters.

This paper is structured as follows. In Section 2, we describe our approach, providing details of the dataset used, the evolutionary model of cluster galaxies, and the technique adopted to compare the two. Our results are given in Section 3. In Section 4, we discuss some possible variations in our modelling scheme and provide an interpretation of our results in the context of the study of M23, with a summary of the cited study given in Section 5. Throughout this paper, we adopt a flat ΛCDM cosmology with Ω_m, 0 = 0.3 and H₀ = 70 km s⁻¹ Mpc⁻¹, and a Kroupa (2001) initial mass function (IMF). Magnitudes and colours are expressed in mag units.

2. Method

Our method is based on the modelling of the photometric and morphological properties of the cluster galaxy population at different redshifts. The model, detailed in Section 2.2, follows the build-up of cluster galaxies as a function of redshift and is engineered to predict the B, V, and I-band magnitudes and the main morphology types (spiral–Spi, S0, and elliptical–Ell) of cluster galaxies. The parameters of the model are constrained using observational data, detailed in Section 2.1. In Fig. 1, we provide a flowchart to guide the reader through our method.

Fig. 1. Flowchart illustrating our method. White blocks show pre-existing dataset or software packages. Yellow blocks show the new products built for this study. The chart indicates also the Sections of the main text that provide details on a given part of the method.

2.1. Data

The two cluster samples used in this study, OmegaWINGS (Fasano et al. 2006; Gullieuszik et al. 2015; Moretti et al. 2014, 2017) and EDisCS (White et al. 2005), have three key properties that make them the optimal choice for our investigation. First, they sample two redshift bins that feature a substantial difference in the morphology distribution of the cluster galaxy population (e.g. Fasano et al. 2000; Vulcani et al. 2011a,b, 2014a), supporting the need for a significant cluster-driven evolution over the last several Gyr. Second, as shown in detail in Section 2.2.1, the typical cluster properties are such that the EDisCS clusters can be considered as progenitor of the OmegaWINGS clusters (Vulcani et al. 2014a). Third, the two samples have compatible morphological classification scheme, as discussed below and in Vulcani et al. (2011a, hereafter V11). These properties ensure that the OmegaWINGS and EDisCS data provide homogeneous and highly significant information to constrain our evolutionary models.

OmegaWINGS is a survey of 46 nearby (0.04 < z < 0.07) galaxy clusters and of their peripheries, imaged in the B and V bands with OmegaCAM (Kuijken 2011) at the VLT Survey Telescope (VLS; Capaccioli & Schipani 2011). Each cluster is covered beyond its virial radius R₂₀₀, extending the environment range covered by the parent WINGS Survey (Fasano et al. 2006) and allowing for the imaging of the surrounding infall region. Spectroscopic follow-ups for 36 OmegaWINGS clusters were taken with the AAOmega spectrograph (Smith et al. 2004; Sharp et al. 2006) at the Australian Astronomical Observatory (AAT) for redshift measurements and assignment of cluster memberships (Moretti et al. 2017). The morphological classification of OmegaWINGS galaxies was performed using the MORPHOT package (Fasano et al. 2012) by Vulcani et al. (2023), based on 21 morphological diagnostics, computed directly from the V-band images, to provide two independent classifications (one based on a maximum likelihood technique and the other one on a neural network machine). These, in turn, were used to infer the morphological type, T_M. Here, we focus on spectroscopically confirmed galaxies within 1 × R₂₀₀ cluster-centric distance and with robust T_M measurements and we assigned a simple morphology label to each galaxy (Ell: −5.5 < T_M < −4.25, S0: −4.25 ≤ T_M ≤ 0, Spi¹: T_M > 0). Our final OmegaWINGS sample is made of 7162 galaxies (2358 Spi, 3201 S0, 1603 Ell), complete down to m_V = 20 and with a median redshift of 0.055, which we took as the characteristic redshift of this sample. We stress that the morphological mix remains steady across the redshift range spanned by the sample, in spite of the fact that the spatial resolution varied by a factor of ∼1.7 from z = 0.04 to z = 0.07. The M_★ estimates were taken from Vulcani et al. (2022), and are based on the relation between M_★/L_B and the rest-frame (B − V) colour from Bell & de Jong (2001), converted to a Kroupa (2001) IMF.

EDisCS is a multi-wavelength photometric and spectroscopic survey of galaxies in clusters at 0.4 < z < 1, based on 20 fields selected from Las Campanas Distant Cluster Survey (LCDCS) catalogue (Gonzalez et al. 2001). For all 20 fields, deep optical multi-band photometry obtained with FORS2/VLT (White et al. 2005) and near-infrared (NIR) photometry obtained with SOFI/NTT was available. FORS2/VLT was additionally used to obtain spectroscopy for 19 of the fields (Halliday et al. 2004; Milvang-Jensen et al. 2008; Vulcani et al. 2012), covering up to R₂₀₀ for all clusters except one. For consistency with WINGS/OmegaWINGS, M_★ estimates were re-derived by V11 using the Bell & de Jong (2001) approach, with total magnitudes inferred from photo-z fits (Pelló et al. 2009) and rest-frame luminosities determined with the method of Rudnick et al. (2009). The EDisCS morphological classification is based on visual inspection of the ACS/HST mosaic F814W imaging, acquired for ten of the highest redshift clusters (Desai et al. 2007). The resulting sample was made up of 254 galaxies (156 Spi, 45 S0, and 53 Ell) with m_I ≲ 22 and a median redshift of ∼0.7, which we used as a characteristic redshift for this sample. We found no significant differences in the morphological mix of EDisCS galaxies above and below z = 0.7, in spite of the F814W filter sampling different rest-frame wavebands across the sample.

Our investigation requires the high- and low-z galaxy samples to have compatible morphology estimates. This was tested in V11, where three of the authors originally performing the EDisCS classification provided a visual classification for a subsample of the WINGS galaxies. They found a very good agreement between the visual and the MORPHOT morphological types, with discrepancies that were compatible with those emerging from comparisons of the three visual classifications amongst each other. We notice that at the typical z of the EDisCS galaxies, the F814W band samples the the rest-frame V-band and provides a spatial resolution (in kpc) comparable with the OmegaCAM images of OmegaWINGS clusters, which is key for a consistent classification between the two samples.

For this study, we used the EDisCS and OmegaWINGS datasets to derive the distribution of Spi and S0 galaxies in a colour-M_★ plane, using the observed-frame (B − V) colour for OmegaWINGS and the (V − I) colour for EDisCS. We focussed on a mass range which is expected to be complete for OmegaWINGS (M_★ ≳ 10^9.8) but we anticipated that in our comparison with the models, we would end up normalising the predicted M_★ distributions for Spi and S0 to those observed in order to minimise issues due to mass incompleteness (see Section 2.3). We accounted for spectroscopic incompleteness in OmegaWINGS and EDisCS using the available correction weights w associated with each cluster member², although we stress that this correction does not impact our results. The resulting 2D distributions are shown in Fig. 4 and discussed in Section 3. They were used to constrain our evolutionary model, as we describe in detail below.

2.2. Model

2.2.1. Cluster growth

The first step in our modelling procedure was to determine cluster mass growth rates as a function of time. This is achieved using the halo growth model of Correa et al. (2015)³, based on the extended Press & Schechter (1974) formalism. This model provides a very good match to the halo accretion histories inferred from more detailed semi-analytical models of galaxy evolution, such as the Millennium (McBride et al. 2009) and the Bolshoi (van den Bosch et al. 2014) suites. The model is fully constrained by assuming a specific halo mass at given redshift: we focussed on a single cluster with a mass M_cl of 5 × 10¹⁴ M_⊙ at z = 0.055, which are the median mass and redshift of OmegaWINGS clusters⁴.

We computed the growth history of this halo from z = 2 to z = 0.055, displaying the results in the top panel of Fig. 2. M_cl has grown by a factor of ∼2 from z = 0.7 to z = 0.055. We inferred M_cl for the EDisCS and OmegaWINGS clusters from their R₂₀₀ estimates (Vulcani et al. 2013; Biviano et al. 2017) and plotted the median value of their distributions in Fig. 2, together with the difference between the 84th and 16th percentiles (blue square and red circle with error-bars). Clearly, there is excellent agreement between our growth model and the typical EDisCS and OmegaWINGS cluster masses: this allows us to treat our high- and low-redshift data within a common evolutionary framework (see also Vulcani et al. 2014b). Additionally, Fig. 2 shows that M_cl at z = 2 was about a factor of 5 and 10 lower than at z = 0.7 and z = 0.055, respectively. This is a key aspect of our method because it ensures that the modelled cluster galaxy population at z = 0.7 and z = 0.055 will be dominated by systems acquired from the field and evolved according to the recipes implemented; rather than by the initial (proto)cluster population at z ≃ 2, whose properties are largely unconstrained. In other words, at the redshifts where the model is compared with the observations, it will have ‘lost memory’ of its (poorly constrained) initial conditions.

Fig. 2. Evolutionary properties of our model cluster. Top panel: Mass (solid black curve) and accretion rate (dashed grey curve) vs redshift for a cluster with Mcl = 5 × 1014 M⊙ at z = 0.055. The blue square and the red circle show the median Mcl and redshifts of the EDisCS and OmegaWINGS clusters, respectively, with the error bars corresponding to half the difference between the 84th and 16th percentiles of the galaxy distribution. Bottom panel: Morphology fractions for field galaxies that join the cluster as a function of z, assumed to be independent of the galaxy M★.

As a next step, we populated the halo with galaxies. These were injected within the cluster from the field in time-steps of 0.1 Gyr, following the redshift-dependent galaxy stellar mass function of McLeod et al. (2021). We assumed that the cluster galaxy mass growth rate is proportional to the halo mass growth rate; that is, if in the time interval from t to t + Δt, the halo mass has increased by $\mathrm{\Delta }{M}_{\mathrm{cl}}\equiv {\dot{M}}_{\mathrm{cl}}\mathrm{\Delta }t$, an additional ΔM_★ ≡ f_★ΔM_cl will be added to the cluster, randomly sampled from the assumed galaxy mass function at redshift z(t). In principle, f_★ can be chosen so that the final number of cluster galaxies above a given M_★ (or below a given magnitude) is similar to that of the typical OmegaWINGS cluster. However, as we were not interested in modelling individual clusters, but, rather, in acquiring the properties (mass, colour, and morphological mix) of the mean galaxy population, we can arbitrarily inflate f_★ to ensure sufficiently large statistics. In particular, we adopted f_★ = 0.6, which for our best-fit model (Section 3.1) provides ∼3 × 10³ disc galaxies (Spi+S0s) with m_I < 22 at z = 0.7, and ∼10⁴ with m_V < 20 at z = 0.055. These are a factor of 5–10 greater than the number of discs in the EDisCS and OmegaWINGS datasets, ensuring a robust statistical characterisation of the model galaxy population.

2.2.2. Field morphology

Our model requires us to assign an initial morphology to all galaxies that join the cluster from the field. Ideally, we would need the morphology distribution of field galaxies as a function of M_★ and z, which we would use to randomly assign an initial type to all galaxies in the model. This is a challenging request: while observational constraints are numerous both in the Local Universe (e.g. Nair & Abraham 2010; Wilman & Erwin 2012; Moffett et al. 2016; Domínguez Sánchez et al. 2018) and at high redshift (e.g. van der Wel et al. 2007; Oesch et al. 2010; Huertas-Company et al. 2016; Cavanagh et al. 2023; Kartaltepe et al. 2023), they are somewhat inhomogeneous and frequently in tension with one another. Our strategy is to focus on field morphology estimates that are compatible with those adopted for our cluster galaxies. We relied on the study by Calvi et al. (2012), who analysed a sample of ∼2000 systems from the Padova–Millennium Galaxy and Group Catalogue (PM2GC; Calvi et al. 2011), using MORPHOT (as in OmegaWINGS) to assign morphological types. For their mass-limited sample (M_★ > 10^10.2 M_⊙), Calvi et al. (2012) found that the morphological mix shows only a weak dependence on M_★ and on the environment (with the exception of clusters). This resulted in 27.0% of Ell, 28.7% of S0s, and 44.3% of Spi in the so-called ‘general field’ environment (PM2-GF), which collects galaxies in groups, binary systems, and in isolation. The general field is built to represent the typical galaxy population that accretes onto nearby clusters; as most galaxies in the PM2-GF belong to groups, we would expect pre-processing effects to be accounted for. We used the percentages quoted above, valid for all M_★, to characterise the field population in our lowest redshift bin. For simplicity, we did not use the estimated uncertainties on these percentages, which are small in any case (1.3 − 1.5%).

The extrapolation of these fractions to higher redshifts requires additional assumptions. As our model does not treat the morphological evolution of ellipticals, the final Ell fraction will depend solely on that accreted from the field. There is observational evidence for an approximately constant ∼30% Ell fraction in clusters at least up to z ≃ 0.6 (Fasano et al. 2000; Vulcani et al. 2023, V11). To reproduce this result, we assumed that the Ell fraction adopted for the low-z field is valid at all z, which guarantees a 27% of Ell in the model cluster at all epochs. A convenient outcome of this assumption is that we can safely factor out the ellipticals from our modelling scheme and focus solely on the evolution of the disc (Spi+S0) population. Conversely, several studies are in agreement on the tendency for the field Spi fraction to increase with z (Dressler et al. 1997; Bundy et al. 2005, 2006; Conselice 2014); however, there is no consensus on the specific trend behind it. Here, we follow the results of Cavanagh et al. (2023, see their Fig. 7) and set the Spi fraction to increase linearly with z so that, by z ≃ 0.9, the disc (Spi+S0) field population is fully dominated by spirals. As this assumption has visible impacts on our results, but only marginally supported by the observational data, we discuss alternative scenarios in Section 4.4.

The bottom panel of Fig. 2 shows the assumed evolution of the field morphological fractions. The field fractions at z = 2 are also used to define the initial morphological mix within the (proto)cluster. This has virtually no impact on our results because, as discussed in Section 2.2.1, the cluster initial conditions are rapidly diluted during the evolution.

2.2.3. Star formation histories and model photometry

The detailed modelling of galaxy photometry and, in particular, of galaxy colours, is key to anchoring our models to the data. In order to have complete control over the photometric evolution of the model galaxies, we have chosen to fully characterise their spectral energy distribution (SED) using the Bayesian Analysis of Galaxies for Physical Inference and Parameter EStimation (BAGPIPES) software package (Carnall et al. 2018), adopting a series of assumptions that we discuss below.

The main ingredient to build the SED of a galaxy is its star formation history (SFH). We constrain the SFH of field spirals by assuming that, until the time at which they join the cluster, they have being lying on the (redshift-evolving) main sequence (MS) of star formation. The reference MS model that we adopt here is that of Popesso et al. (2023); however, in Appendix B, we show that the MS of Speagle et al. (2014) leads to compatible results. Following Gladders et al. (2013), we parametrised the galaxy SFHs using a log-normal function, which is known to aptly describe the evolution of the specific SFR of galaxies with time and the overall star formation rate density of the Universe. In Appendix A, we describe how to derive a relation between the present-day M_★ of a galaxy and the two free parameters of its log-normal SFH, so that the galaxy remains approximately on the MS for 0 < z < 2. As shown in Appendix A, galaxy downsizing arises naturally under these assumptions.

BAGPIPES uses the Bruzual & Charlot (2003) stellar population synthesis (SPS) models (2016 version) with the Kroupa & Boily (2002) initial mass function. The M_★ estimates include stellar remnants, but exclude the gas lost by stellar winds and supernovae. We assigned a unique value for the gas-phase metallicity of any galaxy based on its M_★ and SFR using the fundamental metallicity relation (FMR) of Curti et al. (2020), assumed to be redshift-independent (e.g. Mannucci et al. 2010). We adopted the dust attenuation model of Charlot & Fall (2000), using a value of 0.7 for the power law slope of the attenuation law and a multiplicative factor of 3 for the A_V associated with the birth clouds (which in BAGPIPES have a lifetime of 10 Myr). A visual extinction A_V was assigned to each field spiral following the prescription of Garn & Best (2010), which is based on the galaxy M_★ alone. BAGPIPES was used to infer observed- and rest-frame apparent and absolute magnitudes in the B, V and I bands, using the appropriate OmegaCAM and FORS2 broad-band filters to best match the observed photometric measurements.

2.2.4. RPS and cluster-driven quenching

One of the main photometric features of the cluster galaxy population is the presence of a well-populated red sequence, built of galaxies of various morphological types whose star formation processes have been halted. We modelled environmental quenching in a scenario where its main physical driver is RPS from the ICM, which removes the galaxy’s ISM,thereby inhibiting the star formation processes. As in M23, the stripping proceeds outside-in at an exponential rate, reducing the radius of the ISM disc as

$$\begin{array}{c}\hfill R(t)=R_{\max }exp(-t/t_{\mathrm{s}}),\end{array}$$(1)

where t is the time since the galaxy has joined the cluster, R_max is the disc radius at t = 0, and t_s is the stripping timescale. This exponential form is meant to capture qualitatively the increasing stripping efficiency affecting a galaxy in acceleration within a progressively denser ICM, and allows for a direct comparison with the results of M23, discussed in Section 4.3. We can compute analytically how the SFR drops in this stripping scenario under the assumptions of pure exponential discs and that the SFR density radial profile is proportional to the gas surface density profile (which implies a unitary slope in the Kennicutt-Schmidt relation, as suggested by Leroy et al. 2013). We find

$$\begin{array}{c}\hfill \frac{\mathrm{SFR}(t)}{\mathrm{SFR}(0)}=\frac{1-e^{-\delta (t)}[1+\delta (t)]}{1-e^{-\delta (0)}[1+\delta (0)]},\end{array}$$(2)

where δ(t)≡R(t)/R_d, with R(t) given by Equation (1) and R_d being the exponential scale length of the disc. We can remove the dependence on R_d by imposing an initial disc radius proportional to the scale length; specifically, we set R_max = 5 × R_d, a size that encloses ∼95% of the mass for an exponential profile. Consequently, the drop in SFR modelled by Equation (2) is regulated by (t/t_s) alone. As a reference, from the moment the galaxy joins the cluster, the SFR drops to ∼60% of its initial value after t_s, to ∼15% after 2 × t_s, and to ∼3% after 3 × t_s.

The stripping timescale is one of the main parameter of our model. We assumed a scenario along the lines of that described by Gunn et al. (1972), where the stripping efficiency is regulated by the competition between the gravitational restoring force of the galaxy. In that case, it depends (at least in the proximity of the disc) on M_★, and the cluster gas density and velocity dispersion, both of which depend on M_cl. In this scenario, we can parametrise t_s as a power-law function of the galaxy-to-cluster mass ratio,

$$\begin{array}{c}\hfill {log}_{10}\left(\frac{t_{\mathrm{s}}}{\mathrm{Gyr}}\right)={log}_{10}\left(\frac{{\tau }_{\mathrm{s}}}{\mathrm{Gyr}}\right)+{\alpha }_{\mathrm{s}}[{log}_{10}\left(\frac{M_{★}}{M_{\mathrm{cl}}}\right)+4.7],\end{array}$$(3)

where τ_s and α_s are free parameters. The 4.7 addend on the right-hand side of Equation (3) is chosen so that τ_s is the stripping timescale for a galaxy with log₁₀(M_★/ M_⊙) = 10 in a cluster with log₁₀(M_cl/ M_⊙) = 14.7, the typical M_cl of the OmegaWINGS clusters.

We stress that Equation (3) does not come from a rigorous treatment of the RPS physics in galaxy clusters, but it is a rather simplified and convenient way to outline the trends between t_s, M_★ and M_cl based on our physical intuition. We notice that the RPS efficiency is expected to depend also on the cluster concentration, c_cl. However, in the framework assumed, c_cl can be expressed as a power-law function of M_cl using the halo concentration-mass relations (e.g. Dutton & Macciò 2014); thus, such a dependency is effectively reabsorbed in the α_s parameter of Equation (3).

We tracked the SFH of galaxies subject to this cluster-driven quenching and update their model photometry using the same approach of Section 2.2.3. Similarly to M23, we included the effect of shock-driven dust destruction by lowering the extinction A_V by a factor exp[−t/(t_s + t_dust)], with t_dust = 100  Myr (Jones et al. 1994, 1996). Dust destruction makes galaxies bluer, marginally slowing down their pathway to the red sequence. Also, we considered the stellar metallicity frozen to the value the galaxy had when it joined the cluster. This assumption simplifies our calculations and is justified by the mild (∼0.1 dex) increase in (light-weighted) stellar metallicity that quiescent galaxies show with respect to active galaxies of similar M_★ and redshift (e.g. Gallazzi et al. 2026).

In Fig. 3, we illustrate the case of a representative MS galaxy with initial M_★ ≃ 2 × 10¹⁰ M_⊙ that joins the cluster at z = 1. We follow the mass build-up and the photometric evolution of this system down to z = 0 for three different scenarios: fast stripping with t_s = 0.5 Gyr, slow stripping with t_s = 2.0 Gyr, and no stripping at all. The top panel of Fig. 3 shows that, even in the no-stripping scenario, the galaxy SFR has dropped by a factor ∼8 since z = 1, in agreement with the evolution of the MS and of the cosmic SFR density (e.g. Madau & Dickinson 2014; Enia et al. 2022). Regardless of the scenario considered, the galaxy rest-frame (B − V) increases with time (bottom panel of Fig. 3) due to the combined effect of the ageing of the stellar populations and of the overall decline of the SFR. The no-stripping scenario is the one that provides the bluest and most massive system at z = 0, as expected from galaxies that evolve on the main sequence of star formation. The fast stripping model produces the reddest and least massive system; in fact, it is the only scenario that brings the galaxy on the red sequence, whereas the slow stripping leads to masses and colours that are in between the other two. These illustrative experiments already indicate that short stripping timescales are needed to build a well populated red-sequence of cluster galaxies.

Fig. 3. Photometric evolution and mass build-up of an idealised main-sequence galaxy with initial M★ ≃ 2 × 1010 M⊙ that joins our cluster model at z = 1, for stripping timescales of 0.5 and 2.0 Gyr (dark-red and orange curves, respectively), and for a case without stripping (blue curve). Top panel: SFR as a function of the look-back time. For times larger than ∼8 Gyr (z > 1), the three models share a common SFH shown by the solid black curve. Bottom panel: Galaxy evolution in the rest-frame (B − V) colour vs M★ plane. The white square shows the initial values at z = 1, the coloured circles show the final (z = 0) values for the three scenarios considered. The green dashed-line separates the red sequence from the blue cloud (from Vulcani et al. 2022).

The S0s that join the cluster from the field are treated differently, under the assumption that they are already quenched galaxies and will evolve passively within the cluster environment. We assumed that they were MS galaxies in the past and that their quenching had begun exactly Δt_S0 = 1 Gyr before joining the cluster, occurring on very short time scales equivalently to t_s = 0.1 Gyr in our model. Although observed early type galaxies show a broad variety in their quenching epochs and timescales (e.g. Tacchella et al. 2022), the idealised quenching scenario adopted here is meant to produce galaxies that (by the time they enter the cluster) are already ‘red-and-dead’ regardless of their specific quenching channel. Increasing Δt_S0 and/or further reducing t_s have negligible impacts on the colours of our field S0 population, which is the main observable used in our analysis.

2.2.5. Cluster-driven morphological transformation

Finally, we introduce the cluster-driven morphological transformation in our models. Conversely to what we have assumed for the quenching, the approach that we adopt here is agnostic of the exact mechanism responsible for the transformation: our goal is to determine whether a transformative mechanism is needed to begin with and, if so, what is its characteristic timescale.

We treat the transformation with a statistical approach, where spirals that enter the cluster will progressively increase their probability p_S0 of evolving into an S0 up to a maximum probability value p_max. We parametrise p_S0 as

$$\begin{array}{c}\hfill p_{\mathrm{S}0}(t)=p_{\max }[1-exp(-t/t_{\mathrm{m}})]\end{array}$$(4)

where t_m is the morphological transformation timescale. We notice that the evolution is fully suppressed for p_max = 0. We impose p_S0 = 1 for all S0 field galaxies, for which we assume no morphological evolution within the cluster environment.

Although the parametric form chosen for p_S0 mimics the trends found by M23 for their spectrophotometric evolution of quenched galaxies (see, for instance, their Fig. 5), it is general enough to be applicable to any evolutionary channel. We assume a unique value of p_max for all galaxies at all redshifts, and consider it as a free parameter of the model. By analogy with our treatment of the stripping timescale (Section 2.2.3), we assume that t_m is a power-law function of M_★/M_cl,

$$\begin{array}{c}\hfill {log}_{10}\left(\frac{t_{\mathrm{m}}}{\mathrm{Gyr}}\right)={log}_{10}\left(\frac{{\tau }_{\mathrm{m}}}{\mathrm{Gyr}}\right)+{\alpha }_{\mathrm{m}}[{log}_{10}\left(\frac{M_{★}}{M_{\mathrm{cl}}}\right)+4.7],\end{array}$$(5)

where τ_m and α_m are additional free parameters of the model. In practice, in our model p_max can be considered as a regulator of the global efficiency of the morphological transformation, irrespectively of the galaxy or the cluster masses, which instead regulate the transformation timescales.

We notice that in our agnostic approach, we consider t_s and t_m as separate, independent timescales. However, in a scenario where the morphological transition follows from the galaxy quenching, we expect that t_m > t_s. We prefer not to force this inequality in our models, and expect it to arise naturally as a data-driven result.

2.3. Fitting technique

Our model uses five free parameters, listed in Table 1. Two regulate the cluster-driven stripping timescales (τ_s, α_s), and three regulate the cluster-driven morphological evolution (τ_s, α_s, p_max). These are constrained using the Spi and S0 distribution in two separate colour-M_★ planes (see Section 2.1 and Fig. 4), which our model can predict at any z for any combination of its free parameters, using the p_S0 (Equation (4)) associated with each galaxy as a weight to infer the S0 fraction.

Fig. 4. Comparison of our best-fit cluster model with the EDisCS data at z = 0.7 (top six panels), and with the OmegaWINGS data at z = 0.055 (bottom six panels). The data and the model are shown in the first and second row of each panel set, respectively. Panels in the leftmost (central) column show the distribution of spirals (S0s) in a colour-M★ plane, where the colour is the observed-frame (V − I) for EDisCS and the (B − V) for OmegaWINGS. The colour palette indicate the number of galaxies per bin. Cyan contours are drawn at 2, 5, 10, and 15 galaxies per bin in EDisCS and in the model at z = 0.7, and at 2, 10, 30, and 100 galaxies per bin in OmegaWINGS and in the model at z = 0.055. The total fraction of Spi and S0 galaxies with respect to the total disc (Spi+S0) population is indicated in the top-right corner of each panel. Panels in the rightmost column show the S0 fraction as a function of colour and M★, for bins that have at least five galaxies in the data.

These distributions are built in the same way for the model and the data: we consider the mass range 9.75 < log₁₀(M_★/M_⊙) < 11.25 in bins of 0.25 dex, the (observed-frame) colour range 0.3 < (B − V) < 1.5 in bins of 0.1 mag for the model at z = 0.055 and for the OmegaWINGS data, and the colour range 0.8 < (V − I) < 3.0 in bins of 0.2 mag for the model at z = 0.7 and for the EDisCS data. The use of larger colour bins for EDisCS is justified by the lower number of galaxies available compared with OmegaWINGS. For simplicity, in our models we focus solely on the two time ‘snapshots’ that best match the median redshifts of the two surveys, but we have verified that the use of broader redshift distributions that approximately match those observed does not alter our results.

Before building the 2D distributions, some adjustments to the model were needed in order to make the comparison with the data meaningful. First, we considered only galaxies in our model with m_V < 20 at z = 0.055 and with m_I < 22 at z = 0.7, to account for the magnitude completeness of the observed samples. Second, we weighted our model galaxies so that the disc (Spi+S0) galaxy M_★ function in the model matched the observed ones. This means that, by construction, the total number of model discs galaxies in any bin of M_★ will match the data perfectly, but the colour and morphological mix will still depend on the model parameters. This step allowed us to circumvent potential issues associated with mass incompleteness in the data and/or with an incorrect mass function in the model, which is required to implement the metric discussed below. Finally, we must account for measurement errors that broaden the observed galaxy distribution in the colour-M_★ plane. This is done by introducing random (Gaussian) fluctuations in the colour and M_★ of the model galaxy population. We assumed a 0.3 dex scatter in M_★ based on a comparison between the Bell & de Jong (2001) method and M_★ estimates from full SED modelling (e.g. V11). Following Marasco et al. (2025), we assumed a minimum uncertainty of 5% on all flux measurements, which is added (in quadrature) to the nominal photometric uncertainty leading to a final scatter between 0.05 and 0.07 mag (depending on the band) on our model photometry.

The metric adopted to quantitatively evaluate the differences between the data and the model, both of which are in the form of galaxy number counts in bins of colour and M_★, is based on Poisson statistics. Following Bonamente et al. (2025), we express the Poisson likelihood as

$$\begin{array}{c}\hfill ln\mathcal{L}(\mathbf{x})=k{\displaystyle \sum _{i=1}^N}[D_{i}ln{M}_i(\mathbf{x})-M_i(\mathbf{x})-ln(D_i!)],\end{array}$$(6)

where x are the free parameters of the model, D_i and M_i are data and model galaxy number count in the i-th bin of the colour-M_★ plane, the sum is extended to all N bins considered, and 0 < k ≤ 1 (see below). We computed Equation (6) four times (for the two redshift bins and for the Spi and S0 populations separately) and summed the resulting likelihoods together to obtain the total one.

In an idealised scenario where x can be adjusted so that the difference between the model and the data are driven by Poissonian errors alone, k = 1 would be recommended in Equation (6). However, considering all the assumptions and simplifications done in our approach, this scenario is not realistic and the typical differences between data and model would end up being larger. Lowering k allowed us to compensate for this, artificially boosting the uncertainties in the data and avoiding unrealistically small error-bars in the model parameters. We set k to an ad hoc value of 0.03, chosen so that models whose parameters differ by ±1σ rms from the best-fit values lead to visually appreciable differences in the 2D distributions.

We fit the model to the data in the context of Bayesian statistics using the Nautilus sampler (Lange 2023) to infer the posteriors, using uniform priors in the ranges reported in Table 1 and the likelihood discussed above. We distinguished between the ‘best-fit’ parameters, which are those that maximise our likelihood, and the ‘fiducial’ parameters, given by the median of the marginalised posterior distributions. Upper and lower uncertainties on the fiducial parameters are provided by the 84th and 16th percentiles of the marginalised posteriors.

3. Results

In this section, we show the comparison between the data and our best-fit model (Section 3.1), present our predictions for the RPS and morphological transformation timescales (Section 3.2), and for the evolution of the morphological mix in clusters in the latest ∼6 Gyr (Section 3.3).

3.1. Data versus best-fit model

Fig. 4 shows the comparison between the distribution of cluster galaxies in the colour-M_★ plane for the data and for our best-fit evolutionary model. We first discuss the comparison with EDisCS at z = 0.7, shown in the top-half of Fig. 4. The data show a dominance of Spi over the S0 population (78 ± 3% vs 22 ± 3% of the total disc population), in agreement with what previously found for the same sample by Desai et al. (2007) and V11, and with other studies at similar redshifts (e.g. Dressler et al. 1994; Smail et al. 1997; Fasano et al. 2000). The S0 population is built preferentially by high-M_★ galaxies with (V − I)≳2, while Spi galaxies occupy also the bottom-left corner of the plane at lower M_★ and bluer colours. Interestingly, the data indicate a somewhat bimodal distribution of Spi galaxies, with an over-density around 10.5 < log₁₀(M_★/M_⊙) < 10.75 and 2.0 < (V − I) < 2.5 (which is similar to the region occupied by the S0s) and a separate bluer sequence at log₁₀(M_★/M_⊙) < 10.5 and (B − I)≃1.5.

Compared to the data, our best-fit model predicts a marginally more abundant S0 population at this redshift (66% Spi, 34% S0), but it reproduces the observed trends quite well. In particular, it does feature separate blue and red sequences for the Spi population, although it slightly over-predicts the number of red Spi hence providing a less populated blue sequence. S0s in our model occupy a well defined region of the colour-M_★ plane, in excellent agreement with the data. The rightmost panels in Fig. 4 show the S0 fraction (S0/(S0+Spi)) as a function of the galaxy colour and M_★. In this panel we have considered only bins containing at least 5 disc galaxies, in order to have a minimum statistics. As expected, the data show a trend for the S0 fraction to increase as a function of M_★ and (V − I), reaching a maximum value of ∼0.5. Our model reproduces this trend quite faithfully, reaching similar maximum fractions.

We go on to focus on the bottom-half of Fig. 4, which shows the comparison with OmegaWINGS at z = 0.055. Compared with the high-z population, the low-z data tell a different story. The majority of disc galaxies are S0s (34 ± 1% Spi, 66 ± 1% S0). The S0 population is distributed across the entire M_★ range, with typical (B − V) colours that are a slowly increasing function of M_★: this trend is the well-known ‘ridge’ of the red sequence (e.g. Baldry et al. 2004; Vulcani et al. 2022). The Spi population appears to be distributed similarly, but features a wider tail towards bluer colours for all M_★ (the ‘blue cloud’). The combination of the observed trends leads to an S0 fraction that grows progressively from the blue cloud to the red sequence, reaching maximum values close to unity. The model provides an excellent match to all the observed trends: we clearly see the ridge of the red sequence and tail of the blue Spi building the blue cloud. The overall Spi and S0 fractions predicted by the model match the data almost perfectly. However, there is a minor mismatch between the observed and predicted S0 fraction distribution, with the model featuring higher values within the blue cloud and lower values within the red sequence. We argue that such mismatch could be amended with a minor (< 0.1 mag) increase in the (B − V) colour of the model S0 population, which would increase (decrease) the S0 fraction in the red sequence (blue cloud). We also notice that the model predicts a red sequence that is slightly shallower than that observed. We discuss this further in Section 4.2.

Overall, we find it remarkable that, in spite of all the assumptions and simplifications we made, our model provides a cluster galaxy population that is realistic in terms of mass, colour, and morphology distribution. This is a necessary condition for the investigation of the physics of the cluster-driven evolution, which we examine below.

3.2. RPS and Spi-to-S0 transformation timescales

The last two columns of Table 1 show the fiducial and the best-fit parameters of our model, while corner-plots showing the full posterior distributions are shown in Fig. B.1. We find p_max ≃ 0.8, that is, there is a strong evidence for the morphological evolution of the cluster spiral population: in order to reproduce the data, it is required that ∼80% of spirals that join the cluster evolve into S0 on timescales that we discuss below. We notice that the best-fit parameters of our model are within the quoted uncertainties of the fiducial ones, with the exceptions of τ_m and p_max. The latter takes of a value of 1 for the model that maximises the likelihood. We have attempted to re-fit our model by fixing p_max to 1, finding fiducial values compatible with those reported in Table 1. The same applies for models that use the MS of Speagle et al. (2014), instead of Popesso et al. (2023), as reported in Table B.1.

We can combine the posterior sampling provided by Nautilus with Equations (3) and (5) to investigate the ‘fiducial’ relation between t_s, t_m and M_★/M_cl emerging from our model, showing our results in Fig. 5. We find that t_m > t_s for all values of M_★/M_cl, implying that the morphological transformation occurs after ram-pressure stripping has removed the galaxy’s gas reservoir, in agreement with our expectations. The stripping timescale shows the strongest trend with the mass ratio, ranging from ∼1 Gyr for log₁₀(M_★/M_cl)≃ − 5 and rapidly declining down to ∼100 Myr at log₁₀(M_★/M_cl)≃ − 4. This implies that, for a fixed cluster mass of 5 × 10¹⁴ M_⊙, galaxies with M_★ ≃ 5 × 10⁹ M_⊙ would have t_s = 1 Gyr, and galaxies with M_★ ≃ 5 × 10¹⁰ M_⊙ would have t_s = 100 Myr. Overall, the range of t_s found here is largely compatible with those emerging from previous observational (e.g. Boselli et al. 2006; Jaffé et al. 2018) and theoretical (e.g. Roediger & Brüggen 2007; Tonnesen & Bryan 2009; Marasco et al. 2016) studies, and supports RPS as the most efficient quenching mechanisms in clusters. However, the anti-correlation between t_s and M_★/M_cl, caused by the negative α_s found in our fitting, is counter-intuitive as one would expect most massive galaxies to be less influenced by the cluster environment, at a fixed M_cl. In Section 4.2 we discuss the specific features in the data that drive this trend, and provide a possible interpretation in Section 4.5.

Fig. 5. Characteristic timescales for the ram-pressure stripping (blue) and for the cluster-driven morphological transformation (red) as a function of the galaxy-to-cluster mass ratio (M★/Mcl) derived with our analysis (Equations 3 and 5). The solid lines show the median trend, with the shaded regions showing the difference between the 16th and 84th percentiles. The dashed lines show the best-fit model. As a reference, in the horizontal axis on top we report the M★ scale for a fixed log10(Mcl/ M⊙) = 14.7.

Conversely, t_m shows a weaker correlation with M_★/M_cl, with typical values ranging from 1 to 4 Gyr, which qualitatively agree with the evolutionary timescales found by M23. A more detailed comparison between our findings and those of M23 is provided in Section 4.3. We notice that there is a factor ∼2 difference between the fiducial and the best-fit t_m (solid and dashed red lines in Fig. 5). This is likely driven by the (weak) degeneracy between p_max and τ_m, visible in the corner-plots of Fig. B.1, such that models with a more efficient morphological transformation (larger p_max) but longer timescales (larger τ_m) are similar to those with a less efficient transformation but shorter timescales.

3.3. Evolution of the morphological mix since z ≃ 0.6

We now focus on the evolution of the morphological mix in galaxy clusters. A key observational study in this context is that of Fasano et al. (2000), who investigated the morphological content of clusters over the redshift range 0 < z < 0.6. They used an heterogeneous sample made of 25 clusters at intermediate redshifts observed either with HST/WFPC2 or with high-quality ground-based imaging, complemented by a local comparison set of 22 nearby clusters drawn from Dressler (1980) catalogue. All galaxies were morphologically classified into Ell, S0s, and Spi, and particular care was devoted to ensuring consistency across the different datasets via cross-checks between HST and ground-based classifications.

In Fig. 6, we compare the evolution of the morphological fractions found by Fasano et al. (2000), complete for galaxies down to M_V = −20, with the prediction of our best fit model. Predictions that use the same magnitude completeness of the data (thick solid lines) are shown separately from those that include all galaxies with M_★ > 10⁹ M_⊙ available in the model, in order to highlight the importance of the photometric selection. The predictions for the model obtained with the fiducial parameters are very similar and are not shown here for brevity. The data show a clear trend with redshift, featuring a steady increase of the S0 fraction at the expenses of the Spi, with the fraction of Ell that remains approximately constant. Our model follows very similar trends, especially when the magnitude cut is applied, but with two main differences. First, our Ell fraction is fixed by construction to 27% (Section 2.2.2), which is valid for the local field population but is somewhat a lower limit for the cluster population. Second, the Spi-to-S0 evolution predicted by the model is slightly milder than that shown by the data, which, compared to the model, feature a higher (lower) S0 (Spi) fraction in local galaxy clusters. However, as our model reproduces with remarkable accuracy the Spi and S0 fractions of OmegaWINGS clusters, this issue appears to be associated either with limitations in the dataset chosen to constrain the model, which cover only the initial and final stages of the observed evolution, or with a substantial scatter in the cluster morphological mix, such that the pronounced trend visible in the data of Fasano et al. (2000) is caused by stochastic fluctuations in cluster sample.

Fig. 6. Evolution of the morphological fraction of cluster galaxies for z < 0.6. The black circles with error-bars show the measurements of Fasano et al. (2000), complete down to MV < −20. The thick solid curves show the predictions from our best-fit model using the same photometric completeness, while the thick dashed curves show the predictions using all galaxies available. The thin dashed curves show the assumed field fractions. As discussed in Section 2.2.2, the Ell fraction in the model is set to a constant 27% by construction.

In general, our results agree with those of V11, who found that the EDisCS-to-OmegaWINGS evolution is driven by the combined effect of infall, star formation, and morphological evolution of late-types into early-types. However, our analysis would benefit from tighter observational constrains from novel observations of field and cluster galaxies at various redshift. This will be possible in the near future thanks to the CHileAN Cluster galaxy Evolution Survey (CHANCES; Haines et al. 2023; Sifón et al. 2025), a spectroscopic survey of galaxies in clusters and their surroundings up to z ≃ 0.45 with the 4MOST instrument at the VISTA telescope (de Jong et al. 2019), and to new observational facilities such as the Euclid telescope (Euclid Collaboration: Mellier et al. 2025) or the Vera Rubin observatory (Ivezić et al. 2019).

4. Discussion

4.1. Model limitations

Our models are designed to describe the complex interplay between galaxies and the cluster environment in a simplified (yet flexible) way, and with an affordable computational cost. Apart from the simplified treatment of the environmental mechanisms, perhaps their main limitation is the treatment of galaxy clusters as unique, homogeneous entities, within which galaxies that are accreted from the field (which is also considered homogeneous) evolve according to parametrised recipes. We do not treat effects that are associated with the cluster-centric distance or with the local density, although these are observed to have an impact on the properties of galaxies in clusters (e.g. Pérez-Millán et al. 2023) and in groups (e.g. George et al. 2013). This simplification bypasses detailed treatments for the ICM properties, allowing us to explore a broader parameter space more efficiently. However, it prevents us from investigating the smooth transition of galaxy properties from the field to the cluster environments, which occurs in a transition region that extends beyond the cluster virial radius (e.g. Vulcani et al. 2023).

A related concern is the lack of environment classes aside from clusters and the field, whereas observations and simulations show that 30 − 50% of galaxies reaching clusters today arrive as part of accreted groups rather than individually (e.g. McGee et al. 2009; De Lucia et al. 2012; Benavides et al. 2020). Although we do not discuss specifically the role of pre-processing (e.g. Zabludoff & Mulchaey 1998; Fujita 2004; Sharonova et al. 2025), we stress that the field morphologies adopted at z = 0.055 (Section 2.2.2) are compatible with the actual morphological mix observed around 2 × R₂₀₀ in OmegaWINGS clusters (Fig. 2 of Vulcani et al. 2023). This indicates that, in our model, the properties of the galaxy population infalling onto local clusters are realistic. We discuss further experiments with variations in the field morphology mix in Section 4.4.

Previous studies that modelled the evolution of satellite galaxies in clusters and groups have introduced the so called ‘delayed-than-rapid’ quenching scenario, where the host environment begins to affect galaxies only past a given delay time t_del after their infall (e.g. Wetzel et al. 2013; Oman & Hudson 2016; Fossati et al. 2017). We have verified that the introduction of a delay time parameter in our model - a unique value valid both for the RPS and for the morphology transformation - does not improve the overall quality of the fit, leading to very short t_del (best-fit value of 0.2 Gyr, fiducial value of 0.6 ± 0.4 Gyr) and practically the same likelihood and characteristic timescales of our previous model. A t_del of a few Gyr helps to increase the overall Spi fraction at z = 0.7 due to the delayed transformation, but leads to too many blue galaxies at that redshift due to their delayed quenching. Similarly, it helps to mitigate the excess of low-mass red spirals at z = 0.055 (see Section 4.2), but also produces too many spirals at that redshift. These issues could be alleviated by invoking separate delay times for the morphological and the colour evolutions, or a redshift-dependent t_del. However, exploring these additional parameters goes beyond the purpose of this study.

Finally, our models account only for the ‘negative’ effect of ram pressure on star formation, but both observations (Merluzzi et al. 2013; Vulcani et al. 2018, 2020) and theoretical models (Fujita & Nagashima 1999; Tonnesen & Bryan 2009; Bekki 2014) indicate that, prior to complete gas removal, the star formation activity can be enhanced due to gas compression (e.g. Moretti et al. 2020b,a). Observations suggest that this enhancement can be particularly strong in gas-rich dwarf galaxies (Grishin et al. 2021), and could contribute to populate the sequence of low-M_★ blue spirals which in our model, especially at z = 0.7, is somewhat deficient.

4.2. Variations in the model parameters

To better assess the performance of our approach, we now compare our best-fit model with three variant models where we have modified individual values of the best-fit parameters. In the first variant model we have fixed the stripping timescale to 300 Myr and removed any trend with mass (α_s = 0). In the second, we have increased t_s to an arbitrarily large value, effectively removing the cluster-driven quenching. In the third, we have removed the cluster-driven morphological transformation by imposing p_max = 0.

The results of these experiments are presented in Fig. 7, where we compare the residuals (model-data) in the (B − V)-M_★ plane at the redshift of OmegaWINGS for the four different models. Clearly, the best-fit model provides the smallest residuals and the highest likelihood value. As already mention in Section 3.1, our model appears to slightly underestimate the colour of the S0 population, leading to the offset visible in the bottom-left panel of Fig. 7. It appears that this limitation is intrinsic to our scheme and we cannot amend it with any choice of the model parameters. Redder colours in our model could be obtained by increasing the metallicity of the quenched galaxies, or by using different SPS models. Interestingly, the model with t_s = 300 Myr appears almost identical to our best-fit one, with a single exception: it overestimates the number of red spirals in the low-M_★ regime. The red spiral excess, highlighted with a green ellipse in Fig. 7, is very severe in the lowest bin of M_★, but remains visible up to ∼10^10.3 M_⊙, well above the mass completeness limit in OmegaWINGS. This is a critical feature that pushes the model towards α_s < 0: low-M_★ spirals must have longer stripping timescales (as we have seen in Section 3.2) in order to avoid overpopulating the red sequence. This result holds even when we exclude the lowest M_★ bin from our fit.

Fig. 7. Comparison between the residuals (model-data) in the (B − V)-M★ plane of four different models: our best-fit model (first column), a model with constant ts of 300 Myr (second column), a model without RPS from the ICM (third column), and a model without morphological evolution (fourth column). Red (blue) colours indicate an excess (a deficiency) of galaxies in the model compared to the data. The top and bottom rows show the Spi and S0 populations, respectively. The green ellipse highlights the excess of low-mass red spirals in the model with constant ts. For each model we also list the value of the log likelihood (Equation (6)). For brevity, we only show the results at the redshift of the OmegaWINGS clusters.

The models that exclude either RPS or the morphological transformation provide significantly worse residuals and much lower likelihood values. Excluding RPS leads to overpopulate the blue cloud with both spirals and S0s (and, consequently, to under-populate the red sequence) due to the fact that, while the transformation from Spi to S0 proceeds regularly, galaxies remain on the MS of star formation unless they have joined the cluster as S0 to begin with. Finally, excluding the morphological evolution leads to drastically overestimate the number of spirals, and to underestimate the number of S0s, everywhere in the colour-M_★ plane. We stress that this result holds even if we had fixed the field morphological fractions at all redshifts to the values of Calvi et al. (2012), valid at z ≃ 0, in order to minimise the efficiency of the morphological evolution. Unless the field morphology fractions are significantly different from those adopted here (see Section 4.4), cluster-driven morphological transformation is necessary to explain the data.

4.3. Is RPS as the possible the main driver of the morphological evolution

Although the modelling adopted for the morphological transformation of cluster galaxies described in Section 2.2.5 is inspired by the results of M23, its mathematical implementation is rather general, as the transformation channel is not explicitly specified. The channel can be inferred a posteriori by comparing the resulting transformation timescales with those appropriate for specific channels. However, the right-hand panel of Fig. 5 shows quite some scatter in t_m, with values ranging from less than 1 to several Gyr, which complicates the interpretation of our results.

We show a more detailed comparison with the spectrophotometric evolution results of M23 in Fig. 8. Here, we focus on the evolution of the so called ‘control field’ (CF) galaxies (Vulcani et al. 2018), a subsample of GASP made of 15 field galaxies with typical log₁₀(M_★/M_⊙)≃10.5, which provided the most reliable evolutionary trends due to their regularity and absence of substructures (see Section 3.3 in M23 for additional details). In Fig. 8, we compare the evolution of p_S0 determined for these galaxies by M23 (followed up to 5 Gyr from the onset of the RPS) with that inferred via Equation (4) with our model for log₁₀(M_★/M_⊙) = 10.5 and log₁₀(M_cl/M_⊙) = 14.7. For the CF galaxies of M23 we show the trends for two stripping scenarios, instantaneous and gradual stripping with t_s = 0.5 Gyr, which brackets the t_s range predicted by our model for the M_★ and M_cl considered (left-hand panel of Fig. 5).

Fig. 8. Spiral-to-S0 transformation probability as a function of the time from the galaxy’s entry into the cluster. The blue curves show the results of M23 for the spectrophotometric evolution of a sample of 15 field galaxies with typical log10(M★/M⊙)≃10.5, assuming either instantaneous stripping (solid) or ts = 0.5 Gyr (dashed). The solid black curve shows our fiducial model for log10(M★/M⊙) = 10.5 and log10(Mcl/M⊙) = 14.7, with yellow shaded region showing its scatter given by the difference between the 84th and 16th percentiles. The dotted horizontal line with the grey band on top shows our fiducial pmax and its scatter.

Clearly, the trends predicted by our evolutionary model appears to be consistent with those of the spectrophotometric model for the initial ∼2 Gyr. Beyond that point, the morphology of galaxies ceases to be affected by the ageing of their stellar population: their spectrophotometric evolution is completed, reaching an asymptotic p_max ≃ 0.5. In our evolutionary model, instead, p_S0 keeps on growing with time, flattening to a slightly larger p_max ≃ 0.8. One the one hand, this may point to additional evolutionary channels that keep shaping the morphology of cluster galaxies on longer timescales, such as secular processes (Sellwood & Carlberg 1984; Fujii et al. 2011) or, for the low-M_★ regime, harassment (Moore et al. 1996; Smith et al. 2010a). On the other hand, we stress that the trends resulting from the analysis of M23 are very conservative due to their neglect of stellar kinematics, and in particular of disc rotation, which would smear out the substructures present the galaxy and boost the transition to the S0 type. Additionally, the blue curves in Fig. 8 represent mean trends valid for the whole CF population, but individual galaxies can feature larger p_max. These considerations do not enable us to identify a unique transformation channel, but we can safely state that our results are consistent with the spectrophotometric ageing being a key channel for the spiral-to-S0 transformation observed in cluster galaxies at z ≲ 1.

4.4. Exploring ‘extreme’ field morphologies

One of the main ingredient of our model is the morphological fraction in the field, which we assumed to vary with redshift following the trends shown in the bottom panel of Fig. 2. As these trends are poorly constrained observationally, alternative scenarios deserve to be discussed, at least from a qualitative standpoint.

We consider here two extreme scenarios. In the first one, the field morphological mix does not vary with z and stays constant to the local values found by Calvi et al. (2012). We have refit the data under such conditions, finding that morphological transformation is still required to reproduce the OmegaWINGS data, and with t_m compatible with those of our fiducial model. However, such model largely overestimates the S0 fraction in EDisCS. Trivially, this is caused by the Spi fraction in EDisCS clusters being larger than that in the local field, which itself is an indirect evidence for morphological evolution of the field population.

In the second scenario, we assume that the cluster-driven Spi-to-S0 evolution is suppressed and it is the field morphology that is adjusted to match the cluster data. In practice, in this scenario, we consider pre-processing as the only responsible for the observed morphology mix in clusters. A simple way to investigate this scenario is to ‘reverse engineering’ our best-fit model: at all redshifts, we determine the field fractions that would be required to match the morphological mix of the best-fit cluster model in the absence of internal evolution. The fractions inferred in this way are shown in Fig. 9. Clearly, matching the cluster data requires a field dominated by S0s for z < 1, and with an S0 fraction that steadily increases at later times, in strike contrast with the morphologies observed in the local field and at the periphery of OmegaWINGS clusters (Vulcani et al. 2023). Paradoxically, this approach demands negative spiral fractions at later times, caused by the fact that, in the model, the total number of spiral cluster galaxies decrease with time as the morphological transformation proceeds at a rate faster than that at which new spirals are injected from the field. These findings exclude preprocessing as the main responsible for the morphology mix observed in local clusters, corroborates our original choice for the field morphology fractions, and strengthen the overall results of our analysis.

Fig. 9. Comparison between the field morphological fractions used in our main analysis (dashed curves) and those inferred for an extreme model where the field is fully responsible for the morphological evolution of cluster galaxies (solid line). In both cases the Ell fraction is fixed to 0.27 at all redshifts.

4.5. Stripping and orbit anisotropy

In Section 3.2 we have shown that, in our models, t_s and M_★/M_cl are anti-correlated: for a fixed cluster mass, the stripping proceeds more efficiently in high-M_★ galaxies, which is unexpected. One possible interpretation of this puzzling result is that galaxies with different M_★ join the cluster following substantially different trajectories, with the least (most) massive systems having a stronger tangential (radial) component in their initial velocity.

Strongly radial orbits lead galaxies closer to the pericenter and with a larger speed, maximising the RPS efficiency. We notice that a high RP does not necessarily imply a large occurrence of jellyfish systems, given that low-mass galaxies (which are the most numerous) will be stripped on very short timescales, limiting the possibility to observe the mechanism in action (e.g. Marasco et al. 2016). Observations of jellyfish galaxies in different mass regimes appear to support this, indicating preferentially radial (tangential) orbits for high- (low-) mass jellyfishes (Grishin et al. 2021; Biviano et al. 2024).

Several studies have investigated the orbit anisotropy in galaxy clusters as a function of galaxy mass, morphology and cluster-centric distance, using either dynamical modelling of the observed cluster kinematics (e.g. Adami et al. 1998; Annunziatella et al. 2016; Biviano et al. 2021) or large-scale cosmological simulations (e.g. Lotz et al. 2019). The general picture emerging from these studies is that the anisotropy parameter β declines from positive values (indicating preferentially radial orbits) in the cluster outskirts, where star-forming galaxies reside, to negative values (tangential orbits) in the cluster core, which is dominated by passive, massive systems. This supports a scenario where low-mass, star-forming systems joining the cluster from the field with β > 0 are tidally destroyed, biasing the surviving population within the core towards β < 0.

While the results discussed above are valid for the cluster population as a whole, our findings concern the infalling population alone. The existence of a correlation between M_★ and β for the galaxies that join the cluster from the field is yet to be proven, and it may be challenging to probe it observationally. Instead, it should be straightforward to investigate it using cosmological hydrodynamical simulations. However, this goes beyond the purpose of our study.

5. Conclusions

Galaxy clusters are formidable laboratories for studying the effect of the environment on the evolution of their members. Compelling observational evidence indicate a rapid evolution in the properties of the disc galaxy population in clusters over the past ≃7 Gyr, with passive, lenticular (S0) systems dominating in local clusters, while more actively star forming spirals (Spi) dominate the population at z ≃ 0.7 and beyond. This evolution does not match that observed in the field, implying an active role played by the cluster environment in driving the observed trends (e.g. Guglielmo et al. 2015).

In this study, we investigated whether the observed evolution can be explained by simple cluster evolutionary models in the ΛCDM framework, where the main mechanisms driving the evolution are ram pressure stripping (RPS) from the intra-cluster medium and a gradual transformation of cluster spirals into S0s due to an (initially) unspecified process. We assume that the main model parameters, the stripping timescale, t_s, and the Spi-to-S0 transformation timescale, t_m, are independent power-law functions of the galaxy-to-cluster mass ratio, M_★/M_cl, and constrain their values using the observed spiral and S0 distributions in the colour-M_★ plane from the EDisCS survey at z ≃ 0.7 and the OmegaWINGS survey at z ≃ 0.055. Our results are summarised below.

• Our best-fit model reproduces the observed distribution of spirals and S0 cluster galaxies in the colour-M_★ plane very well, especially at low redshift (Fig. 4), where the model has vastly lost memory of its initial conditions. However, compared to the data, our model slightly overestimates the S0 fraction at z = 0.7, and underestimate the (B − V) colour of the S0 population at z = 0.055 by ∼0.1 dex (bottom-left panel of Fig. 7).

• We find typical t_s values between 100 Myr and 1 Gyr, in agreement with previous estimates (Tonnesen & Bryan 2009; Jaffé et al. 2018; Salinas et al. 2024). Typical t_m values are instead of a few Gyr, albeit with a large scatter (Fig. 5). We interpret t_m > t_s in terms of cause and effect: the morphological transformation follows the quenching of star formation caused by RPS.

• The morphological evolution inferred from our model is compatible with that resulting from the simple ageing of the stellar populations in outside-in quenched galaxies, as proposed by M23. Although we cannot exclude the fact that secular evolution plays a role in shaping galaxy morphology on longer timescales, our results confirm spectrophotometric ageing as a key channel for the Spi-to-S0 transition in galaxy clusters.

An unexpected and somewhat controversial finding is the strong anti-correlation between t_s and M_★/M_cl, driven by the relatively blue colour of low-M_★ spirals in the OmegaWINGS clusters: assuming a constant t_s at all masses leads to overestimate the number of low-M_★ red spirals compared to the data (Fig. 7, second panel in the top row). One possible interpretation is provided in terms of orbit anisotropy, discussed in Section 4.5, but a detailed comparison with large-scale cosmological simulations is required to validate this scenario.

Our analysis is based on detailed information on the cluster galaxy populations available in two specific redshift bins. We expect forthcoming observational datasets from the CHANCES spectroscopic survey and from new facilities, such as the Euclid telescope and the Vera Rubin Observatory, to provide tighter constraints to our model over a broader redshift range, significantly boosting our understanding of the role played by dense environments in the formation and evolution of galaxies.

Acknowledgments

The authors thank an anonymous referee for a thoughtful and constructive report. AM acknowledges funding from the INAF Mini Grant 2023 program ‘The quest for gas accretion: modelling the dynamics of extra-planar gas in nearby galaxies’. BV acknowledges support from the INAF GO grant 2023 “Identifying ram pressure induced unwinding arms in cluster spirals” (P.I. Vulcani).

References

Appendix A: Relation between star formation history and main sequence

In this appendix, we provide details on the SFH models adopted for our main-sequence field galaxies at the various redshifts. Following Gladders et al. (2013), we use log-normal SFH models that are described as

$$\begin{array}{c}\hfill \mathrm{SFR}(t,t_0,\tau )=\frac{1}{t\sqrt{2\pi {\tau }^2}}exp[-\frac{{(lnt-t_0)}^2}{2{\tau }^2}],\end{array}$$(A.1)

where t is the elapsed time since the Big Bang, t₀ is the logarithmic time delay, and τ sets the rise and decay timescale.

Our goal is to find optimised relations between the two parameters of Equation (A.1), t₀ and τ, and the M_★ of a z = 0 galaxy, so that such galaxy stays (approximately) on the MS of Popesso et al. (2023) up to z ≃ 2. This is achieved numerically with the following approach. First, we build an initial, large (∼2 × 10⁵) family of model SFHs with varying values of t₀, τ, and M_f, 0, the total gas mass converted into stars by z = 0 given by the integral of Equation (A.1). For each model we store also total gas mass converted into stars at any z, M_f, z. Next, we convert our M_f, z to M_★ estimates, which exclude the gas lost by stellar winds and supernovae (but include stellar remnants), using the returned gas fractions provided by the Bruzual & Charlot (2003) SPS model (2016 version) available in BAGPIPES. These depend mainly on the galaxy SFH and only to a minor extent on the stellar metallicity, which is assumed to be Solar for the purpose of this calculation.

We then exclude from our family all models that provide M_★ and SFRs incompatible with the MS. This is done by evaluating the typical ‘distance’ between our models and the MS at four representative redshifts (0, 0.5, 1, 2), and then excluding models with distance larger than a fixed threshold. This filters out the vast majority of models, leaving only those that are compatible with the MS. As we show in panels (a) and (b) of Fig (A.1), the surviving models show clear correlations between M_★ and the parameters of Equation (A.1). We regularise these correlation using fifth-degree polynomials, which optimally describe the observed trends. Finally, the polynomial fits are used to build our fiducial SFH library for MS galaxies, which we utilise in our cluster evolutionary model.

Fig. A.1. Panels (a) and (b): Relation between τ (a) and t0 (b), the two parameters of our SFH model (Equation A.1), and the galaxy M★ at z = 0, constrained by the Popesso et al. (2023) main sequence (MS) of star formation. The blue dots show the models in our library that minimise the distance to the MS at various z. The red-dashed curves show our regularisation using fifth degree polynomials. Panel (c): SFR-M★ relations at four different redshifts. The solid curves show the MS of Popesso et al. (2023), the filled circles show the sequence emerging from our fiducial SFH library. Panel (d): Our fiducial SFH library adopted for MS galaxies, colour-coded by the galaxy M★ at z = 0. Our models support the downsizing scenario of M★ build-up.

Panel (c) of Fig. (A.1) shows the comparison between the MS of Popesso et al. (2023) and those predicted by our fiducial SFH library. Clearly, these can reproduce very well the ‘bending’ of the MS visible at low z while still providing a good match to the high-z MS, although their accuracy gets somewhat worse around z ≃ 2. The resulting sequence of SFHs, shown in panel (d) of Fig. (A.1), presents clean distinction between low- and high-M_★ galaxies, with the former showing gently rising SFHs, while the latter feature a sharp peak in their SFR at earlier epochs, followed by a gradual decline. These trends agree with the so-called ‘downsizing’ scenario (e.g. Cowie et al. 1996), where more massive galaxies have assembled their M_★ faster and at earlier epochs than the low-mass ones.

The same technique adopted here can be applied to any MS of star formation, provided that its evolution with z is known. We applied our method to the MS of Speagle et al. (2014), finding analogous correlations between M_★, t₀, and τ, and qualitatively similar SFHs supporting the downsizing scenario.

Appendix B: Supplementary material

Fig. B.1 shows the correlation between the various parameters for the model investigated in the main text. The posteriors are all well behaved, but some degeneracy is present between τ_m and p_max. Table B.1 is equivalent to Table 1, but is valid for models that use the MS of Speagle et al. (2014) instead of Popesso et al. (2023). The two MSs provide perfectly compatible fiducial parameters.

Fig. B.1. Corner plots showing the correlation between the various parameters (shaded regions, with contours at arbitrary iso-density levels) along with their marginalised probability distribution (histograms on top) for our model using the MS of Popesso et al. (2019).

All Tables

All Figures

	Fig. 1. Flowchart illustrating our method. White blocks show pre-existing dataset or software packages. Yellow blocks show the new products built for this study. The chart indicates also the Sections of the main text that provide details on a given part of the method.
In the text

Fig. 2. Evolutionary properties of our model cluster. Top panel: Mass (solid black curve) and accretion rate (dashed grey curve) vs redshift for a cluster with Mcl = 5 × 1014 M⊙ at z = 0.055. The blue square and the red circle show the median Mcl and redshifts of the EDisCS and OmegaWINGS clusters, respectively, with the error bars corresponding to half the difference between the 84th and 16th percentiles of the galaxy distribution. Bottom panel: Morphology fractions for field galaxies that join the cluster as a function of z, assumed to be independent of the galaxy M★.

In the text

Fig. 3. Photometric evolution and mass build-up of an idealised main-sequence galaxy with initial M★ ≃ 2 × 1010 M⊙ that joins our cluster model at z = 1, for stripping timescales of 0.5 and 2.0 Gyr (dark-red and orange curves, respectively), and for a case without stripping (blue curve). Top panel: SFR as a function of the look-back time. For times larger than ∼8 Gyr (z > 1), the three models share a common SFH shown by the solid black curve. Bottom panel: Galaxy evolution in the rest-frame (B − V) colour vs M★ plane. The white square shows the initial values at z = 1, the coloured circles show the final (z = 0) values for the three scenarios considered. The green dashed-line separates the red sequence from the blue cloud (from Vulcani et al. 2022).

In the text

Fig. 4. Comparison of our best-fit cluster model with the EDisCS data at z = 0.7 (top six panels), and with the OmegaWINGS data at z = 0.055 (bottom six panels). The data and the model are shown in the first and second row of each panel set, respectively. Panels in the leftmost (central) column show the distribution of spirals (S0s) in a colour-M★ plane, where the colour is the observed-frame (V − I) for EDisCS and the (B − V) for OmegaWINGS. The colour palette indicate the number of galaxies per bin. Cyan contours are drawn at 2, 5, 10, and 15 galaxies per bin in EDisCS and in the model at z = 0.7, and at 2, 10, 30, and 100 galaxies per bin in OmegaWINGS and in the model at z = 0.055. The total fraction of Spi and S0 galaxies with respect to the total disc (Spi+S0) population is indicated in the top-right corner of each panel. Panels in the rightmost column show the S0 fraction as a function of colour and M★, for bins that have at least five galaxies in the data.

In the text

Fig. 5. Characteristic timescales for the ram-pressure stripping (blue) and for the cluster-driven morphological transformation (red) as a function of the galaxy-to-cluster mass ratio (M★/Mcl) derived with our analysis (Equations 3 and 5). The solid lines show the median trend, with the shaded regions showing the difference between the 16th and 84th percentiles. The dashed lines show the best-fit model. As a reference, in the horizontal axis on top we report the M★ scale for a fixed log10(Mcl/ M⊙) = 14.7.

In the text

Fig. 6. Evolution of the morphological fraction of cluster galaxies for z < 0.6. The black circles with error-bars show the measurements of Fasano et al. (2000), complete down to MV < −20. The thick solid curves show the predictions from our best-fit model using the same photometric completeness, while the thick dashed curves show the predictions using all galaxies available. The thin dashed curves show the assumed field fractions. As discussed in Section 2.2.2, the Ell fraction in the model is set to a constant 27% by construction.

In the text

Fig. 7. Comparison between the residuals (model-data) in the (B − V)-M★ plane of four different models: our best-fit model (first column), a model with constant ts of 300 Myr (second column), a model without RPS from the ICM (third column), and a model without morphological evolution (fourth column). Red (blue) colours indicate an excess (a deficiency) of galaxies in the model compared to the data. The top and bottom rows show the Spi and S0 populations, respectively. The green ellipse highlights the excess of low-mass red spirals in the model with constant ts. For each model we also list the value of the log likelihood (Equation (6)). For brevity, we only show the results at the redshift of the OmegaWINGS clusters.

In the text

Fig. 8. Spiral-to-S0 transformation probability as a function of the time from the galaxy’s entry into the cluster. The blue curves show the results of M23 for the spectrophotometric evolution of a sample of 15 field galaxies with typical log10(M★/M⊙)≃10.5, assuming either instantaneous stripping (solid) or ts = 0.5 Gyr (dashed). The solid black curve shows our fiducial model for log10(M★/M⊙) = 10.5 and log10(Mcl/M⊙) = 14.7, with yellow shaded region showing its scatter given by the difference between the 84th and 16th percentiles. The dotted horizontal line with the grey band on top shows our fiducial pmax and its scatter.

In the text

	Fig. 9. Comparison between the field morphological fractions used in our main analysis (dashed curves) and those inferred for an extreme model where the field is fully responsible for the morphological evolution of cluster galaxies (solid line). In both cases the Ell fraction is fixed to 0.27 at all redshifts.
In the text

Fig. A.1. Panels (a) and (b): Relation between τ (a) and t0 (b), the two parameters of our SFH model (Equation A.1), and the galaxy M★ at z = 0, constrained by the Popesso et al. (2023) main sequence (MS) of star formation. The blue dots show the models in our library that minimise the distance to the MS at various z. The red-dashed curves show our regularisation using fifth degree polynomials. Panel (c): SFR-M★ relations at four different redshifts. The solid curves show the MS of Popesso et al. (2023), the filled circles show the sequence emerging from our fiducial SFH library. Panel (d): Our fiducial SFH library adopted for MS galaxies, colour-coded by the galaxy M★ at z = 0. Our models support the downsizing scenario of M★ build-up.

In the text

	Fig. B.1. Corner plots showing the correlation between the various parameters (shaded regions, with contours at arbitrary iso-density levels) along with their marginalised probability distribution (histograms on top) for our model using the MS of Popesso et al. (2019).
In the text