# On the origin of the Σ1-M★ quenching boundary

A. Cattaneo, P. Dimauro, I. Koutsouridou

► To cite this version:

A. Cattaneo, P. Dimauro, I. Koutsouridou. On the origin of the Σ1-M★ quenching boundary. Monthly Notices of the Royal Astronomical Society, 2025, 537, pp.3929-3942. 10.1093/mnras/staf253. insu-05101681

HAL Id: insu-05101681

https://insu.hal.science/insu-05101681v1

Submitted on 8 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the origin of the $\Sigma_1-M_\star$ quenching boundary

A. Cattaneo,¹,²* P. Dimauro³ and I. Koutsouridou⁴

¹Observatoire de Paris, LUX, PSL University, 61 avenue de l'Observatoire, F-75014 Paris, France

²Institut d'Astrophysique de Paris, CNRS, 98bis Boulevard Arago, F-75014 Paris, France

³INAF – Osservatorio Astronomico di Roma, via di Frascati 33, I-00078 Monte Porzio Catone, Italy

⁴Dipartimento di Fisica e Astronomia, Università degli Studi di Firenze, Via G. Sansone 1, I-50019 Sesto Fiorentino, Italy

Accepted 2025 February 6. Received 2025 February 6; in original form 2024 July 10

## ABSTRACT

We have considered a phenomenologically motivated model in which galaxies are quenched when the energy output of the central black hole exceeds a hundred times the gravitational binding energy of the baryons in the host halo. The model reproduces the mass functions of star-forming and quiescent galaxies at $0 < z < 2.5$ and the quenching boundary on a $\Sigma_1-M_\star$ diagram. The quenching boundary arises because of the colour–morphology relation. The stellar surface density $\Sigma_1$ in the central kiloparsec is a morphological indicator. Galaxies becomes redder as $\Sigma_1$ increases until they cross the quenching boundary and enter the passive population. Mergers drive the growth of supermassive black holes and the morphological evolution that accompany the migration to the red sequence. That is the origin of the population of high-mass passive galaxies. At lower masses, passive galaxies are mainly satellites that ceased to form stars because of environmental effects.

**Key words:** galaxies: evolution – galaxies: formation – galaxies: star formation.

## 1 INTRODUCTION

The stellar surface density $\Sigma_1$ within the central 1 kpc of a galaxy increases with the stellar mass $M_\star$: $\Sigma_1 \propto M_\star^{0.9}$ for star-forming galaxies and $\Sigma_1 \propto M_\star^{0.7}$ for passive galaxies (Barro et al. 2017). For a given $M_\star$, passive galaxies have a higher $\Sigma_1$ than star-forming galaxies. The loci of passive and star-forming galaxies on the $\Sigma_1-M_\star$ diagram are separated by a ‘quenching boundary’ $\Sigma_1 \propto M_\star^{0.66}$ (Chen et al. 2020).

By arguing that $\Sigma_1$ is a proxy for the mass $M_\star$ of the central black hole (BH) and that $M_\star$ is a proxy for the virial mass $M_{vir}$ and the virial velocity $v_{vir}$ of the dark matter (DM) halo, Chen et al. (2020) proposed that the quenching boundary corresponds to a critical BH mass $M_{\star,quench} \propto M_{vir}v_{vir}^2$. Quenching occurs when the energy deposited by the BH into the surrounding gas ($\propto M_\star$) is larger than a fraction or a multiple of the gravitational binding energy of the baryons within the halo ($\propto M_{vir}v_{vir}^2$). The argument has the strength of being entirely empirical. The only theoretical assumption is the cold DM cosmology. Its weakness is that it relies on indirect proxy arguments.

In Koutsouridou & Cattaneo (2022), we used the GalICS 2.2 semi-analytic model (SAM) of galaxy formation to investigate the consistency of this scenario with other observational data. We assumed that feedback from the central BH (quasar feedback) blows out all the gas within the inner parts of a galaxy (i.e. the central starburst, the bulge and the bar, but not the disc) as soon as $M_\star$ is large enough to satisfy the quenching criterion:

$$ \epsilon_{\text{eff}} M_{\star} c^2 > \frac{1}{2} f_{\text{b}} M_{\text{vir}} v_{\text{vir}}^2 \quad (1) $$

where c is the speed of light, $f_b = 0.16$ is the universal baryon fraction, and $\epsilon_{\text{eff}} = 0.00115$ is a parameter of the SAM calibrated to reproduce the quenching boundary that separates SF and passive galaxies on the $M_\star-M_\star$ diagram (Terrazas et al. 2016).

Massive early-type galaxies have old stellar populations (e.g. Thomas et al. 2005). The physical processes that prevent gas accretion from the circumgalactic medium (CGM) from reactivating star formation today constitute a different problem (the so-called ‘maintenance problem’) with respect to what caused the initial quenching of star formation at high redshift. However, in Koutsouridou & Cattaneo (2022) and this article, we assume that, after the initial quenching, gas accretion is permanently shut down. The physics of maintenance are beyond the scope of our phenomenologically motivated analysis. We can nevertheless give two arguments for the plausibility of our assumption. The first is that quasar feedback heats the CGM to high entropy, so that its cooling time becomes long. The second is that quasar feedback is followed by the onset of a self-regulation mechanism preventing effective cooling (e.g. Cattaneo & Teyssier 2007).

With the assumption in equation (1), Koutsouridou & Cattaneo (2022) could reproduce the evolution of the mass function of galaxies since $z = 2.5$ as well as the fraction of passive galaxies at $z = 0$ as a function of $M_\star$ for central and satellite galaxies separately. Our SAM was also in good agreement with the $M_\star-M_\star$ relation and with the morphological properties of galaxies in the local Universe (i.e. bulge-to-total mass ratios). The current article improves our previous work in three ways.

(1) In Koutsouridou & Cattaneo (2022), the comparison with the observations was limited to $z = 0$ except for the evolution of the total mass function of galaxies. Here, we consider the evolution of the mass functions of star-forming and passive galaxies separately

<footnote>*E-mail: andrea.cattaneo@obspm.fr</footnote>
and we extend the comparison to other observables such as the radii of discs and bulges. Radii are important for the calculation of $\Sigma_1$ (see below).

(2) Chen et al. (2020) started from the $\Sigma_1-M_*$ relation and inferred a quenching criterion of the form in equation (1). Here, we close the circle by demonstrating that a model based on the quenching criterion in equation (1) reproduces Chen et al. (2020)'s quenching boundary on the $\Sigma_1-M_*$ diagram. The test is not trivial because equation (1) contains no condition on $\Sigma_1$.

(3) For Chen et al. (2020), the correlation between $M_*$ and $\Sigma_1$ was purely empirical (see also Sahu, Graham & Davis 2022 for direct observational evidence from galaxies with dynamical BH mass estimates). With GALICS 2.3, we can explore its astrophysical origin.

The structure of the article is as follows. In Section 2, we introduce our SAM and the novelties of the GALICS 2.3 version used for this article. In Sections 3 and 4, we present our results for the stellar mass function of galaxies and $\Sigma_1$, respectively. Section 5 summarizes our conclusions. In Appendix A, we discuss the impact that the mass resolution of the N-body simulation used to construct the DM merger trees could have on our results.

## 2 THE MODEL

A systematic presentation of our SAM would take several pages and detract from the focus of the article. Here, we point to the relevant literature and focus on the differences with previous versions.

Cattaneo et al. (2020) gave a systematic presentation of the GALICS 2.1 version (cosmology, merger trees, DM substructures, merger rates, accretion of gas on to galaxies and haloes, supernova feedback, galactic fountain) while referring to Cattaneo et al. (2017) for the internal structure of galaxies (density distributions, radii and characteristic speeds of discs and bulges), star formation rates (SFRs), and morphological transformations, where there had not been any changes since GALICS 2.0 (i.e. Cattaneo et al. 2017). Koutsouridou & Cattaneo (2019) improved GALICS 2.0 by introducing a model for ram pressure and tidal stripping.

### 2.1 Quenching

The main novelty of GALICS 2.2 (Koutsouridou & Cattaneo 2022) with respect to the previous versions 2.0 and 2.1 was the introduction of a quenching mechanism to suppress star formation in massive galaxies. In Koutsouridou & Cattaneo (2022), this mechanism could take two forms: halo quenching, where gas accretion on to galaxies is shut down above a critical halo mass (model A), and BH quenching based on equation (1) (model B).

