⌀$\diameter$

11institutetext: National Centre for Nuclear Research, Pasteura 7, 02-093
Warsaw, Poland
11email<EMAIL_ADDRESS>22institutetext: SISSA, Via Bonomea 265,
34136 Trieste, Italy 33institutetext: IFPU - Institute for fundamental physics
of the Universe, Via Beirut 2, 34014 Trieste, Italy 44institutetext:
Department of Astronomy and Physics, Saint Mary’s University, 923 Robie
Street, Halifax, NS B3H 3C3, Canada 55institutetext: Canada Research Chair in
Astronomy and Astrophysics, Tier II 66institutetext: Department of Physics &
Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver,
BC V6T 1Z1, Canada 77institutetext: INAF - Osservatorio di Astrofisica e
Scienza dello Spazio, Via Gobetti 93/3, I-40129, Bologna, Italy
88institutetext: INAF, OAPD, Vicolo dell’Osservatorio, 5, 35122 Padova, Italy
99institutetext: Department of Astronomy, University of Florida, 211 Bryant
Space Sciences Center, Gainesville, FL, USA 1010institutetext: Cosmic Dawn
Center (DAWN), Copenhagen, Denmark 1111institutetext: Institute for Astronomy,
Royal Observatory, University of Edinburgh, Edinburgh EH9 3HJ, UK
1212institutetext: Department of Physics & Astronomy, University of the
Western Cape, Robert Sobukwe Rd, Bellville, 7535, South Africa
1313institutetext: Steward Observatory, University of Arizona, Tucson, AZ
85721, USA 1414institutetext: Astronomy Program, Department of Physics and
Astronomy, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul 08826,
Republic of Korea 1515institutetext: Smithsonian Astrophysical Observatory, 60
Garden Street, Cambridge, MA 02138, USA 1616institutetext: CEA, CNRS,
Université Paris-Saclay, 91191, Gif-sur-Yvette, France 1717institutetext: Max-
Planck-Institute für Astrophysik, Karl-Schwarzschild-Strasse 1, D-85740
Garching, Germany

# In pursuit of giants:

II. Evolution of dusty quiescent galaxies over the last six billion years from
the hCOSMOS survey

D. Donevski 112233 I. Damjanov 4455 A. Nanni 11 A. Man 66 M. Giulietti 2277 M.
Romano 1188 A. Lapi 2233 D. Narayanan 991010 R. Davé 11111212 I. Shivaei 1313
J. Sohn 14141515 Junais 11 L. Pantoni 1616 Q. Li 1717

(Received ; accepted )

The physical mechanisms that link the termination of star formation in
quiescent galaxies and the evolution of their baryonic components, stars and
the interstellar medium (ISM; dust, gas and metals), are poorly constrained
beyond the local universe. In this work, we characterise the evolution of the
dust content in 548 quiescent galaxies observed at $0.1<z<0.6$ as part of the
hCOSMOS spectroscopic redshift survey. This is the largest sample of quiescent
galaxies at intermediate redshifts, for which the dust, stellar and metal
abundances are consistently estimated to date. We analyse how the crucial
markers of a galaxy dust life-cycle, such as specific dust mass ($M_{\rm
dust}$/$M_{\rm\star}$), evolve with different physical parameters, i.e., gas-
phase metallicity ($Z_{\rm gas}$), time since quenching ($t_{\rm quench}$),
stellar mass ($M_{\star}$) and stellar population age.

We find morphology to be important factor of a large scatter ($\sim 2$ orders
of magnitude) in $M_{\rm dust}/M_{\rm\star}$. Quiescent spirals exhibit strong
evolutionary trends of specific dust mass with $M_{\star}$, stellar age and
galaxy size, in contrast to little to no evolution in ellipticals. When
transitioning from solar to super-solar metallicities ($8.7\lesssim\rm
12+\log(O/H)\lesssim 9.1$), quiescent spirals undergo reversal of $M_{\rm
dust}/M_{\rm\star}$, indicative of a change in dust production efficiency.
Through modelling the star formation histories of our objects we unveil a
broad dynamical range of post-quenching timescales ($60\>\rm Myr<t_{\rm
quench}<3.2\>\rm Gyr$). We show that $M_{\rm dust}/M_{\rm\star}$ is the
highest in recently quenched systems ($t_{\rm quench}<500$ Myr), but its
further evolution is non-monotonic as a consequence of different pathways for
dust formation, growth, or removal on various timescales.

Our data are best described by simulations that include dust growth in the
ISM. While this process is prevalent in the majority of galaxies, for $\sim
15\%$ of objects we find evidence of additional dust content acquired
externally, most likely via minor mergers. Altogether, our results strongly
suggest that prolonged dust production on a timescale $0.5-1\>\rm Gyr$ since
quenching may be common in dusty quiescent galaxies at intermediate redshifts,
even if their gas reservoirs are heavily exhausted (i.e., cold gas fraction
$<1-5\%$).

###### Key Words.:

galaxies: evolution – galaxies: formation–galaxies: ISM

##  1 Introduction

The general population of passively evolved galaxies (the so-called quiescent
galaxies, QGs) are known to have no or little ongoing star-formation activity,
which usually places them below the star-forming main-sequence (MS; e.g.,
Daddi et al. 2005, Toft et al. 2007, van Dokkum et al. 2010, Schreiber et al.
2015). Whatever physical mechanism is in place, one of the important outcomes
of quenching star-formation is the evolution of its agents - gas and dust in
the cold interstellar medium (ISM).

The abundance of dust in the ISM is regulated by a complex network of dust
formation and destruction channels (see e.g., Zhukovska et al. 2016 and
references therein). Even though many details of the involved physical and
chemical processes are not fully understood, the study of dust in galaxies is
essential to understand their evolution for a number of reasons. On one hand,
dust is critical in the thermal balance of gas and in shielding the dense
clouds from ultraviolet (UV) radiation, supporting the star-formation (Cuppen
et al. 2017). Dust further affects the spectral energy distributions (SED) of
massive galaxies to the extent that at shorter wavelengths, stellar light is
absorbed by dust and re-emitted in the far-infrared (FIR). On the other hand,
the interplay of dust with metals and gas in the ISM constitutes a vital
component of the baryon cycle. Metals in the gas phase are proposed to play a
role in the growth of dust particles (Asano et al. 2013, Slavin et al. 2015,
Zhukovska et al. 2016, Hirashita & Nozawa 2017, Popping et al. 2017, Pantoni
et al. 2019), and dust itself is an important constituent of the cold ISM,
often considered its reliable tracer not only in star-forming galaxies (e.g.,
Scoville et al. 2016), but also in QGs (Magdis et al. 2021).

Recent advent of IR/sub-mm instruments such as Herschel, NOEMA and ALMA
allowed us to identify QGs with substantial masses of dust and molecular gas
($M_{\rm dust}>10^{7},M_{\odot}$, $M_{\rm gas}>10^{9}M_{\odot}$, Smith et al.
2014, Gobat et al. 2018, Belli et al. 2021). These findings challenged the
historical picture of QGs containing very little ISM material in proportion to
their stellar mass. The vast majority of known studies have focused on
studying the dust and gas in the QGs in the local Universe ($z\sim 0$) either
through galaxy FIR continuum emission (e.g., Smith et al. 2014, Rowlands et
al. 2014, Michałowski et al. 2019), and/or via low transition CO lines (e.g.
Davis et al. 2013, De Vis et al. 2017, Li et al. 2019b, Sansom et al. 2019).
Moreover, several recent works use extensive ALMA and NOEMA follow-ups to
identify small number of dusty QGs in the distant Universe ($z>1-3$, Williams
et al. 2021, Whitaker et al. 2021b). Nevertheless, there is no consensus view
on whether the observed dust in QGs is solely of internal origin (due to past
star-formation), or whether external contributions (i.e., from mergers) are
influential as well. It is also debatable how long the dust and gas in the ISM
will survive after the quenching. Some studies report significant dust and gas
content and prolonged timescales for ISM removal ($\gtrsim 1$ Gyr; Rowlands et
al. 2014, Rudnick et al. 2017, Gobat et al. 2018, Michałowski et al. 2019,
Woodrum et al. 2022), whereas other report rapid depletion of dust and gas
($\lesssim 100-500$ Myr; Sargent et al. 2015, Williams et al. 2021, Whitaker
et al. 2021b).

Quantifying changes in dust abundance as galaxies evolve is of particular
interest when studying dusty QGs. To that end, specific dust masses ($M_{\rm
dust}/M_{\star}$) have been proposed as a useful marker of dust life-cycle
against multiple destruction processes (Calura et al. 2017, Donevski et al.
2020). Michałowski et al. (2019) use the anti-correlation of $M_{\rm
dust}/M_{\star}$ with stellar age to derive the time for dust removal in QGs
and to understand its connection with the shutdown of star-formation. When
$M_{\rm gas}$ is unavailable, the $M_{\rm dust}/M_{\star}$ is often used as an
observational proxy for molecular gas fraction (Rémy-Ruyer et al. 2014, Magdis
et al. 2021). Furthermore, as dust depletes metals from the gas-phase ISM,
adding the information of gas metallicity ($Z_{\rm gas}$) to the relation
between $M_{\rm dust}$ and $M_{\star}$ allows us to better constrain the
dominant dust formation mechanisms (Feldmann 2015, De Vis et al. 2017). Beyond
the local universe, such connection was only studied in the dusty star-forming
galaxies (DSFGs) observed with ALMA at $z\sim 2.3$ (Shivaei et al. 2022).
Despite its importance, the evolution of specific dust masses as a function of
both redshift and physical quantities (in particular, $Z_{\rm gas}$, stellar
age and time since quenching) has not yet been constrained for a statistical
sample of QGs at intermediate redshifts. The task has multiple challenges: (1)
To properly determine dust quantities (i.e., dust luminosity ($L_{\rm IR}$)
and $M_{\rm dust}$) well-sampled IR SEDs towards the Rayleigh-Jeans (RJ) tail
are required; (2) QGs have considerably fainter IR continuum emission than
DSFGs, and careful de-blending of fluxes obtained with single-dish instruments
is required to model the IR SEDs (e.g., Man et al. 2016, Galliano et al.
2021); (3) It has been demonstrated that obtaining direct gas-metallicity
measurements for QGs is a demanding task (Leethochawalit et al. 2019, Kumari
et al. 2021). To partially overcome those issues, Magdis et al. (2021) applied
multi-wavelength stacking analysis on massive
($\log(M_{\star}/M_{\odot})>10.8$) QGs from the COSMOS field up to $z\sim
1.5$. They reported a sharp decline of $M_{\rm dust}/M_{\star}$ from $z\sim 1$
to $z=0$, interpreting it as a consequence of rapid exhaustion of molecular
gas. However, stacking does not allow the investigation of physical markers of
dust evolution in individual sources. This prevents us from answering
important questions, such as: What processes of dust production/removal are
dominant in QGs during their evolution?

From the theoretical standpoint, there is a growing number of studies that
have modelled the dust evolution in different classes of cosmological
simulations (McKinnon et al. 2017, Davé et al. 2019, Hou et al. 2019, Graziani
et al. 2020, Aoyama et al. 2019), semi-analytic models (Lacey et al. 2016,
Popping et al. 2017, Cousin et al. 2019, Vijayan et al. 2019, Lagos et al.
2019, Pantoni et al. 2019) and chemical evolution models (Asano et al. 2013,
Nanni et al. 2020). While some of them show success in reproducing dust
related properties of star-forming galaxies, they still lack observational
constraints from statistical samples of QGs at $z>0$. On top of this, there is
an active debate regarding the dominant dust producers in massive galaxies.
One group of models advocates for major contribution from stellar sources
(i.e., either from asymptotic giant branch (AGB) stars (e.g., Clemens et al.
2010, Valiante et al. 2011) or supernovae (SNe; Gall & Hjorth 2018)), while
others predict prevalence of the dust grain growth in the ISM (e.g., Asano et
al. 2013, Hirashita et al. 2015). The latter channel is thought to be
important at late cosmic epochs ($z<2$), but it is strongly dependent on
”critical” gas-metallicity beyond which it takes over stellar production
(Triani et al. 2020). Thus, the observed co-evolution between the $M_{\rm
dust},M_{\star}$ and $Z_{\rm gas}$ in QGs is crucial to challenge (and
eventually rule out) some of the prescriptions related to the growth and
destruction of dust.

The present paper addresses these challenges. The first major question we aim
to address is how a specific dust mass evolves since quenching in individually
detected QGs. We assembled a statistical data set of spectroscopically
selected dusty QGs at $0.1<z<0.6$ to achieve this goal. We complement deep
multi-wavelength catalogue with the IR fluxes of the carefully de-blended
sources and apply physically motivated SED modelling to self-consistently
derive sources’ physical properties. By considering independent spectroscopic
information such as $Z_{\rm gas}$ and stellar age index ($D_{\rm n}4000$), we
study how the $M_{\rm dust}/M_{\rm\star}$ changes with metallicity, stellar
mass and time since quenching. We then address the second big question,
namely, how the observed evolution of dust properties in QGs can be understood
within the framework of galaxy formation and evolution.

The paper is organised as follows. In Section 2, we describe the data analysed
in this work. In Section 3, we explain the SED fitting methodology and provide
statistical properties for our sample. In Section 4, we present how the
$M_{\rm dust}/M_{\rm\star}$ in the general population of dusty QGs scale with
the galaxy redshift and sSFR. In Section 5, we investigate the morphological
impact on the evolution of the $M_{\rm dust}/M_{\rm\star}$ with physical
parameters such as $Z_{\rm gas}$ and $M_{\star}$. The interpretation of
results with state-of-the-art simulations is presented in Section 6,, while
our main conclusions are outlined in Section 7. Throughout the paper, we
assume a Planck Collaboration et al. (2016) cosmology and the Chabrier initial
mass function (IMF; Chabrier 2003).

##  2 Data and sample selection

The starting point of our sample is the hCOSMOS spectroscopic survey catalogue
(for comprehensive description of the survey, target selection and
observational design, see Damjanov et al. 2018). Briefly, the hCOSMOS is a
dense spectroscopic survey, designed for providing spectra of galaxies from
the COSMOS UltraVISTA catalogue (Muzzin et al. 2013). The hCOSMOS survey is
the magnitude limited survey with no colour selection, and targets UltraVISTA
galaxies that have $r-$band magnitude within the range of $17.8<r<21.3$. The
bright limit of $r=17.8$ is set to match the limiting magnitude of the Sloan
Digital Sky Survey (SDSS) main galaxy sample. The hCOSMOS survey was done with
the multi-fiber fed spectrograph Hectospec, mounted on the 6.5 m MMT
(Fabricant et al. 2005). Hectospec is a 300-fiber optical spectrograph with an
$\sim 1\rm\>deg^{2}$ field of view (FoV) and a fiber diameter of
$1.5^{\prime\prime}$. In total, hCOSMOS survey observed $\sim 1.5\rm\>deg^{2}$
covering the wavelength range between 370 nm and 910 nm at a resolution of
R$\sim 1500$. The hCOSMOS catalogue includes 4362 galaxies with a science
quality spectra. The galaxies are identified across the redshift range of
$0.01<z<0.7$.

We supplement each source from the hCOSMOS sample with rich panchromatic
photometry from publicly available photometric catalogues. Namely, for SED
modelling of the data we adopt homogeneously calibrated multi-wavelength
catalogues released by the Herschel Extragalactic Legacy Project (Małek et al.
2018, Shirley et al. 2019). The catalogue is created with use of the state-of-
the-art de-blending method and cross-matching procedure of individual galaxies
across broad wavelengths. In particular, the de-blending of confusion limited
PACS (100, 160$\mu$m), and SPIRE (250, 350, and 500 $\mu$m) maps is performed
with the probabilistic de-blending code XID+ (Hurley et al. 2017). The code
exploits prior information on the redshift and source position available from
the higher-resolution maps observed at 3.6 $\mu$m, and extracts PACS and SPIRE
fluxes beyond the conventional $\mathit{Herschel}$ confusion limit, achieving
$1\sigma$ cut depth of $\sim 1$ mJy at 250 $\mu$m. The detailed de-blending
procedure are explained in Pearson et al. (2018) and Shirley et al. (2019),
who demonstrate that such approach takes the big advantage of fluctuations
within the confused maps and place strong constraints on the peak of sources’
FIR SEDs even in IR-faint galaxies. All hCOSMOS sources have wealthy
multiwavelength coverage that consists of at least 15 photometric bands: along
with their information in the IR part of the spectra, the optical to mid-IR
data come from Subaru Suprime-Cam (six broadband filters, $V,g,r,i,z$ and
$Y$), VISTA ($J,H,K_{s}$ bands), Spitzer IRAC (four bands) and Spitzer MIPS.

