Revealing impacts of stellar mass and environment on galaxy quenching

Aims. Galaxy quenching is a critical step in galaxy evolution. In this work, we present a statistical study of galaxy quenching in 17 cluster candidates at 0.5 < z < 1.0 in the COSMOS ﬁeld. Methods. We selected cluster members with a wide range of stellar mass and environment to study their mass and environment dependence. Member galaxies are classiﬁed into star-forming, quiescent and recently-quenched galaxies (RQG) using the rest-frame UVJ diagram. We further separated fast and slow quenching RQGs by model evolutionary tracks on the UVJ diagram. We deﬁned the quenching e ﬃ ciency as the ratio of RQGs over star-forming galaxies and the quenching stage as the ratio of RQGs over quiescent galaxies to quantify the quenching processes. Results. We found quenching e ﬃ ciency is enhanced by both higher stellar mass and denser environment. Massive or dense environment galaxies quench earlier. Slow quenching is more dominant for massive galaxies and at lower redshifts, but no clear dependence on the environment is found. Our results suggest that low-mass galaxies in dense environments are likely quenched through a short-timescale process such as ram pressure stripping, while massive galaxies in a sparse environment are mostly quenched by a longer-timescale process. Using the line strength of H δ and [OII], we conﬁrmed that our UVJ method to select RQGs agrees with high S / N DEIMOS spectra. However, we caution that the visibility time (duration of a galaxy’s stay in the RQG region on the UVJ diagram) may also depend on mass or environment. The method introduced in this work can be applied to RQG candidates for future statistical RQG spectroscopic surveys. The systematic spectroscopic RQG study will disentangle the degeneracy between visibility time and quenching properties.