Model A and B use the same shock stability criterion, described in Cattaneo et al. (2020) and further discussed in Tollet et al. (2022), to discriminate whether gas is accreted in the cold mode or the hot mode. The only difference with respect to accretion is the prescription used to decide at which point gas accretion on to galaxies is shut down. Once that happens, in both models: all the gas in the cold CGM is heated to the virial temperature; gas that accretes on to the halo is automatically incorporated into the hot CGM (the intracluster medium in the case of groups and clusters); and the hot CGM is no longer allowed to cool.

Model A and B also use the same prescriptions to grow supermassive BHs. The only difference is that model A contains no BH feedback. In model A, the growth of BHs is purely limited by how much gas is supplied to them. The quenching model assumed in the current GALICS 2.3 version is the same as in model B of Koutsouridou & Cattaneo (2022).

### 2.2 Morphological transformations in mergers

The other novelty of GALICS 2.2 was the presence of two options for morphological transformations in galaxy mergers. We referred to them as models 1 and 2.

Model 1 was the same as in previous versions. It assumes a very sharp distinction between major and minor mergers (mergers with mass ratios $\mu$ greater and smaller than 1:4, respectively). Only major mergers transform discs into bulges, trigger starbursts, and feed the growth of supermassive BHs.

Model 2 is based on hydrodynamic simulations by Kannan et al. (2015) and assumes a more continuous transition. Mergers transfer a mass fraction $\mu$ of the disc stars to the bulge and a mass fraction $(1 - f_{gas, disc})\mu$ of the gas in the disc to the central starburst ($f_{gas, disc}$ is the gas fraction in the disc). Hence, major mergers do not destroy discs entirely (especially gaseous discs), and minor mergers, too, can contribute the growth of bulges and BHs. The only difference between major and minor mergers is that major mergers are assumed to trigger starbursts throughout the discs and not only in the central regions of galaxies (Koutsouridou & Cattaneo 2022 for details). Model 2 is the default option in GALICS 2.3.

Other minor differences between GALICS 2.2 and the previous versions 2.0 and 2.1 were a new model for disc instabilities (based on numerical simulations by Devergne et al. 2020), a recalibration of the parameters for supernova feedback, and the star formation time-scale in discs. The GALICS 2.3 version used for this article is almost identical to the B2 variant of GALICS 2.2, which uses model B for quenching and model 2 for mergers. The only differences are in the density profiles and sizes of bulges (Section 2.3) and in the star formation law for discs (Section 2.4). We have also improved our model for tidal stripping by discovering and correcting a mistake by a factor of two (Section 2.5).

### 2.3 Density profile and sizes of bulges

Computing $\Sigma_1$ requires a model for the distribution of stars in galaxies, especially for the bulge component, which dominates the central region. Observationally, bulges are fitted with a Sérsic (1963) three-parameter profile. Many bulges exhibit central light excesses (cusps) or deficits (cores) with respect to a Sérsic fit (Faber et al. 1997; Kormendy et al. 2009a). The central light excesses of cuspy ellipticals are assumed to originate from gas that falls into the central regions of galaxies during dissipative mergers and triggers central starbursts (e.g. Faber et al. 1997; Hopkins et al. 2009a).

In our SAM, galaxies comprise four components: the disc, the bar, the bulge, and the cusp, which is present only if the last merger was dissipative. The cusp contains the stars formed in the central starburst at the last major mergers. The bulge corresponds to the pre-existing stellar population.

One of the main difficulties concerning the structural properties of galaxies is to determine the density distributions of the bulge and the cusp. We have explored a number of different models including one with two Sérsic profiles, one for the bulge and one for the central cusp (in addition to the exponential profile used to model the disc). However, introducing more free parameters than constraints just brings further uncertainty into the model. In the end, the model that works best at reproducing the sizes of bulges is the simplest one. In the current GALICS 2.3 version, we model the bulge and the cusp with a single Hernquist (1990) profile, so that the cusp is simply the innermost part of the bulge. The only purpose of maintaining a separate cusp component is to keep track of the stars formed during the last merger.
The Hernquist profile has two parameters. Hence, the mass $M_b$ and the energy $E_b$ of a bulge entirely determine its density distribution. The scale-length $a$ of the Hernquist profile is related to $M_b$ and $E_b$ by the equation (Hernquist 1990):

$$E_b = -\frac{GM_b^2}{12a}. \qquad (2)$$

Equation (2) follows from the virial theorem (the total energy is half the gravitational potential energy) and the assumption that there is not much DM or disc matter inside the bulge, so that one can consider the bulge self-gravitating and compute its gravitational potential energy from its density distribution alone.

The main interest of the Hernquist profile, besides its simple analytic form, is that its two-dimensional projection closely approximates a de Vaucouleurs (1948) profile with effective radius $R_e \simeq 1.8153a$ (the de Vaucouleurs profile is a Sérsic profile with index $n=4$). Once we find $E_b$, it is straightforward to compute

$$R_e = -0.151 \frac{GM_b^2}{E_b}. \qquad (3)$$

and thus the surface density profile of the bulge.

We compute the mechanical energy $E_b$ of a bulge formed after a merger event by assuming that the baryons that end up in the bulge conserve their mechanical energy:

$$E_b = \epsilon_{b1}E_{b1} + \epsilon_{d1}E_{d1} + \epsilon_{b2}E_{b2} + \epsilon_{d2}E_{d2}. \qquad (4)$$

Here $E_{b1}$, $E_{b2}$, $E_{d1}$, $E_{d2}$ are the mechanical energies of the bulges and the discs of the merging galaxies. The subscripts 1 and 2 refer to the larger and the smaller galaxy, respectively. The coefficients $\epsilon_{b1}$, $\epsilon_{b2}$, $\epsilon_{d1}$, and $\epsilon_{d2}$ correspond to the fractions of the baryons in the bulge of galaxy 1, the bulge of galaxy 2, the disc of galaxy 1, and the disc of galaxy 2 that end up contributing to the stellar mass of the final bulge, respectively. We assume that all the stars in the bulge of larger galaxy remain in the bulge component ($\epsilon_{b1} = 1$). The other three coefficients are computed following the prescriptions in Koutsouridou & Cattaneo (2022). Equation (4) does not contain any interaction term because we have assumed that the merging galaxies fall on to each other while starting from rest at infinity (as in the previous versions of our SAM).

The mechanical energies $E_{b1}$ and $E_{b2}$ of the bulges of the merging galaxies are computed from their masses and sizes with equation (2). The mechanical energies $E_{d1}$ and $E_{d2}$ of their discs are computed from the virial theorem by assuming that the mechanical energy $E_d$ of a disc is the opposite of its kinetic energy:

$$E_d = -2\pi \int_0^\infty \frac{1}{2} v_c^2(r) \Sigma_d(r) r dr, \qquad (5)$$

where $\Sigma_d(r)$ is the disc's surface density profile (a decreasing exponential) and $v_c(r)$ is the circular velocity (computed as in Cattaneo et al. 2017).

Rigorously, even in dissipationless ('dry') mergers, the conservation of energy and the virial theorem apply to the whole system and not to its individual components (e.g. the stars or the DM). To understand the impact that this may have on our results, we have looked at the GalMer data base of merger simulations (Chilingarian et al. 2010). In dissipationless simulations, the stars transfer to the DM 2 per cent of their energy in the worst case scenario. The accuracy of applying the virial theorem to the stars alone is measured by $2T_*/|U_*|$, where $T_*$ is the kinetic energy of the stars and $U_*$ is their gravitational potential energy, which is due not only to the stars' self-gravity but also to the DM's gravitational potential. Equation (5) assumes $2T_*/|U_*| = 1$. In GalMer, $0.80 < 2T_*/|U_*| < 0.95$. Our

approximations do introduce errors at the 20 per cent level but such errors are immaterial when considering the large scatter in both observed and modelled sizes (Fig. 1).

In dissipative ('wet') mergers, there is the additional complication that the gas that ends up in the cusp radiates before it forms stars. Hence, equation (4) overestimates $E_b$ and the scale-length $a$ of the Hernquist profile. Hopkins et al. (2009b) run hydrodynamic simulations of the formation of bulges in major mergers. They found that the presence of gas reduced the bulge size by a factor of $(1 + f_c/f_0)^{-1}$, where $f_c$ is the fraction of stellar mass formed during the mergers (i.e. the ratio of the cusp mass to the total bulge mass) and $f_0 \sim 0.3$.

Fig. 1 compares the mass–size relations for discs (blue points) and bulges (orange in points) in the CANDELS data (left), in GALICS 2.3 without dissipation (centre) and in GALICS 2.3 with $f_0 = 0.3$ (right). The blue and red lines are fits to the observational data for discs and bulges, respectively. They have been reproduced on the central column and the right column to assist the comparison with the observations, but one should note the important scatter of the data points around the mean relations.