We define the sample of QGs on the basis of a spectral indicator $D_{\rm
n}4000$, which is a ratio of flux in the 4000-4100 $\AA$ and 3850-3950 $\AA$
bands. The $D_{\rm n}4000$ measures the strength of the 4000 $\AA$ break
produced by a large number of absorption lines where ionised metals are the
main contributors to the opacity. Because the strength of the 4000 $\AA$ break
increases with the population age, it is always used as an indicator of a
quiescence. It is worth noting that $D_{\rm n}4000$ is a strong evolutionary
marker because it is insensitive to reddening and does not require a
K-correction, opposite to galaxy colours. From the full hCOSMOS catalogue of
4362 sources, we select QGs as those with $D_{\rm n}4000>1.5$ (1737 sources in
total). This criterion is commonly used for selecting QGs. The choice of using
$D_{\rm n}4000>1.5$ to select QGs is motivated by findings from large
spectroscopic surveys demonstrating that $D_{\rm n}4000$ is strongly bimodal,
and offers a robust division between emission line SF galaxies and QGs at
$D_{\rm n}4000\sim 1.5$ (e.g., Zahid et al. 2016, Haines et al. 2017, Damjanov
et al. 2018, Wu et al. 2018, Utsumi et al. 2020, Wu et al. 2021).
Additionally, it has been shown that selections based on broadband colours
alone are sometimes inefficient in identifying the recently quenched galaxies
(Vergani et al. 2018, Wu et al. 2020). The selection based on $D_{\rm n}4000$
is also proven to be very robust and effective when combined with the
equivalent width $\rm H{\delta}$ absorption line (EW($\rm H\delta$)). On
average, the S/N of individual hCOSMOS spectra is not high enough for
measuring EW($\rm H\delta$) for the full parent mass-complete sample, but we
note that recent deep spectroscopy studies of QGs at $z\sim 0.8$ from the
LEGA-C survey unveil the tight relation between EW($\rm H\delta$) and $D_{\rm
n}4000$ at $D_{\rm n}4000\gtrsim 1.5$ (Wu et al. 2018, Zhang et al. 2022).
Nevertheless, all QG selection criteria produce samples with some levels of
contamination from SF outliers. For the full hCOSMOS catalogue, Damjanov et
al. (2018) compared the identification of QGs and SF galaxies on the basis of
their rest-frame UVJ colours with the identification based on $D_{\rm n}4000$,
and found only $\sim 13\%$ of contaminants (galaxies that have $D_{\rm
n}4000>1.5$, but their UVJ colours are consistent with being SF).
Interestingly, the reported fraction is lower than for any of the colour-based
classifiers tested in Moresco et al. (2013), who applied a variety of methods
to select QGs up to $z\sim 0.5$.

Finally, from all selected QGs, we build the sample of ”dusty QGs” suitable
for our analysis, by imposing source detections in the IR with $\rm S/N\geq 3$
in at least four photometric bands in the mid-IR-to-FIR range ($8\>\mu$m ¡
$\lambda$ ¡$500\>\mu$m).222When available, for the purpose of our SED fitting
we also included upper limits based on $850\mu$m non-detections from the
SCUBA-2 Cosmology Legacy Survey (S2CLS, Geach et al. 2017).. These
requirements are important for achieving the robustness to dust-related
physical parameters estimated from SED fitting, in particular $M_{\rm dust}$
(see e.g., Berta et al. 2016). From the parent sample of 1737 QGs, we find
detectable dust in 618 galaxies ($35\%$ of total). Because our goal is to
explore the interplay between the stars, dust and metals in individually
detected objects, we further chose the sources with reliable measurements of
gas-phase metallicity ($Z_{\rm gas}$, expressed as $12+\log(\rm O/H)$). Gas
phase metallicities of hCOSMOS galaxies are derived by following the method
described in Zahid et al. (2013) and Sohn et al. (2019). The main details of
the procedure are introduced in Appendix A. In this way, we slightly narrowed
the list to 603 dusty QGs for the further SED analysis. We thus leave the
analysis of unselected sources for our ongoing work (Lorenzon et al., in
prep). All 603 galaxies selected for our SED analysis have available
spectroscopic redshifts ($z_{\rm spec}$), distributed over a wide range
($0.1<z<0.6$).

##  3 Panchromatic SED modelling of the data

###  3.1 Tools: CIGALE