Introduction
Deep surveys have transformed our view of the distant Universe, allowing us to constrain the process of galaxy formation and evolution over the past 13 billion years of cosmic history. The cosmic star formation rate density comes to a peak at z ∼ 2 − 3, and then it declines as the Universe ages until the present day (Lilly et al. 1996;Karim et al. 2011;Burgarella et al. 2013;Madau & Dickinson 2014;Khostovan et al. 2015). This decline in star formation is called galaxy quenching, an essential step in galaxy evolution. Although there have been many efforts to understand quenching mechanisms, there is still a lack of conclusive results.
Quenching is a process where star-forming galaxies (SFGs) transform into quiescent galaxies (QGs). The transitional population between SFGs and QGs are recently quenched galaxies (RQGs), which are ideal for studying the quenching process . The RQGs are still in the quenching process or have just finished (Newberry et al. 1990;Wild et al. 2009;French et al. 2018). This nature makes them capable of delivering information about the physical processes during quenching. RQGs have quenched their star formation in the past few ten million years to a few billion years (Couch & Sharples 1987;Poggianti et al. 2009), but they still contain a large fraction of relatively young A-type stars. This particular stellar composition adds strong Balmer absorption lines to the passive galaxy spectrum (Dressler & Gunn 1983;Couch & Sharples 1987;Wild et al. 2009). In previous studies, this population is mostly referred to as post-starburst galaxies (PSBs), which means galaxies that have experienced a starburst episode followed by a fast quenching. However, this starburst phase is not essential for quenching, and the quenching timescale varies for each galaxy (Poggianti et al. 2009). Therefore, we use recently quenched galaxy in this work as a more precise name for this population.
There are two key parameters for describing galaxy quenching, mass and environment, as the red fraction of galaxies is known to correlate strongly with these two parameters (Peng et al. 2010;Sobral et al. 2011;Muzzin et al. 2012;Darvish et al. 2016). The physical mechanisms of quenching are correspondingly categorised into mass-dependent and environmentdependent types.
Among the stellar mass-driven quenching, the most studied mechanism is active galactic nucleus (AGN) feedback (Best et al. 2005;Hopkins et al. 2005;Kaviraj et al. 2007;Diamond-Stanic et al. 2012;Belli et al. 2021), which is also studied in simulations (e.g. Dubois et al. 2013;Beckmann et al. 2017;Weinberger et al. 2017;Piotrowska et al. 2022;Shi et al. 2022;Xu Article number, page 1 of 20 Article published by EDP Sciences, to be cited as https://doi.org/10.1051/0004-6361/202243733 A&A proofs: manuscript no. output et al. 2022). The quasar host AGNs show evidence of having strong outflows in observation (Dai et al. 2008;Tombesi et al. 2010). AGN outflows cause expulsion of the gas from the gravitational bound of the system, and hence reduce the subsequent star formation activity in their host galaxies (Debuhr et al. 2012;Combes 2017). AGNs can also heat up the cool gas of the interstellar medium, hence preventing the galaxy from forming new stars (Croton et al. 2006;Man & Belli 2018;Zinger et al. 2020). Supernovae feedback is another primary mechanism of mass-dependent quenching that mainly acts on low-mass galaxies. High-energy ejections and superwinds take away the cold gas, stopping further star formation activity (Dekel & Silk 1986;Dekel & Woo 2003).
On the other hand, quenching mechanisms related to interactions between galaxies and their environment, for example other galaxies, the intracluster medium (ICM), and the intergalactic medium (IGM), are considered environmental quenching. A widely established mechanism of environmental quenching is strangulation, also known as starvation (Larson et al. 1980). The virial shocks from the halo heat up the cool gas, ceasing the galaxies' gas supply. The star formation will then be quenched due to the lack of fuel (Bekki et al. 2002;Man & Belli 2018). Galaxy mergers can be another mechanism of environmental quenching. In particular, major mergers can trigger starbursts, which are followed by rapid consumption of all remaining cool gas (Mihos & Hernquist 1996;Gabor et al. 2010;Man & Belli 2018). In addition, when galaxies encounter each other at high speed, their potential wells will interact, causing gravitational shocks. Even though the effect of individual encounters may not be large, the repeated process in high-density environments would make a significant impact on galaxies. This process is known as galaxy harassment (Moore et al. 1996(Moore et al. , 1998. Galaxy harassment leads to instability of galaxy morphology and change in the distribution of material, driving gas into central regions (Moore et al. 1999;Boselli & Gavazzi 2006). Galaxies will be quenched after a starburst consumes their gas. Ram pressure stripping is also an effective environmental quenching mechanism (Ma et al. 2008;Bekki 2009;Roediger et al. 2014). As galaxies move through the ICM in cluster cores, the wind from the ICM will strip the gas from them (Gunn & Gott 1972;Balogh et al. 2000;Steinhauser et al. 2016).
Many mechanisms that probably lead to quenching have been proposed, and surely we need some parameters to discriminate between these mechanisms. One of such critical parameters we are looking for is the quenching timescale, which varies from one mechanism to another. Most of these timescales are obtained from simulations (Wright et al. 2019;Wetzel et al. 2013;Walters et al. 2022). Among mass-dependent quenching mechanisms, the supernova feedback is a strong physical process that quickly quenches star formation in ∼0.1 Gyr (Ceverino & Klypin 2009). On the other hand, the AGN feedback timescale is still under debate. Quasar mode AGN feedback is a strong process. Quasars' lifetimes are not long (<0.1 Gyr, Hopkins et al. 2006), and the gas outflow driven by quasars can quench star formation at million-year-level timescales (Smethurst et al. 2021). On the other hand, the radio mode of AGN feedback is a slow process (Best et al. 2005); it takes up to and beyond 1 Gyr to establish a balance between cooling and heating to reach a low star formation rate phase (Fabian 2012). In Schawinski et al. (2014) the authors maintain that the timescale of AGN feedback leading to a quenching can be related to the host galaxy type. Furthermore, Hirschmann et al. (2017) simulated the star formation history (SFH) of galaxies with AGN and found that the quenching timescale can have a wide range, from a few hundred Myr to a few Gyr. For environmental quenching mechanisms, strangulation is a long-term process that lasts a few billion years (Peng et al. 2015). Merger-driven quenching does not happen directly after galaxy merger events; a median delay time of 1.5 Gyr is expected, and the timescale varies over a wide range (Rodríguez Montero et al. 2019). Ram pressure stripping is a rapid type of quenching, with a timescale of ∼0.2 Gyr (Steinhauser et al. 2016). Thus, we can use the quenching timescale to determine some of the physical mechanisms involved.
To study the quenching process, the main drawback of current RQG works is the lack of observations. Due to their intrinsic short-lived nature, RQGs are rare objects (Tran et al. 2004;Le Borgne et al. 2006;Kaviraj et al. 2007), and we still lack a large spectroscopic sample to carry out systematic studies. Although RQG spectra can provide direct information and strict constraints on their physical properties and SFH, only a limited number of RQG spectra exist at high redshifts (e.g. Pracy et al. 2010;Maltby et al. 2016;D'Eugenio et al. 2020). Increasing the spectroscopic sample is critical, and for this purpose we first need to make a large photometric sample of RQG candidates to feed spectroscopic follow-up observations. Moreover, such a statistical photometrically selected sample of RQG candidates is essential for us to explore the quenching processes prior to spectroscopic studies. As the quenching process is strongly correlated with mass and environment, a systematic study of RQGs must cover a wide range of masses and environments. Therefore, we require a sample of RQGs extending from cluster cores to general fields in order to understand environmental quenching. There are some previous works that study the environment-dependent quenching using the post-starburst galaxy (PSB) population, which belongs to the RQG population (e.g. Socolovsky et al. 2018;Paccagnella et al. 2017Paccagnella et al. , 2019. They study the environment-dependent quenching by estimating the PSB fraction or number in different environments (e.g. cluster, group, or field; distance to cluster centre). These studies prove that the PSBs tend to reside in denser environments. Their results also suggest that there should be various quenching mechanisms in different environments (Paccagnella et al. 2017). In clusters, fast quenching processes, such as ram pressure stripping, are enhanced. In sparse environments, galaxy mergers and interactions produce PSBs (Socolovsky et al. 2018;Paccagnella et al. 2019). There are also works on the mass dependence of galaxy quenching. McNab et al. (2021) studied the mass dependence of quenching for transitional populations (similar to RQGs) using spectroscopic data from the GOGREEN survey (Balogh et al. 2017). Their results support the scenario of preprocessing in the group or proto-cluster environment for massive galaxies and the fast quenching during infall for lower mass galaxies. The GOGREEN survey contains spectra of many cluster galaxies with secure redshift at 0.8 < z < 1.5, which will improve the study of transitional populations in the future.
In this work we present a systematic study of RQGs based on multi-band photometric data in 17 cluster candidates in the COSMOS field at 0.5 < z < 1.0. Our aim is to develop an efficient method of selecting RQGs from photometric surveys (e.g. HSC-SSP, Aihara et al. 2019) to prepare for the coming spectroscopic RQG survey using the Prime Focus Spectrograph (PFS, Tamura et al. 2018) on Subaru. We selected our RQG sample with the rest-frame UVJ colour-colour diagram and confirmed some of the proposed quenching mechanisms using this sample. With the current photometrically selected RQGs, we conducted basic analyses on mass and environment dependence of galaxy quenching and obtained some preliminary results. This paper is structured as follows. In Sect. 2 we present the photometric data we used and the cluster galaxy sample we selected. In Sect. 3 we present the method we use to classify RQGs and correct for possible projection effects. In Sect. 4 we display our results of quenching efficiency, quenching stage, and quenching timescale. In Sect. 5 we discuss the possible contamination and incompleteness of our RQG selection, compare our results with previous works, and then attempt to interpret the results we obtained from the current sample. An H 0 = 69.6kms −1 Mpc, Ω M = 0.286 and Ω Λ = 0.714 cosmology is used. All cosmology calculations are conducted by the cosmology calculator from Wright (2006). All magnitudes presented in this paper are defined in the AB magnitude system.

Dataset and sample
In Sect. 2.1 we describe the photometric data we use in this study. In Sect. 2.2 we describe how we selected our sample galaxies to study the quenching processes.

The COSMOS2015 catalogue
The photometric data that we use in this work is part of the COSMOS2015 catalogue (Laigle et al. 2016). This catalogue mainly combines 26 bands from optical to infrared wavelengths of galaxies in the COSMOS field. It offers deep (m AB ∼25-26 mag) data covering a relatively large area of 2 deg 2 . As we require near-infrared (NIR) photometry to estimate the rest-frame J-band magnitude, the data we use is limited to the central ∼ 1.4deg 2 of the COSMOS field. In this work we use B, V, r, i + , z ++ taken with Suprime-Cam/Subaru; Y, J, H, K s bands taken with VIRCAM/VISTA; and 3.6µm, 4.5µm taken with IRAC/Spitzer. Broad-band B, V, r, i + data are from the previous data release (Capak et al. 2007;Ilbert et al. 2009). Broad z ++ band data are taken with new upgraded CCDs. For each optical band, image PSFs are homogenised so that tile-to-tile variations can reach minimisation (Capak et al. 2007). NIR Y, J, H, K s band data are part of the UltraVISTA survey, and the data combined with the COSMOS2015 catalogue was released in UltraVISTA-DR2, which is much deeper in ultra-deep stripes compared to DR1 (McCracken et al. 2012). The IRAC 3.6µm, 4.5µm data are a combination of several programmes observing the COSMOS field using Spitzer (Capak 2016). We corrected for foreground galactic extinction according to the catalogue included reddening value EBV (=E(B−V)) for each object and foreground extinction factor F f (Allen 1976) according to Laigle et al. (2016): In this work we use photometric redshift and stellar mass estimation results from Laigle et al. (2016). In that work spectral energy distribution (SED) fitting and photometric redshift estimation are performed with the LEPHARE code (Arnouts et al. 2002;Ilbert et al. 2006). They use a set of 31 templates, including spiral and elliptical galaxies from Polletta et al. (2007) together with 12 templates of young blue SFGs generated by the Bruzual & Charlot (2003) models. Dust attenuation is added as a free parameter; 3" aperture magnitudes are adopted since they provide slightly better photometric redshift. The code performs χ 2 analysis between template predicted and observed fluxes. Stellar masses are then derived by the method of Ilbert et al. (2015), using a spectral library generated by the Bruzual & Charlot (2003) model. The Chabrier (2003) initial mass function is assumed.
A 97% confidence-level spectroscopic redshift sample is obtained from Lilly et al. (2007) to examine the reliability of photometric redshift. The precision of photometric redshift is estimated using normalised median absolute deviation (NMAD, Hoaglin et al. 1983), which is defined as 1.48 × median(|z p − z s |/(1 + z s )), where z p refers to photometric redshift and z s refers to spectroscopic redshift. The dispersion of NMAD of the spectroscopically confirmed sample is σ=0.021. Since the width of the redshift slices we used in this study is 0.1, this sigma level is precise enough for our further study. This work focuses on the redshift range z ∼ 0.5 − 1.0, where we have more confirmation of photometric redshift and higher mass completeness (see Sect. 2.2). Moreover, in this range, the UVJ selection criteria are not affected by redshift; hence, we can apply the same UVJ selection region for the whole sample (more information in Appendix A).

Sample selection
In this section we describe how we select galaxy cluster candidates and their member galaxies.
We followed the steps in Laigle et al. (2016) to obtain mass completeness for this sample. For each galaxy the mass required for it to be observed (M lim ) is determined based on K s magnitude and stellar mass M according to the following equation: (2) For each galaxy, the 5σ limiting magnitude (K s,lim ) is used to estimate the stellar mass required for the galaxy to be observed (M lim ). Then in each redshift bin the mass completeness is set to the M lim value that contains 90% of the galaxies. The qualitative conclusions we draw remain unchanged when we apply a slightly higher mass cut, such as 95%. In the overall sample, we remove galaxies without a stellar mass estimation and adopt the most conservative mass limit in the highest redshift bin at M lim = 10 9.8 M . All further analyses are carried out for galaxies with stellar masses higher than the mass limit based on the SED fitting.
We then select galaxies with type flag=0 (galaxy type). Since we use five optical bands from Subaru and four UVista nearinfrared bands to calculate rest-frame UVJ colour, we require our target galaxy to cover at least two optical and two NIR bands to ensure the accuracy of rest-frame colour estimation. After these selections, we have 31983 galaxies. We further select our analysis sample and field sample from these galaxies.
In order to study the environmental dependence of galaxy properties in an unbiased way, we select our sample galaxies from the cluster cores to the outskirts. Thus, our first step is to search for galaxy cluster candidates at 0.5 < z < 1.

Overdensity identification
In this study we select overdensity regions based on the aperture density method, computing the aperture density (i.e. the number density of galaxies in a certain aperture) across our field using 1 Mpc radius apertures. Since there is no complete cluster catalogue at this redshift range, we choose to identify our own cluster candidates to remain consistent throughout this work in the selection and analysis of our cluster and field samples.
We first slice the overall sample into five redshift bins with a width of ∆z = 0.1 each. Then we calculate the aperture density distribution for each slice and make the aperture density maps (e.g. Fig. 1) by placing 1 Mpc radius apertures on 1 Mpc × 1 Mpc grids that cover the whole COSMOS field. We adopt this  (Chiang et al. 2013). Then we define σ as the scatter of densities in all apertures. Regions with > 3σ overdensity (i.e. the density is at least 3σ higher than the average density) are selected as possible cluster candidates. In some cases, the redshift distribution of cluster members is too close to the edge of the redshift slice under scrutiny. Thus, we shift the redshift bin by z∼0.05 to include the majority of the members (see more in Sect. 2.2.3). We selected 37 overdensity regions in all redshift bins, 8 of which are the same overdensity regions at different redshifts. As a result, we found 29 different overdensity regions in 0.5 < z < 1.0.

Member galaxy selection
We adopted 1 Mpc radius apertures to select overdensities, but they cannot trace fine structures and local environments within cluster candidates. Thus, we used the local density, Σ 5 , to probe the local environment; Σ 5 is defined as the average density within a circle centred at a given position with a radius of the distance between the central point and its fifth nearest neighbour. We note that we should not use average Σ 5 as the overall field average density due to the uncertainty of photometric redshift. In redshift confirmed overdensities (see Sect. 2.2.3), the nearest neighbours have lower probabilities of being foreground or background galaxies. However, in low-density field areas, the impact of projection effects is strong and cannot be corrected.
To determine the centre of the cluster candidates, the Σ 5 values are calculated on 0.01 Mpc steps covering the 1 Mpc × 1 Mpc square region centring the peak of aperture density in each cluster candidate. The Σ 5 peaks within the mentioned square region are determined as the centres. The information of these cluster candidate centres is displayed in Table 1. We then put a square aperture at the cluster candidate's centres and selected all galaxies in this redshift slice to be mem- Fig. 1. Aperture density map of a redshift slice at z∼0.5-0.6. The colours show aperture density in the logarithm scale. The red contours indicate the edge of 3σ peaks, and the blue contours the edge of 1σ regions. The area outside the 1σ contours is identified as field. This kind of density map is plotted for all redshift slices to select cluster candidates and construct the field sample, which is used to make a statistical subtraction of the field average density. Stars indicate the confirmed cluster candidates analysed in this work. bers of this cluster candidate. The actual cluster size depends on mass, redshift, and cluster model. Here we apply fixed 2.8 Mpc × 2.8 Mpc apertures for simplicity because most of our clusters do not have a measured virial radius. We choose a relatively large aperture size to capture more members of clusters, especially massive clusters. The aperture size we choose is still consistent with the typical cluster radius ∼ 1.1 -1.6 Mpc at z ∼ 0.8 (Chiang et al. 2013). The field galaxies will not be a huge contamination because we checked the redshift distribution of cluster members (see Sect. 2.2.3) and subtracted the field galaxies to correct for the projection effect (see Sect. 3.3). Since we discuss the environmental effect with local density Σ 5 , the square aperture selection (instead of careful cluster member identification) is good enough for our analysis. Small changes in the dimensions (2.5 -3 Mpc) and the shape of the aperture do not qualitatively change our results. We checked the overlap between cluster candidate regions. There are two cases: cluster candidates 02 and 03 have an overlap of 2.3 Mpc 2 and cluster candidates 14 and 15 have an overlap of 0.56 Mpc 2 . We did not apply any special treatment for the overlapping regions; galaxies in that region are calculated twice as a member of each cluster candidate. Considering that the overlapping region is small, it does not affect our results.

Redshift distribution confirmation
For the overdensity estimation we included all galaxies within ∆z = 0.1, which is larger than the real cluster redshift range. However, the photometric redshift has some uncertainty, and we do not have kinematic information on galaxies. Therefore, it could be possible that some of these overdensities are created by projection effects. Thus, we must confirm the distribution of the potential cluster candidate members in the redshift space.
A real cluster should have Gaussian-like photometric redshift distribution (top panel of Fig. 2), while fake clusters, due to the projection effect, have a flat or a multi-peak distribution Flat distribution, indicating that this overdensity is not a physically connected structure. Bottom panel: Distribution of a real structure, but its peak is close to the edge of our selection redshift slice, missing a large fraction of its galaxy members.
(middle panel of Fig. 2). To remove the contamination from the fake clusters, we calculated the scatter σ z of the photometric redshift distribution of cluster candidate members. A scatter value σ z < 0.025 is required for an overdensity to be selected as a cluster candidate. In some other cases, the redshift distribution shows one side of Gaussian with a peak towards the edge of the redshift bin (bottom panel of Fig. 2). This type of redshift distribution may belong to a real cluster at a different redshift. When |z − z slice edge | < 0.03, we shift the redshift bin by z∼0.05, then check the overdensity and galaxy members' redshift distribution to determine if they belong to the final cluster candidate sample.
The redshift information of all 17 accepted cluster candidates is listed in Table 1. Finally, we have a galaxy sample of 1396 galaxies.

Comparison with other cluster catalogues
We compare the final cluster candidate sample with several other cluster or group catalogues following different selection methods (e.g. overdensity, weak lensing, cosmic web structure). First, we compared it with the Chandra COSMOS X-ray catalogue (Gozaliasl et al. 2019). This work used extended X-ray sources detected on a 32-128 arcsec scale. The authors searched for the cluster or group members within the extended X-ray source area. And the redshift is mainly decided from spectroscopic redshift in Hasinger et al. (2018). We find three matches out of 14 clusters or groups in Gozaliasl et al. (2019) Table 2. Apart from the three cluster candidates mentioned above, another cluster in Gozaliasl et al. (2019) was found by our method, but it does not have NIR so we do not put it in our cluster candidate catalogue; two others were selected by overdensity mapping but were ruled out due to flat photometric redshift distribution. Considering the size of the extended X-ray source is much smaller than the structure we are looking for, we also tested using the smaller 0.5 Mpc radius aperture to search for overdensities again and reproduce almost all clusters or groups in Gozaliasl et al. (2019) except the least massive one. However, these smaller size structures are more likely to be groups with low masses (<10 14 M ) rather than clusters, so we did not use them for our analysis.
Then we compared our sample with optically selected cluster catalogues. There are three matches in Bellagamba et al. (2011). In this work the clusters are selected by weak lensing data. Our cluster candidates 05, 06, and 07 match with clusters 17, 19, and 24, respectively, in their work. The information for these three clusters is listed in Table 3. Since the weak lensing method detects clusters through dark matter spatial distribution rather than through galaxy overdensities, the difference is expected.
We also compared our cluster candidates with the COSMOS environment catalogue (Darvish et al. 2017), where the authors use photometric data from Laigle et al. (2016) to construct a density field. They identify the environment (cluster, filament, field) from the density field for each galaxy without giving cluster information. We slice the galaxy redshifts as in this work, plot cluster galaxies in each slice and compare them with our cluster candidates. Our cluster candidates (01, 02, 04, 05, 06, 07, 08, 10, 14, 15, and 17) match well with dense regions of cluster galaxy distribution. The other cluster candidates also contain cluster galaxies in the aperture, but their centres are offset from the densest regions.
In conclusion, 12 of our 17 cluster candidates have good matches with the clusters discovered in previous studies. In previous works, all 17 cluster candidates have a cluster match within 1 Mpc.

Methods
In Sect. 3.1 we describe how we classify RQGs, SFGs, and QGs. In Sect. 3.2 we explain how we separate fast-and slowquenching galaxies on the UVJ diagram. In Sect. 3.3 we explain how we remove the impact of possible projection effects.
A&A proofs: manuscript no. output Table 2. Properties of X-ray detected clusters taken from Gozaliasl et al. (2019). In the first column, in brackets are our own cluster IDs listed in Table 1. The last column shows the distance between the X-ray detection and our definition of the cluster centre.

Rest frame colour index estimation
We used the EAZY code (Brammer et al. 2008), which fits the SED by maximising the likelihood between the observed and template predicted SED. The templates of EAZY are set purely on stellar population synthesis models, which suit deep optical-NIR surveys. A non-negative matrix factorisation (NMF) algorithm (Blanton & Roweis 2007) is adopted to reduce the number of templates while keeping the ability to reproduce observations. For input of NMF, the library of PEGASE models described in Grazian et al. (2006) is used. This library includes ∼3000 models with ages from 1 Myr to 20 Gyr and several models for SFHs. These templates also consider dust extinction using the Calzetti et al. (2000) law and include emission lines. The EAZY code fits SED using a linear combination of these templates. Through the SED fitting, rest-frame UVJ colours can also be derived. In this work, we use rest-frame U and V filter information from Maíz Apellániz (2006) and rest-frame J band from 2MASS survey (Skrutskie et al. 2006). We fix the redshifts of galaxies to the photometric redshifts estimated in the COSMOS2015 catalogue (Laigle et al. 2016). As a result, EAZY provides the interpolated colour index L X (X=U, V, J). Then the U−V and V−J colour indices are calculated from the U, V, and J colour indices with the form

Rest-frame UVJ diagram
The rest-frame UVJ diagram is a common tool for separating QGs and SFGs (e.g. Daddi et al. 2004;Patel et al. 2012;Arnouts et al. 2013). QGs are redder in U−V colour for a given V−J colour, due to their strong Balmer break feature. Since RQGs are transitional galaxies between these two populations, a natural thought is that they should also reside in a region between these two populations on the rest-frame UVJ diagram.
The location of RQGs on the UVJ diagram ( Fig. 3) has been discussed in many works, such as Almaini et al. (2017) colours derived from galaxy photometry that correlates with galaxy spectral features, mainly D4000 and absorption line), the specific star formation rate (sSFR) obtained from SED fitting, Keck/MOSFIRE spectra, and LEGA-C spectra, respectively. It is noted that the boundaries between different populations are arbitrary. Taking advantage of the existing comparisons, we set our own selection boundaries to include as many RQGs as possible, while keeping a relatively low contamination level (see Sect. 5.1).
In this paragraph we briefly introduce how our boundaries are decided (Fig. 3). We modified the boundaries from Schreiber et al. (2018) (solid black lines), which separates SFGs (region 2, 3, 6, and 7) and QGs (region 1, 4, and 5). In this work dusty and non-dusty galaxies are divided by a diagonal line. We restrict it within the quiescent region to separate the post-starburst galaxies (region 1) and fully quenched QGs (region 4 and 5) as in Wild et al. (2014) and Almaini et al. (2017). However, some works (Almaini et al. 2017;Schreiber et al. 2018;Wu et al. 2020) identified RQGs outside the solid boundaries in both star-forming and quiescent regions by different methods (e.g. PCA method, SED fitting, spectral features). Therefore, we shifted the original solid horizontal boundary to U-V=1.1 (grey dashed horizontal line) and shifted the diagonal boundary by ±∆(U − V) = 0.1 (grey dashed diagonal lines) to add regions 2, 3, and 4 into the original post-starburst region 1. The final RQG sample (regions 1, 2, 3, and 4) should fulfil all the following conditions: (4) Figure 3 depicts all the boundaries used to classify these three galaxy populations. The SFGs (blue points) lie in region 7 and show a sequence depending on dust extinction. The QGs (red points) are in regions 5 and 6, and the dusty QGs are in region 6 (Almaini et al. 2017). The RQGs lie between the SFG and QG populations, in the middle of the evolutionary paths from SFGs to QGs, as shown by the solid curves (see Sect. 3.2). The RQG region (1, 2, 3, and 4) is shown with purple and orange shading. There is a concern that these arbitrary boundaries may lead to some uncertain results. However, we shift the boundaries between QG and RQG by 0.1 magnitudes and find that it does not affect our qualitative results.

Fast-and slow-quenching galaxy separation
The quenching timescale is an important factor in understanding its physical mechanisms. To put a strict constraint on the quenching timescale, it is necessary to obtain high-resolution, high S/N rest-frame optical spectra to fit the galaxy model. However, the number of our sample galaxies with spectra is limited. To give an idea of how the quenching timescale depends on different parameters, we take a photometric approach and apply model colour tracks to separate galaxies with different quenching timescales.
The SFH model we adopt has a constant star formation period followed by an exponentially declining star formation quenching phase, as illustrated in Fig. 4. The star formation rate (SFR) evolution with age is set by Eq 5 where τ is the quenching timescale.
To understand how galaxies move in the rest-frame UVJ space, we used the CIGALE code (Boquien et al. 2019) to make galaxy models and estimate their rest-frame colours. In the estimation of rest-frame UVJ colours, we used exactly the same filter as in Sect. 3.1.1. In this work we assume the galaxies have Chabrier (2003) initial mass function, solar metallicity, and Calzetti et al. (2000) dust attenuation curve with E(B−V)=0.15. Since the RQGs have some remaining star formation, there should be a certain amount of dust. And we arbitrarily choose the value of E(B−V)=0.15 for galaxies in the RQG phase. More details of the effects of dust attenuation on the evolutionary paths are discussed in Appendix A. For a fast quenching galaxy we assume a quenching timescale of 0.1 Gyr, and for a slowquenching galaxy we assume a quenching timescale of 1 Gyr. In Appendix A we discuss how different parameters affect evolutionary pathways in the UVJ space. The evolutionary pathways of the two model galaxies are plotted in Fig. 3. We can see there is a clear difference between the fast and the slow quenching tracks. To separate the fast (purple) and slow (orange) quenching regions, we adopt two lines. For regions 1 and 4, we adopt the solid diagonal line from Schreiber et al. (2018). For regions 2 and 3, we draw a solid red line with the same slope as the previous one where the slow quenching track enters the RQG region.
CIGALE also allows us to estimate the time when the model galaxies stay in the RQG selection area, and we define this time as the visibility timescale. When the quenching timescale is 0.1, 0.3, 1, and 2 Gyr, the visibility timescale is 0.43, 0.58, 1.5, and 2.6 Gyr, respectively. Generally, galaxies with a longer quenching timescale have a longer visibility timescale. This trend will help us understand our results better. However, since the galaxy models here are simplified, the value of the visibility timescale will not be used for quantitative analysis.

Field subtraction
In order to select galaxy cluster candidates we used photometric redshifts adopting redshift intervals (full width) of ∆z = 0.1 at 0.5 < z < 1. However, the intrinsic scatter of a cluster in redshift space is much smaller than ∆z = 0.1. Therefore, there should be many galaxies that are not physically associated with our cluster candidates in each redshift slice. Thus, we need to apply a statistical field subtraction to correct for such a projection effect. To achieve this we first define the field regions following the same method applied during the cluster candidate search. We use aperture density maps for each redshift slice and define the field area as the regions with densities lower than the 1σ threshold, which corresponds to the regions outside the cyan contours in Fig. 1. Information of the field samples in each redshift slice is shown in Table 4.
Before performing the field subtraction, we divide our sample of cluster and field galaxies into several stellar mass bins, which will be used to carry out our environmental analysis in Sect. 4. Finally, we compute the density of each population (SFG, RQG, QG) within each stellar mass bin and subtract the field contribution from each population using the equation where N subtracted is the number of galaxies after the subtraction, N is the number of galaxies before the subtraction, Σ k is the field density of each population, and N cluster is the total number of the Table 4. Information of the field sample in each redshift slice. The comoving area and the number of each population are given for each slice. Fast RQG and slow RQG are RQGs located in fast (purple) and slow (orange) quenching regions in the UVJ diagram (Fig. 3). cluster candidates in this redshift bin. While studying the environmental dependence, galaxies are also binned in local density (Σ 5 ) during our analysis. Thus, we need to subtract possible field galaxies from each local density bin as well. After computing the field density of each population, we then subtract it by applying where Σ k is the field density of each population and Σ f ield is the field density. The first item in brackets is the probability that a cluster candidate member galaxy is a field galaxy; the second item is its probability of being an SFG, an RQG, or a QG. After calculating this probability for each cluster candidate member, we then sum up this value for all galaxies in the local density bin.

Results
In this section, we present the results of our analysis in our selected cluster candidate sample (Table 1). The average field densities are displayed in Table 4. All errors in this section are propagated based on the Poisson statistics for both the original data and the field subtraction using a standard method. The information of all the bins analysed in this section is listed in the tables in Appendix B.

Quenching efficiency
Galaxy quenching drives SFGs to evolve into RQGs. Thus, the ratio of these two populations can represent the efficiency of the quenching process. Hence, we define the quenching efficiency as Figure 5 presents how Q.E. depends on stellar mass. The left panel also shows a comparison between the cluster candidates and the general field. Q.E. shows a strong dependence on the stellar mass, especially in the most massive bin. The Q.E. is always higher in the cluster candidates, showing that the cluster candidate environment is enhancing it. The gap is much more prominent in the most massive bin, indicating that the mass quenching is also stronger in the cluster candidates. To see how the local galaxy density affects the mass-dependent quenching, we divide our sample into high-and low-density subsamples at 10 1.3 /M pc 2 (i.e. 20/M pc 2 ), each of which contains one-half of the galaxy sample (middle and right panels in Fig. 5). It turns out that the local environment does not have a substantial impact on mass-dependent quenching, and that there is a significant rise in the Q.E. with stellar mass for both sparse and dense local environments.
In the left panel of Fig. 6 the Q.E. in the cluster candidates shows a dependence on the local density, which supports the scenario of environmental quenching. In the field the Q.E. is slightly lower than the lowest local density bin in the cluster candidates, indicating that the outskirts of the cluster candidates have a similar enhancement in galaxy quenching to that in the field. The sample is divided into low-mass and high-mass subsamples at M * = 10 10.42 M for analysis. In both of these subsamples, Q.E. increases with local density. However, considering the errorbars, we can only say that the environmental dependence in low-mass systems (the middle panel of Fig. 6) is significant. If we compare the Q.E. in the two subsamples, high-mass galaxies have higher Q.E., which again confirms the presence of massdependent quenching.

Quenching stage
When the quenching process is finished, RQGs become QGs. Therefore, we define quenching stage as the number ratio of these two populations where the Q.S. value indicates how much the quenching process has proceeded for a given galaxy sample. When the Q.S. value is high, many galaxies are in the RQG phase and are undergoing quenching, indicating a relatively early stage of quenching. On the other hand, if the Q.S. value is low, it means that most of the galaxies are already quenched and have become old QGs. As can be seen in the left panel of Fig. 7, the Q.S. decreases with stellar mass, especially for the least massive bin. The Q.S. in the field is higher than in the cluster candidates, showing they are in an earlier phase of quenching. The lowest mass galaxies are at an earlier stage of quenching, which is consistent with the downsizing scenario in which massive galaxies form and quench earlier than low-mass ones.
We then try to reveal how the local environment affects this relation (Fig. 7, centre and right diagrams). The full sample was divided into a high and a low local density subsample using Σ 5 at 20/M pc 2 as a threshold, and ensuring that each subsample contains half of all the galaxies. The trend between the Q.S. and stellar mass is only significant for the denser local environment. As in Fig. 8, the Q.S. clearly decreases with local density. If we separate massive and less massive subsamples at M * = 10 10.42 M (Fig. 8), the trends are still seen in both subsamples, but are less    significant due to decreased sample sizes. The Q.S. in the field is lower than the least dense cluster bin. However, the current field sample size is large, and the errorbar will increase if it has a similar size to that in the cluster. Therefore, this result is not significant and we do not discuss it further here. We also note that mass dependence shows different properties in two local density bins, and environment dependence shows different values in two mass bins. Hence, the mass dependence of quenching also has a dependence on the environment, and the environment dependence of quenching also has a dependence on stellar mass.

Mass, environment, and redshift dependence
As we have emphasised, the quenching timescale is an important factor for discriminating among the different physical mechanisms of quenching. Although we separate fast-and slowquenching galaxies only based on the rest-frame UVJ diagram, we can still use the number ratio of these two populations to illustrate the qualitative dependence of the quenching timescale on mass and environment.
Article number, page 9 of 20 A&A proofs: manuscript no. output  Fig. 6). The number to the left of each data point is the number of RQGs in this density bin, while the number to the right is the number of QGs. The total number of QGs or RQGs in the middle and the right panels is the same as the number of cluster candidates in the left panel. Errors are propagated based on the Poisson statistics for each population in each bin.  In Fig. 9 we show that the number ratio of slow-quenching galaxies over fast-quenching galaxies increases with the stellar mass, which means that low-mass galaxies have shorter quenching timescales compared to those of massive galaxies. This suggests that the quenching of galaxies with different stellar mass is dominated by different physical mechanisms. In Fig. 10 there is a possible dependence of the quenching timescale on the local density. The ratio seems to have a peak in the intermediate density bin with large uncertainty. However, we cannot say whether it is related to any physical trend since our sample size is still limited. Figure 11 shows a clear trend that fast quenching is more important at higher redshifts. This suggests that the main quenching mechanism changes with redshift, and the short timescale quenching becomes more important at higher redshifts, which is consistent with the results in Cucciati et al. (2012) and Rowlands et al. (2018).