In both the dissipationless and the dissipational model, the evolution of $R_e(M_*)$ with $z$ is weaker in GALICS 2.3 than in CANDELS. The dissipationless model (centre) is in reasonably good agreement with the data at all redshifts. Only the most massive galaxies have bulges that are too large at $z > 1.5$ and too small at $z < 1$. Adding dissipation (Fig. 1, right) reduces the sizes of the bulges formed in high-$z$ gas-rich mergers and brings GALICS 2.3 in better agreement with the CANDELS data at $z > 1.5$. Mergers at low $z$ have lower gas fractions. Hence, one would expect that their sizes should be less affected but that is not true because smaller progenitors result in smaller descendants. Dissipation improves the agreement with the data at $z > 1.5$ but spoils it at lower $z$ (right column). The dissipationless model is the one that reproduces the sizes of galaxies better across the entire redshift range $0 < z < 2.5$ and will thus be our reference model in the rest of the article.

There are two possible explanations why the crude model without dissipation works better. First, in real mergers a lot of the gas is in dense molecular clouds that hydrodynamic simulations cannot resolve and that behave like dissipationless N-body particles. Second, by assuming a Hernquist profile we have implicitly assumed that all bulges have the same Sérsic index $n=4$, while we know that is not true and that there are systematic trends of $n$ with $M_*$ (Kormendy et al. 2009b).

## 2.4 Star formation time-scales

Our SAM distinguishes between quiescent star formation in discs and bursty star formation in cusps and bars. This distinction admits one exception. In remnants of major mergers, all the gas is starbursting, including that in the disc component.<fnref n="1" /> Merger-driven starbursts continue until they have depleted all the gas present in the merging galaxies before the mergers. The assumption of bursty star formation does not apply to disc material accreted after major mergers, which forms stars quiescently. Star formation in cusps and bars is always bursty.

In Koutsouridou & Cattaneo (2022), the star formation timescales for quiescent and bursty star formation were 1 and 0.2 Gyr, respectively. That assumption had the advantages of being simple and robust to errors in disc sizes but is not the most physical description of

<footnote n="1">¹In GALICS 2.3, major mergers do not destroy discs entirely (Section 2.2).</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Stellar mass–size relation for discs (blue points) and bulges (orange points) in CANDELS (left), GALICS 2.3 without dissipation (centre), and GALICS 2.3 with $f_0 = 0.3$ (right). The effective radii $R_e$ are half-mass radii on a 2D projection (face-on for discs). The blue and red lines show the mean relations in CANDELS (Dimauro et al. 2019) and are the same on all three columns at a given $z$.</figcaption></figure>
star formation in discs. The star formation time-scale in the starburst mode is less important because it has no impact on the statistical properties of galaxies considered in this article as long as it is short compared with the age of the Universe. Here, we use a starburst star-formation time-scale of 0.24 Gyr (consistent with Kennicutt & De Los Reyes 2021) for both central and extended starbursts (i.e. starbursts in cusps and discs, respectively).

In this article, we revert to the original model of Cattaneo et al. (2017) for the quiescent mode and compute the SFR by assuming that:

$$SFR = \epsilon_{sf} \frac{M_{gas}}{t_{dyn}}, \qquad (6)$$

where $M_{gas}$ is the mass of the gas in the disc, $\epsilon_{sf}$ is the efficiency of quiescent star formation, and

$$t_{dyn} = \frac{2\pi r}{v_c(r)} \qquad (7)$$

is the orbital time at some characteristic radius $r$ (in our SAM, the exponential scale-length). The only difference with GALICS 2.0 (Cattaneo et al. 2017) is that GALICS 2.3 does not contain any gas surface density threshold for star formation. That is based on de los Reyes & Kennicutt (2019), who argue that, unlike the Schmidt law $\Sigma_{SFR} \propto \Sigma_{gas}^k$, the Silk–Elmegreen law $\Sigma_{SFR} \propto \Sigma_{gas}/t_{dyn}$ has no clear turnover at low-gas densities ($\Sigma_{SFR}$ and $\Sigma_{gas}$ are the SFR surface density and the gas surface density, respectively). Equation (6) shows that $\epsilon_{sf}$ is degenerate with respect to $M_{gas}$. That is the reason why our predictions for stellar masses and SFRs are remarkably insensitive to small variations of $\epsilon_{sf}$. If we increase $\epsilon_{sf}$, we consume gas more rapidly, but the higher star formation efficiency compensates for the lower mass of the remaining gas, so that the final SFR is the same. The only way to break the degeneracy is to look at the gas content of galaxies. That is only possible in the local Universe because there are no H I data beyond $z \sim 0.1$.

In Cattaneo et al. (2017), we had calibrated $\epsilon_{sf}$ on the data of Boselli, Cortese & Boquien (2014) and found $\epsilon_{sf} \sim 1/25$. The main difficulty with that calibration was that the rotation speeds in Boselli et al. (2014) were measured from H I half line-widths rather than resolved rotation curves. Hence, it is difficult to know at what radius they correspond. A simple solution is to assume a flat rotation curve. We have also considered an alternative calibration based on los Reyes & Kennicutt (2019) but here there is the additional complication that los Reyes & Kennicutt (2019) computed $t_{dyn}$ at the outer radius of the star formation region, i.e. the radius $r$ that contains $\sim 95$ per cent of the H $\alpha$ flux. The star formation radius can be as small as half the exponential scale-length or as large as three exponential scale-lengths. For an average star formation radius of 1.6 exponential scale-lengths, the data of los Reyes & Kennicutt (2019) suggest $\epsilon_{sf} \sim 1/15$. The two calibrations give an average time-scale for quiescent star-formation at $z = 0$ of 3.5 and 2.1 Gyr, respectively. We have run GALICS 2.3 for both $\epsilon_{sf} = 1/25$ and $\epsilon_{sf} = 1/15$. The results are quite similar but slightly better for $\epsilon_{sf} = 1/25$, especially for the H I content of local galaxies (Fig. 2). We therefore retain $\epsilon_{sf} = 1/25$ as our default model, also for consistency with our previous work.

We have verified that these changes in our model for quiescent star formation have no impact on our previous results (Koutsouridou & Cattaneo 2022). The reason is the one discussed above. Moderate changes of the star formation time-scale have almost no effect on stellar masses and SFRs. They are detectable only through their impact on the gas content of galaxies.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Gas-to-stellar mass ratio versus stellar mass in GALICS 2.3 (black point clouds) and the observations (data points with error bars: Garnett 2002; Swaters & Balcells 2002; Noordermeer et al. 2005; Zhang et al. 2009; Boselli et al. 2014; Calette et al. 2021). The data points from Boselli et al. (2014) are higher than those from the other authors because they include the molecular gas. The model shown here includes tidal stripping, but the figure without tidal stripping is almost identical. The black point cloud is for star-forming galaxies (sSFR $> 10^{-11}$ yr$^{-1}$) with $M_{gas} > 10^8$ M$_\odot$. The red line shows the mean relation in GALICS 2.3 with this selection criterion.</figcaption></figure>

## 2.5 Tidal stripping

We model tidal stripping as in Koutsouridou & Cattaneo (2019) and Tollet et al. (2017). The latter had followed a semi-empirical approach based on abundance matching rather than the semi-analytic one presented here. We find the tidal radius $r_t$ of a satellite galaxy by solving the equation:

$$\frac{M_s(r_t)}{r_t^3} = \epsilon_{ts} \frac{M_h(R)}{R^3} \qquad (8)$$

where $M_s(r)$ is the satellite’s mass within radius $r$ from the centre of the satellite, $R$ is the distance of the satellite from the centre of the host, $M_h(R)$ is the mass of the host system within radius $R$ from its centre, and $\epsilon_{ts}$ is the efficiency of tidal stripping. $M_s(r_t)/r_t^3$ is a decreasing function of $r_t$. The higher $\epsilon_{ts}$, the smaller $r_t$. The tidal radius is minimum when the satellite galaxy is at its orbital pericentre within the host system ($R = R_p$).

We have recently discovered that Tollet et al. (2017) had made a mistake of a factor of two in the calculation of $\epsilon_{ts}$ while passing from their equation (16) to their equation (17). The correct formula is:

$$\epsilon_{ts} = \frac{\alpha^2}{-2\phi_s(r_t)/v_c^2(r_t) - 1} \left[ \frac{V_c(R_p)}{V_p} \right]^2, \qquad (9)$$