We apply full SED (UV+IR) modelling for our QG using the newest release of
Code Investigating GALaxy Emission (CIGALE; Noll et al. 2009, Boquien et al.
2019111https://gitlab.lam.fr/cigale/cigale), a state-of-the-art SED modelling
and fitting code that combines UV-optical stellar SED with an IR component.
For each parameter, CIGALE makes a probability distribution function (PDF)
analysis, providing the output value as the likelihood-weighted mean of the
PDF (with corresponding error as a likelihood-weighted standard deviation). In
the following, we briefly summarise the choice of modules and parameters
presented in Table 1.

#### Stellar component

We construct the stellar component of our SED model adopting Bruzual $\&$
Charlot stellar population synthesis model (Bruzual & Charlot 2003, BC03) with
a Chabrier 2003 IMF. We offer a range of three stellar metallicities closest
to that inferred from the mass-metallicity relation for each object.

Because our primary goal is to infer physical characteristics of QGs, we adopt
the flexible star-formation history (SFH), which is composed of a delayed
component with an additional flexible time at which the star formation is
instantaneously affected. Our choice was motivated by recent studies
demonstrating that the addition of an extra flexibility in the recent SFH is
needed to better recover physical properties (in particular SFR) in a broad
range of QGs (e.g., Ciesla et al. 2016, Hunt et al. 2019, Ciesla et al. 2021,
Suess et al. 2022). The SFH is parametrised as:

$\rm SFR=\begin{cases}\mathit{t}\times e^{-t/\tau_{main}},&\text{when $t\leq
t_{\rm trunc}$}\\\ r\times\rm SFR(\mathit{t}),&\text{when $t>t_{\rm
trunc}$}\end{cases}$ (1)

where $\tau_{main}$ represents the e-folding time of the main stellar
population, while $t_{\rm trunc}$ represents the time at which star formation
is instantaneously affected (quenched). Parameter $r$ is the ratio between
SFRs after quenching and at the moment of quenching, respectively.

Table 1: Parameters used for modelling the SEDs with CIGALE. All ages and
times are given in Gyr.

Parameter | Values |  | Description
---|---|---|---
|  | Star Formation History |
$\tau_{main}$ | 1.0, 1.8, 3.0 |  | e-folding time (main)
Age | 8.5, 9.0, 10.0, 10.5, 11.0, 11.5, 12.0, 12.5, 13.0 |  | Population age (main)
Age of quench | [0.01; 3], 10 values log-sampled |  | Age of the late quenching event
$r$ | [0.001, 1], 10 values log-sampled |  | Ratio of SFR after and prior to quenching
|  | Stellar emission |
IMF | Chabrier 2003 |  | Initial mass function
Z | 0.008, 0.02, 0.05 |  | Stellar metallicity (0.02) in Solar
Separation age | 0.01 |  | Age difference between old and young population
|  | Dust attenuation |
$A^{BC}_{v}$ | 0.3, 0.8, 1.2, 3.3, 3.8 |  | V-band attenuation
slope BC | -0.7 |  | Power law slope of BC attenuation
BC to ISM factor | 0.3, 0.5, 0.8, 1.0 |  | Ratio of the BC-to-ISM attenuation
slope ISM | -0.7 |  | ISM attenuation power law slope
|  | Dust emission |
$q_{PAH}$ | 0.47, 1.12, 3.9, 4.6 |  | Mass fraction of PAH
$U_{\min}$ | 1.0, 2.0, 5.0, 10.0, 15.0, 25.0 |  | Minimum radiation field
$\alpha$ | 2.0 |  | Dust emission power law slope
$\gamma$ | 0.001, 0.01, 0.1 |  | Illuminated fraction
|  | AGN emission |
$r_{ratio}$ | 60. |  | Maximum to minimum radii of the dust torus
$\tau$ | 1.0, 6.0 |  | Optical depth at 9.7$\mu$m
$\beta$ | -0.5 |  | Radial dust distribution within the torus
$\gamma$ | 0.0 |  | Angular dust distribution within the torus
Opening angle | $100^{\circ}$ |  |  $\gamma$ Opening angle of the torus
$\psi$ |  $0^{\circ}$, $90^{\circ}$ |  | Angle between eq.axis and line of sight
$f_{\rm AGN}$ | 0.0, 0.1, 0.25, 0.5, 0.8 |  | AGN fraction

#### Dust attenuation

We adopted Charlot & Fall (2000) double power-law attenuation (CF00), which
assumes that birth clouds (BCs) and the ISM each attenuate light according to
fixed power-law attenuation curves. The formalism is based on age-dependent
differential attenuations between young (age $<10^{7}\>\rm yr$) and old (age
$>10^{7}\>\rm yr$) stars. In this way, the attenuation law intrinsically takes
into account variations in the attenuation curves with time. Both attenuations
are modelled by a power-law function, providing the amount of attenuation in
the $V$ band defined as $\mu=A^{\rm ISM}_{\rm V}/(A^{\rm ISM}_{\rm V}+A^{\rm
BC}_{\rm V})$. Following Battisti et al. (2020) and Ciesla et al. (2021) we
chose to keep both power-law slopes (BC and ISM) of the attenuation fixed at
-0.7, and the parameter $\mu$ fixed to 0.3. Because stars older than
$<10^{7}\>\rm yr$ are only attenuated by $A^{\rm ISM}$, we keep parameter
$A^{\rm ISM}$ varying between 0.3 and 3, which is demonstrated a good choice
for the fitting of quenched objects with detectable dust emission, as recently
pointed out by Ciesla et al. (2021).

#### Dust emission

For modelling the galaxies’ IR SEDs we adopt the physically motivated dust
library of Draine & Li 2007 (DL07 hereafter). In DL07, IR SEDs are calculated
for dust grains heated by starlight for various distributions of intensities.
The majority of the dust is heated by a radiation field with constant
intensity from the diffuse ISM, while muchj smaller fraction of dust
($\gamma$) is exposed to starlight with interstellar radiation field (ISRF)
intensity in a range between $U_{\rm min}$ to $U_{\rm max}$ following a power-
law distribution. We fixed the maximum radiation field intensity to $U_{\rm
max}=10^{6}$ and sample different $U_{\rm min}$, keeping the emission slope
fixed at $\beta=2$. The fraction of the $M_{\rm dust}$ in the form of PAH
grains ($q_{\mathrm{PAH}}$) was sampled within the range
($0.47<q_{\mathrm{PAH}}<4.6$). We log-sample illumination fractions ($\gamma$)
between $0.001$ and $0.1$ with the step of 0.01. As showed in Berta et al.
(2016), sampling towards larger $\gamma$ may increase the chance of
overestimating $M_{\rm dust}$. In our modelling, $L_{\rm IR}$ is an integral
of a SED over the rest-frame wavelength range of $\lambda=8-1000\>\mu$m, while
$M_{\rm dust}$ is derived by fitting and normalising the IR photometry to the
DL07 library.

Figure 1: Distributions of the physical properties estimated for our candidate
quiescent galaxies from the SED fitting with CIGALE. Clockwise, from top left
to bottom right: Goodness of fit expressed as reduced $\chi^{2}$; galaxy
redshift; stellar mass; star-formation rate (SFR), dust mass and linear offset
of the galaxy’s observed SFR to the SFR expected from the modelled MS
($\Delta_{\rm MS}$ in $\log$ scale) as a function of redshift.

We also chose to derive the fractional contribution of the AGN, defined as the
relative impact of the dusty torus of the AGN to the $L_{\rm IR}$ (”AGN
fraction”). We adopt AGN templates presented in Fritz et al. 2006 (see also
Feltre et al. 2012). The parameters in the AGN model were matched to those
from Donevski et al. (2020). Due to computational reasons we somewhat reduce
the number of input options, and model the two extreme values for inclination
angle ($0^{\circ}$ and $90^{\circ}$).

#### Constraints on the parameters

We fit the full datasets with the models defined in previous section. To do
so, we perform the full SED (UV-to-FIR) modelling, but chose the option in
CIGALE to predict $D_{\rm n}4000$. We did several runs of fits until we obtain
a good match between the modelled and the observed $D_{\rm n}4000$. This
procedure is important and ensures that our choice of modelling parameters
with the inclusion of dust correctly matches the galaxy ages and SFH, which is
important to overcome possible degeneracies. Before using our SED-derived
quantities for the science analysis, we confirm that all fitted SEDs are of a
good quality, indicated with a small average $\chi^{2}$ (median is find to be
$\chi^{2}=1.12\pm 0.33$). We also assign the modelling option available within
CIGALE to produce mock catalogues, then following the approach implemented by
Małek et al. (2018) and Donevski et al. (2020) to verify that our SED fitting
procedure does not introduce significant systematics to our dust-related
measurements, in particular $M_{\rm dust}$ (see Appendix B).

###  3.2 Statistical properties of our sample

From the further analysis we exclude AGN powered objects based on their SED
output (55/603 sources with $f_{\rm AGN}>0.25$, or $9\%$ of the total
sample).222We also double-check for additional X-ray-bright AGNs in the COSMOS
(Civano et al. 2016) and find that $\sim 30\%$ of AGNs we identify via SED
modelling also have excess X/radio emission. After this step, the remaining
548 sources are used for our final analysis. We stress here that none of QGs
from our final sample is classified as ”star-forming” based on BPT
diagnostics, while very minor (18/548) number of objects are classified as
”composite”.333Recall that these classifications are based on 1.5” diameter
fiber spectra. Thus, the ”non-classifiable” sample likely includes galaxies
which do not have star formation in the central bulge but may have residual
star formation in the disk.

In Fig. 1, we show the distribution of SED derived properties. We infer the
median redshift of $z=0.37$ with the corresponding 16-84th percentile range
($z=0.26-0.48$). As expected for evolved objects, QGs from our hCOSMOS sample
are very massive systems ($\log(\rm
M_{\star}/M_{\odot})=10.92^{+0.24}_{-0.32}$), The median spectral age
indicator is $D_{\rm n}4000=1.73$, while the median and interquartile range of
IR luminosity is $\log(L_{\rm IR}/L_{\odot})=10.28^{+0.30}_{-0.42}$, which is
two orders of magnitude lower than in DSFGs at a given redshift and stellar
mass range. We apply the functional form of the MS defined by Speagle et al.
2014 (their ’best-fit’, provided with Eq. 28) to check the position of our QGs
with respect to the MS. With this check we unveil that $88\%$ of dusty QGs
from our sample reside well bellow the MS, with the linear offset to the SFR
expected for MS objects being $\log(\Delta_{\rm MS})\leq-0.6$. Remaining $\sim
12\%$ of sources can be considered galaxies approaching ”green valley” or MS.

We also check if the results presented in the rest of the text hold if
alternative selection criteria for QGs are applied. In Appendix C we confront
different tracers for galaxy quiescence (both on the basis of rest-frame
colours and spectroscopy) and demonstrate that selection have little or no
impact on findings presented in the next sections.

Figure 2: Observed redshift evolution of $M_{\rm dust}/M_{\star}$ in dusty QGs
from this work. Individual values are displayed with circles, coloured with
corresponding $D_{\rm n}4000$. Binned medians and associated uncertainties
(16th-84th percentile range) are shown with sandy brown circles and lines,
respectively. For comparison, we show the best fit from the stacking analysis
of Magdis et al. (2021) (red dashed line). The dark-grey line and shaded area
describe the modelled evolution of ALMA detected MS galaxies with a functional
form of $M_{\rm dust}/M_{\star}\propto k\times(1+z)^{2.5}$ (Donevski et al.
2020), while the dark-grey dashed line shows the shift of 1 dex of this MS
scaling relation.

##  4 ISM diagnostics of the general population of QGs from hCOSMOS

###  4.1 Evolution of $M_{\rm dust}/M_{\star}$ with redshift

In order to gain insight into evolution of galaxy dust content in our QG, we
apply $M_{\rm dust}/M_{\star}$ as a tool to assess the efficiency of the
specific dust production and destruction mechanisms. From our multi-band SED
fitting, we find the median of $M_{\rm dust}/M_{\star}=1.5^{+0.6}_{-0.7}\times
10^{-4}$. This is more than an order of magnitude lower than what has been
estimated for MS DSFGs within the same redshift and stellar mass range (e.g.,
da Cunha et al. 2010, Driver et al. 2018), and is two orders of magnitude
lower than in dusty starbursts (da Cunha et al. 2015, Donevski et al. 2020,
Pantoni et al. 2021). As seen in Fig.2,, where we show the redshift evolution
of the $M_{\rm dust}/M_{\star}$ for the full sample, some of our dusty QGs are
approaching the main-sequence (MS) relation modelled for ALMA identified dusty
galaxies at $z<1$ (see the Section 4.4 in Donevski et al. 2020 for details).
That said, a handful of QGs in our sample have $M_{\rm dust}/M_{\star}$
comparable to those of MS DSFGs albeit having lower $L_{\rm IR}$ on average.
The median values of our QGs agrees reasonably well with the modelled cosmic
evolution of $M_{\rm dust}/M_{\star}$ in QGs at $0<z<1$, obtained from
stacking IR-to-radio maps of colour-selected, massive
($M_{\star}>10^{10}M_{\odot}$) galaxies in the COSMOS field (Magdis et al.
2021). The average $M_{\rm dust}/M_{\star}$ from Magdis et al. (2021) are
0.2-0.5 dex bellow the medians from this work, which is expected due to the
stacking technique they used to reach fainter $L_{\rm dust}$. Nevertheless,
our data agree within $1\sigma$, yet suggesting slightly milder decline with
$z$. As radio-mode AGN feedback is expected to dominate the faster decline of
$M_{\rm gas}$ (and supposedly $M_{\rm dust}$; Feldmann 2015, Wu et al. 2018),
the exclusion of candidate AGNs from our final sample likely contribute to the
shallower evolutionary trend with respect to those inferred by Magdis et al.
(2021). In the accompanying paper (Lorenzon et al. in prep.), we stacked the
FIR maps of all unselected QGs from the hCOSMOS parent sample, and obtain the
upper limits of $M_{\rm dust}/M_{\star}=1.4\times 10^{-5}$ and $M_{\rm
dust}/M_{\star}=3\times 10^{-5}$, for $z=0.25$ and $z=0.37$ respectively.
These values, as well as their redshift evolution, are very similar to that of
Magdis et al. (2021), implying a less dramatic cosmic evolution of the ISM in
dusty QGs in comparison to the general population of QGs.

In Fig. 2 we see that galaxies with higher $D_{\rm n}4000$ generally attain
lower specific dust mass than those with smaller $D_{\rm n}4000$. Because
$D_{\rm n}4000$ is a proxy for stellar population age, more evolved systems
with older stellar populations are naturally expected to have lower $M_{\rm
dust}/M_{\star}$ (e.g., Hjorth et al. 2014, Pantoni et al. 2019, Michałowski
et al. 2019). However, Fig. 2 reveals that even for objects observed at the
same redshift, there is a significant spread in $M_{\rm dust}/M_{\star}$ which
does not always coincide with $D_{\rm n}4000$. This could suggest diverse ISM
conditions in our QGs, (i.e., large variability of molecular gas fractions,
recent merger activities etc.) which may result in different dust
production/destruction timescales. This slightly challenges the interpretation
by Magdis et al. (2021), who argue that prompt decline of $M_{\rm
dust}/M_{\star}$ must be due to quick gas removal with no further
replenishment (i.e., due to minor mergers or accretion from the cosmic web).

###  4.2 Evolution of $M_{\rm dust}/M_{\star}$ with sSFR

Figure 3: Evolution with sSFR for $M_{\rm dust}/M_{\star}$ (left) and $M_{\rm
dust}/\rm SFR$ (right). Left panel:The binned medians and corresponding
16th–84th percentiles of this sample are presented with black circles and with
dark grey shaded area, respectively. Open crosses indicate the subsample of
quiescent objects drawn from the Dustpedia archive (Galliano et al. 2021).
Filled crosses show the compilation of dusty QGs and pSB galaxies at $z<0.1$
detected in H-ATLAS (Dariush et al. 2016). Displayed with salmon-coloured
symbols are estimates for the stacked and individual QGs at higher-redshifts:
stacks are represented with diamonds and pentagons (Gobat et al. 2018, Magdis
et al. 2021), while the individual detections are annotated with triangles
(Whitaker et al. 2021a,). In both panels the cyan-shaded region illustrates
the parameter space where majority of dusty pSB galaxies detected with ALMA
reside (Li et al. 2019b); Right panel: The symbols are of the same meaning as
in the left panel. Dashed lines roughly track the position of sources within
different $T_{\rm dust}$, calculated by scaling the strength of the ISRF
intensity as derived from CIGALE.

In star-forming galaxies, dust mass is expected to be a good tracer of the
$M_{\rm gas}$ (da Cunha et al. 2010, Santini et al. 2014, Rowlands et al.
2014, Scoville et al. 2017, Kirkpatrick et al. 2017, Aoyama et al. 2019),
while $M_{\rm gas}$ and SFR are linked through the known Kennicutt–Schmidt
relation (KS, Schmidt 1959, Kennicutt 1998, Sargent et al. 2014). It is thus
important to unveil how $M_{\rm dust}$, $M_{\star}$ and SFR are related in
dusty QGs at intermediate redshifts.

In Fig. 3, we show two ISM diagnostics as a function of sSFR: $M_{\rm
dust}/M_{\star}$ (left panel), and $M_{\rm dust}/\rm SFR$ (right panel). Both
relations extend over two orders of magnitude, with standard deviations of
$\sim 0.4-0.5\>\rm dex$ in each sSFR bin. This hints to complexity of ISM
conditions in dusty QGs, similar to those seen in high-$z$ DSFGs (e.g., da
Cunha et al. 2015, Calura et al. 2017, Strandet et al. 2017, Dudzevičiūtė et
al. 2020, Donevski et al. 2020). Such significant spread in $M_{\rm
dust}/M_{\star}$ and $M_{\rm dust}/\rm SFR$ is interesting to understand, as
QGs are expected to share similar, redshift independent SEDs, consequence of
almost constant dust temperature ($T_{\rm dust}\sim 20$ K) and less turbulent
ISM as compared to DSFGs (Andreani et al. 2018, Magdis et al. 2021). If we
adopt our SED derived $\langle U\rangle_{\rm min}$ and convert it to $T_{\rm
dust}$ via calibration given by Schreiber et al. (2018), we obtain values
ranging from 17K to 33K, with the median and absolute median deviation of
$24\pm 2.3$ K. This $T_{\rm dust}$ is $\sim 4$ K warmer than the average value
found by Magdis et al. (2021). Variations in $T_{\rm dust}$ are seen in each
sSFR bin, in agreement with Martis et al. (2019) and Nersesian et al. (2019)
who reveal the similarly wide range of values among the dusty galaxies with
low sSFR. Differences in $T_{\rm dust}$ can be due to the dust grains being
exposed to different ISRFs (Nersesian et al. 2019), but can also be related to
intrinsic properties of dust, such that galaxies with the higher fraction of
larger (smaller) grains produce cooler (warmer) SEDs (Relaño et al. 2022,
Nishida et al. 2022). This, again, refers to diverse channels for the ISM
evolution in observed QGs.

The relations shown in Fig. 3 can provide insight into the evolutionary stages
of our dusty QGs. For example, the trend of $M_{\rm dust}/M_{\star}$ with $\rm
sSFR$ can be partially interpreted as age-evolutionary sequence (Calura et al.
2017), such that QGs from the upper-right side of the diagram could be
dominated by objects that have younger stellar age and higher gas fraction
($f_{\rm gas}=M_{\rm gas}/(M_{\star}+M_{\rm gas})$). The subsequent decrease
of sSFR is than mostly due to exhaustion of their gas reservoirs, reflecting
the efficiency of ISM removal (Gobat et al. 2020). The anti-correlation of
$M_{\rm dust}/\rm SFR$ and sSFR is often viewed either as a proxy for the
metallicity-dependent gas depletion timescale ($M_{\rm dust}/\rm
SFR\propto\tau_{\rm dep}(Z/Z_{\odot}$), e.g., Béthermin et al. 2015, Magdis et
al. 2021), or inverse of the mean radiation field ($\left\langle U_{\rm
mean}\right\rangle$, e.g., Hunt et al. 2014, Martis et al. 2019). The former
implies that galaxies with smaller $M_{\rm dust}/\rm SFR$ and high sSFR can be
compatible with shorter gas depletion timescales (i.e., $\tau_{\rm dep}\sim
500$ Myr -1 Gyr), indicative of MS DSFGs (Béthermin et al. 2015, Scoville et
al. 2016, Liu et al. 2019). Consequently, objects that populate upper left
part of $M_{\rm dust}/\rm SFR$-sSFR plane are expected to have longer
depletion, order of few Gyr.

Interestingly, the literature lacks systematically selected statistical
datasets of dusty QGs that would be fully suitable for comparison with our
sample. Therefore, we chose to compare our data to the widely used benchmark
samples of dusty galaxies in the nearby universe, such as the H-ATLAS sample
of optically red, dusty galaxies (Dariush et al. 2016, hereafter D16), and the
Dustpedia archive (Galliano et al. 2021, hereafter G21; see also Lianou et al.
2019, Nersesian et al. 2019). The D16 sample comprises of 78 galaxies at
$z<0.1$ that are selected based on $>5\sigma$ Herschel detection at
$250\>\mu$m, and for which SDSS counterparts have red colours based on $NUV-r$
diagnostics. The sample from G21 consists of $\sim 770$ very local galaxies
with angular sizes of $D_{25}>1^{\prime}$ and similarly significant Herschel
detections. The majority of Dustpedia sources are SF galaxies, and for
comparison we only chose those sources (180 in total) within the same range of
sSFR and $M_{\star}$ as our QGs ($\log(\rm sSFR/\rm yr^{-1})<-10$ and
$M_{\star}>10^{10}M_{\odot}$). However, this only allows very rough comparison
because Dustpedia galaxies are not specifically selected as quiescent, i.e.
based on their $D_{\rm n}4000$. On top of this, the median stellar mass of G21
sample is $\log(M_{\star}/M_{\odot})=10.4$, which is $\sim 0.5$ dex below the
median of our QGs. Therefore, we stress that any conclusion resulting from
comparisons with these two samples should be taken with caution. Nevertheless,
as can be derived from both panels of Fig. 3, our estimates follow the similar
evolutionary shape with sSFR as galaxies from G21 and D16, while the derived
dust related properties differ among the samples. The points from G21 depart
from our median $M_{\rm dust}/M_{\star}$ and $M_{\rm dust}/\rm SFR$ at higher
sSFR, while those from D16 lie systematically above our QGs over the entire
sSFR range. Although all three works applied different fitting methods and
model assumptions444In G21 the CIGALE was used for obtaining $M_{\star}$ while
$M_{\rm dust}$ is inferred by the code HerBIE (Galliano et al. 2018) which
adopts the framework of THEMIS dust model (Jones et al. 2017). The properties
of D16 are derived by MAGPHYS (da Cunha et al. 2008), it is more likely that
diversity in values is due to physics and selection bias. The $T_{\rm dust}$
derived in G21 has a median of $\sim 19\rm K$, which is $\sim 5$K colder than
average $T_{\rm dust}$ of our sample. Therefore, the difference with G21 seen
towards higher sSFR likely arises due to G21 galaxies being less massive and
more colder than our QGs. The average $M_{\star}$ of D16 galaxies, on the
other hand, is similar to ours ($\log(M_{\star}/M_{\odot})=10.8$), but their
$M_{\rm dust}$ are systematically higher. This can be a consequence of
selection bias and noise effects. The H-ATLAS survey covers an area that is
$\sim 100$ times larger but $\sim 5$ times shallower than the COSMOS field.
Hence, $>5\sigma$ selection would return the most extreme objects that are
much brighter and ”dustier” than our QGs. In addition, noise effects can
impact $M_{\rm dust}$ estimation (i.e., Berta et al. 2016). As a consequence,
most objects from D16 enter the region of ALMA-observed post-SB galaxies (Li
et al. 2019b, marked as the cyan shaded area in Fig. 3).

The median trend of our specific dust masses agree within $1\sigma$ with those
derived for QGs at higher-$z$ ($0.5<z<2$; Gobat et al. 2018, Magdis et al.
2021, Whitaker et al. 2021a). These studies reach somewhat contrasting
conclusions regarding the physical properties of detected objects. Magdis et
al. (2021) argue of long $\tau_{\rm dep}$ (¿Gyr) and roughly constant $T_{\rm
dust}\sim 20$ K, whereas the two QGs from Whitaker et al. (2021a) are
characterised by shorter depletion timescales ($\tau_{\rm dep}\sim 100$ Myr),
and very different $T_{\rm dust}$ (14 K and 34 K, respectively), and $f_{\rm
gas}$ ($4.6\pm 0.5\%$ and $0.6\pm 0.1\%$, respectively). The dusty QG
identified by Gobat et al. (2018) is found to host a larger cold gas fraction
($5-10\%$) under dust temperature of $T_{\rm dust}=20$ K. Whitaker et al.
(2021a) suggested that a large range in specific dust masses hints to diverse
evolutionary routes to quiescence. Therefore, it is of vital importance to
investigate what causes the large dispersion in $M_{\rm dust}/M_{\star}$ seen
in QGs. We thoroughly explore this question in the next Sections.

Figure 4: Structural and SED-derived physical properties of QGs from hCOSMOS.
Normalised distributions of values for various physical quantities are
displayed for two morphological categories; quiescent spirals and ellipticals
(sQGs and eQGs; yellow and dark cyan colours, respectively). The first row
delineates properties obtained from the Hectospec spectroscopy (redshift,
$D_{\rm n}4000$ and gas-metallicity (as oxygen abundance) expressed as $\rm
12+\log(\rm O/H)$). The second row displays some of the main structural
properties derived from HST images observed via F814W filter (Gini index,
concentration and effective radius in kpc) from the catalogue of Cassata et
al. 2007). The third and fourth row show various physical properties estimated
from multiwavelength SED modelling of our QGs with the code CIGALE. These are,
from upper left to lower right, respectively: $\log(M_{\star}/M_{\odot})$,
$\log(\rm Age/yr)$, $\log(\rm t_{\rm quench}/yr)$, $\log(M_{\rm
dust}/M_{\odot})$, $T_{\rm dust}$ and $\log(L_{\rm dust}/L_{\odot})$. Figure
5: Upper panel: Evolution of $M_{\rm dust}/M_{\star}$ as a function of gas-
phase metallicity, colour-coded for the stellar mass surface densities
(expressed on a logarithmic scale). Galaxies are divided in two mass bins:
$\log(M_{\star}/M_{\odot})<11$ (left), and $\log(M_{\star}/M_{\odot})>11$
(right). Brown markers and shaded regions represent the binned medians and
16th-84th percentile range of the resulting distribution for morphologically
classified sub-samples of dusty QGs: those classified as elliptical-type QGs
(eQGs; squares) and spiral QGs (sQGs; circles). The typical error in gas-
metallicity estimates for solar and super-solar values are displayed with
horizontal lines, respectively.

##  5 Morphological impact on the evolution of dust abundance in QGs

Galaxy morphology is often directly linked to some of the dust-related
properties of galaxies (e.g., Smith et al. 2014, Nersesian et al. 2019,
Casasola et al. 2020). Our goal in this Section is to explore if the observed
evolution of specific dust mass at intermediate redshifts is driven by the
morphological type of our galaxies. For this reason we adopt the morphological
classification by Tasca et al. (2009). They infer morphological types of
galaxies in COSMOS by applying the automated method on parameters such as
Concentration ($C_{\rm n}$), Gini index ($G$), Asymmetry and M20 moment,
calculated from the HST F814W imaging (Cassata et al. 2007). With this we
split the full sample of our QGs into two sub-classes: elliptical QGs (eQGs;
385 objects out of 548; $70\%)$ and spiral QGs (sQGs; 163 objects out of 548;
$30\%$)555We note that morphological classification remains almost unchanged
when using other methods, i.e., Zurich Estimator of Structural Types (ZEST) by
Scarlata et al. (2007). The distributions of their main structural and
physical properties are presented in Fig. 4. Through the rest of the paper, we
extensively investigate evolution of specific dust mass with different
parameters in these two morphological groups.

###  5.1 The relationship between $M_{\rm dust}/M_{\star}$ and gas-phase
metallicity

Gas-phase metallicity is considered one of the crucial parameters influencing
the dust life cycle in galaxies (e.g., Asano et al. 2013, Pantoni et al. 2019,
Hou et al. 2019, Triani et al. 2020). The unique aspect of our hCOSMOS dataset
is a combination of self-consistently SED-derived $M_{\rm dust}$ and
$M_{\star}$, and independent measurements of $Z_{\rm gas}$. This enables us to
investigate, for the first time, how the $M_{\rm dust}/M_{\star}$ evolves as a
function of $Z_{\rm gas}$ in two different morphological categories of QGs
beyond the local Universe. The majority of our objects (76%) have ”super-
solar” $Z_{\rm gas}$ with the median that reaches twice the solar value
($12+\log(\rm O/H)=9.01^{+0.07}_{-0.13}$). Both distributions and medians of
$Z_{\rm gas}$ are similar between the two morphological groups (see Fig. 5).
In Fig. 5. we showcase the relation of $M_{\rm dust}/M_{\star}$ with $Z_{\rm
gas}$ for two stellar mass bins, colour-coding their stellar mass surface
densities (defined as $\Sigma_{\rm M}=M_{\star}/2\pi R_{\rm eff}^{2}$, where
$R_{\rm eff}$ is effective radius in kpc based on the measurements from
Cassata et al. 2007).

