High dust content of a quiescent galaxy at z~2 revealed by deep ALMA observation

We report the detection of cold dust in an apparently quiescent massive galaxy ($\log({M_{\star}/M_{\odot}})\approx11$) at $z\sim2$ (G4). The source is identified as a serendipitous 2 mm continuum source in a deep ALMA observation within the field of Q2343-BX610, a $z=2.21$ massive star-forming disk galaxy. Available multi-band photometry of G4 suggests redshift of $z\sim2$ and a low specific star-formation rate (sSFR), $\log(SFR/M_{\star}) [yr^{-1}] \approx -10.2$, corresponding to $\approx1.2$ dex below the $z=2$ main sequence (MS). G4 appears to be a peculiar dust-rich quiescent galaxy for its stellar mass ($\log({M_{\rm dust}/M_{\star}}) = -2.71 \pm 0.26$), with its estimated mass-weighted age ($\sim$ 1-2 Gyr). We compile $z\gtrsim1$ quiescent galaxies in the literature and discuss their age-$\Delta$MS and $\log({M_{\rm dust}/M_{\star}})$-age relations to investigate passive evolution and dust depletion scale. A long dust depletion time and its morphology suggest morphological quenching along with less efficient feedback that could have acted on G4. The estimated dust yield for G4 further supports this idea, requiring efficient survival of dust and/or grain growth, and rejuvenation (or additional accretion). Follow-up observations probing the stellar light and cold dust peak are necessary to understand the implication of these findings in the broader context of galaxy evolutionary studies and quenching in the early universe.


INTRODUCTION
The general picture drawn over the last two decades in galaxy evolutionary study is that galaxies show bimodality in their colors, morphologies, and star formation rates (SFRs).Star-forming galaxies are characterized by bluer colors and disk-like morphology, having rich cold interstellar medium (ISM).In comparison, quiescent galaxies have redder colors, spheroidal morphology (or higher fraction of bulge component) with less gas (e.g., Strateva et al. 2001;Kauffmann et al. 2003;Baldry et al. 2004;Balogh et al. 2004;Bell et al. 2004;Wuyts et al. 2011;Saintonge et al. 2011;Young et al. 2011).In line with this picture, extensive surveys of cold ISM in star-forming galaxies up to  ∼ 2 have shown that the cosmic star-formation activity (e.g., Madau & Dickinson 2014) is primarily governed by the available gas within (and outside that falls into) the galaxies (e.g., Tacconi et al. 2020;Walter et al. 2020 and references therein).
Therefore, it has been generally considered that a key aspect of quenching star formation is the reduction of cold ISM reservoirs.
Theories and simulations proposed that it can be achieved by expulsion via outflows or depletion via star formation (e.g., Silk & Rees 1998;Di Matteo et al. 2005;Springel et al. 2005).Active galactic nuclei (AGN) have been regarded as a preferable agent of outflows in the high mass (log( ★ / ⊙ ) ≳ 10.7) regime, which is also observationally indicated at  ∼ 1 − 2 (e.g., Förster Schreiber et al. 2018 and references therein;Concas et al. 2022).Supernova feedback (outflow and shock) would be more efficient in sweeping the gas out of the galaxy in a low mass regime (e.g., Dekel & Silk 1986;Efstathiou 2000).These two mechanisms were taken to explain the mismatch between the observed stellar mass function and halo mass function high-mass (AGN) and low-mass (supernovae) regimes early on in galaxy evolution studies (e.g., White & Rees 1978;White & Frenk 1991;Balogh et al. 2001;Springel & Hernquist 2003; see also reviews of Somerville & Davé 2015;Naab & Ostriker 2017).A remarkable consistency in the cosmic evolution of the black hole accretion activity further implies a possible role of AGN in the interplay of cosmic gas and star formation evolution.Additional mechanisms that may contribute to quenching without exhausting the available gas reservoir include heating of the halo gas (via virial shock or kinetic mode feedback from AGN) preventing further infall onto the galaxy by counteracting cooling (e.g., Rees & Ostriker 1977;Blumenthal et al. 1984;Birnboim & Dekel 2003;Bower et al. 2006;Croton et al. 2006;Somerville et al. 2008) or stabilization by a massive central mass component (e.g., Martig et al. 2009;Ceverino et al. 2010).
However, at  ≳ 1, when the most massive quiescent galaxies start to form (e.g., Thomas et al. 2005Thomas et al. , 2010)), observations provided inconclusive results to explain the connection between the removal of cold ISM and these potential quenching mechanisms.Stacking analysis studies found that there is a non-negligible amount of cold dust and gas in quiescent galaxies at  = 1 − 3 (Gobat et al. 2017;Magdis et al. 2021;Blánquez-Sesé et al. 2023), while individual observations on quiescent galaxies did not offer a convergence with cold gas and dust content in the high redshift quiescent galaxies (Sargent et al. 2015;Hayashi et al. 2018;Belli et al. 2021;Whitaker et al. 2021;Caliendo et al. 2021;Suzuki et al. 2022;Kalita et al. 2021;Williams et al. 2021;Morishita et al. 2022;Gobat et al. 2022).The detection and non-detection of cold gas in post-starburst galaxies at  > 0.5 also raised a similar question germane to quenching at high redshift (Suess et al. 2017;Spilker et al. 2018;Woodrum et al. 2022;Bezanson et al. 2022).Overall, there is a large scatter of the observed dust and gas content in the individual quiescent and post-starburst galaxies at high redshift.These observational results perhaps indicate that quenching star formation may not occur from a single mechanism or cannot be explained with simplified bimodal origins, making the full picture of galaxy evolution more complicated.
Whatever the feedback mechanisms are, recent studies advocated the need for bimodal quenching paths in the time scales -slow and fast -both in observations (e.g., Barro et al. 2013Barro et al. , 2016;;Schawinski et al. 2014;Yesuf et al. 2014;Wu et al. 2018;Belli et al. 2019;Carnall et al. 2019;Tacchella et al. 2022) and simulations (e.g., Rodríguez Montero et al. 2019).At  ≳ 2, spectroscopically confirmed massive, quiescent galaxies corroborated the fast quenching mode (e.g., Glazebrook et al. 2017;Forrest et al. 2020;Stockmann et al. 2020;Valentino et al. 2020).How do dust and cold ISM respond to this bimodal quenching time scale and especially fast quenching paths at high redshift?And what can we learn from the constraints of observed dust in high redshift quiescent galaxies to understand galaxy evolution?
This work builds upon a very deep observation with ALMA (∼ 20 hr in total) -one of the deepest images ever made with ALMA at this redshift -targeting a main-sequence galaxy at  = 2.21 (Q2343-BX610).We detect six continuum sources above peak signal-to-noise ratio (/ peak ) of four, five of which were made serendipitously.This Paper delivers the first analysis of these sources, focusing on the discovery of a candidate dust-rich quiescent galaxy at  ∼ 2. Comparing its properties with other high redshift quiescent galaxies available, we allude to ideas to answer the questions above.
The Paper is structured as follows.Section 2 describes the ALMA observation, data analysis, and flux measurement of the continuum sources.Section 3 illustrates the multi-wavelength (optical and nearinfrared, opt/NIR) data ( §3.1), and SED analysis of the quiescent galaxy candidate.It includes the photometric redshift estimate ( §3.2), and spectral energy distribution (SED) fitting ( §3.3).Section 4 describes overall SED fitting results, SFR constraints ( §4.1), UVJ color and stellar age ( §4.2), dust mass measurements ( §4.3).Section 5 delves into the details of the identified quiescent galaxy, by comparing it with respect to high redshift quiescent galaxies in terms of age and offset from the main-sequence ( §5.1) and dust-to-stellar mass ratio at given age.( §5.2).We present a discussion of the potential quenching mechanisms ( §5.3) of G4, followed by the potential origin of the observed dust at  ≳ 1 ( §5.4).Finally, we conclude and summarize our results in Section 6.We adopt a standard ΛCDM cosmology with  0 = 70 km s −1 Mpc −1 and Ω  = 0.3 and Chabrier initial mass function (IMF; Chabrier 2003).

