Tracing the evolutionary pathways of dust and cold gas in high-z quiescent galaxies with SIMBA

G. Lorenzon; D. Donevski; K. Lisiecki; C. Lovell; M. Romano; D. Narayanan; R. Davé; A. Man; K. E. Whitaker; A. Nanni; A. Long; M. M. Lee; Junais Junais; K. Małek; G. Rodighiero; Q. Li

doi:10.1051/0004-6361/202450393

Home

All issues

Volume 693 (January 2025)

A&A, 693 (2025) A118

Full HTML

Open Access

Issue		A&A Volume 693, January 2025


Article Number		A118
Number of page(s)		21
Section		Extragalactic astronomy
DOI		https://doi.org/10.1051/0004-6361/202450393
Published online		09 January 2025

A&A, 693, A118 (2025)

Tracing the evolutionary pathways of dust and cold gas in high-z quiescent galaxies with `SIMBA`

G. Lorenzon¹^⋆, D. Donevski¹^,2, K. Lisiecki¹, C. Lovell³, M. Romano¹^,4, D. Narayanan⁵^,6, R. Davé⁷^,8, A. Man⁹, K. E. Whitaker¹⁰, A. Nanni¹^,11, A. Long¹²^,13, M. M. Lee¹⁴, Junais¹⁷^,18, K. Małek¹, G. Rodighiero¹⁵^,4 and Q. Li¹⁶

¹ National centre for Nuclear Research, Pasteura 7, 02-093 Warsaw, Poland
² SISSA, Via Bonomea 265, 34136 Trieste, Italy
³ Institute of Cosmology and Gravitation, University of Portsmouth, Burnaby Road, Portsmouth PO1 3FX, UK
⁴ INAF, OAPD, Vicolo dell’Osservatorio, 5, 35122 Padova, Italy
⁵ Department of Astronomy, University of Florida, 211 Bryant Space Sciences centre, Gainesville, FL, USA
⁶ Cosmic Dawn centre (DAWN), Copenhagen, Denmark
⁷ Institute for Astronomy, Royal Observatory, University of Edinburgh Edinburgh EH9 3HJ, UK
⁸ Department of Physics & Astronomy, University of the Western Cape, Robert Sobukwe Rd, Bellville 7535, South Africa
⁹ Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada
¹⁰ Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA
¹¹ INAF – Osservatorio astronomico d’Abruzzo, Via Maggini SNC, 64100 Teramo, Italy
¹² Department of Astronomy, University of Washington, Seattle WA 98195-1700, USA
¹³ Department of Astronomy, The University of Texas at Austin, 2515 Speedway Boulevard Stop C1400, Austin TX 78712, USA
¹⁴ DTU-Space, Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark
¹⁵ Dipartimento di Fisica e Astronomia, Università di Padova, Vicolo dell’Osservatorio, 3, 35122 Padova, Italy
¹⁶ Max-Planck-Institute für Astrophysik, Karl-Schwarzschild-Strasse 1, 85740 Garching, Germany
¹⁷ Instituto de Astrofísica de Canarias, Vía Láctea S/N, E-38205 La Laguna, Spain
¹⁸ Departamento de Astrofísica, Universidad de La Laguna, E-38206 La Laguna, Spain

^⋆ Corresponding author; giuliano.lorenzon@ncbj.gov.pl

Received: 15 April 2024
Accepted: 28 May 2024

Abstract

Recent discoveries of copious amounts of dust in quiescent galaxies (QGs) at high redshifts (z ≳ 1 − 2) challenge the conventional view that these objects have a negligible interstellar medium (ISM) in proportion to their stellar mass. We made use of the SIMBA hydrodynamic cosmological simulation to explore how dust and cold gas evolve in QGs and are linked to the quenching processes affecting them. We applied a novel method for tracking the changes in the ISM dust abundance across the evolutionary history of QGs identified at 0 < z ≲ 2 in both cluster and field environments. The QGs transition from a diversity of quenching pathways, both rapidly and slowly, and they exhibit a wide range of times that elapsed between the quenching event and cold gas removal (from ∼650 Myr to ∼8 Gyr). Contrary to some claims, we find that quenching modes attributed to the feedback from active galactic nuclei (AGNs) do not affect dust and cold gas within the same timescales. Remarkably, QGs may replenish their dust content in the quenched phase primarily due to internal processes and marginally by external factors such as minor mergers. Prolonged grain growth on gas-phase metals appears to be the key mechanism for dust re-formation, which is effective within ∼100 Myr after the quenching event and rapidly increases the dust-to-gas mass ratio in QGs above the standard values (δ_DGR ≳ 1/100). Consequently, despite heavily depleted cold gas reservoirs, roughly half of QGs maintain little evolution of their ISM dust with stellar age within the first 2 Gyr following the quenching. Overall, we predict that relatively dusty QGs (M_dust/M_⋆ ≳ 10⁻³ − 10⁻⁴) arise from both fast and slow quenchers, and they are prevalent in quenched systems of intermediate and low stellar masses (9 < log(M_⋆/M_⊙) < 10.5). This strong prediction poses an immediate quest for observational synergy between, for example, the James Webb Space Telescope (JWST) and the Atacama Large Millimetre Array (ALMA).

Key words: evolution / galaxies: evolution / galaxies: high-redshift / galaxies: ISM

© The Authors 2025

Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

The passively evolved, so-called quiescent galaxies (QGs) are often characterised by old stellar populations and a very low or almost absent level of ongoing star formation activity. These features place QGs below the star forming main sequence (MS; e.g. Daddi et al. 2005; Toft et al. 2007; van Dokkum et al. 2010; Rodighiero et al. 2011; Schreiber et al. 2015). There are different paths that star forming galaxies (SFGs) can take to reach the quiescent stage as QGs both from a theoretical (e.g. Feldmann & Mayer 2015; Rodríguez Montero et al. 2019; Zheng et al. 2022) and observational (e.g. Faber et al. 2007; Kriek et al. 2008; Peng et al. 2010; Tacchella et al. 2022; Akhshik et al. 2023) perspective. The conditions required to undergo the different evolutionary processes, which seem to be mainly connected with the timescale of the quenching phase, are still poorly understood (e.g. Park et al. 2023; Corcho-Caballero et al. 2023). The quenching time and the global properties of the QGs are affected by internal mechanisms, such as active galactic nuclei (AGNs; e.g. Crenshaw et al. 2003; Ma et al. 2022; Stacey et al. 2022), star formation activity, and external causes such as ram pressure stripping (e.g. Gunn et al. 1972; Boselli et al. 2020; Junais et al. 2022), starvation (e.g. van de Voort et al. 2017; Maier et al. 2019), and harassment (e.g. Merritt 1983; Moore et al. 1998).

Regardless of the physical mechanism in place, one of the important consequences of quenching the star formation is the evolution of its agents – gas and dust in the cold interstellar medium (ISM). While cold H₂ gas is the fundamental ingredient for star formation, dust is critical in the thermal balance of the gas and in the formation of H₂ molecules (e.g. Cuppen et al. 2017). Metals in the gas phase are thought to support the growth of dust particles (Zhukovska et al. 2016), in turn regulating both the star formation and its quenching. Dust further affects the spectral energy distributions (SEDs) of galaxies, by absorbing the stellar light at shorter wavelengths, and re-emitting it in the far-infrared (FIR). After the star formation ceases, both environmental and internal processes continue to affect ISM evolution, but the quiescent phase of a galaxy has often been associated with a period of simple and predictable evolution. However, recent observations suggest that several processes can still happen within the remaining cold ISM. Namely, various studies unveil anomalously large dust and/or gas content (Gobat et al. 2018; Spilker et al. 2018; Magdis et al. 2021; Belli et al. 2021; Kalita et al. 2021; Morishita et al. 2022; Donevski et al. 2023; Lee et al. 2024; Kakimoto et al. 2024; Setton et al. 2024), which is incompatible with a purely passive evolution. Consequently, quantifying changes in the ISM dust after the cessation of star formation has become one of the key open questions of galaxy formation and evolution.

Despite its importance, very little is known about the evolution of dust and cold gas in QGs, in stark contrast to SFGs, for which a reasonably good census on the dust and cold gas content has been established up to very high redshifts (z ∼ 6; Liu et al. 2019; Donevski et al. 2020; Hodge & da Cunha 2020; Gruppioni et al. 2020; Kokorev et al. 2021; Inami et al. 2022). The main challenge in studying the ISM component of passive galaxies is the faint FIR emission of QGs, and the difficulty in spectroscopically inferring their cold gas masses (M_H₂) and gas phase metallicities (Z_gas), mainly due to the absence of strong emission lines. However, recent observational results have shown great promise, either by stacking the dust continuum of QGs (Man et al. 2016a; Gobat et al. 2020; Millard et al. 2020; Magdis et al. 2021; Blánquez-Sesé et al. 2023), or by large statistical studies of individually detected QGs in deep fields (Donevski et al. 2023). Moreover, a handful of recent works utilised extensive Atacama Large millimetre Array (ALMA) sub-millimetre follow-up of individual QGs in order to attempt to measure M_H₂ of QGs in the distant Universe (z > 1 − 3). This was done either through dust continuum (Whitaker et al. 2021b; Gobat et al. 2022; Blánquez-Sesé et al. 2023; Lee et al. 2024), CO lines (Williams et al. 2021; Belli et al. 2021; Morishita et al. 2022) or atomic carbon emission at 158 μm ([CII], D’Eugenio et al. 2023a).

These studies reported that a non-negligible fraction of massive (10 < log(M_⋆/M_⊙) < 11) QGs at 0.1 < z < 2.5 may be fairly dusty, with average specific dust masses f_dust = M_dust/M_⋆ ≳ 10⁻⁴, systematically lower than in dusty star forming galaxies (DSFGs) at these redshifts. However, these studies also reported few conflicting conclusions regarding dust evolution in QGs. (1) While stacking studies anticipate a strong anti-correlation between dust content and stellar age, statistical studies based on individual detections found a shallower trend accompanied by a large scatter (> 2 orders of magnitude) in f_dust at a given stellar age, M_⋆ and/or Z_gas (Li et al. 2019b; Morishita et al. 2022; Donevski et al. 2023; Lee et al. 2024). The latter finding points towards complex and evolving ISM conditions, and diversity in the quenching mechanisms and/or dust formation pathways. (2) Due to the difficulty in constraining both M_dust and M_H₂, it remains unclear whether QGs obey the canonical dust-to-gas mass ratio (δ_DGR ∼ 1/100) seen in local SFGs, or if they cover a much wider range of values (e.g. Whitaker et al. 2021a; Michałowski et al. 2024). (3) There is no consensus on whether the dust content observed in QGs is of internal origin (owing to past star formation), or if external processes (i.e. mergers, accretion of satellites), have an influential role. (4) Lastly, it is debatable how long the dust and cold H₂ gas will persist after the quenching. Recent studies suggest either a quick depletion of dust and gas (≲100 − 500 Myr; Sargent et al. 2015; Williams et al. 2021; Whitaker et al. 2021b), whereas others report a possibility that timescales for ISM removal may be relatively long (≳1 Gyr; Gobat et al. 2018; Michałowski et al. 2019; Woodrum et al. 2022).

The advent of the James Webb Space Telescope (JWST) has led to the identification of QGs up to very high redshifts (z ∼ 3 − 7, Carnall et al. 2023; Looser et al. 2024), including few candidate dusty QGs (e.g. Rodighiero et al. 2023; Setton et al. 2024), suggesting new and interesting scenarios for the rapid evolution of galaxies in the young Universe. It is then essential to identify and study QGs soon after the quenching episode, as we can then test directly the relative contribution of different quenching processes on the evolution of the cold dust and gas. Because of the many observational challenges outlined above, cosmological simulations of galaxy evolution provide a unique tool for addressing these questions (e.g. Wright & Lagos 2019; Appleby et al. 2020; Donnari et al. 2021). In this work, we made use of the state-of-the-art cosmological simulation SIMBA, which tracks the dust life cycle in a self-consistent way (Davé et al. 2019), meaning that dust grains form, grow and get destroyed accounting for the conditions of the ISM within the simulation itself, evolving together with the galaxy. Our goal is to provide the framework for quantifying changes in the dust and cold gas in QGs during the different stages of their evolution. To this aim, the present study plans to address several key open questions, including (1) understanding what physical mechanisms help replenish the dust reservoirs in QGs after the cessation of star formation; (2) exploring the link between the co-evolution of the dust and cold gas in QGs and the quenching processes affecting them; (3) investigating what is the role of the environment in shaping the fate of ISM in QGs.

We organise this work as follows: in Sect. 2 we present the SIMBA simulation, focusing on the prescriptions for dust evolution and stellar and AGN feedback. Section 3 presents the criteria for the selection of passive galaxies and their environment, defines different representative points of the galaxy evolution and describes the main parameters explored in the analysis. Section 4 describes the predicted global evolution of dust-related properties in QGs as well as their dependence on the environment. In Sect. 5 we provide extensive discussion on the scenarios responsible for the observed diverse evolutionary pathways in QGs. A summary of this work, and possible implications for observational studies, are presented in Sect. 6.

Throughout this work we adopt the same Λ-CDM cosmology of the SIMBA simulation, that is Ω_m = 0.3, Ω_Λ = 0.7, Ω_b = 0.048, H₀ = 68 km s⁻¹ Mpc⁻¹. The same cosmology used in Speagle et al. (2014) is instead used when calculating the star forming main sequence (MS therein), for comparison with our data sample. We assume a Chabrier initial mass function (IMF; Chabrier 2003) consistent with SIMBA.

2. The SIMBA simulation

SIMBA is a state-of-the-art cosmological hydrodynamic simulation that evolves together dark and baryonic matter in a self-consistent way. SIMBA is perfectly suited for our science goal due to its careful treatment of stellar and AGN feedback on sub-galactic scales and active dust modelling. SIMBA allows for a self-consistent evaluation of both dust destruction and dust growth rates at each evolutionary step, together with a rich chemical evolution of the ISM, and a more realistic treatment of AGN feedback with respect to previous iterations, such as MUFASA (Davé et al. 2016).

We made use of the full-physics (m100n1024) simulation that consists of a box of side 100 h⁻¹Mpc containing 2 × 1024³ particles (meaning 1024³ dark matter particles and 1024³ baryon particles). The minimum stellar mass of a resolved galaxy is 5.8 × 10⁸ M_⊙, while the minimum initial mass for a gas particle is 1.82 × 10⁷ M_⊙ and it is 9.6 × 10⁷ M_⊙ for a dark matter particle. The minimum gravitational softening length is 0.5 h⁻¹Mpc. The entire run covers the redshift range 0 < z < 249. Below we describe the major recipes used to model the feedback and dust physics in SIMBA. For a more exhaustive description of SIMBA we refer the reader to Davé et al. (2019).

2.1. Star formation

Star formation in SIMBA is based on the molecular gas fraction following Schmidt (1959) law with a fixed star formation efficiency ϵ_⋆ = 0.02 from Kennicutt (1998). The SFR is proportional to the gas density and molecular gas fraction and inversely scales with the local dynamical time. The molecular gas fraction at a sub-grid level is estimated using the Krumholz et al. (2009a) in order to compensate for the limitations of the simulation’s resolution. The estimation assumes a simple model for the photodissociation of molecular gas in finite clouds (Krumholz et al. 2009b). In the model, the ability of a cloud of molecular gas to self-shield against the interstellar radiation field is strongly dependent on the gas column density and weakly on the metallicity.

2.2. Feedback physics

There are two main categories of feedback in SIMBA: stellar-driven and AGN-driven. For the stellar-driven mechanism, we distinguish between instantaneous feedback from young massive stars and Type II supernovae (SNe) and delayed feedback from evolved stars and Type Ia SNe. Black holes (BHs) are seeded and grown during the simulation using torque-limited (momentum-driven infall in the BH through an accretion disk) and Bondi (isotropic infall) accretion modes. The accretion energy drives feedback that acts to quench galaxies. The AGN-driven mechanism, coming from the black hole accretion, is tuned to mimic the observed dichotomy in black hole growth mode (Heckman & Best 2014) and the observed mode of energy injection into the large-scale gas component. In SIMBA this is obtained by dividing the kinetic feedback into a radiative mode and a jet mode, characterised by different regimes of the Eddington ratio f_Edd = Ṁ_BH/Ṁ_Edd, where Ṁ_BH is the total BH accretion rate (Bondi plus Torque-limited) and Ṁ_Edd is the Eddington accretion rate (accretion rate for which the BH emit at Eddington luminosity). Values of f_Edd > 0.2 correspond to a radiative mode, where mass-loaded winds are ejected radially with velocity of the order of 10³ km s⁻¹. For f_Edd < 0.02 the feedback is a fully jet mode, with hot gas expelled through strictly collimated bipolar jets with zero initial aperture and velocities capped at 7 × 10³ km s⁻¹. The transition regime is characterised by winds of increasing velocity as f_Edd decreases to jet mode. A lower mass limit for the triggering of the jet mode is set to M_BH = 10^7.5 M_⊙ based on considerations from radio observations (Barišić et al. 2017). An additional X-ray feedback heats the gas surrounding the AGN following Choi et al. (2012) when the jet mode is active. The heating is applied directly to the non-ISM gas, deposited as both heat and a radial outward kick for the ISM gas.

2.3. ISM dust in SIMBA

In SIMBA the ISM is treated self-consistently, arising from a direct evolution of the initial conditions which accounts for production, destruction and removal or infall of the different ISM components through feedback mechanisms. The complete description of the dust framework that models the production, growth and destruction of dust grains in SIMBA is introduced in Li et al. (2019a). The baryonic component of a galaxy is initially in the gas phase, and then it is turned into stellar component providing energy injection, represented by a uniform photon field background and a kinetic component that arises from unresolved SNe shocks. The action of both SNe and winds from asymptotic giant branch (AGB) stars are used to update the metallicity value of nearby gas particles. This tracks 9 elements, with good approximation of the whole metal budget: C, N, O, Ne, Mg, Si, S, Ca, Fe. The yields adopted for Type II SNe are from Iwamoto et al. (1999), while for Type Ia SNe and AGB stars, those from Oppenheimer & Davé (2008) are used.

The dust model within SIMBA assumes that dust is produced by condensation of a fraction of metals from Type II SNe and AGB ejecta. The condensation efficiencies for AGB stars and Type II SNe are based on the theoretical models of Ferrarotti & Gail (2006) and Bianchi & Schneider (2007), respectively. The model does not include a contribution from Type Ia SNe, as opposed to some other studies that adopt the same condensation efficiency between Type Ia SNe and Type II SNe (e.g. Popping et al. 2017). Metals can be locked into dust grains on the fly or expelled from the ISM by means of metal-loaded stellar winds. A self-consistent treatment of dust allows grains to grow and be destroyed depending on the energy injection from the unresolved SNe activity and AGN activity. For simplicity, dust grains are assumed to have the same initial size at formation, that is of radius r = 0.1 μm and to follow directly the gas dynamics.

Dust growth is based on a modified version of the Dwek (1998) prescription, deriving the mass of produced dust from the mass of different metals and accounting for the chemical composition and type of the stellar population seeding the ISM. Dust mass can then increase due to accretion in the ISM, based on a two-body collision of grains. This is computed on-the-fly, based on mass, density and temperature of metal-enriched gas hosting the dust (McKinnon et al. 2017).

Destruction of grains can happen by thermal sputtering (Tsai & Mathews 1995), violent action of SNe shocks (Hopkins 2015) or due to astration (consumption by star formation). In the first case, the rate of destruction depends on the temperature, while in the second case, it depends on the mass of the shocked gas and on a parameter ϵ representing the destruction efficiency, since shocks are not resolved in SIMBA. Moreover, dust grains can be destroyed instantaneously if they reside in regions of gas impacted by the X-ray-mode AGN feedback, jets or hot winds, while outflow mechanisms such as the radiative mode AGN feedback and the stellar feedback can heat up and move dust out of the galaxy.

We restricted our results to the redshift range 0 ≲ z ≲ 2 in order to compare it with observationally derived statistics around the cosmic noon. It is worth noting that the SIMBA flagship run results in M_dust functions in a good agreement with observations at z < 2 (Li et al. 2019a), consistently reproduces the observed sub-millimetre number counts up to z ∼ 2 (Donevski et al. 2020; Lovell et al. 2021) and even the observed dust fractions in QGs at z < 1 (Donevski et al. 2023).