Two specific cases
Here we highlight two specified cases in quenching that show the largest difference in the slow and fast quenching ratio. These cases probably correspond most directly to the mass and environment-dependent quenching, respectively.
In Table 5 the first case is the quenching of massive galaxies in a low-density environment. High mass refers to the higher 50% of mass (i.e. > 10 10.42 M ). Low local density refers to the lower 50% of local density (i.e.< 10 1.3 /M pc 2 ). This case should be mainly affected by mass-dependent quenching. The N(slow)/N(fast) ratio has a high value, giving a hint that the mechanism working in this case has a long timescale. This is likely to be, for example, AGN feedback, given that AGNs are Table 5. N(slow)/N(fast) values and their errors in the two specified cases. The high-mass low-density case may be more affected by massdependent quenching, and the low-mass high-density case may be more affected by environmental quenching. Details are discussed in Sect.

4.3.2.
Description N(slow)/N(fast) high mass + low local density 6.500 ± 2.264 low mass + high local density 1.308 ± 0.343 exclusively hosted by massive galaxies, and the AGN feedback can be a slow quenching process with a timescale of a few Gyr (Hirschmann et al. 2017). The second case is the quenching of low-mass galaxies in a dense environment, which is expected to be more affected by the environment. The N(slow)/N(fast) ratio is much lower compared to the first case, showing the impact of a short timescale mechanism. It is likely to be, for example, ram pressure stripping, which strips gas from low-mass galaxies more easily due to their shallow gravitational wells. Moreover, ram pressure stripping can only occur in dense cluster cores on a relatively short timescale (Steinhauser et al. 2016). However, short timescale AGN feedback quenchings also exist, and they can play a role in this case as well.
We note that we cannot reject other possibilities since the information we are able to obtain with photometric data alone is still limited.

