The First Quenched Galaxies: When and How?

Many quiescent galaxies discovered in the early Universe by JWST raise fundamental questions on when and how these galaxies became and stayed quenched. Making use of the latest version of the semianalytic model GAEA that provides good agreement with the observed quenched fractions up to z ∼ 3, we make predictions for the expected fractions of quiescent galaxies up to z ∼ 7 and analyze the main quenching mechanism. We find that in a simulated box of 685 Mpc on a side, the first quenched massive (M ⋆ ∼ 1011 M ⊙), Milky Way–mass, and low-mass (M ⋆ ∼ 109.5 M ⊙) galaxies appear at z ∼ 4.5, z ∼ 6.2, and before z = 7, respectively. Most quenched galaxies identified at early redshifts remain quenched for more than 1 Gyr. Independently of galaxy stellar mass, the dominant quenching mechanism at high redshift is accretion disk feedback (quasar winds) from a central massive black hole, which is triggered by mergers in massive and Milky Way–mass galaxies and by disk instabilities in low-mass galaxies. Environmental stripping becomes increasingly more important at lower redshift.


Introduction
The cessation of star formation in galaxies has drawn considerable attention in recent years, especially given the large numbers of quiescent massive galaxies that have been found in the early Universe (Glazebrook et al. 2017;Schreiber et al. 2018;Girelli et al. 2019;Merlin et al. 2019;Nanayakkara et al. 2024) when the timescale available to assemble and quench these systems is short.Spectroscopic studies suggest that these galaxies experience short periods of intense star formation, grow up to a stellar mass of 10 11 M e in the first 1 or 2 billion yr of the Universe, and then stop forming stars within a few tens of million years (Forrest et al. 2020;Valentino et al. 2020;Carnall et al. 2023a;Kakimoto et al. 2023).This rapid assembly and quenching process might challenge our current understanding of galaxy formation (Finkelstein et al. 2023).
The number densities of quenched massive galaxies M å > 10 10.5 M e increase rapidly from ∼10 −6 Mpc −3 at z ∼ 5 to as much as a factor of 10 times higher at z ∼ 3 (Marsan et al. 2022) although the measured number densities have a relatively large scatter due to different selection criteria and cosmic variance (Valentino et al. 2023).The classical UVJ color selection of quenched galaxies (Wuyts et al. 2008) is found to be incomplete and underestimates the number of quenched galaxies at z > 3 (Schreiber et al. 2018).Some studies propose a modified UVJ selection (Belli et al. 2019), others favor an NUVrJ (Ilbert et al. 2013) or ugi color selection (Antwi-Danso et al. 2023) to identify galaxies that have been quenched recently, which is important for galaxies at z > 3 (Gould et al. 2023;Kubo et al. 2024).
Despite a large scatter in observational measurements, it is a solid conclusion that most theoretical models underpredict the number densities of quenched galaxies (Cecchi et al. 2019;Girelli et al. 2019;Gould et al. 2023) at z > 4 by about an order of magnitude.Weaver et al. (2023) found that quenched galaxies in the SHARK model (Lagos et al. 2023) and IllustrisTNG simulation (Pillepich et al. 2018) at 3.5 < z < 4.5 are not as massive as the observed ones.Either creating enough massive galaxies at early cosmic epochs or quenching them on a short timescale remains a challenge for current galaxyformation models.
Various physical mechanisms have been proposed to explain the rapid assembly of mass in the early Universe, including weaker feedback (Dekel et al. 2023), enhanced star formation efficiencies (Wang et al. 2023), and a top-heavy initial mass function (IMF; Trinca et al. 2023).The physical mechanisms driving quenching also remain unclear.Quenching could be caused by internal feedback from active galactic nuclei (AGN) and supernova (SN) feedback or by external physical processes including galaxy-galaxy interactions and environmental stripping.With a minimal physical model, Gelli et al. (2023) suggests that SN feedback is not powerful enough to quench galaxies of ∼10 8 M e at high redshift.The fact that many high-z quenched galaxies are found to host luminous AGN (Ito et al. 2022;Carnall et al. 2023b;Belli et al. 2023;D'Eugenio et al. 2023;Davies et al. 2024;Shimakawa et al. 2024) suggests an important contribution to quenching from feedback from their central supermassive black holes (SMBHs).This appears to be confirmed in recent theoretical works: Kurinchi-Vendhan et al.
(2023) and Kimmig et al. (2023) analyze the massive quenched galaxies in IllustrisTNG and Magneticum at z ∼ 3 and show that these galaxies are indeed quenched by AGN feedback.Lovell et al. (2023) found that AGN feedback is the dominant quenching mechanism for galaxies above 10 9 M e at z ∼ 5. Qin et al. (2017) use semianalytic models to identify analogs of quenched galaxies observed at z ∼ 5 and show that these have grown through significant mergers and host the most massive black holes at their redshifts.Some quenched galaxies are found in overdense environments (Kubo et al. 2021;McConachie et al. 2022;Alberts et al. 2023;Tanaka et al. 2023), suggesting environmental quenching may also contribute as early as z ∼ 4.
In our recent work (De Lucia et al. 2024), we present the latest version of our Galaxy Evolution and Assembly (GAEA) model (hereafter GAEA2023) and show that it can correctly reproduce the observed quenched fractions up to redshift ∼3, as well as the number densities of quenched galaxies up to redshift ∼5, better than many state-of-the-art models and simulations.The good agreement with observations makes it a perfect tool for studying the physical origin of quenched galaxies at high redshift.
In this work, we extend the analysis presented in De Lucia et al. (2024) to the quenched fractions and their quenching mechanisms since z ∼ 7.In Section 2, we introduce the semianalytic model and our methodology.In Sections 3 and 4, we present the results and give our conclusions.

Model and Simulation
GAEA2023 (De Lucia et al. 2014) now combines independent versions of the model including an improved treatment for the SN feedback (Hirschmann et al. 2016) of the multiphase cold gas (Xie et al. 2017), environmental effects (Xie et al. 2020), and AGN accretion and feedback (Fontanot et al. 2020).In particular, GAEA2023 implements treatment for tidal stripping and ram pressure stripping that gradually removes hot gas, as well as ram pressure stripping of cold gas for satellite galaxies.These implementations help to solve the overquenching of low-mass satellite galaxies and to improve the model predictions in terms of gas fractions (Xie et al. 2020).GAEA2023 also implements updated modeling for cold-gas accretion onto SMBHs.Mergers and disk instability cause a fraction of cold gas to lose angular momentum and flow toward the center, where it forms an accretion disk around the SMBH.This material is then accreted onto the SMBH on a viscous timescale: accretion can heat the surrounding gas and cause an outflow (for details about the model, we refer to Fontanot et al. 2020).Below, we refer to this process as accretion disk feedback.GAEA2023 has been calibrated to reproduce the galaxy stellar mass function up to z ∼ 3, H I mass function and quenched fractions at z ∼ 0, as well as AGN luminosity function up to z ∼ 4. In De Lucia et al. (2024) we show that this model version reproduces well the observed quenched fraction, the stellar mass function of the quenched population up to z ∼ 3, as well as the number densities of quenched massive galaxies up to z ∼ 5.
The model is run on the Millennium Simulation (Springel et al. 2005) with a box size of 685 Mpc based on a Wilkinson Microwave Anisotropy Probe 1 yr cosmology (Spergel et al. 2003) with Ω m = 0.25, σ b = 0.045, σ 8 = 0.9, and h = 0.73.
In the following, we will compare GAEA2023 results with predictions from TNG100 and TNG300 of the IllustrisTNG project (Marinacci et al. 2018;Naiman et al. 2018;Nelson et al. 2018;Pillepich et al. 2018;Springel et al. 2018).The TNG project is a suite of cosmological magnetohydrodynamical simulations, adopting the Planck cosmology (Planck Collaboration et al. 2016) with Ω m = 0.3089, Ω b = 0.0486, σ 8 = 0.8159, and h = 0.6774.The TNG100 and TNG300 simulate cubic boxes of side lengths approximately 100 and 300 Mpc.TNG considers two modes of AGN feedback: for high accretion rates, the surrounding gas is heated by thermal feedback from AGN.When the accretion rates are low, gas instead receives a kinetic "kick" that causes gas outflows.In this work, we use the publicly available database 16 to retrieve the simulated stellar mass, instantaneous star formation rate, and black hole mass within twice the stellar half-mass-radius.

Quenched Fraction
Figure 1 shows the evolution of quenched fractions as predicted by GAEA2023 for low-mass (2-4 × 10 9 M e ), Milky Way-mass (2-4 × 10 10 M e ), and massive galaxies (0.8-1.5 × 10 11 M e ).We consider different definitions for quenched galaxies: first, we select a sample including specific star formation rate sSFR = SFR/M å < 0.3/t H ,17 where t H is the Hubble time at given redshift (Franx et al. 2008).It is important to stress that the star formation rates from GAEA2023 are averaged over the time interval between two snapshots, which is ∼80 Myr at z ∼ 7 and increases to ∼300 Myr at z ∼ 0. Therefore, our sample of model quenched galaxies does not include those that are only temporarily quenched by, e.g., bursty periods of star formation.To get a fair comparison with observational data, we also present quenched fractions based on a UVJ (synthetic) color selection.The magnitudes are computed by convolving the star formation and chemical evolution history with photometric tables from Bruzual & Charlot (2003) and accounting for dust attenuation (De Lucia & Blaizot 2007).We tried different cuts commonly adopted in the literature (Williams et al. 2009;Whitaker et al. 2011;Muzzin et al. 2013;Martis et al. 2016) and plot the scatter obtained as shaded regions in Figure 1.For massive and Milky Way-mass galaxies at z > 1, we also used the diagonal selection cut proposed by Belli et al. (2019), which is designed to select recently quenched galaxies.The error boxes show observational measurements of quenched fractions for galaxies from the UltraVISTA Data Release 1 (DR1) and 3D-HST surveys (Martis et al. 2016).In the redshift range 0 < z < 3, there is a good agreement between GAEA2023 and data at all mass scales.The quenched fractions defined using the synthetic UVJ photometry are similar to those defined by sSFR at z < 2 but decrease more rapidly at z > 2. As for GAEA2023, the UVJ color selection underestimates the fraction of quiescent galaxies at z > 2. For TNG, the quenched fractions of galaxies with M å ∼ 10 11 M e are larger than observational measurement.
Moving to higher redshifts, the predicted quenched fractions decrease.In the framework of GAEA2023, the quenched fraction of low-mass galaxies is 0.2% at z ∼ 7.3, i.e., 6 out of 2880 galaxies are quenched.The quenched fraction remains below 1% until z ∼ 3. We traced 42 quenched low-mass galaxies at z ∼ 6.2 forward in time and found that 11 of them moved back to the main sequence within 0.5 Gyr.Most of the high-redshift quenched low-mass galaxies, however, remain quiescent for more than 1 Gyr.
The first quenched galaxies with mass similar to the Milky Way appear at z ∼ 6.2, i.e., 2 out of 116 Milky Way-like galaxies are quenched.One of these two returns to the main sequence after ∼0.5 Gyr.The other is a satellite that remains quenched until it merges with another galaxy.The quenched fraction, for galaxies of this mass, grows quickly to 10% between 4 < z < 5.
The first massive quenched galaxies (with M å ∼ 10 11 M e ) are found at z ∼ 4.5.Two out of 35 galaxies in this mass bin are quenched at this redshift.In their subsequent evolution, these two galaxies continue to have a low star formation rate for most of the time.
Predictions from TNG are quite different from those based on GAEA2023, with systematically lower quenched fractions at high redshift and no quenched galaxy at z > 4.This is in clear tension with the existence of spectroscopically confirmed quenched massive galaxies at z > 4 (e.g., Carnall et al. 2023b).A similar result is reported in Merlin et al. (2019), where a lower sSFR cut 10 −11 yr −1 was used.Though most massive simulated galaxies at z ∼ 4 have a central massive black hole of 10 8 M e , only a minor fraction of them are quenched.In TNG, kinetic feedback from AGN represents the most efficient quenching mechanism for massive galaxies (Kurinchi-Vendhan et al. 2023).However, the accretion rates at high redshift are so large that the assumed mode for AGN feedback in TNG comes in the form of thermal feedback, which is not efficient enough to quench galaxies at z > 3.

Quenching Mechanisms
Broadly speaking, we can consider two quenching scenarios: internal quenching, i.e., AGN feedback and SN feedback, and external quenching, i.e., environmental stripping and galaxygalaxy interactions.In this section, we analyze the relative importance of these quenching mechanisms at different redshifts in the GAEA framework.
First of all, we traced the history of high-redshift quenched galaxies in the three stellar mass ranges considered.Figure 2 shows the evolution histories of three representative galaxies selected from each stellar mass bin.All quenched model galaxies have experienced high-rate black hole accretion and suffered accretion disk feedback right before quenching, suggesting that this mechanism is the main quenching channel at this redshift.This is confirmed by the dotted-dashed lines in Figure 1, showing predictions for Milky Way-mass galaxies from the version of GAEA (Xie et al. 2020) that does not include disk accretion feedback and that significantly underpredicts the fractions of quenched galaxies at high redshift (see also in De Lucia et al. 2024).
Since large accretion rates give rise to luminous quasars, in the top panel of Figure 3 we compare the fraction of quasarhost galaxies in recently quenched galaxies (dashed) and the entire population (solid).The recently quenched galaxies are those that got quenched since the last snapshots, namely quenched in the previous ∼80 Myr at z ∼ 7 and ∼300 Myr at z = 0.The black hole (BH) accretion rate that is used to calculate the bolometric luminosity is also averaged on the same timescale.
About 60% and 30% of massive and Milky Way-mass galaxies host AGN with bolometric luminosity brighter than 10 44 erg s −1 at z > 2. These fractions rise to 100% for recently quenched galaxies at z > 2. Surprisingly, more than 70% of low-mass recently quenched galaxies at z > 4 also host luminous AGN, whereas the fraction is only 10% for all galaxies in this mass range.The elevated fraction of luminous AGN for quenched galaxies confirms that the disk accretion feedback from SMBH is the dominant quenching mechanism for high-redshift galaxies in the framework of GAEA2023.The fraction of luminous AGN decreases at lower redshift, and the differences between all galaxies and quenched samples are also reduced: at low redshift, accretion disk feedback is less important for quenching.
When comparing to observations, the model-predicted AGN fractions are larger than the observational measurements for X-ray-selected AGN at z < 3 (Aird et al. 2022;Ji et al. 2022).However, this is not surprising given that the model-predicted AGN fractions represent upper limits to the actual AGN population, as many of these AGN events may have expired within the timescale we use to estimate the average BH accretion rate.This is especially true for low redshift, where the timescale is much longer and leads to an overestimate of AGN hosts concerning observational measurements.Additionally, the observational measurements should be intended as lower limits and can be significantly affected by selection, obscuration (Hickox & Alexander 2018).We also note that there have been many recent studies reporting discoveries of quenched galaxies hosting AGN (Ito et al. 2022;Carnall et al. 2023a;Bluck et al. 2024), as well as several cases where AGN emission is absent (Jin et al. 2024;Nanayakkara et al. 2024), leaving this issue under debate.
The disk accretion feedback is triggered by both disk instabilities and mergers.While tracing the evolution of individual galaxies, we find the large accretion rates are associated with merger events for Milky Way-mass and massive galaxies.In the middle panel of Figure 3, we show the fraction of galaxies that just experienced mergers between the recently quenched and the entire population.All mergers with a mass ratio larger than 1:100 are considered, motivated by the rapid growth of black holes driven by multiple minor  ), the fraction of galaxies that have experienced recent mergers with a mass ratio larger than 1/100, and the fraction of satellites, respectively.Solid, dotted, and dashed lines correspond to the total, quenched, and newly quenched galaxy populations.Color code is the same as in Figure 1.mergers or even very minor mergers, which we find to be common for high-z quiescent galaxies.Compared to the entire sample of Milky Way-mass and massive galaxies, newly quenched galaxies have much higher merger rates, especially at high redshift.Therefore, mergers represent the main channel for black hole accretion at these mass ranges in the GAEA framework.
For low-mass galaxies, more than half of the recently quenched galaxies have not experienced any mergers around the time of quenching.We thus conclude that their SMBH accretion is not primarily driven by mergers.We find that lowmass galaxies are more likely to have unstable disks at high redshift, and it is this disk instability that triggers efficient black hole accretion.
The connection between the quenching process and black hole accretion becomes weaker at lower redshift, where environmental effects become increasingly important.The bottom panel of Figure 3 shows the fraction of galaxies that are satellites.In GAEA2023, satellite galaxies lose hot gas and cold gas gradually by tidal stripping and ram pressure stripping, whereas central galaxies are unaffected.For low-mass galaxies, a larger fraction of satellite galaxies are quenched as early as z ∼ 6, so the dependence of quenching on the environment starts at very early epochs in the framework of our model.The phenomenon of environmental quenching starts since z ∼ 3 for Milky Way-mass galaxies.The difference in satellite fraction between quenched and all massive galaxies is negligible, which is consistent with previous results that massive galaxies are mainly quenched by AGN feedback (Qin et al. 2017;Xie et al. 2020;Kimmig et al. 2023;Lovell et al. 2023).