3. Methods

3.1. Identification of quenched galaxies

For identifying QGs in the range 0 ≲ z ≲ 2 we adopted a redshift-dependent parametrisation from Pacifici et al. (2016) which defines both a threshold for the star forming phase and a threshold for the quenched phase. Such thresholds are based on the instantaneous specific SFR (sSFR; sSFR ≡ SFR/M_⋆), with SFGs having sSFR > 1/τ(z) and QGs having sSFR < 0.2/τ(z), with τ(z) being the age of the Universe at redshift z. The thresholds have an empirical origin, coming from the observed evolution of the normalisation of the MS (i.e. Whitaker et al. 2012; Fumagalli et al. 2014)¹. We ensure that this method provides a great agreement with a selection based on the distance to the MS (Speagle et al. 2014), as validated in previous SIMBA works on QGs (Appleby et al. 2020; Akins et al. 2022).

In Fig. 1 we exemplify six representative star formation histories (SFH) of SIMBA galaxies that quench at z ≲ 2. The region between the intersection of the SFH with the two thresholds is highlighted in green. Fig. 1 outlines the variety of SFHs for the displayed QGs. It is noticeable that SFHs are not always smooth, and can exhibit large oscillations between subsequent evolutionary points. Due to this irregular behaviour, the sSFR of a galaxy can temporarily fall below the quenching threshold and immediately rise again above it, remaining in a state of low-but-consistent star formation for a long time. Such events must not be accounted for as quenching events since they are ambiguous in nature. We then follow the prescription of Rodríguez Montero et al. (2019) to avoid such cases by interpolating the SFH with a cubic B-Spline (Dierckx 1975). The point at which the SFH drops below the SF and the quenching thresholds is calculated as the intersection of the functions sSFR ≡ 1/τ(z) and sSFR ≡ 0.2/τ(z) respectively with the interpolated SFH. We further request that the sSFR remains below the star forming threshold for a minimum time of 0.2 × τ(z_q), where τ(z_q) is the cosmic time at the redshift z_q of the candidate quenching event.

Fig. 1.

