High Molecular-Gas to Dust Mass Ratios Predicted in Most Quiescent Galaxies

Observations of cold molecular gas reservoirs are critical for understanding the shutdown of star formation in massive galaxies. While dust continuum is an efficient and affordable tracer, this method relies upon the assumption of a"normal"molecular-gas to dust mass ratio, $\delta_{\mathrm{GDR}}$, typically of order one hundred. Recent null detections of quiescent galaxies in deep dust continuum observations support a picture where the cold gas and dust has been rapidly depleted or expelled. In this work, we present another viable explanation: a significant fraction of galaxies with low star formation per unit stellar mass are predicted to have extreme $\delta_{\mathrm{GDR}}$ ratios. We show that simulated massive quiescent galaxies at $0$4 orders of magnitude. The dust in most simulated quiescent galaxies is destroyed significantly more rapidly than the molecular gas depletes, and cannot be replenished. The transition from star-forming to quiescent halts dust formation via star formation processes, with dust subsequently destroyed by supernova shocks and thermal sputtering of dust grains embedded in hot plasma. After this point, the dust growth rate in the models is not sufficient to overcome the loss of $>$3 orders of magnitude in dust mass to return to normal values of $\delta_{\mathrm{GDR}}$ despite having high metallicity. Our results indicate that it is not straight forward to use a single observational indicator to robustly pre-select exotic versus normal ratios. These simulations make strong predictions that can be tested with millimeter facilities.


INTRODUCTION
After over a decade and thousands of hours of observations, we now have a reasonably good census of the dust and cold gas content in 'normal' star forming galaxies out to z ∼ 2, as well as other galaxy properties that correlate (see reviews by, e.g., Carilli & Walter 2013;Hodge & da Cunha 2020;Tacconi et al. 2020). However, the dust and cold gas content of quiescent galaxies remain uncertain, especially towards higher redshift. Spectroscopic studies measuring the cold gas content of quiescent galaxies at z > 0.5 are limited, as robust constraints are observationally expensive. First results from intermediate redshift studies using CO emission to trace metal-rich molecular hydrogen are conflicting: some yield surprisingly large gas reservoirs (e.g., Suess et al. 2017;Rudnick et al. 2017;Hayashi et al. 2018;Belli et al. 2021), while others result in strong upper limits suggesting rapid gas depletion (e.g., Sargent et al. 2015;Bezanson et al. 2019;Williams et al. 2021), and one study results in a mixture (Spilker et al. 2018). The diversity in M H2 may be due to the range in specific star formation rates (sSFR≡SFR/M ) of the samples themselves, as low sSFRs are hard to constrain (e.g., Leja et al. 2019).
By stacking low-resolution far-infrared to submillimeter imaging, the first studies of statistically meaningful samples of quiescent galaxies at z ∼ 1 − 2 find moderate dust and inferred gas content of order f H2 = M H2 /M ∼ 5 − 10% (Gobat et al. 2018;Magdis et al. 2021). These results are within 2σ of an extrapolation of the Tacconi et al. (2018) scaling relations of sSFR and f H2 . However, deblending low-resolution data is non-trivial and systematic uncertainties can be large (Viero et al. 2012). While these analyses perform a comprehensive modeling of the full infrared spectral energy distribution (SED), local studies often leverage correlations between the total gas mass and a single far-infrared band. That said, debate remains regarding which wavelengths provide the strongest correlation. Some studies find the highest correlation of the molecular-gas mass with dust mass as traced by the longest wavelengths * NHFP Hubble Fellow (e.g., Bourne et al. 2013), whereas others find that the molecular-gas mass is best correlated with the obscured SFR at the peak of the infrared SED (rest-frame 100-160µm; Groves et al. 2015). Moreover, there is not yet consensus whether the gas traced by the dust continuum is primarily molecular or a combination with neutral hydrogen (Janowiecki et al. 2018). Scoville et al. (2016) calibrate the Rayleigh-Jeans dust continuum as a tracer of molecular-gas mass, using a sample of star-forming galaxies at z ∼ 0 and z ∼ 2 with both CO(1-0) and dust continuum fluxes. Combining this methodology based on optically-thin dust continuum with sensitive millimeter telescopes may be one of the most efficient ways to observe gas in distant quiescent galaxies. Uncertainties due to variations in dust temperature are thought to be minimal at λ rest > 250 µm (e.g. Scoville et al. 2014), though some tension exists at high redshift (Harrington et al. 2021). The major driver of the scatter in the conversion of dust continuum to gas mass is instead the ratio of the molecular-gas mass to dust mass, δ GDR =M H2 /M dust (Privon et al. 2018). Dust continuum is demonstrated to be a robust tracer of the molecular-gas mass for 'normal' star-forming galaxies (Privon et al. 2018;Liang et al. 2018;Kaasinen et al. 2019), but calibrations are not yet tested for quiescent galaxies. Inherent to the Scoville et al. (2016) methodology is the assumption that δ GDR =150, though it is bundled together with the assumed dust emissivity per unit mass. High δ GDR ratios in nearby quiescent galaxies are thought to be the result of thermal sputtering resulting from the impact of dust grains with hot plasma (e.g., Galliano et al. 2018;Smercina et al. 2018).
Even when leveraging strong gravitational lensing, Whitaker et al. (2021) detect cold dust in only two out of a sample of six quiescent galaxies at z ∼ 2. When assuming δ GDR of 100, these data impose strong upper limits of f H2 < 1%. The most obvious conclusion from the null detections is that the galaxies harbor very small dust and molecular-gas reservoirs. This is consistent with other studies using independent tracers of molecular hydrogen that find similarly low gas fractions among a sample of massive quiescent galaxies (e.g., Williams et al. 2021). However, sample sizes remain small and we cannot rule out that the assumed δ GDR ratio is actu- Figure 1. Dust masses are dramatically lower than the molecular-gas mass for model galaxies with low sSFRs (red), assuming the standard relative abundance of δGDR=MH2/M dust =100 (solid line). The trend becomes most pronounced by z=1, as the quiescent population grows. The models predict that dust is rapidly destroyed when star formation tapers.
ally much higher at low sSFR, implying under-predicted inferred gas fractions.
In this paper, we seek to gain physical intuition on δ GDR ratios at low sSFR by studying populations of simulated galaxies ranging from z=0 to z=3 in the simba galaxy formation simulation that separately tracks H 2 and dust. We present details of the simulated galaxy sample in §2, and show results in §3 including a comparison to the literature. In §4, we discuss the implications of a theoretical model tuned to reproduce dust-to-metal ratios at z=0 in the context of the existing astronomical literature. simba adopts a ΛCDM cosmology with Ω M = 0.3 and H 0 = 68 km s −1 Mpc −1 , and a Chabrier (2003) initial mass function.