The Fig. 5 reveals several interesting features: (1) At fixed $Z_{\rm gas}$
the observed $M_{\rm dust}/M_{\star}$ varies amongst morphological types. In
sQGs $M_{\rm dust}/M_{\star}$ is systematically higher than in eQGs, and it
evolves in a more complex way. When transiting from (sub-)solar to the super-
solar gas metallicities, sQGs undergo turnover of the relation between $M_{\rm
dust}/M_{\star}$ and $Z_{\rm gas}$. Such behaviour is visible in both mass
bins, being more prominent for $10<\log(M_{\star}/M_{\odot})<11$. We note that
the trend observed in sQGs should be consider rather tentative, as typical
uncertainties in $Z_{\rm gas}$ measurements are quite large ($\sim 0.4$) for
$12+\log(\rm O/H)\leq 8.8$. Irrespectively of this, eQGs maintain remarkably
flat $M_{\rm dust}/M_{\star}$ over range of metallicities. To our knowledge,
this is the first time that the morphological dependence in $M_{\rm
dust}/M_{\star}$ with $Z_{\rm gas}$ has been observed in QGs. This result
strongly indicates that morphology-type impact on the dust content extends at
least to $z\sim 0.6$, complementing the conclusions from known studies of QGs
at $z\sim 0$ (Smith et al. (2014), Beeston et al. 2018, Nersesian et al. 2019,
Casasola et al. 2022); (2) For $12+\log(\rm O/H)\sim 9$ there is a range of
specific dust masses that extends over $\sim 2$ orders of magnitude and
exceeds the median deviation in both morphological types. This implies that
our QGs are subject to complex interplay of processes contributing to the dust
life cycle.

In the first place, dust particles can be internally produced via core-
collapse supernovae (CC SNe) and in outflows of AGB stars, and can further
grow via collisions and accretion of free metals in the ISM. The absence of a
clear trend in Fig. 5 suggests that balance between formation and destruction
of dust grains is altered at super-solar $Z_{\rm gas}$. Specific dust masses
in eQGs saturates around $M_{\rm dust}/M_{\star}\sim 10^{-4}$ over range of
metallicities. Interpreting such behavior is not trivial. For evolved
galaxies, the $M_{\rm dust}/M_{\star}$ is expected to arise from the
”competition” between the dust-to-gas ratio ($\delta_{\rm DGR}$) increasing
and $f_{\rm gas}$ decreasing (i.e., Asano et al. 2013, Béthermin et al. 2015).
By using resolved observations of M101 galaxy, Vílchez et al. (2019) find
constant $M_{\rm dust}/M_{\star}$ with $Z_{\rm gas}$ in the outer part
characterised by sub-solar $Z_{\rm gas}$. They interpret such flat trend as a
consequence of constant yield ratio dominated by stellar sources under
constant gas fraction ($f_{\rm gas}\sim 0.1$). Nevertheless, such scenario
seems unlikely for our sample, as it is favoured for sub-solar environments,
but not for the metal-rich ISM. Furthermore, cold gas fraction in our QGs is
expected to vary, consequence of the cosmic evolution. If the dust yield is
solely by stellar sources and subsequently removed by efficient outflows, one
would expect strong anti-correlation between $M_{\rm dust}/M_{\star}$ and
$Z_{\rm gas}$, as known in most of local, DSFGs (i.e, Casasola et al. 2022).
This contradicts our data and suggests that prolonged dust productions of
varying efficiencies may be required to balance destructive processes
responsible for gas (and subsequently, dust) decline.

In this regard, growing number of theoretical studies propose the grain growth
of dust in ISM to play important role in metal-rich galaxies (e.g., Asano et
al. 2013, Rémy-Ruyer et al. 2014, Hirashita et al. 2015, Zhukovska et al.
2016, De Vis et al. 2019, Aoyama et al. 2019). In such scenario, if $Z_{\rm
gas}$ in a galaxy exceeds a certain critical value which strongly depends on
galaxy star-formation history, the grain growth becomes active and the $M_{\rm
dust}$ rapidly increases until metals are depleted from the ISM. The evolution
of specific dust mass in our sQGs qualitatively agrees with the predictions
from the grain growth scenario. $M_{\rm dust}/M_{\star}$ first rises up to
certain $Z_{\rm gas}$, and levels off when most of the available metals are
depleted onto the dust grains. At the same time, global dust destruction
should decrease available cold gas, reducing the material for the further dust
production. This interpretation was proposed in Casasola et al. (2020) who
explore morphological impact on the scatter in dust scaling relations with
$Z_{\rm gas}$ in Dustpedia galaxies. They invoke effect of grain growth to
explain observed specific dust masses ($M_{\rm dust}/M_{\star}>0.0001$) and
high average dust-to-gas ratio ($\delta_{\rm DGR}\gtrsim 1/100)$) in objects
with super-solar $Z_{\rm gas}$. In principle, expected accretion time for
grain growth scales to the inverse of the cold gas fraction, being very rapid
($\sim 10^{7}-10^{8}$ yr) in QGs with $f_{\rm gas}>0.1$ (Hirashita et al.
2015, Hirashita & Nozawa 2017). With this production channel being dominant,
the observed scatter in $M_{\rm dust}/M_{\star}$ at a given $Z_{\rm gas}$ may
reflect the variation of growth efficiencies. The significance of such process
can be traced via dust depletion to metals (quantified as ”dust-to-metal
ratio”, $\delta_{\mathrm{DTM}}$, De Vis et al. 2019), which we explicitly
investigate with models in Section 6.4.

Our QGs span a wide range of optical sizes, with stellar mass surface
densities (defined as $\Sigma_{\rm M}=M_{\star}/2\pi R_{\rm eff}^{2}$)
extending over $\gtrsim 2$ orders of magnitude. From Fig. 5 we see that over
range of metallicities, dusty QGs with lower surface densities tend to have
larger $M_{\rm dust}/M_{\star}$. This may suggests that processes driving the
size growth in the QGs may leave observational imprints on their relative dust
abundance. Because the presence of a companion galaxy can affect the dust
scaling relations, several works argue that i.e., extended sizes of dusty QGs
can support the idea that some fraction of dust and metals is acquired via
external sources, most likely via mergers with dust-rich satellites (Rowlands
et al. 2012, Martini et al. 2013, Lianou et al. 2016, Dariush et al. 2016,
Kokusho et al. 2019). We will return to the connection with surface densities
in the sections where we examine implications related to the channels of dust
production and dust removal, and their link with post-quenching timescales
(Section 5.4 and Section 5.5). The comprehensive theoretical interpretation
based on the models will be presented in Section 6.

###  5.2 Dissecting the age-dependent evolution of $M_{\rm dust}/M_{\star}$

We now focus on analysing the $M_{\rm dust}/M_{\star}$ under the age-
evolutionary framework. To achieve this goal, we adopt the mass-weighted
stellar ages ($t_{\star}$), derived from our self-consistent SED fitting. In
Fig. 6 we show the trends of $M_{\rm dust}/M_{\star}$ with the mass-weighted
stellar ages ($t_{\star}$) for two morphological subgroups. From Fig. 6 we see
that $M_{\rm dust}/M_{\star}$ in sQGs has the overall decreasing trend with
increasing stellar age, opposite to eQGs that maintain very shallow evolution.
The median stellar age of sQGs is, on average, $\sim 1$ Gyr younger than in
eQGs ($\log(\rm Age/yr)=9.65$ and $\log(\rm Age/yr)=9.79$, respectively). The
large scatter is especially pronounced at older ages (e.g. $t_{\star}>5$ Gyr),
while the difference between the two sub-groups almost vanishes towards the
largest stellar ages.

Figure 6: Evolution of $M_{\rm dust}/M_{\star}$ as a function of mass-weighted
stellar age, colour-coded for the stellar-mass surface density. Brown symbols
have the same meaning as in the previous Figure. Dotted line is the best fit
(exponentially declining function) to the sQGs.

Inferred mass-weighted stellar ages inform that the majority of stellar
content was produced in early times in both morphological groups. There is a
tail towards younger ages visible for the sQGs, while the same in not observed
for the eQGs. The large spread in $t_{\star}$ as seen for our sQGs is in
agreement with Tojeiro et al. (2013) and Zhou et al. (2021), who found that
some fraction of sQGs at low-$z$ have exceptionally younger ages (¡0.5 Gyr)
than the general population of eQGa and sQGs. By investigating the mass-size
growth in the full sample of hCOSMOS QGs, Damjanov et al. (2022) show that the
non-negligible number of extended QGs with younger stellar populations (as
seen in our dustiest sQGs) are recently quenched galaxies. Interestingly, they
found the evolution of the mass-size relation for those recently quenched QGs
to be mostly a consequence of progenitor bias rather than of minor mergers.

Negative correlation between $M_{\rm dust}/M_{\star}$ and stellar age (as seen
in our sQGs) can be interpreted as age-evolutionary sequence (Michałowski et
al. 2019, Burgarella et al. 2020). Within this framework, the level of
specific dust decline with $\tau_{\star}$ can be used for quantifying the
removal time for dust. For the small sample of H-ATLAS red spirals and
ellipticals at $z\sim 0.1$, Michałowski et al. (2019) found the exponentially
declining function to the $M_{\rm dust}/M_{\star}$ with age, inferring long
dust removal timescales ($\sim 2.25$ Gyr). If we apply the same method to our
sQGs, we would obtain $\tau_{\rm d}=2.75\pm 0.5$ Gyr, which agrees within
$1\sigma$ with estimates from Michałowski et al. (2019), despite some
noticeable differences between the two samples.666The sample from Michałowski
et al. (2019) is composed of H-ATLAS objects that are described to have very
bright SPIRE fluxes. In addition, their sample also contains ”bluer” objects
with $D_{\rm n}4000$ smaller than the cut we applied to systematically select
our QGs ($D_{\rm n}4000$¿1.5). Contrary to sQGs, the evolution of $M_{\rm
dust}/M_{\star}$ in eQGs is almost independent of stellar age, indicative of
very slow dust removal, which likely be prolonged over several Gyr. Weak
evolution of $M_{\rm dust}/M_{\star}$ at high stellar ages, and especially in
eQGs, can strengthen the idea presented in previous Section that dust
abundance in our QGs could be of mixed origin. For example, flattening of
specific dust mass with age is expected if dust is continuously produced by
AGB stars or is supported externally via mergers (see e.g., Smith et al. 2014,
Richtler et al. 2018, Kokusho et al. 2019). An additional clue on this can be
gained from mass surface densities. Many older QGs are generally expected to
have undergone the size growth via minor mergers, which would increase the
$M_{\star}$, but decrease the surface density. This would introduce the
scatter in $\Sigma_{\rm M}$ per given age, in line to what is seen at Fig. 6.
Optical studies of QGs at low and intermediate redshifts (i.e., Beverage et
al. 2021, Barone et al. 2022) pointed out that the scatter in surface
densities per given age reflects the diverse conditions of the Universe when a
galaxy becomes quiescent.

###  5.3 The evolution of $M_{\rm dust}/M_{\star}$ with $M_{\star}$

Our data indicate a general anti-correlation between the $M_{\rm
dust}/M_{\star}$ with $M_{\star}$. As evident from Fig.7, for given
$M_{\star}$, there is a distinction between the sQGs and eQGs at
$\log(M_{\star}/M_{\odot})\lesssim 10.75$, with very little overlap between
the two. Strong difference vanishes above this stellar mass with many sQGs
entering the region populated by eQGs ($M_{\rm dust}/M_{\star}<10^{-4}$). The
anti-correlation with $M_{\star}$ in sQG and eQGs can be best described with
following fitting functions:

$\log(M_{\rm dust}/M_{\star})=-0.58\times M_{\star}+2.71$ (2)

for sQGs, and

$\log(M_{\rm dust}/M_{\star})=-0.26\times M_{\star}-1.11$ (3)

for eQGs. The Spearman rank coefficient is find to be -0.62 for sQGs, and
-0.32 for eQGs, while for both morphological groups we obtain a small
probability (corresponding to the high significance of $\sim 5.5\sigma$ and
$4.2\sigma$, respectively) that there is no correlation. Recent works on local
dusty galaxies reported morphological impact on the relation between $M_{\rm
dust}/M_{\star}$ with $M_{\star}$ (Nersesian et al. 2019, Casasola et al.
2020). Interestingly, the evolutionary slope we infer for our eQGs is
identical to the one found in studies of nearby SF galaxies (Casasola et al.
2020), and in ALMA detected MS DSFGs at $z\sim 1$ (Donevski et al. 2020). This
suggests that DSFGs and eQGs likely share some common processes in their dust
evolution, i.e., removal and/or quenching (Lapi et al. 2018).

Figure 7: Upper panel: The evolution of $M_{\rm dust}/M_{\star}$ with
$M_{\star}$ in our QGs, colour-coded as a function of galaxy size (represented
with effective radius, $R_{\rm eff}$, in kpc). Median values in different
stellar mass bins are displayed with brown circles and squares for sQGs and
eQGs, respectively. The typical errors are shown as brown crosses. The
straight lines are the best linear fits describing the data for sQGs and eQGs.
Line symbols and fitting functions are indicated in the legend, while the
shaded regions represent $95\%$ confidence interval of the fit; Lower panel:
The evolution of $M_{\rm dust}/M_{\star}$ with $M_{\star}$ in sQGs (left) and
eQGs (right), in the connection with time after quenching ($t_{\rm quench}$).
Coloured lines are describing evolution of specific dust mass (as median
values in three stellar mass bins) for different post-quenching intervals. The
meaning of coloured lines is indicated in the legend.

Negative correlation of $M_{\rm dust}/M_{\star}$ with $M_{\star}$ observed in
our QGs is in agreement with known observational studies of dusty galaxies in
different environments (e.g., Bourne et al. 2012, Casasola et al. 2020,
Donevski et al. 2020). It is usually interpreted as a natural reflection of
the dust life-cycle: $M_{\star}$ grows with time as galaxy evolve, while dust
grains decrease from the budget being incorporated into the stellar mass. By
applying the chemical galaxy model on proto-spheroidal galaxies, Calura et al.
(2017) find that $M_{\rm dust}/M_{\star}$ is larger in galaxies with
$10<\log(M_{\star}/M_{\odot})<11$ because these are characterised both by a
lower specific destruction rate and by a larger growth rate than in more
massive galaxies. This may explain why our QGs with relatively high specific
dust mass ($\log(M_{\rm dust}/M_{\star})>-3$) are identified exclusively at
$\log(M_{\star}/M_{\odot})<11)$. Casasola et al. (2022) recently analyse the
spatially resolved Dustpedia galaxies and show that variations in $M_{\rm
dust}/M_{\star}$ exist even for galaxies within the same morphological group,
attributing it to the variations in dust growth efficiencies.

There is another important aspect of our analysis presented in Fig.7. Namely,
the uppermost part of the plane at $\log(M_{\star}/M_{\odot})\lesssim 11$ is
dominated by sQGs that are characterised with larger sizes and stellar
populations younger than in eQGs with same $M_{\star}$. Further increase in
optical sizes towards $\log(M_{\star}/M_{\odot})\gtrsim 11$ is followed by
reduced difference between the specific dust masses of eQGs and sQGs. That
said, the growth in galaxy size with $M_{\star}$ (as confirmed by many works
on the mass-size relation) does not correspond to the similar increase in
$M_{\rm dust}$. Damjanov et al. (2022) demonstrate that sizes of QGs from this
work increase as $\log(R_{\rm eff}/\rm kpc)\sim
0.75\times\log(M_{\star}/M_{\odot})$. While this evolution is find to be
almost redshift independent, it emerges from different processes, both
internal and external. Such dichotomy can modulate gas fraction and quenching
efficiency which in turn would impact evolution of a dust content in galaxies.

###  5.4 How the dust evolution and post-quenching timescales are connected?