where $\phi(r)$ is the gravitational potential of the satellite, $v_c(r)$ is the circular velocity of the satellite, $V_c(R)$ is the circular velocity of the host and $\alpha$ is a parameter of the model. Equation (9) was derived assuming that the density distributions of the host and the satellite are power laws of $R$ and $r$, respectively, and that they have the same exponent $\alpha$. We made that simplification so that we could tackle the problem analytically. For the profile of Navarro, Frenk & White (1997; NFW), the exponent varies from $\alpha = -3$ at large radii to $\alpha = -1$ in the inner region. Numerical experiments where the host and the satellite follow NFW profiles are in reasonably good agreement with simple analytic estimates for $\alpha = -3$ (Tollet et al. 2017, appendix A), which, incidentally, corresponds to the classical
<figure><img src="image_4.png" /><figcaption>**Figure 3.** The mass functions of all galaxies (black), star-forming galaxies (blue), and passive galaxies (red). The curves show the predictions of GALICS 2.3. The data points are the observations by Tomczak et al. (2014) at $0.2 < z < 3.0$ (circles) and Weigel, Schawinski & Bruderer (2016) at $0.02 < z < 0.06$ (squares). The grey shaded area in the panel at lowest $z$ shows the possible range for the mass function in the SDSS according to Bernardi et al. (2017). The SAM has a formal mass resolution limit of $M_* \sim 10^7 M_\odot$ (the stellar mass corresponding to the minimum halo mass that we can resolve based on the $M_*-M_{vir}$ relation in Fig. 5).</figcaption></figure>

Jacobi limit. Hence, we choose $\alpha = -3$ as our default value for models with tidal stripping.

Tollet et al. (2017) thought they were computing tidal stripping for $\alpha = -3$ but they were effectively using a lower value (by a factor of $\sqrt{2}$) because of the mistake in their equation (17). The model was in reasonably good agreement with the mass functions of central and satellite galaxies, but, on a closer look, the green curves in fig. 13 of Tollet et al. (2017) overestimated the stellar masses of galaxies by at least 0.1 dex, which is logical, since their stripping was not strong enough.

All figures in this article are for $\alpha = -3$. We have also run our model with $\alpha = 0$ (no tidal stripping) for comparison. We found very little difference at $z > 0.5$. At $z < 0.5$, the model without tidal stripping displays an excess of galaxies with $M_* \sim 10^{12} M_\odot$. They are brightest-cluster galaxies (BCGs) that grow to very large masses by cannibalising satellite galaxies. Tidal stripping transfers some of the stellar mass of satellite galaxies to the intracluster light (ICL) before satellites merge with the BCGs of their host holes. It thereby reduces the masses of BCGs by 0.1 dex on average, in agreement with previous findings by Tollet et al. (2017). Montes (2022) finds the ICL of groups and clusters accounts for 50 to 70 per cent of the total light of the ICL and the BCG combined. Assuming that stellar mass is proportional to luminosity, our finding is consistent with hers if 10 to 25 per cent of the observed ICL light comes from satellite galaxies that no longer exist because they have merged with the BCG (10 and 25 per cent correspond to an ICL fraction of 50 and 70 per cent, respectively). One should note, however, that the fundamental picture is the same with or without tidal stripping, and that none of the results of this article is sensitive to our model for tidal stripping.

## 3 MASS FUNCTIONS

In Koutsouridou & Cattaneo (2022), we had already shown that our SAM reproduces the mass function of galaxies at $0.5 < z < 2.5$ (Ilbert et al. 2013; Muzzin et al. 2013; Tomczak et al. 2014) as well as the local data by Yang, Mo & van den Bosch (2009), Baldry et al. (2012), Bernardi et al. (2013) and Moustakas et al. (2013). The novelty here is that we look at the mass functions of star-forming and passive galaxies, separately. To compare our predictions with their data, we must consider how they separated the two populations in order to adopt a specific star formation rate (sSFR) threshold consistent with theirs.

