⌀$\diameter$

11institutetext: National Center for Nuclear Research, Pasteura 7, 02-093
Warsaw, Poland
11email<EMAIL_ADDRESS>22institutetext: SISSA, Via Bonomea
265, 34136 Trieste, Italy 33institutetext: Institute of Cosmology and
Gravitation, University of Portsmouth, Burnaby Road, Portsmouth PO1 3FX, UK
44institutetext: INAF, OAPD, Vicolo dell’Osservatorio, 5, 35122 Padova, Italy
55institutetext: Department of Astronomy, University of Florida, 211 Bryant
Space Sciences Center, Gainesville, FL, USA 66institutetext: Cosmic Dawn
Center (DAWN), Copenhagen, Denmark 77institutetext: Institute for Astronomy,
Royal Observatory, University of Edinburgh, Edinburgh EH9 3HJ, UK
88institutetext: Department of Physics & Astronomy, University of the Western
Cape, Robert Sobukwe Rd, Bellville, 7535, South Africa 99institutetext:
Department of Physics & Astronomy, University of British Columbia, 6224
Agricultural Road, Vancouver, BC V6T 1Z1, Canada 1010institutetext: Department
of Astronomy, University of Massachusetts, Amherst, MA 01003, USA
1111institutetext: INAF - Osservatorio astronomico d’Abruzzo, Via Maggini SNC,
64100, Teramo, Italy 1212institutetext: Department of Astronomy, University of
Washington, Seattle, WA 98195-1700, USA 1313institutetext: Department of
Astronomy, The University of Texas at Austin, 2515 Speedway Boulevard Stop
C1400, Austin, TX 78712, USA 1414institutetext: DTU-Space, Technical
University of Denmark, Elektrovej 327, DK2800 Kgs. Lyngby, Denmark
1515institutetext: Dipartimento di Fisica e Astronomia, Università di Padova,
Vicolo dell’Osservatorio, 3, I-35122 Padova, Italy 1616institutetext: Max-
Planck-Institute für Astrophysik, Karl-Schwarzschild-Strasse 1, 85740,
Garching, Germany

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

G. Lorenzon 11 D. Donevski 1122 K. Lisiecki 11 C. Lovell 33 M. Romano 1144 D.
Narayanan 5566 R. Davé 7788 A. Man 99 K. E. Whitaker 1010 A. Nanni 111111 A.
Long 12121313 M. M. Lee 1414 Junais 11 K. Małek 11 G. Rodighiero 151544 Q. Li
1616

Recent discoveries of copious amounts of dust in quiescent galaxies (QGs) at
high redshifts ($z\gtrsim 1-2$) challenge the conventional view that these
objects have negligible interstellar medium (ISM) in proportion to their
stellar mass. We make use of the SIMBA hydrodynamic cosmological simulation to
explore how dust and cold gas evolve in QGs in relation to the quenching
processes affecting them. We apply a novel method for tracking the changes in
the ISM dust abundance across the evolutionary history of QGs identified at
$0<z\lesssim 2$ in both cluster and field environments. The QGs transition
from a diversity of quenching pathways, both rapid and slow, and exhibit a
wide range of times elapsed between the quenching event and cold gas removal
(from $\sim 650$ Myr to $\sim 8$ Gyr). Contrary to some claims, we find that
quenching modes attributed to the feedback from active galactic nuclei (AGN)
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 as the key mechanism for
dust re-formation, being effective within $\sim 100$ Myr after the quenching
event, and rapidly increasing the dust-to-gas mass ratio in QGs above the
standard values ($\delta_{\rm DGR}\gtrsim 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_{\rm
dust}/M_{\star}\gtrsim 10^{-3}-10^{-4}$) arise from both fast and slow
quenchers, and are prevalent in quenched systems of intermediate and low
stellar masses ($9<\log(M_{\star}/M_{\odot})<10.5$). This strong prediction
poses an immediate quest for observational synergy between e.g., James Webb
Space Telescope (JWST) and the Atacama Large Millimeter Array (ALMA).

###### Key Words.:

Galaxies: ISM, Galaxies: evolution, (ISM:) evolution

##  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 (AGN; 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 & Gott 1972, Boselli et al. 2020, Junais et al. (2022)),
starvation (e.g., van de Voort et al. 2016, Maier et al. 2019), and harassment
(e.g., Merritt 1983, Moore et al. 1998).

