1. Introduction
In modern galaxy evolution studies, dust is crucial in how galaxies evolve and change across cosmic time. The abundance of dust in galaxies is directly connected with galaxy evolution (Santini et al. Reference Santini2014), as stars would not form effectively without it, as dust catalyses the formation of molecules (Wakelam et al. Reference Wakelam2017; Chen et al. Reference Chen2018) and enables the fragmentation of gas clouds (Omukai et al. Reference Omukai2005; Schneider et al. Reference Schneider2006). However, the exact abundance of this dust in galaxies is unclear, especially at earlier epochs when the age of the Universe was comparable to the typical dust formation time scales (Todini & Ferrara Reference Todini and Ferrara2001; Bianchi & Schneider Reference Bianchi and Schneider2007; Lésniewska, Aleksandra & Michałowski, Michał Jerzy Reference Lésniewska2019).
Modern galaxy evolution models have attempted to match the different dust abundances across galaxy populations at different epochs (McKinnon et al. Reference McKinnon2016; McKinnon et al. Reference McKinnon2017; Popping et al. Reference Popping2017; Hou et al. Reference Hou2019; Triani et al. Reference Triani2020; Vijayan et al. Reference Vijayan2019; Parente et al. Reference Parente2022; Parente et al. Reference Parente2023). Although significant advancements have been made in modelling the life cycle of dust, the mechanisms behind the creation, growth, and destruction of dust, particularly concerning dust content, remain a topic of debate in the modern era (Hensley & Draine Reference Hensley and Draine2023; Ragone-Figueroa et al. Reference Ragone-Figueroa2024; Yates et al. Reference Yates2024). Enhancing our models through direct comparisons of simulated and observed data is crucial for deepening our understanding of dust’s life cycle and its contents in galaxies and, consequently, galaxy evolution. Recent results suggest a substantial evolution of the dust content of galaxies since
$z\sim1$
(Parente et al. Reference Parente2023; Eales & Ward Reference Eales and Ward2024).
Despite these current efforts, most comparisons currently lack a sufficiently large sample size to reveal significant discrepancies, particularly in earlier epochs where surveys are lacking. Most modern comparisons with models and observations (Li et al. Reference Li2019; McKinnon et al. Reference McKinnon2017; Popping et al. Reference Popping2017, e.g.) rely on low sample size observational datasets nearly a decade old (Eales et al. Reference Eales2009; Dunne et al. Reference Dunne2011; Clemens et al. Reference Clemens2013) that focus on later epochs (
$0 \lt z \lt 1$
) containing 82, 1867, and 234 sources, respectively. The only outlier most studies utilise is Beeston et al. (Reference Beeston2018), which includes 15 750 sources at a redshift
$\lt 0.1$
due to its utilisation of more modern surveys. However, this is only at a small redshift range at later epochs. To properly understand the properties that govern the abundance of dust in galaxies, more complete samples at earlier epochs are necessary for enhancing our comparisons.
The current aim of this study is to compare the simulated and observed datasets with a larger sample size and later redshifts than previously accomplished to help remedy this issue. We strive to identify the discrepancies between the dust mass generated by modern galaxy evolution models and measured observational data. Using an extensive data set to identify these discrepancies, this paper serves as a first-effort basis for improving galaxy evolution models across epochs up to
$z=1.5$
.
The complete datasets of GAMA, G10-COSMOS, 3D-HST with stellar mass (left) and dust mass (right) up till the intermediate Universe (
$z\approx 1.5$
). The more yellow the data displayed on the figures, the greater the data contained in that area, and the bluer the less. The clumps observed correspond to different surveys. GAMA covers up to the near Universe (
$z \lt 0.5$
) and corresponds to the high-density clump at the top left of the plots. G10-COSMOS and 3D-HST cover up to the intermediate Universe (
$z \lt 1.75$
) and are harder to distinguish as there is little separation between the data.
We, therefore, utilise the dust mass data produced by the hydrodynamical model Simba and compare it using various methods with a large, homogeneous, and detailed dust mass dataset, as described in Driver et al. (Reference Driver2018). The observed dataset includes data from Galaxy and Mass Assembly (GAMA), an offshoot of GAMA, the Cosmic Evolution Survey (G10-COSMOS), and the Hubble Space Telescope 3D project (3D-HST). We note that the data received from G10-COSMOS and Simba were obtained via private communication and are not publicly available.
The paper is organised as follows: in Section 2, we summarise our observational datasets: GAMA, G10-COSMOS, and 3D-HST. In Section 3, we briefly describe the dust evolution model present in Simba and the simulations used in our paper. In Section 4, we explain our galaxy selection progress for both the observational and simulated datasets. In Section 5, we explore the results of our comparisons between the two datasets, and in Section 6, we discuss further the implications of these results and put them into context with other works. Finally, we finish with our conclusions in Section 7.
2. Data
We combine the three datasets outlined in Driver et al. (Reference Driver2018): GAMA (Driver et al. Reference Driver2011; Liske et al. Reference Liske2015), G10-COSMOS (Davies et al. Reference Davies2015; Andrews et al. Reference Andrews2017), and 3D-HST (Momcheva et al. Reference Momcheva2016). These studies contain data from the ultraviolet (UV), mid-infrared (MIR), and far-IR (FIR) wavelengths. The GAMA and G10-COSMOS studies also contain constraints from Herschel Space Observatory’s SPIRE (Griffin et al. Reference Griffin2010; Poglitsch et al. Reference Poglitsch2010) and PACS instruments (Eales et al. Reference Eales2010; Oliver et al. Reference Oliver2012), which allow for robust dust mass measurements. The three studies GAMA, G10-COSMOS, and 3D-HST, extend to the nearby (
$z \leq 0.5$
), intermediate (
$z \lt 1.75$
), and high-z Universe (
$z \lt 5.0$
), respectively. Each dataset comprises approximately 200 000 galaxies that sample a broad range in stellar mass, dust mass, and look-back time. Each dataset was subject to data cuts based on flux limitations set in Driver et al. (Reference Driver2018) and active galactic nuclei (AGN) contamination removal. Each sample was processed using magphys, a spectral energy distribution fitting (SED) code, to provide estimates of dust mass, stellar mass, and star-formation rates (SFR) based on the flux limitations (see Sections 2.1, 2.2, and 2.3). Refer to Figure 1 for an overview of the stellar and dust mass distributions of these combined datasets up till the intermediate Universe (
$z\approx 1.5$
).
2.1. GAMA
The GAMA survey (Driver et al. Reference Driver2009; Driver et al. Reference Driver2011; Liske et al. Reference Liske2015) is a complete (98 % to
$r \lt 19.8$
mag) spectroscopic survey up to the near universe (
$z \leq 0.5$
). The survey consists of five regions G02 (
$\sim 56$
deg
$^2$
), G09 (
$\sim 60$
deg
$^2$
), G12 (
$\sim 60$
deg
$^2$
), G15 (
$\sim 60$
deg
$^2$
), and G23 (
$\sim 51$
deg
$^2$
) for a total area of roughly 180 deg
$^2$
with about 70 000 galaxies per region. Driver et al. Reference Driver2018 uses LAMBDARCatv01, a catalogue for GAMA made with in-house software (lambdar; Wright et al. Reference Wright2016) which provides flux limits with their errors, upper-limits, and flags. These objects were filtered from an initial value of 200 246 to 128 568 based on specific quality criteria, with values subsequently adjusted to be compatible with magphys. They also implemented an additional cut step to eliminate high stellar mass outliers. We use the magphys SED fits data products with the high stellar mass cut (Driver et al. Reference Driver2018) for the described 128 568 objects in this study.
2.2. G10-COSMOS
The G10-COSMOS survey (Davies et al. Reference Davies2015; Andrews et al. Reference Andrews2017) is a subset of HST COSMOS survey that covers a 1 deg
$^2$
region of the survey (Scoville et al. Reference Scoville2007b, a) and is 100 percent redshift complete up to a specified flux limit of
$i \lt 25$
mag. The survey includes wavelengths from UV to FIR. The survey uses the data collected by the Galaxy Evolution Explorer (GALEX; Martin et al. Reference Martin2005; UV), Canada–France–Hawaii Telescope (CFHT; Cuillandre et al. Reference Cuillandre, Peck, Seaman and Comeron2012; optical), Subaru (Taniguchi et al. Reference Taniguchi2007; optical), HST (Scoville et al. Reference Scoville2007b; optical), Visible and Infrared Survey Telescope (VISTA; Emerson & Sutherland Reference Emerson, Sutherland, Tyson and Wolff2002; Jarvis et al. Reference Jarvis2013; NIR), Spitzer (Werner et al. Reference Werner2004; Sanders et al. Reference Sanders2007; MIR), and some Herschel (Pilbratt Reference Pilbratt and Mather2003; Shirley et al. Reference Shirley2021; FIR).
Similar to the GAMA dataset, G10-COSMOS utilises the lambdar G10CosmosLAMBDARCatv06 catalog (Andrews et al. Reference Andrews2017a; Wright et al. Reference Wright2016). This catalogue also similarly underwent AGN cuts based on different criteria outlined (Donely et al. Reference Donely2012; Seymour et al. Reference Seymour2008; Laigle et al. Reference Laigle2016), as well as the removal of high stellar mass outliers. It includes a total of 142 260 objects found in the intermediate Universe (
$z \lt 1.75$
). We use the magphys SED fits data products with the AGN cut (Driver et al. Reference Driver2018) for the 142 260 objects described in this study.
2.3. 3D-HST
The 3D-HST survey (Brammer et al. Reference Brammer2012; Momcheva et al. Reference Momcheva2016) is an almost complete (up to 85 percent at
$F160W = 26.0$
mag) photometric, grating-prism (GPRISM), and spectroscopic survey up to the high-z universe (
$z \lt 5$
) that covers a
$\sim0.2$
deg
$^2$
area of the sky. The fields 3D-HST covers are sub-regions of the All-Wavelength Extended Groth Strip International Survey (AEGIS) (Davis et al. Reference Davis2007), Cosmic Evolution Survey (COSMOS) (Scoville et al. Reference Scoville2007b), Great Observatories Origins Deep Survey (GOODS-South) (Nonino et al. Reference Nonino2009), GOODS-North (Barger et al. Reference Barger2008), UKIRT Infrared Deep Sky Survey (UKIDSSUDS) Almaini et al. (Reference Almaini, Metcalfe and Shanks2007), and Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS) (Grogin et al. Reference Grogin2011; Koekemoer et al. Reference Koekemoer2011). Unlike GAMA and G10-COSMOS, 3D-HST does not include FIR coverage from Herschel, which is currently under development (Hurley et al. Reference Hurley2017). The 3D-HST catalog consists of 204 294 galaxies and AGN, with a redshift estimate of 3 839 from spectroscopic, 15 518 from GRISM, and 185 843 from photometric data. We use the magphys SED fits data products from the fair AGN cut, which follows the criteria set by Donely et al. (Reference Donely2012), as well as the exclusion of high stellar mass outliers (Driver et al. Reference Driver2018), for the 188 235 objects discussed in this study.
2.4. magphys
To calculate the physical properties for dust mass, stellar mass, and SFR used in this study, the three datasets containing 128 568, 142 260, and 188 235 objects, respectively, were parsed through the SED-fitting code magphys. For further details on the data preparation process and magphys, we point the reader to Driver et al. Reference Driver2018 and da Cunha et al. Reference da Cunha2008. This parsing through magphys standardises the assumptions and systematics of dust mass, stellar mass, and SFR estimates between the datasets, which minimises the errors between the surveys and allows us to combine the datasets and make proper comparisons with the simulated galaxies in Simba.
3. The Simba simulation
In this study, we utilise the cosmological simulation Simba (Davé et al. Reference Davé2019), which is the successor to mufasa and is built upon the gizmo hydrodynamic code (Hopkins & Lee Reference Hopkins and Lee2015), for our comparisons. Here, we will present a summary of the relevant physics in Simba to help explain our findings in Section 5.
Simba presents a self-consistent model for dust production, growth, and destruction (Li et al. Reference Li2019; Davé et al. Reference Davé2019). This model tracks eleven elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe) across cosmic time, with enrichment sourced from Type II SNe, Type Ia SNe, and AGB stars. In Simba, dust is fully coupled with gas flows, a treatment considered accurate due to the simulation’s under-resolved drag force and radiative pressure. Dust grains are also modeled with a consistent size of
$0.1 \mu \text{m}$
and a density of
$2.4 \, \text{g cm}^{-3}$
(Draine Reference Draine2003).
Dust production in Simba is determined by taking fixed fractions of metals from Type II SNe and AGB star condensations, based on the work of Popping et al. (Reference Popping2017), which updates earlier findings by Dwek (Reference Dwek1998). Below are the equations (outlined in Li et al. Reference Li2019; Davé et al. Reference Davé2019) that describe how dust mass is determined using AGB stars and Type II SNe where
$\text{M}_{i,d}^{\text{j}}$
refers to the ith element (C, O, Mg, Si, S, Ca, and Fe) produced by the jth stellar process (AGB stars or SNeII) and where
$\text{ M}_{i,ej}^{\text{j}}$
refers to the mass ejected from the jth process.
The mass of dust produced by AGB stars with a carbon-to-oxygen ratio of greater than one (
$C/O \gt 1$
) is expressed as
\begin{equation} \text{M}_{i,d}^{\text{AGB}} = \begin{cases} \delta_C^{\text{AGB}} \left( \text{M}_{C,ej}^{\text{AGB}} - 0.75 \text{M}_{O,ej}^{\text{AGB}} \right), & i = C \\[5pt] 0, & \text{otherwise}, \end{cases}\end{equation}
where
$\delta_i^{\text{AGB}}$
is the fixed condensation efficiency of element i in AGB stars based on Ferrarotti & Gail Reference Ferrarotti and Gail2006. The dust mass produced by AGB stars with carbon-to-oxygen ratios less than one (
$C/O \lt 1$
) is expressed as
\begin{equation} \text{M}_{i,d}^{\text{AGB}} = \begin{cases} 0, & i = C \\[5pt] 16 \sum\limits_{i=\text{Mg,Si,S,Ca,Fe}} \delta_i^{\text{AGB}} \text{M}_{i,ej}^{\text{AGB}}, & i = O \\[14pt] \delta_i^{\text{AGB}} \text{M}_{i,ej}^{\text{AGB}}, & \text{otherwise}, \end{cases}\end{equation}
where
$\mu_i$
is the mass of element i. Finally, the mass of dust produced by Type II SNe is described as
\begin{equation} \text{M}_{i,d}^{\text{SNII}} = \begin{cases} \delta_C^{\text{SNII}} \text{M}_{C,ej}^{\text{SNII}}, & i = C \\[5pt] 16 \sum\limits_{i=\text{Mg,Si,S,Ca,Fe}} \delta_i^{\text{SNII}} \text{M}_{i,ej}^{\text{SNII}}, & i = O \\[14pt] \delta_i^{\text{SNII}} \text{M}_{i,ej}^{\text{SNII}}, & \text{otherwise}, \end{cases}\end{equation}
where
$\delta_i^{\text{SNII}}$
is the fixed condensation efficiency of element i for SNe II based on Bianchi & Schneider Reference Bianchi and Schneider2007.
Results of our sample selections. Data with prevalent volume limitations have been removed and labelled as “
$\ldots$
”.
Dust grains are then formed through the accretion of local gaseous metals following Dwek (Reference Dwek1998).
\begin{equation} \left(\frac{d\text{M}_d}{dt}\right)_{grow} = \left(1 - \frac{\text{M}_d}{\text{ M}_{metal}}\right)\left(\frac{\text{M}_d}{\tau_{accr}}\right),\end{equation}
where
$\text{M}_{metal}$
is the total mass of dust and local gas-phase metals. The accretion timescale
$\tau_{accr}$
is then found following Hirashita (Reference Hirashita2000) and Asano et al. (Reference Asano2013) and is
\begin{equation} \tau_{accr} = \tau_{ref}\left(\frac{\rho_{ref}}{\rho_g}\right)\left(\frac{T_{ref}}{T_g}\right) \left(\frac{Z_{\text{M}_\odot}}{Z_g}\right),\end{equation}
where
$p_g$
,
$T_g$
, and
$Z_g$
are the local gas density, temperature, and metallicity, respectively, and the others are reference values.
$p_{ref}=100$
H atoms cm
$^{-3}$
,
$T_{ref} = 20$
K, and
$\tau_{ref} = 10$
Myr.
These dust grains are eventually destroyed by thermal sputtering following Popping et al. (Reference Popping2017) and McKinnon et al. (Reference McKinnon2017), with the timescale expressed as
\begin{align} \tau_{sp} &= a\left|\frac{da}{dt}\right|^{-1} \sim (0.17\text{Gyr})\left(\frac{a}{0.1\mu\text{m}}\right)\left(\frac{10^{-27}\text{g cm}^{-3}}{\rho_g}\right) \notag \\[5pt] &\times \left[\left(\frac{T_0}{T_g}\right)^{w} + 1\right].\end{align}
where
$w=2.5$
controls the low-temperature scaling of the sputtering and
$T_0=2\times10^6\text{K}$
is the temperature above where the sputtering curve starts to flatten. The growth rate of the dust due to this sputtering is then calculated by
\begin{equation} \left(\frac{d\text{M}_d}{dt}\right)_{sp}=-\frac{\text{M}_d}{\tau_{sp}/3}.\end{equation}
SN blasts are not resolved in the simulations and are implemented by a subgrid model for dust destruction by SN shocks following Dwek & Scalo (Reference Dwek and Scalo1980), Seab & Shull (Reference Seab and Shull1983), and McKee et al. (Reference McKee1987); McKee (Reference McKee1989). The timescale
$\tau_{de}$
is
\begin{equation} \tau_{de} = \frac{\text{M}_g}{\varepsilon\gamma \text{M}_s},\end{equation}
where
$\text{M}_g$
is the local gas mass,
$\varepsilon=0.3$
is the efficiency of local gain destruction,
$\gamma$
is the local SNe II rate, and
$\text{M}_g$
is the mass of the local gas shocked at
$100\text{km s}^{-1}$
. Finally, a solar abundance of
$Z_{\odot} = 0.0134$
is assumed for the star formation and grain growth models taken from Asplund et al. (Reference Asplund2009).
The parameters governing these processes are listed in Table 1 (Li et al. Reference Li2019). It is important to note that these parameters are adjusted to align with observations at low redshifts; therefore, the
$z=0$
simulation matches observations by construction rather than prediction.
We utilise the Simba m100n1024 simulations at redshifts
$z = 0.0, 0.1, 0.5, 1.0, 1.5$
for our comparisons. This simulated cosmological cube has a side length of
$100h^{-1}$
Mpc and contains
$1024^3$
dark matter and gas elements. The mass resolution limits of Simba are
$9.6 \times 10^7 \text{M}_\odot$
for dark matter particles and
$1.82 \times 10^7 \text{M}_\odot$
for gas elements. This simulation adheres to the Planck16 concordant cosmology method (Planck Collaboration et al. Reference Collaboration2016).
4. Sample Selection
To start our comparisons between the observations and Simba, we must ensure that our combined observational dataset is comprehensive and addresses any limitations inherent to each dataset. To accomplish this, we applied specific restrictions regarding redshift, dust mass, stellar mass, and SFR to our datasets. These measures helped facilitate fair and accurate comparisons between observations and simulations.
Our implementation of redshift restrictions relies primarily on the snapshots provided by Simba, as described in Section 3. Accordingly, we also must limit our observational datasets to these ranges. A cut of
$z=0.025$
was then applied at each redshift on either side of the initial redshift to ensure consistent comparisons across cosmic time, resulting in a total bin width of
$\Delta z = 0.05$
. Since no data is available before
$z = 0.0$
, we have limited the right-hand side of the data to the specified redshift bin while maintaining the same range.
These datasets also have limitations concerning volume and sensitivity (see Section 6.2). Given these constraints, along with the stellar mass and dust mass resolution limits of Simba described in Section 3, we chose to limit our dust mass to
$10^{6} \lt \text{M}_D [\text{M}_\odot] \lt 10^{9}$
, stellar mass to
$10^{8.59} \lt \text{M}_\odot \lt 10^{11.5}$
, and SFR to
$10^{-2} \lt \text{M}_{\odot/\text{yr}} \lt 10^{2}$
. These selected ranges proved to be the most effective in minimising the limitations of each dataset and ensuring a complete sample across the datasets. We also excluded GAMA at
$z=0.5$
, G10-COSMOS at
$z=0.0$
, and 3D-HST for
$z\lt0.5$
because these volume limitations at these redshifts negatively influenced our results.
Please refer to Table 1 and Figure 2 for a comprehensive overview of our redshift selections and a quantitative comparison of the number of galaxies between observations and simulations.
The dust mass ranges selected for this study. We plot GAMA/G10-COSMOS/3D-HST data as black dots. We plot the selected data with the ranges applied in green dots with a box surrounding them. The red dotted line represents the dust mass volume limits of the surveys. Note that we exclude GAMA data at
$z=0.5$
due to volume-related issues in the selections.
5. Results
After selecting and constraining our samples, we compare the dust content and properties between GAMA/G10-COSMOS/3D-HST and Simba, following the selections outlined in Section 4. We acknowledge that the significant differences in sample size (Table 1) may introduce some bias in our results and could account for some observed differences between the observations and the model.
DMFs from observations and simulations at
$z = 0-1.5$
. Eales et al. (Reference Eales2009) is plotted from data in the range
$0.6 \lt z \lt 1.0$
and Dunne et al. (Reference Dunne2011) and Beeston et al. (Reference Beeston2018) is plotted from data in the range
$0.0 \lt z \lt 0.1$
. Our results are not standardised to the cosmological parameters described in Li et al. (Reference Li2019).
5.1 Dust mass functions
Figure 3 shows the redshift evolution of the DMFs, comparing ours, the simulation and other observational data consisting of data from Dunne et al. (Reference Dunne2011), Clemens et al. (Reference Clemens2013), and Beeston et al. (Reference Beeston2018) at
$z=0.0$
and Eales et al. (Reference Eales2009) at
$z=1.0$
. We use these findings, presented in Li et al. (Reference Li2019), as an additional complement to our results. We also do not standardise our data to their assumed dust mass absorption coefficient and cosmological parameters, as the difference between the observations and simulations is trivial. This may introduce some slight discrepancies, but it should be sufficient for our comparisons.
At
$z=0.0$
, our data agrees well with Simba and broadly follows the other observed data. colour blackSimba data underestimates the DMF in the medium-mass range, where the observational dataset contains about
$\sim36\%$
of the total galaxies observed. Otherwise, it matches quite well in the low-mass end. The model at
$z=0.1$
exhibits similar characteristics with a more significant underestimate in the medium-mass range where a substantial portion of the observational dataset lies (
$\sim50\%$
) and a slight underestimation near the high-mass end, but this range consists only of
$\sim5\%$
of the total galaxies in the observations. The model at
$z=0.5$
is where Simba begins to noticeably deviate from the observational data. Simba significantly underestimates the low-mass end where a significant portion of the observational data lies (
$\sim74\%$
) and slightly underestimates the high-mass end, but similar to
$z=0.1$
, it consists of only about
$\sim5\%$
of the total observational total. This pattern continues in the
$z=1.0$
model, where Simba underestimates dust near the low and high ends again, where the observational data constitutes
$\sim52\%$
and
$\sim9\%$
of the total galaxies, respectively. In contrast, the
$z=1.5$
model largely agrees near the low-mass end and similarly underestimates the high end. However, unlike
$z=0.1$
and
$z=0.5$
, the high end of the observational dataset contains about
$\sim25\%$
of the total galaxies. Please refer to Table 2 for the complete counts of the observational data in each redshift, dust mass, and stellar mass bin. We will explore the stellar and dust mass bins further in the next section.
Galaxy counts for GAMA/G10-COSMOS/3D-HST in different stellar mass, dust mass, and SFR bins.
5.2. Dust mass versus stellar mass
To further investigate the inconsistencies shown in Figure 3, we will examine the evolution of dust with stellar mass throughout cosmic time between the two datasets. We illustrate our comparisons in Figure 4 (log(SFR)
$\leq$
0) and 5 (log(SFR)
$\geq$
0) as hex-bin (simulations) and contour (observations) plots colour-coded with the star formation rates of the datasets and normalised stellar mass in Figure 6 to examine the populations explicitly. The graphs (Figure 4 and 5) are divided into log(SFR)
$\leq$
0 and log(SFR)
$\geq$
0 galaxies and catalogues (Figure 6) to highlight where the differences between the datasets are the most concentrated.
We will start with our results from Figure 4 representing low-SFR galaxies. At
$z=0$
, our data generally agrees well with Simba. However, there is a noticeable mismatch between the observations and the simulations. The observations’ low-SFR dust-poor galaxies are found at the lowest stellar masses (
$10^{8.59}-10^{9.5}$
) consisting of about 74% of the total observed galaxies in this SFR range, while in the simulations, they are more prevalent at medium stellar masses (
$10^{9.5}-10^{10.5}$
). However, compared to Simba, the sample size of GAMA (see Table 1) is significantly smaller, which may explain this difference. In contrast to
$z=0.0$
, at
$z=0.1$
, the data matches closely with the model but misses massive dust-poor galaxies and dust-rich galaxies. The
$z=0.5$
model, on the other hand, resembles the
$z=0.0$
model but offers an improved match in the previously identified mismatched region. Regardless, the model misses a small fraction (
$\lt2\%$
) of dust-rich galaxies at this redshift. At
$z=1.0$
, the model matches closely with the observations but similarly misses some dust-rich galaxies (
$\lt5\%$
). Once again, the
$z=1.5$
model closely matches the observations; however, the observations do not reflect the massive dust-poor galaxies predicted by the simulations.
We will now focus on the results from Figure 5 representing high-SFR galaxies. At
$z=0.0$
, the observational data with galaxies near the minimum SFR closely aligns with the simulated galaxies, yet it misses massive high-SFR galaxies. The data shows a trend similar to that in Figure 4 at
$z=0.1$
, with the simulated data failing to account for dust-rich galaxies of all sizes, which consists of
$\sim11\%$
of the total number of galaxies in this SFR range. Additionally, while the concentration of high SFR in the observations is comparable, it is more concentrated in medium-sized galaxies than in massive ones. Once again, the
$z=0.5$
model closely resembles the
$z=0.0$
data and misses massive high-SFR galaxies. At
$z=1.0$
, the model produces results consistent with those shown in Figure 4. At
$z=1.5$
, the general trend of the observational data aligns with the simulations. However, the model fails to predict the clump of small dust-rich galaxies at this redshift, consisting of
$\sim20\%$
of the total observed galaxies at this SFR range.
A relation between dust mass and stellar mass of Simba and GAMA/G10-COSMOS/3D-HST of quenching galaxies (SFR
$\leq$
0) at
$z=0.0-1.5$
. Observations, represented by contours weighted by star formation rates, are drawn at −8, −4, −2, −0.5, 0.5, 2, 4, and 8. The simulation, Simba, is plotted as hex-bins. Both are colour-coded according to star formation rates.
A relation between dust mass and stellar mass of Simba and GAMA/G10-COSMOS/3D-HST of star-forming galaxies (SFR
$\geq$
0) at
$z=0.0-1.5$
. Observations, represented by contours weighted by star formation rates, are drawn at
$-$
8,
$-$
4,
$-$
2,
$-$
0.5, 0.5, 2, 4, and 8. The simulation, Simba, is plotted as hex-bins. Both are colour-coded according to star formation rates.
Normalised counts of stellar mass from Simba and GAMA/G10-COSMOS/3D-HST. The observational dataset is separated into individual surveys to highlight the distinctions between them.
Figure 6 coincides with a lot of the discrepancies we are seeing in these low- and high-SFR galaxies. At
$z=0.0$
, we can clearly see a peak of the galaxy population in the Simba data around
$10^9-10^{10} \text{M}_\odot$
that aligns perfectly with the large distribution in this range in low-SFR galaxies in Figure 4. This peak continues to
$z=0.5$
and vanishes from
$z=1.0$
onward, matching again what we see in the low-SFR figure. The massive dust-poor and rich low-SFR galaxies missing in the Simba data can also be seen here, as the observational data favours high-mass galaxies, while Simba favours low- to intermediate-mass galaxies. However, this high mass favour in the observations is likely due to the sensitivity limits in GAMA (see Section 6.2). At redshifts
$z=0.5$
to
$z=1.5$
, the population of galaxies reveals a distinct difference between Simba and the observations, as the observations appear to favour and peak at low masses, whereas Simba peaks at more intermediate masses. This aligns perfectly with what we see in Figures 4 and 5.
6. Discussion
Through our analysis comparing observations with simulations and using dust masses derived from Driver et al. (Reference Driver2018) and simulated by Davé et al. (Reference Davé2019), we have identified several significant discrepancies in the modern simulation model Simba. Notably, we find a tension between the dust-rich and dust-poor galaxies modelled in Simba and those observed. In the following sections, we will attempt to address follow-up questions related to our conclusions, including discussing notable limitations of our datasets, comparing our work with similar analyses conducted by others, and outlining potential future research that could enhance the outcome of our results and improve our models.
6.1. Comparison to other works
Our understanding of how dust evolves across cosmic time has evolved rapidly over the past decade, thanks in part to the hard work done to simulate the physical processes that govern its creation, growth, and destruction. To ensure the precision of these simulations, it is crucial to compare leading models, such as Simba, directly with observational data and each other. This comparison helps us determine how effectively each model reproduces the observed data and where future models can improve. In this context, we compare our findings with previous works, focusing on those that have performed analyses similar to ours.
Since our work builds directly on the findings of Li et al. (Reference Li2019), a paper that also compares with the Simba model, it is encouraging that our results align well with theirs. Especially at
$z=0.0$
, we find that our findings agree well with the previous observations, also included in our paper, and also follow a trend similar to the simulations provided by Simba. We also find agreement at
$z=1.0$
with their findings and previous observations from Eales et al. (Reference Eales2009) that Simba underestimates the observed DMF, especially at the higher end. With our additional data, we also found that Simba underestimates a large fraction of dust-poor galaxies at this redshift and similarly at
$z=0.5$
. We believe that this underestimate is probably a byproduct of an inherent resolution limit (Zheng et al. Reference Zheng2021) in Simba, as this is consistent between
$z=0.5$
and
$z=1.5$
.
The M16 (McKinnon et al. Reference McKinnon2016) and L25n256 (McKinnon et al. Reference McKinnon2017) models are hydrodynamical dust models that track dust production, growth, and destruction up to the early universe, similar to Simba. We note that the L25n256 model is the successor of the M16 model that includes thermal sputtering and uses different dust growth parameters following the work by Hirashita (Reference Hirashita2000). At a redshift
$z=1.0$
, the high-mass end of the DMF predicted by the M16 model aligns well with our observations, though it overestimates it at
$z=0.0$
. The M16 model appears to excel at assembling massive dust-rich galaxies, more so than Simba, which may explain its overestimation at
$z=0.0$
. On the other hand, the L25n256 model presents the opposite and matches the high-mass end at
$z=0.0$
but significantly underestimates the DMF at
$z=1.0$
, similar to Simba.
The fiducial dust model in Popping et al. (Reference Popping2017) is a semi-analytical model (SAM) that includes new recipes to track the production and destruction of dust up to the early universe. Despite the difference in how Simba and Popping et al. (Reference Popping2017) model dust across cosmic time, both present similar dust mass functions that generally agree with each other. Due to this agreement, the fiducial model similarly underestimates the DMF at
$z=1.0$
. Unlike Simba, the model presents a slightly worse match at
$z=0.0$
as it overestimates the DMF at the high end compared to our observations. This work also considers the dust-to-stellar mass relation as we do in Figures 4 and 5. However, Popping et al. (Reference Popping2017) presents their comparisons as functions, and although different, their overall shape seems to agree well with our observations at
$z=0.0$
and
$z=1.0$
.
Donevski et al. (Reference Donevski2020) also compares this dust-to-stellar mass relation and uses Simba, similar to ours. However, they focus on earlier redshifts outside the scope of this paper, so we are unable to compare their results directly.
Similarly to the work done by Popping et al. (Reference Popping2017), Triani et al. (Reference Triani2020) uses a SAM model incorporated with their Dusty SAGE model. Their model effectively matches the dust mass functions compared to ours at
$z=0.0$
and shows good agreement with Simba and may even be a better match. They also present a dust mass function at
$z=1.386$
, which underestimates the function compared to our similar redshift results at
$z=1.0$
and
$z=1.5$
. Additionally, their model includes a dust-to-stellar mass relation; however, like Popping et al. (Reference Popping2017), their comparison is based on functions, preventing a direct comparison. However, the shape of their data appears to follow a trend similar to ours.
The Hou et al. (Reference Hou2019) model is a hydrodynamical simulation with a dust enrichment model that takes into account two different grain sizes and accounts for stellar dust production and interstellar dust processing. Unlike these other models, the proposed model disagrees with our observations at all epochs. At both
$z=0.0$
and
$z=1.0$
, the model overproduces dust-rich galaxies and does not align with our observed DMFs in shape. It is possible that they produce too many massive galaxies at all redshifts in this case.
6.2. Limitations
Here, we discuss the limitations of our datasets to attempt to explain the discrepancies observed between Simba and our observations in Section 4.
6.2.1 Observational limitations
As noted in Sections 2.1-2.3, our datasets vary significantly in the FIR coverage, directly affecting the dust mass estimation performed by magphys. The GAMA dataset provides full FIR coverage, while G10-COSMOS offers partial coverage, and the 3D-HST dataset lacks any coverage. In cases where FIR coverage is present, dust mass estimations in G10-COSMOS depend on constraints from total IR emissions (Jin et al. Reference Jin2018; Magnelli et al. Reference Magnelli2024). Meanwhile, dust mass estimations for the 3D-HST dataset depend solely on extinction measurements due to the absence of FIR data in magphys (da Cunha et al. Reference da Cunha2008). These indirect methods for estimating the dust mass can lead to systematic underestimations for G10-COSMOS and 3D-HST. As noted in Driver et al. (Reference Driver2018), the absence of FIR coverage in these surveys adds systematic error to the dust mass estimates produced by magphys, particularly in the earlier redshifts where G10-COSMOS and 3D-HST coverage is prevalent (see Table 6). Future research should concentrate on improving our far-infrared coverage to address this limitation.
Beyond this, these datasets are volume-limited down to specific dust mass, stellar mass, and star formation rate at low redshifts and similarly sensitivity-limited at high redshifts due to the inherent limitations of each survey (see Driver et al. Reference Driver2018, Table 3). These volume and sensitivity limits result in significant systematic uncertainties and limit our sample sizes (see Table 1). For example, volume limits result in a minimum dust mass that correlates with the redshift for
$z \geq 1.5$
, adding restrictions to our sample (see Figure 2). Furthermore, in Figure 6, we can see that GAMA and G10-COSMOS are sensitivity limited at
$z=0.1$
and
$z=1.5$
, respectively, while G10-COSMOS and 3D-HST are volume limited at
$z=0.1$
and
$z=0.5$
, respectively. These limitations also introduce random errors, such as cosmic variance and Poisson errors, as well as errors from Eddington Bias (see Driver et al. Reference Driver2018, Tables 3–6). Improving survey volume and sensitivity is imperative to overcome this limitation, especially in earlier epochs.
Finally, while not as significant as the other limitations, AGN contamination in the data can lead to extreme estimates of stellar masses and star formation rates, with some minor impact on dust estimates. Therefore, it was essential to remove AGN contaminants (described in Sections 2.1, 2.2, and 2.3) to avoid these effects. However, as noted in Driver et al. (Reference Driver2018), the removal process introduces errors (see Tables 4 and 5).
6.2.2 Model caveats
Although robust, the Simba model also has some caveats that affect our comparisons. We also summarise these caveats in three main points. First, the parameters that govern dust production, growth, and destruction are not well constrained, which introduces uncertainties in the model predictions (Li et al. Reference Li2019, Table 1). A combination of these different free parameters from what Li et al. (Reference Li2019) chose may lead to a better match with the observations. Second, the model does not fully incorporate dust physics. It only models the ISM by depleting gas-phase metals and includes inactive dust particles that are only coupled to gas particles and do not change in size over time. However, the latter is more of an issue in a higher resolution simulation. Finally, the simulation of
$100h^{-1}$
Mpc lacks sufficient resolution to resolve the smallest scales of a multiphase ISM accurately. To overcome this resolution limitation, the model was adjusted such that
$\tau_{\text{ref}}$
, the growth timescale of dust in the Simba, was modified to increase the effective gas density and the parameters governing dust destruction and condensation were fixed. Popping & Péroux (Reference Popping and Péroux2022), a paper that compares multiple simulations, including Simba, also suggests that this time scale is too short and that modelling these effective yields requires more parameters than our simulations cannot yet resolve. Thus, constraining the free parameters in a fully resolved multiphase ISM in a higher resolution simulation will allow us to understand small-scale dust physics better and help improve the match to observational data, ushering in the need for better modelling and observations to contain those models. However, despite the presence of these dust-related caveats, we see no indication that these caveats are having a significant effect on our results.
7. Conclusion
We have presented dust comparisons between our observational data and the simulation model Simba. The Simba model exhibits rough agreement with our observations in DMFs and dust mass versus stellar mass across all epochs, but its limitations are apparent. From our comparisons, we concluded the following:
-
(i) Simba misses dust-rich galaxies for all epochs above
$z=0.0$
(Figure 3–5). This is consistent in low and high SFR and is likely a byproduct of the limitations of the Simba model at later redshifts and a byproduct of survey limitations at earlier redshifts (see Sections 6.2.1 and 6.2.2). -
(ii) Simba has a higher concentration of intermediate-mass low-SFR galaxies when compared to observations (Figure 4–6). This is consistent in redshifts up to
$z=0.5$
. At these redshifts, Simba may not be able to accumulate sufficient stellar mass to form more massive dust-rich galaxies as Simba galaxies seem to clump around
$10^{9}-10^{10}\text{M}_\odot$
. -
(iii) Simba does not accurately model low-dust mass galaxies at earlier redshifts, specifically at
$z=0.5$
and
$z=1.0$
. We believe that this is a limitation in Simba presenting at these redshifts (Figure 6).
Overall, with our current comparisons, we do not see specific indications that Simba has problems related to its implementations of dust physics, despite the limitations of the dust model itself. Rather, we predominantly see that issues arise from differences in galaxy populations, leading to the observed discrepancies in dust mass. Some of these issues could be solved with a higher resolution simulation, especially between
$z=0.5$
and
$z=1.5$
, but later redshifts indicate a problem within the Simba model itself. As for the observational data, significant survey limitations impact some of our results, so future observations should focus on improving survey volume and FIR coverage. Several surveys promise to improve statistics here, such as WAVES for spectroscopic redshifts, Rubin Observatory’s Deep Drilling Fields combined with Euclid Deep Field observations, and the ALMA-wide survey on COSMOS, which covers only the brightest galaxies and will improve statistics there.
Acknowledgements
This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. Reference Collaboration2013; Astropy Collaboration et al. Reference Collaboration2018).
Open access