Tomczak et al. (2014)'s classification is based on Whitaker et al. 2011's redshift-dependent colour-colour (CC) criterion (Fig. 4, left panels). The central column shows the *loci* of red and blue galaxies (according to Whitaker et al. 2011's CC criterion) from a CANDELS sample with SFR measurements on an sSFR-$M_*$ diagram. The sSFR distribution in CANDELS is bimodal. One can clearly distinguish a mainly blue high sSFR population and a mainly red low sSFR population. The intermediate region with lower point density shifts to lower sSFRs at low $z$.

Our approach is to identify an sSFR boundary that separates star-forming and passive galaxies in CANDELS and to apply that boundary to separate star-forming and passive galaxies in our SAM. That approach makes sense if the main sequence of star-forming galaxies occupies the same sSFR range in GALICS 2.3 and CANDELS. If the sSFRs in GALICS 2.3 were systematically lower than those in CANDELS, then a boundary based on CANDELS may run through our SAM's main sequence, so that many galaxies on our SAM's main sequence would be misclassified as passive.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Whitaker et al. (2011)'s CC selection criterion in four bins of redshift (left) and the location of red and blue galaxies from CANDELS (classified according to this criterion) on a SSFR-stellar mass diagram (centre). The right column shows the sSFR–M* in GALICS 2.3 (grey point cloud), which contains many more galaxies than CANDELS at low M*. We have randomly selected a subsample of the grey points so that there is the same number of objects per interval of M* in GALICS 2.3 as in CANDELS and we have coloured this subsample using four colours: blue for unquenched central galaxies, green for unquenched satellite galaxies, light red for quenched central galaxies, and dark red for quenched satellite galaxies. The green dashes show the critical sSFR used to separate star-forming and passive galaxies. The cyan curves show the mean sSFR–M* of Popesso et al. (2023) ±0.45 dex. All lines are the same in the central and the right column to assist the comparison. The green shaded rectangles on the right panels show that, for most galaxies, the transition from unquenched to quenched and from star-forming to passive occurs in the stellar mass range 10¹⁰.² M☉ ≲ M* ≲ 10¹¹.² M☉.</figcaption></figure>
<table><caption>**Table 1.** Critical sSFRs used to separate star-forming and passive galaxies.</caption><thead><tr><th>Redshift</th><th>sSFR<sub>crit</sub> (yr<sup>-1</sup>)</th></tr></thead><tbody><tr><td>z &lt; 0.5</td><td>10<sup>-11</sup></td></tr><tr><td>0.5 &lt; z &lt; 1.0</td><td>10<sup>-10</sup></td></tr><tr><td>1.0 &lt; z &lt; 1.5</td><td>1.8 10<sup>-10</sup></td></tr><tr><td>1.5 &lt; z &lt; 2.0</td><td>2.5 10<sup>-10</sup></td></tr><tr><td>z &gt; 2.0</td><td>5.0 10<sup>-10</sup></td></tr></tbody></table>

Popesso et al. (2023) combined the data from 27 publications and fitted the main sequence with the cyan solid curves. Expectedly, most of the star-forming in the CANDELS sample lie within ±0.45 dex from Popesso et al. (2023)'s main sequence (i.e. within the cyan dashed curves; Fig. 4). If they did not, it would mean that our CANDELS sample is at odds with other observational studies. The right panels shows how galaxies populate the sSFR–M* diagram in GALICS 2.3. The colour of the symbols (discussed later) is not important for our current argument. The cyan lines are the same as in the middle panels. They have been replicated in the right panel to aid the comparison. Reassuringly, the bulk of the star-forming population in GALICS 2.3 lies within the cyan dashed curves, too. That means that we can indeed apply a sSFR criterion derived from CANDELS to separate star-forming and passive galaxies in GALICS 2.3.

We have explored two ways of separating star-forming and quiescent galaxies on an sSFR–M* diagram: horizontal cuts at redshift-dependent critical sSFRs or cuts parallel the ridge line of the main sequence of star-forming galaxies (i.e. to the cyan curves). We found that cuts at the sSFRs in Table 1 (corresponding to the green horizontal dashed lines in Fig. 4) produced the best fit to the observational data (Fig. 3).

A word of caution is that the comparison with the mass function of Tomczak et al. (2014) is necessarily imperfect. The CANDELS data show the presence of red galaxies in the middle of the blue points on the main sequence of star-forming galaxies, so that no sSFR criterion can be entirely successful at reproducing the effects of a colour classification. Table 2 analyses this problem quantitatively. For 86 per cent (at high z) to 90 per cent (at low z) of the objects in the CANDELS sample, a classification based on Whitaker et al. 2011's CC criterion and one based on sSFR return the same outcome. In most cases, blue galaxies are star-forming and red galaxies are passive. If the number of red galaxies incorrectly classified as passive and the number of blue classified as star-forming even though they have sSFR < sSFRcrit were equal, then these errors would cancel each other out and they would have no impact on the number densities of star-forming and passive galaxies. Table 2 does show systematic differences, but the difference between the number of blue galaxies with sSFR < sSFRcrit and the number of red galaxies with sSFR > sSFRcrit (between the third and the second column) constitutes a small fraction of the total galaxy population at any given redshift: 8 per cent at 2 < z < 2.5 and only 2 per cent at 0.5 < z < 1. Hence, the error that we commit by comparing predictions based on a sSFR criterion with observations based on a CC criterion is unlikely to be significant.

The above conclusion is vindicated by comparing the mass functions of Tomczak et al. (2014) at 0.2 < z < 0.5 with those by Weigel et al. (2016) at 0.02 < z < 0.06 (Fig. 3). Weigel et al. (2016) did not use a CC criterion. They had SFR measurements and they separated star-forming and passive galaxies at sSFRcrit = 10⁻¹¹ (yr⁻¹) as we do in our model at z < 0.5. If an sSFR criterion and a CC criterion gave substantially different results, there should have been substantial differences between red/blue squares at 0.02 < z < 0.06 and the red/blue circles at 0.2 < z < 0.5. We see none besides the normal evolution of the mass functions with redshift.

Fig. 3 shows that GALICS 2.3 is in excellent agreement with the observations at intermediate redshifts (0.5 < z < 2.0). Its only weakness is a slight tendency to underestimate the number densities of passive galaxies at intermediate masses (panels at 0.5 < z < 0.75, 0.75 < z < 1, and 1.5 < z < 2). At z > 2, GALICS 2.3 tends to overestimate the number densities of passive galaxies at high masses (M* ≳ 10¹¹ M☉) and to underestimate them at lower masses. The discrepancy at M* ≳ 10¹¹ M☉ could be easily fixed by lowering sSFRcrit at z > 2 but that would come at the price of an ever larger tension with the data at lower masses. At low z, there is the opposite problem. At M* ≳ 10¹¹.⁵ M☉, GALICS 2.3 overpredicts the number densities of both star-forming and passive galaxies but especially of the former. The problem manifests clearly at 0.02 < z < 0.06, where the blue curve ends with a sudden upturn.

Bernardi et al. (2017) have shown that the mass function of galaxies at z ~ 0 is much more uncertain than suggested by a purely statistical analysis because there are also systematic errors from the stellar population model and from the procedure used to fit the light profile. One cannot exclude that the mass function at 0.02 < z < 0.06 may be as high at the upper boundary of the grey shaded area (Fig. 3). In that case, the discrepancy with the data points at M* ≳ 10¹¹.⁵ M☉ would derive from Weigel et al. (2016)'s underestimating the stellar mass function at high M*.

Still, our SAM predicts the existence of star-forming galaxies with M* ~ 10¹² M☉, for which there are no observational counterparts. We have focused on the most massive galaxy in our computational volume, which happens to be star-forming, to investigate how that could be possible. We found that it has an anomalous growth history due to our model for orphan galaxies and ghost subhaloes (Appendix A).

Fig. 5 shows the local stellar–halo mass relation in GALICS 2.3. The blue dots are ‘unquenched’ galaxies that have never satisfied equation (1). The red dots are quenched passive galaxies (galaxies with sSFR < 10⁻¹¹ yr⁻¹). The few green dots (corresponding to a number density of 0.675 × 10⁻³ Mpc⁻³) are quenched galaxies with residual star formation. GALICS 2.3 allows the existence of such galaxies because in our SAM BH blows away all the gas in the central starburst but the gaseous disc is not blown away (Section 1). We have studied these objects in some detail.

Quenched star-forming galaxies have a Gaussian mass distribution with M* = 10¹⁰.⁶±⁰.³ M☉ (Fig. 6, left). Their sSFR distribution has a peak at 10⁻¹⁰.² yr⁻¹ (Fig. 6, centre) corresponding to the ridge line of the main sequence of star-forming galaxies at z = 0 for M* =

<table><caption>**Table 2.** Breakdown of the CANDELS sample by colour and sSFR in five redshift intervals (numbers of galaxies).</caption><thead><tr><th>Redshift</th><th>Blue with sSFR&gt; sSFR<sub>crit</sub></th><th>Red with sSFR&gt; sSFR<sub>crit</sub></th><th>Blue with sSFR&lt; sSFR<sub>crit</sub></th><th>Red with sSFR&lt; sSFR<sub>crit</sub></th></tr></thead><tbody><tr><td>0.5 &lt; z &lt; 1.0</td><td>4384</td><td>252</td><td>387</td><td>1327</td></tr><tr><td>1.0 &lt; z &lt; 1.5</td><td>2852</td><td>110</td><td>325</td><td>537</td></tr><tr><td>1.5 &lt; z &lt; 2.0</td><td>1443</td><td>76</td><td>153</td><td>247</td></tr><tr><td>2.0 &lt; z &lt; 2.5</td><td>597</td><td>24</td><td>86</td><td>64</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** M<sub>⋆</sub>−M<sub>vir</sub> relation in GALICS 2.3 at z = 0.05: scatter plot (blue: unquenched galaxies; green: quenched star-forming galaxies; red: quenched passive galaxies), and median (red curve). The symbols show the data points of Reyes et al. (2012, cyan squares) and Wojtak & Mamon (2013) for late-type (blue triangles) and early-type (red circles) galaxies. To avoid overcrowding, we have shown only a subsample of the model galaxies.</figcaption></figure>

10¹⁰.⁶ M☉ but is clearly skewed towards low sSFRs. A little more than half of the quenched star-forming galaxies are still on the main sequence but many are entering or have entered the green valley. Very few have high sSFRs. They correspond to gas-rich major mergers, in which BH feedback has kicked in and all the gas is star bursting, whether in the cusp or the disc (Section 2.2).

The distribution for the time of quenching *t*_q_ Fig. 6, right) shows that many star-forming galaxies experienced BH quenching several gigayears ago. To appreciate the significance of this finding, we compare *t*_q_ with the time-scale on which a disc runs out of gas after it has ceased to accrete gas from the environment. In the local Universe, a galaxy with sSFR~10⁻¹⁰.² yr⁻¹ becomes passive when its sSFR drops below 10⁻¹¹ yr⁻¹ (Table 2). For a gas fraction of 20 per cent (the gas fraction that corresponds to M* = 10¹⁰.⁶ M☉; Fig. 2), the sSFR drops by 0.8 dex in 1.7 star-formation time-scales, that is, ~6 Gyr for a star-formation time-scale of 3.5 Gyr (Section 2.2). That is an upper limit because, in our SAM, BHs grow because of mergers, and mergers always transfer some gas from the disc to the central cusp, where star formation is more efficient. In fact, in major mergers, all the gas becomes starbursting with a depletion time of 0.24 Gyr (Section 2.2). Hence, the long delay of 6 Gyr between the quasar

phase and the actual shutdown of star formation is relevant only to galaxies that were quenched by a minor merger and that have not experienced any major merger thereafter.

Galaxies that were quenched more than 6 Gyr ago (at $z > 0.6$) cannot have kept making stars for such a long time without an external gas supply. As we do not allow gas accretion on to quenched galaxies (Section 1), mergers are the only possible supply mechanism in our SAM. Star-forming galaxies with $t_q \lesssim 10$ Gyr are, in their majority, galaxies where mergers with gas-rich satellites have temporarily reactivated star formation.

Fig. 5 shows that quenched passive galaxies dominate the M*−M*vir relation at high masses (the most massive galaxy clearly stands out as an outlier). The red curve shows the median M* as a function of M*vir in our SAM. GALICS 2.3 is in reasonably good agreement with the lensing data of Reyes et al. (2012; cyan squares) for local disc galaxies and with Wojtak & Mamon (2013)'s mass estimates from satellite kinematics for late-type (blue triangles) and late-type (red circles) galaxies.

GALICS 2.3 reproduces the upturn of the mass functions of passive galaxies at low M*, which is clearly visible in the data points at 0.25 < z < 0.75. This feature was already present in the luminosity function of early-type galaxies from the Two-degree Field Galaxy Redshift Survey (Folkes et al. 1999), which was close to a Gaussian at high luminosities but rose again at low luminosities. More than 20 years ago, Cattaneo (2001) had remarked that merger remnants account for the Gaussian part at high luminosities but cannot explain the low-mass early-type population, for which one needs a different explanation.

Ilbert et al. (2010), Tomczak et al. (2014), and Huertas-Company et al. (2015) have later confirmed the presence of two early-type populations: a massive bulge-dominated population with a Gaussian mass function centred on $M_* = 10^{10.7\pm0.5}$ M☉, which already existed at $z \sim 2-3$, and a lower mass population that began to emerge much later. Many of the galaxies in this second population are red discs, which one may visually classify as lenticular because of the lack of spiral arms (Huertas-Company et al. 2015).