One of the main questions that arise from presented results is: How can we use
the observed evolution of dust content to learn about evolution of QGs
following their quenching? To answer this question, from star-formation
histories modelled in CIGALE, we adopt post-quenching time ($t_{\rm quench}$),
which is the time passed since the star formation is quenched, and is defined
as age of the galaxy minus $t_{\rm trunc}$ from Eq.1. Our modelling procedure
reveals wide distribution of $t_{\rm quench}$, from 60 Myr to 3.2 Gyr, which
ensures that QGs from this work probe broad dynamical ranges in terms of
quenching age. The median $t_{\rm quench}$ is found to be $436\pm 215$ Myr for
sQGs, and $704\pm 282$ Myr for eQGs. Lower panel in Fig. 7 shows that QGs
generally experience trend of $M_{\rm dust}/M_{\star}$ declining with post-
quenching time increasing, which applies for both morphological groups.
Morphology impact on $M_{\rm dust}/M_{\star}$ is evident by comparing the two
sub-populations within the same range of $M_{\star}$ and $t_{\rm quench}$. The
difference in $M_{\rm dust}/M_{\star}$ between the sQGs and eQGs is the
largest in ”recently-quenched” QGs ($t_{\rm quench}<500$ Myr), and smallest
for ”early-quenched” galaxies ($t_{\rm quench}>\rm 1Gyr$). Nevertheless, the
average vertical drop in specific dust mass per given $M_{\star}$ is
relatively small within the morphological groups: it is $\sim 0.5\>\rm dex$
for sQGs and even lower ($\sim 0.25\>\rm dex$) for eQGs. Therefore, the
observed trend of $M_{\rm dust}/M_{\star}$ with $M_{\star}$ is mainly driven
by increase of $M_{\star}$ with $t_{\rm quench}$, while $M_{\rm dust}$
exhibits remarkably shallow decline within the first $\sim 1$ Gyr following
quenching.

Considering above, we can draw several conclusions: (1) Our data provide
strong evidence for morphology playing a role in observed $M_{\rm
dust}/M_{\star}$. Under the assumption that the level of SFR cessation
correlates with $t_{\rm quench}$, we can deduce that for the majority of our
QGs the removal of cold gas and dust are interlinked. However, our data
strongly suggest the link is non-monotonic, which is manifested as flattening
of trend between $M_{\rm dust}/M_{\star}$ and $M_{\star}$ within certain
$t_{\rm quench}$; (2) Such shallow evolution in $M_{\rm dust}/M_{\star}$ with
increasing $M_{\star}$ cannot be explained with age-evolutionary scenario, and
is rather indicative of ”dust abundance recovery” counterbalancing the drop of
star-formation over different timescales. The ”dust abundance recovery” in
sQGs happens on timescales shorter than in eQGs ($t_{\rm quench}<500$ Myr vs.
$t_{\rm quench}>500$ Myr, respectively). The difference may be caused by
initial cold gas available after quenching which would affect the accretion
timescales (see Section 5.5). Evidence of flattening can be seen even in the
most massive QGs that quenched $>1Gyr$ ago, but these objects have lower than
average $M_{\rm dust}/M_{\star}$, which is likely a consequence of ISM removal
mechanisms (e.g. gas heating, Smercina et al. 2018) starting to be visible
against the prolonged dust production.

Dusty QGs from this work may serve as a benchmark for better understanding the
connection between their structural properties and the ISM evolution. For
instance, if $M_{\rm dust}/M_{\star}$ mirrors the change of molecular gas
fraction (e.g., Imara et al. 2018, Whitaker et al. 2021a, Magdis et al. 2021),
results from Fig. 7 implies that $f_{\rm gas}$ is the highest soon after the
quenching episode, and that sQGs have on average a higher cold gas content
than eQGs. Within both morphological categories, more extended galaxies have
larger specific dust mass. Chen et al. (2020) and Barone et al. (2022) argue
that extended QGs with smaller surface densities undergo less extreme kinetic
mode feedback (and consequently slower quenching), than similarly massive
compact sources, which would leave them with more residual molecular gas.

Furthermore, the effective presence of dust in our QGs along with inferred
$t_{\rm quench}$ may provide important clues on quenching mechanisms. As
evolution in $M_{\rm dust}/M_{\star}$ changes its pace (Fig. 7 reveals phases
that are both accelerated and slowed down within the same range of $t_{\rm
quench}$), it is likely no one mechanism is responsible for quenching star
formation in our QGs. This is consistent with observational studies of passive
spirals (Fraser-McKelvie et al. 2018) and passive ellipticals (Zhou et al.
2021). It is also consistent with the TNG-50 simulation’s study which reported
variety of quenching timescales and time since quiescence, ranging from 1-4
Gyr in eQGs and sQGs (Park et al. 2022). Similarly, Mahajan et al. (2020)
found that local dusty QGs from GAMA survey are primarily affected by
quenching acting on longer timescales (order of several Gyr) due to a lack of
morphological transformation associated with the transition in optical colour.
Ciesla et al. (2021) demonstrate how variations in dust properties (i.e.
$L_{\rm dust}$) per given $t_{\rm quench}$ resemble the ”strength” of
quenching processes in a small sample of local dusty QGs from the Herschel
reference survey (HRS). They proposed how the dust abundance observed in
objects that quenched $\sim 1$ Gyr ago can be explained if quenching is caused
by ram-pressure stripping, which is known for affecting mostly gas, but not
stars and dust in the central regions. The same mechanism is recently found to
be responsible for galaxy size increase (Grishin et al. 2021) and some of our
most extended sQGs with the highest specific dust masses appear to be
compatible with this scenario. Interestingly, the average $M_{\rm
dust}/M_{\star}$ as well as the large spread with stellar age, are comparable
to those of dusty post-starbursts for which Li et al. 2019b reported offset of
several hundred Myr between the SFR quenching and decline in gas-fraction and
dust mass.

###  5.5 Implications for the production of dust in QGs

Results presented throughout the Section 5 impose important constraints on the
timelines and mechanisms of dust production and dust removal in QGs. Below we
summarise possible scenarios, and continue with analysis with models in
Section 6.

Implications for internal dust growth: In galaxies dominated by older stellar
populations, there will be time for AGB stars to contribute significantly to
the dust mass budget. The total mass of dust formed by AGB intermediate mass
stars ($1\>M_{\odot}\lesssim M_{\star}<8\>M_{\odot}$) depends on the initial
mass of the star, with the maximal yield proposed to lay within $\sim 4\times
10^{-4}$ and $\sim 3\times 10^{-3}\>M_{\odot}$ for the case of solar and
super-solar metallicities, respectively (Zhukovska et al. 2008, Ventura et al.
2020). By imposing maximal production in $\sim 10^{10}$ AGB stars, we would
obtain $M_{\rm dust}$ between $4\times 10^{6}M_{\odot}$ and $3\times
10^{7}M_{\odot}$. However, the average production efficiency by the AGB stars
is at least 10-20 times (depending on the stellar model) lower than the
maximal one (Schneider et al. 2014), resulting in dust masses $\lesssim
3\times 10^{6}M_{\odot}$. This is still below the average value of eQGs,
suggesting insufficient dust yield solely produced by AGB stars to fully match
our data.777It is worth noting here that this does not rule out contribution
from AGB stars. Under certain conditions, they may be important suppliers of
prolonged dust production (especially in eQGs), but it is unlikely for them to
be the dominant dust production channel. On the other hand, average value in
sQGs is even higher ($M_{\rm dust}>2.6\times 10^{7}M_{\odot}$), difficult to
explain even with maximally efficient SNe dust production in the early stages
of galaxy history, as demonstrated by several works (i.e., Michałowski et al.
2019). Instead, our data favour dust grain growth as a relatively fast
production mechanism to overcome the dust loss. As mentioned in Section 5.1,
in metal-rich galaxies this process can increase $M_{\rm dust}$ until metals
are depleted from the ISM. Timescales for dust growth in QGs are proposed to
be relatively rapid, up to $\sim 50-250$ Myr, if the lifetime of the cold gas
is long enough ($\sim 10^{7}$ yr, Hirashita et al. 2015). If we stay
conservative and account for the factor of $\sim 3$ uncertainty due to
expected variations in the $f_{\rm gas}$ of our QGs, accretion timescales
would remain shorter than $t_{\rm quench}$ for most of our galaxies. This may
explain why the dust emission in our QGs is detectable in timescales longer
than those expected for dust destruction due to gas heating ($\sim 10^{8}$
yr).

Implications for external sources of dust: The most massive of our eQGs lack
the clear evolution of $M_{\rm dust}/M_{\star}$ with stellar ages, metallicity
and stellar mass. This qualitatively agrees with scenarios proposing external
contribution to $M_{\rm dust}$ from minor mergers (Kartaltepe et al. 2010,
Rowlands et al. 2012, Martini et al. 2013, Lianou et al. 2016), Richtler et
al. 2018, Kokusho et al. 2019). Our sample contains only 17 QGs ($4\%$) with
”disturbance” signatures (extended objects with morphologies that are on the
edge between disks and irregulars). However, this should be considered a lower
limit since there is a possibility that mergers occur long time ago enough to
erase the clear disturbance signature. We thus follow the approach by Nevin et
al. (2019) who investigate simulated mergers (both gas-poor and gas-rich)
within a wide range of mass ratios. They leverage the strength of
concentration coefficient and Gini index as the most sensitive predictors for
minor mergers in their later stages, since these parameters generally enhance
as the merger progresses. If we adopt the same criteria proposed by Nevin et
al. (2019) (their Fig. 5), we would end up with 86 merger candidates or $16\%$
of our sample. Such criteria correspond to galaxies in their final, post-
coalescence stage, $t>3.5$ Gyr since the merging started. Interestingly, the
same authors find that the difference in $f_{\rm gas}$ between gas-poor and
gas-rich mergers does not produce a significant distinction in the
concentration of the remnant. Thus, we caution that candidate merging QGs in
our sample can be a mix of both, gas-rich and gas-poor. In principle, galaxies
that show strong positive excess in their $M_{\rm dust}/M_{\star}$ are likely
to be compatible with the gas-rich scenario. The majority of our minor merger
candidates (75/86) are eQGs of older ages ($>7-10$ Gyr) and higher stellar
masses ($\log(M_{\star}/M_{\odot})>10.9$). Only 11/86 candidates are sQGs. As
a result, we anticipate that internal mechanisms will dominate dust content in
QGs, with minor mergers contributing more to eQGs (75/386, or $\sim 19\%$)
than sQGs (11/162, or $6\%$) We check the catalogue of merging Spitzer
galaxies in the COSMOS field (Kartaltepe et al. 2010) and find 14 objects from
our sample that are present in this catalogue, and are classified as minor
mergers. These QGs have $M_{\rm dust}/M_{\star}$ fluctuating within $\pm 0.3$
dex of the median of our full sample. In addition, there is a tentative
evidence in our data that candidate mergers have $Z_{\rm gas}$ $\sim 0.1$ dex
below the median, consistent with suggestions from the literature of mergers
causing the scatter in mass-metallicity relation (Bustamante et al. 2018,
Griffith et al. 2019). This all support the idea that merging environment can
have visible impact on the evolution of $M_{\rm dust}/M_{\star}$ with
$M_{\star}$, $Z_{\rm gas}$ and stellar age in QGs. Recent observational tests
conducted with ALMA supported such possibility through identifying ISM
complexity in a small sample of the brightest dusty eQGs at $z\sim 0.05$
(Sansom et al. 2019). This study reveals examples of massive molecular gas
discs, but also objects with dusty companions contributing to $>50\%$ of the
FIR emission.

###  5.6 Implications for the survival of dust in massive QGs

The lack of evolution in specific dust mass of eQGs with $Z_{\rm gas}$, age
and $M_{\star}$ is compatible with very slow dust removal. This can be a
result of weak outflows, long destruction timescales due to SN shocks and
inefficient thermal sputtering. Below we briefly examine some of the
possibilities.

  * •

Thermal sputtering
Galaxies in massive halos and clusters are expected to undergo significant
thermal destruction of their grains. As presented in Donevski et a. (in prep),
$\sim 15\%$ of QGs from this work reside in galaxy clusters. The observed
$M_{\rm dust}/M_{\star}$ in these sources would be difficult to interpret if
the sputtering timescale due to the exposition of dust grains to a hot ISM is
extremely short. Timescales within which dust grains are efficiently destroyed
in QGs are given as $10^{5}(1+(10^{6}K/T_{\rm gas})^{3})n_{e}^{-1}$ yr for
typical 0.1 $\mu$m-sized grains (Hirashita & Nozawa 2017). While the typical
values are quite short (1-10 Myr, Vogelsberger et al. 2019), it would be
possible to obtain longer destruction timescales ($t_{\rm sput}>1$ Gyr) if
grains are larger and/or gas is of lower density (Smercina et al. 2018). Such
longer sputtering timescales of $t_{\rm quench}>0.5-1$Gyr are another
indication for prolonged dust growth as it supports the production of large
grains which are expected to be reduced relatively slowly in hot halos (Relaño
et al. 2022, Priestley et al. 2022). On top of this, along with the sputtering
in hot halos, dust cooling may be operating through the emission of IR
radiation, as proposed by recent hydrodynamical simulations (Vogelsberger et
al. 2019). However, this solution would likely produce excess in the sub-mm
wavelengths, which we did not observe in existing SCUBA-2 data.

  * •

Destruction due to SN shocks
The total rate of dust mass destruction due to SN shocks is given by
$\dot{M}_{\rm dest}\propto M_{\rm dust}/\tau_{\rm destr}$, where dust
destruction timescale is usually approximated as (Slavin et al. 2015):

$\tau_{\rm destr}=\frac{\Sigma M_{\rm gas}}{f_{\rm ISM}R_{\rm SN}M_{\rm
cl}}=\frac{\tau_{\rm SN}M_{\rm gas}}{M_{\rm cl}}$ (4)

Here $\Sigma M_{\rm gas}$ is the surface density of molecular gas mass,
$f_{\rm ISM}$ is the value that accounts for the effects of correlated SNe,
$R_{\rm SN}$ is the SNe rate, $\tau_{\rm SN}$ is the mean interval between
supernovae in the Galaxy (the inverse of the rate) and $M_{\rm cl}$ is the
total ISM mass swept-up by a SN event. Here is important to note that $M_{\rm
cl}$ varies with the ambient gas density and metallicity, and as metals being
efficient cooling channel in the ISM, higher $Z_{\rm gas}$ would result in
smaller swept mass (see Asano et al. 2013, Hou et al. 2019). We can roughly
approximate range of destraction timescales by using estimated $Z_{\rm gas}$,
and applying conservative approach for converting $M_{\rm dust}$ to $M_{\rm
gas}$ through metallicity-dependent $\delta_{\mathrm{DGR}}$ (Rémy-Ruyer et al.
2014). Following Eq 4, we infer range of $\tau_{\rm destr}$ from $0.11\>\rm
Gyr$ to $9.5\>\rm Gyr$, with the median of $1.28\>\rm Gyr$ for sQGs, and
$1.75\>\rm Gyr$ for eQGs. This informs that destruction due to SNe shocks in
not expected to be crucial in our sample.

Overall, the brief analysis from this section hints at an expected small
influence of dust destruction due to SNe shocks, while effect of thermal
sputtering can be significantly prolonged if grains are growing inside QGs. On
top of this, effect of outflows (e.g., due to stellar and X-ray feedback) is
expected to be more important (i.e., De Looze et al. 2016). To investigate
this, in the next Section we analyse outputs from different simulations of
dust formation that have these effects included.

##  6 Interpretation of results with models

To evaluate our results within the framework of dusty galaxy formation, in
this section we inspect analytic models and the state-of-the-art cosmological
simulations that track the dust life cycle in a self-consistent way.

###  6.1 Models

We analyse the predictions from the: (I) cosmological galaxy formation
simulation with self-consistent dust growth and feedback (SIMBA, Davé et al.
2019); and (II) chemical model of Nanni et al. (2020).

####  6.1.1 SIMBA cosmological simulation

The cosmological galaxy formation simulation, SIMBA, utilizes mesh-free finite
mass hydrodynamics, and is successor of MUFASA simulation (Davé et al. 2016)
with improvements to the sub-grid prescriptions for star formation, AGN
feedback and X-ray feedback, as well as implementation of dust physics. We
refer the reader to Davé et al. (2019) for extensive description of the
simulation. Self-consistent dust framework that models the production, growth
and destruction of dust grains in SIMBA is introduced in Li et al. (2019a).
Overall, SIMBA accounts for dust grains produced from stellar populations (SNe
and AGB stars), which can further grow via accreting gas-phase metals. Dust
can be destroyed via different processes such as thermal sputtering, SNe
shocks and astration (consumption by star-formation). Dust can instantaneously
be destroyed in gas impacted by jet-mode AGN feedback, while outflow processes
such as radiative-mode Eddington AGN feedback and stellar feedback can heat up
and transport dust out of the galaxy.

The net rate of dust production/destruction in SIMBA can be generalised as:

$\Sigma\dot{M}_{\rm dust}\propto\dot{M}_{\rm dust}^{\rm SNe}+\dot{M}_{\rm
dust}^{\rm ISM}-\dot{M}_{\rm dust}^{\rm destr}-\dot{M}_{\rm dust}^{\rm
SF}+\dot{M}_{\rm dust}^{\rm inf}-\dot{M}_{\rm dust}^{\rm out}.$ (5)