OBSERVATIONS AND SOURCE DETECTION
The original program aimed to study one of the well-studied, typical main-sequence galaxies at  = 2.21, Q2343-BX610 (hereafter BX610), at an angular resolution of 0 ′′ .2(ID: ADS/JAO.ALMA#2019.1.01362.S, PI: R. Herrera-Camus).Band 4 receivers were used to cover both CO (4 − 3) and [C I] ( 3  1 − 3  0 ) (hereafter, [C I]) emissions to study the cold interstellar medium properties and gas kinematics to resolve the inner structure of the galaxy down to ∼ 1.5 kpc scale, close to the Toomre length at this redshift.Therefore, antenna configuration was chosen to achieve ≈0 ′′ .2resolution, with baseline lengths ranging between 15 m and 3.7 km using 44-48 antennas.On-source time was 13.6 hrs (≃20 hrs in total including the overheads), which enables the characterization of the cold molecular gas kinematics including radial gas inflow in addition to the global ordered rotational motions for the target galaxy BX610.The detailed analyses on the main target, BX610, are presented by Genzel et al. (2023) and S. Arriagada and R. Herrera-Camus et al. 2023 (in preparation).This deep integration then allowed us to detect other sources within the ALMA field of view (FoV) accordingly.
Four spectral windows (SPWs) are used, two of each placed in the upper and the lower sideband, respectively, to cover the lines and the underlying continuum.Two SPWs are set in the Frequency-Division Mode (FDM) to detect the redshifted CO (4 − 3) and [C I] with a channel width of 3.906MHz (∼ 8.2 km s −1 ) covering 1.875 GHz bandwidth.The remaining two SPWs are observed in the Time-Division Mode (TDM) with a 2.0-GHz bandwidth at the 15.6-MHz resolution to cover the dust continuum at 2 mm.
We used CASA (McMullin et al. 2007) version 6.1.1.15for the calibration of visibility data and imaging.For visibility calibration, we used the pipeline script provided by the ALMA Regional Center staff.Continuum and cube images are produced by CASA task, tclean, and deconvolved down to 2 noise level, which is initially measured from the dirty map.
The continuum image is obtained using the line-free spectral windows observed in the TDM mode.The synthesized beam is 0 ′′ .23×0′′ .22,0 ′′ .22×0′′ .20 and 0 ′′ .22×0′′ .21for CO (4 − 3), [C I] and 2 mm continuum with natural weighting, respectively.Tapered images are also created by setting uvtaper parameters of 0 ′′ .15 and 0 ′′ .35 in the tclean task to check the presence of extended emissions.We then circularize the beam for these tapered maps, by specifying the restoringbeam parameter as 0 ′′ .29 and 0 ′′ .48respectively; so the final beams for the tapered images are 0 ′′ .29×0′′ .29 and 0 ′′ .48×0′′ .48.With the natural weighting, the typical continuum noise level around the phase center (BX610) is ≈ 3 Jy beam −1 and it goes up to ≈ 7 Jy beam −1 around the area close to the farthest continuum detected source from the center.
We identified six dust continuum sources with a signal-to-noise ratio of the peak flux greater than four in the ALMA map as shown in Figure 1.The ALMA photometry is performed based on a fixed aperture and cross-checked by the result with the 2D Gaussian fitting using imfit and the growth curve.We took the aperture size of 1 ′′ .0,encompassing most of the emissions from the sources, which are compact.We list the continuum fluxes for the continuum detected sources in Table 1.We checked the publicly available data on the ALMA archive observed at shallower depths and different configurations (#2013.1.00059S, 2017.1.00856.S, 2017.1.01045.S) in Band 4. We confirmed that BX610's flux values in those programs are consistent with those reported here within the errors1 .We note that the BX610's 2 mm continuum flux measured in Brisbin et al. (2019) is a factor ∼ 2 larger than our measurement at the same frequency.Although it is unclear whether this difference can be attributed to calibration errors in the flux in the Plateau de Bure Interferometer (PdBI) observations, our flux estimate is the conservative estimate of the total flux.This ensures our later discussion on the dust mass and dust-to-stellar mass ratio for G4.
We also checked the ALMA archive if there is additional data available at higher frequencies, finding that two ALMA project (#2015.1.00250.S and #2019.1.00853.S) covers three galaxies (BX610, G3, and G4) at Band 6 (1.2 mm) within the field of view.G3 and G4 are on the edge of the coverage where the sensitivity is degraded by 50-60% with respect to the central position.We find ≈ 5 and ≈ 3 peak features in G3 and G4, respectively using #2015.1.00250.S.The resolution of #2019.1.00853.S's Band 6 program is 0 ′′ .087×0.′′ 076 for natural weighting, where we do not find a noticeable signal at the position of G4 but G3.Photometry of Band 6 is also listed in Table 1 using the same method applied in the 2 mm map.We note, however, that the uncertainty of the photometry is larger due to the low signal-to-noise ratios compared to the 2 mm detection which needs deeper follow-up observations.This paper will focus on the analyses of G4 tification is performed on all objects described in the next section.
We defer the full discussions on other continuum and line-detected sources from the ALMA program to another paper (M. Lee et al. in preparation).Instead, we present a summary of the spectroscopic redshift of the 2 mm-continuum detected galaxies as follows and in Table 1.Among six ALMA continuum sources, CO and [C I], lines are detected in BX610 and G1, and one line for another (G5).The redshift of G1 is 2.21 based on CO (4 − 3) and [C I] line detections, confirming it being at the same redshift as BX610.With a single line and the photometric redshift described below, G5 is less likely associated with the BX610-G1 pair system and is rather a background source at  ≈ 3.01.

Counterpart identification and photometry at shorter wavelengths
We gathered all available data from optical and near-/mid-infrared (NIR/MIR) observations from the ground (  , , R  , , ,   ) and space (HST; F814W/F110W/F140W, Spitzer; IRAC Channel 2, MIPS 24m) and archival data from ALMA.Due to the very different point-spread functions (PSFs, ranging from ≈0 ′′ .15 to ≈5 ′′ .5) in different wavelengths, we choose to use different, reasonably large aperture sizes and correct it to recover the total fluxes (so that color corrections have less effect on the photometry for the same aperture correction) described as follows.
, , and R  band images were obtained by the Palomar telescope using the Large Format Camera (LFC) in 2001.The description of the data and analysis of the observations (Q2343 field) is available in Steidel et al. (2004) and the references therein.The   images are also from the Palomar 5m telescope, taken with the Wide Field Infrared Camera (WIRC) in 2003 and 2004.We refer the readers to Erb et al. (2006) for the details of the observations and data analyses.We used  and  band images taken from Magellan/FourStar, obtained   in 2013.The  image is a stack of the 2 and 3 intermediate bands, and the  band is a stack of 1 and 2.All of the Magellan data were taken by Gwen Rudie.All of the above photometry was performed using SExtractor, after matching the point-spread function (PSF) in all bands (PSF ∼ 1 ′′ ) using the   (for most objects) or  band (for very red objects) image as the detection band in "dual-band" mode.This strategy is to optimize the source detection and the corresponding isophotal aperture given that UV color-selected galaxies, the original main target of the observations, are better detected in the optical band than in the NIR, while the red objects, like G4, are conversely detected in the NIR bands with higher S/N.The magnitudes were corrected from isophotal to "total" using the SExtractor mag_auto in the relevant detection band.
For HST, we obtained the imaging data from the MAST archive processed with Grizli2 pipeline (Brammer & Matharu 2021;Brammer et al. 2022).The original work using the WFC3/F110W image was presented in Tacchella et al. (2015).We used the Spitzer data from the golfir3 project (Brammer 2022) which generated the Spitzer mosaics from the pipeline processed data.The maps were corrected for the astrometry.
For the HST bands and IRAC 4.5m photometry, we used SEP4 (version 1.2.0) to extract sources.For the same reason of different PSF, the source extraction and photometry measurement were done independently between the HST and Spitzer maps.For the HST bands, we used the longest wavelength (F110W (G1) or F140W (other galaxies)) available for source detection.We used an aperture of 0 ′′ .7 for the photometry in all bands and corrected it to "total" values within an elliptical Kron aperture (Kron 1980).The latter was calculated based on the longest HST bands (F140W or F110W), and the same amount of aperture correction is applied for all the other bands.Given the aperture size is large enough the color correction is negligible.For the IRAC photometry, we used a larger aperture of 3 ′′ .6 and applied the aperture correction which ranges between 10-18%.We also checked the curve of growth of the detected galaxies in the IRAC map, confirming that the aperture correction and thus the estimated flux is reasonable.
The MIPS 24m data was taken in 2006 and the description of the data reduction is available in Reddy et al. (2010).The MIPS photometry was done using an aperture of 3 ′′ centered at the 2 mm sources, and an aperture correction factor of 3.097 was applied based on the MIPS Instrument Handbook (v3.2.1)Table .4.13.We obtained the errors based on the median value of standard deviation values given the same aperture taken from random, emission-free positions around the source.
Counterpart identification was performed based on the position of the 2 mm emission.We identified the optical/NIR/MIR counterpart at 0 ′′ .10-0 ′′ .24offset from the 2 mm peak positions except for G2.Considering the beam and the point-spread function (PSF) sizes of ALMA and optical/NIR observations, and the signal-to-noise ratio, we regard the offset as negligible to pin down the counterparts.The offsets would rather be originated from the internal structures and different levels of dust extinction across the galaxy.Finally, we note that there is no additional counterpart identified in the HST imaging within the IRAC and MIPS apertures that could contribute to the photometry.
Figure 2 shows the three HST filter and Spitzer/IRAC Ch2 (4.5m) maps overlaying the 2 mm ALMA continuum contours.G2 does not have clear counterparts in any of the bands which may be partly owing to the poor sensitivity (i.e., on the edge of the coverage) around the source.G3 starts to appear in the J, F140W maps and the longer wavelengths, which may imply either an obscured (indicated by the dust detection) or a Lyman/Balmer-break galaxy at higher redshift.G2 and G3 have insufficient data points to perform the full SED fitting and get a good photometric redshift estimate.Followup observations of these two are necessary to characterize these two optically dark/faint galaxies.Hereafter, we consider only G4 for further analyses.

Photometric redshift
We use eazy-py5 (Brammer et al. 2008) to get a first handle on the photo- estimate.Full details concerning eazy-py and the template set can be found in Brammer et al. (2008); Kokorev et al. (2022); see also Blanton & Roweis (2007) for the template creation algorithm.Briefly, the best-fit photo- is estimated from a linear combination of a set of 13 templates from the Flexible Stellar Populations Synthesis code (fsps, Conroy & Gunn 2010) and one additional template recently obtained from  = 8.5 -ID4590 from Carnall et al. (2022) (namely, carnall_sfhz_13 template in eazy-py) 6 .The first thirteen templates are constructed based on redshift-dependent SFHs (a combination of four star-formation histories) and 2-3 dust models.Star-formation histories that start earlier than the age of the Universe are excluded at a given redshift.We note that there is a limitation of this set-up for eazy-py that there will be -SFH degeneracies if there are an insufficient number of bands.We explore different SFHs using different SED fitting codes in Appendix A which does not alter our conclusion on G4.The set of thirteen templates is chosen based on an algorithm to reasonably cover the rest-frame color space following the philosophy of Brammer et al. (2008).The additional template can better explain the extremely strong emission lines seen in early JWST spectra of z>6 galaxies.Templates are constructed on the Chabrier (2003) IMF.Kriek & Conroy (2013) dust attenuation law (dust index  = -0.1, V = 3.1) is adopted and its maximum attenuation is set to be redshift-dependent.On the fsps thirteen templates, a fixed grid of nebular emission lines and continuum from cloudy (v13.03)models are added (metallicity: log (/ ⊙ ) ∈ [−1.2, 0], ionization parameter log () ∈ [−1.64, −2]7 ; Byler 2018).Given the unconstrained nature of fixed parameters used in eazy-py, we will only use the photo- information from eazy-py.Finally, a correction for the effect of dust in the Milky Way ( ( − ) = 0.0327) is applied to the templates within eazy-py, pulling the Galactic dust map by Schlafly & Finkbeiner (2011) from dustmaps (Green 2018).For the fitting process, we set minimum and maximum redshifts to 0.01 and 12, respectively, and the redshift step (z_step) to 0.01.Referring to Brammer et al. (2008), we use the extended K-band prior (prior_K_extend.dat) to lower the weight of the low redshift solution.
To gain more credentials of photo- estimate from eazy-py, we also run bagpipes and cigale using different SFH.The details of the SED results and the set-up to fit redshift are summarized in Appendix A. In general, they agree well with each other, giving the best-fit photometric redshift of  ≃ 2 for G4.We overlay the final results in Figure 3 for the () distribution from eazy-py, bagpipes, and cigale photo-z constraint.Hereafter, we consider G4 to be at  ∼ 2.

SED fitting
Given the consistency between the estimated photo- from different SED codes, we use bagpipes8 (Carnall et al. 2018), fast++9 (Kriek et al. 2009;Schreiber et al. 2018) and cigale 10 (Boquien et al. 2019) by fixing the redshift based on the best-fit photometric redshift from eazy-py ( = 2.13).The photometric data used for G4's SED fitting is summarized in Table 2.We use three different codes to see if these codes are able to fit the G4's photometry and thus strengthen the reliability of the overall results.We regard the results with fixed redshift as the best estimate from SED, which is summarized in Table 3.In Appexdix A, we explore the impact of fitting the photometry with the redshift left free to vary.The comparison indicates reasonable agreement between parameters derived with and without fixing the redshift, and does not alter the overall conclusion that G4 is a massive quiescent galaxy at  ∼ 2 with low specific star-formation rate log( /yr −1 ) ≲ −12.5.

bagpipes
We take into account two different types of parametric star-formation histories (SFH) to fit the data points: double-power law SFH ( () ∝ (/)  + (/) − −1 ) and non-parametric SFH.These two sets of SFH were chosen because an exponentially declining SFH (-model) model did not deliver a good fit for G4 using bagpipes giving high  2 without reasonable convergence.Although the -model (e.g., Schmidt 1959;Conroy 2013) is the most popular model that has been extensively used in studies based on large surveys (e.g., Ilbert et al. 2013;Skelton et al. 2014;Stefanon et al. 2017) recent studies advocated the need for more complex or different models to explain the star-formation histories of galaxies (e.g., Maraston et al. 2010;Wuyts et al. 2011;Pacifici et al. 2016;Carnall et al. 2018;Leja et al. 2019a).
The double-power law SFH has more flexibility with an additional parameter (a total of 6 parameters) compared to a -model.We vary  and  values between 0.01 and 10000, with uniform priors in log space ("log_10").The age of the Universe the peak of star formation () is sampled between 0.1 and 13.5 Gyr and the metallicity range is between 0 and 4.0  ⊙ .We adopt the reddening law of Calzetti  2000) with the -band attenuation (  ) between 0 and 4 magnitudes.Nebular emission is not taken into account for G4; we tested nebular emission with G4 and the nebular emission did not change our analysis.We allow the stellar mass formed by the assumed SFH in 7.0 < log( ★ / ⊙ ) < 13.5.
We also test with non-parametric SFH with continuity prior using bagpipes to test if non-parametric SFH can produce better results.By its construction, the non-parametric continuity SFH priors can yield a better fit and may overfit because we have more free parameters to fit.To set the SFH bin, we follow the methodology described in Leja et al. (2019a) with six bins.The first two recent bins are set to 0-30, and 30-100 Myr, and the remaining four bins are equally spaced in logarithmic time.Leja et al. (2019a) showed that the number of bins in the mock analysis is insensitive as long as  bins > 4. We also confirm that the number of bins has a negligible impact on stellar mass, age, and  V , and SFR for G4 for 4 <  bins < 9. Therefore, we adopt the bin number  bins = 6, which is computationally less expensive without excessively overfitting given the number of data points.For G4, we get the consistent stellar mass, and SFR, placing a galaxy far below the star-forming main sequence at  = 2, as summarized in Table 3.

fast++
We run fast++ with the following assumptions: Bruzual & Charlot ( 2003) library (bc03) assuming Chabrier (Chabrier 2003) initial mass function (IMF).We assume a delayed- model for the sake of simplicity and to allow more flexibility than the  model.The minimum -folding time is set to log  min = 8.0 (100 Myr) and maximum of log  max = 10.1 (12.5 Gyr) in steps of log  = 0.1 The minimum time since the onset of star formation is 50 Myr and the maximum is the age of the universe.The log(age) grid is made in steps of 0.05.Solar metallicity and the Calzetti dust attenuation law (Calzetti et al. 2000) are adopted with visual extinctions in the range 0 <   < 4.

cigale
We use cigale (Boquien et al. 2019) in order to take into account the ALMA data points (2 mm and 1.2 mm) to avoid the potential underestimation of the SFR from optical to mid-IR photometry alone, especially for an obscured star-forming case.We adopt the delayed- model for the SFH.The other set-up is similar to those used in the other SED fitting process that we used Bruzual & Charlot (2003) library assuming Chabrier (Chabrier 2003) initial mass function (IMF) and use the Calzetti dust attenuation law (Calzetti et al. 2000).As for the dust emission, we use the Draine et al. (2014) model constructed by varying the PAH mass fraction ( PAH = 1.12 − 3.90), minimum radiation field ( min = 0.5 − 10), and power law slope ( = 1.0 − 2.0), which describes the distribution of the radiation field per mass.We also set the metallicity and nebular emission to vary.cigale fits G4's SED well with reasonable  2 ≈ 11.3 (Figure 4).

Consistency in the SED fitting results and SFR from other methods
G4 is a massive galaxy with log ( ★ / ⊙ ) ≈ 11.0 based on the SED fitting.Best-fit SED results of G4 are summarized in Table 3.The stellar masses between different codes are overall in good agreement.The agreement in the stellar mass between different models is consistent with earlier studies that stellar masses are the most robust parameter estimated from SED fitting for "normal" (and not so blue) galaxies (e.g., Papovich et al. 2001;Shapley et al. 2001;Conroy 2013;Pacifici et al. 2023).
All codes yield broadly consistent results favoring the quiescent nature of G4 though there are different levels of inactivity from the SED fitting.The multiple SED codes indicate that G4 is characterized as a quiescent galaxy with a low level of star-formation activity with log( /yr −1 ) ≲ −12.5, which is -3.5 dex lower than the starforming main sequence at  = 2.The cigale fit, where we include the ALMA data points for the fitting with energy-balance assumption also gives a low level of star-formation and a quiescent solution for G4 with a reasonable  2 value.The fitting results do not change with and without the 1.2 mm data point for G4.We show the best fit from cigale in Figure 4 and the other best-fit SED models in Appendix A. However, the SFR from the SED could introduce artificially low SFR, especially for the parametric SFH models by its construction.In this regard, the attenuation corrected R  -band (i.e., rest-frame 2000 Å) magnitude suggests the SFR of ∼ 8  ⊙ yr −1 and log( / −1 ) ∼-10.3 instead: we applied the implied extinction correction ( V ∼ 0.5) assuming the Calzetti law, i.e., (2000Å)∼ 9× (B − V) ∼ 1.5 mag, which gives  ∼ 24.5 mag and converted this into SFR attributing the emission to UV continuum of young stars.
Given the different degrees of low SFR from the SED best fit and inferred SFR from the R  -band photometry, we estimate the SFR from three other different ways using 24 m and 2 mm fluxes.These are more reliable tracers of SFR because they are not affected by dust obscuration.The estimated dust extinction ( V ) is moderate, which can be attributed mainly to the old stellar population (≈ 0.4 in  V ), but it is still good to check whether other different SFR measurements agree with each other.A summary of underlying estimates is presented in Table 4 If we use the modified black body (MBB) and take the 2 mm flux, the obscured star formation can range ≈ 2 − 136  ⊙ yr −1 assuming  d = 17-35 K and  = 1.5-2.2,using Kennicutt (1998) to convert  IR to SFR and dividing it by 1.6 for Chabrier IMF (Madau & Dickinson 2014).Here, we assume no contribution from old stellar populations heating the dust.In general, however, there is a non-negligible contribution of the old stellar population in dust emission for quiescent galaxies, which could play another source of uncertainty constraining the SFR of G4 from  IR and in effect the SFR can be lower than this.Further, without constraining higher frequency near the dust SED peak and the dust temperature, there is considerable uncertainty in inferring the actual star-formation rate (Table 4).ALMA Band 4 covers the rest-frame 660m at  ∼ 2, and should mainly trace the cold dust component rather than the SFR.Low S/N (≈ 1.5 for total flux, Table 1) of the 1.2 mm data point does not provide a strong indication of extremely high star formation, giving a similar range of 2 − 72  ⊙ yr −1 .Future highfrequency observations at ALMA Band 9, and 10 tracing the dust peak will help to pin down the SFR.Reddy et al. (2010) used the MIPS 24m flux for ⟨⟩ ∼ 2 galaxies to constrain the relationship between the rest-frame 8m luminosity, which traces the emission from polycyclic aromatic hydrocarbons (PAHs), and SFR.Using our 24m flux in Table 2 and equation (2) in Reddy et al. (2010), we get SFR of ∼ 30  ⊙ yr −1 (total dustcorrected SFR).The scatter of this fit is about 0.24 dex, giving a range of 16-55  ⊙ yr −1 from the 24m flux.However, the relationship between rest-frame 8m and SFR is constructed based on galaxies forming stars at a higher rate, which might be a source of uncertainty playing a role in this conversion.
The third approach is to estimate SFR by scaling the 24m flux with respect to Q2343-BX610, which has a better constraint with more photometric data points from UV-to-FIR, and use the latest SFR measurements of Q2343-BX610.Brisbin et al. (2019) measured the star-formation rate of BX610 to be 140  ⊙ yr −1 from their SED fitting including a few FIR data points.They estimated slightly higher than attenuation corrected SFR from UV-to-NIR (60  ⊙ yr −1 ) and H (115  ⊙ yr −1 ) (Förster Schreiber et al. 2014Schreiber et al. , 2018)).By taking the MIPS 24m photometry for BX610 and G4, we estimate the G4's SFR is ∼ 8-10  ⊙ yr −1 .
SFR from different tracers can change by two orders of magnitude (or more; Schreiber et al. 2018;Belli et al. 2019;Man et al. 2021) between   SED ,   IR+UV , and   [OII] or H for high redshift quiescent galaxies.Based on multiple approaches to estimate the SFR, we suggest that the SFR of G4 from the scaling from the 24m point (the third method, SFR∼10  ⊙ yr −1 ) as the most reasonable (or robust) estimate with available data.Inferring SFR with only a single data point from 2 mm needs heavy interpolation and the conversion between rest-8m to SFR might have unknown systematics for low luminosity regime.This also minimizes our concern of potential underestimation of SFR from the SED fitting.This also agrees with the dust attenuation corrected SFR from the rest-frame UV magnitude.Finally, the adopted SFR of 10  ⊙ yr −1 also places the galaxy to be in (somewhat) better agreement with the observations and simulations in the age-ΔMS plane (Section 5.1).
With the SFR value of ∼ 10  ⊙ yr −1 , G4 is a quiescent galaxy with lower sSFR (log ( /yr −1 ) ≈ −10.2) at least by ≈ 1.2 dex from the  = 2 main-sequence, compared to main-sequence galaxies at  ∼ 2. We note that log ( /yr −1 ) = −10.0(−10.5) at  ≈ 2 corresponds to the mass-doubling time being more than 3 (10) times the age of the universe at that redshift.One of the definitions of quiescent galaxies used in the literature is sSFR smaller than 0.3/ Hubble (Franx et al. 2008), and G4 satisfies this.We use this sSFR (log ( /yr −1 ) = −10.2) as a proxy of the galaxy's quiescence.
We note that available data can not completely reject the possibility of a higher star-formation rate than currently best guessed.Future observations at higher frequencies at submm covering the dust SED peak are needed to highlight potential obscured star formation that we may still miss.
As for the other parameters, we take the bagpipes results using non-parametric SFH in the following discussions, because the results are close to the median values of four different SED fitting results.Also, the SFR of G4 from bagpipes non-parametric gives one of the highest star-formation rates that we use as the upper limit constraint of the SFR (from the SED).

UVJ color and stellar age
We investigate whether G4 also satisfies the selection criteria of quiescent galaxies using rest-frame , ,  magnitudes, based on the best fit from bagpipes (non-parametric SFH).The  −  and  −  colors (hereafter   colors) are widely used to distinguish quiescent galaxies from reddened dusty galaxies and star-forming galaxies as originally demonstrated by Labbé et al. (2005) and Wuyts et al. (2007); star-forming galaxies have bluer  −  colors than the other two classes, and dusty star-forming galaxies tend to have redder  −  colors than their quiescent counterpart.This idea has been further supported by simulations for  ∼ 2 galaxies (e.g., Donnari et al. 2019, but see Akins et al. 2022).The   diagnostics become less reliable for  > 3 galaxies with decreasing purity both indicated by observations (e.g., Merlin et al. 2018;D'Eugenio et al. 2020;Forrest et al. 2020) and simulations (Lustig et al. 2021(Lustig et al. , 2023)).Nevertheless, the   color selection works fairly well to select quiescent galaxies at  ∼ 2 observationally.
We adopt the   color selection of Whitaker et al. (2011) for 2 <  < 3.5 galaxies which are shown as a black diagonal line and a dashed line in Figure 5.The diagonal black line separates the plane into two regions of quiescent and star-forming galaxies from their definition.Additional criteria of  −  > 1.2 and  −  < 1.4 (dashed line) distinguish unobscured and dusty star-forming galaxies, respectively.G4 is classified as a quiescent galaxy from these criteria.The   colors based on eazy-py (filled square in Figure 5) also show a good agreement.In Appendix A, we also list the best-fit   colors obtained from the best-fit SED without fixing redshift.All but one SED best-fit results give the rest-frame   colors of the quiescent population.One exceptional case is in the bagpipes run with a flat redshift prior using double-power law SFH, where the convergence was bad to reliably constrain the UVJ colors with large error bars (Table A1).The   colors of G4 also satisfy two other different quiescent galaxy criteria proposed by Williams et al. (2009) and Leja et al. (2019b).Williams et al. (2009) proposed a slightly different definition as a quiescence ridge line,  −  > 1.3,  −  < 1.6 and ( − ) > 0.88 × ( − ) + 0.49 for  ∼ 2 galaxies (a dashed-dotted line in Figure 5).The solid red line in Figure 5 is from Leja et al. (2019b) for log( /yr −1 ) = −10.5, in which the sample purity of quiescent galaxies is high enough (≳ 92%) in the   diagram using the 3D-HST galaxies.There is a tendency for galaxies with lower sSFR to be placed on the top left corners, more distantly perpendicular to the QG-SFG dividing line (thick solid black line in Figure 5).However, the sSFR distribution saturates beyond this sSFR limit such that   colors no longer give any further constraints.
Given that G4's   colors satisfy the quiescent criteria, we use this to infer the stellar age of the galaxy and compare it with the SED best fit.Calibrated with their deep spectroscopic observations, Belli et al. (2019) found a good correlation with the (median) stellar age (for > 300 Myr) and the   colors for 1.5 <  < 2.5 quiescent galaxies and the age distribution is along the diagonal quiescent ridge line.Belli et al. (2019) introduced a rotated coordinate system to calculate the age,   = 0.75( − ) + 0.66( − )   = −0.66(− ) + 0.75( − ), (1) and the following median stellar age ( 50 ) can be estimated by log( 50 / yr) = 7.03 + 1.12 •   .This age-color trend was also corroborated by (Carnall et al. 2020) based on the broad-band photometric data points for 2 <  < 5 quiescent galaxies.
The mass-weighted age from the SED matches well with the one from the empirical relation.The color of G4's symbol in Figure 5 represents the mass-weighted age (∼ 1.9 Gyr) based on the SED fitting using bagpipes (non-parametric SFH), as listed in Table 3.The small grey dots in the background show log( ★ / ⊙ ) > 10.5 galaxies at the redshift range of 2.0 <  < 2.5 from the COSMOS2020 catalog (Weaver et al. 2022) by taking the rest-frame   colors from their eazy-py fit.The number density of them is also plotted as contours.For visual comparison, we color-code the empirical agecolor relation proposed in Belli et al. (2019) using the COSMOS2020 galaxies within the quiescent region.The region marked in a blue box indicates the age range between 300 < age/Myr < 800, where post-starburst galaxies are expected to place according to Belli et al. (2019).Although we need spectroscopic confirmation to better constrain the stellar age, the SED-based mass-weighted ages (Table 3, ∼ 1 − 2 Gyr) and empirically derived age (≈ 1.2 Gyr) are broadly consistent within the errorbars.
Typical ages of quiescent galaxies with spectroscopic observations at  ∼ 2 are observed to be 1−2 Gyr, with a relatively large age spread (Kriek et al. 2006(Kriek et al. , 2009;;Whitaker et al. 2013;van de Sande et al. 2013;Mendel et al. 2015;Onodera et al. 2015;Fumagalli et al. 2016;Belli et al. 2019;Stockmann et al. 2020;Estrada-Carpenter et al. 2020).This is also noticeable, although perhaps being much crude, by the inferred age and density distribution of COSMOS2020 galaxies in Figure 5; the density peak of the "red cloud" is at ∼ 1 Gyr but there are many galaxies younger and older than this.Considering the age distribution in the   diagram, the estimated stellar age of the G4 (= 1 − 2 Gyr) appears to be reasonable and common at  ∼ 2.

Dust-to-stellar mass ratio
We estimate the dust mass of G4 based on a single measurement at 2 mm, assuming a modified black body with where  is the dust emissivity index, and  abs is the grain absorption cross section per unit mass, following the discussion in Bianchi (2013) and Berta et al. (2016).The lower dust temperature of a quiescent galaxy is indicated by Magdis et al. (2021), which performed the stacking analysis 11 of the SED and obtained the best fit using Draine & Li (2007) templates.The adopted value of 23.5 K is the value for  ∼ 2 quiescent galaxies 12 .The estimated dust mass of G4 is log( d / ⊙ ) = 8.51 ± 0.24 where the error takes into account the photometric error only.
We note two caveats in our dust mass estimate: first, the inferred 11 We note that they excluded the 24m detection sources above ∼45 Jy (3) and G4 is below the limit, satisfying their selection criteria. 12Magdis et al. ( 2021) noted a potential evolution in the dust temperature giving higher dust temperature with increasing redshift.If we take the mean value of 21 K estimated for  = 0 − 2, the dust mass of G4 becomes slightly higher by 0.08 dex, giving log(  d / ⊙ ) = 8.59 ± 0.24.2021) for massive (log(  ★ / ⊙ ) > 10.8) quiescent galaxies at  ≲ 2 are shown as grey squares.Grey shaded region with a dotted curve is based on the functional form for the best-fit at  < 1,  gas (=  gas / ★ ) ∝ (1 + ) 5 , assuming a fixed gas-to-dust mass ratio of 92 and normalized at  = 1 with  gas = 8%.For comparison to star-forming galaxies of similar mass, we calculated dust-to-stellar mass ratios using the scaling relation suggested by Tacconi et al. (2018), assuming a gas-to-dust mass ratio (GDR) of 100 for log (  ★ / ⊙ ) = 11.0 and 10.5 (red and blue shaded bands, respectively; more massive galaxies have smaller gas (and dust) at fixed GDR.
dust mass is a function of the assumed dust temperature, with higher  d corresponding to lower dust masses given an observed flux density.
The dust-to-stellar mass ratio for G4, assuming  d = 23.5 K, is log( d / * ) = −2.71± 0.26; if one instead assumes  d = 25 K, often adopted for normal star-forming galaxies (Scoville et al. 2016), the dust-to-stellar mass ratio reduces to log( d / * ) = −2.75-still at the high end of the range measured for quiescent galaxies (< −3); see also Figure 6.Second, the modified black body assumption we have used is a conservative estimate of dust mass.It is known that it gives systematically lower values compared to estimates based on a fully-sampled SED fitted with Draine & Li (2007) models.The current assumption of the  and the corresponding  values mitigates this discrepancy (Bianchi 2013;Berta et al. 2016).The dust mass of G4 could become even higher if we adopt Draine & Li (2007) and once the SED is well sampled; if this is the case, it will strengthen the dust-rich nature of G4.As further support of our estimated dust mass, we note that the dust mass derived from cigale, where we assumed the Draine et al. (2014) model (an updated model from Draine & Li 2007), is log( d / ⊙ ) = 8.62 ± 0.40, consistent with our measurement assuming the modified black-body.
In Figure 6, we show the dust-to-stellar mass ratio of G4 (log( d / ★ ) = −2.71± 0.26) based on the above calculation and stellar mass from the SED fitting using bagpipes (non-parametric SFH).For comparison, we plot the stacking results of quiescent galaxies with ⟨log( ★ / ⊙ )⟩ ≈ 11.0 and 0 <  < 2 from Magdis et al. (2021) (grey shaded curve and square data points after factoring in ≈ 0.21 dex for converting Salpeter (Salpeter 1955) IMF to Chabrier IMF in stellar mass).
The  d / ★ value for G4 is comparable with those based on stacking analysis (Man et al. 2016;Hayashi et al. 2018;Gobat et al. 2018;Magdis et al. 2021;Blánquez-Sesé et al. 2023) and higher than most of the individual gravitationally-lensed QGs (Whitaker et al. 2021;Caliendo et al. 2021) (Figure 6).Whitaker et al. (2021) put an order of magnitude lower limit of dust-to-stellar mass ratio (< −4 dex) for non-detected sources, although one must be cautious interpreting non-detections for more spatially-extended lensed sources (Gobat et al. 2022).Individual dust measurements and upper limits of non-lensed sources from Hayashi et al. (2018); Suzuki et al. (2022); Kalita et al. (2021) are also shown using the same modified blackbody model assumptions for ,  abs , and  d in Eq.( 2).Most individual data points have dust-to-stellar mass ratios at least a factor of 3 lower than G4; the exceptions are ALMA.14( = 1.4;Hayashi et al. 2018) and ZF-COS-19589 ( = 3.7; Suzuki et al. 2022), which have inferred dust-to-stellar mass ratios comparable to or even higher than G4 (labeled in Figure 6.) ALMA.14 is located at the outer edge of the core region of a galaxy cluster at  = 1.413 .It is detected in both CO (2 − 1) and [O II] (narrow-band); while it has a low SFR ∼ 3 M ⊙ yr −1 based on an SED fit to the optical/near-IR photometry (see Table B1), this estimate may significantly under-estimate the actual SFR, as we have discussed (Section 4.1).We recall that G4's best-fit SFR using the same star formation history parametrization ( model) and SED fitting code (fast++) suggested SFR∼ 0.01 M ⊙ yr −1 .Given the CO and [O II] detections of ALMA.14, its status as a QG may be questionable.
ZF-COS-19589 has an inferred SFR comparable to our best estimate SFR of G4 (Section 4.1) based on a dust-corrected H luminosity: 10 +28.5 −12.9  ⊙ yr −1 ; Schreiber et al. (2018).Its SFR inferred from an SED fit is 0.00 +4.41  −0.00  ⊙ yr −1 .Suzuki et al. ( 2021) estimated the   IR , extrapolated from a single continuum detection at an observed wavelength of 870m giving a range of 6 − 33  ⊙ yr −1 depending on the assumed dust temperatures (between 20 K and 40 K).These inferences suggest that ZF-COS-19589 may have properties similar to G4, with M * smaller by ∼ 0.3 − 0.4 dex but with a comparable deviation (−1.2 dex) from the star-forming main-sequence for its redshift.
In Figure 6, we also plot the dust-to-stellar mass constraints of the star-forming main-sequence galaxies for different stellar mass bins of log ( ★ / ⊙ ) = 11.0 and 10.5.We used Tacconi et al. (2018) relation to get gas fraction and converted it to dust-to-stellar mass ratio assuming a fixed gas-to-dust mass ratio (GDR) of 100.The shaded region corresponds to ±0.3 dex scatter of the main sequence.The adopted GDR is reasonable for these two massive populations with high metallicity e.g., Rémy-Ruyer et al. (2014).For a higher GDR, the shaded bands would go lower and vice versa.
If the Kennicutt-Schmidt relation is also applicable to this galaxy, we expect the gas-to-dust mass ratio of G4 should be higher than 100 to explain the observed dust-to-stellar mass ratio.A high gas-to-dust mass ratio for lensed (less massive) quiescent galaxies at  ∼ 2−3 was also advocated in Whitaker et al. (2021) based on simba simulation.To put it in a different way, the galaxy's depletion time scale is long (≳2.5 Gyr) at  ∼ 2 with the estimated SFR and observed dust at fixed gas-to-dust mass ratio.We will explore the dust-rich nature of the galaxy in the next section (Section 5.2).
Taking these all into account, G4 is highly dust-rich with a dustto-stellar mass ratio of log( d / ★ ) = −2.71± 0.26, compared to other high redshift quiescent galaxies reported in the literature, except for two and the stacked results.Considering its current starformation rate at a low level (log( /yr) ≲ −10.2), unless there is a process igniting star formation at later times, it is likely G4 will remain gas(dust)-rich while being "quenched" for a long period of time (≳2.5 Gyr).

Morphology
The HST composite (Figure 2) image shows that G4 has a central core-like component and an extended disk.Here, we quantify G4's morphology by fitting the 2D light distribution with the Sérsic profile based on the F140W map for later discussion.
We use statmorph (Rodriguez-Gomez et al. 2019) 14 and galfit (ver3.0;Peng et al. 2002Peng et al. , 2010)).We fit the magnitude (), half-light radius (  ), Sérsic index (), axis ratio (), position angle (PA), and the central position.The initial guess parameters in the galfit are taken from the statmorph result.The fit gives  = 1.66 ± 0.03 (statmorph), where the error is estimated based on the bootstrapping after 1000 realization, and  = 1.66 ± 0.11 (galfit), respectively.The visual inspection of the residual images for both fits suggests the fit is reasonably good.
We also obtained Gini (Abraham et al. 2003) and  20 (Lotz et al. 2004) non-parametric values from statmorph.They are 0.52 ± 0.01 and −1.85 ± 0.03, respectively.These values locate the galaxy close to the borderline between Sb/Sbc and E/S0/Sa classification based on Lotz et al. (2008).

A DUST-RICH QUIESCENT GALAXY IN DISTANT UNIVERSE
The available opt/NIR/MIR photometry and best-fit   colors (although these two are not independent measurements) suggest that G4 is a quiescent galaxy with low sSFR and mass-weighted age of ≈ 1 − 2 Gyr at  ≈ 2. Meanwhile, ALMA observations reveal that the galaxy has a high dust-to-stellar mass ratio, comparable to the stacked results which are (at least threefold) higher than most of the quiescent galaxies individually studied.In this section, we take advantage of these findings and compile available measurements in the literature to address the implication of G4's properties in the context of quenching at high redshift.

Age-ΔMS relation: passive evolution at 𝑧 > 1?
If galaxies evolve passively after the quenching, a negative correlation between the age and the distance from the main sequence is expected.
With this idea, we first search for a signature of passive evolution at  ≳ 1 in the age-ΔMS relation.We gather the measurements of  ≳ 1 quiescent galaxies in the literature from Estrada-Carpenter  2020); Kubo et al. (2021) for  > 3 sources, all of which are plotted in Figure 7. First, we investigate  ∼ 1 quiescent galaxies based on Estrada-Carpenter et al. ( 2020) and Belli et al. (2019) which provide 177 quiescent galaxies in total.The age ( 50 ), stellar mass, and SFR measurements in these studies are estimated in a similar way, alleviating the systematic errors between the two.We choose a subset of goodquality data sets with the following criteria: galaxies having  ★ errors within 0.2 dex with the stellar mass of log( ★ / ⊙ ) > 10.5 and SFR errors within 0.6 dex.The adopted errors would sufficiently incorporate the potential systematic uncertainties in the different SED codes, as demonstrated by Pacifici et al. (2023).The final number used in the fitting is 147 in total.Galaxies with error bars in Figure 7 show the subset and are used in the following fitting procedure.
With the following functional form, we perform the Bayesian Markov chain Monte Carlo fitting using PyMC (Salvatier et al. 2016) The Spearman coefficient (with the good quality data) is −0.43 (thus negative) with -value of 7-8.Therefore, we confirm that the anticorrelation of mass-weighted stellar age and the deviation from the main sequence exists at  ∼ 1, indicating passive evolution.The best fit is shown as a dashed blue line in Figure 7.
Compared to the situation at  ∼ 1, confirming the existence of age-ΔMS relation at  ≳ 2 is largely limited by the data.At  ∼ 2, Stockmann et al. (2020) provided mostly the upper limit constraints of SFR from the SED fitting.At least for those detected in 24m (table 2 in their paper), they are located close to G4 (i.e., log( ) ≃ −10.2 and log(age) of ≃9.2).Higher redshift quiescent galaxies at  ≳ 3 (Schreiber et al. 2018;Valentino et al. 2020;Forrest et al. 2020;D'Eugenio et al. 2020;Kubo et al. 2021) are typically younger than G4 because the age of the universe gets younger.While there is a hint of negative correlation with lower normalization in  ∼ 3 quiescent galaxies from the visual inspection, substantiating the existence of the relation  ≳ 2 (as evidence of passive evolution) requires more observational data.
For completeness, in Figure 7, we mark high- quiescent galaxies with dust continuum and/or CO/[CI] line measurements by compiling all available information ( ★ , age, SFR, dust/CO/[C I]constraints) from the following papers: Onodera et al. ( 2012 (2023).We provide a summary in Table B1.Although it is tentative, it seems that many galaxies detected in dust and/or CO/[C I] are more scattered with respect to age-ΔMS relation.We investigate this more by connecting the dust-to-stellar mass ratio together with age and sSFR in Section 5.2.
The anti-correlation at  ∼ 1 is also discernible in the cosmological simulation in simba simulations (Davé et al. 2019 15 ; Figure 8),  2020) good quality samples (see Section 5.1).Highz ( > 1) quiescent galaxies with dust/CO/[C I] observations are also plotted, taken from the various literature (see Section 5.2 and Table B1).
Figure 8. Age-ΔMS relation in simba simulations at  = 2.0 (left) and  = 1.0 (right).Data points are for massive (log  ★ > 10.5) quiescent galaxies (defined as log(ΔMS) < −0.3) color-coded by the dust-to-stellar mass ratio.The distribution of  = 0 quiescent galaxies are also shown in grey contours on the right panel.We overlay G4's position in star symbols for the  = 2.0 snapshot.The dashed line on the right panel is the best fit obtained from the observations of  ∼ 1 quiescent galaxies shown in Figure 7.
Our best-fit result of  ∼ 1 quiescent galaxies (dashed line) is shown in the  ∼ 1 snapshot (Figure 8 right), which is comparable to the simulation result.The consistency in the observation and simulation at  = 1.0 and the evolution down to  = 0 (shown in grey contours) suggest that ΔMS-age relation is likely well established up to  ∼ 1 for quiescent galaxies.
The correlation at  = 2.0 is not clearly recognizable in the simulation (but a tail exists).Mendel et al. (2015), with their spectroscopic campaigns (VIRIAL survey), demonstrated that stellar ages of  ≲ 1 quiescent galaxies are reasonably explained by passive evolution while at  ≳ 2, galaxies' pathways after the quenching (or while being quenched) and the evolution may not be simply a "passive" evolution.The lack of age-ΔMS correlation at  ∼ 2 in the simulation perhaps corroborates the finding from the VIRIAL survey.
Meanwhile, simba simulations hardly find galaxies displaying similar properties of G4 for its given ΔMS and age and considering the associated uncertainties (see black contours in the left panel of Figure 8).The distribution of G4 and available  ∼ 2 quiescent galaxies highly encourages more observations to check if there is any difference between the observation and simulation at  ≳ 2 and whether G4 is instead a rare population in the ΔMS-age relation.Finally, the SFR from 24m shows a better agreement with the distribution of simulated galaxies, providing indirect support for our choice of the SFR.

Long depletion time scale
If galaxy quenching requires the removal/consumption of gas and dust, we can conjecture an additional correlation (that may not necessarily be linear) with the gas and dust content with the deviation from the main sequence (or sSFR) and age.Initial reduction/consumption of gas and dust (just after quenching) would have an imprint on the strength of the initial quenching mechanism, while the following evolution would give us a hint on additional mechanisms to remove/consume/supply gas and dust out of/to the galaxy.Such a possibility is hinted at in the simba galaxies as shown in Figure 8 where the colors represent the dust-to-stellar mass ratio at given ΔMS and age; we find a decrease of  d / ★ along the "passive evolutionary" sequence identified in Section 5.1.Observations also found such a trend between  d / ★ and sSFR in the local galaxies (e.g., Cortese et al. 2012;De Vis et al. 2017).All these above lead us to connect the observed sSFR and age with the observed dust.

Dust depletion time scale
With the data set compiled in Table B1, we show the log( d / ★ ) as a function of log(Age) in Figure 9. Individual data points are from this work (G4), Gobat et al. (2022) 17 , Suzuki et al. (2022), and G4 is plotted as a star color-coded with its sSFR.We also plot other  > 1 quiescent galaxies with dust constraints from the literature (detection for circles, 3 upper limit with upside-down triangles for non-detection; Gobat et al. 2022;Whitaker et al. 2021;Suzuki et al. 2022;Kalita et al. 2021).The three lines are inspired by Michałowski et al. (2019): dashed line is the fit based on the local Herschel selected early-type galaxies, while the other two lines are tweaked and rescaled assuming  = 0.7 Gyr (dashed-dotted) and 2.5 Gyr (dotted) to guide our eye for dust-poor high-redshift quiescent galaxies.The colors represent the specific star-formation rate of the galaxies.Kalita et al. (2021), giving nine sources in total.From this figure, the dust-richness of G4 at a given stellar age is a noticeable feature.It is also conceivable that there are two different subgroups in this plane: five out of nine galaxies display coherent decreasing sSFR with their age, while the remaining four do not with higher dust-to-stellar mass ratios including G4.In the following discussion, we call the former dust-poor quiescent galaxy and the latter dust-rich quiescent galaxy.
We compare G4 with these high redshift quiescent galaxies by taking the functional form (  d  ★ =  exp(−age/)) described in Michałowski et al. (2019).Michałowski et al. (2019) used the Herschel-detected early-type galaxies (ETG) in the local universe ( < 0.1), and measured the -folding time of  = (2.5 ± 0.4) Gyr.This is shown as a dashed line in Figure 9.
G4 is well aligned with the measurement of the local dust-rich ETG.The other (three) dust-rich quiescent galaxies are also close to this relation.They are MRG-2129 ( = 2.1; Whitaker et al. 2021;Man et al. 2021;Gobat et al. 2022;Morishita et al. 2022) -known to be rotation dominated (Toft et al. 2017;Newman et al. 2018b) and host AGN (Man et al. 2021) -, Galaxy-D (z=2.9;Kalita et al. 2021) -located in a group environment -and ZF-COS-19589 ( = 3.7 but uncertain redshift; Schreiber et al. 2018;Suzuki et al. 2022).The SFR levels are different by more than three orders of magnitude in these galaxies.We also note that except for G4, dust emissions are captured at shorter wavelengths (< 400 m) in the rest frame which is a subject of a larger uncertainty to get the true dust mass.Having this cautionary note in mind, the two most dust-rich quiescent galaxies (G4 and ZF-COS-19589) are aligned with the local relation measured in Michałowski et al. (2019).Intriguingly, the stacked data from Blánquez-Sesé et al. (2023) at ⟨⟩ ∼ 1.5 is located closer to G4 if we use the inferred stellar age using Belli et al. (2019) from Nonetheless, G4's dust-to-stellar mass ratio would still be higher than any of those sources by ≳ 0.3 dex.
their rest-frame   colors ( −  = 1.73,  −  = 1.05, weighted mean, giving log(age)=9.2), also being consistent with the trend of Michałowski et al. (2019) with longer depletion time.Donevski et al. (2023) explored 548 dusty quiescent galaxies in the hCOSMOS survey at 0.01 <  < 0.7, finding a similar evolutionary trend in the spiral quiescent galaxies as observed in Michałowski et al. (2019); elliptical quiescent galaxies show an almost flat trend with age, demonstrating that morphological type is an important factor of the scatter in  d / ★ .We recall the G4's morphology is a composite of disk and bulge structures (Section 4.4).Perhaps, G4 and other dust-rich quiescent galaxies could be the progenitors of dusty (spiral) quiescent galaxies which constitute ∼ 10% of the total quiescent population in Donevski et al. (2023), and of those in Michałowski et al. (2019).We may be witnessing the emergence of such (relatively rare) galaxy population observed at lower redshift.
If we rescale and match the functional form by hand for the remaining dust-poor quiescent galaxies, we estimate the dust-depletion time scale of  ≲ 0.7 Gyr as shown as a dash-dotted line in Figure 7. Whitaker et al. (2021) explored the nature of high redshift quiescent galaxies without dust detection.They showed a substantial drop (or a huge scatter) in the  d / ★ ratio at the (light-weighted)18 age of ∼ 0.3 − 0.5 Gyr for  ∼ 2 galaxies (Figure 5 in their paper).This corresponds to the (slightly longer) dust destruction time scale of theoretical expectations (0.1 − 0.4 Gyr, e.g., Draine & Salpeter 1979;Jones et al. 1994).They attributed the dust non-detection as an indication of extremely high gas-to-dust mass ratio ( GDR ) and of the mechanisms requiring strong dust destruction and shutting down star formation -and this may also hamper efficient grain growths.
Altogether, G4 may have experienced less efficient feedback (see also Section 5.4) that would be in the form of dust destruction by the supernova (SN) reverse shock (e.g., Nozawa et al. 2006), thermal sputtering by hot electrons (e.g., De Vis et al. 2017;Galliano et al. 2021) and astration (e.g., Clark et al. 2015) compared to dust-poor ones.Further, a relatively larger scattered distribution of dust-rich quiescent galaxies in Figure 7 could also mean some additional processes that play a role by postponing the full depletion of the dust in G4.

Quenching path and mechanism to maintain quiescence
The dust-richness with only a moderate SFR of G4 is puzzling because of the large inferred amounts of gas (from dust) to potentially (re-)fuel star formation if it cools fast enough.Unless there is a mechanism that hampers gas to form stars e.g., morphological quenching (Martig et al. 2009), it is difficult to explain simultaneously the quiescence and the dust-rich nature of G4.We describe in the following that G4 cannot be simply explained with a fast-quenching mode unless there is an additional mechanism which aids sustaining the dust-rich nature.
We revisit the toy model presented in Belli et al. (2019) using the python interface of fsps (Conroy et al. 2009), python-fsps (Johnson et al. 2022).Synthetic galactic spectral energy distributions and the   trajectories are constructed using the exponentially declining SFH model for two different -folding times of 100 Myr and 1 Gyr.The isochrones are based on MIST (Choi et al. 2016;Dotter 2016) and the spectral library is from MILES (Falcón-Barroso et al. 2011).We used the dust attenuation law of Calzetti et al. (2000) and the old stellar population attenuation (dust2) being 0.4 (which corresponds to roughly  V ≈ 0.4), which would also incorporate the  V values from the SED fitting (Table 3).We fixed the metallicity as solar metallicity.Two trajectories are shown in Figure 5.We also show the expected offset with dust which we indicate as the arrow of  V = 0.4.The trajectories would depend on SFH and dust attenuation but show the overall trends of the different quenching time scales.The information of G4's SFH before the quenching remains unconstrained (i.e., the best fit can be obtained in many parametric models as well as non-parametric SFH, see Section 3.3).Nevertheless, G4's colors and the expected age support that the galaxy has experienced rather a faster quenching episode based on the fsps model trajectories.
The fraction of rapid quenching mode increases from ∼ 4% at  = 1.5 to ∼ 23% at  = 2.2 (Park et al. 2022).The number would go higher if one adds a less bursty post-starburst track.Park et al. (2022) also found in the TNG100 simulation that rapid quenching in these galaxies is driven by AGN.For half of the cases, gas-rich major mergers trigger the central starburst and make a compact remnant, although TNG100 underpredicts the fraction of fast-quenching phase galaxies than the observations.
In a classical picture where AGN plays a role, it accompanies the efficient removal/destruction of gas and dust.If the fast quenching mode was driven by AGN, G4 is not aligned with this picture because it still has a considerable amount of dust at a given stellar mass.It can be a moderate starburst or AGN that does not efficiently or entirely remove or destroy dust, and the remaining dust (and gas) is not cooled and collapsed fast enough to form stars, perhaps by the gravitational potential of a bulge component.
G4's dust-richness at a given sSFR and age and long dust-depletion times may indicate a morphological quenching together with inefficient (or failed) feedback.G4's morphology (Section 4.4) is fitted with a moderately high Sérsic index,  ≈ 1.6.Gobat et al. (2017) measured high dust content for ⟨⟩ = 1.7 quiescent galaxies with high (= 3.5 ± 0.1) based on stacking analysis.Although G4 has a somewhat lower Sérsic index, our results still hint at a compact bulge-like structure that could stabilize itself and hamper the fragmentation of cold gas to form stars (see also discussions in Hjorth et al. 2014 which also proposed a need of morphological quenching for galaxies with excessive  d at low SFR regime).
We infer the Toomre  parameter of G4 foreseeing future followup observations.For a stable disk, the ratio of rotation velocity and intrinsic gas dispersion is connected to a gas fraction (  gas ) and Toomre Q parameter (Förster Schreiber et al. 2006;Genzel et al. 2011): MRG-2129 is one of the dust-rich quiescent galaxies that do not follow  d / ★ -age -sSFR relation (see Section 5.2.1, Figure 9).This galaxy is known to be rotation dominated (Toft et al. 2017;Newman et al. 2018b) with "high"  rot / 0 (stellar) value of ∼ 2−3.If we assume that G4's disk is also rotationally supported (/ = 2 − 3 with  = √ 2, though the gas distribution can be different from the stellar component), its gas fraction (  gas ≈ 0.2(20%), converted from the observed  d assuming GDR of 100) would suggest  ∼ 2 − 3. Future stellar and gas kinematics studies will verify if G4 also has a stable disk further supporting the morphological quenching.

High dust yield
In this last section, we briefly discuss the possible origins of the observed dust by quantifying dust yields.The discussion also aims and Cas A  ⊙ ).The gray shaded region indicates the range of the dust yield expected from theory for SN II without destruction, while the blue area represents a more realistic yield range due to dust destruction (≪ 0.1  ⊙ ).The pink line indicates the maximum theoretical yield possible for the AGB stars (0.04  ⊙ ) but the observational results claim yields lower by more than an order of magnitude (see details in the main text).G4 requires almost the maximum dust yields from both SN II and AGB stars if grain growth is not taken into account.
to provide support for inefficient feedback acted on G4 and the need for additional mechanisms to explain the observed dust.
We follow the calculation of Michałowski et al. (2010) but with slightly different mass ranges when calculating the number of stars,  ★ = 8 − 40  ⊙ for Type II SN (hereafter SN II), and  ★ = 1.5 − 8  ⊙ for AGB stars.With our estimated stellar mass from bagpipes (non-parametric SFH), one SN II should produce ≈ 0.13  ⊙ of dust, while ≈ 0.013  ⊙ for AGB star, to explain the observed dust mass (3.2 × 10 8  ⊙ ).
Both values are close to the high end of the dust yield required for SN II and AGB stars.For example, the maximum dust yields of the observed SN II are reported to be 0.45 − 0.8  ⊙ for SN1987A (Dwek & Arendt 2015;Matsuura et al. 2015), 0.3 − 0.5  ⊙ for Cas A (De Looze et al. 2017), and 0.03 − 0.23  ⊙ for Crab (Gomez et al. 2012;Temim & Dwek 2013;De Looze et al. 2019).However, these high values are for freshly formed dust.Theoretical yield in AGB stars is 0.04  ⊙ (e.g., Morgan & Edmunds 2003;Ferrarotti & Gail 2006;Ventura et al. 2012;Nanni et al. 2013Nanni et al. , 2014;;Schneider et al. 2014), which is only valid for a very narrow range of the stellar mass.The net yields of either origin are expected to be at least an order of magnitude lower if we take into account the destruction of the reverse shock for the SNe (e.g., Todini & Ferrara 2001;Nozawa et al. 2003Nozawa et al. , 2006;;Bocchio et al. 2016; blue shaded area in Figure 10) and the observational result of the AGB stars (e.g., Rowlands et al. 2012).Donevski et al. (2023) explored simba simulation and the chemical model of Nanni et al. (2020), and offered the need for prolonged dust grain growth in dusty quiescent galaxies at intermediate redshifts.
The mechanism allowing efficient grain growth for high redshift sources is yet fully understood, with a limited amount of data.Exploring this is beyond the scope of this paper and the grain growth remains a possibility of the observed dust in quiescent galaxies.
On the other hand, Gobat et al. (2017) and Belli et al. (2017) opened a possibility of the external origin of dust and gas for high redshift quiescent galaxies.Considering the downsizing in the concordance model, the majority of massive quiescent galaxies in the early universe are expected to be located in dense environments.At the time when the accretion rate from the large-scale structure is higher ( = 1 − 3), there may be also a higher chance of a rejuvenation event for a (temporarily) "quenched galaxy" no matter what quenching mechanism was playing a role initially.In this regard, Tacchella et al. (2022) found an indication higher fraction of rejuvenated galaxies (≈ 30%) in a massive halo (log( h / ⊙ ) > 13.0) compared to a lower halo (≈ 4%) for  ∼ 0.8 quiescent galaxies.
Based on the stellar-to-halo mass relation (Moster et al. 2010), the massive nature of G4 hints at the halo mass is log( ℎ / ⊙ ) ∼ 12.9, which indicates the galaxy is located in a "group-like" environment.Confirming the overdensity and the association of G4 with starforming galaxies (BX610 and G1) will further strengthen this idea 19 .We estimate up to a third of the observed dust can be originated externally based on a crude calculation following the major merger rates for galaxies in the CANDELS field (Duncan et al. (2019)), which is 0.3-0.4Gyr −1 for  ∼ 2 log ( ★ / ⊙ ) > 10.3 galaxies.Still, it would require a high dust survival rate with no destruction, likely assisted by efficient grain growth.Follow-up observation of stellar and gas kinematics and their alignment could give us a hint of the external origin, as indicated in local fast rotators (Davis et al. 2011(Davis et al. , 2013)).
Taking these all together, the amount of dust in G4 suggests an efficient survival of dust (or inefficient feedback on dust) before and during initial quenching, efficient grain growth afterward, and perhaps a rejuvenation event during/after its quenching.We performed the same calculation with other high-redshift galaxies (summarized in Figure 10), and two dust-rich quiescent galaxies (ALMA.14 and ZF-COS-19589) require similarly high dust yields.Finally, as discussed in Michałowski (2015), it is also worth noting that Salpeter (Salpeter 1955) and top heavier IMF would require a comparable amount of and even higher dust yield; Chabrier IMF is the most conservative estimate.

CONCLUSION
We reported the detection of cold dust in a massive quiescent galaxy (log( ★ / ⊙ ) ≈ 11), G4, at  ∼ 2. The galaxy is one of the six 2 mm continuum sources in the deep ALMA observations targeting a massive main-sequence galaxy at  = 2.21.We identified the optical/near-infrared counterparts for all but one galaxy.This paper focused on one among them, G4, whose photometric redshift is estimated at  ∼ 2, close to the redshift of the Q2343-BX610 and G1 which are spectroscopically confirmed.
The avilable photometry suggest low SFR well below the  = 2 main sequence (≈ 1.2 dex) with log( / ★ ) [ −1 ] ≈ −10.2.The quiescent nature of the source is also supported by the   colors based on the best-fit SED.
We compiled the data in the literature for  ≳ 1 quiescent galaxies and compared it with G4.We discovered the existence of a negative correlation between age and ΔMS in quiescent galaxies at  ∼ 1, suggesting a passive evolution up to  ∼ 1. Observational data at  ∼ 2 − 3 including G4 did not provide a strong indication of the negative correlation (and hence the passive evolution) with limited data points.simba simulations showed a comparable negative relation at  ∼ 1 but at  ∼ 2 G4's age and ΔMS is rare in their simulation box.More observational data is needed to determine whether G4 is unique (and rare) or whether the simulation underproduces such a population.
G4 exhibits unique features considering its estimated massweighted age (≈ 1 − 2 Gyr), sSFR, and dust-to-stellar mass ratio compared to other dust-poor quiescent galaxies: the galaxy is dustrich for its stellar mass (log( d / ★ ) = −2.71± 0.26) and its age (1-2 Gyr).Based on the dust-to-stellar mass ratio and age relation, a longer dust depletion time is suggested, similar to dusty quiescent galaxies observed in local and intermediate redshift ( = 2.50 − 2.75 Gyr).This may indicate a potential evolutionary connection with each other.
The   color trajectories imply that the galaxy has experienced a fast-quenching mode.However, morphological quenching along with inefficient (initial) feedback would be required to explain its dust-rich nature and G4's morphology supports this.The high dust yield of G4 supports an inefficient quenching mode requiring efficient survival of dust and rejuvenation, perhaps assisted by efficient grain growth.
Our discovery encourages many observational programs for follow-up.Future rest-frame optical and high-frequency submillimeter observations to constrain the residual star-forming activity (if any), resolved distribution of stellar metallicity and age, stellar and gas kinematics, and deeper CO/[C I]/[C II] observations to constrain gas-to-dust ratio will give us more hints on the origin of the dust observed in the dust-rich quiescent galaxies and quenching mechanisms in the early universe.Finding a non-negligible number of dust-rich galaxies (4/9) in the literature also pushes us to build a large sample of such, which will further help us understand the formation of these galaxies in the broader context of galaxy evolutionary studies and quenching in high redshift.
The deep ALMA Band 4 integration on a single field with a wealth of data allowed us to uncover a dust-rich quiescent candidate at  = 2.This offered a new and complementary window to study high redshift quiescent galaxies and understand the connection between the available gas/dust and galaxy quenching.Based on this, we envisage deep ALMA observations of individual quiescent galaxies at high redshift (in a few selected areas of the sky) will also provide another insightful guidance to galaxy quenching at high redshift, together with statistical studies of the quiescent population from a large survey.

APPENDIX A: PHOTOMETRIC REDSHIFT
In addition to eazy-py run, we explore two other different SED fitting codes, bagpipes (Carnall et al. 2018) and cigale (Boquien et al. 2019), to cross-check the photometric constraints of G4.The free parameters and assumptions except for SFH are the same as described in Section 3.3.For Bagpipes, we choose to use a doublepower law for simultaneously fitting the photo-, because the fit with delayed- and  model were not well converged with high  2 (≫ 30).The inspection of the corner plots of the fitted parameters also convinced us to adopt the assumption of SFH to double-power law.
For bagpipes, we take Gaussian and uniform priors to fit the redshift.For the Gaussian priors, the center is set at the best photo- from eazy-py with relatively broad width (redshift_prior_sigma ( zprior ) = 2.0).For cigale, we put uniform priors of redshift between 0 and 6.
Table A1 shows the best fit without fixing the redshift using bagpipes and cigale for G4.All SED codes give a redshift solution of  ∼ 2 for G4.The uniform prior from bagpipes did not provide a good convergence giving higher uncertainty of the   colors.We also note that at any redshift the estimated star-formation rate from bagpipes is close to zero, making the galaxy quiescent.Finally, Figure A1 shows the best fit listed in Table 3 fixing the redshift at  = 2.13 to illustrate that all SED codes fit the opt/NIR data points well.

Figure 1 .
Figure 1.The ALMA Band 4 continuum image near BX610 where we detect six continuum sources at 2 mm.In the left panel, we show the scale bar to indicate 100 kpc (physical) at  = 2.2, which corresponds to the redshift of BX610 and G1.The right panel size for each galaxy is 2 ′′ ×2 ′′ .The contours are shown from 3 in steps of 2, and the beam (0. ′′ 22×0.′′ 21) is shown on the bottom left for the natural weighting CLEAN map.-3 contours would be shown as a dashed line, which we do not find around the sources.
continuum flux is measured using the aperture of 1 ′′ .00 in the natural weighting map either based on the aperture photometry and 2D Gaussian fitting using imfit function in CASA.b: For 1.2 mm, we used the ALMA archival data (project code: 2015.1.00250.S) with 0 ′′ .32×0 ′′ .29 resolution for G4 and G3.For Q2343-BX610, we used a tapered map at 0 ′′ .40×0 ′′ .37 based on the investigation of the growth curve.c: G1's redshift is based on CO (4 − 3)and [C I], detections and more details of line emitting sources will be presented in M. Lee et al. (in preparation).d: The  = 2.21 solution is not rejected.See Section 3.2 and Figure 3.

Figure 2 .
Figure 2. Cutout images of HST/F814W, HST/F110W, HST/F140W and Spitzer/IRAC/ch2 and false color images (red: F140W, green: F110W, blue: F814W) for all galaxies detected in ALMA 2 mm continuum.The centers of individual images are all at the 2 mm peaks.The Spitzer image is zoomed out to afford the larger pixel size.A scale bar indicating 5 kpc at  = 2.21 is shown.Overlaid contours are ALMA 2 mm continuum detection (starting from 3 in steps of 2) and the bottom left shows the beam size of the ALMA image (0 ′′ .22×0′′ .21).G4 is highlighted in a red box.

Figure 3 .
Figure 3. Photometric redshift estimation from SED fitting for G4.The best-fit SED template fit from eazy-py is shown on its left panel, and the best-fit redshift is indicated as green lines.The grey line is the best fit at a fixed redshift of  = 2.21, the spectroscopic redshift of BX610 and G1.Black-filled circles show the observed flux and open squares show the best-fit model convolved with the filter response.On the right panel,  () distribution of the eazy-py (yellow shaded) and bagpieps (dark blue, uniform redshift prior), and cigale (red shaded) are shown.See also Table A1 in Appendix A.

Figure 5 .
Figure 5.The rest-frame   colors of G4 based on the bagpipes (non-parametric SFH) best-fit (star symbol), color-coded by the massweighted age obtained from the SED fitting.The overlaying data points are log(  ★ / ⊙ ) > 10.5 galaxies at 2.0 <  < 2.5 in the COSMOS2020 catalog (Weaver et al. 2022) and their number density distributions in contours.eazy-py best-fit   colors of G4 are shown as a square for a fair comparison with COSMOS2020.We define a quiescent region based on Whitaker et al. (2011).A thick solid diagonal line represents the definition by Whitaker et al. (2011) with additional cut in the  −  and  −  colors shown as a dashed line.COSOMOS2020 galaxies in the quiescent region are color-coded with their expected stellar age based on the empirical age-color relation derived by Belli et al. (2019) (for age > 300 Myr).This shows a good match between age from the empirical relation and the one from the SED fitting for G4, giving a (mass-weighted) stellar age of 1 − 2 Gyr.A region corresponding to the stellar age between 300 and 800 Myr where most post-starburst galaxies would fall, is marked according toBelli et al. (2019).Two other different definitions of the quiescent region in the   color space are plotted: the dashed-dotted line is the definition suggested byWilliams et al. (2009), and the red solid line represents the criteria byLeja et al. (2019b) for log( /yr −1 ) = −10.5.The dust attenuation direction and the amount are indicated as the arrow for   = 0.4.We show two simple toy models of   color trajectory assuming the -model SFH, for different -folding times, 100 Myr (dark blue solid, fast track) and 1 Gyr (dark red dotted, slow track) using fsps-py followingBelli et al. (2019) (see Section 5.3 for details).

Figure 6 .
Figure6.Constraints on dust-to-stellar mass ratios as a function of redshift for  > 1 quiescent galaxies.G4 is indicated by a star symbol.Sources from the literature are indicated as circles if they are detected in dust emission; otherwise, inverted triangles indicate upper limits (3).ALMA.14 and ZF-COS-19589, two galaxies with high dust-to-stellar mass ratios comparable to G4, are also labeled.The color of each source indicates the stellar mass of the galaxy.The stacking results byMagdis et al. (2021) for massive (log(  ★ / ⊙ ) > 10.8) quiescent galaxies at  ≲ 2 are shown as grey squares.Grey shaded region with a dotted curve is based on the functional form for the best-fit at  < 1,  gas (=  gas / ★ ) ∝ (1 + ) 5 , assuming a fixed gas-to-dust mass ratio of 92 and normalized at  = 1 with  gas = 8%.For comparison to star-forming galaxies of similar mass, we calculated dust-to-stellar mass ratios using the scaling relation suggested byTacconi et al. (2018), assuming a gas-to-dust mass ratio (GDR) of 100 for log (  ★ / ⊙ ) = 11.0 and 10.5 (red and blue shaded bands, respectively; more massive galaxies have smaller gas (and dust) at fixed GDR.

Figure 7 .
Figure 7. Age (mass-weighted) distribution as a function of the offset from the main sequence.G4 is plotted as star symbols, using different SFR estimates (24m flux scaling and SED); the fainter one is the best-fit SED using bagpipes (non-parametric SFH) and see Section 4.1 for detailed discussions.We use the main-sequence definition from Speagle et al. (2014).Literature values are from Schreiber et al. (2018); Belli et al. (2019); Estrada-Carpenter et al. (2020); Stockmann et al. (2020); D'Eugenio et al. (2020); Valentino et al. (2020); Kubo et al. (2021) for  ∼ 1 − 4 galaxies with constraints of the mass-weighted age, stellar mass, and star-formation rate.Here, we use SFR values from the SED fitting, except for two cases (ZF-COS-19589, D'Eugenio et al. 2020); for ZF-COS-19589, the 870m flux is used to calculate SFR(IR) assuming  d = 40 K taken from Suzuki et al. 2021; D'Eugenio et al. (2020) performed stacking analysis and we took the dust-corrected SFR from the (stacked) [OII] line flux.The thick blue dashed line shows our best fit of age-ΔMS relation for  ∼ 1 quiescent galaxies using Belli et al. (2019) and Estrada-Carpenter et al. (2020) good quality samples (see Section 5.1).Highz ( > 1) quiescent galaxies with dust/CO/[C I] observations are also plotted, taken from the various literature (see Section 5.2 and TableB1).

Figure 9 .
Figure9.Dust-to-stellar mass ratio as a function of mass-weighted age.G4 is plotted as a star color-coded with its sSFR.We also plot other  > 1 quiescent galaxies with dust constraints from the literature (detection for circles, 3 upper limit with upside-down triangles for non-detection;Gobat et al. 2022;Whitaker et al. 2021;Suzuki et al. 2022;Kalita et al. 2021).The three lines are inspired byMichałowski et al. (2019): dashed line is the fit based on the local Herschel selected early-type galaxies, while the other two lines are tweaked and rescaled assuming  = 0.7 Gyr (dashed-dotted) and 2.5 Gyr (dotted) to guide our eye for dust-poor high-redshift quiescent galaxies.The colors represent the specific star-formation rate of the galaxies.

Figure 10 .
Figure10.Dust yield required for SNII and AGB stars to explain the observed dust in high redshift ( > 1) quiescent galaxies including G4.The dust yield of SN II is in yellow star and for AGB in magenta circle for dust detected sources.Upper limits are shown as inverted triangles of each color.Two dashed lines indicate the maximum dust yield range observed in SN 1987A and Cas A  ⊙ ).The gray shaded region indicates the range of the dust yield expected from theory for SN II without destruction, while the blue area represents a more realistic yield range due to dust destruction (≪ 0.1  ⊙ ).The pink line indicates the maximum theoretical yield possible for the AGB stars (0.04  ⊙ ) but the observational results claim yields lower by more than an order of magnitude (see details in the main text).G4 requires almost the maximum dust yields from both SN II and AGB stars if grain growth is not taken into account.