Example sSFR versus look-back time for six representative QGs from this study. The sources’ SFHs are illustrated with blue lines. Panels display QGs along with their corresponding names, quenching times (t_quench), the redshift of the quenching event (z_q), the formation redshift (z_form) and the stellar mass (log(M_⋆/M_⊙). Look-back time is shown on the bottom x-axis in Gyr, while the top x-axis shows the corresponding redshift. The sSFR is illustrated on the y-axis as log-scaled. The SFHs are superimposed with orange and green lines, which mark the thresholds for the start and end of the quenching event, respectively (see the description in Sect. 3.1). The green vertical band delimits the duration of the quenching phase, while red stars are the four critical points we use to determine the sample evolution as described in Sect. 3.2.

A galaxy can also have multiple quenching events if it is rejuvenated after crossing the quenching threshold. Rejuvenation can happen for different reasons, usually involving the formation of a cold gas region through inflows of molecular hydrogen (Young et al. 2014) because of morphological changes, environmental effects or merging events (Rowlands et al. 2018; Woodrum et al. 2022; Chauke et al. 2018; Zhang et al. 2023). It is then important to evaluate these different quenching phases in order to connect them to their respective possible causes. However, as noted in Sect. 5.1 the fraction of rejuvenated galaxies found in our sample appears to be fairly low, so we can neglect the effects of multiple quenching in our statistical treatment. We selected in the specific case of multiple quenching only the first quenching event. The quenching time is calculated as the cosmic time interval that elapsed between the quenching event (at z_q) and the crossing of the star forming phase threshold previous to that (at z_sf):

$\begin{matrix} t_{quench} = τ (z_{q}) - τ (z_{sf}) . \end{matrix}$ $\begin{aligned} t_{\rm quench} = \tau (z_{\rm q}) - \tau (z_{\rm sf}). \end{aligned}$ (1)

As we studied the redshift range 0 < z ≲ 2, we excluded from our analysis all galaxies that do not cross the quenching threshold in this range. In our accompanying work (Lorenzon et al., in prep.) we present the analysis of the effect of formation redshifts to global ISM properties for all SIMBA galaxies.

3.2. Evolutionary phases and timescales

The use of a cosmological simulation allows a simultaneous study from both a statistical and single-case perspective. We linked these two aspects by comparing different galaxies during a few significant phases of their evolution. In this way, it is possible to determine the different SFH that galaxies experience to reach the same stage in their evolution, and to monitor which properties are influencing such paths.

We identified four main evolutionary points that we used in order to monitor the changes in ISM properties of our sample. These points can be understood as critical phases in the galaxy’s history since their formation, and are displayed in Fig. 1. After the formation phase which is the first appearance of the galaxy in the simulation, we identified (Fig. 1, from right to left, respectively) a peak phase (P1) which is the global peak of the galaxy star formation history; quenching phase (P2) which is the start of the quenching process as identified in Sect. 3.1; quenched phase (P3), which marks the end of the quenching phase; gas-removed phase (P4), the point at which the specific gas mass (f_gas = M_gas/M_⋆) drops below the arbitrary threshold of 0.1% that we use to characterise the gas-dry stage. We adopted this limit as such amount of gas is practically negligible and close to the resolution limit of SIMBA for intermediate mass galaxies in our sample. This limit is connected to observations as it is an order of magnitude below the f_gas measured with ALMA for massive (M_⋆ > 10^11.5 M_⊙) QGs at z > 1 (Williams et al. 2021), QGs identified at z ∼ 0.8 in the LEGA-C survey (Spilker et al. 2018), and in local ellipticals (Davis et al. 2016). We are interested in studying how the ISM properties are linked to quenching timescales as well as to stellar properties such as age, stellar mass and distance to the MS. These properties can be measured from observations and can be used to predict the state of the ISM in a given QG. The use of the evolutionary points allows us to compare these global properties consistently for the whole sample.

3.3. Selection of cluster and field galaxies

The analysis is performed using the full-physics SIMBA catalogue of side 100 h⁻¹Mpc, and contains 55 609 galaxies at z = 0. We select clusters using the dark halo mass at z = 0, requiring a virial mass M_halo > 10¹⁴ M_⊙. Clusters are then defined on the basis of their virial radius, which, despite not being representative of the entire structure, selects the region with the largest density contrast with respect to the field. We find a total of 25 clusters for a final sample of 3321 galaxies.

Field galaxies are selected as the SIMBA galaxies at z = 0 that are outside the halos which we identified as clusters. The merger tree of each galaxy from the cluster sample has been reconstructed up to z ∼ 2 using the python library of CAESAR (Thompson et al. 2014), optimised for the handling of SIMBA data. Galaxies in the cluster sample at higher redshift may not be part of a virialised structure (Chiang et al. 2013; Lovell et al. 2018), but they will eventually be part of the dense environment we defined at z = 0. This allows us to ask if the evolution of QGs that end up in massive and dense structures differs from purely field galaxies. The choice of not tracking directly the environment at each epoch avoids the ambiguity in defining high redshift large structures. For the field we used the public information available in the SIMBA catalogues repository.

We selected the final sample of passive galaxies as described in Section 3.1 obtaining a total of 2144 cluster QGs and 7436 field QGs. The galaxies for which it is possible to calculate a quenching time in the redshift range 0 ≲ z ≲ 2 are selected for both cluster and field samples. Due to resolution effects, we limited the stellar mass range of the total SIMBA sample to log(M_⋆/M_⊙) > 9. In this way, we also ensured our sample to be mass-complete.

3.4. Identification of major and minor mergers

We traced major and minor mergers along the galaxy’s history. Mergers are identified using the relative position of galaxies in the SIMBA catalogues. For each galaxy (primary) we explored a volume of ∼0.5 cMpc in radius around it and traced the position of all other galaxies (secondary) in that region at each snapshot. We measured the average maximum proper velocity of galaxies in clusters to be ∼2000 km s⁻¹, and we approximated the typical timescale between snapshots in the redshift range 0 < z < 2 as ∼1 × 10⁸ yr, meaning that we assumed the same galaxy to move of proper motion inside a volume of radius ∼200 ckpc. We traced in this way each galaxy to its next position in the next snapshot. When a secondary galaxy enters a region of radius ∼10 times the radius of the primary galaxy, it is classified as a potential merger. The secondary galaxy can then disappear, reappear or multiple smaller galaxies can temporarily form in that region due to the complex interaction and tidal forces. Each case is considered and evaluated in order to establish the validity of the merger. If the stellar mass ratio of the merging galaxies is R = M_secondary/M_primary > 1:4 then the merger is major merger, while if it is R < 1:4 the merger is minor (Kaviraj 2014; Rodríguez Montero et al. 2019). If multiple minor mergers are happening at the same time, their combined stellar mass ratio is compared to the major merger threshold R > 1:4. Mergers exceeding this threshold are considered a single major merger undergoing intense morphological changes due to tidal forces.

4. The global ISM evolution in cluster and field QGs

To gain insight into the overall evolution of galaxy dust content in our QGs, we used f_dust = M_dust/M_⋆ and δ_DGR = M_dust/M_H₂, which are well-suited proxies to assess the efficiency of the specific dust production and destruction mechanisms (Dunne et al. 2011; Rowlands et al. 2014; Tan et al. 2014; Béthermin et al. 2015; Calura et al. 2017; Michałowski et al. 2019; Burgarella et al. 2020). A fundamental step in the evolution of these ISM properties is the quenching phase of the galaxy. The mechanism responsible for the quenching is considered to act on the molecular gas suppressing the SF. If the dust and molecular gas components are assumed to evolve together with the stellar component (Michałowski et al. 2019), dust is also impacted by the quenching process. If thus the duration of the quenching process is assumed to be connected to its efficiency and, possibly, to entirely different mechanisms, studying t_quench will provide information on both the possible mechanism of quenching and its efficiency in affecting the ISM content.

4.1. Evolution of the quenching mode with redshift

Fig. 2 shows the distribution of the quenching times as measured in Sect. 3.1, normalised by the cosmic time τ(z_q) at the redshift of quenching z_q. Both cluster (red) and field (blue) samples follow a clearly bimodal distribution. We performed a Gaussian fit for each of the two peaks in cluster (field), inferring t_quench, 1/τ(z_q) = − 2.04(−2.16)±0.22(0.22) for the first peak and t_quench, 2/τ(z_q) = − 1.16(−1.07)±0.28(0.28) for the second, respectively. In terms of non-normalised timescales, the peak values correspond to t_quench,1 ∼ 53.7(40.5) Myr and t_quench,2 ∼ 407.2(499.9) Myr at z ≡ 1 and t_quench,1 ∼ 30.1(22.74) Myr and t_quench,2 ∼ 228.4(280.5) Myr at z ≡ 2 for cluster(field). We defined two different regions of the distribution by calculating the intersection point of the Gaussian used to fit each peak. The obtained thresholds are t_{quench, c, lim}/τ(z_q) = − 1.64 and t_{quench, f, lim}/τ(z_q) = − 1.67 for cluster and field galaxies respectively. The bimodal distribution of normalised t_quench has been found in SIMBA for a similar range in stellar masses and redshift (Zheng et al. 2022) and up to z ∼ 4 (Rodríguez Montero et al. 2019) with a similar position of the peaks and a separation threshold (t_{quench, lim}/τ(z_q) = − 1.5), roughly corresponding to t_quench ∼ 175 Myr.

Fig. 2.

PDF distribution of the normalised quenching time for cluster (dark red) and field (blue) galaxies. The quenching time t_quench is divided by the cosmic time τ(z_q) at the redshift z_q at which the galaxy starts quenching. The vertical dotted dark red(blue) line represents the log(t_quench/τ(z_q)) = −1.64(−1.67) separation between fast and slow quenching galaxies in cluster(field).

Galaxies with associated t_quench/τ(z_q) lower than the respective threshold are defined as fast QGs, while larger values identify what we call slow QGs. This definition considers slow QGs as galaxies that quench in a timescale comparable with the cosmic age at the time of quenching. The normalisation by τ(z_q) ensures that this definition remains redshift-independent so that a galaxy that quenches in a short timescale in the local Universe, can be considered slow quenching if its star formation ceases with the same timescale but in a much younger state of the Universe.

A comparison with the literature highlights how SIMBA predicts a peak of much shorter t_quench with respect to other simulations, with a separation threshold which is also generally smaller. Wright & Lagos (2019) finds timescales of the order of some Gyr, with a division around t_quench ∼ 1.5 Gyr from the EAGLE simulations. The characteristic t_quench dividing fast and slow QGs as obtained from the IllustrisTNG simulation is around t_quench ∼ 1 Gyr (Walters et al. 2022). Furthermore, the duration of quenching events we infer in SIMBA qualitatively agrees with the fast and slow quenching timescales reported in observational studies of QGs at intermediate and high-z (Socolovsky et al. 2018; Wu et al. 2018; Belli et al. 2019; Wu et al. 2021; Tacchella et al. 2022; Hamadouche et al. 2023; Park et al. 2023). Observations of massive QGs at z ∼ 0.8 showed a broad range of quenching timescales, indicative of a possible interplay of different quenching mechanisms in shaping the distribution of the sample (Tacchella et al. 2022). A more environment-focused observational study found quenching timescales for clusters at 0.8 < z < 1.35 to be around t_quench = 1.0 ± 0.3 Gyr (Foltz et al. 2018), which is consistent with the lower-end tail (near the separation threshold) of the slow quenching broad peak of the timescale distribution in SIMBA (see also Zheng et al. 2022 Fig. 2). These values, however, do not provide a clear categorisation of the sample in fast or slow quenchers, rather inserting it into an undetermined category in our classification. In Sect. 5, we provide comprehensive analysis on how the quenching routes and their underlying physical mechanisms affect the ISM gas and dust in QGs.

To better have a sense of the evolution of t_quench with the cosmic time, we show in Fig. 3 the fraction of fast and slow QGs at different redshifts with respect to the total number (cluster and field) of QGs in the sample. Normalising by the fixed size of the total sample allows to see how blindingly searching for QGs that quenched in the range 0 ≲ z ≲ 2 would systematically lead to a larger probability of finding a QG in the field. The number of fast and slow quenchers for the cluster sample increases steadily towards z = 0, with slow quenchers dominating at all redshifts. In the field, the situation is analogous for galaxies that quench in the range 0.8 ≲ z ≲ 2.0, but for z ≲ 0.8 the number of fast quenchers rapidly decreases. This scenario is in line with both observations and semi-analytical model predictions of galaxies up to z ∼ 3 (Pandya et al. 2017) showing that at low-z galaxies preferentially evolve through a slow quenching mode, while the fast quenching mode increase in importance at higher redshifts.

Fig. 3.

Fraction of fast (cyan) and slow (red) QGs with redshift for the cluster (continuous lines) and field (dashed lines).

4.2. The δ_DGR-metallicity plane for quiescent galaxies

Gas-phase metallicity is one of the crucial parameters influencing the dust life cycle in galaxies (e.g. Asano et al. 2013; Wiseman et al. 2017; Popping & Péroux 2022). The unique dust model implemented in SIMBA enables us to investigate how the different markers of galaxies’ dust-life cycles (i.e. δ_DGR and f_dust) evolve as a function of Z_gas (expressed as 12 + log(O/H)).

In Fig. 4 we show the distribution on the dust-metallicity plane for cluster(field) QGs. The panels represent from top to bottom the samples observed at the start of their quenching phase (P2), at the end of the quenching phase (P3), and at the last evolutionary point (gas-removed phase, P4), respectively. Background contours show the 50th and 90th percentile level of the distribution at the peak of the SFH (P1). As described in Sect. 3.2 we defined the evolutionary point P4 using an arbitrary threshold for the cold gas fraction, namely f_H₂ = M_H₂/M_⋆ = 10⁻⁴. If a galaxy does not cross this threshold then it is not counted in the plane. The dashed vertical line marks the solar oxygen abundance 12 + log(O/H) = 8.69 (Asplund et al. 2009), while the horizontal dot-dashed line is the literature reference value (δ_DGR ∼ 1/100, i.e. Lisenfeld et al. 2000; Leroy et al. 2011; Magdis et al. 2011; Sandstrom et al. 2013; Rémy-Ruyer et al. 2014; Whitaker et al. 2021a). Arrows show the position of the galaxies presented in Fig 1. We show the empirical fit δ_DGR = 10.54 − 0.99 × (12 + log(O/H) from Magdis et al. (2012) as a grey line, and we shade in grey the 0.15 dex scatter of the relation. The fit was done on a sample of star forming galaxies.

Fig. 4.

δ_DGR-metallicity plane for cluster QGs (upper panel) and field QGs (lower panel). Positions of QGs within the plane are shown for different evolutionary phases, as indicated in the legend. Contours define positions typical for the peak star formation phase. The intensity of colours refers to the number of sources. The vertical dashed line marks the solar value for the oxygen abundance (Z_⊙), while the horizontal dashed-dotted line shows the reference value for the dust-to-gas ratio in SFGs from the local Universe δ_DGR = 1/100. The positions of galaxies introduced in Fig. 1 are indicated with arrows. The scaling relation between Zgas and δ_DGR for local SFGs (Magdis et al. 2012) is shown as a solid line with equation indicated in the legend. Histograms on the y-axis show the PDF distribution of log(δ_DGR) for both samples, with small ticks marking the 16th and 84th percentile.

Our sample recovers the same global trend observed for the entire population of galaxies at z = 0 in Davé et al. (2019). The distribution has a characteristic, environment-independent ‘boomerang shape’. It reconciles a linear increase of the δ_DGR with the oxygen abundance, combined with a large scatter towards super-solar Z_gas (10⁻¹ ≤ δ_DGR ≤ 10⁻⁵).

Surprisingly, soon after quenching, ∼50% (∼40%) of cluster (field) QGs attain f_dust ≳ 10⁻³ and thus roughly follow the trend expected for SF galaxies (i.e. best-fit relation from Magdis et al. 2012). The scattered region which deviates from the fit is populated by ∼30% (50%) of cluster(field) galaxies, respectively. As the galaxies evolve, the ISM reservoir sustains and has measurable Z_gas, M_dust and M_H₂, up to the quenched phase (P3) in ∼95% of QGs in the field sample. For cluster QGs, the fraction falls to ∼80%. In the remaining galaxies, the cold gas content is below the resolution of SIMBA, so that the properties of dust, metallicity and SFR are not measurable. The linear trend of the distribution is significantly populated even in the dry phase (P4; ∼30% for clusters and 20% for field QGs), showing no indication of evolution with cosmic time. This behaviour resembles the one observed in DSFGs, for which a redshift independent relation of δ_DGR with Z_gas has been reported both at low redshifts (z ∼ 0, e.g. De Vis et al. 2019) and at high redshifts (1.5 < z < 2.5; e.g. Shapley et al. 2020; Popping et al. 2023).

The ’boomerang-shaped’ dust-metallicity plane from Fig. 4 reveals several interesting features: (1) The bulk of our QGs populates the same part of the plane as SIMBA SFGs (see Fig. 3 in Li et al. 2019b). Consequently, these QGs also follow the trend expected from the empirical relation between the δ_DGR and Z_gas obtained for MS DSFGs (Magdis et al. 2012). Such a well-behaved trend is mostly pronounced for galaxies of intermediate masses (see Appendix A). Intriguingly, a non-negligible number of QGs has δ_DGR that is compatible or even higher, than observed for MS DSFGs. As an example, in Fig. 4 we also highlight that four out of six QGs from Fig. 1 maintain (or even increase), their δ_DGR during post-quenching evolution. This is clearly visible in Appendix B, showing that for the same 4 galaxies (G1, G3, G4, G5) f_dust remains relatively constant, with changes of less than an order of magnitude during their evolution up to P4. (2) By confronting the upper and lower panels of Fig. 4 it is evident that the environment does not introduce the change in the shape of δ_DGR and Z_gas plane, but reduces the dust abundance more efficiently in the cluster than in the field. Namely, both samples have a strong peak around the typical literature value of δ_DGR ∼ 1/100, but field QGs cover the plane more uniformly along the region of large scatter. The difference between the distributions is strongest towards the last evolutionary point (P4) that describes the gas-dry phase. The fraction of QGs in this phase is rather similar, with only 45%(47%) of galaxies retaining an amount of ISM above the resolution limit of SIMBA in the field and cluster environment, respectively.

For ∼10% of the field QGs, the gas fraction remains above our arbitrary threshold (f_gas > 10⁻⁴) over the entire evolution following the quenching event. Such QGs do not satisfy the ‘dry-phase’ criteria, therefore their P4 point is not calculated. The actual fraction of field QGs with a negligible trace of ISM is thus around 30–40%. This number is only ∼4% for cluster QGs. For both cluster and field QGs the bulk of the distribution remains centred around δ_DGR ≳ 1/100. These QGs tend to maintain a large fraction of dust during their last phases of passive evolution where the cold gas content drops significantly; (3) At 12 + log(O/H) ≳ 8.8 there is a range of δ_DGR that extends over ∼4 orders of magnitude independently of evolutionary phase or environment. The absence of a clear trend in Fig. 4 suggests that the balance between the formation and destruction of dust grains is strongly altered at super-solar Z_gas. The large scatter appearing for metal-rich QGs has been observed locally in Casasola et al. (2020) and subsequently at intermediate redshifts in Donevski et al. (2023). A number of processes can contribute to this scatter, that is, ISM-rich mergers, an increase in f_gas by environment-dependent inflows (Aoyama et al. 2022), and efficient re-growth of dust grains that can compensate the destruction due to thermal sputtering (Zhukovska et al. 2016; Pantoni et al. 2019; Gillman et al. 2021). We explore possible scenarios in Sect. 5. We note that throughout the rest of this work, we show the analysis for the sample of field QGs since we clarified in Fig. 4 that the field and cluster samples have the same ISM properties even if differently distributed.

4.3. ISM dust-gas abundance within the MS

To better understand the ISM evolution of our quenched sample, we trace the change in f_dust with respect to the offset from the MS (Δ_MS), as defined by Speagle et al. (2014)². In Fig. 5 we illustrate this relation along the four key evolutionary stages introduced in Sect. 3.2.

Fig. 5.

Evolution of f_dust with Δ_MS for the field QGs. The Δ_MS is derived following the best-fit relation from Speagle et al. (2014). Each column represents a different phase of the galaxy evolution as indicated in the legend. These are, from left to right: sSFR-peak phase (P1), quenching phase (P2), quenched phase (P3), and gas-removed phase (P4). Each row corresponds to different stellar mass bins calculated at the quenching time. Stellar masses increase from top to bottom, while the sample evolves from left to right. The grey shaded region defines the 0.3 dex scatter of the MS. Points are colour-coded for their δ_DGR calculated at the quenched phase (P3). Points with the same colours in different columns represent the same galaxy across different evolutionary stages. The horizontal arrow marks the direction of evolution with respect to the MS, while the vertical arrow shows the direction expected for increasing dust growth. For each panel, we show the fraction of fast(slow) QGs in the column corresponding to P3. To guide the eye, we show a dashed horizontal line (f_dust = 10⁻⁴), which roughly corresponds to a cold gas fraction of 1% and marks an approximate observational limit of known dusty QGs detected with ALMA (see Section 5.2.4).

As QGs evolve, they get further away from the MS, but we observe little to no changes in f_dust during the first two phases of their evolution. This behaviour is noticeable for all mass bins. Nevertheless, the later evolution of ISM dust and gas is different for QGs of different M_⋆. Namely, for QGs of lower masses (log(M_⋆/M_⊙) < 9.75) the dust content remains remarkably stable across the entire evolutionary history. More complex behaviour emerges for more massive QGs (log(M_⋆/M_⊙)≳9.75), for which we observe a larger scatter in f_dust after the quenched phase. The QGs of these stellar masses show no clear evolution with Δ_MS. Our data indicate a general trend going from the top to the bottom row, with higher masses having on average lower f_dust at the gas-removed phase. This implies that M_⋆ has an important impact on the late evolution of ISM dust and gas in QGs, such that less massive QGs are characterised by a lower dust destruction rate and/or by a larger growth rate than in more massive galaxies. The significant spread in f_dust at the given MS-offset is interesting to explore, as QGs are generally expected to share similar, redshift independent SEDs, a consequence of a less complex and less turbulent ISM as compared to DSFGs (e.g. Magdis et al. 2021). Therefore, it has been suggested that the relation between the f_dust and Δ_MS can be interpreted as an ‘age-evolutionary sequence’, such that QGs from the upper-right side of an f_dust − Δ_MS diagram are dominated by objects with younger stellar age and higher gas fraction. The later decrease in star formation rate (SFR) is likely due to the exhaustion of their gas reservoirs, reflecting efficient gas removal (Gobat et al. 2020). However, the large scatter and absence of a clear trend during the later phases of evolution, suggests a complexity of processes affecting the ISM dust abundance. Whitaker et al. (2021b) suggests that a large range in f_dust hints at diverse evolutionary routes to quiescence, a question we thoroughly explore in the subsequent sections.

Furthermore, Fig. 5 illustrates that δ_DGR and f_dust follow each other until the start of the quenched phase. Hence we pay particular attention to the last column of the plot in Fig. 5. This evolutionary point is set around a specific value of f_H₂, that is f_H₂ ∼ 10⁻⁴ (see Sect. 3.2), meaning that the main ISM property changing between the different galaxies inside a mass bin is f_dust. At the start of the quenched phase, the scatter in f_dust presents a smooth gradient in δ_DGR, with low δ_DGR being associated to low f_dust. This trend is not maintained up to the final, gas-removed phase. The last column of Fig. 5 showcases QGs with low f_dust and δ_DGR in the quenched phase, attaining higher values of f_dust even with f_H₂ ∼ 10⁻⁴. The vertical arrow indicates the direction along which the dust growth could have an effect. If there is a mechanism preventing dust grains to be destroyed efficiently, the major effect would be to increase both f_dust and δ_DGR. From the column corresponding to the P4 point in Fig. 5 we see that the bulk of QGs with f_dust > 10⁻⁴ at the gas-removed phase is for low mass galaxies (9 < log(M_⋆/M_⊙) < 9.75), but a non-negligible amount is found for galaxies in the range 9.75 < log(M_⋆/M_⊙) < 10.5 and few galaxies are found for even larger masses. This suggests that even for QGs with virtually no cold ISM gas (that is the case in P4 by construction), copious amounts of dust can be found. The number of galaxies with large f_dust for the intermediate and large mass ranges is likely a consequence of the underlying mass distribution as shown in the mass division of the δ_DGR-metallicity plane (Appendix A).

5. Exploring the mechanisms influencing the dust-gas evolution in high-z QGs

To reveal the key sources influencing a complex dust-gas interplay presented above, in this Section we provide an in-depth analysis of both external and internal mechanisms that may be responsible for the overall evolution of the ISM in QGs. The focus will be set in particular on exploring the origin of the large scatter in f_dust. This allows us to characterise the properties of QGs with copious amounts of dust in the ISM and to provide a qualitative evolution of their ISM content.

5.1. External channels affecting ISM dust evolution

5.1.1. Cluster environment

As seen in Fig. 4, the fraction of QGs with measurable ISM dust in the super-solar Z_gas regime is reduced more significantly for the cluster QGs than for the field ones. As a consequence, a large number of QGs with super-solar Z_gas has δ_DGR that extends over > 4 orders of magnitude, including exotically low values (δ_DGR ≲ 10⁻⁴) that are first reported in SIMBA work by Whitaker et al. (2021a). Such low values are similar to the low upper limit (δ_DGR ∼ 10⁻⁴) obtained from Planck observations of galaxy clusters in the local Universe (Planck Collaboration XLIII 2016). It is generally considered that thermal sputtering on dust grains embedded in a hot medium must be relatively rapid in massive halos typical for galaxy clusters (e.g. Galliano et al. 2021). The timescale for dust destruction by hot gas sputtering depends on the temperature T_gas and density n_e of the gas as $t_{sp} = 10^{5} (1 + {(10^{6} K / T_{gas})}^{3}) n_{e}^{- 1} yr$ $t_{\mathrm{sp}} = 10^{5}(1+(10^{6} K/T_{\mathrm{gas}})^{3})n_{\mathrm{e}}^{-1}\, \rm yr$ (Hirashita et al. 2015). This means that at the typical conditions of the hot halos of massive QGs (T_gas ∼ 1 × 10⁷ K, n_e ∼ 1 × 10⁻³ cm³), the hot gas is expected to support rapid destruction of dust, with a timescale of the order of t_sp ∼ 100 Myr. Since the sputtering is efficient for T_gas > 1 × 10⁶ K and decreases rapidly with t_gas < 10⁶ K (Nozawa et al. 2006), the effect is expected to be more significant in the cluster with respect to the field galaxies. The fraction of QGs in SIMBA during the dry, gas-removed phase is different for the two environments. We illustrate this in Appendix A, showing how the bulk of cluster QGs have almost entirely removed their ISM dust during the final evolutionary stage, in contrast to field QGs that keep the low, but measurable level of f_dust.

Vogelsberger et al. (2019) uses the AREPO cosmological simulation (Springel 2010) and finds a typical timescale of 10 Myr to 1 Gyr for the dust destruction within the cluster core (around 100 kpc to 1 Mpc from the centre). We calculated the timescale for SIMBA using the temperature of the circumgalactic medium (CGM), which in the simulation catalogues tops at around 0.5 − 1 × 10⁵ K. For such temperatures the timescale is expected to be larger than for the typical scales of T ∼ 10⁶ − 10⁷ K, namely t_sp ∼ 100 Myr. Despite this, a large population (∼30%) of galaxies in SIMBA clusters show dust fractions values well above the measurements of Planck Collaboration XLIII (2016), and consistent with MS galaxies (Magdis et al. 2012). This suggests a physical process that is able to alleviate the grain destruction due to thermal sputtering, most likely by prolonging the timescales by one order of magnitude at least. This is in line with recent works (Priestley et al. 2022; Relaño et al. 2022) showing that longer sputtering timescales of t_sp > 0.5 − 1 Gyr support the ongoing production of larger dust grains, which are expected to be reduced relatively slowly in hot halos. The possible dust re-growth scenario has also been proposed by Hirashita et al. (2015) and Hirashita & Nozawa (2017), as an explanation for the observed high f_dust in local elliptical QGs. These studies concludes that the variation of the dust growth efficiency is caused by large variations observed in δ_DGR, similar to those seen in our QGs. Whitaker et al. (2021a) conduces several numerical experiments with higher-resolution SIMBA runs, finding that both thermal sputtering and dust growth are needed to explain the observed δ_DGR in galaxies with low sSFR. We return to this point in Sect. 5.3.

5.1.2. Mergers

Merging events (in particular minor mergers) are important contributors to the size growth for massive QGs during their passive evolution (Bundy et al. 2009; Kaviraj 2014; Matharu et al. 2019). Generally, mergers can form galaxies with different amount of dust depending on the dust content of the progenitors. Several studies propose that mergers with dust-rich galaxies can significantly enhance ISM dust abundance in QGs, potentially acting as the main channel for dust enrichment (Rowlands et al. 2012; Martini et al. 2013; Lianou et al. 2016; Dariush et al. 2016; Kokusho et al. 2019).

We explored the occurrence of merging events among our QGs. From a complete track of the merger tree of QGs up to z ∼ 2 we inferred that ∼20% of QGs experience at least one major merger in this redshift range and ∼15% of QGs experiencing multiple minor mergers. The major merger fraction is ∼10% during the gas-removal phase (from P2 to P4), which can last up to several Gyr and occurs mostly for z < 1. This all confirms a non-negligible fraction of merging events in SIMBA galaxies as previously reported in Rodríguez Montero et al. (2019). Our merger fractions are consistent with those reported in the observational study of Man et al. (2016b). They use a large catalogue of mergers observed up to z ∼ 3, and for major and minor mergers within the same redshift range as in this study, they infer ∼20% and ∼15%, respectively. We further explored whether mergers are the main drivers of the scatter in dust-related quantities as observed for the quenched phase in Fig. 4 and Fig. 5. We inferred a relatively low fraction of mergers with respect to the fraction of QGs populating the δ_DGR scatter, which we measured to be ∼45%(∼20%) for QGs in the upper region of the scatter (log(δ_DGR) > − 3 and 12 + log(O/H) > 8.69), and ∼55%(∼36%) considering the entire scattered region (no constrain on δ_DGR) for P3(P4), respectively. We also checked and confirmed that the highest δ_DGR and f_dust are not prevalent in mergers. Consequently, we conclude that the observed amount of dust in QGs cannot be explained by considering the stochastic action of mergers only. The conclusion applies to both cluster and field QGs.

Our conclusion is in line with recent observational results reporting that dust in QGs is not preferentially of external origin (e.g. Michałowski et al. 2019; Blánquez-Sesé et al. 2023). It also qualitatively agrees with Rodríguez Montero et al. (2019). They trace the evolution of the merging and quenching rates in SIMBA, and report insufficient merger rate to explain the observed quenching rate at z ≤ 1.5. The same study reports that fast quenchers are dominating the mass range around 1–3 × 10¹⁰ M_⊙, which corresponds to the second row of Fig. 5, where the scatter appears more prominent in both P3 and P4. Moreover, they do not find a significant physical connection between mergers and fast or slow quenchers. Since the scatter appears during the quenching phase, the conclusions of Rodríguez Montero et al. (2019) suggest that a physical connection between the scatter and the merger process must be weak.

5.1.3. Rejuvenations

As explored in Rodríguez Montero et al. (2019), rejuvenation events can occur in SIMBA mostly in relation to minor mergers with gas-rich satellites. Bursts in the sSFR can be triggered by mergers and galaxy fly-bys and can thus be associated with a redistribution of the ISM in the galaxy and a burst in dust formation due to stellar activity. We make use of the empirical thresholds for star formation in order to identify rejuvenation events. We required a galaxy to be quenched for at least 0.2 × τ(z_q), with τ(z_q) being the cosmic time at quenching. A rejuvenation happens if the sSFR of a quenched galaxy is sSFR(z)≥1/τ(z), that is if the galaxy returns to the MS. The resulting fraction of rejuvenated galaxies is ∼3% for the clusters and ∼2% in field galaxies. Some examples of these objects are clearly seen in Fig. 5, as they re-enter the MS region during the late stages of evolution. We further checked and found that our rejuvenated galaxies are not limited to a particular section of δ_DGR-metallicity plane, or f_dust-MS plane, but populate its entirety. This all suggests that rejuvenations are not prevalent for producing a large scatter in f_dust and δ_DGR in our QGs.

Seemingly low rejuvenation fractions have recently been obtained from an independent work of Szpila et al. (2024). They use Simba-C (Hough et al. 2023) QGs, and obtain a rapid downfall of rejuvenating fraction at z < 2 to the point that < 5% of galaxies at z = 0.5 experience rejuvenation before z = 0. The trend reverses at z > 2, where Simba-C predicts a fraction > 30%, along with a population of high redshift QGs that reaches z ∼ 5. The low fraction of rejuvenations for low-z galaxies can be misleadingly attributed to a lower amount of cosmic time left until z = 0 for the rejuvenation to have a large chance of occurrence. However, observations of the contribution of O-type stars show that rejuvenation generally occurs during the first < 1 Gyr after the quenching event (Zhang et al. 2023).

The lack of low-z rejuvenated galaxies seems to be physically supported rather than an effect of redshift selection. Interestingly, some observations seem to point to an opposite trend, with the local Universe showing that 10–30% of z < 0.1 QGs have signs of ongoing rejuvenation (Treu et al. 2005; Schawinski et al. 2007; Donas et al. 2007; Thomas et al. 2010), while this number decreases to an uncertain 1–16% for z > 1 (Belli et al. 2017; Chauke et al. 2019; Tacchella et al. 2022; Woodrum et al. 2022; Akhshik et al. 2021). Such discrepancies are due to the lack of consensus on the definition of rejuvenated galaxy, of which an accurate comparison between simulated and observational data must carefully take into account, as highlighted by Alarcon et al. (2023). Here we explored a simple, but constraining definition for SIMBA rejuvenated galaxies, requesting the overall SF of a galaxy to be comparable with the low end of the MS, which can result in a selection of only the strongest and most extended onsets of new star formation. We also verified that lowering the SF threshold for the selection of rejuvenation to sSFR = 0.5/τ(z) can increase the fraction to the observed 10-20% at z < 2. However, this choice is not physically motivated. As described in Szpila et al. (2024), further constraints on the duration of the rejuvenation event would result in an even lower fraction.

5.2. Internal mechanisms

Analysis from Sect. 5.1 reveals that external mechanisms are not capable of producing the major effect on the observed dust-gas evolution in QGs. We therefore turn our attention to investigate in detail the internal processes affecting them. For simplicity, throughout the following sections, we present the trends for the field QGs only but ensure that these do not change if we consider objects that reside in clusters.

5.2.1. Impact of AGN feedback

The presence of an AGN can have a significant impact on the evolution of the ISM. The energy from the AGN is injected into the ISM through the ejection of winds and high-energy radiation. Winds can sweep ISM material from the nuclear region of the galaxy into the CGM. If the accretion of matter into the central black hole is rapid enough, bipolar highly collimated jets may form and effectively remove the ISM material, increasing its temperature. Comparing the Illustris, Illustris TNG and SIMBA simulations, Ma et al. (2022) concludes that the main effect of internal quenching has to be found in the cumulative energy released during the central black hole growth. This energy affects the formation of cold gas, leading to the cessation of star formation. Using an analogous diagnostic, we explore the effect of the AGN energy injection on the ISM.

The jet mode of energy injection in the ISM from the AGN is triggered when f_Edd < 0.02 and M_BH > 10^7.5 M_⊙Davé et al. (2019). Figure 6 illustrates the effect of feedback variations (quantified with the Eddington ratio and M_BH) on the ISM dust abundance in QGs. We present the results for field QGs as a reference for comparison. All properties of the sample are measured at the quenched phase of the galaxies (P3), where the effect of star formation is minimal. The majority (∼90%) of fast-quenchers exclusively reside in the jet-dominated region of the plot. Slow quenchers are more evenly divided, with ∼26% of QGs being influenced by the radiative mode feedback, and jet-dominated AGNs influencing the rest. We also find that the AGN jet feedback is most commonly present in the stellar mass range of 1 − 5 × 10¹⁰ M_⊙, confirming the earlier results of Rodríguez Montero et al. (2019), Appleby et al. (2020) and Yang et al. (2024).

Fig. 6.

Distribution of the f_Edd over the central black hole mass M_BH of the sample of fast (left) and slow (right) field QGs at the quenched stage. The colour coding shows the median dust fraction f_dust computed for the quenched stage (that is, at the P3 evolutionary point). The vertical dashed line represents the lower mass limit for the activation of jet mode in SIMBA (M_BH = 10^7.5). The dashed horizontal line delimits the starting of the transition from a radiative feedback mode to a jet mode (f_Edd = 0.2), while the dashed-dotted line delimits the region with fully set jet mode (f_Edd < 0.02). The red and blue contours represent the 50 and 90 percentile levels for the slow and fast quenching sub-samples respectively.

A net drop of ∼2 orders of magnitude in f_dust is visible going from the radiative mode to the jet mode dominated region in the sample of slow QGs. Moreover, QGs in the jet-dominated region of slow QGs have a larger scatter in their f_dust with respect to the QGs experiencing the radiative mode feedback, as visible in Fig. 6. However, by looking at the impact of quenching routes on f_dust no explicit division can be found. Namely, the majority of both fast and slow quenchers appears dust-rich at the quenched phase, with median log(f_dust)≳ − 4. Interestingly, fast quenchers have predominantly high f_dust. Their median log(f_dust) is high in both feedback variants ( $- 3 . 09_{0.57}^{0.27}$ $-3.09^{0.27}_{0.57}$ in the jet mode, and $- 3 . 06_{0.39}^{0.25}$ $-3.06^{0.25}_{0.39}$ in the radiative mode, respectively). These values are 2 − 5× higher than those of slow quenchers ( $- 4 . 81_{0.83}^{1.18}$ $-4.81^{1.18}_{0.83}$ and $- 3 . 25_{0.45}^{0.26}$ $-3.25^{0.26}_{0.45}$ , respectively). In other words, while AGN jet feedback is primarily responsible for generating a population of rapidly quenched QGs, it does not correspond to an immediate removal of dust during the quenched phase.

This all points towards an important complexity between the timescales for the cessation of star formation, gas removal and dust evolution in the late phases of QGs. The conclusion is qualitatively consistent with recent observational studies (e.g. Greene et al. 2020) suggesting that the peak of AGN activity may occur during an earlier phase of the quenching process but it is possible that the AGN jet mode feedback has not fully cleared or destroyed the cold gas reservoir. Such a possibility is further sustained by Herschel observations of large amount of dust in outflows from NGC 3079 in Veilleux et al. (2021). They found that the sputtering timescale in the outflowing filament is consistent with the timescale of AGN activity required to lift the filament from the ISM. From observation of star formation-driven outflows in dwarf galaxies McCormick et al. (2018) found dust entrained by galactic-scale winds, suggesting the possibility of shielded regions in which dust can survive thermal sputtering and collisions. Dust surviving in the CGM is also observed by Romano et al. (2024), likely carried out of the ISM by outflows. This provides a scenario in which the dust is expelled from the ISM, but not efficiently destroyed, and can then rain back down in the galaxy, as observed with the filament of NGC 3079 reconnecting to the galactic plane. The H₂ gas, on the other hand, is observed to be deposited in the CGM, where its destiny is uncertain and will depend on the AGN activity. This underlines possible different evolution pathways for dust and gas in the ISM. Further evidence comes from recent observations of dust-obscured quasars at z ∼ 2.5 showing a correlation between the reddening of the source and the presence of radiation-driven molecular outflows (Stacey et al. 2022). The study finds that dust and molecular gas can survive after being removed from the galaxy even in the most luminous quasar in a process that effectively and rapidly quench the galaxy nuclear region. This scenario has been supported by high resolution simulations of galactic winds as AREPO-RT (Kannan et al. 2021) and Cholla (Schneider & Robertson 2015) for different grain sizes and distances from the galactic plane (Richie et al. 2024). Dust survival depends on the conditions of the ISM, with SNe-heated gas being more efficient in destroying dust than gas affected by the AGN. After being deposited on the CGM, the surviving material can eventually rain down on the galaxy bringing new gas and dust to the ISM. The amount of ISM replenished in this way would depend on the potential well of the galaxy, the energy injected on the outflowing material, and the conditions of the CGM. This delicate balance is not yet fully understood. In Sect. 5.2.3, we further investigate this possibility by analysing the properties of fast and slow QGs and their influence on the f_dust evolution.

5.2.2. X-ray feedback

AGNs in SIMBA contribute to the energy injection into the ISM and CGM also through X-ray radiation. This effect is triggered during the full-velocity jet-phase of the AGN and it depends on f_gas since the gas surrounding the super-massive black hole (SMBH) can absorb and radiate away the X-ray energy. The main effect of the X-ray contribution is on the heating of the gas component, on the non-ISM gas and a combination of heating and momentum transfer for the ISM gas around the SMBH. This leads to two main consequences: on the one hand, massive halos have larger T_gas, especially in the central region, as observed in X-ray profiles of SIMBA by Robson & Davé (2020); on the other hand, gas is efficiently removed from low-mass halos and is pushed outwards from the central region of galaxies, preventing formation of new stars, as found in Appleby et al. (2020).

At the start of the quenching phase, ∼30% of both field and cluster QGs from our sample satisfy the condition for X-ray heating activation. This fraction, however, increases differently towards the end of a quenching phase, reaching ∼60% for the field QGs, and only ∼40% for the cluster QGs. This is likely connected to a slightly larger fraction of radiative mode dominated QGs in clusters, since galaxies that do not experience jet mode do not trigger X-ray heating in SIMBA. The lack of strong AGN feedback manifesting through both collimated jets and powerful X-ray radiation can help explain the extended tail of slowly-quenched QGs containing substantial ISM dust abundance. Only a handfull of such QGs is observed amongst the fast quenchers.

5.2.3. Impact of quenching routes and timescales of H₂ gas removal on the dust content of QGs

One of the main questions arising from the presented results is how the dust in QGs is impacted by the different quenching routes and the consequent variety of timescales for the cold gas removal that follows the quenching process. The quenching timescale (as introduced in Sect. 3.2) reflects the transition period leading a galaxy to quiescence, and it is linked to the quenching mechanism itself. If, after the cessation of the SF, the changes of cold gas and dust tightly follow each other, one can conclude that the evolution in f_dust is solely dependent on the quenching route (see e.g. Gobat et al. 2020). This assumption is at the basis of the age-evolutionary framework, which proposes that once quenched, QGs undergo consumption (or destruction) of the dust and cold gas without further replenishment. Such evolutionary picture requires a strong anticorrelation of f_dust with stellar age (Michałowski et al. 2019; Caliendo et al. 2021; Whitaker et al. 2021a; Blánquez-Sesé et al. 2023).

To probe this scenario, we explored the overall evolution of f_dust with the mass-weighted stellar age. Figure 7 showcases this relation for the same four critical points of SFHs as presented in Fig. 5. We colour-code with the timescales for cold gas removal, defined as time elapsing between the beginning of a quenching phase (P2) and the beginning of the subsequent gas-removed phase (P4) as:

$\begin{matrix} t_{removal}^{H_{2}} = τ (z_{P 4}) - τ (z_{P 2}), \end{matrix}$ $\begin{aligned} t^{\mathrm{H}_2}_{\rm removal} = \tau (z_{\rm P4}) - \tau (z_{\rm P2}), \end{aligned}$ (2)

Fig. 7.

Evolution of f_dust with stellar age. We track the same evolutionary stages as introduced in Sect. 3.2. Points are colour-coded with the mean quenched timescale, that is the time interval between the P2 and P4. grey dashed and dashed-dotted lines show the exponential decrease of dust abundance assuming a timescale of τ = 2.5 Gyr and τ = 0.7 Gyr respectively and are motivated by Michałowski et al. (2019) and Lee et al. (2024). White squares denote the binned medians from the sample of dusty QGs observed in the COSMOS field (Donevski et al. 2023). Other markers show various QGs followed-up with ALMA: empty star (z = 2, Lee et al. 2024), empty circles (z ∼ 1.5 − 2, Gobat et al. 2022), and empty diamond (z = 2.149, Morishita et al. 2022). Cyan(red) contours define the 50th and 90th percentile of the distribution inferred for the fast(slow) QGs from this study, respectively. In the bottom left panel, the three main peaks of the distribution are labelled A, B, and C.

where τ(z) is the cosmic time at z. Cyan(red) contours represents the 50th and 90th percentile of the distribution of the fast(slow) QGs. Fig. 7 reveals that QGs in SIMBA start their quenching phase with a wide distribution in stellar ages and with a relatively narrow range of f_dust. While moving from the peak of sSFR to the start of the quenching event, stellar ageing in QGs does not introduce significant scatter in f_dust. The situation drastically changes when QGs transition to the quenched stage. At the given stellar age, this evolutionary point unveils not only a significant scatter in f_dust, but also a division between the underlying quenching routes.

The lower left panel of Fig. 7 reveals three prominent regions. The first one is reserved for fast quenchers that occupy the upper left side of the plane, encompassing dustier objects with younger stellar populations (A). This region only partly overlaps with those of slow quenchers, who are characterised by generally older stellar populations (log(Age/yr) ≳ 9.4). Slow quenchers consist of two distinct clouds: one of slow, relatively dust-rich QGs (B; log(f_dust)≳ − 4), and one of slow dust-poor QGs (C; log(f_dust)≲ − 5). The positions of these three clouds change when QGs reach their final phase, characterised by a significant deficit of cold H₂ gas. While both group A and C end-up dust-poor and populate the lower-right side of the diagram, group B maintain an extremely shallow or almost absent evolution in f_dust, contrasting the expectations from age-evolutionary models. This trend can be seen in Appendix B, where galaxy G1, G4, G5 and G3 experience little evolution of f_dust across the evolution from P1 to P4 and opposite to the jet mode dominated G2 and G6, which experience a sudden decrease in dust content during the quenching phase ending up in group C.

To understand the diverse pathways for dust evolution amongst the three distinct QG groups mentioned above, we first checked their median quenching times and found them to be $t_{quench, A} \equiv 0 . 05_{- 0.02}^{+ 0.05}$ $t_{\mathrm{quench, A}} \equiv 0.05^{+0.05}_{-0.02}$ Gyr for fast quenchers, and $t_{quench, B} \equiv 0 . 82_{- 0.48}^{+ 0.98}$ $t_{\mathrm{quench, B}} \equiv 0.82^{+0.98}_{-0.48}$ Gyr and $t_{quench, C} \equiv 0 . 61_{- 0.28}^{+ 0.70}$ $t_{\mathrm{quench, C}} \equiv 0.61^{+0.70}_{-0.28}$ Gyr, for dust-rich and dust-poor slow quenchers, respectively. Values are presented in Table 1 together with other characterising properties of the sub-samples. Errors are assigned as the 16th and 84th percentile of each sub-sample distribution. The scale of the error is a direct consequence of the fast and slow QGs definition, which is not directly based on t_quench, but accounts for normalisation by the cosmic time at the epoch of quenching as explained in Sect. 4.1. This means that slow QGs can have both short and long t_quench if they are quenched at low-z or high-z respectively, spreading out the distribution.

Table 1.

Median roperties of SIMBA field galaxies at quenching.

If we focus on the median of the t_quench it appears that the division between fast and slow QGs is much clearer than a division between the groups B and C that compose the slow quenchers. Therefore, migration of the slow QGs to low f_dust does not depend solely on the duration of the quenching event but is largely affected by the dominating quenching mechanism. Both the A and C groups are mostly affected by a jet-dominated AGN, while the B group is mainly influenced by radiative mode feedback, as shown in Fig. 6.

Akins et al. (2022) couples SIMBA with POWDERDAY radiative transfer code (Narayanan et al. 2021), and tracks the main evolutionary routes among galaxies transitioning through the quenching phase and those that are recently quenched after an intense episode of starburst (so-called post-starbursts galaxies, PSBs). Interestingly, they find further movement across the UVJ colour plane independent of quenching timescales. As a result the quenching mode itself seems to be insufficient to characterise the ISM content of QGs that quenched at z ∼ 1 − 2 (see also Sun & Egami 2022; French et al. 2023). The C group from our sample show analogous characteristics to the PSBs described. Such galaxies have slightly higher f_H₂ with respect to other subsamples, $log (f_{H_{2}, B}) = - 1 . 37_{- 0.17}^{+ 0.17}$ $\log(f_{\mathrm{H_2, B}}) = -1.37^{+0.17}_{-0.17}$ , $log (f_{H_{2}, A}) = - 1 . 37_{- 0.55}^{+ 0.27}$ $\log(f_{\mathrm{H_2, A}}) =-1.37^{+0.27}_{-0.55}$ , and $log (f_{H_{2}, C}) = - 1 . 20_{- 0.14}^{+ 0.13}$ $\log(f_{\mathrm{H_2, C}}) = -1.20^{+0.13}_{-0.14}$ respectively.

A higher incidence rate of AGN feedback in QGs (and/or PSB galaxies) with younger stellar ages has recently been confirmed by several studies (Greene et al. 2020; Wu et al. 2023). Additionally, Lammers et al. (2023) explored MaNGA (Bundy et al. 2015) galaxies, and found that the radio-mode AGN feedback that significantly affects central SF may not always be efficient in expelling cold gas in a galaxy-wide fashion. They infer a long-lasting timescale for gas removal of ∼5 − 6 Gyr. The cold gas removal timescales obtained for rapidly quenched dusty QGs in this study are similarly long (see Table 1). This may support the idea that even for rapidly quenched QGs that are being affected by jet mode feedback, AGN efficiency is not directly translated to further dust evolution.

5.2.4. Dissecting the age-evolution of dusty QGs

The degeneracy in the overall evolution and in the large scatter observed in f_dust at the quenched phase is not resolved by considering the mass-weighted stellar age as presented in Fig. 7. Within the age-evolutionary scenario, both the dust and cold gas content of QGs are expected to follow each other and decline exponentially as stellar population ages, as proposed in several studies (e.g. Michałowski et al. 2019; Li et al. 2019b; Nanni et al. 2020; Lee et al. 2024). This evolution takes the functional form:

$\begin{matrix} f_{dust} = \frac{M_{dust}}{M_{⋆}} ≃ A \cdot \exp (- \frac{Age}{τ}), \end{matrix}$ $\begin{aligned} f_{\rm dust} = \frac{M_{\text{dust}}}{M_{\star }} \simeq A \cdot \mathrm{exp} \left(-\frac{\mathrm{Age}}{\tau }\right), \end{aligned}$ (3)

where A is a normalisation constant and τ is the e-folding time. The typical timescale of dust removal is still uncertain, with some observations proposing values of τ > 1 Gyr (Michałowski et al. 2019; Gobat et al. 2018; Lee et al. 2024) and others stating much faster timescales τ < 100 − 500 Myr (Williams et al. 2021; Whitaker et al. 2021b). For comparison, we show in Fig. 7 the evolution of f_dust with age adopting a timescale of τ ≡ 0.7 Gyr and one of τ ≡ 2.5 Gyr. While part of the distribution of SIMBA QGs falls into the region spanned by the expected timescales, a cloud of old QGs (log(age_⋆/yr) ≳ 9.4) maintains relatively large dust content (−3 < log f_dust) < − 4). These values are 3 − 5× larger than expected from the longest dust destruction timescale proposed by Michałowski et al. (2019). On top of this, the evolution of QGs that become dust-poor (log(f_dust)≲ − 5) after the quenching phase, appears to be much steeper than empirically predicted. These galaxies move almost vertically on the plane, from the quenching phase (P1) to the quenched phase (P2), in the time t_quench. This behaviour is clearly visible for galaxy G6 and G2 in Appendix B. The two galaxies in the example host a jet-dominated AGN during the quenching phase and their f_dust decreases rapidly. However, they differ substantially for t_removal^H₂, which in the case of G2 is two orders of magnitude longer and results in almost a two orders of magnitude lower f_dust at the point P4. Contours show the sample distribution for the fast and slow QGs, highlighting how the scatter generated during the quenching phase is intrinsic to the slow QGs. The evolution along the dust-age plane for the specific examples galaxies of Fig. 1 are shown and discussed in Appendix B.

The fact that the H₂ removal timescale can be short even for relatively dusty QGs, challenges the common idea that the dust and H₂ components of the ISM evolve alike. From Table 1 we see that mean values of f_H₂ are similar for all three sub-samples at the quenched stage, while the dust amount in QGs from group C is systematically reduced, resulting in average δ_DGR and f_dust being > 2 orders of magnitude lower than in other QGs. This implies that caution should be taken when inferring M_gas from M_dust, as highlighted in Whitaker et al. (2021a).

Finally, we compare the distribution of simulated galaxies with observations from Morishita et al. (2022), Gobat et al. (2022), Donevski et al. (2023) and Lee et al. (2024). The QGs from the first three studies reside at z ∼ 1.5 − 2 while QGs from the later study are at intermediate redshifts of 0.2 < z < 0.7. The spread of the observational points agrees with the large scatter in f_dust predicted by our SIMBA analysis. From the two lower panels in Fig. 7 we can reveal the possible nature of f_dust as well as the evolutionary stage of observed QGs. Individually detected QGs with younger stellar ages are caught after the completion of their quenching phase and are mostly consistent with the fast mode of quenching. The particularly prominent object from Fig. 7 is the ‘dustiest’ QG identified at z ∼ 2 by Lee et al. (2024). The position of this source within the f_dust-age plane coincides with group A QGs which suggests the role of AGN jet feedback. Despite this, the f_dust is still considerably high due to a brief and efficient episode of dust growth, while the cold gas still preserves detectable amounts over ∼ Gyr timescales. This interpretation matches remarkably well the observational characterisation inferred by Lee et al. (2024), that is, evidence for a fast quenching episode but fairly slow ISM removal (see also Kakimoto et al. 2024). The same authors point out that the excess of dust observed for their source is difficult to reproduce assuming stellar dust production, even when accounting for maximal SNe yields with no destruction. This implies a support of some additional dust re-formation processes. Switching to other observations of QGs with similarly young stellar populations we can notice substantially lower f_dust as compared to the QG from Lee et al. (2024). The QG from Morishita et al. (2022) has a dust temperature (∼ 28 K) warmer than in typical QGs. This is a consequence of an extremely strong AGN which is found to power the galaxy-scale outflows yet leaving a substantial content of warm gas long after quenching (Man et al. 2021). Such conditions may suppress the efficient growth which may further explain the drop in f_dust. Interestingly, the two most massive QGs (M_⋆ > 10¹¹ M_⊙) have the lowest f_dust and come from the sample of strongly lensed QGs analysed by Gobat et al. (2022). They are on a track consistent with the fast dust depletion, in agreement with the suggestion from Whitaker et al. (2021b). Nevertheless, caution should be taken for these objects, as their M_dust may be subject to larger uncertainties due to their extended sizes, as discussed in Gobat et al. (2022).

The majority of QGs from Donevski et al. (2023) occupy the region dominated by slow quenchers, but are not exclusively attributed to the quenched phase. As can be seen from the bottom-right panel of Fig. 7, the oldest QGs may be even entering the dry-phase devoid of significant cold gas. The expected gas-removal timescale of $t_{removal}^{H_{2}} ≲ 1 Gyr$ $t^{\mathrm{H}_2}_{\mathrm{removal}} \lesssim 1\, \rm Gyr$ predicted by SIMBA is in agreement with what Donevski et al. (2023) derived for early-type QGs as the amount of time elapsed between the quenching and the time of epoch of the observation. Their sample has been selected excluding the presence of X-ray bright AGN candidates, which supports the dominance of low-efficiency radiative feedback as seen for the galaxies in group B in this study.

5.3. Emergence of dusty QGs via prolonged dust growth: a new pathway for the late ISM evolution

The analysis from previous sections reveals the existence of a considerable number of QGs for which the dust evolution following the quenching event occurs independently of gas-removal timescales. The most striking examples are dust-rich QGs for which f_dust remains high even in a late phase characterised by a cold gas deficit. Following the Eq. (2), we infer the median values for the gas-removal timescales and find them to be $t_{removal, A}^{H_{2}} = 5 . 73_{- 3.00}^{+ 2.18}$ $t^{\mathrm{H}_2}_{\mathrm{removal, A}} = 5.73^{+2.18}_{-3.00}$ Gyr, $t_{removal, B}^{H_{2}} = 2 . 81_{- 1.66}^{+ 4.22}$ $t^{\mathrm{H}_2}_{\mathrm{removal, B}} = 2.81^{+4.22}_{-1.66}$ Gyr, $t_{removal, C}^{H_{2}} = 6 . 42_{- 2.54}^{+ 1.80}$ $t^{\mathrm{H}_2}_{\mathrm{removal, C}} = 6.42^{+1.80}_{-2.54}$ Gyr. From this we can deduce that in dusty QGs the gas-removal timescales are not tightly related to the dominant quenching mode, introducing new complexity in the pathways for the dust and gas (co-)evolution.

In QGs dominated by older stellar populations, there would be time for AGB stars to contribute to the dust mass budget. However, as pointed out in several works (Nanni et al. 2013; Hirashita et al. 2015; Ventura et al. 2020), the dust mass production efficiency decreases in metal-rich environments similar to our QGs. On top of this, the AGB condensation efficiency adopted in SIMBA models is quite low (< 20%), which would return the maximal M_dust order of ≲10⁶ M_⊙, unable to capture the higher values seen in our dusty QGs. This suggests an insufficient dust yield solely produced by AGB stars which prevents them from being a dominant channel of dust re-formation in QGs, in line with conclusions from several independent works (e.g. Hirashita et al. 2015; Michałowski et al. 2019).

To explain the excess of dust content of QGs, a handful of studies pointed to the importance of dust growth that would allow the dust grains to increase in size, thus counteracting the joint actions of AGN and thermal sputtering (Mattsson et al. 2014; Hirashita & Nozawa 2017; Richings & Faucher-Giguère 2018; Whitaker et al. 2021a; Donevski et al. 2023). This process is expected to be viable in metal-rich environments once the heated gas cools down. To probe this scenario in a quantitative way, we first checked a dust re-growth timescale (τ_acc) in QGs which we roughly derive following Li et al. (2019a) and using the same prescription used for SIMBA:

$\begin{matrix} \begin{matrix} τ_{acc} = & 1 \times 10^{7} yr \\ \times {(\frac{a}{0.1 μ m})}^{- 1} {(\frac{n_{H}}{100 {cm}^{- 3}})}^{- 1} \\ \times {(\frac{T_{gas}}{20 K})}^{- 1 / 2} {(\frac{Z_{gas}}{0.02})}^{- 1} [yr], \end{matrix} \end{matrix}$ $\begin{aligned} \begin{aligned} \tau _{\rm acc} =&\, 1\times 10^{7} \, \text{ yr} \\&\times \left(\frac{ a }{0.1\, \upmu \mathrm{m}}\right)^{-1} \left(\frac{ n _{\rm H}}{100\, \mathrm{cm}^{-3}}\right)^{-1} \\&\times \left(\frac{ T _{\rm gas}}{20\, \mathrm{K}}\right)^{-1/2} \left(\frac{ Z _{\rm gas}}{0.02}\right)^{-1} [\mathrm {yr}], \end{aligned} \end{aligned}$ (4)

where the first term defines the reference timing of growth activation, a is a dust grain size with a typical size of ∼0.1 μm (Li et al. 2019a). The n_H and T_gas are the number density and temperature of the ISM gas respectively, and Z_gas is a gas-phase metallicity³. Molecular gas number density is calculated from the surface gas density following Popping et al. (2017). We calculated the gas surface density using the gas mass and the half-mass gas radius from SIMBA together with the metallicity. Grain size is assumed to be the size of newly formed grain in SIMBA, that is a ≡ 0.1 μm. We assume T_gas ≡ 30 K and we use Z_gas as provided by the simulation. We find that dust re-growth time is the most rapid amongst the slowly quenched dust-rich QGs of group B ( $τ_{acc} = 70 . 64_{- 42.2}^{+ 66.98}$ $\tau_{\mathrm{acc}}=70.64^{+66.98}_{-42.2}$ Myr). For the fast QGs of group A and the slow dust-poor QGs of group C, typical re-growth timescales are slightly longer on average.

The re-growth timescales obtained for all three groups of QGs are still sufficiently quick to compete against the sputtering timescales estimated for our dusty QGs (≳ 10⁸ yr). It is feasible for dust grains to re-form after the swept-up ISM has cooled. As pointed out by Hirashita & Nozawa (2017), as long as sputtering timescales are longer than the AGN cycle, and with short-enough re-growth timescales, there would be enough residual material to reform the dust on ISM metals. Donevski et al. (2023) propose the grain growth mechanism as a plausible explanation for the observed different pathways of early-type QGs, for which more than 500 Myr elapsed after the quenching. Timescales for such process counteracting the dust destruction are expected to be short, of the order of 50 − 250 Myr (Hirashita et al. 2015). Since the timescales for H₂ gas removal in group B are similar to that for dust destruction, the presence of large f_dust in such galaxies at the dry-phase strongly supports the idea of dust re-growth happening in < 500 μyr. Ni et al. (2023) investigates the relation between the observed fraction of X-dominated AGNs and their stellar age in the COSMOS field. They find a strong anti-correlation, with the oldest objects having lower AGN accretion rates. They conclude that a significant fraction of the accretion on the AGN is due to hot recycled gas, which cools and affects the AGN on a timescale of ∼1 Gyr. This means that, after the reduction of H₂ gas from the quenching, the ISM would have enough time before the re-activation of an intense X-ray AGN feedback to undergo dust re-growth, in qualitative agreement with SIMBA predictions.

It is important to note that prolonged dust growth is not the sole process influencing ISM dust abundance and it does not rule out the contribution of AGB stars. Dust growth through accretion on ISM metals requires seed grains to be deposited into the ISM for which AGB stars are likely sources (see e.g. Nanni et al. 2013). An alternative possibility, as proposed in several observational works, suggests dust mixing via the interaction of AGN winds with star formation-driven outflows. The latter effect seems to be particularly relevant for lower-mass galaxies (see e.g. Romano et al. 2023). Investigating this possibility is out of the scope of this study.

We recover the same dust-poor but gas-rich QGs (group C) characterised with exotically low δ_DGR that are first identified by Whitaker et al. (2021a). The presence of such QGs can be interpreted by means of dilution and inefficient dust grain growth, as already proposed by Hirashita & Nozawa (2017). In this scenario, during long-lasting (∼5 − 6 Gyr) timescales, radiative mode AGN feedback steadily removes the H₂ gas, and dust grains become efficiently destroyed by thermal sputtering and do not have a way to withstand the erosion by growing through condensation of metals. On top of this, inflows of dust-poor cold gas would decrease the δ_DGR, keeping the group C sources in the lower part of Fig. 5. It is worth noting that for the case of inefficient dust growth, their δ_DGR becomes highly sensitive to the M_⋆, and strongly affects most massive QGs, as we shown in Section 4.3.

A direct test of grain growth in SIMBA can be performed by re-running the simulation without the prescription for grain growth. However, the exclusion of this mechanisms would significantly impact the entire ISM physics even during the SF phase of the galaxies, thus producing entirely different evolutionary paths and resulting QGs, making the comparison significantly imprecise. A similar test is performed in Whitaker et al. (2021a), where the removal of dust growth in SIMBA results in a shift of the entire galaxy sample towards lower δ_DGR. SFGs shows a larger effect, with δ_DGR decreasing by > 1 dex, while for QGs it shifts by ∼ 0.5 dex on average. We leave a careful investigation of the dust growth in an accompanying work (Lorenzon et al. in prep.) by exploring the high-resolution run of the simulation and accounting for the full dust cycle connecting ISM and CGM. We choose here to monitor the changes of the ISM in the quenched phase by means of the most direct marker available for observations, that is the dust-to-metal-ratio (δ_DTM), which is proven to be a great proxy for the dust re-growth (see e.g. De Vis et al. 2019).

We follow De Vis et al. (2019) and calculate δ_DTM ≡ M_d/M_Z where M_Z = f_z × M_g + M_d is the total mass of metals, M_d is the dust mass, M_g is the gas mass and $f_{Z} = 27.36 \times 10^{12 + \log_{10} (O / H) - 12}$ $f_{\mathrm{Z}}=27.36 \times 10^{12+\mathrm{log}_{10}(\mathrm{O/H})-12}$ is the mass fraction of metals. Values are listed in Table 1. We find that at the end of the quenching phase, both groups A and B have large values of log(δ_DTM) > − 0.8, while for group C the value is at least twice as big. The values for groups A and B are comparable with the literature reference log(δ_DTM)∼ − 0.52 (i.e. Camps et al. 2015; Clark et al. 2016) and the values from DustPedia galaxies (log(δ_DTM)∼ − 0.67; De Vis et al. 2019) measured for QGs. At the gas-removed evolutionary point the scenario is rather different, with median of $log (δ_{DTM, A}) = - 2 . 82_{- 0.81}^{+ 2.46}$ $\log(\delta_{\mathrm{DTM, A}}) = -2.82^{+2.46}_{-0.81}$ , $log (δ_{DTM, B}) = - 0 . 27_{- 1.16}^{+ 0.1}$ $\log(\delta_{\mathrm{DTM, B}}) = -0.27^{+0.1}_{-1.16}$ , $log (δ_{DTM, C}) = - 3 . 14_{- 0.69}^{+ 0.89}$ $\log(\delta_{\mathrm{DTM, C}}) = -3.14^{+0.89}_{-0.69}$ . While there is a significant decrease in δ_DTM for both groups A and C, the value for group B stays rather high up to the dry-phase, indicating the presence of dust growth in a large part of the sample. We qualitatively summarise our results in Fig. 8, where the properties of the groups A, B, and C are shown in relation to f_dust and stellar age, together with the direction of evolution and dust growth.

Fig. 8.

Qualitative sketch explaining the evolution of f_dust with stellar age within the quenched phase, as presented in the bottom panel of Fig. 7. The cyan cloud illustrates the position of fast quenchers, while the two red clouds illustrate the position of dust-poor and dust-rich slow quenchers. The evolutionary directions (marked with arrows) for each sample are shown along with the main physical mechanisms involved in their dust evolution. The meaning of lines in the large panel is the same as in Fig. 7. The dark-grey shaded regions in the two subplots from the right panel delineate the part of the age-dust plane devoid of the given QG.

5.4. Observational implications for studying high-z dusty QGs

While it has been commonly assumed that QGs are dust-poor, the presence of significant dust content in SIMBA QGs predicted by this study imposes several important implications. Firstly, the substantial ISM dust abundance can severally affect the intrinsic SEDs of QGs and selection techniques based on rest-frame colours (see e.g. Akins et al. 2022; Long et al. 2023; Setton et al. 2024). Secondly, we show that observed dust and cold gas tightly follow each other until the end of a quenching phase. After the galaxy becomes fully quenched, the observed f_dust may be resulting from a dust growth, with re-formation timescales which may be offsetting from those attributed to H₂ gas exhaustion. For that reason, we urge caution when using f_dust to convert to f_gas and ultimately to M_gas in QGs. This seems to be particularly important for planning the ALMA observations to probe the cold gas content via commonly used tracers such as CO transitions. At a given MS-offset (or stellar age), older and metal-richer QGs with enhanced f_dust would require longer integration times for the detection of molecular gas than those simply extrapolated from the typical δ_DGR ∼ 1/100. This is a consequence of enhanced δ_DGR and subsequently lower M_gas. Indeed, recent studies pointed out alternative tracers such as [CII] 158 μm emission line to be more efficient than CO for probing a gaseous component of QGs at z > 2 (D’Eugenio et al. 2023a). We thoroughly investigate the two above-mentioned challenges in our accompanying study (Lorenzon et al., in prep.). Finally, our findings suggest the existence of dust-rich QGs at z ∼ 1 − 2 that were quenched > 5 Gyr ago (z_q ≳ 3 − 5) through the rapid and intense AGN activity, similar to those that were confirmed observationally with JWST (Stacey et al. 2022; Kocevski et al. 2023; D’Eugenio et al. 2023b) and with ALMA (Salak et al. 2024). Joint observational strategies between JWST and ALMA are required to test this specific prediction. If confirmed, such rapidly quenched but dust-rich galaxies will provide a unique insight into the SMBH-ISM co-evolution in the evolved objects caught in the first 1–2 Gyr of the Universe’s history.

6. Conclusions

We use the full suite of SIMBA cosmological simulation and provide the theoretical framework for the pathways of dust and cold gas evolution in QGs at high-z. We identify a mass-complete sample (M_⋆ > 10⁹ M_⊙) of QGs that are quenched over a wide redshift range (0 < z ≲ 2) and study their ISM dust and gas abundance at different critical points in evolutionary history. We account for their large-scale environments, exploring separately the QGs residing in the field and in clusters.

Our main results are summarised as follows:

Quenching timescales in SIMBA are strongly bimodal, and no prevalence of a specific quenching mode has been found for a given environment. The bimodality does not translate directly into dust-rich and dust-poor galaxies, suggesting that the quenching mode alone does not have a prevalent effect on the ISM evolution.
The δ_DGR-Z_gas plane for QGs consists of two distinct regions. Roughly half of QGs follow the relation expected for SFGs and the other half exhibits the absence of a clear trend accompanied by a large scatter (∼4 orders of magnitude) at super-solar Z_gas. Such behaviour suggests the diversity of efficiencies for the production and removal of dust and gas in QGs. Although the overall shape remains the same for the field and cluster QGs, we find that dust destruction is enhanced by sputtering in cluster environments.
The stellar masses in the quenching phase have a profound impact on further dust or gas evolution in QGs. As they depart from the MS towards the final gas-dry phase (f_gas < 0.1%), more than a half of QGs with intermediate masses (log(M_⋆/M_⊙)≲10.5) sustain high ISM dust abundance (f_dust > 10⁻³ − 10⁻⁴; δ_DGR > 1/100). This is in contrast to massive QGs (10.5 < log(M_⋆/M_⊙) < 11.25) that experience a big drop (> 1–2 orders of magnitude) in their f_dust and δ_DGR.
Enhanced f_dust and δ_DGR are mostly a consequence of internal processes, and only marginally of external effects such as mergers and rejuvenations. The dominant dust re-formation mechanism is prolonged dust growth in metal-rich ISM. This key process is effective in the first 100 − 500 Myr of a quenching phase, allowing QGs to extend the dust destruction timescales expected for thermal gas heating.
We successfully reproduce the observational relation of f_dust with stellar age in QGs identified at 0 < z ≲ 2. We interpret the emergence of these sources as a mix of both fast and slow quenchers, with fast quenchers predominantly having younger stellar ages (< 3 Gyr) and slow quenchers dominating the large scatter towards older stellar ages (> 3–10 Gyr). During the quenching phases, the bulk of QGs shows little to no evolution in their f_dust with age. This implies long dust depletion timescales (> 2–3 Gyr), a direct consequence of prolonged grain growth.
Overall, we predict that for a significant number of QGs, removal timescales for the cold H₂ gas and dust are mutually different, owing to a complex interplay of quenching pathways attributed to the AGN, and various dust re-growth efficiencies. This prediction may have important implications for characterising the SMBH-ISM co-evolution in QGs.

We stress that the aforementioned predictions urge immediate observational investigations of dusty QGs. In particular, to better understand the strength of the AGN activity and its coupling to the ISM dust and gas, our study suggests a systematic search and analysis of ISM in high-z QGs of low and intermediate masses (9 < log(M_⋆/M_⊙) < 10.5). Although challenging, this observational quest is achievable with the current generation of instruments such as JWST, ALMA and Euclid. The JWST MIRI, in particular, can be used to probe the warm dust and molecular gas through MIR diagnostics in individual dusty QGs at various cosmic epochs up to z ∼ 2 − 3, while its combination with NIRspec will help access the level of AGN activity and better constrain gas-phase metallicity. Synergy with sub-mm instruments (ALMA, NOEMA) and near-IR survey telescopes (Euclid) is needed to infer the census on the cold dust and large-scale environments of selected QGs, respectively. In our complementary work, we provide the recipes for such observational tests and present the key differences in SEDs of dust-rich versus dust-poor QGs that were formed at different cosmic epochs (Lorenzon et al., in prep.).

¹

The SF threshold roughly corresponds to the evolution of the level of the sSFR in the MS minus 0.3 dex, which is generally considered to be the lower limit of the MS scatter. The quenching threshold roughly corresponds to the average observed distance of at least 0.4 dex below the MS, ensuring that the galaxy can be considered quenched independently of its mass.

²

For this task, we homogenise the MS functional form from Speagle et al. (2014) to the IMF by Chabrier (2003).

³

Note that given the SIMBA’s approximate kiloparsec resolution, a multi-phase galaxy ISM is not resolved, which prevent us from knowing the exact dependence of the modelled dust content to the gas surface density.

Acknowledgments

G.L. and D.D. acknowledge support from the National Science centre (NCN) grant SONATA (UMO-2020/39/D/ST9/00720). K.L acknowledges support of the Polish Ministry of Education and Science through the grant PN/01/0034/2022 under ‘Perły Nauki’ programme. AWSM acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) through grant reference number RGPIN-2021-03046. J is grateful for support from the Polish National Science Centre via grant UMO-2018/30/E/ST9/00082. A.N, and M.R. acknowledge support from the Narodowe Centrum Nauki (NCN), Poland, through the SONATA BIS grant UMO2020/38/E/ST9/00077. J. is funded by the European Union (MSCA EDUCADO, GA 101119830 and WIDERA ExGal-Twin, GA 101158446). M.R. acknowledges support from the Foundation for Polish Science.