In GALICS 2.3, these populations correspond to two pathways from the star-forming to the passive population: BH quenching and environmental strangulation. To show the importance of these two mechanisms, we have coloured the model galaxies in the right panels of Fig. 4 according to both quenching status and their status as central or satellite galaxies: we have used blue for unquenched central galaxies, green for unquenched satellite galaxies, light red for quenched central, and dark red for quenched satellite galaxies. Many galaxies are shown in grey, although they must fall in one of the four categories, because, to have the same numbers of coloured

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Distribution of stellar mass M* (left), sSFR (centre), and quenching time t_q for quenched star-forming galaxies at z = 0.05; t_q is the cosmic time (time since the big bang) at which quenching occurs. N is the number of galaxies in each bin (out of 675).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Σ₁−M* relation in GALICS 2.3 at z ~ 0.05 (left) and z = 0.64 (right). The small circles are galaxies that have not been quenched (based on the quenching criterion in equation 1). The large circles are quenched passive galaxies. The stars are galaxies that have been quenched but have not become passive yet because of residual star formation in the disc component. Colour shows the bulge-to-total mass ratio *B/T*. The red line at z = 0.05 corresponds to the ridge line of quiescent galaxies in ZENS (Tacchella et al. 2017). The green line at z = 0.64 is Chen et al. (2020)'s quenching boundary at that redshift.</figcaption></figure>

symbols for the model and the observations, we have coloured only a subset of all model galaxies.

In the first pathway (BH quenching), galaxies grow along the main sequence (cyan curves) until they reach the mass range where quenching occurs (10$^{10.2}$ M$_{\odot}$ < M$_{*}$ < 10$^{11.2}$ M$_{\odot}$; Ilbert et al. 2010; Huang et al. 2013; Barro et al. 2017). That is where the red symbols progressively replace the blue and green ones. Once quenched, galaxies drop out of the main sequence through the green rectangles, which correspond to the mass range above. This is the route through which the high-mass passive population was formed. This population is mainly composed of central galaxies (light red) although it also contains satellites galaxies (dark red).

The second pathway (environmental effects such as ram pressure and tides) is more important at $M_{*} < 10^{11.2}$ M$_{\odot}$, where most of the galaxies with $sSFR < sSFR_{crit}$ are unquenched satellites (green circles). In GALICS 2.3, satellite galaxies fall systematically below the main sequence of star-forming galaxies. At $z > 1.5$, most of the green circles lie between the lower cyan dashed curve and the green horizontal dashed line (below the lower boundary of the main sequence but still within the star-forming population). At $z < 1.5$, however, a growing number of satellite galaxies start slipping below the green dashes, which are now running through the middle of the green symbols. Interestingly, in CANDELS, too, there seem to be two populations of blue galaxies, one on the main sequence and one just below it (blue circles in the middle panels of Fig. 4; the presence of two populations is more obvious at $z > 1.5$, where there are fewer data points). We have no information on the environment of these galaxies but we speculate that the blue CANDELS population with lower sSFRs may be the population of satellite galaxies, many of which are still making stars.

The finding that passive galaxies of low mass ceased to make stars because of environmental effects is consistent with chemical evidence that in low-mass galaxies star formation shut down slowly over several gigayears (Peng, Maiolino & Cochrane 2015). Also see Goubert et al. 2024 for evidence that quiescence in low-mass satellites correlates with the environment more than it does with BH mass.

Two considerations suggest that GALICS 2.3 may not fully account for the importance of this second route. One is the deficit of passive

galaxies in our SAM at $M_{*} \sim 10^{10}$ M$_{\odot}$ (Fig. 3; most prominent at 0.2 < z < 0.5), where environmental strangulation should be the main pathway to quiescence. The other is the excess of star-forming galaxies at high masses (most prominent at 0.02 < z < 0.06). Most of these massive star-forming galaxies were quenched by BH feedback but have been replenished with gas by mergers with gas-rich satellite (Figs 5 and 6). Perhaps the reason there are too many of them is that, in our SAM, ram pressure and tidal stripping are not as effective as they should be.

## 4 SURFACE DENSITIES

In both CANDELS and GALICS 2.3, discs and bulges follow different mass–size relations. For a same $M_{*}$, bulges have smaller effective radii than discs (Fig. 1) and thus higher surface densities. Fig. 7 shows that, in our SAM, galaxies with higher bulge-to-total stellar mass ratios *B/T* have statistically higher Σ₁. Unquenched galaxies (small circles) follow a different, steeper relation than quenched ones (large symbols).

Fig. 7 at $z = 0.05$ is analogous to fig. 2 of Tacchella et al. (2017), who plotted Σ₁ versus $M_{*}$ for galaxies in the Zurich ENVironmental Study (ZENS, 0.05 < z < 0.0585) and found a steeper relation for star-forming galaxies. The red line is the Σ₁−M* relation for passive galaxies in ZENS and the SDSS (Fang et al. 2013). Fig. 7 shows that at $z = 0.05$ quenched galaxies in GALICS 2.3 follow the SDSS/ZENS relation for passive galaxies. The transition between the steeper relation for star-forming and the shallower relation for quenched galaxies occurs at the same mass scale as in Tacchella et al. (2017).

The difference, which we have already seen in Fig. 5, is that in GALICS 2.3 not all quenched galaxies are passive. In Fig. 7, we have used stars to denote quenched galaxies with residual star formation and distinguish them from passive quenched galaxies (large circles). Once again, most quenched galaxies are passive, but Fig. 7 does show seven star-forming galaxies on the Σ₁−M* relation for passive galaxies at $M_{*} > 10^{11}$ M$_{\odot}$. In Tacchella et al. (2017), there is only one despite the fact that our Fig. 5 and fig. 2 of Tacchella et al. (2017) contain a similar amount of objects at $M_{*} > 10^{11}$ M$_{\odot}$. With the exception of the most massive galaxy, which is anomalous and the
<figure><img src="image_9.png" /><figcaption>**Figure 8.** $\Sigma_1-M_\star$ relation in GALICS 2.3 at $0.5 < z < 2.5$. Galaxies are shown with blue or red circles according to whether they are star-forming or quiescent, respectively. Large circles correspond to quenched galaxies, small circles to galaxies that have not been quenched. The green lines are the quenching boundaries of Chen et al. (2020) in the mass range where they have been established observationally. The population that ceased to make stars because of environmental effects (small red circles) does not appear on this figure because we have shown much fewer galaxies here than in Fig. 4. If we had shown the same number of objects, the red sequence would be saturated by large circles.</figcaption></figure>

only unquenched galaxy on the diagram at $M_\star > 10^{11} M_\odot$, the origin of the problem is that the discs of quenched galaxies keep forming stars for too long. Quenched ellipticals are all passive. The fraction of

quenched galaxies with residual star formation is larger when there is a significant disc component. In both ZENS and GALICS 2.3, almost all star-forming galaxies lie below the red line.
Tacchella et al. (2017) coloured their data points by $B-I$. We have coloured our model galaxies by $B/T$. Their similarity is due to the colour–morphology relation. Galaxies become redder as their morphologies evolves towards earlier types. Morphological evolution drives $\Sigma_1$ up. When $\Sigma_1$ crosses the quenching boundary, they leave the main sequence of star-forming galaxies and migrate to the red sequence of passive galaxies. In ZENS, passive galaxies above the ridge line of the red sequence are systematically redder than those below it. In their vast majorities, passive galaxies with $B-I \gtrsim 1.4$ are above the red line and passive galaxies with $B-I \lesssim 1.2$ are below it. Fig. 7 shows an analogous behaviour with respect to morphology. At a given $M_*$, the average $B/T$ increases with $\Sigma_1$. Most elliptical galaxies ($B/T > 0.8$) are above the red line and they all correspond to quenched galaxies (large symbols). Most of the galaxies in which star formation has been quenched are objects with $B/T > 0.5$. Quenched galaxies with $0.1 < B/T < 0.3$ do exist (large blue symbols). Their $\Sigma_1$ are lower than those of quenched galaxies with higher $B/T$, in the same way that not-so-red quiescent galaxies in ZENS have lower $\Sigma_1$ than redder objects. Nearly all disc-dominated galaxies ($B/T < 0.5$) are below the red line, independently of their quenching status.

In Fig. 8, we extend our analysis of the $\Sigma_1-M_*$ relation to the six redshift bins in Fig. 1. Galaxies are classified as star-forming or passive and shown in blue or red based on sSFR (Table 1). The green lines are the quenching boundaries of Chen et al. (2020). We have shown them in green as a visual reminder that they correspond to a sort of ‘green valley’ and that they should not be confused with the ridge line of quiescent galaxies in Fig. 7, which we have shown in red because it corresponds to the red sequence.