The first term on the right side of Eq. 7 describes the dust produced by
condensation of a fraction of metals from the ejecta of $\rm SNe$ and $\rm
AGB$ stars; the second term describes the dust by accretion in the ISM; the
third term describes the dust destructed by $\rm SNe$ shock waves and thermal
sputtering; the fourth term is the consumption of dust due to astration; the
fifth term is an additional dust production by gas infall; the sixth term
describes the dust mass that is ejected out of the ISM due to feedback
processes. The latter two mechanisms are responsible for heating up and
removal of gas into the DM halo, or even further out. The dust model within
SIMBA does not include contribution from $\rm Ia\>SNe$ which is opposite to
some models that proposed the same condensation efficiency between Type $\rm
Ia\>SNe$ and Type II SNe (Dwek 1998, McKinnon et al. 2017, Popping et al.
2017). The condensation efficiencies for AGB stars and and $\rm CC\>SNe$ are
updated based on the theoretical models of Ferrarotti & Gail (2006) and
Bianchi & Schneider (2007). The evolution of dust grains in ISM is simulated
with single-sized dust model where the grains are assumed to all have the same
initial radius and density ($a=0.1\>\mu$m and $\sigma=2.4\>g/cm^{3}$,
respectively). SIMBA was run on a number of volumes with different
resolutions. Here we use the largest one which has $1024^{3}$ dark matter
particles and $1024^{3}$ gas elements in a cube of $140\>{\rm Mpc}h^{-1}$ side
length. This simulation run includes the full feedback and dust physics, which
is ideal for our goal of studying the statistical sample of dusty QGs. We note
that this SIMBA run results in dust mass functions and dust-to-metal ratios in
a good agreement with observations at $z<1$ (Li et al. 2019a), and is able to
explain the observed sub-mm number counts up to $z\sim 2$ (Lovell et al. 2021)
and the dust content of QGs observed at $z\sim 1-2$ (Williams et al. 2021).

####  6.1.2 The Nanni et al. (2020) model

The chemical model of Nanni et al. (2020, hereafter N20) introduces a new set
of analytical solutions for the dust evolution in galaxies. The model assumes
that dust formed in supernovae remnants, around evolved AGB stars and in the
ISM is composed predominantly by silicates (olivine and pyroxene), amorphous
carbon dust and metallic iron. Similar to SIMBA, Type Ia SNe are not
considered important contributors to the dust enrichment. In N20 different
physical processes are considered for affecting the time evolution of $M_{\rm
dust}$: astration of gas and metals due to formation of stars; dust
destruction by SN shock waves that propagate in the ISM, and dust removal by
galactic outflows. The model does not include thermal sputtering of dust in
hot halos. The following theoretical yields are adopted: Kobayashi et al.
(2006) for AGB stars, Cristallo et al. (2015) for Type II SNe, and Iwamoto et
al. (1999) for Type Ia SNe Consistently with the calculations of gas and metal
evolution, N20 accounts for the inflow material being composed of pristine gas
that does not change the total content of dust. At each time-step the galactic
outflow is considered to be regulated through the ”mass-loading factor” (ML)
and assumes the outflow to be due to the stellar feedback in the ISM. The ML
parameter is quantitatively defined as the ratio of outflow rate and SFR ($\rm
ML=\dot{M}_{\rm out}/\rm SFR$) and is related to initial baryon mass in the
galaxy, that is fully composed of gas at the beginning of simulation and
converted to stars afterwords. In N20 model, $M_{\rm gas}$ is set as a
multiple of $M_{\star}$ normalised to 1$\>M_{\odot}$ at 13 Gyr. While the
fiducial N20 model excludes explicit growth of dust particles in ISM, here we
test the importance of this mechanism by switching it on and off in different
simulation runs. It is worth noting that N20 model showed success in
reproducing total dust mass and star formation rate in nearby galaxies and
their Lyman-Break galaxy (LBG) analogues at high-$z$ (Burgarella et al. 2020).

###  6.2 Comparison between models and data

We now confront our observational findings to the models described above. To
ensure consistency between observed and simulated data, we impose the same
range of modelled $M_{\star}$, $M_{\rm dust}$ and $\rm sSFR$ as in our
observations ($\log(M_{\star}/M_{\odot})>10$, $\log(M_{\rm dust}/M_{\odot})>6$
and $\log(\rm sSFR/\rm yr^{-1})<-10$). In addition, we follow Damjanov et al.
(2018) and infer the $90\%$ stellar-mass completeness limit of our hCOSMOS
data by fitting the function $M_{\star}=10.6-\log(0.4/z)-0.2$. We find that 46
dusty QGs reside below the completeness limit defined in this way.
Furthermore, only sources identified within the $90\%$ completeness limit are
taken into account when comparing to models. We note that same IMF (Chabrier
2003) is adopted both in simulations, and for the SED fitting of observed
data.

Figure 8: Evolution of $M_{\rm dust}$/$M_{\star}$ as a function of $\rm sSFR$
(left) and $Z_{\rm gas}$ (right) modelled with the state-of-the-art
cosmological simulation SIMBA (Davé et al. 2019, blue lines) and chemical
evolution model of Nanni et al. (2020) (light and dark violet lines). In the
panels, darker (lighter) curves related to N20 model are simulations that
include (exclude) dust grain growth in the ISM, respectively. These simulation
runs include different values of mass-loading factors and the grain
condensation efficiencies, as indicated in the legend. Solid blue line shows
the SIMBA model prediction which is the result of the flagship run (100 Mpc/h
box) that includes ISM dust growth on metals and all feedback variants. In the
right panel, along with the median trend for all selected dusty QGs in SIMBA,
we also display average tracks for QGs with different $\delta_{\rm DGR}$ (blue
dashed line for $\delta_{\rm DGR}<1/100$, and blue dotted-dashed line for
$1/100<\delta_{\rm DGR}<1/2000$). Data from this work (contributing to the
$\sim 90\%$ completeness) are displayed as binned 2D histograms, with colorbar
showing the number of observed QGs in each cell.

In Fig. 8 we show how the $M_{\rm dust}/M_{\star}$ in models evolves as a
function of sSFR and $Z_{\rm gas}$. Predictions from SIMBA match our data very
well as the simulation shows a great potential of simultaneously reproducing
observed trends of $M_{\rm dust}/M_{\star}$ with sSFR and $Z_{\rm gas}$.

For the N20 model we showcase predictions of different simulation runs. These
account for different outflow rates (ML=0.5, 0.55 and 0.6) and different dust
condensation fractions by Type II SNe (low fraction of $25\%$ and high
fraction of $50\%$). The total ISM mass swept in each SNe event is adopted to
be $M_{\rm sw}=6800M_{\odot}$, which is a standard value proposed in the
literature (Dwek et al. 2014). We also explore cases where we switch on/off
modelled dust growth in ISM. We find that the evolution of $M_{\rm
dust}/M_{\star}$ with sSFR in N20 is well reproduced when the dust grain
growth in the ISM is included, and is accompanied with the moderate outflow
and higher dust condensation fraction ($\sim 50\%$). By performing controlled
simulation experiments with N20 we find the evolution of specific dust mass is
extremely sensitive to the changes of ML factor. Namely, the higher outflow
rates ($\rm ML>0.6$) would rapidly reduce $M_{\rm dust}/M_{\star}$ to the
values lower than $10^{-5}$ in older stellar populations ($>5-10$ Gyr). This
would produce swift decline of specific dust mass with sSFR. In addition, for
weaker outflows $M_{\rm dust}/M_{\star}$ remains too flat for decreasing
values of the sSFR, indicating that dust astration and destruction from SNe
are not sufficient to explain the anti-correlation of $M_{\rm dust}/M_{\star}$
with sSFR, even if we impose maximum efficiency of dust destruction. While the
best-matching model of N20 reproduces observed trend of $M_{\rm
dust}/M_{\star}$ with sSFR, at the same time it falls relatively short to
reproduce the relation with $Z_{\rm gas}$. As seen from Fig. 8, $M_{\rm
dust}/M_{\star}$ are 4-5 times too low for the values we infer in super-solar
$Z_{\rm gas}$, while decline in dust abundance is slightly steeper as compared
to SIMBA. At (sub-)solar $Z_{\rm gas}$ both models predict higher $M_{\rm
dust}/M_{\star}$ than inferred by our data. As we note in Section 5.1, the
difference may be just an artefact due to large uncertainties on estimated
$Z_{\rm gas}$. Otherwise, lower $M_{\rm dust}/M_{\star}$ than seen in models
for (sub-)solar $Z_{\rm gas}$ may indicate interesting case of objects missed
by simulations, e.g. galaxies that accrete gaseous material from the cosmic
web, but still did not have enough time to significantly increase their dust
abundance.

###  6.3 Sources of dust in QGs at intermediate redshifts: insights from
models

Confronting the different theoretical models and simulations allows us to
better understand the contributors for dust abundance in QGs. Results from
both observations and models presented throughout this paper strongly support
the scenario where dust growth due to grain collisions in the metal-rich ISM
may play a vital role in dusty QGs. This conclusion is in general agreement
with other observational studies exploring dust content in QGs (Michałowski et
al. 2019, Richtler et al. 2018), and with recent theoretical predictions
(Hirashita et al. 2015, Hirashita & Nozawa 2017, Vijayan et al. 2019 and
Triani et al. 2021). Indeed, it has been demonstrated that excluding ISM grain
growth in SIMBA yields $M_{\rm dust}$ that are $\sim 0.5$ dex lower than in
QGs for which ISM dust growth is included (e.g., Whitaker et al. 2021b).
Similarly, by turning off the effect of ISM grain growth in N20, simulated
$M_{\rm dust}/M_{\star}$ are lower by $\gtrsim 0.4$ dex as compared to cases
when dust growth is included.888Note that various combinations of galactic
outflows and condensation efficiencies nay introduce some degeneracies in the
models, making the predictions strongly dependent on how one regulates the
outflow (see Nanni et al. 2020 for the discussion. We note that the shortfall
to data in N20 cannot be cured by imposing lower ML factors which would result
in an overall flat trend with sSFR, inconsistent with our data (left panel of
Fig. 6). Interestingly, Nanni et al. (2020) found that dust growth has
irrelevant role in the dust content in LBGs and local dwarf galaxies. This
clearly demonstrates that different mechanisms for the dust build-up are in
place in later phases of massive galaxy evolution.

Since the timescale for dust growth in the ISM changes as a function of gas
surface density for given $Z_{\rm gas}$, good agreement of our results with
SIMBA implies that dusty QGs in hCOSMOS have enough available cold gas to
support dust growth over quite short timescales ($\sim 100-250$ Myr). SIMBA
predicts median gas masses of $M_{\rm gas}\sim 2\times 10^{9}\>M_{\odot}$, but
higher $f_{\rm gas}$ in sQGs than in eQGs ($f_{\rm gas}=0.11\pm 0.04$ vs.
$f_{\rm gas}=0.06\pm 0.03$ respectively). In addition, while checking
simulation results, we also notice that sQGs in SIMBA tend to have more
extended distributions of cold gas mass as compared to the stellar mass, which
is not the case in eQGs. This can be explained either with higher ratios of
accreting gas from the outskirts after minor mergers, or that gas is mostly
distributed in the disc regions, while the stellar component is very compact.
Such difference in $f_{\rm gas}$ support our earlier qualitative conclusion
from Section 5.4 and Section 5.5, and may be the at the heart of why observed
$M_{\rm dust}/M_{\star}$ is larger in sQGs than in eQGs.

Quenching timescales in SIMBA are strongly bimodal, with the rapid ones
($<250-500$ Myr) being overly effective within very specific mass range
($10<\log(M_{\star}/M_{\odot})<10.5$) due to black hole feedback mechanism
(Rodríguez Montero et al. 2019, Zheng et al. 2022). The vast majority of dusty
QGs that quench at higher stellar masses ($\log(M_{\star}/M_{\odot})>10.5$)
are expected to be preferentially influenced by the ”slow quenching” mode
(order of $\gtrsim 0.5-2$ Gyr). This mass range corresponds to the $\sim 75\%$
of our sources. Observed shallower decline and large spread of $M_{\rm
dust}/M_{\star}$ towards older ages and higher $M_{\star}$ are consistent with
this scenario. While black hole feedback in SIMBA is responsible for heating
halo gas, it explicitly does not affect most of the cold gas in the ISM since
it is hydrodynamically decoupled. In other words, even if there is no star-
formation, there still would be a material to support the dust growth within
$\sim 10^{7}-10^{8}$ yr after the quenching. Consequently, prolonged growth of
dust would not be possible in the case of X-ray feedback as it has been showed
that in SIMBA it efficiently pushes central gas outwards over very short
timescales (Appleby et al. 2020).

Requirement that ISM growth is needed to simultaneously explain the evolution
of specific dust mass with sSFR and $Z_{\rm gas}$ also supports predictions of
semi-analytical model by Triani et al. (2021). They show that lack of strong
SNe destruction events in dust-rich quenched galaxies supports the super-solar
$Z_{\rm gas}$ and prolonged build-up of dust. While abundance of metals in the
ISM provides material for galaxies to accumulate dust mass, stating that grain
growth is the only process responsible for observed $M_{\rm dust}$ in QGs can
be misleading. In the previous section we offer different lines of evidences
that $\sim 15\%$ of our QGs obey signs of (past) merging activities. Accretion
of a gas-rich satellite would sustain prolonged dust grain growth in the
accreted cold ISM if $Z_{\rm gas}$ is above certain threshold999The exact
value of $Z_{\rm gas}$ where dust growth becomes dominant over stellar
production is proposed to vary. For example, recent analytic models suggest
$\rm 12+\log(O/H)>8.5$ for $M_{\rm gas}>10^{9}M_{\odot}$, e.g., Triani et al.
(2020). We however expect that minor mergers play secondary role in the dust
build up in dusty QGs in SIMBA. Rodríguez Montero et al. (2019) show that
there is a small fraction ($<10\%$) of low-level rejuvenation events mostly
related to minor mergers with gas-rich satellites at $z<1$. These processes
contribute to the young discs and more cold gas in the outskirts, and are
followed by fast ($\sim 100$ Myr) quenching events.

Aside from the two favourable mechanisms outlined above, there are other
possibilities that might be responsible for the dust content observed in our
QGs. Nanni et al. (2020) claims that using a non-standard ”top-heavy” IMF
increases dust yields by a factor of $\times\>3-4$ if a non-standard ’top-
heavy’ IMF is applied (but, see also discussion by McKinnon et al. 2017). At
the same time, we expect a non trivial connection between dusty QGs and their
large-scale environments. The detailed analysis of this interplay is out of
the scope of this paper, but is addressed in our ongoing work (Donevski &
Kraljic, in prep.).

###  6.4 To which extent does the specific dust mass mirror the molecular gas
fraction in QGs?

Finally, relying on a good agreement with SIMBA, we link our observed and
simulated data to inform future observations on expected interplay between the
dust, gas and metals in QGs. In previous sections we show that observed
$M_{\rm dust}/M_{\star}$ can be resulting from a dust growth, for which
timescales may be offsetting from those attributed to the decline of star-
formation and gas exhaustion. We can question: To which extent the observed
evolution of $M_{\rm dust}/M_{\star}$ mirrors the evolution of molecular gas?
To give a sense of how the evolution of $M_{\rm dust}/M_{\star}$ with $Z_{\rm
gas}$ can be used to deduce the cold gas-mass properties, we follow Tan et al.
(2014) and rewrite the $M_{\rm dust}/M_{\star}$ as:

$\frac{M_{\rm dust}}{M_{\star}}\propto\frac{M_{\rm gas}}{M_{\star}}\times
Z_{\rm gas}\times\delta_{\mathrm{DTM}},$ (6)

where $M_{\rm gas}/M_{\star}$ is molecular gas ratio, and
$\delta_{\mathrm{DTM}}$ is dust-to-metal ratio.

Figure 9: Upper panel: Modelled $M_{\rm dust}/M_{\star}$ as a function of
$Z_{\rm gas}$ in SIMBA. Circles represent the modelled QGs within the
restricted range of dust-to-gas mass ratios
($1/20<\delta_{\mathrm{DGR}}<1/100$). Grey contours reveal the distribution of
all QGs selected in the simulations with the criteria described in Section
6.1. The points are colour-coded for the $\delta_{\rm{DTM}}$, with circle
sizes scaling to $\rm f_{\rm gas}$. Inset sketch qualitatively describes the
behaviour of $M_{\rm dust}/M_{\star}$: up to the certain $Z_{\rm gas}$ it
mostly resembles change of molecular gas fraction, while at $\rm
12+\log(O/H)\gtrsim 8.9$ it reflects the level of dust growth through the
gradual increase in $\delta_{\rm{DTM}}$; bbcbbcLower panel: Evolution of
$M_{\rm dust}/M_{\star}$ with $f_{\rm gas}$ for the same simulated sources
plotted in the upper panel. The blue and orange dashed lines show median
trends for two regimes of $Z_{\rm gas}$ as indicated in the legend. The effect
of increased dust growth in QGs is most pronounced in sources with heavily
exhausted gas reservoirs, as indicated by lower molecular gas fractions
($f_{\rm gas}<1-5\%$). To highlight this effect through the rise in
$\delta_{\rm{DTM}}$, we sketch the shaded area between the two tracks and
colour-code it to roughly match the colorbar convention from the upper panel.

There are two possible approaches to solve this equation. One is to follow the
method applied by Magdis et al. (2021), who assumed constant, metallicity
dependent $\delta_{\mathrm{DGR}}$ for all QGs (they adopt
$\delta_{\mathrm{DGR}}=1/92(1/35))$ for solar and super-solar $Z_{\rm gas}$,
respectively). However, this solution implies a constant dust-to-metal ratio,
in which case the evolution of $M_{\rm dust}/M_{\star}$ entirely mirrors the
evolution of molecular gas fraction. While this offers an appealing empirical
explanation for the observed evolution of the ISM properties of QGs, it does
not account for the possible variations in dust-to-metal ratio. Indeed,
$\delta_{\mathrm{DTM}}$ would remain constant if dust and metals are created
from stars at same rates throughout galaxy evolution. Otherwise, with metals
supporting the dust grain growth, $\delta_{\mathrm{DTM}}$ would change and
becomes important marker of the growth efficiency (i.e., Feldmann 2015, De Vis
et al. 2019).