NUMERICAL METHODS
We turn to hydrodynamic cosmological simulations that include predictive modeling of the gas and dust physics to better understand this problem. Li et al. (2019) present a self-consistent model for the dust-togas and dust-to-metal ratios in galaxies in the simba cosmological hydrodynamic galaxy formation simulation , tracking the production, growth, and destruction of dust grains over time, broadly following McKinnon et al. (2017) with notable improvements (Li et al. 2019). The dust content is primarily governed by (1) formation of dust in stellar ejecta (e.g., Type II supernovae and asymptotic giant branch stars), (2) growth of dust via the accretion of metals, and (3) destruction of dust via thermal sputtering, consumption by star formation, or supernovae (SNe) shocks. Molecular-gas is computed within dense (n H > 0.13 cm −3 ) gas using a metallicity-dependent prescription following Krumholz & Gnedin (2011). We refer to Davé et al. (2019), Li et al. (2019), and Narayanan et al. (2020) for full details on Simba.
Herein, we define f H2 =M H2 /M from the simulations and ignore neutral hydrogen. The Scoville et al. (2016) methodology and CO measurements (Bolatto et al. 2013) are both calibrated to trace molecular hydrogen alone, under the assumption that the fraction of dust associated with molecular hydrogen is greater than that of neutral hydrogen. This is not necessarily true for quiescent galaxies; if dust traces both molecular and neutral hydrogen, with significant neutral reservoirs expected (e.g., Zhang et al. 2019), dust continuum serves as upper limits for M H2 . However, defining the gas mass to include neutral hydrogen (i.e., M gas =M HI +M H2 ) only exaggerates the implied exotic values of δ GDR . Thus, the molecular hydrogen is summed for all particles within 30 kpc physical spherical aperture. Dust is fully coupled to the gas and its mass is calculated in the same way. This physical aperture is important for removing contributions from nearby small, unresolved galaxies, which can Figure 2. Histograms of the molecular-gas to dust ratio, δGDR, for simba model galaxies from z=0 (top left) to z=3 (bottom right). Star-forming galaxies peak at δGDR ∼ 100, whereas we see the buildup of the quiescent population spanning ∼4 orders of magnitude and peaking at δGDR ∼ 10 4 − 10 5 . Once the dust is destroyed, rarely can a galaxy regrow it. Galaxy populations are separated into star-forming and quiescent based on their distance from the log(SFR)log(M ) relation. be significant for a quenched galaxy that has very little molecular-gas itself.

