The Dust-to-Gas and Dust-to-Metals Ratio in Galaxies from z=0-6

TL;DR

This study models the evolution of dust-to-gas and dust-to-metal ratios in galaxies from redshift 0 to 6, using simulations and machine learning to understand their relationships with galaxy properties.

Contribution

It introduces a comprehensive dust evolution model integrated into cosmological simulations and provides a publicly available code for estimating dust ratios based on galaxy properties.

Findings

The simulation reproduces the local dust mass function but underestimates high-redshift values.

DGR correlates strongly with gas-phase metallicity, weakly with other properties.

Little evolution in DGR-metallicity relationship from z=0 to 6.

Abstract

We present predictions for the evolution of the galaxy dust-to-gas (DGR) and dust-to-metal (DTM) ratios from z=0 to 6, using a model for the production, growth, and destruction of dust grains implemented into the \simba\ cosmological hydrodynamic galaxy formation simulation. In our model, dust forms in stellar ejecta, grows by the accretion of metals, and is destroyed by thermal sputtering and supernovae. Our simulation reproduces the observed dust mass function at z=0, but modestly under-predicts the mass function by ~x3 at z ~ 1-2. The z=0 DGR vs metallicity relationship shows a tight positive correlation for star-forming galaxies, while it is uncorrelated for quenched systems. There is little evolution in the DGR-metallicity relationship between z=0-6. We use machine learning techniques to search for the galaxy physical properties that best correlate with the DGR and DTM. We find that the DGR is primarily correlated with the gas-phase metallicity, though correlations with the depletion timescale, stellar mass and gas fraction are non-negligible. We provide a crude fitting relationship for DGR and DTM vs. the gas-phase metallicity, along with a public code package that estimates the DGR and DTM given a set of galaxy physical properties.

Tables

Table 1. Table 1: Simulation Free Parameters

Parameter	Description	Value	Range accepted by literatures
Thermal sputtering
$a$	Grain radius ($\mu m$)	0.1	—
$\sigma$	Density of solid matters within grains (g cm-3)	2.4	2.2(graphite), 3.3(silicate)a
Production
${\delta }_{i,\mathrm{dust}}^{\mathrm{AGB},C/O>1}$	Condensation efficiency	0.2 for $i$ = C	0.2 – 1.0b
		0 otherwise	0
		0 for $i$ = O	0
		0.2 otherwise	0.2 – 0.8b
${\delta }_{i,\mathrm{dust}}^{\mathrm{SNII}}$		0.15 for $i$ = C	0.15 – 1.0b
		0.15 otherwise	0.15 – 0.8b
Growth
${\rho }^{\mathrm{ref}}$	Reference density (g cm-3)	$2.3\times {10}^{-22}$	—
$T^{\mathrm{ref}}$	Reference temperature (K)	20	—
${\tau }_g^{\mathrm{ref}}$	Growth time-scale with $T=T^{\mathrm{ref}}$ and $\rho ={\rho }^{\mathrm{ref}}$ (Myr)	10	2 – 500c
Destruction (SNe Shock)
$E_{\mathrm{SN},51}$	Energy per SN (${10}^{51}$ erg)	1.0	—
$ϵ$	The efficiency of destruction by SN shocks	0.3	0.1 – 0.5d

Equations