Conclusion
In this work we use the semianalytic model GAEA2023 to study the quenched fractions predicted for massive (M å ∼ 10 11 M e ), Milky Way-mass (M å ∼ 10 10.5 M e ), and lowmass galaxies (M å ∼ 10 9.5 M e ) since z = 7. GAEA2023 predictions are in good agreement with the observed quenched fraction measured from UltraVista and 3D-HST in the redshift range 0 < z < 3.
The quenched fractions defined by UVJ color are consistent with those defined by sSFR up to z ∼ 3.At higher redshift, the quenched fractions are underestimated by UVJ color selection.When adopting an sSFR selection, about 5% of massive galaxies are found to be quenched at z ∼ 4.5, about 2% of Milky Way-mass galaxies are first found to be quenched at z ∼ 6.2, and the quenched fraction of low-mass galaxies is 0.2% at z ∼ 7.More than half of galaxies maintain a low star formation rate for over ∼1 Gyr.
All recently quenched Milky Way-mass and massive galaxies at z > 2 and more than 60% of low-mass newly quenched galaxies at z > 4 host luminous quasars (with bolometric luminosity brighter than 10 44 erg s −1 ).This suggests that accretion disk feedback from SMBHs is the main reason for quenching at high redshift, as confirmed by analyzing predictions from an alternative model where this physical process is switched off.We find that disk accretion feedback responsible for quenching is driven by galaxy mergers for massive and Milky Way-mass galaxies and by disk instabilities for lower-mass galaxies.Environmental effects become increasingly important for low-mass galaxies at z < 6 and for Milky Way-mass galaxies at z < 2. Massive galaxies are not quenched by environmental processes in the framework of GAEA (Hirschmann et al. 2016;De Lucia et al. 2019).
The earliest quenched galaxy so far is at z ∼ 7 (Looser et al. 2023), which has a stellar mass of M å ∼ 5 × 10 8 M e .At z ∼ 5, most quenched galaxies discovered are massive galaxies.Based on our model predictions, we expect to find nonnegligible numbers of quenched galaxies with stellar mass ∼10 9.5 M e at z ∼ 7 or even higher redshift.A large fraction of these galaxies are expected to host luminous quasars.