MOLECULAR-GAS TO DUST MASS RATIOS IN THE SIMBA COSMOLOGICAL SIMULATION
The relation between the molecular-gas mass and dust mass of simba model galaxy populations ranging from z=0 (left) to z=3 (right) is found in Figure 1. The solid line is the standard value of δ GDR =100, with dotted lines offset in increments of 1 dex. Quiescent galaxies are broadly identifiable on the red end of the color spectrum. We show that star-forming galaxies (blue points) follow a relatively tight relation between dust mass and molecular-gas mass consistent with δ GDR ∼ 100, but we see a large intrinsic scatter among the low sSFR population. In particular, δ GDR for low sSFR galaxies spans six orders of magnitude and most objects have significantly higher ratios than standard assumptions of δ GDR ∼100-200 (e.g., Tacconi et al. 2020). This dramatic change in δ GDR is the result of rapidly decreasing dust masses. Figure 2 shows the redshift evolution of δ GDR when separating star-forming (blue) and quiescent (red) galaxies into bimodal galaxy populations. In order to do this in a meaningful way, we match the log(SFR)log(M ) relations parameterized by Whitaker et al. (2014) to the mean ridge-line of the model galaxies. The polynomial fit to the observations at z=1-1.5, z=1.5-2 and z=2-2.5 are matched to running mean of the models at z=1, 2, and 3, respectively. Shifting the z=0.5-1 fit down by 0.5 dex aligns with the model relation at z=0. These shifts remove known offsets between observed UV+IR SFRs and simulated datasets (Nelson et al. 2021). Our goal is to identify galaxies at a given epoch that are >0.7 dex below the star-forming sequence, labeling this population as quiescent.
Through this exercise, we learn something important: star-forming galaxies consistently peak at δ GDR ∼ 100 with σ=0.4 dex, whereas the quiescent population has a wider range of δ GDR . This conclusion is robust to perturbations of the star-forming/quiescent boundary. In the models, the dust mass is governed by dust formation (Type II SNe, AGB stars), dust growth (accretion of metals), and dust destruction by thermal sputtering and SNe shocks. Here, thermal sputtering in the hot ISM and SNe shocks both contribute to the dust destruction process 1 . To verify this, we conduct a series of controlled numerical experiments in which we run comparable resolution (but smaller box) (25/h Mpc) 3 simulations with 256 3 particles, and systematically turn off dust growth, all destruction processes (except by star formation itself -this is not posssible as the dust is tied to gas particles in simba), and just thermal sputtering (see Figure 3). If we turn off dust growth altogether in the models, this has the most dramatic effect on the peak location of δ GDR for star-forming galaxies, shifting 2 dex higher. Turning off dust destruction mechanisms does not completely remove the bimodality, indicating that dust destruction in star-forming regions also plays an important role. A key aspect here is that once the dust is destroyed, rarely can the galaxy regrow it, as indicated by the build-up of the quiescent population at δ GDR ∼ 10 4 − 10 5 . This bimodality within the quiescent population at z = 0 is partially driven by the star-forming/quiescent definition, where some ambiguity exists.
We color-code model galaxies by the median δ GDR value in Figure 4 to demonstrate the impact of the large intrinsic scatter in the factor transcribing f dust (bottom) into f H2 (top). This wide range in δ GDR produces a moderately tight relation between f H2 and log(sSFR), despite the variations observed in f dust . Again, this is because the intrinsic scatter of δ GDR is driven by variations in M dust . Figure 4 includes a compilation of CO Figure 3. Comparison of (25/h Mpc) 3 box simulation runs for the full model relative to (a) no dust growth, (b) turning off dust destruction (labeled 'No Destruction') by SN shocks and thermal sputtering, but not astration -this is not possible to turn off, and (c) no thermal sputtering. The top panels demonstrate the redshift evolution of the median molecular-gas to dust mass ratio (left) and dust mass (right) for star-forming (blue) and quiescent (red) galaxies. Only bins with greater than 10 galaxies are shown. The bottom panels show a comparison of the full model distribution at z = 1 (gray, no separation of quiescent/star-forming) relative to each of the cases separated into quiescent/star-forming populations. Dotted lines represent the mode.
The lack of dust continuum detections at extremely low dust fractions may not be physical, but rather a selection bias. While the local early-type galaxies and red spirals from Rowlands et al. (2012) (compilation adopted from Micha lowski et al. 2019) shown in Figure 4 appear to have significantly larger dust masses than predicted by the models, this sample is blindly selected at sub-millimeter wavelengths and thus may not be representative. Values of f dust as low as 10 −5 -10 −6 have been observed in nearby early-type galaxies (Smith et al. 2012). However, placing such low limits on dust mass becomes technically challenging, even in the nearby universe.
If we take these predicted exotic molecular-gas to dust mass ratios at face value, even with extremely sensitive limits of millimeter observations, there may be no hope to detect the lowest sSFR galaxies in dust continuum in the near term. In some simulated metal-rich quiescent galaxies, dust regrowth is not enough to significantly increase the dust mass to counteract the destruction processes. Such model predictions have important implications which we explore in the following section.