Visibility timescale
In addition to the reasons we analysed in the previous sections, we should note that the visibility timescale of galaxies can also play a role in our results. The visibility timescale is a galaxy's timescale of staying in the RQG region on the UVJ diagram. If the visibility timescale also depends on stellar mass or local density, it would degenerate with the Q.E. and Q.S. In Sect. 3.2 we show that galaxies with longer quenching timescales generally spend a longer time in the RQG region on the UVJ diagram. However, with photometric data solely, we cannot constrain the visibility time with a small enough uncertainty to completely break the degeneracy between the mass and environment dependence of the visibility time and quenching properties. Hence we only separate the fast quenching (short visibility time) and the slow quenching (long visibility time) bins, and estimate their contributions to determine how the visibility timescale affects our results.
In Fig. 12, we show how the fast-quenching (short visibility time) galaxies and the slow-quenching (long visibility time) galaxies contribute to the mass and the environment dependence of Q.E. and Q.S . In the mass dependence of the Q.E. panel (top left panel), we can see that the long visibility time galaxies have a strong mass dependence, while the short visibility time galaxies do not. The figure shows that the visibility time affects the Q.E.; the increasing trend of Q.E. is partly due to the mass dependence of visibility time, and the massive galaxies tend to stay in the RQG region for more time. In the local density dependence of the Q.E. panel (top right panel), both fast-and slowquenching galaxies show an increasing trend, indicating that the visibility time does not have a strong impact on the environment dependence of Q.E. This is consistent with Sect. 4.3, where we show that the quenching timescale has a dependence on the stellar mass but not on the local density. In the Q.S. panels (bot-tom two panels), both fast-and slow-quenching galaxies show a decreasing trend (similar to the overall trend), showing that the visibility time does not strongly degenerate with the mass or environment dependence of the Q.S.
In Sect. 4.3.1 the degeneracy between the quenching timescale and visibility timescale exists as well. However, there is no proper way to break this degeneracy with only galaxy photometry. We should keep this degeneracy in mind when interpreting the physical scenario of quenching.
In conclusion, while analysing the mass dependence of Q.E., we should be careful of the influence of RQG visibility time. We leave this issue to our future work. When we have the spectra of RQGs, we can disentangle the degeneracy between the quenching properties and the visibility time.