References

Akhshik, M., Whitaker, K. E., Leja, J., et al. 2021, ApJ, 907, L8 [NASA ADS] [CrossRef] [Google Scholar]
Akhshik, M., Whitaker, K. E., Leja, J., et al. 2023, ApJ, 943, 179 [NASA ADS] [CrossRef] [Google Scholar]
Akins, H. B., Narayanan, D., Whitaker, K. E., et al. 2022, ApJ, 929, 94 [NASA ADS] [CrossRef] [Google Scholar]
Alarcon, A., Hearin, A. P., Becker, M. R., & Chaves-Montero, J. 2023, MNRAS, 518, 562 [Google Scholar]
Aoyama, K., Kodama, T., Suzuki, T. L., et al. 2022, ApJ, 924, 74 [NASA ADS] [CrossRef] [Google Scholar]
Appleby, S., Davé, R., Kraljic, K., Anglés-Alcázar, D., & Narayanan, D. 2020, MNRAS, 494, 6053 [NASA ADS] [CrossRef] [Google Scholar]
Asano, R. S., Takeuchi, T. T., Hirashita, H., & Inoue, A. K. 2013, Earth Planets Space, 65, 213 [NASA ADS] [CrossRef] [Google Scholar]
Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
Barišić, I., van der Wel, A., Bezanson, R., et al. 2017, ApJ, 847, 72 [Google Scholar]
Belli, S., Genzel, R., Förster Schreiber, N. M., et al. 2017, ApJ, 841, L6 [NASA ADS] [CrossRef] [Google Scholar]
Belli, S., Newman, A. B., & Ellis, R. S. 2019, ApJ, 874, 17 [Google Scholar]
Belli, S., Contursi, A., Genzel, R., et al. 2021, ApJ, 909, L11 [NASA ADS] [CrossRef] [Google Scholar]
Béthermin, M., Daddi, E., Magdis, G., et al. 2015, A&A, 573, A113 [Google Scholar]
Bianchi, S., & Schneider, R. 2007, MNRAS, 378, 973 [NASA ADS] [CrossRef] [Google Scholar]
Blánquez-Sesé, D., Gómez-Guijarro, C., Magdis, G. E., et al. 2023, A&A, 674, A166 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Boselli, A., Fossati, M., Longobardi, A., et al. 2020, A&A, 634, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bundy, K., Fukugita, M., Ellis, R. S., et al. 2009, ApJ, 697, 1369 [NASA ADS] [CrossRef] [Google Scholar]
Bundy, K., Bershady, M. A., Law, D. R., et al. 2015, ApJ, 798, 7 [Google Scholar]
Burgarella, D., Nanni, A., Hirashita, H., et al. 2020, A&A, 637, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Caliendo, J. N., Whitaker, K. E., Akhshik, M., et al. 2021, ApJ, 910, L7 [NASA ADS] [CrossRef] [Google Scholar]
Calura, F., Pozzi, F., Cresci, G., et al. 2017, MNRAS, 465, 54 [Google Scholar]
Camps, P., Misselt, K., Bianchi, S., et al. 2015, A&A, 580, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Carnall, A. C., McLeod, D. J., McLure, R. J., et al. 2023, MNRAS, 520, 3974 [NASA ADS] [CrossRef] [Google Scholar]
Casasola, V., Bianchi, S., De Vis, P., et al. 2020, A&A, 633, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
Chauke, P., van der Wel, A., Pacifici, C., et al. 2018, ApJ, 861, 13 [NASA ADS] [CrossRef] [Google Scholar]
Chauke, P., van der Wel, A., Pacifici, C., et al. 2019, ApJ, 877, 48 [NASA ADS] [CrossRef] [Google Scholar]
Chiang, Y.-K., Overzier, R., & Gebhardt, K. 2013, ApJ, 779, 127 [Google Scholar]
Choi, E., Ostriker, J. P., Naab, T., & Johansson, P. H. 2012, ApJ, 754, 125 [NASA ADS] [CrossRef] [Google Scholar]
Clark, C. J. R., Schofield, S. P., Gomez, H. L., & Davies, J. I. 2016, MNRAS, 459, 1646 [NASA ADS] [CrossRef] [Google Scholar]
Corcho-Caballero, P., Ascasibar, Y., Sánchez, S. F., & López-Sánchez, Á. R. 2023, MNRAS, 520, 193 [NASA ADS] [CrossRef] [Google Scholar]
Crenshaw, D. M., Kraemer, S. B., & George, I. M. 2003, ARA&A, 41, 117 [NASA ADS] [CrossRef] [Google Scholar]
Cuppen, H. M., Walsh, C., Lamberts, T., et al. 2017, Space Sci. Rev., 212, 1 [Google Scholar]
Daddi, E., Renzini, A., Pirzkal, N., et al. 2005, ApJ, 626, 680 [NASA ADS] [CrossRef] [Google Scholar]
Dariush, A., Dib, S., Hony, S., et al. 2016, MNRAS, 456, 2221 [NASA ADS] [CrossRef] [Google Scholar]
Davé, R., Thompson, R., & Hopkins, P. F. 2016, MNRAS, 462, 3265 [Google Scholar]
Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827 [Google Scholar]
Davis, T. A., Greene, J., Ma, C.-P., et al. 2016, MNRAS, 455, 214 [Google Scholar]
De Vis, P., Jones, A., Viaene, S., et al. 2019, A&A, 623, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
D’Eugenio, C., Daddi, E., Liu, D., & Gobat, R. 2023a, A&A, 678, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
D’Eugenio, F., Perez-Gonzalez, P., Maiolino, R., et al. 2023b, Nat. Astron., submitted, [arXiv:2308.06317] [Google Scholar]
Dierckx, P. 1975, J. Comput. Appl. Math., 1, 165 [CrossRef] [Google Scholar]
Donas, J., Deharveng, J.-M., Rich, R. M., et al. 2007, ApJS, 173, 597 [NASA ADS] [CrossRef] [Google Scholar]
Donevski, D., Lapi, A., Małek, K., et al. 2020, A&A, 644, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Donevski, D., Damjanov, I., Nanni, A., et al. 2023, A&A, 678, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Donnari, M., Pillepich, A., Joshi, G. D., et al. 2021, MNRAS, 500, 4004 [Google Scholar]
Dunne, L., Gomez, H. L., da Cunha, E., et al. 2011, MNRAS, 417, 1510 [NASA ADS] [CrossRef] [Google Scholar]
Dwek, E. 1998, ApJ, 501, 643 [NASA ADS] [CrossRef] [Google Scholar]
Faber, S. M., Willmer, C. N. A., Wolf, C., et al. 2007, ApJ, 665, 265 [Google Scholar]
Feldmann, R., & Mayer, L. 2015, MNRAS, 446, 1939 [Google Scholar]
Ferrarotti, A. S., & Gail, H. P. 2006, A&A, 447, 553 [CrossRef] [EDP Sciences] [Google Scholar]
Foltz, R., Wilson, G., Muzzin, A., et al. 2018, ApJ, 866, 136 [NASA ADS] [CrossRef] [Google Scholar]
French, K. D., Smercina, A., Rowlands, K., et al. 2023, ApJ, 942, 25 [NASA ADS] [CrossRef] [Google Scholar]
Fumagalli, M., Labbé, I., Patel, S. G., et al. 2014, ApJ, 796, 35 [Google Scholar]
Galliano, F., Nersesian, A., Bianchi, S., et al. 2021, A&A, 649, A18 [EDP Sciences] [Google Scholar]
Ghodsi, L., Man, A. W. S., Donevski, D., et al. 2024, MNRAS, 528, 4393 [NASA ADS] [CrossRef] [Google Scholar]
Gillman, S., Tiley, A. L., Swinbank, A. M., et al. 2021, MNRAS, 500, 4229 [NASA ADS] [Google Scholar]
Gobat, R., Daddi, E., Magdis, G., et al. 2018, Nat. Astron., 2, 239 [NASA ADS] [CrossRef] [Google Scholar]
Gobat, R., Magdis, G., D’Eugenio, C., & Valentino, F. 2020, A&A, 644, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Gobat, R., D’Eugenio, C., Liu, D., et al. 2022, A&A, 668, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Greene, J. E., Setton, D., Bezanson, R., et al. 2020, ApJ, 899, L9 [NASA ADS] [CrossRef] [Google Scholar]
Gruppioni, C., Béthermin, M., Loiacono, F., et al. 2020, A&A, 643, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Gunn, J. E., Gott, J., & Richard, I. 1972, ApJ, 176, 1 [Google Scholar]
Hamadouche, M. L., Carnall, A. C., McLure, R. J., et al. 2023, MNRAS, 521, 5400 [NASA ADS] [CrossRef] [Google Scholar]
Heckman, T. M., & Best, P. N. 2014, ARA&A, 52, 589 [Google Scholar]
Hirashita, H., & Nozawa, T. 2017, Planet Space Sci., 149, 45 [NASA ADS] [CrossRef] [Google Scholar]
Hirashita, H., Nozawa, T., Villaume, A., & Srinivasan, S. 2015, MNRAS, 454, 1620 [NASA ADS] [CrossRef] [Google Scholar]
Hodge, J. A., & da Cunha, E. 2020, Roy. Soc. Open Sci., 7, 200556 [NASA ADS] [CrossRef] [Google Scholar]
Hopkins, P. F. 2015, MNRAS, 450, 53 [Google Scholar]
Hough, R. T., Rennehan, D., Kobayashi, C., et al. 2023, MNRAS, 525, 1061 [NASA ADS] [CrossRef] [Google Scholar]
Inami, H., Algera, H. S. B., Schouws, S., et al. 2022, MNRAS, 515, 3126 [NASA ADS] [CrossRef] [Google Scholar]
Iwamoto, K., Brachwitz, F., Nomoto, K., et al. 1999, ApJS, 125, 439 [NASA ADS] [CrossRef] [Google Scholar]
Junais, Boissier, S., Boselli, A., et al. 2022, A&A, 667, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Kakimoto, T., Tanaka, M., Onodera, M., et al. 2024, ApJ, 963, 49 [NASA ADS] [CrossRef] [Google Scholar]
Kalita, B. S., Daddi, E., D’Eugenio, C., et al. 2021, ApJ, 917, L17 [NASA ADS] [CrossRef] [Google Scholar]
Kannan, R., Vogelsberger, M., Marinacci, F., et al. 2021, MNRAS, 503, 336 [NASA ADS] [CrossRef] [Google Scholar]
Kaviraj, S. 2014, MNRAS, 437, L41 [NASA ADS] [CrossRef] [Google Scholar]
Kennicutt, R. C., Jr 1998, ApJ, 498, 541 [Google Scholar]
Kocevski, D. D., Onoue, M., Inayoshi, K., et al. 2023, ApJ, 954, L4 [NASA ADS] [CrossRef] [Google Scholar]
Kokorev, V. I., Magdis, G. E., Davidzon, I., et al. 2021, ApJ, 921, 40 [NASA ADS] [CrossRef] [Google Scholar]
Kokusho, T., Kaneda, H., Bureau, M., et al. 2019, A&A, 622, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Kriek, M., van der Wel, A., van Dokkum, P. G., Franx, M., & Illingworth, G. D. 2008, ApJ, 682, 896 [NASA ADS] [CrossRef] [Google Scholar]
Krumholz, M. R., McKee, C. F., & Tumlinson, J. 2009a, ApJ, 693, 216 [NASA ADS] [CrossRef] [Google Scholar]
Krumholz, M. R., McKee, C. F., & Tumlinson, J. 2009b, ApJ, 699, 850 [NASA ADS] [CrossRef] [Google Scholar]
Lammers, C., Iyer, K. G., Ibarra-Medel, H., et al. 2023, ApJ, 953, 26 [NASA ADS] [CrossRef] [Google Scholar]
Lee, M. M., Steidel, C. C., Brammer, G., et al. 2024, MNRAS, 527, 9529 [Google Scholar]
Leroy, A. K., Bolatto, A., Gordon, K., et al. 2011, ApJ, 737, 12 [Google Scholar]
Li, Q., Narayanan, D., & Davé, R. 2019a, MNRAS, 490, 1425 [CrossRef] [Google Scholar]
Li, Z., French, K. D., Zabludoff, A. I., & Ho, L. C. 2019b, ApJ, 879, 131 [NASA ADS] [CrossRef] [Google Scholar]
Lianou, S., Xilouris, E., Madden, S. C., & Barmby, P. 2016, MNRAS, 461, 2856 [NASA ADS] [CrossRef] [Google Scholar]
Lisenfeld, U., Isaak, K. G., & Hills, R. 2000, MNRAS, 312, 433 [CrossRef] [Google Scholar]
Liu, D., Schinnerer, E., Groves, B., et al. 2019, ApJ, 887, 235 [Google Scholar]
Long, A. S., Antwi-Danso, J., Lambrides, E. L., et al. 2023, ApJ, submitted, [arXiv:2305.04662] [Google Scholar]
Looser, T. J., D’Eugenio, F., Maiolino, R., et al. 2024, Nature, 629, 53 [Google Scholar]
Lovell, C. C., Thomas, P. A., & Wilkins, S. M. 2018, MNRAS, 474, 4612 [NASA ADS] [CrossRef] [Google Scholar]
Lovell, C. C., Geach, J. E., Davé, R., Narayanan, D., & Li, Q. 2021, MNRAS, 502, 772 [NASA ADS] [CrossRef] [Google Scholar]
Ma, W., Liu, K., Guo, H., et al. 2022, ApJ, 941, 205 [NASA ADS] [CrossRef] [Google Scholar]
Magdis, G. E., Daddi, E., Elbaz, D., et al. 2011, ApJ, 740, L15 [NASA ADS] [CrossRef] [Google Scholar]
Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012, ApJ, 760, 6 [NASA ADS] [CrossRef] [Google Scholar]
Magdis, G. E., Gobat, R., Valentino, F., et al. 2021, A&A, 647, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Maier, C., Ziegler, B. L., Haines, C. P., & Smith, G. P. 2019, A&A, 621, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Man, A. W. S., Greve, T. R., Toft, S., et al. 2016a, ApJ, 820, 11 [NASA ADS] [CrossRef] [Google Scholar]
Man, A. W. S., Zirm, A. W., & Toft, S. 2016b, ApJ, 830, 89 [Google Scholar]
Man, A. W. S., Zabl, J., Brammer, G. B., et al. 2021, ApJ, 919, 20 [NASA ADS] [CrossRef] [Google Scholar]
Martini, P., Dicken, D., & Storchi-Bergmann, T. 2013, ApJ, 766, 121 [NASA ADS] [CrossRef] [Google Scholar]
Matharu, J., Muzzin, A., Brammer, G. B., et al. 2019, MNRAS, 484, 595 [Google Scholar]
Mattsson, L., Gomez, H. L., Andersen, A. C., et al. 2014, MNRAS, 444, 797 [Google Scholar]
McCormick, A., Veilleux, S., Meléndez, M., et al. 2018, MNRAS, 477, 699 [NASA ADS] [CrossRef] [Google Scholar]
McKinnon, R., Torrey, P., Vogelsberger, M., Hayward, C. C., & Marinacci, F. 2017, MNRAS, 468, 1505 [NASA ADS] [CrossRef] [Google Scholar]
Merritt, D. 1983, ApJ, 264, 24 [NASA ADS] [CrossRef] [Google Scholar]
Michałowski, M. J., Hjorth, J., Gall, C., et al. 2019, A&A, 632, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Michałowski, M. J., Gall, C., Hjorth, J., et al. 2024, ApJ, 964, 129 [CrossRef] [Google Scholar]
Millard, J. S., Eales, S. A., Smith, M. W. L., et al. 2020, MNRAS, 494, 293 [NASA ADS] [CrossRef] [Google Scholar]
Moore, B., Lake, G., & Katz, N. 1998, ApJ, 495, 139 [Google Scholar]
Morishita, T., Abdurro’uf, Hirashita, H., et al. 2022, ApJ, 938, 144 [Google Scholar]
Nanni, A., Bressan, A., Marigo, P., & Girardi, L. 2013, MNRAS, 434, 2390 [NASA ADS] [CrossRef] [Google Scholar]
Nanni, A., Burgarella, D., Theulé, P., Côté, B., & Hirashita, H. 2020, A&A, 641, A168 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Narayanan, D., Turk, M. J., Robitaille, T., et al. 2021, ApJS, 252, 12 [NASA ADS] [CrossRef] [Google Scholar]
Ni, Q., Aird, J., Merloni, A., et al. 2023, MNRAS, 524, 4778 [CrossRef] [Google Scholar]
Nozawa, T., Kozasa, T., & Habe, A. 2006, ApJ, 648, 435 [NASA ADS] [CrossRef] [Google Scholar]
Oppenheimer, B. D., & Davé, R. 2008, MNRAS, 387, 577 [NASA ADS] [CrossRef] [Google Scholar]
Pacifici, C., Kassin, S. A., Weiner, B. J., et al. 2016, ApJ, 832, 79 [Google Scholar]
Pandya, V., Brennan, R., Somerville, R. S., et al. 2017, MNRAS, 472, 2054 [NASA ADS] [CrossRef] [Google Scholar]
Pantoni, L., Lapi, A., Massardi, M., Goswami, S., & Danese, L. 2019, ApJ, 880, 129 [NASA ADS] [CrossRef] [Google Scholar]
Park, M., Belli, S., Conroy, C., et al. 2023, ApJ, 953, 119 [NASA ADS] [CrossRef] [Google Scholar]
Peng, Y.-J., Lilly, S. J., Kovač, K., et al. 2010, ApJ, 721, 193 [Google Scholar]
Planck Collaboration XLIII. 2016, A&A, 596, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Popping, G., & Péroux, C. 2022, MNRAS, 513, 1531 [CrossRef] [Google Scholar]
Popping, G., Somerville, R. S., & Galametz, M. 2017, MNRAS, 471, 3152 [NASA ADS] [CrossRef] [Google Scholar]
Popping, G., Shivaei, I., Sanders, R. L., et al. 2023, A&A, 670, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Priestley, F. D., Chawner, H., Barlow, M. J., et al. 2022, MNRAS, 516, 2314 [NASA ADS] [CrossRef] [Google Scholar]
Relaño, M., De Looze, I., Saintonge, A., et al. 2022, MNRAS, 515, 5306 [CrossRef] [Google Scholar]
Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, A31 [Google Scholar]
Richie, H. M., Schneider, E. E., Abruzzo, M. W., & Torrey, P. 2024, ApJ, 974, 81 [NASA ADS] [CrossRef] [Google Scholar]
Richings, A. J., & Faucher-Giguère, C.-A. 2018, MNRAS, 474, 3673 [Google Scholar]
Robson, D., & Davé, R. 2020, MNRAS, 498, 3061 [NASA ADS] [CrossRef] [Google Scholar]
Rodighiero, G., Daddi, E., Baronchelli, I., et al. 2011, ApJ, 739, L40 [Google Scholar]
Rodighiero, G., Bisigello, L., Iani, E., et al. 2023, MNRAS, 518, L19 [Google Scholar]
Rodríguez Montero, F., Davé, R., Wild, V., Anglés-Alcázar, D., & Narayanan, D. 2019, MNRAS, 490, 2139 [CrossRef] [Google Scholar]
Romano, M., Nanni, A., Donevski, D., et al. 2023, A&A, 677, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Romano, M., Donevski, D., Junais, et al. 2024, A&A, 683, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Rowlands, K., Dunne, L., Maddox, S., et al. 2012, MNRAS, 419, 2545 [NASA ADS] [CrossRef] [Google Scholar]
Rowlands, K., Gomez, H. L., Dunne, L., et al. 2014, MNRAS, 441, 1040 [Google Scholar]
Rowlands, K., Wild, V., Bourne, N., et al. 2018, MNRAS, 473, 1168 [NASA ADS] [CrossRef] [Google Scholar]
Salak, D., Hashimoto, T., Inoue, A. K., et al. 2024, ApJ, 962, 1 [NASA ADS] [CrossRef] [Google Scholar]
Sandstrom, K. M., Leroy, A. K., Walter, F., et al. 2013, ApJ, 777, 5 [Google Scholar]
Sargent, M. T., Daddi, E., Bournaud, F., et al. 2015, ApJ, 806, L20 [NASA ADS] [CrossRef] [Google Scholar]
Schawinski, K., Kaviraj, S., Khochfar, S., et al. 2007, ApJS, 173, 512 [NASA ADS] [CrossRef] [Google Scholar]
Schmidt, M. 1959, ApJ, 129, 243 [NASA ADS] [CrossRef] [Google Scholar]
Schneider, E. E., & Robertson, B. E. 2015, ApJS, 217, 24 [NASA ADS] [CrossRef] [Google Scholar]
Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Setton, D. J., Khullar, G., & Miller, T. B., et al. 2024, ApJ, submitted [arXiv:2402.05664] [Google Scholar]
Shapley, A. E., Cullen, F., Dunlop, J. S., et al. 2020, ApJ, 903, L16 [NASA ADS] [CrossRef] [Google Scholar]
Socolovsky, M., Almaini, O., Hatch, N. A., et al. 2018, MNRAS, 476, 1242 [NASA ADS] [CrossRef] [Google Scholar]
Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15 [Google Scholar]
Spilker, J., Bezanson, R., Barišić, I., et al. 2018, ApJ, 860, 103 [NASA ADS] [CrossRef] [Google Scholar]
Springel, V. 2010, MNRAS, 401, 791 [Google Scholar]
Stacey, H. R., Costa, T., McKean, J. P., et al. 2022, MNRAS, 517, 3377 [NASA ADS] [CrossRef] [Google Scholar]
Sun, F., & Egami, E. 2022, MNRAS, 517, L126 [NASA ADS] [CrossRef] [Google Scholar]
Szpila, J., Davé, R., Rennehan, D., Cui, W., & Hough, R. 2024, MNRAS, submitted, [arXiv:2402.08729] [Google Scholar]
Tacchella, S., Conroy, C., Faber, S. M., et al. 2022, ApJ, 926, 134 [NASA ADS] [CrossRef] [Google Scholar]
Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179 [Google Scholar]
Tan, Q., Daddi, E., Magdis, G., et al. 2014, A&A, 569, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Thomas, D., Maraston, C., Schawinski, K., Sarzi, M., & Silk, J. 2010, MNRAS, 404, 1775 [NASA ADS] [Google Scholar]
Thompson, R., Nagamine, K., Jaacks, J., & Choi, J.-H. 2014, ApJ, 780, 145 [Google Scholar]
Toft, S., van Dokkum, P., Franx, M., et al. 2007, ApJ, 671, 285 [CrossRef] [Google Scholar]
Treu, T., Ellis, R. S., Liao, T. X., et al. 2005, ApJ, 633, 174 [NASA ADS] [CrossRef] [Google Scholar]
Tsai, J. C., & Mathews, W. G. 1995, ApJ, 448, 84 [CrossRef] [Google Scholar]
van de Voort, F., Bahé, Y. M., Bower, R. G., et al. 2017, MNRAS, 466, 3460 [NASA ADS] [CrossRef] [Google Scholar]
van Dokkum, P. G., Whitaker, K. E., Brammer, G., et al. 2010, ApJ, 709, 1018 [Google Scholar]
Veilleux, S., Meléndez, M., Stone, M., et al. 2021, MNRAS, 508, 4902 [NASA ADS] [CrossRef] [Google Scholar]
Ventura, P., Dell’Agli, F., Lugaro, M., et al. 2020, A&A, 641, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Vogelsberger, M., McKinnon, R., O’Neil, S., et al. 2019, MNRAS, 487, 4870 [NASA ADS] [CrossRef] [Google Scholar]
Walters, D., Woo, J., & Ellison, S. L. 2022, MNRAS, 511, 6126 [NASA ADS] [CrossRef] [Google Scholar]
Whitaker, K. E., van Dokkum, P. G., Brammer, G., & Franx, M. 2012, ApJ, 754, L29 [Google Scholar]
Whitaker, K. E., Narayanan, D., Williams, C. C., et al. 2021a, ApJ, 922, L30 [NASA ADS] [CrossRef] [Google Scholar]
Whitaker, K. E., Williams, C. C., Mowla, L., et al. 2021b, Nature, 597, 485 [NASA ADS] [CrossRef] [Google Scholar]
Williams, C. C., Spilker, J. S., Whitaker, K. E., et al. 2021, ApJ, 908, 54 [NASA ADS] [CrossRef] [Google Scholar]
Wiseman, P., Schady, P., Bolmer, J., et al. 2017, A&A, 599, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Woodrum, C., Williams, C. C., Rieke, M., et al. 2022, ApJ, 940, 39 [NASA ADS] [CrossRef] [Google Scholar]
Wright, R. J., Lagos, C. D. P., Davies, L. J. M., et al. 2019, MNRAS, 487, 3740 [NASA ADS] [CrossRef] [Google Scholar]
Wu, P.-F., van der Wel, A., Bezanson, R., et al. 2018, ApJ, 868, 37 [NASA ADS] [CrossRef] [Google Scholar]
Wu, P.-F., Nelson, D., van der Wel, A., et al. 2021, AJ, 162, 201 [NASA ADS] [CrossRef] [Google Scholar]
Wu, P.-F., Bezanson, R., D’Eugenio, F., et al. 2023, ApJ, 955, 75 [NASA ADS] [CrossRef] [Google Scholar]
Yang, T., Davé, R., Cui, W., et al. 2024, MNRAS, 527, 1612 [Google Scholar]
Young, L. M., Scott, N., Serra, P., et al. 2014, MNRAS, 444, 3408 [NASA ADS] [CrossRef] [Google Scholar]
Zheng, Y., Dave, R., Wild, V., & Montero, F. R. 2022, MNRAS, 513, 27 [NASA ADS] [CrossRef] [Google Scholar]
Zhang, J., Li, Y., Leja, J., et al. 2023, ApJ, 952, 6 [NASA ADS] [CrossRef] [Google Scholar]
Zhukovska, S., Dobbs, C., Jenkins, E. B., & Klessen, R. S. 2016, ApJ, 831, 147 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Evolution of the δ_DGR-metallicity plane with stellar mass