We thus apply another approach which benefits from the good agreement with
SIMBA. We explore what lies behind the scatter in $M_{\rm gas}/M_{\star}$ at a
given $Z_{\rm gas}$ in simulations. To illustrate how $f_{\rm gas}$ and
$\delta_{\mathrm{DTM}}$ relate to $M_{\rm dust}/M_{\star}$, in Fig.9 we
display simulated QGs that have $\delta_{\mathrm{DGR}}$ within the restricted
range of values that are usually adopted in observations to convert observed
$M_{\rm dust}$ to $M_{\rm gas}$ ($1/20<\delta_{\mathrm{DGR}}<1/100$; e.g.,
Leroy et al. 2011, Schreiber et al. 2018, Dudzevičiūtė et al. 2020, Whitaker
et al. 2021a, Magdis et al. 2021). From the upper panel of Fig.9 we see that
even for such a relatively small difference in $\delta_{\mathrm{DGR}}$,
specific dust masses display wider spread in values ($\sim 8\times
10^{-5}<M_{\rm dust}/M_{\star}<3\times 10^{-3}$).

Most importantly, Fig.9 reveals that gas fraction is not the faithful tracer
of $M_{\rm dust}/M_{\star}$ over the entire range of $Z_{\rm gas}$. SIMBA
predicts that up to $\rm 12+\log(O/H)\lesssim 8.85-8.9$, observed $M_{\rm
dust}/M_{\star}$ mainly traces $f_{\rm gas}$, but this correlation is weaker
for $12+\log(\rm O/H)\gtrsim 8.85-8.9$, when increase in
$\delta_{\mathrm{DTM}}$ becomes more important. In other words, it is gradual
evolution of $\delta_{\mathrm{DTM}}$, rather than of gas fraction, that are
responsible for the spread in observed specific dust mass in super-solar dusty
QGs. Therefore, we should expect that some of very metal-rich QGs with
$\delta_{\mathrm{DTM}}>0.5$ may still have relatively high specific dust mass
($10^{-3}<M_{\rm dust}/M_{\star}<10^{-4}$) despite heavily exhausted gas
reservoirs ($f_{\rm gas}<1-5\%$). This is illustrated in the lower panel of
Fig.9. The effect of prolonged dust growth on metals causes different pathways
(slower and faster) for the decline in $M_{\rm dust}/M_{\star}$. The
corresponding offset in $M_{\rm dust}/M_{\star}$ starts to increase below
certain $f_{\rm gas}$. In galaxies that show positive excess, thermal
sputtering acts on longer timescales ($>1$ Gyr) due to inefficiency in swiftly
destroying larger grains produced via grain growth. If this is not the case,
QGs would end up with extremely low $\delta_{\mathrm{DGR}}$ as reported in
some dusty QGs observed in the local Universe (e.g., Lianou et al. 2016,
Casasola et al. 2020). For that reason, any use of $M_{\rm dust}/M_{\star}$ to
convert to $f_{\rm gas}$ and ultimately to $M_{\rm gas}$ must be carefully
treated when QGs have super-solar $Z_{\rm gas}$. This is particularly
important for planning the observations (i.e., with ALMA) to detect the
molecular gas via CO[2-1]. At a given $M_{\rm dust}/M_{\star}$, galaxies with
the highest $Z_{\rm gas}$ would require significantly longer integration time
than those with $Z_{\rm gas}\sim Z_{\odot}$ due to their higher
$\delta_{\mathrm{DGR}}$ and subsequently lower $M_{\rm gas}$. This suggests
that care should be taken when interpreting the nature of dust evolution in
metal-rich QGs.

##  7 Conclusions

We present the first statistical study of the evolution of dust-related
properties in QGs observed at intermediate redshifts ($0.1<z<0.6$) as part of
the hCOSMOS spectroscopic survey. We analyse 548 massive
($M_{\star}>10^{10}\>M_{\odot}$) QGs that are systematically selected based on
their spectral age indices ($D_{\rm n}4000>1.5$). We combine spectroscopic
information with a self-consistent, multi-band SED fitting method to explore
trends of $M_{\rm dust}/M_{\star}$ with various physical parameters, including
independently measured gas-phase metallicity. This is the largest sample of
QGs for which the evolution of dust, metal and stellar content have been
estimated to date. We fully evaluate our findings with state-of-the-art
simulations of dusty galaxy formation. Our main results are summarised as
follows:

  * •

Dusty QGs at intermediate redshifts experience the cosmic evolution in $M_{\rm
dust}/M_{\star}$, but we find it to be shallower than in known studies on
stacked samples. A large spread ($>2$ orders of magnitude) in $M_{\rm
dust}$/$M_{\rm\star}$ suggests non-uniform ISM conditions in these massive
galaxies.

  * •

Morphological and structural parameters are important factors in the scatter
in specific dust masses. The $M_{\rm dust}/M_{\rm\star}$ in sQGs is $\gtrsim
4$ times greater than in eQGs of a similar $M_{\star}$. At
$M_{\star}<10^{11}M_{\odot}$, higher specific dust masses are associated with
lower stellar mass surface densities and younger ages of stellar populations,
whereas the same relationship is less obvious at $M_{\star}>10^{11}M_{\odot}$.

  * •

The sQGs exhibit evolutionary trends of $M_{\rm dust}/M_{\rm\star}$ with
$M_{\star}$, stellar age and galaxy size, in contrast to little to no
evolution in eQGs.

  * •

Although both sQGs and eQGs have median super-solar $Z_{\rm gas}$
($12+\log(\rm O/H)\sim 9$), sQGs have a reversal of $M_{\rm
dust}/M_{\rm\star}$ towards the highest $Z_{\rm gas}$, whereas eQGs have a
flat trend. We interpret this as a change in efficiency of dust production and
removal since quenching.

  * •

We derive a broad dynamical range of post-quenching timescales in QGs
($60\>\rm Myr<t_{\rm quench}<3.2\>\rm Gyr$). In general, $M_{\rm
dust}/M_{\rm\star}$ is the highest in recently quenched systems ($t_{\rm
quench}<500$ Myr), but its further evolution is non-monotonic, which is strong
evidence for different pathways for prolonged dust growth, or removal on
various timescales. The moderately shallow evolution of $M_{\rm
dust}/M_{\rm\star}$ imply that quenching mechanisms usually attributed to
quick removal of gas and dust (i.e., powerful outflows and X-ray feedback) are
not overly influential within our sample.

  * •

State-of-the art cosmological simulation SIMBA and the chemical model of Nanni
et al. (2020) are both capable of producing the dustiest QGs observed in this
work ($M_{\rm dust}/M_{\star}\gtrsim 10^{-3}$), with the difference that
chemical models fall short of fully explaining the observed specific dust
masses at super-solar $Z_{\rm gas}$. The prolonged grain growth on ISM metals
is strongly required in simulations to account for the observed dust content
in QGs. Without this assumption, the models have a larger discrepancy with the
data, even when adopting substantial stellar yields with high condensation
efficiencies and moderate outflows.

  * •

Our results strongly suggest that the observed $M_{\rm dust}$/$M_{\rm\star}$
in dusty QGs is twofold: (1) The dominant channel is prolonged dust growth in
metal-rich ISM. Such a process is viable in the first 1 Gyr since quenching,
and it helps dusty QGs significantly prolong the destruction timescales
expected for gas heating; (2) For $\sim 15\%$ of sources we find evidence that
the dust content is supported or acquired externally, most likely via minor
mergers.

  * •

Cosmological simulation SIMBA predicts that enhanced $M_{\rm
dust}/M_{\rm\star}$ in many QGs with super-solar $Z_{\rm gas}$ mirrors the
enhancement in dust-to-metal ratio rather than rise in molecular gas fraction.
This prediction has an important observational consequence in that QGs can
still be observed as dust-rich ($M_{\rm dust}/M_{\rm\star}\simeq 10^{-3}$)
even if their gas reservoirs are heavily exhausted (with the gas fraction
$<1-5\%$).

The complexity of galaxies’ dust life-cycles revealed by our study motivates
future observational investigations of dusty QGs. The model predictions
described in Section 6. (in particular the diversity in $\delta_{\rm{DTM}}$)
are rooted in the mechanisms dominating dust production and removal. Many
observational tests of these predictions are now directly accessible via JWST
instruments. The JWST MIRI, in particular, can be used for probing the warm
dust and molecular gas through mid-IR diagnostics in individual dusty QGs at
various cosmic epochs up to $z\sim 2-3$. Along with independent metallicity
measurements, the JWST MIRI can be used to place strong constraints on the
fractions and sizes of PAH grains through multiple PAH features with a large
wavelength separation. At the same time, it can help unveiling the impact of
dust-hidden AGN activity that would influence the dust removal. Furthermore,
silicate strengths and a rising mid-IR emission longwards of $\geq 10\mu$m
would offer insights on dust grain growth. There are numerous synergistic
opportunities with high-resolution sub-mm observations (e.g., with ALMA and
NOEMA) that would provide information on the cold dust and gas in QGs.
Finally, in our ongoing works, we examine how the large-scale environments
influence the interplay of dust and metals in QGs (Donevski & Kraljic, in
prep.) and present the key differences in SEDs of dusty vs. non-dusty QGs
(Lorenzon et al., in prep.).

###### Acknowledgements.

We are thankful to the anonymous referee for very constructive comments which
improved the initial version of the paper. We would like to acknowledge Anna
Feltre, Kasia Malek, Marcella Massardi, Lorenzo Zanisi, Katarina Kraljic and
Alessandro Bressan for useful discussions and comments D.D acknowledges
support from the National Science Center (NCN) grant SONATA
(UMO-2020/39/D/ST9/00720). I.D. acknowledges the support of the Canada
Research Chair Program and the Natural Sciences and Engineering Research
Council of Canada (NSERC; funding reference No. RGPIN-2018-05425). A.M.
acknowledges the support of the Natural Sciences and Engineering Research
Council of Canada (NSERC) through grant reference number RGPIN-2021-03046.
M.R. and A.N. acknowledge support from the Narodowe Centrum Nauki
(UMO-2020/38/E/ST9/00077). D.N. acknowledges support from the NSF via
AST-1909153.