RQG comoving number density
As mentioned in Sect. 3.3, the RQGs, fast-quenching RQGs, and slow-quenching RQGs all show relations with redshift. To better understand this dependence, we calculated the comoving number densities in each redshift bin. In this section not only do we use the galaxies in our 17 cluster candidates, but also all the galaxies in the COSMOS field with M * > 10 10 M .
In Fig. 13 we present the number densities of fast-quenching galaxies, slow-quenching RQGs, and total RQG sample in different redshift bins (including both field and cluster candidate regions). We bin all galaxies into two groups to increase the statistics. The low-z bin range is 0.5 < z < 0.75, while the high-z bin refers to 0.75 < z < 1.0. Both slow-and fastquenching RQGs show an increasing trend with redshift, and the fast-quenching RQGs show a steeper slope. The increasing trend of fast-quenching RQGs is consistent with Belli et al. (2019) and Wild et al. (2016) who reveal that the number density of poststarburst galaxies (corresponding to fast-quenching RQGs in this work) increases with redshift. Comparing the slope of slow-and fast-quenching galaxies, we can tell that, similar to the trend in cluster candidate regions, fast quenching is more important at higher redshifts (see Sect. 4.3.1).

Discussions
In Sect. 5.1 we discuss the contamination and completeness of our photometric RQG selection method. In Sect. 5.2 we discuss the physical implications of our results.