$$m_{i,d}^{\mathrm{AGB}}=\{\begin{array}{ll}{\delta }_{\mathrm{C}}^{\mathrm{AGB}}(m_{C,\mathrm{ej}}^{\mathrm{AGB}}-0.75m_{O,\mathrm{ej}}^{\mathrm{AGB}}),& i=\mathrm{C}\\ 0,& \mathrm{otherwise},\end{array}$$

$$m_{i,d}^{\mathrm{AGB}}=\{\begin{array}{ll}0,& i=\mathrm{C}\\ 16\sum _{i=\mathrm{Mg},\mathrm{Si},\mathrm{S},\mathrm{Ca},\mathrm{Fe}}{\delta }_i^{\mathrm{AGB}}{m}_{i,\mathrm{ej}}^{\mathrm{AGB}},& i=\mathrm{O}\\ {\delta }_i^{\mathrm{AGB}}{m}_{i,\mathrm{ej}}^{\mathrm{AGB}},& \mathrm{otherwise},\end{array}$$

$$m_{i,d}^{\mathrm{SNII}}=\{\begin{array}{ll}{\delta }_{\mathrm{C}}^{\mathrm{SNII}}{m}_{\mathrm{C},\mathrm{ej}}^{\mathrm{SNII}},& i=\mathrm{C}\\ 16\sum _{i=\mathrm{Mg},\mathrm{Si},\mathrm{S},\mathrm{Ca},\mathrm{Fe}}{\delta }_i^{\mathrm{SNII}}{m}_{i,\mathrm{ej}}^{\mathrm{SNII}},& i=\mathrm{O}\\ {\delta }_i^{\mathrm{SNII}}{m}_{i,\mathrm{ej}}^{\mathrm{SNII}},& \mathrm{otherwise},\end{array}$$

$${\left(\frac{\mathrm{d}{M}_d}{\mathrm{d}t}\right)}_{\mathrm{grow}}=\left(1-\frac{M_d}{M_{\mathrm{metal}}}\right)\left(\frac{M_d}{{\tau }_{\mathrm{accr}}}\right),$$

$${\tau }_{\mathrm{accr}}={\tau }_{\mathrm{ref}}\left(\frac{{\rho }_{\mathrm{ref}}}{{\rho }_g}\right)\left(\frac{T_{\mathrm{ref}}}{T_g}\right)\left(\frac{Z_{\odot }}{Z_g}\right).$$

$${\tau }_{\mathrm{sp}}$$

$$\sim (0.17\mathrm{Gyr})\left(\frac{a}{0.1\mu m}\right)\left(\frac{10^{-27}\mathrm{g}\mathrm{c}{\mathrm{m}}^{-3}}{{\rho }_g}\right)\left[{\left(\frac{T_0}{T_g}\right)}^{\omega }+1\right],$$

$${\left(\frac{\mathrm{d}{M}_d}{\mathrm{d}t}\right)}_{\mathrm{sp}}=-\frac{M_d}{{\tau }_{\mathrm{sp}}/3}$$

$${\tau }_{\mathrm{de}}=\frac{M_g}{ϵ\gamma M_s},$$

$$\log \text{DGR}=(2.445\pm 0.006)\log \left(\frac{Z}{Z_{\odot }}\right)-(2.029\pm 0.003).$$

Full text

††pubyear: 2019††pagerange: The Dust-to-Gas and Dust-to-Metal Ratio in Galaxies from $z=0-6$–References

1 Introduction

Dust plays a critical role in the physics of the interstellar medium (ISM) and galaxy evolution. The surfaces of dust grains catalyze a range of chemical reactions that influence the structure of ISM and star formation (Hollenbach & Salpeter, 1971; Mathis, 1990; Weingartner & Draine, 2001; Draine, 2003; Wolfire et al., 2008; Hollenbach et al., 2012; Gong et al., 2017), including the formation of molecular hydrogen and grain-catalyzed recombinations of H+ and C+. The ejection of dust from galaxies can contribute to metal abundances in the intergalactic medium and offers an additional cooling channel (Ostriker & Silk, 1973; Bouché et al., 2007; Peeples et al., 2014; Ménard et al., 2010; Peek et al., 2015; Vogelsberger et al., 2018), while dust absorption of far ultraviolet and optical photons can shape the temperature structure of the neutral ISM (Goldsmith, 2001; Krumholz et al., 2011; Narayanan et al., 2011, 2012; Narayanan & Davé, 2012).

A complex set of physical processes contributes to the evolving dust content of the Universe. It can be produced via condensation of dust grains from the gas-phase metals in the ejecta of asymptotic giant branch (AGB) stars and supernovae (SNe; Gehrz 1989; Todini & Ferrara 2001; Nozawa et al. 2003; Ferrarotti & Gail 2006; Nozawa et al. 2007; Zhukovska et al. 2007; Nanni et al. 2013; Schneider et al. 2014), after which it can grow in the ISM via accretion of gas-phase metals (Dominik & Tielens, 1997; Dwek, 1998; Hirashita & Kuo, 2011; Zhukovska, 2014). It can be destroyed via enhanced non-thermal sputtering in SN blast waves, thermal sputtering, and via grain-grain collisions (Draine & Salpeter, 1979a, b; Seab & Shull, 1983; McKee et al., 1987; Jones et al., 1996; Bianchi & Ferrara, 2005; Nozawa et al., 2007).

Dust properties in galaxies have been intensively studied through statistics and scaling relations, of which three particularly interesting are dust mass functions (DMFs) (Dunne et al., 2003; Vlahakis et al., 2005; Eales et al., 2009; Dunne et al., 2011; Clemens et al., 2013; Beeston et al., 2018), dust-to-gas mass ratios (DGRs) and dust-to-metal mass ratios (DTMs) as a function of galaxy metallicity or stellar mass (Issa et al., 1990; Lisenfeld & Ferrara, 1998; Hirashita et al., 2002; Draine et al., 2007; Galametz et al., 2011; Rémy-Ruyer et al., 2014; Giannetti et al., 2017; Chiang et al., 2018; Kahre et al., 2018; De Vis et al., 2019; De Cia et al., 2013, 2016; Zafar & Watson, 2013; Sparre et al., 2014; Wiseman et al., 2017). These relationships provide a convenient method for determining gas masses in galaxies, as well as providing constraints on the baryon cycle that governs galaxy evolution at low and high-redshifts (e.g. Magdis et al., 2012).

Theorists have commonly used an assumed constant dust-to-metal ratio in galaxies in order to model the evolving dust content in hydrodynamic or semi-analytic models of galaxy formation (e.g. Silva et al., 1998; Granato et al., 2000; Baugh et al., 2005; Lacey et al., 2010; Narayanan et al., 2010, 2015; Narayanan et al., 2018a, b; Fontanot & Somerville, 2011; Niemi et al., 2012; Somerville et al., 2012; Hayward et al., 2013; Cowley et al., 2017; Katz et al., 2019; Ma et al., 2019). However, there is growing evidence from both integrated and resolved far-infrared studies of galaxies at both low and high-redshift that the dust-to-gas and dust-to-metal ratios in galaxies are not constant, and may not even be straight-forwardly modeled by a simple linear relationship with a galaxy physical property (such as metallicity). For example, while the gas to dust ratio appears to scale with the metallicity of galaxies in the local Universe (Dwek, 1998; Draine et al., 2007; Bendo et al., 2010), there may be deviations from this trend at the lowest metallicities (e.g. Galliano et al., 2005; Galametz et al., 2011; Rémy-Ruyer et al., 2014; De Vis et al., 2019). Similarly, the DTG measured by damped Lyman-alpha (DLA) and gamma-ray burst (GRB) absorbers (e.g. De Cia et al., 2013, 2016; Wiseman et al., 2017) from $z=0.1$ to $z=6.3$ are similar to those in the Local group, though drop at metallicities lower than 0.05 $Z_{\odot}$. Hence more sophisticated theoretical modeling of galaxy dust content and its evolution is needed.

In recent years, galaxy evolution models have progressed from treating dust as a simple scale factor of the metal mass (see Somerville et al., 2012, and references therein) to including the physics of dust formation, growth and destruction in galaxies as they evolve. The first generation of these sorts of simulations treated galaxies as one-zone models (e.g. Issa et al., 1990; Dwek, 1998; Inoue, 2003; Morgan & Edmunds, 2003; Calura et al., 2007; Zhukovska et al., 2007; Hirashita & Yan, 2009; Asano et al., 2013; Calura et al., 2014; Rowlands et al., 2014; Zhukovska, 2014; Feldmann, 2015; De Vis et al., 2017), though more recently a number of groups have begun to incorporate self-consistent dust physics on-the-fly into bona fide hydrodynamic models of galaxy formation and evolution. Bekki (2015), McKinnon et al. (2016), and Aoyama et al. (2017) established some of the initial frameworks for including dust in hydrodynamic galaxy formation simulations to study the evolution of dust properties in individual galaxies. Building on this, McKinnon et al. (2017) performed full-volume cosmological simulations using moving mesh code arepo to study dust properties across galaxies over cosmic time, but were unable to successfully reproduce the DGR – metallicity relation. Vogelsberger et al. (2018) extended this framework by implementing high temperature dust cooling channels to study dust in galaxy clusters and its impact on the intergalactic medium. Aoyama et al. (2017) and Aoyama et al. (2018) developed a $2$-grain size model into a SPH cosmological simulation, where they studied overall dust properties in a whole cosmological volume and IGM, while Hou et al. (2019) built on this to add a phenomenological Active Galactic Nuclei (AGN) feedback model. Finally, Popping et al. (2017) and Vijayan et al. (2019) have implemented the physics of dust formation, growth and destruction into semi-analytic galaxy formation models. The growing interest in modeling dust evolution highlights its importance in more accurately modeling the observed properties of galaxies.

What has been missing thus far is a predictive self-consistent model for the dust-to-gas and dust-to-metal ratios in galaxies across cosmic time in a large-volume cosmological galaxy formation simulation. In this paper, we aim to develop this model. To do this, we incorporate into state-of-the-art cosmological hydrodynamic simulation Simba (Davé et al., 2019) a model to track on-the-fly dust formation and evolution, broadly following the McKinnon et al. (2017) passive scalar dust algorithm. Here, passive refers to the dust being advected with the gas, and scalar refers to the dust having a fixed grain size distribution.

We include dust production from Type II SNe and Asymptotic Giant Branch (AGB) stars, and further growth via accretion of metals, while destruction can occur from sputtering, consumption by star formation, or SN shocks. We explore the evolution of the galaxy dust mass function and the scaling relations of the DGR and DTM with metallicity over cosmic time. We then build on this, and investigate the physical drivers of the DGR and DTM using a machine-learning framework trained by our simulated dataset to understand the scatter in the DGR/DTM-metallicity relation. We use these tools to develop an algorithm (that we release publicly) for the dust mass from galaxies without the assumption of an overly simplistic dust-to-gas or dust-to-metal ratio. We additionally provide a simple scaling relation for the DGRs in galaxies.

This paper is organized as follows. In §2, we summarize the Simba simulation suite, with a particular focus on the model for dust formation and evolution. We present the dust mass functions and scaling relations between the DGR/DTM and gas phase metallicities in §3. In §3.4, we model the underlying physical drivers of the DGR and DTM, and establish a connection between the DGR/DTM and various physical properties of galaxies. We then discuss our results, compare them to other theoretical work, and discuss potential caveats in §4, and conclude in §5.

2 Methodology

2.1 Cosmological Simulations

This work utilizes the Simba cosmological hydrodynamic simulation . We refer the reader to Davé et al. (2019) for full details, and we summarize the salient points here.

The primary Simba simulation we use here has $1024^{3}$ dark matter particles and $1024^{3}$ gas elements in a cube of $100h^{-1}{\rm Mpc}$ side length, and is run from $z=249$ down to $z=0$. We assume a Planck16 (Planck Collaboration et al., 2016) concordant cosmology of $\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$, $\Omega_{b}=0.048$, $H_{0}=68\;{\rm km}\,{\rm s}^{-1}\;{\rm Mpc}^{-1}$, $\sigma_{8}=0.82$, and $n_{s}=0.97$. Our Simba run has a minimum gravitational softening length $\epsilon_{\rm min}=0.5h^{-1}{\rm kpc}$, mass resolution $9.6\times 10^{7}\;{\rm M}_{\odot}$ for dark matter particles and $1.82\times 10^{7}\;{\rm M}_{\odot}$ for gas elements. The system is evolved using a forked version of the Gizmo cosmological gravity plus hydrodynamics solver (Hopkins, 2015), in its Meshless Finite Mass (MFM) version. This code, modified from Gadget-3 (Springel, 2005), evolves dark matter and gas elements together including gravity and pressure forces, handling shocks via a Riemann solver with no artificial viscosity.

Radiative cooling and photoionisation heating are modeled using the Grackle-3.1 library (Smith et al., 2017), including metal cooling and non-equilibrium evolution of primordial elements. An H2-based star formation rate is used, where the H2 fraction is computed based on the sub-grid model of Krumholz et al. (2009) based on the metallicity and local column density, with minor modifications as described in Davé et al. (2016) to account for variations in numerical resolution. The star formation rate is given by the H2 density divided by the dynamical time: SFR$=\epsilon_{*}\rho_{\rm H2}/t_{\rm dyn}$, where we use $\epsilon_{*}=0.02$ (Kennicutt, 1998). These stars drive winds in the interstellar medium. This form of feedback is modeled as a two-phase decoupled wind, with 30% of wind particles ejected hot, i.e. with a temperature set by the supernova energy minus the wind kinetic energy. The modeled winds have an ejection probability that scales with the the galaxy circular velocity and stellar mass (calculated on the fly via fast friends-of-friends galaxy identification). The nature of these scaling relations follow the results from higher-resolution studies in the Feedback In Realistic Environments zoom simulation campaign (e.g. Muratov et al., 2015; Anglés-Alcázar et al., 2017b; Hopkins et al., 2014, 2018).

The chemical enrichment model tracks eleven elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe) during the simulation, with enrichment tracked from Type II supernovae (SNe), Type Ia SNe, and Asymptotic Giant Branch (AGB) stars. The yield tables employed are: Nomoto et al. (2006) for SNII yields, Iwamoto et al. (1999) for SNIa yields, and AGB star enrichment following Oppenheimer & Davé (2006). Type Ia SNe and AGB wind heating are also included, along with ISM pressurisation at a minimum level as required to resolve the Jeans mass in star-forming gas as described in Davé et al. (2016).

Simba incorporates black hole physics. Black holes are seeded and grown during the simulation via two-mode accretion. The first mode closely follows the torque-limited accretion model presented in Anglés-Alcázar et al. (2017a), and the second mode uses Bondi accretion, but solely from the hot gas component. The accretion energy is used to drive feedback that serves to quench galaxies, including a kinetic subgrid model for black hole feedback, along with X-ray energy feedback. Simba additionally includes a dust physics module to track the lifecycle of cosmic dust, which we describe in the following section.

2.2 Modeling the Dust Lifecycle

In our implementation, dust is fully coupled with gas flows. This treatment is essentially accurate, as the drift caused by the gas-dust drag force and the radiative pressure is under-resolved in our simulations. Additionally, dust grains are assumed to have the same physical properties with a constant radius $a\ =\ 0.1\ {\rm\mu m}$ and density $\sigma=2.4\ {\rm g\ cm^{-3}}$ Draine (2003). We ignore dust cooling channels which will be implemented in future work.

Dust is produced by condensation of a fraction of metals from SNe and AGB ejecta. We follow the prescription described by Equation (4) to (7) in Popping et al. (2017) which updates the work of Dwek (1998). In the following, $m_{i,d}^{j}$ refers to the dust mass of the $i$th element (C, O, Mg, Si, S, Ca, Fe) produced by the $j$th stellar process (SNII or AGB stars), whereas $m_{i,{\rm ej}}^{j}$ refers to the mass of ejecta from the $j$th process.

The mass of dust produced by AGB stars with a carbon-to-oxygen mass ratio C/O $>$ 1 is expressed as

[TABLE]

where $\delta_{i}^{\rm AGB}$ is the condensation efficiency of element $i$ for AGB stars. The mass of dust produced by AGB stars with C/O $<$ 1 is expressed as

[TABLE]

where $\mu_{i}$ is the mass of element $i$ in atomic mass units. The mass of dust produced by Type II SNe is described as

[TABLE]

where $\delta_{i}^{\rm SNII}$ is the condensation efficiency of element $i$ for SNII.

We choose a fixed dust condensation efficiency $\delta^{\rm AGB}_{i,\rm dust}=0.2$ based on the computation of Ferrarotti & Gail (2006), and $\delta^{\rm SNII}_{i,\rm dust}=0.15$ guided by Bianchi & Schneider (2007) to match the low-metallicity end of the observed $z=0$ relation between the DGR and gas-phase metallicities (Rémy-Ruyer et al., 2014)111Note - because our condensation efficiency for Type II SNe is tuned to match the low-$Z$ end of the local DGR-metallicity relation, the $z=0$ version of this relationship in the remainder of this paper should be treated as matching observations by construction, and not as a bona fide prediction.. We omit the condensation of Type Ia SNe ejecta, as recent work suggests that Type Ia SNe are not significant sources of dust production (see Nozawa et al., 2011; Dwek, 2016; Gioannini et al., 2017). This is different from McKinnon et al. (2016) and Popping et al. (2017) where Type Ia SNe are assumed to have the same condensation efficiency as Type II SNe.

Once dust grains are seeded, they grow by accreting gas-phase metals. Following Dwek (1998), the growth rate is expressed as

[TABLE]

where $M_{\rm metal}$ is the total mass of dust and local gas-phase metals. Following Hirashita (2000) and Asano et al. (2013), the accretion time scale $\tau_{\rm accr}$ is

[TABLE]

where $\rho_{g}$, $T_{g}$ and $Z_{g}$ are the local gas density, temperature and metallicity, respectively. $\rho_{\rm ref}$, $T_{\rm ref}$ and $Z_{\rm ref}$ are the reference values correspondingly. We take $\rho_{\rm ref}=100$ H atoms cm*-3*, $T_{\rm ref}=20$ K and $\tau_{\rm ref}=10$ Myr in this work.

Dust grains can be eroded by colliding with thermally excited gas especially in hot halos (e.g. Barlow 1978; Draine & Salpeter 1979b; Tielens et al. 1994). We adopt the approximation of the thermal sputtering rate of grain radii derived by Tsai & Mathews (1995), following McKinnon et al. 2017 and Popping et al. 2017. The sputtering time scale is expressed as

[TABLE]

where $\omega$ = $2.5$ controls the low-temperature scaling of the sputtering rate and $T_{0}\ =\ 2\times 10^{6}$ K is the temperature above which the sputtering rate flattens. The growth rate of dust mass due to thermal sputtering is then calculated by

[TABLE]

Because SN blast waves are not resolved in our simulations, we implement a subgrid model for dust destruction by SN shocks (Dwek & Scalo, 1980; Seab & Shull, 1983; McKee et al., 1987; McKee, 1989). The characteristic time scale $\tau_{\rm de}$ is

[TABLE]

where $M_{g}$ is the local gas mass, $\epsilon=0.3$ is the efficiency with which grains are destroyed in SNII shocks (McKee, 1989), $\gamma$ is the local SNII rate, and $M_{s}$ is the mass of local gas shocked to at least 100 km/s, calculated using the Sedov-Taylor solution to a homogeneous medium of $n_{\rm H}=0.13$ H atoms cm*-3* (the minimum SF threshold density of our simulations).

We additionally destroy dust completely in hot winds and during star formation and AGN X-ray heating (§2.1). The parameters adopted in this simulation is listed in Table 1.

Finally, we note that for the star formation and grain growth models, we need to provide a total metallicity in solar units. For this, we assume a solar abundance ($Z_{\odot}=0.0134$) taken from Asplund et al. (2009).

2.3 Data Analysis – A Machine-Learning Framework

We seek to accurately quantify how galaxy dust properties, particularly the DGR, trace other global galaxy properties. This represents a regression problem, where from a set of input variables, the prediction for the DGR is desired that most closely follows what is predicted directly by the simulation.

We employ machine learning for this regressor, as is now becoming common for a wide variety of astrophysical applications (e.g. Ball et al., 2007; Fiorentin et al., 2007; Gerdes et al., 2010; Carrasco Kind & Brunner, 2013; Ness et al., 2015; Kamdar et al., 2016; Agarwal et al., 2018; Rafieferantsoa et al., 2019). Taking advantage of the large training set offered by Simba simulation of tens of thousands of galaxies, we use machine learning to relate the galaxy DGR to a set of galaxy properties, i.e. an $N$-dimensional vector $X$, the components of which are the global galaxy properties as detailed in §3.4.

The primary algorithm used in this work is extremely randomized trees (ERTs; Geurts et al. 2006). ERTs build a large ensemble of regression trees, each of which splits the training set – here, an $(N+1)$-dimensional space comprising of data points ($X$,DGR) from 70% of the simulated galaxies – recursively among one randomly-selected subset of the galaxy properties. Each splitting divides the $(N+1)$-dimensional space into two $(N+1)$-dimensional subspaces, and it stops once the resulting subspace only contains one ($X$,DGR) point or the user-defined maximum tree depth is reached, in which case a relation between $X$ and DGR is established. The estimates produced by all the regression trees in the ERT ensemble are averaged to build a final map from $X$ to DGR. We refer readers interested in further details to Geurts et al. 2006 for the details of splitting and randomization in ERTs. For this work, we used the implementation of ERTs in the python package, scikit-learn (Pedregosa et al., 2012).

3 Dust Properties Over Cosmic Time

3.1 Dust Mass Functions

Figure 1 shows the redshift evolution of dust mass function (DMF), comparing against the observational result of Dunne et al. (2003) at $z=2$, Eales et al. (2009) at $z=1$, and Dunne et al. (2011); Clemens et al. (2013) and Beeston et al. (2018) at $z=0$. Unlike the comparison presented in Davé et al. (2019), here we standardized their results to our cosmological parameters (c.f. §2.1) and our assumed dust mass absorption coefficient $\kappa(850\ \mu{\rm m})=0.77\ {\rm cm^{2}\ g^{-1}}$.

At $z=0$, Simba agrees well with observed data. Our simulation underproduces the DMF at the low mass end, due to our mass resolution and the minimum mass of identified galaxies (24 baryonic particles $\approx 4.37\times 10^{8}\;{\rm M}_{\odot}$ baryonic mass). The $z=2$ model dust mass function under-predicts the observational one by a modest factor $\sim 3$. This is still much better than early attempts in this area, where galaxies with $M_{d}\gtrsim 10^{8}M_{\odot}$ are hardly produced (e.g. McKinnon et al., 2017). We note that the observational mass function by Dunne et al. (2003) and Dunne et al. (2011) are from surveys of sub-mm sources with large beam sizes, which could result in multiple objects being blended within one beam therefore overestimating their dust masses (e.g. Narayanan et al., 2010; Hayward et al., 2013; Narayanan et al., 2015). Beyond this, once we fold in the uncertainties in deriving dust masses from sub-millimetre photometry, it is probably premature to use the high-redshift DMF as a strong constraint on models. Overall, the dust mass function predicted by Simba broadly agrees with currently observed determinations, with very good agreement in the overall dust mass function shape. This indicates that Simba viably models dust evolution over cosmic time in galaxies, and sets the stage for examining more detailed dust-related properties.

Figure 2 shows the ratio of the comoving cosmic dust mass density $\rho_{\rm d}$ and the comoving cosmic stellar mass density $\rho_{*}$ as a function of redshift. The cosmic dust (or stellar) mass density is computed by summing the dust masses $M_{d}$ (or stellar masses $M_{*}$) of all gas cells and dividing by the total comoving volume. We get $\rho_{*}=310\rho_{\rm d}$ at $z=0$. For the convenience of comparison, $\rho_{*}$ is normalized such that $\rho_{*}=\rho_{\rm d}$ at $z=0$.

The cosmic dust mass density rapidly increases from $z=6$ by over 1.5 dex. At late times, the dust density flattens as the global star formation rate falls, the amount of metals available to be accreted drops, and a quasi-balance is reached among dust production, growth, destruction, and astration. The evolution of the dust mass density only slightly lags behind the stellar mass density at early epochs and catches up with rapid grain growth.

The comoving dust mass density excluding dust ejected out of galaxies via galactic winds is shown as the cyan line. It increases monotonically at high redshifts, following the total dust mass density. At $z\sim 2$, it starts declining as star formation rates decline on average with the onset of quenching massive galaxies, which slows down the metal enrichment and thus limits the grain growth. Meanwhile the destructive processes remain strong, and are even enhanced in massive galaxies that harbor little cold gas. Comparing this trend to the evolution of total dust mass density, we infer that the destruction of dust ejected to halos, dominated by thermal sputtering, is not strong enough to quickly eliminate dust grains owing to the low gas density. The wind model may need to be modified so that dust can be efficiently destroyed during the decoupled wind phase from galaxies into circum-galactic gas.

3.2 The Dust-to-Gas Ratio

The agreement of our predicted dust mass functions with observations suggests that Simba represents a plausible dust evolution model. We now turn to examining the main target relations of our paper, the dust-to-gas and dust-to-metal ratios in Simba galaxies, in comparison to data.

Figure 3 shows the hexbin plot of the DGR as a function of gas-phase metallicity ($Z_{\rm gas}$) between $z=0-6$, compared against the data points as observed by Rémy-Ruyer et al. (2014) and De Vis et al. (2019). Hexbins are color-coded with specific star formation rates $\text{sSFR}=\text{SFR}/M_{*}$.

At $z=0$, Simba shows a DGR that is in good agreement with observations. However, we emphasise that our dust model was tuned to do so via our choice of the dust condensation efficiencies. These quantities mostly change the $Z<0.2Z_{\odot}$ part of the DGR, without changing the slope much. Hence the slope of the DGR vs. $Z_{\rm gas}$ is a robust prediction, as is our predictions for the redshift evolution.

There are two main regimes in the DGR$-Z_{\rm gas}$ plane. The first regime corresponds to star-forming galaxies, where the DGR increases with $Z_{\rm gas}$. Our models predict a weak evolution of the DGR-$Z_{\rm gas}$ relation between $z=0$ to $z=6$, as is also predicted by Popping et al. (2017) using a semi-analytic model. The evolution of this relation in the DGR–$Z_{\rm gas}$ plane is mainly driven by the metal enrichment of galaxies as more galaxies just move along the sequence to slightly higher $Z_{\rm gas}$ at lower redshifts. The second regime corresponds primarily to quenched galaxies, and shows low DGR values and no correlation with $Z_{\rm gas}$. This is driven by AGN feedback that builds up the quenched galaxy population, in which dust production is stopped but dust destruction is enhanced by sputtering in surrounding hot gas. Even though there are no observed galaxies in this regime, this is likely due to an observational selection effects rather than an actual lack of galaxies (De Vis et al., 2019). Finally, we note that the predicted DGR ratios in low-$Z_{\rm gas}$ galaxies generally lie below the observations. These objects are typically very gas-rich dwarfs. It is unclear whether they are overly gas-rich in Simba, or else they have too little dust production (or both). There may also be observational selection effects that bias in favor of higher dust masses in such small, faint systems.

We determine a best-fit power law to the DGR-$Z_{\rm gas}$ relation in the star-forming regime. We separate the star-forming sequence from quenched galaxies by applying a density-based spatial clustering of applications with noise (DBSCAN) algorithm (Ester et al., 1996) to galaxies in {DGR, $Z_{\rm gas}$, sSFR} space. We then fit the star-forming sequence in the DGR–$Z_{\rm gas}$ plane using a power law:

[TABLE]

This is quite close to the best fit power law to the De Vis et al. (2019)’s data $\log\text{DGR}=(2.45\pm 0.12)\log\left(\frac{Z}{Z_{\odot}}\right)-(2.0\pm 1.4)$, quantitatively confirming the good agreement of the predicted and observed slopes. Still, the scatter of the best fit is large ($\sigma=0.31$ dex), even though the correlation is clear. Better estimates of the DGR might be obtained by incorporating secondary physical parameters in addition to $Z_{\rm gas}$; we explore this in §3.4 using machine learning.

3.3 The Dust-to-Metal Ratio

In Figure 4, we apply a similar analysis to the DTM, plotting the DTM ratio as a function of $Z_{\rm gas}$ between $z=0-6$. This relation is equivalent to the relation between the dust-to-gas ratio and $Z_{\rm gas}$. As discussed above, our simulation shows only a weak evolution of the DTM ratio from $z=6$ to $z=0$ and approximately reproduces the result of Rémy-Ruyer et al. (2014) and De Vis et al. (2019). The DTM ratio increases rapidly as $Z_{\rm gas}$ increases at the low metallicity regime until it is capped at a roughly constant value $\sim 0.8$ when $Z_{\rm gas}>0.5Z_{\odot}$. For the quenched galaxies, the DTM ratio drops off quickly for the same reasons as in the DGR case.

Overplotting the $z\sim 2-4$ data from high redshift observations of Damped Lyman Alpha (DLA) and Gamma Ray Burst (GRB) absorbers (De Cia et al., 2013, 2016; Wiseman et al., 2017) against the simulated data, however, we find a systematic discrepancy at $Z_{\rm gas}<0.5Z_{\odot}$, where the predicted DTM ratios show a much steeper dependence with metallicity than the observations. The source of this discrepancy is unclear. We note that the nature of these absorbers can be significantly different from the physical conditions in galaxy disks. DLAs are thought to arise from the outskirts of gas disks in galaxies and perhaps even from metal-poor gas in the circum-galactic medium (Berry et al., 2014). Moreover, these studies measure metallicities and DTM via abundances acquired from optical/UV absorption-line spectroscopy, which is different from methods typically used for galaxies observed in the local universe (i.e. strong line calibrations for metallicities and infrared emission for the dust mass). This discrepancy thus may reflect a difference in the dust production versus destruction in different environments or scales, combined with the selection effect and methodologies of observations. We defer a more careful comparison of the DGR in these particular types of objects to future work, but note that there is potentially a discrepancy at low-$Z$. This also may be responsible for the too-low DGR at low-$Z$. We note, however, that significantly larger amounts of dust in low-mass galaxies would steepen the dust mass function, which may put our currently viable predictions into conflict with observations.

3.4 Better DGR and DTM prediction via machine learning

We now investigate the physical drivers of the DGR and DTM ratios in galaxies, as well as their scatter. While equation 9 provides a rough fit, it is clear that there is correlated scatter, and hence fitting with more variables should give a tighter relation. To approach this agnostically in terms of the form and input variables, we employ a machine learning approach using Extremely Randomized Trees (ERTs).

We seek to establish a map between a range of physical properties – namely gas-phase metallicity ($Z_{\rm gas}$), gas depletion time scale ($\tau_{\rm depletion}\equiv M_{g}/{\rm SFR}$), stellar mass ($M_{*}$), half-baryonic mass radius ($R_{b,50}$), gas mass fraction ($f_{g}\equiv M_{g}/(M_{g}+M_{*})$), and gas surface density ($\Sigma_{g}$) – onto the galaxy DGR and DTM ratios. We limit our analyses to $z=0$ and concentrate our efforts on fitting the star-forming sequence of galaxies, since quenched galaxies show little dust and no obvious trend with any physical property and our relations show little evolution with redshift. It would be straightforward to apply this methodology to other redshift outputs.

The ERT is set up using $70\%$ of the selected Simba galaxy sample as a training set, with a maximum depth of 20 levels. We then use the remaining 30% for validation. The algorithm then generates a mapping between the inputs and the two desired outputs (DGR and DTM). By using ERT, we also have access to importance levels, which are determined as the relative depth of a given input parameter was used for branching the tree.

Figure 5 shows the fitted relation between the DGR and various galaxy physical properties. Cyan points show the training set, and orange points denote the validation set; taken together, they represent all simulation star-forming galaxies. Observations are shown from De Vis et al. (2019) as black points.

Comparing to observations, Simba reproduces the observed $z=0$ DGR values as a function of various galaxy physical quantities reasonably well. The DGR increases with metallicity and stellar mass, though less tightly so with the latter. The DGR also drops with the gas fraction and depletion time, probably reflecting underlying trends from the stellar mass dependence. The DGR shows no clear dependence on $R_{b,50}$ or $\Sigma_{g}$. For $Z_{\rm gas}$, $\tau_{\rm depletion}$, $M_{*}$ and $f_{g}$, the Spearman’s rank correlation coefficients are $r=0.87,-0.63,0.64,0.81$, respectively, compared to 0.33 for $R_{b,50}$ and 0.14 for $\Sigma_{g}$.

Figure 6 shows the analogous plot for the DTM ratio. The trends are broadly similar, with DTM increasing with $Z_{\rm gas}$ and $M_{*}$, and decreasing with $f_{g}$ and $\tau_{\rm depletion}$. This suggests that $M_{*}$, $f_{g}$, and $\tau_{\rm depletion}$ may provide additional information that will enable tighter prediction of both the DGR and the DTM. Rank correlation coefficients are similar to the DGR case.

Since we know the true values for the 30% validation set, we can examine how well the ERT is able to reproduce these true values. The quality of fitting is shown in the left panels of Figure 7, for the DGR (top panels) and DTM (bottom). We find a very tight relation with a $\sim 0.15$ dex scatter estimated by the mean squared error (MSE) between the predicted DGR and true (simulated) DGR of the cross-validation set. This scatter is significantly reduced from $\sim 0.28$ dex when only $Z_{\rm gas}$ is used, showing that the machine learning is effective at generating better predictions for the DGR. Similarly, for the DTM, the scatter is reduced from 0.27 dex when only $Z_{\rm gas}$ is used, to $0.14$ using the full ERT mapping.

To examine the sensitivity to the ERT tree depth, we show in the right panel of Figure 7 the MSE as a function of tree depth. We see that increasing the tree depth initially greatly improves predictions, but beyond a depth of $\ga 9$ levels, there is essentially no improvement. This is true for both the DGR and the DTM ratio. At this point, given the sample size and the number of parameters, there is no more information contained in additional tree levels. Hence we find an optimal maximum ERT tree depth of 9 levels for this sample.

Finally, we examine importance levels of the input physical parameters, as returned by the ERT algorithm, shown in Figure 8. At the optimal depth of 9 levels, the left panel shows that the DGR appears to be most directly correlated with the metallicity, followed by the gas fraction. The depletion time and stellar mass also assist with the fitting at a lower level. As expected, the half-mass radius and the gas surface density do not contribute significant information.

The trends are broadly similar for the DTM, as shown in the right panel of Figure 8. However, it now appears that the metallicity, depletion time, and gas fractions all show similar levels of importance. The stellar mass still shows lower importance, and $R_{b,50}$ and $\Sigma_{g}$ are again irrelevant.

In summary, our ERT-based machine learning framework is able to significantly improve the predictive power for the DGR and DTM relations. Using only $Z_{\rm gas}$ results in scatters of $\sim 0.3$ dex, while using the ERT-generated mapping reduces the scatter to $\sim 0.15$ dex. The key quantities driving this are the metallicity, gas fraction, depletion time, and (to a lesser extent) stellar mass, while in Simba the dependencies on the baryonic half-mass radius and gas surface density are negligible. The map determined via ERT using Simba can be applied by modelers who usually have no information about dust or do not track dust evolution in a self-consistent way, which will provide a more accurate estimation of dust mass than that based on a simple assumption of DGR (or DTM). Alternatively, observers with information about these global galaxy quantities can utilise these algorithms to estimate the dust-to-gas or dust-to-metals ratios in their galaxies.

4 Discussion

4.1 Physical Underpinnings

We now take a deeper dive in to the details of the trends of the DGR with various physical properties, by examining the correlations of the DGR versus various galaxy physical quantities as shown in Figure 5. We focus our discussion on the DGR, as the trends and interpretations for the DTM ratio are analogous.

In panel (a) of Figure 5 we show the DGR-$Z_{\rm gas}$ relation. At $Z<0.2Z_{\odot}$, the relation is roughly linear, corresponding to an approximately constant but low DTM (c.f Figure 6a). While there is a steep increase from $Z\sim 0.2Z_{\odot}$ to $Z\sim 0.5Z_{\odot}$, the relation is again roughly linear at $Z>0.5Z_{\odot}$, corresponding to an approximately constant DTM $\sim 0.8$ (see Figure 6a). These three different trends shows three different regimes of dust enrichment. At lower $Z_{\rm gas}$, the galaxies are under-evolved and the dust enrichment is dominated by dust production via condensation of ejecta from late-stage stars. In the intermediate $Z_{\rm gas}$ regime, dust growth via accreting gas-phase metals gradually take over the enrichment process. At higher $Z_{\rm gas}$, the growth is extremely strong and the dust mass is mainly determined by the gas-phase metals available for accretion.

In panel (c) of Figure 5, we show the DGR-$M_{*}$ relation. The relation flattens after $M_{*}\ga 10^{9.5}M_{\odot}$. This transition follows the flattening of the metal-metallicity relations above a comparable $M_{*}$ (Tremonti et al. 2004; see Figure 9 in Davé et al. 2019) plus the fact that the DGR primarily traces the galaxy metallicity as Figure 5(a) shows. At higher masses, most galaxies are quenched, so the DGR actually drops and has a large scatter owing to dust destruction.

In Figure 5(b) and Figure 5(e) we show the DGR-$\tau_{\rm depletion}$ and DGR-$f_{g}$ relation respectively. Galaxies with lower $f_{g}$ and $\tau_{\rm depletion}$, which implies that they are quiescent and highly evolved, tend to have a higher DGR. The relations are flatter at low $f_{g}$ and $\tau_{\rm depletion}$, because the rapid grain growth due to abundant metals is countered by enhanced destructive processes, i.e. shock waves from supernovae, thermal sputtering and astration. These trends shows the correlation between DGR and galaxy evolutionary stages which proceeds as star formation deplete gas and build up metals. We note, that there are high-$M_{*}$ (usually high-$Z_{\rm gas}$) galaxies, whose $f_{g}$ is still relatively high and whose star formation rates are not highly suppressed, having a relatively low DGR. A similar situation applies to some low-$M_{*}$ galaxies. This contributes to the scatter in the DGR-$Z_{\rm gas}$ plane.

We conclude that the DGR (and DTM) can be determined by $Z_{\rm gas}$ along with $M*$ which reflects the chemical enrichment history and $\tau_{\rm SF}$ along with $f_{g}$ which indicates the evolutionary stage of galaxies. Physically, this suggests that metal enrichment history reflected by $Z_{\rm gas}$ and $M_{*}$, and evolutionary stages quantified by $\tau_{\rm depletion}$ (see also Asano et al. 2013; Zhukovska 2014; Feldmann 2015) and $f_{g}$ (see also De Vis et al. 2019) are the main drivers of the scatter in the DGR$-Z_{\rm gas}$ plane. Meanwhile, galaxy compactness as quantified by $R_{b,50}$ and $\Sigma_{g}$ does not seem to impact the DGR or DTM, showing that at least in Simba dust content is insensitive to galactic structure – with the caveat that given Simba’s $\sim 1$ kpc resolution, galactic structure may not be faithfully modeled in detail.

4.2 Comparison with other models

The DMF and the DGR-$Z_{\rm gas}$ relation have been studied by cosmological hydrodynamic simulations (McKinnon et al., 2017; Hou et al., 2019) and semi-analytic models (Popping et al., 2017). Like McKinnon et al. (2017) and Popping et al. (2017), our work predicts that the dust mass function increases monotonically from z = 2 to 0 at the high-mass end ($M_{d}\gtrsim 10^{8}M_{\odot}$), and is unable to simultaneously match the $z=0$ to $z=2$ DMF. Nevertheless, our result appears to have the closest match to observations to date, underpredicting the $z=2$ DMF by a factor $\sim 3$, which is significantly better than McKinnon et al. (2017) where galaxies with $M_{d}\gtrsim 10^{8}M_{\odot}$ are hardly produced. On the other hand, Hou et al. (2019) tracks two types of grains (small/large) and grain-grain shattering and coagulation. They also implement a subgrid-model to boost the density of unresolved dense gas, which is not adopted by our model (therefore we use a short $\tau_{\rm ref}$) and a simple AGN feedback model to suppress star forming activity of massive galaxies. They are able to reproduce the non-monotonic trend of the high-mass-end DMF evolution from $z=2$ to $z=0$, which was observed by Dunne et al. (2011). However, they fail to match the observed DMF at $z=0$ (where they overproduce high $M_{d}$ galaxies) to $z=2$ (where they underproduce high $M_{d}$ galaxies).

What steps forward are necessary for simulations to match the observed $z=2$ DMF? Some possible solutions include: (1) implementing dust yields in stellar ejecta as a function of a star’s mass and metallicity (Ferrarotti & Gail, 2006; Bianchi & Schneider, 2007; Zhukovska et al., 2007; Nanni et al., 2013; Schneider et al., 2014), or the local ISM density or temperature, as pointed out by McKinnon et al. (2017); (2) tracking the evolution of grains with different sizes (see e.g. Hou et al. 2019 for two-size grains and McKinnon et al. 2018; Aoyama et al. 2019 for a continuous distribution of grain sizes). Solution #2 will be implemented in an upcoming paper. We expect that a non-monotonic evolution of high-mass-end DMF would result from the intensified grain-grain collisions at lower redshift which will generally lead to an increasing abundance of smaller grains that experience faster destruction.

Turning to the DGR-$Z_{\rm gas}$ relation: The simulations by McKinnon et al. (2017) obtain a rather flat relation between the DGR and gas-phase metallicity from low $Z_{\rm gas}$ to high $Z_{\rm gas}$, mainly because their accretion timescale for dust growth (see Equation (5) of McKinnon et al. 2016) does not vary with the local ISM $Z_{\rm gas}$. We find that the dependency of the accretion timescale on metallicity is essential to reproduce the observed DGR-$Z_{\rm gas}$ relation as is shown in Popping et al. (2017), Hou et al. (2019) and our work, in our case by the results of two test runs in Figure 9. If we assume a metallicity-independent accretion timescale, then we get a flatter and higher relation (blue points) compared to the Rémy-Ruyer et al. (2014) observations. The actual observations are slightly shallower than the predictions from the metallicity-dependent model, suggesting that perhaps our metallicity dependence should be softened somewhat.

Similar to Popping et al. (2017), our results show a weak evolution of the DGR-$Z_{\rm gas}$ relation from $z=0$ to $z=6$, especially for high metallicity galaxies. This is encouraging given that the underlying galaxy formation models (SAM vs hydrodynamic) are rather different from one another. Beyond this, the treatment of dust growth in both models is different: grain growth by accretion in our work does not use the information of any detailed subgrid ISM model, whereas in the Popping et al. (2017) SAM, the accretion time-scale is calculated by inferring the gas density in molecular clouds from SFR laws. The fact that these two very different models arrive at the same conclusion for a relatively modest evolution in the DGR-$Z$ relation likely underscores its robustness. On the other hand, Hou et al. (2019) shows that the DGR at a fixed $Z_{\rm gas}$ builds up significantly a factor of $\sim 5$ from $z=5$ to $z=0$. This discrepancy mainly comes from the different treatment of grain growth and feedback. Grain growth in our simulation is overall stronger, therefore our simulated galaxies are able to reach the quasi-static states within a shorter period of time. Moreover, Simba has a more sophisticated feedback mechanisms that suppress further metal enrichment at lower redshifts, particularly in massive galaxies.

Hou et al. (2019) explores the scaling relations with a variety of galaxy properties other than $Z_{\rm gas}$. Comparing Figure 5 (c) and (e) against their Figure 4 (b) and (d), we find that the trend predicted by both simulations are similar but there are big discrepancies at the high $M_{*}$ ($M_{*}>10^{9.5}M_{\odot}$, the turning point of mass-metallicity relation) and low $f_{g}$ end, where Hou et al. (2019)’s model overproduces the DGR compared against observational constraints. In this regime, our DGR-$M_{*}$ and DGR-$f_{g}$ relations flatten. The sophisticated black hole feedback model in Simba  compared to the simple phenomenological AGN feedback model used in their work which may underestimate the suppressing power, can explain the difference. A similar analysis applies to the low $f_{g}$ regime where most of the galaxies are highly evolved and massive, and black hole feedback mechanisms are influential.

4.3 Caveats

Here we point out some caveats of our simulation. First, we note that the choice of free parameters for dust production, grain growth and grain destruction via SNe shocks are not well constrained (c.f. Table 1) and most likely degenerate. Though grain growth dominates the evolution of dust content, a combination of free parameters different from what is chosen in this work could potentially lead to an equally good (or even slightly better) match to observations, e.g. one with stronger production, weaker grain growth and stronger destruction, or one with stronger grain growth and much stronger destruction. On the other hand, the observed metallicities depend strongly on which strong-line calibration is used. The match to De Vis et al. (2019) would be worse if e.g. calibrations of Pettini & Pagel (2004) are used instead of Pilyugin & Grebel (2016).

Besides, dust plays an important role in cooling and shielding gas, catalyzing important chemical reactions (e.g. formation of molecular hydrogen) and recombination processes, and alternating the interstellar radiation fields. However, dust in our simulation only affects the ISM by depleting gas-phase metals thus reducing the efficiency of gas cooling channels. Work to implement the dust physics in a fully self-consistent manner is under way.

Beyond this, in our current model implementation, dust grains are fully coupled to gas flows. In reality, dust grains can be decoupled from gas due to radiative forces and lack of pressure and experience/apply drag forces from/to gas (Squire & Hopkins, 2018; Hopkins & Squire, 2018a, b). Dust grains in this work are assumed to have the same grain size to capture the major processes that evolve the dust mass. In the future we will implement “active” dust particles that are not strictly coupled to gas to track the evolution of grain size distributions and take grain-grain collisions into account. As mentioned in §4.2, we expect the evolution of the grain size distribution would alter the evolution of the cosmic dust content. The use of dust super-particles to sample the spatial distribution of dust would also save the memory usage, making it computationally feasible to track multiple-size grains in a large volume cosmological simulation, even though we expect the effect on dynamics would be negligible because of the lack of spatial resolution of our large-volume cosmological simulations and strong radiative fields.

Finally, the 100$h^{-1}{\rm Mpc}$ Simba volume lacks the resolution to resolve a multi-phase ISM, which is a common issue for large-volume cosmological simulations. As a result, parameters such as the reference accretion time scale $\tau_{\rm ref}$ have to be tuned such that the effective gas density is boosted. We also assume fixed dust destruction and condensation efficiencies, which may actually be functions of local ISM properties (Seab & Shull, 1983; McKee et al., 1987; Ferrarotti & Gail, 2006; Bianchi & Schneider, 2007; Zhukovska et al., 2007; Yamasawa et al., 2011; Nanni et al., 2013; Schneider et al., 2014; Temim et al., 2015). Calibrating parameters we have used against numerical simulations of resolved ISM is one potential approach to make improvement.

5 Summary

We have developed a self-consistent model for the formation, growth and destruction of dust in the Simba cosmological galaxy formation simulation, and used these to compare to predictions to observed predictions, as well as study the physical drivers of the dust-to-gas and dust-to-metals ratios in galaxies. We also develop a machine learning framework to relate the dust-to-gas ratio (DGR) and dust-to-metal ratio (DTM) to various input global galaxy properties. Our main results are as follows:

Our main results follow:

• •

Simba broadly reproduces the observed dust mass functions across the measured redshifts of $z=0-2$, albeit with modest under-prediction at high-$z$. The low-mass end steepens at high redshift.

• •

We find a relationship between the dust-to-gas ratio of star-forming galaxies and their gas phase metallicity such that lower metallicity galaxies have lower dust-to-gas ratios. This is broadly in accord with observations. There is little evolution with redshift in this relationship. Meanwhile, quenched galaxies show lower DGR and essentially no relationship in this space (Figure 3). This non-constant dust-to-metals ratio with metallicity (eq. 9) has implications for galaxy formation models that historically have assumed a constant dust-to-metals ratio (see e.g. Silva et al., 1998; Granato et al., 2000; Baugh et al., 2005; Lacey et al., 2010; Narayanan et al., 2010, 2015; Narayanan et al., 2018a, b; Fontanot & Somerville, 2011; Niemi et al., 2012; Somerville et al., 2012; Hayward et al., 2013; Cowley et al., 2017; Katz et al., 2019; Ma et al., 2019).

• •

The DTM ratio vs. metallicity relation drops at low metallicity, akin to the DGR relationship (Figure 4). This is consistent with low redshift result, yet may be in tension with observational constraints at low metallicities at high redshifts from GRBs and DLAs. Other trends are qualitatively similar to those of the DGR.

• •

In order to help both modelers and observers estimate more accurate DGR and DTM ratios, we have developed a publicly-available machine learning framework that generates a mapping between the DGR and DTM ratio and a set of galaxy physical properties, showing that it depends significantly on various galaxy properties. The machine learning framework reduces the scatter in the prediction from $\sim 0.3$ dex in the DTM and DGR using a simple fit to $Z_{\rm gas}$ alone, down to $\sim 0.15$ dex using the machine learning framework. This code is available at https://bitbucket.org/lq3552/dust_galaxy_analyzer.

• •

While the DGR and DTM ratios depend most sensitively on the gas phase metallicity in galaxies, we demonstrate that there are important secondary relationships between these ratios and the depletion time scale, stellar mass and gas fraction of galaxies (Figure 5). The DGR and DTM ratio both drop to lower $M_{*}$, and rise to lower $f_{\rm gas}$ and gas depletion times. There is no dependence on gas surface density and baryonic half-mass radius. Hence dust content is governed by both long-term evolutionary processes such as metal content and stellar mass, as well as short-term variations such as varying gas content and (commensurately) depletion times.

Overall, the Simba dust model is at least as successful compared with other current dust models implemented in cosmological simulations. However, there remain various caveats and potential directions for improvement. These include having active dust (not tied to the gas), multiple dust grain sizes, and implementing more sophisticated dust cooling. Furthermore, there are various free parameters that are constrained indirectly by observations, which might be better constrained using high-resolution ISM simulation. By using such a multi-scale approach to combine high resolution simulations and observational constraints into a cosmological galaxy formation model, we are moving towards more comprehensively studying the evolution of galaxy dust on cosmological scales.

Acknowledgements

The authours acknowledge helpful discussions with Ryan McKinnon, Gergo Popping, and Paul Torrey, and thank the anonymous referee for constructive comments. QL acknowledges support from the University of Florida Informatics Institute. QL and DN acknowledge support from the NSF via grant AST-1715206, and the Space Telescope Science Institute via HST AR-15043.0001. RD acknowledges support from the Wolfson Research Merit Award program of the U.K. Royal Society. This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility. The equipment was funded by BEIS capital funding via STFC capital grants ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure. Some of the simulations presented here were performed on the University of Florida HiPerGator2.0 supercomputing facility, and the authours are grateful to the staff that supports the software and hardware on HiPerGator.