Fig. A.1 shows the δ_DGR-metallicity plane at the evolutionary points defined in Sect. 3.2, evolving from SF to completely passive (left to right) and from low to high M_⋆ (top to bottom) for both cluster (red) and field (blue) samples. Galaxies are selected with log(M_⋆/M_⊙) > 9 in order to achieve completeness of the sample (Tacconi et al. 2018, Ghodsi et al. 2024). It is clear that the behaviour of both cluster and field samples is qualitatively similar. As the low mass (log(M_⋆/M_⊙) < 9.75) galaxies evolve, they maintain a quite large δ_DGR. The intermediate mass range (9.75 < log(M_⋆/M_⊙) < 10.5) presents a large scatter of δ_DGR that appears during the quenching of the galaxies (between P2 and P3). Even for super-solar metallicities, galaxies can maintain δ_DGR comparable to the low mass interval. Galaxies with the highest metallicities have systematically lower δ_DGR. For the range of large masses (10.5 < log(M_⋆/M_⊙) < 11.25) metallicities are systematically higher at all evolutionary points. The scatter in P3 is less prominent with respect to the sample of intermediate masses, but there are less galaxies with large δ_DGR. Very few galaxies have measurable metallicity, M_H₂ and M_dust in P4, and they mostly populate the low δ_DGR region of the plane. At this stage, field galaxies more dominantly populate the log(δ_DGR) < − 2 region with respect to cluster galaxies, that are sparsely found with both log(δ_DGR) < − 2 and log(δ_DGR) > − 2.