Spectroscopic confirmation of our RQG selection
Our current result is based on photometric data. However, the UVJ selection is not fully confirmed by galaxy spectra, while the only way to precisely determine a galaxy's RQG properties is through its spectral features. The main limitation of RQG research is that it is a rare population, making it expensive to obtain a large sample of RQG spectra. Fortunately, the COSMOS field already has deep spectroscopic data available for many galaxies, enabling us to match our targets with the existing spectra.
We use spectroscopic data from DEIMOS 10K spectroscopic survey catalogue (Hasinger et al. 2018) and LEGA-C survey data release 3 (van der Wel et al. 2021). The DEIMOS 10K survey obtained the spectra of some of the COSMOS2015 catalogue galaxies. Depending on the galaxy, the spectral resolution can be either R∼ 2000 or ∼ 2700, which are both large enough to measure the equivalent width of absorption lines. The LEGA-C survey obtained galaxy spectra with R=2500. The parent sam- Fig. 12. In this figure, we show results of Q.E. and Q.S. in fast-and slow-quenching galaxies separately. Slow-quenching galaxies (same as long visibility timescale galaxies) are marked with orange points, while fast-quenching galaxies (same as short visibility timescale galaxies) are marked with light blue points. Fig. 13. Comoving number density evolution of all RQGs, slowquenching RQGs, and fast-quenching RQGs. To increase the statistics, all galaxies are put into two bins, 0.5 < z < 0.75 (low redshift) and 0.75 < z < 1.0 (high redshift), consisting of 1468 and 2419 RQGs, respectively.
ple LEGA-C adopts for observation is the UltraVISTA catalogue presented in Muzzin et al. (2013).
The spectral features characteristic of RQGs are strong Balmer absorption lines (Couch & Sharples 1987). We measure the galaxies' Hδ absorption line equivalent width (EW), representing the Balmer absorption. For all the galaxies within our sample, galaxies with other strong Balmer absorption lines (e.g. Hγ, Hβ absorption line) have strong Hδ absorption lines as well. Another feature of RQGs is that their star formation activities are at a low level, so we use [OII] emission line strength to examine the existence or absence of star formation activity in the galaxy in this study. We define a galaxy with EW(Hδ)<-4Å and EW([OII])<5Å (negative value represents absorption) as an RQG galaxy in this section.
The equivalent width of each line feature is calculated using together with Lick spectral indices (Salaris & Cassisi 2005) and the index of the [OII] line from Balogh et al. (1999). Here F I,λ is the index flux, which is obtained from the index band; F C,λ is the continuum flux, which we assume to be constant over the index band. The continuum flux is calculated from the average flux over blue and red continuum bands. We do not remove the infilling emission in the calculation of EW(Hδ). All the spectral index information is listed in Table 6. For the DEIMOS 10K observations, in addition to matching our RQG sample with their spectroscopic catalogue, we also checked the spectra they labelled with the Hδ line feature. We have 105 matches in the DEIMOS 10K catalogue. For the LEGA-C survey, we only matched our RQG sample with their catalogue. We matched the target positions within 1", and we find 23 matched RQGs. However, we find that there are spectra in both surveys that cannot provide a high enough S/N for accurate absorption line EW measurement, for which there are various reasons. In particular, the DEIMOS 10K spectra mostly have exposure times from 1 to 3 hours, which is not long enough for less luminous galaxies. To avoid contamination from poorquality spectra, we select galaxies with EW errors of all four measured spectral lines σ i < 0.8Å (i=Hδ, [OII]). This error cut leaves us with 40 high-quality spectra (28 from DEIMOS 10K and 12 from LEGA-C). Another thing we should note is that the DEIMOS spectrograph has two chips that produces an artificial break in the spectrum at around 7600Å . To avoid this effect, we checked all good-quality spectra and remove those with this break in either the continuum or the index band of any of the four lines. There are RQG spectra without Hδ coverage in both DEIMOS 10K and LEGA-C, and we remove these galaxies as well. In total, we are left with 27 UVJ-selected RQGs with good enough spectra for spectroscopic confirmation. Their information is listed in Appendix C. The positions of these 27 galaxies on the UVJ diagram are plotted in Fig. 14. Among the 27 RQGs, 14 galaxies have strong Balmer absorption lines (EW(Hδ)<-4Å). For the 13 galaxies below this threshold, we checked their colour index errorbars on the UVJ diagram (red marks in Fig. 14). Based on the fluxes and errors given in the COSMOS2015 catalogue, we randomly generate 1000 flux values following Gaussian distribution and then calculate rest-frame colour indices from these fluxes. In this way we can calculate the 3σ errors for all galaxies. Nine galaxies out of 13 are close to the boundary. We also cross-matched these 13 galaxies with Chandra COSMOS Legacy Survey data (Civano et al. 2016), and found two galaxies with Hδ emission Fig. 15. Two checked galaxy spectra (in observed frame). Spectral features are indicated with green dashed lines. An artificial break between two chips of the DEIMOS spectrograph is shown as a grey shaded area. Top panel: Typical spectrum for strong Balmer absorption line galaxy, at z=0.864. Bottom panel: Typical spectrum for weak Balmer absorption line galaxy, at z=0.553. lines (ID 908883 and 635561) are X-ray AGNs. We should also note that all these galaxies with weak Balmer absorption features lie in the slow-quenching region, which confirms that the quenching timescale of this UVJ region is longer, that their stellar populations have already grown old, and thus that the Balmer absorption feature has started to fade. Six of the 13 galaxies (ID 314012, 375916, 391497, 438059, 559691, 810957, 816426, and 837551) have -4Å <EW(Hδ)<-1.5Å; meanwhile, their spectra are passive. This is consistent with what we discuss in Appendix A, Fig. A.2. In Fig. A.2 we show the evolution of EW along the evolutionary paths of model galaxies. The galaxies in slowquenching regions show weak Hδ absorption lines (as weak as EW(Hδ)∼1.5Å). The remaining three galaxies are possible contamination. Therefore, we conclude that the overall contamination rate is not high.
However, in the 14 strong absorption line galaxies, 8 (6 of which are slow-quenching galaxies) also have strong [OII] lines (EW([OII])>5Å). This is consistent with the work of Rudnick et al. (2017), who find [OII] emission in quenched galaxies. The [OII] emission may indicate that the depletion time is indeed longer for slow-quenching galaxies, keeping the galaxy forming stars for a long time. However, we still need further work to fully understand the source of [OII] in RQGs. In Fig. 15, we present typical examples of one strong and one weak absorption galaxy spectra.
In addition to the RQGs, we found two QGs and two SFGs that satisfy the RQG criteria in all high-quality spectra (Fig. 14). Their information is listed in Appendix C. The SFG at the bottom (ID 897743) has clear AGN features (i.e. a broad MgII emission line, broad Hγ and Hβ emission lines). The other SFG and one of the QGs are near the RQG region boundaries. However, for the other QG, which is far away from any boundary, we may have a risk of missing a certain subpopulation of RQGs that, for some reason (e.g. undergoing stronger starburst before quenching) spend more time in the recently quenched phase.
In conclusion, our photometric UVJ selection method still suffers from certain contamination and completeness problems, but judging from the current spectroscopic data, neither the contamination nor the incompleteness is high. Hence, our results from photometric data are still trustworthy. Due to the lack of high-quality spectra for absorption measurements, we do not have a statistical sample to estimate a reliable contamination or completeness rate. Deeper and wider spectroscopic observations of RQGs are required to better assess our photometric method.