Table 1 .
A summary of ALMA 2 mm detected sources

Table 2 .
Opt/NIR Photometric data points for G4 Smith et al. 2012imitation of the SFR from the SED which is likely to introduce an unrealistically low value of SFR and hence ΔMS especially for the parametric SFH models.The photometry (R  -band) and the best-fit dust attenuation ( V ) suggest the SFR of ∼ 8  ⊙ yr −1 and log(ΔMS) ∼ −1.6 (see also the main text in Section 4.1).‡:We quote the absolute  2 value, as the templates employed (as is the case for many SED fitting codes) are not independent of each other, and degrees of freedom are ill-defined (e.g.,Smith et al. 2012). : best-fit photo- from eazy-py. : mass-weighted. : averaged over 100 Myr  : including the ALMA 2 mm data point.

Table 4 .
Brisbin et al. (2019)8)and 24 m data points  d = 17-35 K for modified black body and usingKennicutt (1998)for  IR to SFR conversion.We note that the 1.2 mm data point gives a range of 2 − 72  ⊙ yr −1 , but the uncertainty is much larger if we take into account the photometry uncertainty.Based on the 24 m flux scaling between BX610 and G4 and using the SFR of BX610 from using FörsterSchreiber et al. (2018)(115  ⊙ yr −1 ) andBrisbin et al. (2019)(140  ⊙ yr −1 ).
: Assuming  = 1.5 − 2.2,  : Based onReddy et al. (2010)and the values in the parenthesis are taking the scatter of the relation between   and L(H) (and SFR), not the 24m flux errors. :