The predictions of GALICS 2.3 are in good agreement with the CANDELS data. The green lines mark the lower boundary of the passive population. In Chen et al. (2020), almost all quiescent galaxies lie above the quenching boundaries. In Fig. 8, almost all the red circles are above the green lines. In Chen et al. (2020), most of the star-forming galaxies lie below the quenching boundaries, although there is a minority that lies above it. The same applies to our model.

Quenched galaxies with residual star formation (the large blue circles above the green lines) have lower $\Sigma_1$ than the red circles (passive galaxies). That means that $\Sigma_1$ increases as the star formation fades away. The conversion of gas into stars is not the main cause of the increase of $\Sigma_1$. $\Sigma_1$ increases because of morphological evolution. As we had already seen in Koutsouridou & Cattaneo 2022, fig. 10, panel B2, $B/T$ increases with decreasing SFR not only as galaxies migrate from the star forming to the passive population (i.e. to $sSFR < 10^{-11} \text{ yr}^{-1}$ at $z \sim 0$) but also all the way down to $sSFR = 10^{-13} \text{ yr}^{-1}$. This morphological evolution is due to further mergers between the one that triggered quenching and the final complete shutdown of star formation. The quenched population is mainly composed of massive central galaxies (Fig. 4). Mergers play a major role in the mass assembly of these systems (e.g. Cattaneo et al. 2011). While it is certainly possible and even probable that the delay between quenching and quiescence may be too long in our SAM, we consider this feature to be robust. Galaxies evolve towards earlier morphological types (higher $B/T$ and $\Sigma_1$) as their stellar populations age after the episode that has triggered quenching.

## 5 SUMMARY AND CONCLUSION

We have considered a simple model, in which major and minor mergers funnel gas into the central regions of galaxies. Supermassive BHs feed on this gas until the energy deposited into the ambient gas

is enough to unbind it or to heat it to such high entropy that its accretion on to galaxies is permanently shut down (‘quenching’). In Koutsouridou & Cattaneo (2022), we had focused on the comparison with the local data. In this article, we have shown that our SAM:

(1) is in good agreement with the evolution of the mass functions of star-forming and passive galaxies with redshift (especially at $0.5 < z < 2$, Fig. 3), and

(2) reproduces Chen et al. (2020)’s quenching boundary on the $\Sigma_1-M_*$ diagram (Fig. 8).

The second finding is noteworthy because, although Chen et al. (2020)’s empirical model had provided the initial motivation for our research, nowhere do the structural properties of galaxies enter our quenching criterion (equation 1). Neither were any free parameters tuned to obtain the result in Fig. 8, which is thus a genuine prediction of our model.

Our SAM’s ability to reproduce the observations bolsters our confidence that we can use it to interpret them. Our conclusion is that $\Sigma_1$ is a morphological indicator and that the quenching boundary arises because the growth of supermassive BHs is linked to the formation of bulges. For a given stellar mass, galaxies with a larger $B/T$ have a larger BH mass.

In our interpretation, the quenching boundary on the $\Sigma_1-M_*$ diagram is simply an aspect of the colour–morphology relation. The $\Sigma_1-M_*$ relation is stratified in both colour/stellar age (Tacchella et al. 2017; Luo et al. 2020) and $B/T$ (Luo et al. 2020 and Fig. 7 of this article). Spiral galaxies (low $\Sigma_1$) are mainly star-forming. Elliptical galaxies (high $\Sigma_1$) are almost all passive. The transition from late types to early types is both morphological and spectral. The $\Sigma_1$ quenching boundary corresponds to the critical $B/T$ that separates the star-forming spiral population from the passive S0/elliptical population. In our SAM, the critical $B/T$ is in the range $0.3 < B/T < 0.5$ at $z \lesssim 1.3$ (Fig. 7, right panel) and slightly higher at earlier epochs, consistently with Chen et al. (2020)’s finding that the quenching threshold increases with $z$.

Fig. 7 of this article and fig. 10 of Koutsouridou & Cattaneo (2022) show that $B/T$ keeps increasing after the event that triggered quenching because of further mergers while the SFR fades away. This finding is important to explain why the $\Sigma_1-M_*$ relation displays colour stratification even within the passive population and is in agreement with Cattaneo et al. (2011), who estimated that galaxies with $M_* > 10^{11} M_\odot$ accreted most of their mass through dissipationless mergers.

BH quenching explains the high-mass passive population. The mass function of this population has a maximum at $M_* \sim 10^{10.7} M_\odot$ (Fig. 3). The mass function of passive galaxies rises again towards low masses at $M_* \lesssim 10^9 M_\odot$. This low-mass passive population is almost entirely composed of satellite galaxies and is the product of environmental effects such as strangulation, ram-pressure stripping and tidal stripping. This finding is consistent with Huertas-Company et al. (2015)’s morphological analysis (the lower mass passive population is dominated by red discs/lenticulars rather than elliptical galaxies) and Peng et al. (2015)’s chemical analysis (star formation has depleted the gas reservoir over a few gigayears after accretion from the environment has stopped).

## DATA AVAILABILITY

The results of these articles do not depend on any non-public data.

## REFERENCES

Baldry I. K. et al., 2012, MNRAS, 421, 621

Barro G. et al., 2017, ApJ, 840, 47
Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013a, ApJ, 762, 109

Behroozi P. S., Wechsler R. H., Wu H.-Y., Busha M. T., Klypin A. A., Primack
J. R., 2013b, ApJ, 763, 18

Bernardi M., Meert A., Sheth R. K., Fischer J. L., Huertas-Company M., Maraston C., Shankar F., Vikram V., 2017, MNRAS, 467, 2217

Bernardi M., Meert A., Sheth R. K., Vikram V., Huertas-Company M., Mei S., Shankar F., 2013, MNRAS, 436, 697

Boselli A., Cortese L., Boquien M., 2014, A&A, 564, A65

Calette A. R., Rodríguez-Puebla A., Avila-Reese V., Lagos C. D. P., 2021, MNRAS, 506, 1507

Cattaneo A. et al., 2017, MNRAS, 471, 1401

Cattaneo A., 2001, MNRAS, 324, 128

Cattaneo A., Koutsouridou I., Tollet E., Devriendt J., Dubois Y., 2020, MNRAS, 497, 279

Cattaneo A., Mamon G. A., Warnick K., Knebe A., 2011, A&A, 533, A5

Cattaneo A., Teyssier R., 2007, MNRAS, 376, 1547

Chen Z. et al., 2020, ApJ, 897, 102

Chilingarian I. V., Di Matteo P., Combes F., Melchior A. L., Semelin B., 2010, A&A, 518, A61

de los Reyes M. A. C., Kennicutt R. C. J., 2019, ApJ, 872, 16

de Vaucouleurs G., 1948, L'Astronomie, 62, 344

Devergne T. et al., 2020, A&A, 644, A56

Diemer B., Behroozi P., Mansfield P., 2023, preprint (arXiv:2305.00993)

Dimauro P. et al., 2019, MNRAS, 489, 4135

Faber S. M. et al., 1997, AJ, 114, 1771

Fang J. J., Faber S. M., Koo D. C., Dekel A., 2013, ApJ, 776, 63

Folkes S. et al., 1999, MNRAS, 308, 459

Garnett D. R., 2002, ApJ, 581, 1019

Goubert P. H., Bluck A. F. L., Piotrowska J. M., Maiolino R., 2024, MNRAS, 528, 4891

Hernquist L., 1990, ApJ, 356, 359

Hopkins P. F., Cox T. J., Dutta S. N., Hernquist L., Kormendy J., Lauer T. R., 2009a, ApJS, 181, 135

Hopkins P. F., Lauer T. R., Cox T. J., Hernquist L., Kormendy J., 2009b, ApJ, 181, 486

Huang J. S. et al., 2013, ApJ, 766, 21

Huertas-Company M. et al., 2015, ApJ, 809, 95

Ilbert O. et al., 2010, ApJ, 709, 644

Ilbert O. et al., 2013, A&A, 556, A55

Kannan R., Macciò A. V., Fontanot F., Moster B. P., Karman W., Somerville R. S., 2015, MNRAS, 452, 4347

Kennicutt R. C. J., De Los Reyes M. A. C., 2021, ApJ, 908, 61

Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340

Kormendy J., Fisher D. B., Cornell M. E., Bender R., 2009a, ApJS, 182, 216

Kormendy J., Fisher D. B., Cornell M. E., Bender R., 2009b, ApJS, 182, 216