Physical scenarios of quenching processes
In Sect. 4.1 we show that the Q.E. increases with both the stellar mass and the local density. From the Q.E., we can see that the mass dependence trend is not greatly affected by the local density, but the efficiency is higher in denser local environments, especially in the highest local density bin. It shows that mass quenching happens in all kinds of environments, but the cluster environment (i.e. dense environment) can also enhance the quenching of massive galaxies. The most significant increase in efficiency appears in the most massive bin, indicating that mass-dependent quenching is much stronger for the massive galaxies. The Q.E.'s dependence on local density is not independent. It is also affected by the stellar mass, so that for massive galaxies, this dependence almost disappears (see the middle and right panels of Fig. 6). Therefore, although the environmental effects of quenching also exist, they are only efficient in affecting low-mass systems. Our definition of Q.E. is different from commonly used definitions based on the fraction of QGs (e.g. Kawinwanichakij et al. 2017;van der Burg et al. 2020;Reeves et al. 2021). Since our definition introduces the RQG population, which may still be in the process of quenching, it can convey the efficiency more directly. In Kawinwanichakij et al. (2017) the authors separate mass-dependent quenching efficiency and environmental quenching efficiency. To compare our result with theirs, we also calculated our mass quenching efficiency using their definition (see equation 5 in Kawinwanichakij et al. 2017). In Kawinwanichakij et al. (2017) the mass-dependent quenching efficiency increases with stellar mass; the trend is consistent with our result. In the least massive bin our quenching efficiency is ∼0.3 dex higher, and in other bins our efficiencies are ∼0.1 dex higher than the Kawinwanichakij et al. (2017) result. This might be due to different sample completeness or cosmic variance of our data. The environmental quenching efficiency in Kawinwanichakij et al. (2017) does not show a significant evolution with stellar mass. In our work we only use two stellar mass bins while discussing environmental quenching, but instead we dis-cuss its dependence on local density. In general, low-mass galaxies in our sample have lower efficiency, but this strongly depends on local density. We note that in Kawinwanichakij et al. (2017) the efficiency at the massive end has extremely large errorbars, and therefore we cannot establish a reliable comparison in that mass range. In van der Burg et al. (2020) the authors discuss quenching efficiency as a function of stellar mass; the resulting trends and binned values both agree with our result. In Reeves et al. (2021) the authors discussed how the quiescent galaxy fraction depends on stellar mass. The trends they found agree with our results, except in the least massive bin where our quiescent galaxy value is slightly higher than theirs. This is likely a consequence of the different redshift range studied in their analysis (1.0 < z < 1.5). At that redshift low-mass galaxies are predominantly not quenched, hence their fraction is lower. The authors also pointed out that the environment and mass dependence of quenching are not separable, which is consistent with our result. McNab et al. (2021) studied the mass dependence of quenching by the excess of cluster SFG fraction over the field SFG fraction. The increasing trend of the excess represented by quenching efficiency agrees with our result. Webb et al. (2020) studied the SFR declination rate in field and cluster galaxies by spectroscopy and photometry fitting. Their results show that cluster galaxies (i.e. denser environment galaxies) quench faster, which also agrees with our results.
In Sect. 4.2 we see that the Q.S. decreases with both the stellar mass and the local density. The Q.S. provides us with a glimpse of quenching in the past. It significantly decreases with stellar mass at the low-mass end (similarly for the cluster candidates and the field), showing that more massive galaxies are at a much later stage of quenching. This is consistent with the downsizing scenario (e.g. Cowie et al. 1996;Cattaneo et al. 2008;Webb et al. 2020) that massive galaxies form and quench earlier, and hence they are at a later, more advanced stage of evolution compared to low-mass galaxies. The mass dependence of the quenching stage is strongly affected by the local density. The trend is only significant in high-density environments. For sparse environments, the errorbars are too large for us to say anything. The scenario is consistent with Hatch et al. (2011), where faster evolution in protoclusters is found. This can also be due to the migration of galaxies from the low-density cluster outskirts to the high-density cluster cores where the quenching happens. The Q.S. also shows a stable decreasing trend with the local density, indicating that galaxies in denser environments (cluster cores) are at a later, more advanced stage of quenching. The earlier Q.S. of field galaxies also supports this scenario. This is consistent with the inside-out quenching in clusters, where galaxies in cluster cores quench early, while the galaxies in the outer parts quench progressively later, which is consistent with the scenario in Koyama et al. (2013) and Shimakawa et al. (2018). Another possibility is that during the in-falling process their gas will be stripped off, which directly leads to their quenching (Bekki et al. 2002). In this case the effect should be stronger for low-mass galaxies since their gas components are easier to strip. This environmental dependence is not strongly affected by the stellar mass, which supports the inside-out quenching that has an impact on low-and high-mass galaxies equally. In McNab et al. (2021) the authors use three types of transitional populations: green valley (GV) selected by rest-frame (NUV-V) and (V-J) colours; blue quiescent (BQ) selected by rest-frame (U-V) and (V-J) colours; and post-starburst (PSB) selected by galaxy spectra. They connect the quenching timescale to the transitional population fraction excess in clusters. They found that the lowest mass systems (< 10 10.5 M ) have short quenching timescales (<1 Gyr), which agrees with our result that low-mass galaxies are dominated by fast quenching (<1 Gyr).
The slow-to-fast quenching ratio is an indicator of the quenching timescale. In Sect. 4.3 we find that the massive galaxies are dominated by slow-quenching, and the less massive galaxies are dominated by fast quenching. This result is consistent with Walters et al. (2022), where the authors find that slow quenchers tend to have higher stellar mass. The authors used IllustrisTNG simulation (Nelson et al. 2019) to estimate the quenching timescale. Walters et al. (2022) also used MaNGA IFU survey data (Bundy et al. 2015) and took measurements of the stellar age gradient from Woo & Ellison (2019) to split fast and slow quenchers. In this way the authors confirmed that slow quenchers have higher stellar mass. In Hahn et al. (2017) the authors matched N-body simulation from Wetzel et al. (2013) and their evolutionary model with SDSS DR7 Abazajian et al. (2009) data to investigate and constrain the quenching timescale. They find that satellite galaxies have shorter quenching timescales than the central ones. This result broadly agrees with this work, assuming that satellite galaxies are less massive compared to central galaxies. However, we note that the aperture we use for cluster member galaxy selection in this work is about a factor of 2 larger than R 200 in our clusters (see Table 2 and R 200 in Table 2 of Gozaliasl et al. 2019). Therefore, we cannot exclude the possibility that some galaxies far from the cluster centre could also be considered central galaxies of smaller haloes infalling towards the main cluster core.
In Sect. 4.3.2 two specified cases are discussed. These two cases show clearly different timescales of quenching. The massdependent quenching in sparse regions should be less affected by the environment. It turned out to have a longer timescale, and very likely to be AGN feedback that caused quenching. This is consistent with the observed timescale of quenching for this process (Hopkins et al. 2006) and with the typical stellar mass of host galaxies (Kauffmann et al. 2003). AGN feedback also plays an indirect role in the morphology transformation by suppressing in situ star formation, and AGN feedback is essential for forming ellipticals (Dubois et al. 2016). This could also explain the bimodality in galaxy morphology. The environmental quenching is only significant for low-mass galaxies, and it is surely stronger in denser regions. Quenching, in this case, appears to be a short timescale type, which we reasonably speculate to be caused by ram pressure stripping. As mentioned in Sect. 2.2.4, three cluster candidates are detected in X-rays, and their total masses are measured, which indicates that they are massive enough to allow ram pressure stripping to happen.
In Sects. 4.3 and 4.5 we performed a preliminary exploration of the quenching's dependence on redshift. In both cluster candidates and the general field, the fast-quenching RQGs increase with redshift. In cluster candidates the ratio of slow-quenching RQG to fast-quenching RQG decreases with redshift, indicating that fast quenching becomes more important at higher redshift. In Tinker & Wetzel (2010), Balogh et al. (2016), Fossati et al. (2017), andFoltz et al. (2018), the authors discussed the galaxy quenching timescale's dependence on redshift in clusters and groups. Combining their results, the quenching timescale also decreases with redshift. We estimated the comoving number density of fast-quenching and slow-quenching RQGs in the field. Fast-quenching RQGs show a strong increasing trend. This trend is consistent with the post-starburst analysis in Wild et al. (2016) and Belli et al. (2019). The redshift range they investigate is wider than ours, but both show that the number density of post-starburst galaxies increases with redshift up to z∼2. Wild et al. (2016) analyses the PSBs in the UKIDSS UDS field. To compare the two sets of results, we use the PSB catalogue from Wilkinson et al. (2021) that includes PSBs from the COSMOS field. To make a direct comparison, we limit our RQG sample to the traditional PSB region in Wild et al. (2014) (region 1 in Fig. 3). We also apply the same mass cut for both samples. Then we compare the comoving number density of two samples in Fig. 16. The increasing trend is found in both samples, while the absolute value has ∼0.4 dex difference, especially in higher redshift bins. We note that Fig. 1 in Wilkinson et al. (2021) shows a spectroscopic confirmation of their PSB sample, where a large fraction of PSBs are missed by their PCA selection. The incompleteness of PCA selection may lead to the difference between the two samples.
In this work, slow-quenching RQGs also show an increasing trend, although it is much weaker. This indicates that fast quenching processes become more important at higher redshifts. At lower redshifts, pre-processing, which is also a long timescale mechanism (Taranu et al. 2014), may be an additional quenching mechanism. This may not be as common at high redshifts because not many dense groups would exist. Therefore, the growing influence of pre-processing can result in the rising importance of slow quenching at lower redshift.
In conclusion, our photometric selection method is also a powerful tool in determining target galaxies for large systematic spectroscopic surveys of RQGs. As we apply this method to COSMOS galaxies in cluster candidates at 0.5<z<1.0, we obtain results with high consistency with previous works. However, our statistical analysis still suffers from the limited sample size and the uncertainty introduced by the photometric selection. We plan to increase our sample by including all Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP) DEEP fields (Aihara et al. 2019). If we include IRAC data, we can also expand our sample redshift range from the current 0.5-1.0 to 0.5-2.0, which traces the quenching process back to the cosmic noon. More importantly, spectroscopic confirmation of our sample's RQG features is required for further analysis. The preliminary spectroscopic confirmation shows low contamination and high completeness, but more spectra are required for more accurate confirmation. The high accuracy contamination and completeness rates of our UVJ selection method must be estimated and corrected to get the final results. Moreover, RQG galaxy spectra can be used to constrain quenching timescales through the spectral fitting. In this way, we can better distinguish between different mecha-nisms. We look forward to the Prime Focus Spectrograph (PFS) (Tamura et al. 2018) on Subaru that will come into use in 2023, when we can efficiently obtain the spectra of a large sample of RQGs due to the PFS's large field of view (1.3 deg 2 ) and great multiplicity (2,400 fibres).