Fig. A.1.

Evolution of the δ_DGR-metallicity plane according to the four evolutionary points (from left to right) introduced in Sect. 3.2 and the same mass bins used in Fig. 5 (from top to bottom). Galaxies from the cluster sample are represented in red, while galaxies from the field sample are in blue. Vertical grey dotted lines show the solar value for the metallicity (12 + log(O/H)≡8.69). Horizontal grey dashed lines show the literature reference value of δ_DGR ≡ 1/100.

Appendix B: Evolutionary tracks for representative QGs

We show in Fig. B the evolution of the four representative field galaxies (right panel) from Fig. 1 on the dust-stellar-age plane. Three out of the four galaxies (G1, G3, G4) are in the region of high f_dust despite having different evolutionary paths. Such galaxies start the quenching process with an older stellar population than G2. G1 increases its dust content up to the quenching phase and then it decreases it. The t_quench is the longest in the sample. The f_dust at z = 0 is below the resolution for SIMBA, thus is not present in the figure. G2 is the only galaxy ending in the bulk of low dust fraction, after a long and continuous evolution. G2 has the longest t_removal^H₂, one order of magnitude larger than G1, and the second longest t_quench of the four galaxies. The f_dust of G2 at z = 0 is low, but still above the resolution limit for SIMBA. Moreover, G2 is the only galaxy of the sample experiencing the jet-dominated AGN feedback mode during the quenching phase. G3 is the only fast quencher of the four galaxies. G4 has the shortest t_removal^H₂, one order of magnitude lower than G1 and G3 and two orders of magnitude lower than G2, and its f_dust increases at z = 0. It is important to notice that, despite the longer t_removal^H₂, G2 end up with a larger f_dust than G1 and G3 at z = 0, suggesting that the timescale by itself is not enough to predict the evolution of the dust content.