Koutsouridou I., Cattaneo A., 2019, MNRAS, 490, 5375

Koutsouridou I., Cattaneo A., 2022, MNRAS, 516, 4194

Luo Y. et al., 2020, MNRAS, 493, 1686

Montes M., 2022, Nat. Astron., 6, 308

Moustakas J., Coil A. L., Aird J., Blanton M. R., Cool R. J., Eisenstein D. J., Mendez A. J., Wong K. C., 2013, ApJ, 767, 50

Muzzin A. et al., 2013, ApJ, 777, 18

Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493

Noordermeer E., van der Hulst J. M., Sancisi R., Swaters R. A., van Albada
T. S., 2005, A&A, 442, 137

Peng Y., Maiolino R., Cochrane R., 2015, Nature, 521, 192

Popesso P. et al., 2023, MNRAS, 519, 1526

Reyes R., Mandelbaum R., Gunn J. E., Nakajima R., Seljak U., Hirata C. M., 2012, MNRAS, 425, 2610

Sahu N., Graham A. W., Davis B. L., 2022, ApJ, 927, 67

Sérsic J. L., 1963, Bol. Asociacion Argentina Astron. Plata Argentina, 6, 41

Sheth R. K., Mo H. J., Tormen G., 2001, MNRAS, 323, 1

Swaters R. A., Balcells M., 2002, A&A, 390, 863

Tacchella S., Carollo C. M., Faber S. M., Cibinel A., Dekel A., Koo D. C., Renzini A., Woo J., 2017, ApJ, 844, L1

Terrazas B. A., Bell E. F., Henriques B. M. B., White S. D. M., Cattaneo A., Woo J., 2016, ApJ, 830, L12

Thomas D., Maraston C., Bender R., Mendes de Oliveira C., 2005, ApJ, 621, 673

Tollet E., Cattaneo A., Macciò A. V., Kang X., 2022, MNRAS, 515, 3453

Tollet É., Cattaneo A., Mamon G. A., Moutard T., van den Bosch F. C., 2017, MNRAS, 471, 4170

Tomczak A. R. et al., 2014, ApJ, 783, 85

Tweed D., Devriendt J., Blaizot J., Colombi S., Slyz A., 2009, A&A, 506, 647

Weigel A. K., Schawinski K., Bruderer C., 2016, MNRAS, 459, 2150

Whitaker K. E. et al., 2011, ApJ, 735, 86

Wojtak R., Mamon G. A., 2013, MNRAS, 428, 2407

Yang X., Mo H. J., van den Bosch F. C., 2009, ApJ, 695, 900

Zhang W., Li C., Kauffmann G., Zou H., Catinella B., Shen S., Guo Q., Chang R., 2009, MNRAS, 397, 1243

APPENDIX: MASS RESOLUTION OF THE MERGER TREES

The *N*-body simulation used to construct the merger trees employed 1024³ particles for a volume of (100 Mpc)³. For our cosmology that corresponds to an *N*-body particle mass of 3.7 × 10⁷ M☉. The halo finder (HALOMAKER; Tweed et al. 2009) detects haloes with more than a hundred particles. Previous comparison of the halo mass functions from *N*-body simulations with the analytical fit by Sheth, Mo & Tormen (2001) have shown that our halo catalogues are complete above a mass limit of a few hundred particles (Fig. 1 of Cattaneo et al. 2017). In our case, that corresponds to a halo mass of *M*vir ~ 10¹⁰ M☉ and a stellar mass of *M*⋆ ~ 10⁷ M☉ based on stellar-halo mass relation in Fig. 5 of this article.

Subhaloes are much more difficult to resolve, especially when they
are inside massive systems (i.e. groups and clusters of galaxies). A
subhalo that the halo finder can no longer resolve is deemed to have
merged with its host system. That leads to the numerical overmerging
of satellite galaxies with the central dominant galaxies of their host
haloes.

To obviate this problem, we complete the HALOMAKER catalogues with ghost subhaloes, i.e. subhaloes that the halo finder no longer detects but that should still exist based on the model of Tollet et al. (2017), where we showed that a halo catalogue constructed from simulation with 512³ particles and completed with ghosts is a good as one constructed from a simulation with 1024³ particles without ghosts in terms of its ability to reproduce the conditional mass functions of satellite galaxies in groups and clusters.

To test how numerical resolution affects our results, we have rerun GALICS 2.3 on the merger trees of Cattaneo et al. (2017) and used the results to recompute the mass functions of star-forming and passive galaxies (Fig. A1). Cattaneo et al. (2017) used a simulation with the same cosmology, the same volume and the same initial condition as ours but with only employed 512³ particles and no ghost subhaloes. In Fig. A1, the mass resolution is eight times lower than in Fig. 3 for haloes and far worse for subhaloes. Unsurprisingly, there are differences, which we are about to discuss, but they are not large.

At low masses, the black curves are slightly lower in Fig. A1 than
in Fig. 3. For the red curves, the difference is much more pronounced.
Expectedly, at lower resolution there are fewer dwarf galaxies and
especially fewer passive dwarves, most of which are satellites and
therefore poorly resolved.

At high masses, the low-resolution simulation gives a much higher
fraction of star-forming galaxies. In Fig. A1, the red and blue curves
for 0.02 < z < 0.06 cross at M* ~ 10¹¹.² M☉. In the observations,
the fraction of star-forming galaxies at M* ~ 10¹¹.² M☉ is ~10
rather than ~50 per cent. At high resolution (Fig. 3), the fraction of
star-forming galaxies is still on the high side (around 20 per cent)
<figure><img src="image_10.png" /><figcaption>**Figure A1.** As Fig. 3 but for merger trees: (1) from a simulation with 512³ instead of 1024³ particles and (2) without ghost subhaloes. The model shown here has a stellar mass resolution of $M_* \sim 10^8 M_\odot$ for central galaxies. For satellite galaxies, the mass resolution is worse by at least an order of magnitude.</figcaption></figure>

but the discrepancy with the observations is much smaller. The trend is reassuring and let us presume that the agreement with the observations would further improve if we could use merger trees from an *N*-body simulation with even higher resolution.

We have run several tests to investigate why higher resolution increases the fraction of passive galaxies. At low resolution, galaxies of the same mass have lower *B/T* ratios but higher *ex situ* fractions (stellar mass fractions acquired through mergers). That can happen only if the additional *ex situ* stars come from minor mergers, which are less effective at increasing *B/T* and growing supermassive BHs but can supply gas to the discs of massive galaxies.

Ram pressure and tidal stripping remove gas from satellite galaxies and thus limit the importance of this route to replenish quenched massive central galaxies with gas. However, both ram pressure and tidal stripping are more effective when satellite galaxies are closer to their orbital pericentres. At low resolution many subhaloes merge with their hosts long before reaching their pericentres. We have checked that, for a same set of DM merger trees, *B/T* ratios decrease and *ex situ* fractions increase when these processes are turned off.

There are two reasons why minor mergers become more important at low resolution. Both are connected to the fact that low resolution leads to overmerging (see above) and that degrading the resolution has a greater impact on galaxies with lower masses. First, degrading the resolution increases mergers rates at high redshift, where halo masses are smaller, and that explains why the number densities of massive galaxies at $2 < z < 2.5$ are higher at low resolution (compare Figs 3 and A1). Once a galaxy has grown too much, it

becomes hard for it to find companions of comparable mass later on. Hence, minor mergers become more probable. Second, at all redshifts, the merging rates of the galaxies with lowest masses receive the greatest boost and that contributes to make minor mergers comparatively more important.

The most massive galaxy in Fig. 5 is unquenched because it was able to reach a very large mass without any major merger. We have conducted several tests to investigate how that could possibly happen. We have discovered that this system has a pathological growth due to the algorithm that creates ghost subhaloes. The goal is to prevent premature merging but for this object some mergers are excessively delayed.

To investigate the extent to which pathologies of the merger trees may affect our results, we have rerun GALICS 2.3 on DM merger trees constructed using a different algorithm (Behroozi et al. 2013b; CONSISTENT TREES), a different halo finder (Behroozi, Wechsler & Wu 2013a; ROCKSTAR) and a different *N*-body simulation (Klypin et al. 2016; SMDPL). Clearly, GALICS 2.3 perform better with its native merger trees, on which it was calibrated, and other merger trees have issues of their own (see Diemer, Behroozi & Mansfield 2023). However, the results were not substantially different. That bolsters our confidence that the resolution of the *N*-body simulation and the algorithms used to extract merger trees from the simulation should not be a major cause of concern.

This paper has been typeset from a *TEX/LATEX* file prepared by the author.