Conclusions
In this study we constructed a multiband photometric sample of galaxies from 17 cluster candidates in the COSMOS field at 0.5 < z < 1.0. We then adopted the rest-frame UVJ selection method to classify RQGs, SFGs, and QGs. Moreover, we apply two exponentially declining SFH models with different quenching timescales to separate fast-quenching and slow-quenching galaxies on the UVJ diagram. We then correct for the projection effects introduced by the uncertainties of photometric redshifts and investigate different quenching scenarios as a function of stellar mass and local density. The main results of this study are summarised as follows: 1. We define Q.E. as the ratio of the number of RQGs to SFGs in a given volume. Q.E. increases with both mass and local density, confirming the existence of mass and environment-dependent quenching. Mass-dependent quenching is stronger for high-mass galaxies. 2. We define the Q.S. as the ratio of RQGs to QGs, aiming to quantify the history of quenching in a given volume. We find that this parameter also depends on both mass and local density. Massive galaxies are at a later stage of the quenching process, which is consistent with the downsizing scenario. Galaxies in the cluster core regions are at an earlier stage of quenching, which is a sign of inside-out quenching in cluster environments. 3. We compare the Q.E. and Q.S. of cluster and field galaxies.
We find that mass quenching is enhanced by the cluster environment, and the field galaxies are at an earlier Q.S. compared to the cluster galaxies of the same stellar mass. 4. The ratio of the number of slow-to fast-quenching galaxies is a proxy of the quenching timescale. In cluster environments the quenching timescale depends both on the stellar mass and on the redshift. Low-mass galaxies are dominated by fast quenching, while massive galaxies are dominated by slow quenching. Slow quenching is also more important at lower redshifts. 5. The quenching of massive galaxies in the sparse environment is dominated by the slow quenching process. This is likely due to the relatively long timescale of AGN feedback, as AGNs are hosted exclusively by massive galaxies. 6. Fast quenching is much more important to quench low-mass galaxies in the dense environment. We believe this is likely caused by short timescale processes such as ram pressure stripping, which only happens close to the cluster core (i.e. high-density environments) and affects low-mass galaxies more efficiently. 7. Visibility time of RQGs may degenerate with the quenching properties, especially with Q.E. and timescale. However, we cannot separate well the visibility time and the quenching properties without spectroscopic data. We leave this problem for our future spectroscopic work. 8. Existing spectroscopic data confirms that our photometrybased UVJ selection method is reliable enough. By measuring the EW of Balmer absorption lines and the [OII] emission line from DEIMOS 10K spectroscopic survey data, we can confirm the RQG nature of a limited number of galaxies with high-quality spectra. Although seven UVJ-selected RQGs do not show expected strong absorption lines, two of them are AGN host galaxies and four of them lie within the slow-quenching RQG region and show relatively weak absorption lines. Thus, only one of the 17 galaxies cannot be explained, and is likely contamination. Moreover, only three galaxies outside the RQG regions have RQG spectral features, and two of them are very close to the selection boundary. Combining all these pieces of evidence, the UVJ selection method is a powerful and reliable analysis tool to investigate the RQG population before more spectroscopic data become available.
The results we present in this work give evidence of both mass-dependent quenching and environmental quenching. They also shed light on the mass and environmental dependence of quenching timescales. With future spectroscopic data the physical properties of our sample galaxies can be better constrained, and the uncertainty of current photometric results can be reduced. The UVJ selection method we displayed in this paper can be applied to large photometric datasets (e.g. HSC-SSP data) of optical to NIR. We can efficiently select RQG candidates for further spectroscopic surveys performed by high-multiplicity multiobject instruments such as PFS.

Appendix A: Quenching model tracks on the rest-frame UVJ diagram
Here we investigate the evolution tracks of model galaxies in the rest-frame UVJ space and discuss how they depend on each parameter.
First, we discuss how the redshift might influence the UVJ selection criteria. In Fig. A.1 we show model galaxies with different star formation durations. The time after quenching started is indicated with star symbols in different colours. We can see that the evolution path is almost the same in the RQG and QG regions, indicating that a star formation duration difference of 3 Gyr will not change the selection criteria of RQGs. If we suppose that two galaxies started to form at the same time, but have two different durations of star formation (1 and 4 Gyr), they will show almost the same colours if the times after the start of quenching are the same. The observed epochs of these two RQGs can, for example, correspond to z=1.0 and 0.5, respectively. Therefore, we do not need to shift the selection boundaries depending on the redshifts in this range.
Then we discuss the evolution of spectral features along the model paths. We use the high-resolution mode of CIGALE to output high-resolution SED of model galaxies at each position and calculate the equivalent widths (EWs) of Hδ and [OII] line features. Figure A.2 shows the fast and slow quenching paths with blue and red curves, respectively. The EWs of Hδ (black) and [OII] (light blue) at each point are presented next to the position. A negative EW represents the absorption feature, while a positive EW represents the emission feature. All EWs are calculated using the same method as in Sect. 5.1. We can see that the Hδ absorption line strength first increases and then decreases along the evolution paths. In the RQG region, both the fast-quenching and slow-quenching galaxies show strong Hδ absorption. For the fast-quenching galaxy there is no [OII] emis- sion left. However, the Hδ absorption is weaker for the slowquenching galaxy compared to the fast-quenching galaxy due to its longer quenching timescale, and the [OII] emission is stronger due to remaining star formation. These features are also consistent with the spectral confirmation in Sect. 5.1. In conclusion, our boundary selection is reasonable and consistent with the spectral features of RQGs.
Finally, we discuss how dust attenuation affects the evolution paths. We use different dust parameters in CIGALE and calculate the evolution paths on the UVJ diagram. In Fig. A.3, we show the evolution path with E(B−V)=0.15 and 0.2. The paths are shifted to the top right, and the reddening vector of E(B−V)=0.1 is also plotted in the diagram. We note that for both E(B−V) values, the fast quenching and slow quenching tracks are still in the fast (purple) and slow (orange) quenching regions. Considering that RQGs should have a certain amount of dust as the fuel of the remaining star formation, the dust attenuation will lead to some contamination, but will not affect the separation of fast and slow quenching much.