Fig. B.1.

Evolution of galaxies from Fig. 1 in the dust-stellar-age plane for the cluster (left) and field (right) galaxies. The mass-weighted stellar age is on the x-axis, while the f_dust is shown on the y-axis. Each galaxy is colour-coded with the timescale t_removal^H₂ to reach the dry-phase (P4) from the quenched phase (P3). The properties of each galaxy are evaluated from left to right at the P1, P2, P3, and P4. If present, the red circle symbolises the values at z = 0.

The situation is similar for the cluster galaxies (left panel), where the galaxy experiencing jet mode feedback (G6) reaches lower values of f_dust at the dry-phase. G6 has a very short t_removal^H₂, and a relatively low t_quench. The majority of the dust content is lost during the quenching phase, cold gas is lost soon after and the dust content rapidly follows with no measurable dust at z = 0. The galaxy G5 follows the same trend observed for G1, having also a similar SFH.

All Tables

Table 1.

Median roperties of SIMBA field galaxies at quenching.

In the text

All Figures

Fig. 1.

Example sSFR versus look-back time for six representative QGs from this study. The sources’ SFHs are illustrated with blue lines. Panels display QGs along with their corresponding names, quenching times (t_quench), the redshift of the quenching event (z_q), the formation redshift (z_form) and the stellar mass (log(M_⋆/M_⊙). Look-back time is shown on the bottom x-axis in Gyr, while the top x-axis shows the corresponding redshift. The sSFR is illustrated on the y-axis as log-scaled. The SFHs are superimposed with orange and green lines, which mark the thresholds for the start and end of the quenching event, respectively (see the description in Sect. 3.1). The green vertical band delimits the duration of the quenching phase, while red stars are the four critical points we use to determine the sample evolution as described in Sect. 3.2.

In the text

	Fig. 2. PDF distribution of the normalised quenching time for cluster (dark red) and field (blue) galaxies. The quenching time t_quench is divided by the cosmic time τ(z_q) at the redshift z_q at which the galaxy starts quenching. The vertical dotted dark red(blue) line represents the log(t_quench/τ(z_q)) = −1.64(−1.67) separation between fast and slow quenching galaxies in cluster(field).
In the text

	Fig. 3. Fraction of fast (cyan) and slow (red) QGs with redshift for the cluster (continuous lines) and field (dashed lines).
In the text

Fig. 4.

δ_DGR-metallicity plane for cluster QGs (upper panel) and field QGs (lower panel). Positions of QGs within the plane are shown for different evolutionary phases, as indicated in the legend. Contours define positions typical for the peak star formation phase. The intensity of colours refers to the number of sources. The vertical dashed line marks the solar value for the oxygen abundance (Z_⊙), while the horizontal dashed-dotted line shows the reference value for the dust-to-gas ratio in SFGs from the local Universe δ_DGR = 1/100. The positions of galaxies introduced in Fig. 1 are indicated with arrows. The scaling relation between Zgas and δ_DGR for local SFGs (Magdis et al. 2012) is shown as a solid line with equation indicated in the legend. Histograms on the y-axis show the PDF distribution of log(δ_DGR) for both samples, with small ticks marking the 16th and 84th percentile.

In the text

Fig. 5.

Evolution of f_dust with Δ_MS for the field QGs. The Δ_MS is derived following the best-fit relation from Speagle et al. (2014). Each column represents a different phase of the galaxy evolution as indicated in the legend. These are, from left to right: sSFR-peak phase (P1), quenching phase (P2), quenched phase (P3), and gas-removed phase (P4). Each row corresponds to different stellar mass bins calculated at the quenching time. Stellar masses increase from top to bottom, while the sample evolves from left to right. The grey shaded region defines the 0.3 dex scatter of the MS. Points are colour-coded for their δ_DGR calculated at the quenched phase (P3). Points with the same colours in different columns represent the same galaxy across different evolutionary stages. The horizontal arrow marks the direction of evolution with respect to the MS, while the vertical arrow shows the direction expected for increasing dust growth. For each panel, we show the fraction of fast(slow) QGs in the column corresponding to P3. To guide the eye, we show a dashed horizontal line (f_dust = 10⁻⁴), which roughly corresponds to a cold gas fraction of 1% and marks an approximate observational limit of known dusty QGs detected with ALMA (see Section 5.2.4).

In the text

Fig. 6.

Distribution of the f_Edd over the central black hole mass M_BH of the sample of fast (left) and slow (right) field QGs at the quenched stage. The colour coding shows the median dust fraction f_dust computed for the quenched stage (that is, at the P3 evolutionary point). The vertical dashed line represents the lower mass limit for the activation of jet mode in SIMBA (M_BH = 10^7.5). The dashed horizontal line delimits the starting of the transition from a radiative feedback mode to a jet mode (f_Edd = 0.2), while the dashed-dotted line delimits the region with fully set jet mode (f_Edd < 0.02). The red and blue contours represent the 50 and 90 percentile levels for the slow and fast quenching sub-samples respectively.

In the text

Fig. 7.

Evolution of f_dust with stellar age. We track the same evolutionary stages as introduced in Sect. 3.2. Points are colour-coded with the mean quenched timescale, that is the time interval between the P2 and P4. grey dashed and dashed-dotted lines show the exponential decrease of dust abundance assuming a timescale of τ = 2.5 Gyr and τ = 0.7 Gyr respectively and are motivated by Michałowski et al. (2019) and Lee et al. (2024). White squares denote the binned medians from the sample of dusty QGs observed in the COSMOS field (Donevski et al. 2023). Other markers show various QGs followed-up with ALMA: empty star (z = 2, Lee et al. 2024), empty circles (z ∼ 1.5 − 2, Gobat et al. 2022), and empty diamond (z = 2.149, Morishita et al. 2022). Cyan(red) contours define the 50th and 90th percentile of the distribution inferred for the fast(slow) QGs from this study, respectively. In the bottom left panel, the three main peaks of the distribution are labelled A, B, and C.

In the text

Fig. 8.

Qualitative sketch explaining the evolution of f_dust with stellar age within the quenched phase, as presented in the bottom panel of Fig. 7. The cyan cloud illustrates the position of fast quenchers, while the two red clouds illustrate the position of dust-poor and dust-rich slow quenchers. The evolutionary directions (marked with arrows) for each sample are shown along with the main physical mechanisms involved in their dust evolution. The meaning of lines in the large panel is the same as in Fig. 7. The dark-grey shaded regions in the two subplots from the right panel delineate the part of the age-dust plane devoid of the given QG.

In the text

Fig. A.1.

Evolution of the δ_DGR-metallicity plane according to the four evolutionary points (from left to right) introduced in Sect. 3.2 and the same mass bins used in Fig. 5 (from top to bottom). Galaxies from the cluster sample are represented in red, while galaxies from the field sample are in blue. Vertical grey dotted lines show the solar value for the metallicity (12 + log(O/H)≡8.69). Horizontal grey dashed lines show the literature reference value of δ_DGR ≡ 1/100.

In the text

Fig. B.1.

Evolution of galaxies from Fig. 1 in the dust-stellar-age plane for the cluster (left) and field (right) galaxies. The mass-weighted stellar age is on the x-axis, while the f_dust is shown on the y-axis. Each galaxy is colour-coded with the timescale t_removal^H₂ to reach the dry-phase (P4) from the quenched phase (P3). The properties of each galaxy are evaluated from left to right at the P1, P2, P3, and P4. If present, the red circle symbolises the values at z = 0.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.

[1] Akhshik, M., Whitaker, K. E., Leja, J., et al. 2021, ApJ, 907, L8 [NASA ADS] [CrossRef] [Google Scholar]

[2] Akhshik, M., Whitaker, K. E., Leja, J., et al. 2023, ApJ, 943, 179 [NASA ADS] [CrossRef] [Google Scholar]

[3] Akins, H. B., Narayanan, D., Whitaker, K. E., et al. 2022, ApJ, 929, 94 [NASA ADS] [CrossRef] [Google Scholar]

[4] Alarcon, A., Hearin, A. P., Becker, M. R., & Chaves-Montero, J. 2023, MNRAS, 518, 562 [Google Scholar]

[5] Aoyama, K., Kodama, T., Suzuki, T. L., et al. 2022, ApJ, 924, 74 [NASA ADS] [CrossRef] [Google Scholar]

[6] Appleby, S., Davé, R., Kraljic, K., Anglés-Alcázar, D., & Narayanan, D. 2020, MNRAS, 494, 6053 [NASA ADS] [CrossRef] [Google Scholar]

[7] Asano, R. S., Takeuchi, T. T., Hirashita, H., & Inoue, A. K. 2013, Earth Planets Space, 65, 213 [NASA ADS] [CrossRef] [Google Scholar]

[8] Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]

[9] Barišić, I., van der Wel, A., Bezanson, R., et al. 2017, ApJ, 847, 72 [Google Scholar]

[10] Belli, S., Genzel, R., Förster Schreiber, N. M., et al. 2017, ApJ, 841, L6 [NASA ADS] [CrossRef] [Google Scholar]

[11] Belli, S., Newman, A. B., & Ellis, R. S. 2019, ApJ, 874, 17 [Google Scholar]

[12] Belli, S., Contursi, A., Genzel, R., et al. 2021, ApJ, 909, L11 [NASA ADS] [CrossRef] [Google Scholar]

[13] Béthermin, M., Daddi, E., Magdis, G., et al. 2015, A&A, 573, A113 [Google Scholar]

[14] Bianchi, S., & Schneider, R. 2007, MNRAS, 378, 973 [NASA ADS] [CrossRef] [Google Scholar]

[15] Blánquez-Sesé, D., Gómez-Guijarro, C., Magdis, G. E., et al. 2023, A&A, 674, A166 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[16] Boselli, A., Fossati, M., Longobardi, A., et al. 2020, A&A, 634, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[17] Bundy, K., Fukugita, M., Ellis, R. S., et al. 2009, ApJ, 697, 1369 [NASA ADS] [CrossRef] [Google Scholar]

[18] Bundy, K., Bershady, M. A., Law, D. R., et al. 2015, ApJ, 798, 7 [Google Scholar]

[19] Burgarella, D., Nanni, A., Hirashita, H., et al. 2020, A&A, 637, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[20] Caliendo, J. N., Whitaker, K. E., Akhshik, M., et al. 2021, ApJ, 910, L7 [NASA ADS] [CrossRef] [Google Scholar]

[21] Calura, F., Pozzi, F., Cresci, G., et al. 2017, MNRAS, 465, 54 [Google Scholar]

[22] Camps, P., Misselt, K., Bianchi, S., et al. 2015, A&A, 580, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[23] Carnall, A. C., McLeod, D. J., McLure, R. J., et al. 2023, MNRAS, 520, 3974 [NASA ADS] [CrossRef] [Google Scholar]

[24] Casasola, V., Bianchi, S., De Vis, P., et al. 2020, A&A, 633, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[25] Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]

[26] Chauke, P., van der Wel, A., Pacifici, C., et al. 2018, ApJ, 861, 13 [NASA ADS] [CrossRef] [Google Scholar]

[27] Chauke, P., van der Wel, A., Pacifici, C., et al. 2019, ApJ, 877, 48 [NASA ADS] [CrossRef] [Google Scholar]

[28] Chiang, Y.-K., Overzier, R., & Gebhardt, K. 2013, ApJ, 779, 127 [Google Scholar]

[29] Choi, E., Ostriker, J. P., Naab, T., & Johansson, P. H. 2012, ApJ, 754, 125 [NASA ADS] [CrossRef] [Google Scholar]

[30] Clark, C. J. R., Schofield, S. P., Gomez, H. L., & Davies, J. I. 2016, MNRAS, 459, 1646 [NASA ADS] [CrossRef] [Google Scholar]

[31] Corcho-Caballero, P., Ascasibar, Y., Sánchez, S. F., & López-Sánchez, Á. R. 2023, MNRAS, 520, 193 [NASA ADS] [CrossRef] [Google Scholar]

[32] Crenshaw, D. M., Kraemer, S. B., & George, I. M. 2003, ARA&A, 41, 117 [NASA ADS] [CrossRef] [Google Scholar]

[33] Cuppen, H. M., Walsh, C., Lamberts, T., et al. 2017, Space Sci. Rev., 212, 1 [Google Scholar]

[34] Daddi, E., Renzini, A., Pirzkal, N., et al. 2005, ApJ, 626, 680 [NASA ADS] [CrossRef] [Google Scholar]

[35] Dariush, A., Dib, S., Hony, S., et al. 2016, MNRAS, 456, 2221 [NASA ADS] [CrossRef] [Google Scholar]

[36] Davé, R., Thompson, R., & Hopkins, P. F. 2016, MNRAS, 462, 3265 [Google Scholar]

[37] Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827 [Google Scholar]

[38] Davis, T. A., Greene, J., Ma, C.-P., et al. 2016, MNRAS, 455, 214 [Google Scholar]

[39] De Vis, P., Jones, A., Viaene, S., et al. 2019, A&A, 623, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[40] D’Eugenio, C., Daddi, E., Liu, D., & Gobat, R. 2023a, A&A, 678, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[41] D’Eugenio, F., Perez-Gonzalez, P., Maiolino, R., et al. 2023b, Nat. Astron., submitted, [arXiv:2308.06317] [Google Scholar]

[42] Dierckx, P. 1975, J. Comput. Appl. Math., 1, 165 [CrossRef] [Google Scholar]

[43] Donas, J., Deharveng, J.-M., Rich, R. M., et al. 2007, ApJS, 173, 597 [NASA ADS] [CrossRef] [Google Scholar]

[44] Donevski, D., Lapi, A., Małek, K., et al. 2020, A&A, 644, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[45] Donevski, D., Damjanov, I., Nanni, A., et al. 2023, A&A, 678, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[46] Donnari, M., Pillepich, A., Joshi, G. D., et al. 2021, MNRAS, 500, 4004 [Google Scholar]

[47] Dunne, L., Gomez, H. L., da Cunha, E., et al. 2011, MNRAS, 417, 1510 [NASA ADS] [CrossRef] [Google Scholar]

[48] Dwek, E. 1998, ApJ, 501, 643 [NASA ADS] [CrossRef] [Google Scholar]

[49] Faber, S. M., Willmer, C. N. A., Wolf, C., et al. 2007, ApJ, 665, 265 [Google Scholar]

[50] Feldmann, R., & Mayer, L. 2015, MNRAS, 446, 1939 [Google Scholar]

[51] Ferrarotti, A. S., & Gail, H. P. 2006, A&A, 447, 553 [CrossRef] [EDP Sciences] [Google Scholar]

[52] Foltz, R., Wilson, G., Muzzin, A., et al. 2018, ApJ, 866, 136 [NASA ADS] [CrossRef] [Google Scholar]

[53] French, K. D., Smercina, A., Rowlands, K., et al. 2023, ApJ, 942, 25 [NASA ADS] [CrossRef] [Google Scholar]

[54] Fumagalli, M., Labbé, I., Patel, S. G., et al. 2014, ApJ, 796, 35 [Google Scholar]

[55] Galliano, F., Nersesian, A., Bianchi, S., et al. 2021, A&A, 649, A18 [EDP Sciences] [Google Scholar]

[56] Ghodsi, L., Man, A. W. S., Donevski, D., et al. 2024, MNRAS, 528, 4393 [NASA ADS] [CrossRef] [Google Scholar]

[57] Gillman, S., Tiley, A. L., Swinbank, A. M., et al. 2021, MNRAS, 500, 4229 [NASA ADS] [Google Scholar]

[58] Gobat, R., Daddi, E., Magdis, G., et al. 2018, Nat. Astron., 2, 239 [NASA ADS] [CrossRef] [Google Scholar]

[59] Gobat, R., Magdis, G., D’Eugenio, C., & Valentino, F. 2020, A&A, 644, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[60] Gobat, R., D’Eugenio, C., Liu, D., et al. 2022, A&A, 668, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[61] Greene, J. E., Setton, D., Bezanson, R., et al. 2020, ApJ, 899, L9 [NASA ADS] [CrossRef] [Google Scholar]

[62] Gruppioni, C., Béthermin, M., Loiacono, F., et al. 2020, A&A, 643, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[63] Gunn, J. E., Gott, J., & Richard, I. 1972, ApJ, 176, 1 [Google Scholar]