Figure 1 .
Figure 1.Quenched fractions as a function of redshift.Different colors represent galaxies of different stellar masses.Vertical error bars show the standard deviations obtained for 100 randomly selected subsamples.Shaded regions show the uncertainties in quiescent fractions from slightly different cuts in the UVJ diagram (more details in Section 3.1).The dashed-dotted line shows results from a GAEA run where accretion disk feedback is switched off.Dotted and dashed lines show the quenched fractions measured from TNG300 and TNG100, respectively.Error boxes are observed quenched fractions for UltraVISTA DR1 and 3D-HST survey from Martis et al. (2016).

Figure 2 .
Figure2.Evolution of three representative galaxies.The upper and lower panels show the evolution of sSFR and the SMBH accretion rate associated with the disk accretion mode normalized by stellar mass.The sSFR is indicated by solid and dotted lines when a galaxy is classified as central or satellite, respectively.The gray dashed line in the upper panel is the separation between quenched and star-forming galaxies.Large and small triangles in the lower panel indicate merger events with mass ratios above 0.3 and 0.01.Color code is the same as in Figure1.

Figure 3 .
Figure 3.The top, middle, and bottom panels show the fraction of AGN hosts ( [ ] L log erg s 44 bol Red and blue circles show observationally estimated AGN fractions with [ Milky Way-mass galaxies(Ji et al. 2022).Orange open and filled squares show AGN fractions with [ forming and quenched galaxies(Aird et al. 2022).These studies are based on X-ray-selected AGN fractions and should be intended as lower limits for the overall AGN fractions.The model predictions, however, are the upper limits of the actual AGN population.