DISCUSSION
In this letter, we present predictions from the hydrodynamic cosmological simba simulations showing there is a dramatic increase in the molecular-gas mass to dust mass ratios, δ GDR , when star formation slows and falls below log(sSFR) -10 yr −1 . The scatter in δ GDR spans >4 orders of magnitude, with some quiescent galaxies having "normal" ratios whereas most others have exotic ratios.
The wide range of predicted δ GDR values within the simulations cannot be explained by variations in metallicity as almost all model galaxies at low sSFR have so- The models predict a tight correlation between decreasing fH2 with decreased sSFR. The squares represent the median molecular-gas to dust ratio, with small circles representing individual model galaxies plotted in sparsely populated parameter space. Black symbols represent CO measurements of quiescent galaxies in the literature, with upper limits indicated with an arrow (or left caret at z=0). Bottom: f dust fans out to cover up to 5 orders of magnitude below log(sSFR)<-10 yr −1 , given the predicted high molecular-gas to dust mass ratios. Black symbols represent dust continuum measurements of quiescent galaxies in the literature. lar or super-solar metallicities. Instead, it is likely that there exists rapid dust destruction in some galaxies shutting down star formation -beyond dust consumed in the major episode of star formation itself. The simba simulation includes three processes that destroy dust in low sSFR galaxies: (1) thermal sputtering by hot electrons, and destruction by (2) astration (the incorporation of dust into a stellar interior -a star particle here -during star formation) prior to the quenching (which is mainly due to supermassive black hole feedback) combined with (3) unresolved SN shocks. All three play a role here (e.g., Figure 3), driving a rapid decrease in the predicted dust masses when star formation quenches. Thermal sputtering in particular ramps up as gas density and/or temperature increases, and is invoked to explain high δ GDR values in nearby post-starburst galaxies (e.g., Smercina et al. 2018).
Interestingly, the first cold dust continuum detections of low sSFR galaxies outside the local universe remain unresolved despite the extreme extended stellar light profiles resulting from strong gravitational lensing . Though only a sample of two, these data suggest high molecular gas and dust surface densities. This has also been found in star-forming galaxies at similar redshifts (e.g., Tadaki et al. 2017) and nearby post-starburst galaxies (Smercina et al. 2021). Star formation may be suppressed by significant turbulent heating under these conditions, as internal turbulent pressure is proportional to gas surface density (Smercina et al. 2021). Unfortunately, the simulations cannot resolve effects down to physical scales of the individual clouds. Empirical confirmation of moleculargas surface densities in significant samples of quiescent galaxies would corroborate this idea.
Chemical evolution models presented in Rémy-Ruyer et al. (2014) predict that galaxies with the shortest star formation timescales of 0.5 Gyr sustain high δ GDR ratios for the longest time periods, up until they reach high metallicity, at which point δ GDR drops to more typical values on average. The majority of quiescent galaxies in the simba simulation do not experience such a latestage drop in δ GDR (e.g., Figure 2) despite having high metallicity. This effect in the Rémy-Ruyer et al. (2014) models is the result of dust growth by the accretion of metals, but the dust mass in simba is too small to recover from depletion.
It is also possible for dust production to increase once again as the stellar populations age. A common process to replenish dust reservoirs is via AGB stars. This phase of stellar evolution occurs when stars are around Figure 5. At z=2, there exists a large scatter in f dust at light-weighted ages older than 500 Myr. More observations of dust in quiescent galaxies, even deeper than those augmented by strong lensing (black Caliendo et al. 2021;Whitaker et al. 2021), are required to test the predicted exotic molecular-gas to dust ratios and the steep turnover for light-weighted ages >500 Myr.
∼ 1 Gyr, similar to the ages of quiescent galaxies at cosmic noon. AGB stars produce copious amounts of dust that could replenish dust reservoirs. However, nearby studies of quiescent galaxies find that f dust (Micha lowski et al. 2019) and f H2 (French et al. 2015) both correlate with the age of the stellar populations, with older galaxies having lower dust and molecular-gas masses. In Figure 5, we consider the trend within the simulations at z=2 between f dust and the light-weighted age, colorcoded by δ GDR . The predicted trend in the high redshift universe is not gradual. There is a dramatic turnover for ages >500 Myr. We defer a more detailed analysis using star formation histories to discern between typical and exotic δ GDR in quenched galaxies to future work.
The exotic molecular-gas to dust mass ratios predicted by simba call into question the efficiency with which one can use the observed dust continuum as an indirect measure of the interstellar medium at low sSFRs. To develop a more complete understanding of the efficacy of this technique, we must secure detections of both CO emission and dust continuum (ideally in multiple bands) for the same galaxies. Uncertainties in the temperature and composition of dust in quiescent galaxies introduce important systematics when estimating the moleculargas mass from a single dust continuum measurement (Liang et al. 2019;Privon et al. 2018). Whereas such observations are costly with our current generation of telescopes, the boost from gravitational lensing for spec-tacular sources offers a pathway forward in the near term for small samples. Upcoming upgrades to existing facilities will also be fruitful: the improved sensitivity and wide field of view of the Toltec instrument on the Large Millimeter Telescope will enable similarly sensitive dust continuum studies of mass-representative quiescent samples at high redshift, and the spectroscopic capabilities predicted for the next generation Very Large Array will explore new parameter space for CO in these sources. Until then, the simulations paint a grim picture for the future utility of dust continuum for most low sSFR galaxies.
Predictions from simba encompass an extreme range of δ GDR ratios that are not yet observed. Owing to the scatter in these predicted ratios, the models do not elucidate a clear physical observable (i.e., stellar mass, surface density, sSFR, etc) to distinguish the expected δ GDR (Li et al. 2019). But as a whole these predictions can be empirically tested and validated; future deeper CO observations and/or larger statistical dust continuum samples will help improve our understanding of the redshift evolution and behavior of δ GDR at low sSFR. While current dust continuum sensitivity limits preclude confirmation of the low dust fractions predicted in Figure 4, deeper CO measurements have the power to rule out these exotic ratios. By also measuring the spatial distributions of the molecular-gas, we can confirm the compact nature. If the model outcomes cannot be verified, the assumptions about dust processes must be revisited. There exists a beautiful synergy between these future observations and the model outcomes that will guide both observers and theorists alike.