Whatever physical mechanism is 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_{2}$ gas is the fundamental
ingredient for the star formation, dust is critical in the thermal balance of
the gas and in the formation of $H_{2}$ 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 (SED) 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), 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
of the dust and cold gas content has been established up to very high
redshifts ($z\sim 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_{\rm H_{2}}$) and gas phase metallicities ($Z_{\rm 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. 2016, 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 ALMA sub-millimeter
follow up of individual QGs in order to attempt to measure $M_{H_{2}}$ 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 $\mu$m ($\left[\rm CII\right]$,
D’Eugenio et al. 2023a).

These studies reported that a non-negligible fraction of massive
($10<\log(M_{\star}/M_{\odot})<11$) QGs at $0.1<z<2.5$ may be fairly dusty,
with average specific dust masses $f_{\rm dust}=M_{\rm dust}/M_{\star}\gtrsim
10^{-4}$, 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_{\rm dust}$ at
a given stellar age, $M_{\star}$ and/or $Z_{\rm 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_{\rm dust}$ and $M_{\rm H_{2}}$, it remains
unclear whether QGs obey the canonical dust-to-gas mass ratio ($\delta_{\rm
DGR}\sim 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 $\rm H_{2}$ gas will persist after the
quenching. Recent studies suggest either a quick depletion of dust and gas
($\lesssim 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 ($\gtrsim 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\sim 3-7$, Carnall et al.
2023, Looser et al. 2023), including few candidate dusty QGs (e.g., Rodighiero
et al. 2022, 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 et al.
2019, Appleby et al. 2020, Donnari et al. 2021). In this work we utilise 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 of the dust and cold gas in QGs
during the different stages of their evolution. To this aim, the present study
plans to give viable answers to key open questions, namely (1) What physical
mechanisms help replenish the dust reservoirs in QGs after the cessation of
star formation?; (2) What is the link between the co-evolution of the dust and
cold gas in QGs and the quenching processes affecting them?; (3) 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 with 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 adopt the same $\Lambda$-CDM cosmology of the SIMBA
simulation, that is $\Omega_{\rm m}=0.3$, $\Omega_{\Lambda}=0.7$, $\Omega_{\rm
b}=0.048$, $H_{0}=$
$68\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{Mpc}}^{-1}$.
We adopt the same cosmology used in Speagle et al. 2014 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 make use of the full-physics (m100n1024) simulation that consists of a box
of side $100\>\rm h^{-1}Mpc$ containing $2\times 1024^{3}$ particles (meaning
$1024^{3}$ dark matter particles and $1024^{3}$ baryon particles). The minimum
stellar mass of a resolved galaxy is $5.8\times 10^{8}\>\rm M_{\odot}$, while
the minimum initial mass for a gas particle is $1.82\times 10^{7}\>\rm
M_{\odot}$ and it is $9.6\times 10^{7}\>\rm M_{\odot}$ for a dark matter
particle. The minimum gravitational softening length is $0.5\>\rm h^{-1}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
$\epsilon_{\star}=0.02$ from Kennicutt Jr (1998). The SFR is proportional to
the gas density and molecular gas fraction, and inversely scale 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 a delayed
feedback from evolved stars and Type $\rm 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, characterized by different regimes of the
Eddington ratio $f_{\rm Edd}=\dot{M}_{\rm BH}/\dot{M}_{\rm Edd}$, where
$\dot{M}_{\rm BH}$ is the total BH accretion rate (Bondi plus Torque-limited)
and $\dot{M}_{\rm Edd}$ is the Eddington accretion rate (accretion rate for
which the BH emit at Eddington luminosity). Values of $f_{\rm Edd}>0.2$
correspond to a radiative mode, where mass-loaded winds are ejected radially
with velocity of the order of $10^{3}\>\rm km\>s^{-1}$. For $f_{\rm 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\text{\times}{10}^{3}\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$. The
transition regime is characterized by winds of increasing velocity as $f_{\rm
Edd}$ decreases to jet mode. A lower mass limit for the triggering of the jet
mode is set to $M_{\rm BH}=10^{7.5}M_{\odot}$ 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 outwards kick for the ISM gas.

###  2.3 ISM dust in SIMBA

In SIMBA the ISM evolves naturally from the initial conditions, accounting for
the production, destruction, and transport of its various components through
feedback mechanisms, thus ensuring self-consistence. 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 gas phase, and then it is turned into
stellar component, providing energy injection by modeling the net effect from
the unresolved contribution of SNe shocks, radiation pressure, and stellar
winds. A spatially uniform photoionizing background is used in the context of
gas heating/cooling. 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 $\rm 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 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 $\rm 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\text{\,}\mathrm{\SIUnitSymbolMicro 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 (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 $\epsilon$ 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 restrict our results to the redshift range $0\lesssim z\lesssim 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_{\rm dust}$
functions in a good agreement with observations at $z<2$ (Li et al. 2019a). It
is also consistent with the observed sub-millimeter number counts up to $z\sim
2$ (Donevski et al. 2020, Lovell et al. 2021) and even with 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\lesssim z\lesssim 2$ we adopt 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; $\rm
sSFR\equiv\rm SFR/M_{\star}$), with SFGs having $\rm sSFR>1/\tau(z)$ and QGs
having $\rm sSFR<0.2/\tau(z)$, with $\tau(z)$ being the age of the Universe at
redshift $z$. The thresholds have empirical origin, coming from the observed
evolution of the normalisation of the MS (i.e., Whitaker et al. 2012,
Fumagalli et al. 2014, Whitaker et al. 2012)111The 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.. 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 examplify six representative star formation histories (SFH) of
SIMBA galaxies that quench at $z\lesssim 2$. The region between the
intersection of the SFH with the two thresholds is highlighted in green. Fig.
1 outlines 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 behavior, the sSFR of the galaxy
can temporary 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
quenchign thresholds is calculate as the intersection of the functions $\rm
sSFR\equiv 1/\tau(z)$ and $\rm sSFR\equiv 0.2/\tau(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\times\tau(z_{q})$, where
$\tau(z_{q})$ is the cosmic time at the redshift $\rm z_{q}$ of the candidate
quenching event.

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 select 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_{\rm q}$) and the crossing of the star forming
phase threshold previous to that (at $z_{\rm sf}$):

$t_{\rm quench}=\tau(z_{\rm q})-\tau(z_{\rm sf}).$ (1)

As we study the redshift range $0<z\lesssim 2$, we exclude from our analysis
all galaxies that do not cross 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 link these two aspects by
comparing different galaxies during a few significant phases of their
evolution. In this way, we can determine the different SFH that galaxies
experience to reach the same stage in their evolution, and monitor which
properties are influencing such paths.

We identify four main evolutionary points that we use in order to monitor the
changes in ISM properties of our sample. These points can be understood as
critical phases in the galaxy 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 identify (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_{\rm gas}=M_{\rm gas}/M_{\star}$) drops below the
arbitrarily threshold of 0.1% that we use to characterise the gas-dry stage.
We adopt 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_{\rm gas}$ measured with ALMA for massive
($M_{\star}>10^{11.5}M_{\odot}$) QGs at $z>1$ (Williams et al. 2021), QGs
identified at $z\sim 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.

Figure 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 ($\rm t_{\rm
quench}$), the redshift of the quenching event ($z_{\rm q}$), the formation
redshift ($z_{\rm form}$) and the stellar mass ($\log(M_{\star}/M_{\odot}$).
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.

###  3.3 Selection of cluster and field galaxies

The analysis is performed using the full-physics SIMBA catalog, and contains
55609 galaxies at $z=0$. We select clusters using the dark halo mass at $z=0$,
requiring a virial mass $M_{\rm halo}>10^{14}M_{\odot}$. 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 using a uniform random sampler on the SIMBA
galaxies at $z=0$ that are outside the halos which we identify as clusters.
The merger tree of each galaxy from the cluster sample has been reconstructed
up to $z\sim 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 redshit may not be part of a virialized 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 to ask if the evolution of QGs
that end up in massive nad dense structures differ from purely field galaxies.
The choice of not tracking directly the environment at each epoch avoid the
ambiguity in defining high redshift large structures. For the field we used
the public information available in the SIMBA catalogs repository.

We select 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\lesssim z\lesssim 2$ are selected for both cluster and field samples. Due
to resolution effects, we limit the stellar mass range of the total SIMBA
sample to $\log(M_{\star}/M_{\odot})>9$. In this way, we also ensure our
sample to be mass-complete.

###  3.4 Identification of major and minor mergers

We trace major and minor mergers along the galaxy history. Mergers are
identified using the relative position of galaxies in the SIMBA catalogs. For
each galaxy (primary) we explore a volume of
$\sim$0.5\text{\,}\mathrm{c}\mathrm{M}\mathrm{pc}$$ in radius around it and
trace the position of all other galaxies (secondary) in that region at each
snapshot. We measure the maximum proper velocity of galaxies in clusters to be
$\sim$2000\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$ on average, and we
approximate the typical timescale between snapshots in the redshift range
$0<z<2$ as $\sim$1\text{\times}{10}^{8}\text{\,}\mathrm{yr}$$, meaning that we
assume the same galaxy to move of proper motion inside a volume of radius
$\sim$200\text{\,}\mathrm{c}\mathrm{k}\mathrm{pc}$$. We trace in this way each
galaxy to its next position in the next snapshot. When a secondary galaxy
enters a region of radius $\sim 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 $\rm R=M_{\rm
secondary}/M_{primary}>1:4$ then the merger is major merger, while if it is
$\rm R<1:4$ the merger is minor (Kaviraj 2013, 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 $\rm 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 use $f_{\rm dust}=M_{\rm dust}/M_{\star}$ and $\delta_{\rm DGR}=M_{\rm
H_{2}}/M_{\star}$, 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.
2016, 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_{\rm 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,
normalized by the cosmic time $\tau(z_{\rm q})$ at the redshift of quenching
$z_{\rm q}$. Both cluster (red) and field (blue) samples follows a clearly
bimodal distribution. We perform a Gaussian fit for each of the two peaks in
cluster (field), inferring $t_{\rm quench,1}=-2.04(-2.16)\pm 0.22(0.22)$ for
the first peak and $t_{\rm quench,2}=-1.16(-1.07)\pm 0.28(0.28)$ for the
second, respectively. In terms of non-normalized timescales, the peak values
correspond to $t_{\rm quench,1}\sim 53.7(40.5)\rm Myr$ and $t_{\rm
quench,2}\sim 407.2(499.9)\rm Myr$ at $z\equiv 1$ and $t_{\rm quench,1}\sim
30.1(22.74)\rm Myr$ and $t_{\rm quench,2}\sim 228.4(280.5)\rm Myr$ at $z\equiv
2$ for cluster(field). We define two different regions of the distribution by
calculating the intersection point of the Gaussian used to fit each peak. The
obtained thresholds are $t_{\rm quench,c,lim}/\tau(z_{\rm q})=-1.64$ and
$t_{\rm quench,f,lim}/\tau(z_{\rm q})=-1.67$ for cluster and field galaxies
respectively. The bimodal distribution of normalised $t_{\rm quench}$ has been
found in SIMBA for a similar range in stellar masses and redshift (Zheng et
al. 2022) and up to $z\sim 4$ (Rodríguez Montero et al. 2019) with a similar
position of the peaks and a separation threshold ($t_{\rm
quench,lim}/\tau(z_{\rm q})=-1.5$), roughly corresponding to $t_{\rm
quench}\sim$175\text{\,}\mathrm{Myr}$$.

Galaxies with associated $t_{\rm quench}/\tau(z_{\rm q})$ lower that the
respective threshold are defined as fast QGs, while larger values identify
what we call slow QGs. This definition consider slow QGs as galaxies that
quench in a timescale comparable with the cosmic age at the time of quenching.
The normalization by $\tau(z_{\rm q})$ ensures that this definition remain
redshift-independent so that a galaxy that quench 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 highlight how SIMBA predicts a peak of much
shorter $t_{\rm quench}$ with respect to other simulations, with a separation
threshold which is also generally smaller. Wright et al. (2019) found
timescales of the order of some Gyr, with a division around $t_{\rm
quench}\sim$1.5\text{\,}\mathrm{Gyr}$$ from the EAGLE simulations. The
characteristic $t_{\rm quench}$ dividing fast and slow QGs as obtained from
the IllustrisTNG simulation is around $\rm t_{\rm
quench}\sim$1\text{\,}\mathrm{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\sim$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 $\rm t_{\rm
quench}=1.0\pm$0.3\text{\,}\mathrm{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
categorization 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.

Figure 2: PDF distribution of the normalized quenching time for cluster (dark
red) and field (blue) galaxies. The quenching time $\rm t_{\rm quench}$ is
divided by the cosmic time $\tau(z_{\rm q})$ at the redshift $z_{\rm q}$ at
which the galaxy starts quenching. The vertical dotted dark red(blue) line
represents the $\log(\rm t_{\rm quench}/\tau(z_{\rm q}))=-1.64(-1.67)$
separation between fast and slow quenching galaxies in cluster(field). Figure
3: Fraction of fast (cyan) and slow (red) QGs with redshift for the cluster
(continuous lines) and field (dashed lines).

To better have a sense of the evolution of $t_{\rm 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\lesssim z\lesssim
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 increase
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\lesssim z\lesssim 2.0$, but for $z\lesssim 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\sim 3$ (Pandya et
al. 2017) showing that at low-$z$ galaxies preferentially evolve through a
slow quenching mode, while the fast quenchign mode increase in importance at
higher redshifts.

###  4.2 The $\delta_{\rm 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.,
$\delta_{\rm DGR}$ and $f_{\rm dust}$) evolve as a function of $Z_{\rm gas}$
(expressed as $12+{\rm 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
define the evolutionary point P4 using an arbitrary threshold for the cold gas
fraction, namely $f_{H_{2}}=\rm M_{H_{2}}/\rm M_{\star}=10^{-4}$. 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+{\rm log(O/H)}=8.69$
(Asplund et al. 2009), while the horizontal dot-dashed line is the literature
reference value ($\delta_{\rm DGR}\sim 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 $\delta_{\rm
DGR}=10.54-0.99\times(12+\log(\rm O/H)$ from Magdis et al. (2012) as a gray
line, and we shade in gray the 0.15 dex scatter of the relation. The fit was
done on a sample of star forming galaxies.

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 $\delta_{\rm DGR}$ with the oxygen abundance, combined
with a large scatter towards super-solar $Z_{\rm gas}$
($10^{-1}\leq\delta_{\rm DGR}\leq 10^{-5}$).

Figure 4: The $\delta_{\rm 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 legend. Contours
define positions typical for the peak star-formation phase. Colorbars refer to
the number of sources. The vertical dashed line marks the solar value for the
oxygen abundance ($Z_{\odot}$), while the horizontal dashed-dotted line shows
the reference value for the dust-to-gas ratio in SFGs from the local Universe
$\delta_{\rm DGR}=1/100$. The positions of galaxies introduced in Fig. 1 are
indicated with arrows. The scaling relation between $Z{\rm gas}$ and
$\delta_{\rm 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(\delta_{\rm DGR})$ for both samples, with small
ticks marking the 16th and 84th percentile.

Surprisingly, soon after quenching, $\sim$50% ($\sim$40%) of cluster (field)
QGs attain $f_{\rm dust}\gtrsim 10^{-3}$ 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 $\sim$30%
(50%) of cluster and field galaxies, respectively. As the galaxies evolve, the
ISM reservoir sustains and has measurable $Z_{\rm gas}$, $M_{\rm dust}$ and
$M_{\rm H_{2}}$, up to the quenched phase (P3) in $\sim 95\%$ of QGs in the
field sample. For cluster QGs, the fraction falls to $\sim 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; $\sim$30% for clusters and 20% for flield QGs), showing no
indication of an evolution with cosmic time. This behaviour resembles the one
observed in DSFGs, for which redshift independent relation of $\delta_{\rm
DGR}$ with $Z_{\rm gas}$ has been reported both at low redshifts ($z\sim 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 populate 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
$\delta_{\rm DGR}$ and $Z_{\rm gas}$ obtained for MS DSFGs (Magdis et al.
2012). Such well-behaved trend is mostly pronounced for galaxies of
intermediate masses (see Appendix A). Intriguingly, a non-negligible number of
QGs has $\delta_{\rm DGR}$ that is compatible or even higher, than observed
for MS DSFGs. As an example, in Figure 4 we also highlight that four out of
six QGs from Fig. 1 maintain (or even increase), their $\delta_{\rm DGR}$
during post-quenching evolution. (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 $\delta_{\rm DGR}$ and $Z_{\rm 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 $\delta_{\rm
DGR}\sim 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 a measurable amount of their ISM in the field and cluster
environment, respectively.

For $\sim 10\%$ of the field QGs, gas fraction remains above our arbitrary
threshold ($f_{\rm gas}>10^{-4}$) 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
negligible trace of ISM is thus around 30-40%. This number is only $\sim 4\%$
for cluster QGs. For both cluster and field QGs the bulk of distribution
remains centered around $\delta_{\rm DGR}\gtrsim 1/100$. These QGs tend to
maintain a large fraction of dust during their last phases of passive
evolution where cold gas content drops significantly; (3) At $12+\log(\rm
O/H)\gtrsim 8.8$ there is a range of $\delta_{\rm DGR}$ that extends over
$\sim 4$ orders of magnitude independently of evolutionary phase or
environment. The absence of a clear trend in Fig. 4 suggests that the balance
between formation and destruction of dust grains is strongly altered at super-
solar $Z_{\rm 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, i.e., ISM-rich mergers, an increase in $f_{\rm 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. 2020). 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_{\rm dust}$ with respect to the offset from the MS ($\Delta_{\rm
MS}$), as defined by Speagle et al. 2014.222For this task, we homogenise the
MS functional form from Speagle et al. (2014) to the IMF by Chabrier (2003).
In Fig. 5 we illustrate this relation along the four key evolutionary stages
introduced in Sect. 3.2.

Figure 5: The evolution of $f_{\rm dust}$ with $\Delta_{\rm MS}$ for the field
QGs. The $\Delta_{\rm 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 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 gray shaded region defines the 0.3 dex scatter
of the MS. Points are color-coded for their $\delta_{\rm DGR}$ calculated at
the quenched phase (P3). Points with the same colors 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 an 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_{\rm
dust}=10^{-4}$), which roughly corresponds to 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_{\rm 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_{\star}$. Namely, for
QGs of lower-masses ($\log(M_{\star}/M_{\odot})<9.75$) the dust content
remains remarkably stable across entire evolutionary history. More complex
behavior emerges for more massive QGs ($\log(M_{\star}/M_{\odot})\gtrsim
9.75$), for which we observe a larger scatter in $f_{\rm dust}$ after the
quenched phase. The QGs of these stellar masses show no clear evolution with
$\Delta_{\rm MS}$. Our data indicate a general trend going from the top to the
bottom row, with higher masses having on average lower $f_{\rm dust}$ at the
gas-removed phase. This implies that $M_{\star}$ has important impact of 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_{\rm 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_{\rm dust}$ and
$\Delta_{\rm MS}$ can be interpreted as an ”age-evolutionary sequence”, such
that QGs from the upper-right side of an $f_{\rm dust}-\Delta_{\rm 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
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, suggest a complexity of processes affecting the ISM
dust abundance. Whitaker et al. (2021a) suggested that a large range in
$f_{\rm dust}$ hints at diverse evolutionary routes to quiescence, a question
we thoroughly explore in the subsequent sections.

Furthermore, Fig. 5 illustrates that $\delta_{\rm DGR}$ and $f_{\rm 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_{2}}$, that is
$f_{H_{2}}\sim 10^{-4}$ (see Sect. 3.2), meaning that the main ISM property
changing between the different galaxies inside a mass bin is $f_{\rm dust}$.
At the start of the quenched phase, the scatter in $f_{\rm dust}$ presents a
smooth gradient in $\delta_{\rm DGR}$, with low $\delta_{\rm DGR}$ being
associated to low $f_{\rm dust}$. This trend is not maintained up to the
final, gas-removed phase. The last column of Fig. 5 showcases QGs with low
$f_{\rm dust}$ and $\delta_{\rm DGR}$ in the quenched phase, attaining higher
values of $f_{\rm dust}$ even with $f_{H_{2}}\sim 10^{-4}$. 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_{\rm dust}$ and $\delta_{\rm DGR}$.
From the column corresponding to the P4 point in Fig. 5 we see that the bulk
of QGs with $f_{\rm dust}>10^{-4}$ at the gas-removed phase is for low mass
galaxies ($9<\log(M_{\rm\star}/M_{\odot})<9.75$), but a non-negligible amount
of galaxies in the range $9.75<\log(M_{\rm\star}/M_{\odot})<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_{\rm dust}$ for the intermediate and large mass ranges is likely a
consequence of the underling mass distribution as shown in the mass division
of the $\delta_{\rm DGR}$-metallicity plane (Fig. A).

##  5 What are 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 in-depth analysis of both external and
internal mechanisms that may be responsible for the overall evolution of the
ISM in QGs. The focus would be set in particular on exploring the origin of
the large scatter in $f_{\rm dust}$. This allows to characterize the
properties of QGs with copious amount 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 regime is reduced stronger for the cluster QGs than for the field ones.
As a consequence, a large number of QGs with super-solar $Z_{\rm gas}$ has
$\delta_{\rm DGR}$ that extends over ¿4 orders of magnitude, including
exotically low values ($\delta_{\rm DGR}\lesssim 10^{-4}$) that are first
reported in SIMBA work by Whitaker et al. (2021a). Such low values are similar
to the low upper limit ($\delta_{\rm DGR}\sim 10^{-4}$) obtained from Planck
observations of galaxy clusters in the local Universe (Adam et al. 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_{\rm gas}$ and density
$n_{\rm e}$ of the gas as $t_{\rm sp}=10^{5}(1+(10^{6}\rm K/T_{\rm
gas})^{3})n_{\rm e}^{-1}\rm yr$ (Hirashita et al. 2015). This means that at
the typical conditions of the hot halos of massive QGs ($T_{\rm
gas}\sim$1\text{\times}{10}^{7}\text{\,}\mathrm{K}$$, $n_{\rm
e}\sim$1\text{\times}{10}^{-3}\text{\,}{\mathrm{cm}}^{3}$$), the hot gas is
expected to support rapid destruction of dust, with a timescale of the order
of $t_{\rm sp}=$1\text{\times}{10}^{8}\text{\,}\mathrm{yr}$$. Since the
sputtering is efficient for $T_{\rm
gas}>$1\text{\times}{10}^{6}\text{\,}\mathrm{K}$$ and decreases rapidly with
$T_{\rm gas}<$1\text{\times}{10}^{6}\text{\,}\mathrm{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_{\rm dust}$.

Vogelsberger et al. (2019) used the AREPO cosmological simulation (Springel
2010) and found a typical timescale of 1-10 Myr for the dust destruction
within the cluster core. We calculate the timescale for SIMBA using the
temperature of the circumgalactic medium (CGM), which in the simulation tops
at around $0.5-1\times 10^{5}$ K. For such temperatures the timescale is
expected to be larger than for the typical scales of T$\sim 10^{6}-10^{7}$ K,
namely $t_{\rm sp}\sim 100\>\rm Myr$. Despite this, a large population
($\sim$30%) of galaxies in SIMBA clusters show dust fractions values well
above the measurements of Adam et al. (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_{\rm sp}>0.5-1$ Gyr support 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_{\rm dust}$ in local elliptical QGs.
These studies concluded that the variation of the dust growth efficiency is
caused by large variations observed in $\delta_{\rm DGR}$, similar to those
seen in our QGs. Whitaker et al. (2021a) conducted several numerical
experiments with higher-resolution SIMBA runs and found that both thermal
sputtering and dust growth are needed to explain the observed $\delta_{\rm
DGR}$ in galaxies with low sSFR. We will 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 (2013), 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 explore the occurrence of merging events among our QGs. From a complete
track of the merger tree of QGs up to $z\sim 2$ we infer that $\sim 20\%$ of
QGs experience at least one major merger in this redshift range and $\sim
15\%$ of QGs experiencing multiple minor mergers. The major merger fraction is
$\sim 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 observational study of Man et al. (2016). They use a large
catalog of mergers observed up to $z\sim 3$, and for major and minor mergers
within the same redshift range as in this study, they infer $\sim 20\%$ and
$\sim 15\%$, respectively. We further explore whether mergers are main drivers
of the scatter in dust-related quantities as observed for the quenched phase
in Fig. 4 and Fig. 5. We infer relatively low fraction of mergers with respect
to the fraction of QGs populating the $\delta_{\rm DGR}$ scatter, which we
measure to be $\sim 45\%$($\sim 20\%$) for QGs in the upper region of the
scatter ($\log(\delta_{\rm DGR})>-3$ and $12+\log(\rm O/H)>8.69$), and $\sim
55\%$($\sim 36\%$) considering the entire scattered region (no constrain on
$\delta_{\rm DGR}$) for P3(P4), respectively. We also check and confirm that
the highest $\delta_{\rm DGR}$ and $f_{\rm dust}$ are not prevalent in
mergers. Consequently, we conclude that observed amount of dust in QGs cannot
be explained by considering the stochastic action of mergers only. The
conclusion applies for 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\leq 1.5$. Same study reported that fast quenchers are
dominating the mass range around 1-3$\times 10^{10}M_{\odot}$, 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 require a galaxy
to be quenched for at least $0.2\times\tau(z_{\rm q})$, with $\tau(z_{\rm q})$
being the cosmic time at quenching. A rejuvenation happens if the sSFR of a
quenched galaxy is sSFR(z) $\geq$ 1/$\tau$(z), that is if the galaxy returns
to the MS. The resulting fraction of rejuvenated galaxies is $\sim 3\%$ for
the clusters and $\sim 2\%$ in field galaxies. Some examples of these objects
are clearly seen in Fig. 5, as they re-renter the MS region during the late
stages of evolution. We further check and find that our rejuvenated galaxies
are not limited to a particular section of $\delta_{\rm DGR}$-metallicity
plane, or $f_{\rm dust}$-MS plane, but populate its entirety. This all
suggests that rejuvenations are not prevalent for producing a large scatter in
$f_{\rm dust}$ and $\delta_{\rm 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\sim 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 have shown that rejuvenation generally occurs during the first
$<$1\text{\,}\mathrm{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 explore 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 verify that
lowering the SF threshold for the selection of rejuvenation to $\rm
sSFR=0.5/\tau(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 constrains on the duration of the rejuvenation event would
result in an even lower fraction.

###  5.2 Internal mechanisms

Analysis from Sect. 5.1 reveal 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) concluded
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_{\text{Edd}}<0.02$ and $M_{\text{BH}}>10^{7.5}M_{\odot}$ Davé et al.
(2019). Figure 6 illustrates the effect of feedback variations (quantified
with the Eddington ratio and $M_{\text{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 ($\sim
90\%$) of fast-quenchers exclusively resides in the jet-dominated region of
the plot. Slow quenchers are more evenly divided, with $\sim$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\times 10^{10}M_{\odot}$, confirming the earlier
results of Rodríguez Montero et al. (2019), Appleby et al. (2020) and Yang et
al. (2024).

A net drop of $\sim 2$ orders of magnitude in $f_{\rm 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_{\rm 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_{\rm 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_{\rm dust})\gtrsim-4$.
Interestingly, fast quenchers have predominantly high $f_{\rm dust}$. Their
median $\log(f_{\rm dust})$ is high in both feedback variants
($-3.09^{0.27}_{0.57}$ in the jet mode, and $-3.06^{0.25}_{0.39}$ in the
radiative mode, respectively). These values are $2-5\times$ higher than those
of slow quenchers ($-4.81^{1.18}_{0.83}$ and $-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.

Figure 6: Distribution of the $f_{\rm Edd}$ over the central black hole mass
$M_{\rm BH}$ of the sample of fast (left) and slow (right) field QGs at the
quenched stage. The color coding shows the median dust fraction $f_{\rm 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 ($\rm M_{\rm BH}=10^{7.5}$). The dashed horizontal line delimits
the starting of the transition from a radiative feedback mode to a jet mode
($f_{\rm Edd}=0.2$), while the dashed-dotted line delimits the region with
fully set jet mode ($f_{\rm Edd}<0.02$). The red and blue contours represent
the 50 and 90 percentile levels for the slow and fast quenching sub-samples
respectively.

This all points towards an important complexity between the timescales for
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 the 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 to
thermal sputtering and collisions. Dust surviving in the CGM was also observed
by Romano et al. (2024), likely carried out of the ISM by outflows. This
provides a scenario in which the dust in 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_{2}$
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\sim
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 found 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 will
depend on the conditions of the ISM, with SNe-heated gas being more efficient
in destroying dust than gas affected by the AGN. After being 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 to
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_{\rm 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_{\rm gas}$, since the gas
surrounding SMBHs 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_{\rm 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, $\sim$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 $\sim$60% for the field QGs, and only $\sim$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. Such QGs are not observed amongst the fast quenchers.

####  5.2.3 Impact of quenching routes and timescales of $\rm H_{2}$ 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_{\rm 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 a further
replenishment. Such evolutionary picture requires a strong anticorrelation of
$f_{\rm 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 explore the overall evolution of $f_{\rm 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:

$t^{H_{2}}_{\rm removal}=\tau(z_{P4})-\tau(z_{P2}),$ (2)

where $\tau(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 dusty QGs in SIMBA start their quenching phase with a wide
distribution in stellar ages and with a relatively narrow range of $f_{\rm
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_{\rm
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_{\rm dust}$, but also a division between the
underlying quenching routes.

Figure 7: Evolution of $f_{\rm dust}$ with stellar age. We track the same
evolutionary stages as introduced in Sect. 4.3. Points are colour-coded with
the mean quenched timescale, that is the time interval between the P2 and P4.
Gray dashed and dashed-dotted lines show the exponential decrease of dust
abundance assuming a timescale of $\tau=$2.5\text{\,}\mathrm{Gyr}$$ and
$\tau=$0.7\text{\,}\mathrm{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\sim 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 labeled A, B, and C.

The lower left panel of Fig. 7 reveals three prominent regions. The first one
is reserved for fast quenchers that occupy 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 with
generally older stellar populations ($\log(\rm Age/yr)\gtrsim 9.4$). Slow
quenchers consist of two distinct clouds: one of slow, relatively dust-rich
QGs (B; $\log(f_{\rm dust})\gtrsim-4$), and one of slow dust-poor QGs (C;
$\log(f_{\rm dust})\lesssim-5$). The positions of these three clouds change
when QGs reach their final phase, characterised with significant deficit of
cold $H_{\rm 2}$ gas. While both group A and B end-up dust-poor and populate
the lower-right side of the diagram, group C maintain an extremely shallow or
almost absent evolution in $f_{\rm dust}$, contrasting the expectations from
age-evolutionary models.

To understand the diverse pathways for dust evolution amongst the three
distinct QG groups mentioned above, we first check their median quenching
times and find them to be $t_{\rm quench,A}\equiv 0.05^{+0.05}_{-0.02}$ Gyr
for fast quenchers, and $t_{\rm quench,B}\equiv 0.82^{+0.98}_{-0.48}$ Gyr and
$t_{\rm quench,C}\equiv 0.61^{+0.70}_{-0.28}$ Gyr, for dust-rich and dust-poor
slow quenchers, respectively. Values are presented in Tab 5.3 together with
other characterizing 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_{\rm quench}$, but accounts for a normalization 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_{\rm quench}$ if they
quenched at low-$z$ or high-$z$ respectively, spreading out the distribution.

If we focus on the median of the $t_{\rm 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_{\rm 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 influenced by a radiative-mode feedback, as shown in Fig. 6.

Akins et al. (2022) coupled SIMBA with POWDERDAY radiative transfer code
(Narayanan et al. 2021), and tracked 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 characterize the ISM content of QGs
that quenched at $z\sim 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_{\rm H_{2}}$ with respect to
other subsamples, $\log(f_{\rm H_{2},B})=-1.37^{+0.17}_{-0.17}$, $\log(f_{\rm
H_{2},A})=-1.37^{+0.27}_{-0.55}$, and $\log(f_{\rm
H_{2},C})=-1.20^{+0.13}_{-0.14}$ respectively.

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. 2014) 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 $\sim 5-6$ Gyr. The cold gas removal timescales obtained for
rapidly quenched dusty QGs in this study are similarly long (see Tab. 5.3).
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
the further dust evolution.

####  5.2.4 Dissecting the age-evolution of dusty QGs

The overall evolution and the large scatter observed in $f_{\rm 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:

$f_{\rm dust}=\frac{M_{\text{dust}}}{M_{\star}}\simeq A\cdot{\rm
exp}\left(-\frac{\rm Age}{\tau}\right),$ (3)

where A is a normalization constant and $\tau$ is the e-folding time. The
typical timescale of dust removal is still uncertain, with some observations
proposing values of $\tau>1$ Gyr (Michałowski et al. 2019, Gobat et al. 2018,
Lee et al. 2024) and others stating much faster timescales $\tau<100-500$ Myr
(Williams et al. 2021, Whitaker et al. 2021b). For comparison, we show in Fig.
7 the evolution of $f_{\rm dust}$ with age adopting a timescale of
$\tau\equiv$0.7\text{\,}\mathrm{Gyr}$$ and one of
$\tau\equiv$2.5\text{\,}\mathrm{Gyr}$$. While part of the distribution of
SIMBA QGs falls into the region spanned by the expected timescales, a cloud of
old QGs (log(${\rm age}_{\star}$/yr) $\gtrsim$ 9.4) maintains relatively large
dust content ($-3<\log f_{\rm dust})<-4$). These values are $3-5\times$ 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_{\rm dust})\lesssim-5$) after the quenching phase, appear
to be much steeper than empirically predicted. These galaxies move almost
vertically on the plane, from the quenching phase to the quenched phase, in
the time $t_{\rm quench}$. 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 $\rm H_{2}$ removal timescale can be short even for
relatively dusty QGs, challenges the common idea that the dust and $\rm H_{2}$
components of the ISM evolve alike. From Tab 5.3 we see that mean values of
$f_{\rm H_{2}}$ 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 $\delta_{DGR}$ and $f_{\rm dust}$ being ¿2 orders of magnitude
lower than in other QGs. This implies that caution should be taken when
inferring $M_{\rm rmgas}$ from $M_{\rm dust}$, as highlighted in Whitaker et
al. (2021a).

Figure 8: A qualitative sketch explaining the evolution of $f_{\rm dust}$ with
stellar age within the quenched phase, as presented in the bottom panel of
Fig. 7. Cyan cloud illustrates the position of fast quenchers, while 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.

Finally, we compare the distribution of simulated galaxies with observations
from Morishita et al. (2022), Lee et al. (2024), Gobat et al. (2022), and
Donevski et al. (2023). The QGs from the first three studies reside at $z\sim
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_{\rm dust}$ predicted by our SIMBA analysis. From the two lower
panels in Fig. 7 we can reveal the possible nature of $f_{\rm dust}$ as well
as 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\sim 2$ by Lee et al.
(2024). The position of this source within the $f_{\rm dust}$-age plane
coincides with group A QGs which suggests the role of AGN jet feedback.
Despite this, the $f_{\rm dust}$ is still considerably high due to a brief and
efficient episode of dust growth, while the cold gas still preserves
detectable amount over $\sim\rm 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 pointed 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_{\rm dust}$ as compared to
the QG from Lee et al. (2024). The QG from Morishita et al. (2022) has a dust
temperature ($\sim 28\rm 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_{\rm dust}$. Interestingly, the two most massive QGs
($M_{\star}>10^{11}M_{\odot}$) have the lowest $f_{\rm dust}$ and come from
the sample of strongly lensed QGs analysed by Gobat et al. (2022). They are on
the 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_{\rm 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^{H_{2}}_{\rm 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 analyses from previous sections reveal 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_{\rm dust}$ remains high even in a late phase
characterised with a cold gas deficit. Following the Eq.2, we infer the median
values for the gas-removal timescales, and find them to be $t^{H_{2}}_{\rm
removal,A}=5.73^{+2.18}_{-3.00}$ Gyr, $t^{H_{2}}_{\rm
removal,B}=2.81^{+4.22}_{-1.66}$ Gyr, $t^{H_{2}}_{\rm
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_{\rm
dust}$ order of $\lesssim 10^{6}M_{\odot}$, 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 for 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 check a dust re-growth
timescale ($\tau_{\rm acc}$) in QGs which we roughly derive following Li et
al. (2019a) and using the same prescription used for SIMBA:

$\displaystyle\tau_{\rm acc}=$ $\displaystyle\,1\times 10^{7}\,\text{yr}$ (4)
$\displaystyle\times\left(\frac{\mathit{a}}{0.1\,$\text{\,}\mathrm{\SIUnitSymbolMicro
m}$}\right)^{-1}\left(\frac{\mathit{n}_{\rm
H}}{100\,$\text{\,}\mathrm{c}\mathrm{m}^{-3}$}\right)^{-1}$
$\displaystyle\times\left(\frac{\mathit{T}_{\rm
gas}}{20\,$\text{\,}\mathrm{K}$}\right)^{-1/2}\left(\frac{\mathit{Z}_{\rm
gas}}{0.02}\right)^{-1}[yr],$

where the first term defines the reference timing of growth activation, $a$ is
a dust grain size with a typical size of
$\sim$0.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ (Li et al. 2019a). The
$n_{\rm H}$ and $T_{\rm gas}$ are the number density and temperature of the
ISM gas respectively, and $Z_{\rm gas}$ is a gas-phase metallicity333Note that
given the SIMBA’s $\sim\rm kpc$ resolution, a multiphase galaxy ISM is not
resolved, which prevent us from knowing the exact dependence of the modelled
dust content to the gas surface density. Molecular gas number density is
calculated from the surface gas density following Popping et al. (2017). We
calculate the gas surface density using the gas mass and the half-mass gas
radius from SIMBA together with the metallicity. We find that dust re-growth
time is the most rapid amongst the slowly quenched dust-rich QGs of group B
($\tau_{\rm acc}=70.64^{+66.98}_{-42.2}$ Gyr). 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.

| Fast | Slow
---|---|---
| $\bm{\log(f_{\rm dust})>-4}$ |
Stellar properties
$\Delta_{\rm MS}$ | $-1.4^{+0.24}_{-0.18}$ | $-1.37^{+0.27}_{-0.25}$ | $-1.44^{+0.35}_{-0.2}$
$\log(M_{{\star}}/M_{\odot})$ | $10.23^{+0.4}_{-0.27}$ | $9.8^{+1.08}_{-0.39}$ | $10.45^{+0.31}_{-0.24}$
$\rm Age[Gyr]$ | $1.69^{+1.55}_{-0.56}$ | $3.92^{+1.56}_{-1.7}$ | $2.24^{+1.63}_{-0.76}$
ISM properties
$\log(f_{{\rm dust}})$ | $-3.09^{+0.27}_{-0.55}$ | $-3.29^{+0.27}_{-0.39}$ | $-5.13^{+0.7}_{-0.62}$
$\log(f_{{\rm H_{2}}})$ | $-1.37^{+0.2}_{-0.3}$ | $-1.37^{+0.17}_{-0.17}$ | $-1.2^{+0.13}_{-0.14}$
$\log(\delta_{{\rm DGR}})$ | $-1.72^{+0.37}_{-0.52}$ | $-1.89^{+0.25}_{-0.45}$ | $-3.95^{+0.81}_{-0.64}$
$12+\log(\rm O/H)$ | $8.79^{+0.12}_{-0.12}$ | $8.73^{+0.26}_{-0.18}$ | $8.96^{+0.07}_{-0.08}$
$\log(\delta_{{\rm DTM}})$ | $-0.27^{+0.15}_{-0.39}$ | $-0.31^{+0.09}_{-0.42}$ | $-2.35^{+0.79}_{-0.64}$
Timescales
$t_{{\rm quench}}[\rm Gyr]$ | $0.05^{+0.05}_{-0.02}$ | $0.82^{+0.98}_{-0.48}$ | $0.61^{+0.7}_{-0.28}$
$t^{{H_{2}}}_{{\rm removal}}[\rm Gyr]$ | $5.73^{+2.18}_{-3.0}$ | $2.81^{+4.22}_{-1.66}$ | $6.42^{+1.8}_{-2.54}$
$t_{{\rm acc}}[Myr]$ | $87.75^{+28.7}_{-28.94}$ | $70.64^{+66.98}_{-42.2}$ | $86.66^{+30.82}_{-27.97}$

Table 1: Medians and associated uncertainties (16th–84th percentile range) of
various physical properties related to fast and slow quenching QGs. The
properties are computed for QGs corresponding to the three regions in Fig. 7
at the quenching phase (see also Fig. 8). Note that for the slow quenchers we
made additional division based on $f_{\rm dust}$ to reflect the difference
between the two clouds as illustrated in Fig. 7.

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 ($\gtrsim 10^{8}\rm 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 an 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\>\rm Myr$ (Hirashita et al.
2015). Since the timescales for $\rm H_{2}$ gas removal in group B are similar
to that for dust destruction, the presence of large $f_{\rm dust}$ in such
galaxies at the dry-phase strongly supports the idea of dust re-growth
happening in $<$500\text{\,}\mathrm{Myr}$$. Ni et al. (2023) investigated the
relation between the observed fraction of X-dominated AGNs and their stellar
age in the COSMOS field. They found a strong anti-correlation, with the oldest
objects having lower AGN accretion rates. They concluded 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
$\sim$1\text{\,}\mathrm{G}\mathrm{y}\mathrm{r}$$. This means that, after the
reduction of $H_{2}$ 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 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 source (see e.g., Nanni
et al. 2013). An alternative possibility, as proposed in several observational
works, suggests a dust mixing via interaction of AGN winds with star
formation-driven outflows. The latter effect seem to be particularly relevant
for lower-mass galaxies (see e.g., Romano et al. 2023). Investigating this
possibility is out of scope of this study.

We recover the same dust-poor but gas-rich QGs (group C) characterised with
exotically low $\delta_{\rm DGR}$ that are first identified by Whitaker et al.
(2021a). The presence of such QGs can be interpret by mean of dilution and
inefficient dust grain growth, as already proposed by Hirashita & Nozawa
(2017)). In this scenario, during a long-lasting ($\sim 5-6$ Gyr) timescales,
radiative-mode AGN feedback steadily removes the $\rm H_{2}$ 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 $\delta_{\rm 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 $\delta_{\rm DGR}$ becomes
highly sensitive to the $M_{\star}$, and strongly affects most massive QGs, as
we show in Section 4.3.

To further explore dust growth scenario we analyse the dust-to-metal ratio
($\delta_{\rm 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 $\delta_{\rm DTM}$ as $\delta_{\rm DTM}$ = $M_{\rm Z}$/$M_{\rm d}$
where $M_{\rm Z}$ = $f_{\rm z}\times{\rm M}_{\rm g}+{\rm M}_{\rm d}$ is the
total mass of metals, $M_{\rm d}$ is the dust mass, $M_{\rm g}$ is the gas
mass and $f_{\rm Z}=27.36\times 10^{12+{\rm log}_{10}({\rm O/H})-12}$ is the
mass fraction of metals. Values are listed in Tab. 5.3. We find that at the
end of the quenching phase, both group A and B have large values of
$\log(\delta_{\rm DTM})>-0.8$, while for group C the value is at least twice
as big. The values for group A and B are comparable with the literature
reference $\log(\delta_{\rm DTM})\sim-0.52$ (i.e. Camps et al. 2015, Clark et
al. 2016) and the values from DustPedia galaxies ($\log(\delta_{\rm
DTM})\sim-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(\delta_{\rm DTM,A})=-2.82^{+2.46}_{-0.81}$, $\log(\delta_{\rm
DTM,B})=-0.27^{+0.1}_{-1.16}$, $\log(\delta_{\rm
DTM,C})=-3.14^{+0.89}_{-0.69}$. While there is a significant decrease in
$\delta_{\rm 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 summarize our results in Fig. 8, where
the properties of the groups A, B, and C are shown in relation to $f_{\rm
dust}$ and stellar age, together with the direction of evolution and dust
growth.

###  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 galaxy becomes fully
quenched, the observed $f_{\rm dust}$ may be resulting from a dust growth,
with re-formation timescales which may be offsetting from those attributed to
$\rm H_{2}$ gas exhaustion. For that reason, we urge caution when using
$f_{\rm dust}$ to convert to $f_{\rm gas}$ and ultimately to $M_{\rm 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_{\rm dust}$ would require longer integration times for
the detection of molecular gas than those simply extrapolated from the typical
$\delta_{\mathrm{DGR}}\sim 1/100$. This is a consequence of enhanced
$\delta_{\mathrm{DGR}}$ and subsequently lower $M_{\rm gas}$. Indeed, recent
studies pointed out to alternative tracers such as [CII] 158 $\mu$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 suggests the existence of dust-rich QGs at $z\sim 1-2$
that were quenched $>5$ Gyr ago ($z_{\rm q}\gtrsim 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_{\star}>10^{9}M_{\odot}$)
of QGs that are quenched over a wide redshift range ($0<z\lesssim 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 $\delta_{\rm DGR}$-$Z_{\rm 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 ($\sim 4$ orders of magnitude) at super-solar $Z_{\rm gas}$. Such
behavior 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/gas evolution in QGs. As they depart from the MS towards the final gas-
dry phase ($f_{\rm gas}<0.1\%$), more than a half of QGs with intermediate
masses ($\log(M_{\star}/M_{\odot})\lesssim 10.5$) sustain high ISM dust
abundance ($f_{\rm dust}>10^{-3}-10^{-4}$; $\delta_{\rm DGR}>1/100$). This is
in contrast to massive QGs ($10.5<\log(M_{\star}/M_{\odot})<11.25$) that
experience a big drop (¿1-2 orders of magnitude) in their $f_{\rm dust}$ and
$\delta_{\rm DGR}$.

  * •

Enhanced $f_{\rm dust}$ and $\delta_{\rm 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_{\rm dust}$ with
stellar age in QGs identified at $0<z\lesssim 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 show little to no evolution in their
$f_{\rm dust}$ with age. This implies long dust depletion timescales (¿2-3
Gyr), a direct consequence of prolonged grain growth.

  * •

Overall, we predict that for significant number of QGs, removal timescales for
the cold $\rm H_{2}$ 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_{\star}/M_{\odot})<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\sim
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).

###### Acknowledgements.

G.L.and D.D acknowledge support from the National Science Center (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’ program. 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

## References

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

##  Appendix A Evolution of the $\delta_{\rm DGR}$-metallicity plane with
stellar mass

Fig. 9 shows the $\delta_{\rm 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_{\star}$ (top to bottom) for both cluster (red)
and field (blue) samples. Galaxies are selected with
$\log(M_{\star}/M_{\odot})>9$ in order to achieve completeness of the sample
(Tacconi et al. 2018, Ghodsi et al. 2023). It is clear that the behavior of
both cluster and field samples is qualitative similar. As the low mass
($\log(M_{\star}/M_{\odot})<9.75$) galaxies evolve, they maintain a quite
large $\delta_{\rm DGR}$. The intermediate mass range
($9.75<\log(M_{\star}/M_{\odot})<10.5$) presents a large scatter of