[64] Hamadouche, M. L., Carnall, A. C., McLure, R. J., et al. 2023, MNRAS, 521, 5400 [NASA ADS] [CrossRef] [Google Scholar]

[65] Heckman, T. M., & Best, P. N. 2014, ARA&A, 52, 589 [Google Scholar]

[66] Hirashita, H., & Nozawa, T. 2017, Planet Space Sci., 149, 45 [NASA ADS] [CrossRef] [Google Scholar]

[67] Hirashita, H., Nozawa, T., Villaume, A., & Srinivasan, S. 2015, MNRAS, 454, 1620 [NASA ADS] [CrossRef] [Google Scholar]

[68] Hodge, J. A., & da Cunha, E. 2020, Roy. Soc. Open Sci., 7, 200556 [NASA ADS] [CrossRef] [Google Scholar]

[69] Hopkins, P. F. 2015, MNRAS, 450, 53 [Google Scholar]

[70] Hough, R. T., Rennehan, D., Kobayashi, C., et al. 2023, MNRAS, 525, 1061 [NASA ADS] [CrossRef] [Google Scholar]

[71] Inami, H., Algera, H. S. B., Schouws, S., et al. 2022, MNRAS, 515, 3126 [NASA ADS] [CrossRef] [Google Scholar]

[72] Iwamoto, K., Brachwitz, F., Nomoto, K., et al. 1999, ApJS, 125, 439 [NASA ADS] [CrossRef] [Google Scholar]

[73] Junais, Boissier, S., Boselli, A., et al. 2022, A&A, 667, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[74] Kakimoto, T., Tanaka, M., Onodera, M., et al. 2024, ApJ, 963, 49 [NASA ADS] [CrossRef] [Google Scholar]

[75] Kalita, B. S., Daddi, E., D’Eugenio, C., et al. 2021, ApJ, 917, L17 [NASA ADS] [CrossRef] [Google Scholar]

[76] Kannan, R., Vogelsberger, M., Marinacci, F., et al. 2021, MNRAS, 503, 336 [NASA ADS] [CrossRef] [Google Scholar]

[77] Kaviraj, S. 2014, MNRAS, 437, L41 [NASA ADS] [CrossRef] [Google Scholar]

[78] Kennicutt, R. C., Jr 1998, ApJ, 498, 541 [Google Scholar]

[79] Kocevski, D. D., Onoue, M., Inayoshi, K., et al. 2023, ApJ, 954, L4 [NASA ADS] [CrossRef] [Google Scholar]

[80] Kokorev, V. I., Magdis, G. E., Davidzon, I., et al. 2021, ApJ, 921, 40 [NASA ADS] [CrossRef] [Google Scholar]

[81] Kokusho, T., Kaneda, H., Bureau, M., et al. 2019, A&A, 622, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[82] Kriek, M., van der Wel, A., van Dokkum, P. G., Franx, M., & Illingworth, G. D. 2008, ApJ, 682, 896 [NASA ADS] [CrossRef] [Google Scholar]

[83] Krumholz, M. R., McKee, C. F., & Tumlinson, J. 2009a, ApJ, 693, 216 [NASA ADS] [CrossRef] [Google Scholar]

[84] Krumholz, M. R., McKee, C. F., & Tumlinson, J. 2009b, ApJ, 699, 850 [NASA ADS] [CrossRef] [Google Scholar]

[85] Lammers, C., Iyer, K. G., Ibarra-Medel, H., et al. 2023, ApJ, 953, 26 [NASA ADS] [CrossRef] [Google Scholar]

[86] Lee, M. M., Steidel, C. C., Brammer, G., et al. 2024, MNRAS, 527, 9529 [Google Scholar]

[87] Leroy, A. K., Bolatto, A., Gordon, K., et al. 2011, ApJ, 737, 12 [Google Scholar]

[88] Li, Q., Narayanan, D., & Davé, R. 2019a, MNRAS, 490, 1425 [CrossRef] [Google Scholar]

[89] Li, Z., French, K. D., Zabludoff, A. I., & Ho, L. C. 2019b, ApJ, 879, 131 [NASA ADS] [CrossRef] [Google Scholar]

[90] Lianou, S., Xilouris, E., Madden, S. C., & Barmby, P. 2016, MNRAS, 461, 2856 [NASA ADS] [CrossRef] [Google Scholar]

[91] Lisenfeld, U., Isaak, K. G., & Hills, R. 2000, MNRAS, 312, 433 [CrossRef] [Google Scholar]

[92] Liu, D., Schinnerer, E., Groves, B., et al. 2019, ApJ, 887, 235 [Google Scholar]

[93] Long, A. S., Antwi-Danso, J., Lambrides, E. L., et al. 2023, ApJ, submitted, [arXiv:2305.04662] [Google Scholar]

[94] Looser, T. J., D’Eugenio, F., Maiolino, R., et al. 2024, Nature, 629, 53 [Google Scholar]

[95] Lovell, C. C., Thomas, P. A., & Wilkins, S. M. 2018, MNRAS, 474, 4612 [NASA ADS] [CrossRef] [Google Scholar]

[96] Lovell, C. C., Geach, J. E., Davé, R., Narayanan, D., & Li, Q. 2021, MNRAS, 502, 772 [NASA ADS] [CrossRef] [Google Scholar]

[97] Ma, W., Liu, K., Guo, H., et al. 2022, ApJ, 941, 205 [NASA ADS] [CrossRef] [Google Scholar]

[98] Magdis, G. E., Daddi, E., Elbaz, D., et al. 2011, ApJ, 740, L15 [NASA ADS] [CrossRef] [Google Scholar]

[99] Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012, ApJ, 760, 6 [NASA ADS] [CrossRef] [Google Scholar]

[100] Magdis, G. E., Gobat, R., Valentino, F., et al. 2021, A&A, 647, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[101] Maier, C., Ziegler, B. L., Haines, C. P., & Smith, G. P. 2019, A&A, 621, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[102] Man, A. W. S., Greve, T. R., Toft, S., et al. 2016a, ApJ, 820, 11 [NASA ADS] [CrossRef] [Google Scholar]

[103] Man, A. W. S., Zirm, A. W., & Toft, S. 2016b, ApJ, 830, 89 [Google Scholar]

[104] Man, A. W. S., Zabl, J., Brammer, G. B., et al. 2021, ApJ, 919, 20 [NASA ADS] [CrossRef] [Google Scholar]

[105] Martini, P., Dicken, D., & Storchi-Bergmann, T. 2013, ApJ, 766, 121 [NASA ADS] [CrossRef] [Google Scholar]

[106] Matharu, J., Muzzin, A., Brammer, G. B., et al. 2019, MNRAS, 484, 595 [Google Scholar]

[107] Mattsson, L., Gomez, H. L., Andersen, A. C., et al. 2014, MNRAS, 444, 797 [Google Scholar]

[108] McCormick, A., Veilleux, S., Meléndez, M., et al. 2018, MNRAS, 477, 699 [NASA ADS] [CrossRef] [Google Scholar]

[109] McKinnon, R., Torrey, P., Vogelsberger, M., Hayward, C. C., & Marinacci, F. 2017, MNRAS, 468, 1505 [NASA ADS] [CrossRef] [Google Scholar]

[110] Merritt, D. 1983, ApJ, 264, 24 [NASA ADS] [CrossRef] [Google Scholar]

[111] Michałowski, M. J., Hjorth, J., Gall, C., et al. 2019, A&A, 632, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[112] Michałowski, M. J., Gall, C., Hjorth, J., et al. 2024, ApJ, 964, 129 [CrossRef] [Google Scholar]

[113] Millard, J. S., Eales, S. A., Smith, M. W. L., et al. 2020, MNRAS, 494, 293 [NASA ADS] [CrossRef] [Google Scholar]

[114] Moore, B., Lake, G., & Katz, N. 1998, ApJ, 495, 139 [Google Scholar]

[115] Morishita, T., Abdurro’uf, Hirashita, H., et al. 2022, ApJ, 938, 144 [Google Scholar]

[116] Nanni, A., Bressan, A., Marigo, P., & Girardi, L. 2013, MNRAS, 434, 2390 [NASA ADS] [CrossRef] [Google Scholar]

[117] Nanni, A., Burgarella, D., Theulé, P., Côté, B., & Hirashita, H. 2020, A&A, 641, A168 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[118] Narayanan, D., Turk, M. J., Robitaille, T., et al. 2021, ApJS, 252, 12 [NASA ADS] [CrossRef] [Google Scholar]

[119] Ni, Q., Aird, J., Merloni, A., et al. 2023, MNRAS, 524, 4778 [CrossRef] [Google Scholar]

[120] Nozawa, T., Kozasa, T., & Habe, A. 2006, ApJ, 648, 435 [NASA ADS] [CrossRef] [Google Scholar]

[121] Oppenheimer, B. D., & Davé, R. 2008, MNRAS, 387, 577 [NASA ADS] [CrossRef] [Google Scholar]

[122] Pacifici, C., Kassin, S. A., Weiner, B. J., et al. 2016, ApJ, 832, 79 [Google Scholar]

[123] Pandya, V., Brennan, R., Somerville, R. S., et al. 2017, MNRAS, 472, 2054 [NASA ADS] [CrossRef] [Google Scholar]

[124] Pantoni, L., Lapi, A., Massardi, M., Goswami, S., & Danese, L. 2019, ApJ, 880, 129 [NASA ADS] [CrossRef] [Google Scholar]

[125] Park, M., Belli, S., Conroy, C., et al. 2023, ApJ, 953, 119 [NASA ADS] [CrossRef] [Google Scholar]

[126] Peng, Y.-J., Lilly, S. J., Kovač, K., et al. 2010, ApJ, 721, 193 [Google Scholar]

[127] Planck Collaboration XLIII. 2016, A&A, 596, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[128] Popping, G., & Péroux, C. 2022, MNRAS, 513, 1531 [CrossRef] [Google Scholar]

[129] Popping, G., Somerville, R. S., & Galametz, M. 2017, MNRAS, 471, 3152 [NASA ADS] [CrossRef] [Google Scholar]

[130] Popping, G., Shivaei, I., Sanders, R. L., et al. 2023, A&A, 670, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[131] Priestley, F. D., Chawner, H., Barlow, M. J., et al. 2022, MNRAS, 516, 2314 [NASA ADS] [CrossRef] [Google Scholar]

[132] Relaño, M., De Looze, I., Saintonge, A., et al. 2022, MNRAS, 515, 5306 [CrossRef] [Google Scholar]

[133] Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, A31 [Google Scholar]

[134] Richie, H. M., Schneider, E. E., Abruzzo, M. W., & Torrey, P. 2024, ApJ, 974, 81 [NASA ADS] [CrossRef] [Google Scholar]

[135] Richings, A. J., & Faucher-Giguère, C.-A. 2018, MNRAS, 474, 3673 [Google Scholar]

[136] Robson, D., & Davé, R. 2020, MNRAS, 498, 3061 [NASA ADS] [CrossRef] [Google Scholar]

[137] Rodighiero, G., Daddi, E., Baronchelli, I., et al. 2011, ApJ, 739, L40 [Google Scholar]

[138] Rodighiero, G., Bisigello, L., Iani, E., et al. 2023, MNRAS, 518, L19 [Google Scholar]

[139] Rodríguez Montero, F., Davé, R., Wild, V., Anglés-Alcázar, D., & Narayanan, D. 2019, MNRAS, 490, 2139 [CrossRef] [Google Scholar]

[140] Romano, M., Nanni, A., Donevski, D., et al. 2023, A&A, 677, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[141] Romano, M., Donevski, D., Junais, et al. 2024, A&A, 683, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[142] Rowlands, K., Dunne, L., Maddox, S., et al. 2012, MNRAS, 419, 2545 [NASA ADS] [CrossRef] [Google Scholar]

[143] Rowlands, K., Gomez, H. L., Dunne, L., et al. 2014, MNRAS, 441, 1040 [Google Scholar]

[144] Rowlands, K., Wild, V., Bourne, N., et al. 2018, MNRAS, 473, 1168 [NASA ADS] [CrossRef] [Google Scholar]

[145] Salak, D., Hashimoto, T., Inoue, A. K., et al. 2024, ApJ, 962, 1 [NASA ADS] [CrossRef] [Google Scholar]

[146] Sandstrom, K. M., Leroy, A. K., Walter, F., et al. 2013, ApJ, 777, 5 [Google Scholar]

[147] Sargent, M. T., Daddi, E., Bournaud, F., et al. 2015, ApJ, 806, L20 [NASA ADS] [CrossRef] [Google Scholar]

[148] Schawinski, K., Kaviraj, S., Khochfar, S., et al. 2007, ApJS, 173, 512 [NASA ADS] [CrossRef] [Google Scholar]

[149] Schmidt, M. 1959, ApJ, 129, 243 [NASA ADS] [CrossRef] [Google Scholar]

[150] Schneider, E. E., & Robertson, B. E. 2015, ApJS, 217, 24 [NASA ADS] [CrossRef] [Google Scholar]

[151] Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[152] Setton, D. J., Khullar, G., & Miller, T. B., et al. 2024, ApJ, submitted [arXiv:2402.05664] [Google Scholar]

[153] Shapley, A. E., Cullen, F., Dunlop, J. S., et al. 2020, ApJ, 903, L16 [NASA ADS] [CrossRef] [Google Scholar]

[154] Socolovsky, M., Almaini, O., Hatch, N. A., et al. 2018, MNRAS, 476, 1242 [NASA ADS] [CrossRef] [Google Scholar]

[155] Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15 [Google Scholar]

[156] Spilker, J., Bezanson, R., Barišić, I., et al. 2018, ApJ, 860, 103 [NASA ADS] [CrossRef] [Google Scholar]

[157] Springel, V. 2010, MNRAS, 401, 791 [Google Scholar]

[158] Stacey, H. R., Costa, T., McKean, J. P., et al. 2022, MNRAS, 517, 3377 [NASA ADS] [CrossRef] [Google Scholar]

[159] Sun, F., & Egami, E. 2022, MNRAS, 517, L126 [NASA ADS] [CrossRef] [Google Scholar]

[160] Szpila, J., Davé, R., Rennehan, D., Cui, W., & Hough, R. 2024, MNRAS, submitted, [arXiv:2402.08729] [Google Scholar]

[161] Tacchella, S., Conroy, C., Faber, S. M., et al. 2022, ApJ, 926, 134 [NASA ADS] [CrossRef] [Google Scholar]

[162] Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179 [Google Scholar]

[163] Tan, Q., Daddi, E., Magdis, G., et al. 2014, A&A, 569, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[164] Thomas, D., Maraston, C., Schawinski, K., Sarzi, M., & Silk, J. 2010, MNRAS, 404, 1775 [NASA ADS] [Google Scholar]

[165] Thompson, R., Nagamine, K., Jaacks, J., & Choi, J.-H. 2014, ApJ, 780, 145 [Google Scholar]

[166] Toft, S., van Dokkum, P., Franx, M., et al. 2007, ApJ, 671, 285 [CrossRef] [Google Scholar]

[167] Treu, T., Ellis, R. S., Liao, T. X., et al. 2005, ApJ, 633, 174 [NASA ADS] [CrossRef] [Google Scholar]

[168] Tsai, J. C., & Mathews, W. G. 1995, ApJ, 448, 84 [CrossRef] [Google Scholar]

[169] van de Voort, F., Bahé, Y. M., Bower, R. G., et al. 2017, MNRAS, 466, 3460 [NASA ADS] [CrossRef] [Google Scholar]

[170] van Dokkum, P. G., Whitaker, K. E., Brammer, G., et al. 2010, ApJ, 709, 1018 [Google Scholar]

[171] Veilleux, S., Meléndez, M., Stone, M., et al. 2021, MNRAS, 508, 4902 [NASA ADS] [CrossRef] [Google Scholar]

[172] Ventura, P., Dell’Agli, F., Lugaro, M., et al. 2020, A&A, 641, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[173] Vogelsberger, M., McKinnon, R., O’Neil, S., et al. 2019, MNRAS, 487, 4870 [NASA ADS] [CrossRef] [Google Scholar]

[174] Walters, D., Woo, J., & Ellison, S. L. 2022, MNRAS, 511, 6126 [NASA ADS] [CrossRef] [Google Scholar]

[175] Whitaker, K. E., van Dokkum, P. G., Brammer, G., & Franx, M. 2012, ApJ, 754, L29 [Google Scholar]

[176] Whitaker, K. E., Narayanan, D., Williams, C. C., et al. 2021a, ApJ, 922, L30 [NASA ADS] [CrossRef] [Google Scholar]

[177] Whitaker, K. E., Williams, C. C., Mowla, L., et al. 2021b, Nature, 597, 485 [NASA ADS] [CrossRef] [Google Scholar]

[178] Williams, C. C., Spilker, J. S., Whitaker, K. E., et al. 2021, ApJ, 908, 54 [NASA ADS] [CrossRef] [Google Scholar]

[179] Wiseman, P., Schady, P., Bolmer, J., et al. 2017, A&A, 599, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[180] Woodrum, C., Williams, C. C., Rieke, M., et al. 2022, ApJ, 940, 39 [NASA ADS] [CrossRef] [Google Scholar]

[181] Wright, R. J., Lagos, C. D. P., Davies, L. J. M., et al. 2019, MNRAS, 487, 3740 [NASA ADS] [CrossRef] [Google Scholar]

[182] Wu, P.-F., van der Wel, A., Bezanson, R., et al. 2018, ApJ, 868, 37 [NASA ADS] [CrossRef] [Google Scholar]

[183] Wu, P.-F., Nelson, D., van der Wel, A., et al. 2021, AJ, 162, 201 [NASA ADS] [CrossRef] [Google Scholar]

[184] Wu, P.-F., Bezanson, R., D’Eugenio, F., et al. 2023, ApJ, 955, 75 [NASA ADS] [CrossRef] [Google Scholar]

[185] Yang, T., Davé, R., Cui, W., et al. 2024, MNRAS, 527, 1612 [Google Scholar]

[186] Young, L. M., Scott, N., Serra, P., et al. 2014, MNRAS, 444, 3408 [NASA ADS] [CrossRef] [Google Scholar]

[187] Zheng, Y., Dave, R., Wild, V., & Montero, F. R. 2022, MNRAS, 513, 27 [NASA ADS] [CrossRef] [Google Scholar]

[188] Zhang, J., Li, Y., Leja, J., et al. 2023, ApJ, 952, 6 [NASA ADS] [CrossRef] [Google Scholar]

[189] Zhukovska, S., Dobbs, C., Jenkins, E. B., & Klessen, R. S. 2016, ApJ, 831, 147 [NASA ADS] [CrossRef] [Google Scholar]

Tracing the evolutionary pathways of dust and cold gas in high-z quiescent galaxies with SIMBA

1. Introduction

2. The SIMBA simulation

2.1. Star formation

2.2. Feedback physics

2.3. ISM dust in SIMBA

3. Methods

3.1. Identification of quenched galaxies

3.2. Evolutionary phases and timescales

3.3. Selection of cluster and field galaxies

3.4. Identification of major and minor mergers

4. The global ISM evolution in cluster and field QGs

4.1. Evolution of the quenching mode with redshift

4.2. The δDGR-metallicity plane for quiescent galaxies

4.3. ISM dust-gas abundance within the MS

5. Exploring the mechanisms influencing the dust-gas evolution in high-z QGs

5.1. External channels affecting ISM dust evolution

5.1.1. Cluster environment

5.1.2. Mergers

5.1.3. Rejuvenations

5.2. Internal mechanisms

5.2.1. Impact of AGN feedback

5.2.2. X-ray feedback

5.2.3. Impact of quenching routes and timescales of H2 gas removal on the dust content of QGs

5.2.4. Dissecting the age-evolution of dusty QGs

5.3. Emergence of dusty QGs via prolonged dust growth: a new pathway for the late ISM evolution

5.4. Observational implications for studying high-z dusty QGs

6. Conclusions

Acknowledgments

References

Appendix A: Evolution of the δDGR-metallicity plane with stellar mass

Appendix B: Evolutionary tracks for representative QGs

All Tables

All Figures

Tracing the evolutionary pathways of dust and cold gas in high-z quiescent galaxies with `SIMBA`

4.2. The δ_DGR-metallicity plane for quiescent galaxies

5.2.3. Impact of quenching routes and timescales of H₂ gas removal on the dust content of QGs

Appendix A: Evolution of the δ_DGR-metallicity plane with stellar mass