## References

  * Andreani et al. (2018) Andreani, P., Boselli, A., Ciesla, L., et al. 2018, A&A, 617, A33
  * Aoyama et al. (2019) Aoyama, S., Hirashita, H., & Nagamine, K. 2019, MNRAS, 2834
  * Appleby et al. (2020) Appleby, S., Davé, R., Kraljic, K., Anglés-Alcázar, D., & Narayanan, D. 2020, MNRAS, 494, 6053
  * Asano et al. (2013) Asano, R. S., Takeuchi, T. T., Hirashita, H., & Inoue, A. K. 2013, Earth, Planets, and Space, 65, 213
  * Barone et al. (2022) Barone, T. M., D’Eugenio, F., Scott, N., et al. 2022, MNRAS, 512, 3828
  * Battisti et al. (2020) Battisti, A. J., Cunha, E. d., Shivaei, I., & Calzetti, D. 2020, ApJ, 888, 108
  * Beeston et al. (2018) Beeston, R. A., Wright, A. H., Maddox, S., et al. 2018, MNRAS, 479, 1077
  * Belli et al. (2021) Belli, S., Contursi, A., Genzel, R., et al. 2021, ApJ, 909, L11
  * Berta et al. (2016) Berta, S., Lutz, D., Genzel, R., Förster-Schreiber, N. M., & Tacconi, L. J. 2016, A&A, 587, A73
  * Béthermin et al. (2015) Béthermin, M., Daddi, E., Magdis, G., et al. 2015, A&A, 573, A113
  * Beverage et al. (2021) Beverage, A. G., Kriek, M., Conroy, C., et al. 2021, ApJ, 917, L1
  * Bianchi & Schneider (2007) Bianchi, S. & Schneider, R. 2007, MNRAS, 378, 973
  * Bisigello et al. (2020) Bisigello, L., Kuchner, U., Conselice, C. J., et al. 2020, MNRAS, 494, 2337
  * 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
  * Boquien et al. (2019) Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103
  * Bourne et al. (2012) Bourne, N., Maddox, S. J., Dunne, L., et al. 2012, MNRAS, 421, 3027
  * Brinchmann et al. (2004) Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151
  * Bruzual & Charlot (2003) Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000
  * Burgarella et al. (2020) Burgarella, D., Nanni, A., Hirashita, H., et al. 2020, A&A, 637, A32
  * Bustamante et al. (2018) Bustamante, S., Sparre, M., Springel, V., & Grand, R. J. J. 2018, MNRAS, 479, 3381
  * Calura et al. (2017) Calura, F., Pozzi, F., Cresci, G., et al. 2017, MNRAS, 465, 54
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245
  * Casasola et al. (2020) Casasola, V., Bianchi, S., De Vis, P., et al. 2020, A&A, 633, A100
  * Casasola et al. (2022) Casasola, V., Bianchi, S., Magrini, L., et al. 2022, arXiv e-prints, arXiv:2210.15993
  * Cassata et al. (2007) Cassata, P., Guzzo, L., Franceschini, A., et al. 2007, ApJS, 172, 270
  * Chabrier (2003) Chabrier, G. 2003, Publications of the Astronomical Society of the Pacific, 115, 763
  * Charlot & Fall (2000) Charlot, S. & Fall, S. M. 2000, ApJ, 539, 718
  * Chen et al. (2020) Chen, Z., Faber, S. M., Koo, D. C., et al. 2020, ApJ, 897, 102
  * Ciesla et al. (2016) Ciesla, L., Boselli, A., Elbaz, D., et al. 2016, A&A, 585, A43
  * Ciesla et al. (2021) Ciesla, L., Buat, V., Boquien, M., et al. 2021, A&A, 653, A6
  * Civano et al. (2016) Civano, F., Marchesi, S., Comastri, A., et al. 2016, ApJ, 819, 62
  * Clemens et al. (2010) Clemens, M. S., Jones, A. P., Bressan, A., et al. 2010, A&A, 518, L50
  * Cousin et al. (2019) Cousin, M., Buat, V., Lagache, G., & Bethermin, M. 2019, A&A, 627, A132
  * Cristallo et al. (2015) Cristallo, S., Straniero, O., Piersanti, L., & Gobrecht, D. 2015, ApJS, 219, 40
  * Cuppen et al. (2017) Cuppen, H. M., Walsh, C., Lamberts, T., et al. 2017, Space Sci. Rev., 212, 1
  * da Cunha et al. (2008) da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595
  * da Cunha et al. (2010) da Cunha, E., Eminian, C., Charlot, S., & Blaizot, J. 2010, MNRAS, 403, 1894
  * da Cunha et al. (2015) da Cunha, E., Walter, F., Smail, I. R., et al. 2015, ApJ, 806, 110
  * Daddi et al. (2005) Daddi, E., Renzini, A., Pirzkal, N., et al. 2005, ApJ, 626, 680
  * Damjanov et al. (2022) Damjanov, I., Sohn, J., Geller, M. J., Utsumi, Y., & Dell’Antonio, I. 2022, arXiv e-prints, arXiv:2210.01129
  * Damjanov et al. (2018) Damjanov, I., Zahid, H. J., Geller, M. J., Fabricant, D. G., & Hwang, H. S. 2018, ApJS, 234, 21
  * Dariush et al. (2016) Dariush, A., Dib, S., Hony, S., et al. 2016, MNRAS, 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. (2013) Davis, T. A., Alatalo, K., Bureau, M., et al. 2013, MNRAS, 429, 534
  * De Looze et al. (2016) De Looze, I., Baes, M., Bendo, G. J., et al. 2016, MNRAS, 459, 3900
  * De Vis et al. (2017) De Vis, P., Gomez, H. L., Schofield, S. P., et al. 2017, MNRAS, 471, 1743
  * De Vis et al. (2019) De Vis, P., Jones, A., Viaene, S., et al. 2019, A&A, 623, A5
  * Donevski et al. (2020) Donevski, D., Lapi, A., Małek, K., et al. 2020, A&A, 644, A144
  * Draine & Li (2007) Draine, B. T. & Li, A. 2007, ApJ, 657, 810
  * Driver et al. (2018) Driver, S. P., Andrews, S. K., da Cunha, E., et al. 2018, MNRAS, 475, 2891
  * Dudzevičiūtė et al. (2020) Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2020, MNRAS, 494, 3828
  * Dwek (1998) Dwek, E. 1998, ApJ, 501, 643
  * Dwek et al. (2014) Dwek, E., Staguhn, J., Arendt, R. G., et al. 2014, ApJ, 788, L30
  * Fabricant et al. (2005) Fabricant, D., Fata, R., Roll, J., et al. 2005, PASP, 117, 1411
  * Feldmann (2015) Feldmann, R. 2015, MNRAS, 449, 3274
  * Feltre et al. (2012) Feltre, A., Hatziminaoglou, E., Fritz, J., & Franceschini, A. 2012, MNRAS, 426, 120
  * Ferrarotti & Gail (2006) Ferrarotti, A. S. & Gail, H. P. 2006, A&A, 447, 553
  * Fraser-McKelvie et al. (2018) Fraser-McKelvie, A., Brown, M. J. I., Pimbblet, K., Dolley, T., & Bonne, N. J. 2018, MNRAS, 474, 1909
  * Fritz et al. (2006) Fritz, J., Franceschini, A., & Hatziminaoglou, E. 2006, MNRAS, 366, 767
  * Gall & Hjorth (2018) Gall, C. & Hjorth, J. 2018, ApJ, 868, 62
  * Gallazzi et al. (2005) Gallazzi, A., Charlot, S., Brinchmann, J., White, S. D. M., & Tremonti, C. A. 2005, MNRAS, 362, 41
  * Galliano et al. (2018) Galliano, F., Galametz, M., & Jones, A. P. 2018, ARA&A, 56, 673
  * Galliano et al. (2021) Galliano, F., Nersesian, A., Bianchi, S., et al. 2021, A&A, 649, A18
  * Geach et al. (2017) Geach, J. E., Dunlop, J. S., Halpern, M., et al. 2017, MNRAS, 465, 1789
  * Gobat et al. (2018) Gobat, R., Daddi, E., Magdis, G., et al. 2018, Nature Astronomy, 2, 239
  * Gobat et al. (2020) Gobat, R., Magdis, G., D’Eugenio, C., & Valentino, F. 2020, A&A, 644, L7
  * Graziani et al. (2020) Graziani, L., Schneider, R., Ginolfi, M., et al. 2020, MNRAS, 494, 1071
  * Griffith et al. (2019) Griffith, E., Martini, P., & Conroy, C. 2019, MNRAS, 484, 562
  * Grishin et al. (2021) Grishin, K. A., Chilingarian, I. V., Afanasiev, A. V., et al. 2021, Nature Astronomy, 5, 1308
  * Haines et al. (2017) Haines, C. P., Iovino, A., Krywult, J., et al. 2017, A&A, 605, A4
  * Hirashita & Nozawa (2017) Hirashita, H. & Nozawa, T. 2017, Planet. Space Sci., 149, 45
  * Hirashita et al. (2015) Hirashita, H., Nozawa, T., Villaume, A., & Srinivasan, S. 2015, MNRAS, 454, 1620
  * Hjorth et al. (2014) Hjorth, J., Gall, C., & Michałowski, M. J. 2014, ApJ, 782, L23
  * Hou et al. (2019) Hou, K.-C., Aoyama, S., Hirashita, H., Nagamine, K., & Shimizu, I. 2019, MNRAS, 485, 1727
  * Hunt et al. (2019) Hunt, L. K., De Looze, I., Boquien, M., et al. 2019, A&A, 621, A51
  * Hunt et al. (2014) Hunt, L. K., Palazzi, E., Michałowski, M. J., et al. 2014, A&A, 565, A112
  * Hurley et al. (2017) Hurley, P. D., Oliver, S., Betancourt, M., et al. 2017, MNRAS, 464, 885
  * Imara et al. (2018) Imara, N., Loeb, A., Johnson, B. D., Conroy, C., & Behroozi, P. 2018, ApJ, 854, 36
  * Iwamoto et al. (1999) Iwamoto, K., Brachwitz, F., Nomoto, K., et al. 1999, ApJS, 125, 439
  * Jones et al. (2017) Jones, A. P., Köhler, M., Ysard, N., Bocchio, M., & Verstraete, L. 2017, A&A, 602, A46
  * Kartaltepe et al. (2010) Kartaltepe, J. S., Sanders, D. B., Le Floc’h, E., et al. 2010, ApJ, 721, 98
  * Kennicutt (1998) Kennicutt, Jr., R. C. 1998, ARA&A, 36, 189
  * Kirkpatrick et al. (2017) Kirkpatrick, A., Pope, A., Sajina, A., et al. 2017, ApJ, 843, 71
  * Kobayashi et al. (2006) Kobayashi, C., Umeda, H., Nomoto, K., Tominaga, N., & Ohkubo, T. 2006, ApJ, 653, 1145
  * Kobulnicky & Kewley (2004) Kobulnicky, H. A. & Kewley, L. J. 2004, ApJ, 617, 240
  * Kokusho et al. (2019) Kokusho, T., Kaneda, H., Bureau, M., et al. 2019, A&A, 622, A87
  * Kumari et al. (2021) Kumari, N., Maiolino, R., Trussler, J., et al. 2021, A&A, 656, A140
  * Lacey et al. (2016) Lacey, C. G., Baugh, C. M., Frenk, C. S., et al. 2016, MNRAS, 462, 3854
  * Lagos et al. (2019) Lagos, C. d. P., Robotham, A. S. G., Trayford, J. W., et al. 2019, MNRAS, 489, 4196
  * Lapi et al. (2018) Lapi, A., Pantoni, L., Zanisi, L., et al. 2018, ApJ, 857, 22
  * Leethochawalit et al. (2019) Leethochawalit, N., Kirby, E. N., Ellis, R. S., Moran, S. M., & Treu, T. 2019, ApJ, 885, 100
  * Leja et al. (2019) Leja, J., Tacchella, S., & Conroy, C. 2019, ApJ, 880, L9
  * Leroy et al. (2011) Leroy, A. K., Bolatto, A., Gordon, K., et al. 2011, ApJ, 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. (2019) Lianou, S., Barmby, P., Mosenkov, A. A., Lehnert, M., & Karczewski, O. 2019, A&A, 631, A38
  * Lianou et al. (2016) Lianou, S., Xilouris, E., Madden, S. C., & Barmby, P. 2016, MNRAS, 461, 2856
  * Liu et al. (2019) Liu, D., Schinnerer, E., Groves, B., et al. 2019, ApJ, 887, 235
  * Lovell et al. (2021) Lovell, C. C., Geach, J. E., Davé, R., Narayanan, D., & Li, Q. 2021, MNRAS, 502, 772
  * Magdis et al. (2021) Magdis, G. E., Gobat, R., Valentino, F., et al. 2021, A&A, 647, A33
  * Mahajan et al. (2020) Mahajan, S., Gupta, K. K., Rana, R., et al. 2020, MNRAS, 491, 398
  * Małek et al. (2018) Małek, K., Buat, V., Roehlly, Y., et al. 2018, A&A, 620, A50
  * Man et al. (2016) Man, A. W. S., Greve, T. R., Toft, S., et al. 2016, ApJ, 820, 11
  * Martini et al. (2013) Martini, P., Dicken, D., & Storchi-Bergmann, T. 2013, ApJ, 766, 121
  * Martis et al. (2019) Martis, N. S., Marchesini, D. M., Muzzin, A., et al. 2019, ApJ, 882, 65
  * McKinnon et al. (2017) McKinnon, R., Torrey, P., Vogelsberger, M., Hayward, C. C., & Marinacci, F. 2017, MNRAS, 468, 1505
  * Michałowski et al. (2019) Michałowski, M. J., Hjorth, J., Gall, C., et al. 2019, A&A, 632, A43
  * Moresco et al. (2013) Moresco, M., Pozzetti, L., Cimatti, A., et al. 2013, A&A, 558, A61
  * Muzzin et al. (2013) Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18
  * Nanni et al. (2020) Nanni, A., Burgarella, D., Theulé, P., Côté, B., & Hirashita, H. 2020, A&A, 641, A168
  * Nersesian et al. (2019) Nersesian, A., Xilouris, E. M., Bianchi, S., et al. 2019, A&A, 624, A80
  * Nevin et al. (2019) Nevin, R., Blecha, L., Comerford, J., & Greene, J. 2019, ApJ, 872, 76
  * Nishida et al. (2022) Nishida, K. Y., Takeuchi, T. T., Nagata, T., & Asano, R. S. 2022, MNRAS, 514, 2098
  * Noll et al. (2009) Noll, S., Burgarella, D., Giovannoli, E., et al. 2009, A&A, 507, 1793
  * Osterbrock & Ferland (2006) Osterbrock, D. E. & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei
  * Pantoni et al. (2021) Pantoni, L., Lapi, A., Massardi, M., et al. 2021, MNRAS, 504, 928
  * Pantoni et al. (2019) Pantoni, L., Lapi, A., Massardi, M., Goswami, S., & Danese, L. 2019, ApJ, 880, 129
  * Park et al. (2022) Park, M., Tacchella, S., Nelson, E. J., et al. 2022, MNRAS, 515, 213
  * Pearson et al. (2018) Pearson, W. J., Wang, L., Hurley, P. D., et al. 2018, A&A, 615, A146
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13
  * 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. J., et al. 2022, MNRAS, 516, 2314
  * Relaño et al. (2022) Relaño, M., De Looze, I., Saintonge, A., et al. 2022, MNRAS, 515, 5306
  * Rémy-Ruyer et al. (2014) Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, A31
  * Richtler et al. (2018) Richtler, T., Hilker, M., Voggel, K., et al. 2018, A&A, 620, A133
  * Rodríguez Montero et al. (2019) Rodríguez Montero, F., Davé, R., Wild, V., Anglés-Alcázar, D., & Narayanan, D. 2019, MNRAS, 490, 2139
  * Rowlands et al. (2014) Rowlands, K., Dunne, L., Dye, S., et al. 2014, MNRAS, 441, 1017
  * Rowlands et al. (2012) Rowlands, K., Dunne, L., Maddox, S., et al. 2012, MNRAS, 419, 2545
  * Rudnick et al. (2017) Rudnick, G., Hodge, J., Walter, F., et al. 2017, ApJ, 849, 27
  * Sansom et al. (2019) Sansom, A. E., Glass, D. H. W., Bendo, G. J., et al. 2019, MNRAS, 482, 4617
  * Santini et al. (2014) Santini, P., Maiolino, R., Magnelli, B., et al. 2014, A&A, 562, A30
  * Sargent et al. (2014) Sargent, M. T., Daddi, E., Béthermin, M., et al. 2014, ApJ, 793, 19
  * Sargent et al. (2015) Sargent, M. T., Daddi, E., Bournaud, F., et al. 2015, ApJ, 806, L20
  * Scarlata et al. (2007) Scarlata, C., Carollo, C. M., Lilly, S., et al. 2007, ApJS, 172, 406
  * Schmidt (1959) Schmidt, M. 1959, ApJ, 129, 243
  * Schneider et al. (2014) Schneider, R., Valiante, R., Ventura, P., et al. 2014, MNRAS, 442, 1440
  * Schreiber et al. (2018) Schreiber, C., Elbaz, D., Pannella, M., et al. 2018, A&A, 609, A30
  * Schreiber et al. (2015) Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74
  * Scoville et al. (2017) Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150
  * Scoville et al. (2016) Scoville, N., Sheth, K., Aussel, H., et al. 2016, ApJ, 820, 83
  * Shirley et al. (2019) Shirley, R., Roehlly, Y., Hurley, P. D., et al. 2019, MNRAS, 490, 634
  * Shivaei et al. (2022) Shivaei, I., Popping, G., Rieke, G., et al. 2022, ApJ, 928, 68
  * Siudek et al. (2018) Siudek, M., Małek, K., Pollo, A., et al. 2018, A&A, 617, A70
  * Slavin et al. (2015) Slavin, J. D., Dwek, E., & Jones, A. P. 2015, ApJ, 803, 7
  * Smercina et al. (2018) Smercina, A., Smith, J. D. T., Dale, D. A., et al. 2018, ApJ, 855, 51
  * Smith et al. (2014) Smith, D. J. B., Jarvis, M. J., Hardcastle, M. J., et al. 2014, MNRAS, 445, 2232
  * Sohn et al. (2019) Sohn, J., Geller, M. J., Zahid, H. J., & Fabricant, D. G. 2019, ApJ, 872, 192
  * Speagle et al. (2014) Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15
  * Strandet et al. (2017) Strandet, M. L., Weiss, A., De Breuck, C., et al. 2017, ApJ, 842, L15
  * Suess et al. (2022) Suess, K. A., Leja, J., Johnson, B. D., et al. 2022, ApJ, 935, 146
  * Tan et al. (2014) Tan, Q., Daddi, E., Magdis, G., et al. 2014, A&A, 569, A98
  * Tasca et al. (2009) Tasca, L. A. M., Kneib, J. P., Iovino, A., et al. 2009, A&A, 503, 379
  * Toft et al. (2007) Toft, S., van Dokkum, P., Franx, M., et al. 2007, ApJ, 671, 285
  * Tojeiro et al. (2013) Tojeiro, R., Masters, K. L., Richards, J., et al. 2013, MNRAS, 432, 359
  * Triani et al. (2021) Triani, D. P., Sinha, M., Croton, D. J., Dwek, E., & Pacifici, C. 2021, MNRAS, 503, 1005
  * Triani et al. (2020) Triani, D. P., Sinha, M., Croton, D. J., Pacifici, C., & Dwek, E. 2020, MNRAS
  * Utsumi et al. (2020) Utsumi, Y., Geller, M. J., Zahid, H. J., et al. 2020, ApJ, 900, 50
  * Valiante et al. (2011) Valiante, R., Schneider, R., Salvadori, S., & Bianchi, S. 2011, MNRAS, 416, 1916
  * van Dokkum et al. (2010) van Dokkum, P. G., Whitaker, K. E., Brammer, G., et al. 2010, ApJ, 709, 1018
  * Ventura et al. (2020) Ventura, P., Dell’Agli, F., Lugaro, M., et al. 2020, A&A, 641, A103
  * Vergani et al. (2018) Vergani, D., Garilli, B., Polletta, M., et al. 2018, A&A, 620, A193
  * Vijayan et al. (2019) Vijayan, A. P., Clay, S. J., Thomas, P. A., et al. 2019, MNRAS, 489, 4072
  * Vílchez et al. (2019) Vílchez, J. M., Relaño, M., Kennicutt, R., et al. 2019, MNRAS, 483, 4968
  * Vogelsberger et al. (2019) Vogelsberger, M., McKinnon, R., O’Neil, S., et al. 2019, MNRAS, 487, 4870
  * Whitaker et al. (2021a) Whitaker, K. E., Narayanan, D., Williams, C. C., et al. 2021a, ApJ, 922, L30
  * 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
  * Woodrum et al. (2022) Woodrum, C., Williams, C. C., Rieke, M., et al. 2022, arXiv e-prints, arXiv:2210.03832
  * Wu et al. (2021) Wu, P.-F., Nelson, D., van der Wel, A., et al. 2021, AJ, 162, 201
  * Wu et al. (2020) Wu, P.-F., van der Wel, A., Bezanson, R., et al. 2020, ApJ, 888, 77
  * Wu et al. (2018) Wu, P.-F., van der Wel, A., Bezanson, R., et al. 2018, ApJ, 868, 37
  * Zahid et al. (2016) Zahid, H. J., Baeza Hochmuth, N., Geller, M. J., et al. 2016, ApJ, 831, 146