sígame v3: Gas Fragmentation in Postprocessing of Cosmological Simulations for More Accurate Infrared Line Emission Modeling

We use the 3D radiative transfer code Skirt (Camps & Baes 2020) to calculate the local dust-attenuated ISRF. Skirt uses information on the dust and stellar distribution in the galaxy simulation to inform a Monte Carlo algorithm on how to emulate the relevant physical processes including scattering, absorption, and emission of photons as they pass through the interstellar dust. For the stellar component ("SourceSystem" in Skirt), an IMF similar to that used in the galaxy simulation is chosen, and the stellar emission is set to match the Binary Population and Spectral Synthesis (BPASS v2.2.1) models according to the age and metallicity of each star particle (Eldridge et al. 2017; Stanway & Eldridge 2018).

Skirt requires radiative smoothing lengths around each stellar particle's position to spread out the launch locations of the photon packet around that position. If such a smoothing length is not available in the simulation, sígame will calculate one that scales linearly with stellar mass and spans 100–300 pc. Often, Skirt divides the dust medium into cells of sizes much smaller than 100 pc, but we note that the gravitational smoothing lengths of cosmological simulations are typically much larger (of the order of 1 kpc). Therefore, although SKIRT inserts a lot of cells to ease the radiative transfer calculation, adjacent cells often have similar densities. The choice of scaling the radiative smoothing length to between 100 and 300 pc was made after having tested two extreme cases: (a) stellar smoothing lengths of 10 pc and (b) stellar smoothing lengths of 1000 pc. In case (a) the stellar radiation field ends up being concentrated in point sources and in case (b) it creates artificial "rings" around the original stellar particle position. Hence, we chose the smoothing lengths somewhere in between. ¹⁹

Finally, Skirt requires the initial mass of each stellar particle rather than the current stellar mass, since the spectral energy distributions used by Skirt already take into account the mass evolution of the population based on its age and metallicity. If not available in the simulation, sígame will calculate the initial stellar mass, before mass loss due to stellar evolution, using the publicly available Python-FSPS code, which itself is a Python translation of the Flexible Stellar Population Synthesis code, in order to convert current stellar mass and formation time into initial stellar masses (Conroy et al. 2009; Conroy & Gunn 2010; Foreman-Mackey et al. 2014). For the dust component ("MediumSystem" in Skirt), sígame can either use the dust masses directly, if provided by the galaxy simulation, or calculate them from the metallicity of the gas using a fixed DTM mass ratio. For the dust types and composition, we use the built-in THEMIS model (Jones et al. 2017) for dust in the diffuse ISM that can be invoked in Skirt. This dust mixture contains two families of dust particles: amorphous silicate and amorphous hydrocarbon (we refer to Jones et al. 2017 for more detail). The postprocessing with Skirt returns a spectrum for each cell of each galaxy on an oct-tree adaptive grid. The rest of sígame works on this grid structure, for which gas densities are calculated with the SWIFTsimIO package (Borrow & Borrisov 2020), and all other properties are inherited from the nearest fluid element or gas particle in the galaxy simulation.

In addition to being attenuated by dust grains, interstellar light from stars is also absorbed by gas, in particular in the hydrogen-ionizing regime at photon energies above 13.6 eV. With Skirt a spectrum is generated in each cell that has been attenuated by dust only, which tends to absorb radiation in the FUV at 6–13.6 eV (Gnedin et al. 2008). In order to account for the additional attenuation by gas, we run a suite of Cloudy models for a fixed luminosity source, over a range of total hydrogen column densities. For the luminosity source, an unobscured stellar spectrum of a stellar population of age 1 Myr and solar metallicity, scaled to a luminosity of 10⁶ L_⊙, using the same BPASS v2.2.1 tables as those used in the modeling with Skirt. The stellar spectrum is added to a background continuum spectrum corresponding to the observed local ISM ISRF as available with the "table ism" command in Cloudy. For these models, we keep the metallicity at 1 Z_⊙ and maintain the dust size distribution and abundance at values appropriate for the ISM of the Milky Way as given by the Cloudy "grains ISM" command. The latter includes both a graphitic and a silicate component and generally reproduces the observed overall extinction properties for a ratio of extinction per reddening of RV ≡ A_v /E(B − V) = 3.1 (we refer to the Cloudy manual for details). We use the local ISM abundance table of Cloudy, which is a mean of the warm and cold phases of the ISM from Cowie & Songaila (1986) and Savage & Sembach (1996). In these Cloudy models, and all other Cloudy models in this work, the cosmic microwave background at z = 0 is included as an additional blackbody radiation field. By changing the column density at which Cloudy stops the calculation, a set of transmitted spectra of different shapes are produced containing the effects of different amounts of dust and gas absorption as shown in Figure 2.

Zoom In Zoom Out Reset image size

From the shape of the spectra in Figure 2, we quantify the effect of dust absorption as the ratio between optical-to-near-infrared (OIR; 0.5–2 eV) flux and FUV flux (red and teal shaded regions in Figure 2). We can now match the derived OIR/FUV flux ratio with that of Skirt in each simulation cell to translate the dust extinction into an average column density of gas that the stellar light has passed through, and thereby complete the shape of the incident spectrum (not the normalization) to be used for the final Cloudy grid described in Section 2.4. The inclusion of absorption by gas in this way is not an exact solution because we are essentially approximating the effect of several luminosity sources and column densities taken into account by Skirt with a single luminosity source and one column of gas and dust in Cloudy. We do this because there is no clear way of recovering the exact path of each photon packet ejected by Skirt. Furthermore, this method assumes a constant DTM ratio (fixed at roughly 0.5 here) and uses a single metallicity for all the intervening gas. Finally, we note that the chosen stellar population sets the X-ray portion of the spectra in Cloudy, which is further attenuated by gas and dust as seen in Figure 2. Yet this part of the spectrum could look very different for a stellar population of different age and metallicity, and the intensity of X-rays could potentially affect the emission line flux (see Vallini et al. 2018 for a study on CO emission from high-z active galactic nuclei (AGNs)).

Due to the limited resolution of cosmological simulations, the individual cells typically do not contain densities higher than ∼10 cm⁻³ (see Appendix C for the gas and SFR density distributions in one of the simulated galaxy samples used to test sígame in the following sections). However, all the FIR emission lines considered here can or will only be excited at much higher densities, as can be seen from the critical densities listed in Table 1. In order to perform the dust radiative transfer, Skirt introduces cells of smaller size than those in the source simulation, but this procedure does not increase the density of the gas.

To mitigate the lack of resolution in density, this version of sígame turns toward simulations made at much higher spatial resolution in order to use them as look-up tables that can help describe the fragmentation of the gas on sub-grid scales. In principle, the user can import their simulation of choice, but for the purpose of testing sígame in this paper we chose the high-resolution data from a simulation performed with the Voronoi mesh code arepo (Springel 2010) of an interacting M51-like galaxy model (Tress et al. 2020). The simulation (hereafter described as arepo-M51) reached subparsec resolution in dense gas and allowed for an analysis of the formation and destruction of GMCs, which showed that the evolution of GMCs depends only weakly on galaxy–galaxy interactions.

2.3.1. The arepo-M51 Simulation

2.3.2. Parameterizing the Density PDF

The arepo-M51 simulation contains starbursting regions as well as regions with hardly any ongoing star formation. The question we pose is: how does the gas density PDF change from one region to another and can we parameterize it for use in a cosmological simulation of much lower resolution? To investigate this, we divide the arepo-M51 simulation volume into cubes of 200 pc on a side and calculate the gas density PDF within each cube. The region size of 200 pc was chosen in order to represent typical Skirt cell sizes, and at the same time be large enough to contain enough elements in arepo-M51 to properly sample the gas density PDF. We expect the PDF to move toward higher densities as the volume-averaged gas density of a region increases, but also as SFR density increases (Kainulainen et al. 2009). Figure 3 shows the resulting mean PDFs of cells within 200 pc regions using teal dashed lines. In gray we show the 1σ spread for 12 bins of volume-averaged density, 〈n_H〉_V , and SFR density, 〈n_SFR〉_V . A red vertical dashed line in each panel indicates the critical density, above which the cell data start to be incomplete as some gas cells are converted into sink particles.

Zoom In Zoom Out Reset image size

Due to the spread in PDFs, as indicated with gray areas in Figure 3, we cannot just interpolate between PDFs. Instead, we make a parametric fit, such that the fit parameters depend on 〈n_H〉_V and 〈n_SFR〉_V . For the parametric fit, we build on previous work showing that the density distribution of a collapsing cloud can be approximated by a lognormal with a power-law tail. As discussed in the Introduction, both observations and simulations show that the density field of an isothermal star-forming cloud is well approximated by a lognormal distribution at low density and a power-law distribution at high density. Burkhart (2018) describes the resulting PDF as

$$\begin{eqnarray}{p}_{\mathrm{LN}+\mathrm{PL}}(s)=\left\{\begin{array}{cl}N\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{s}}{e}^{\tfrac{{\left(s-{s}_{0}\right)}^{2}}{2{\sigma }_{s}^{2}}}, & s\lt {s}_{t}\\ {{Ne}}^{-\alpha s}, & s\gt {s}_{t},\end{array}\right.\end{eqnarray} \tag{ 1 }$$

where N is the normalization constant, σ_s is the standard deviation of the lognormal, α is the power-law slope, s₀ is the mean logarithmic overdensity, and s is the logarithmic overdensity

$$\begin{eqnarray}&&s\equiv \mathrm{ln}(\rho /{\rho }_{0}),\end{eqnarray} \tag{ 2 }$$

where ρ₀ is the volume-averaged density of the entire cloud. The transition density, s_t , above which the distribution approximates a power law, was shown by Burkhart (2018) to relate to the lognormal width σ_s and power-law slope α as

$$\begin{eqnarray}&&{s}_{t}=(\alpha -1/2){\sigma }_{s}^{2}.\end{eqnarray} \tag{ 3 }$$

In Figure 3 a fit is made to all mean PDFs by combining a lognormal function with a power law at the high-density end, keeping σ_s and α as free parameters.

However, this fit does not take into account gas mass locked in sink particles. The gas in sink particles where no SN has exploded will be dense and distributed in relatively undisturbed GMCs, while gas in other sink particles with supernova remnants (SNRs) has already been at least partly dispersed in the surrounding ISM due to SN feedback. In order to account for the dense gas mass in sink particles, we divide the sink particles into two categories:

• 1.  
  Category I sink particles contain dense, cold gas that creates stellar populations with time, but no SN explosions. Here, we expect the gas to be relatively undisturbed and we count all gas mass in the particle as dense gas following a power-law density PDF regardless of the sink particle age. The percentage of sink particles in this category is only 0.13%.
• 2.  
  Category II sink particles contain dense, cold gas that has formed or will produce at least one SN. Due to feedback from the massive stars that eventually produce SN explosions, a GMC will not survive for long in this environment. In a study looking for OH(1720 MHz) masers in the inner Galaxy, Hewitt & Yusef-Zadeh (2009) found that dense gas interacts with only 15% of SNRs. Considering that SNRs are visible only during the first ∼0.1 Myr of their lifetime and assuming that SNe occur at a roughly constant rate after the first SNe go off, we can convert the SNR number fraction of 15% into an upper age limit. We count gas in sink particles with SNRs as being dense gas only if the sink particle age since the first SN exploded is less than 15% of the longest SNR lifetime. The first SNe, created by the most massive stars, are expected to begin exploding at a sink particle age of ∼3 Myr (e.g., Yungelson et al. 2008). The percentage of sink particles in this category is 22.7%.

We count the internal gas mass of the remaining sink particles as diffuse gas that would already have been dispersed by feedback if ionization and winds were included. As sink particles accrete more mass than forms into stars, there is a maximum mass limit of ∼ 2 × 10⁵ M_⊙ above which they are not allowed to accrete further and instead form a new sink particle. This retains nominal resolution in the collisionless particles instead of having single particles with anomalously high mass. Figure 4 illustrates the selection criteria and the distribution in sink particle ages and total current gas masses (total sink mass – sink stellar mass). The older the sink particles, the more likely it is that stars have formed and the gas has been used up or dispersed. This outweighs the additional mass from accretion over time. A main branch can be seen of sink particles that start out with the maximum mass and begin to reduce in gas mass after a period of 3 Myr when SNe start to disperse the gas (Category II sink particles by the definition above). Another branch is of sink particles that ran out of gas, many of which are Category I particles that will never produce an SN.

Zoom In Zoom Out Reset image size

The additional sink gas mass is expected to change only the high-density portion of each PDF, because it has already passed the accretion density threshold of ρ_c = 600 cm⁻³. We first attempt to accommodate the sink gas mass by iterating on the power-law slope until enough additional mass has been added. This method results in the orange dashed curves in Figure 3, while the orange solid curves are combined lognormal and power-law fits made to the new distributions.

The resulting fits, shown with orange curves in Figure 3, are used as look-up tables for each cell in the Skirt output, by interpolating in 〈n_H〉_V and 〈n_SFR〉_V . The corresponding sonic Mach numbers, M_s , for each fit can be estimated from σ_s through (Padoan et al. 1997)

$$\begin{eqnarray}&&{\sigma }_{s}^{2}=\mathrm{ln}(1+{b}^{2}{M}_{s}^{2}).\end{eqnarray} \tag{ 4 }$$

Adopting a turbulent forcing parameter of b = 1/3, our density PDFs display Mach numbers from 7 to 27. For cell densities below the minimum grid point in the arepo-M51 look-up tables, we adopt a single, narrow lognormal PDF corresponding to a Mach number of 1 and b = 1/3 (Federrath et al. 2008) without the power-law tail. All PDFs are cropped to the density range from n_H = 10⁻⁴ cm⁻³ to n_H = 10⁷ cm⁻³.

The final step of sígame v3 is to derive the line emission, having gathered all the necessary information by postprocessing the simulation. To this end, we use a library of one-zone Cloudy models that span the parameter ranges listed in Table 2. As mentioned in Section 2.2, the column densities ($\mathrm{log}\,$ N_H) in Table 2 are not actual column densities in the simulations, but a parameter used to set the shape of the spectra. For each value of the F_FUV flux, the cosmic-ray ionization rate ξ is scaled as

$$\begin{eqnarray}&&\xi =({F}_{\mathrm{FUV}}/{G}_{0}){\xi }_{0},\end{eqnarray} \tag{ 5 }$$

where ξ₀ = 10⁻¹⁶ s⁻¹ is taken as the average Milky Way value (Indriolo et al. 2007) and the solar neighborhood FUV flux G₀ = 1.6 × 10⁻³ erg cm⁻² s⁻¹ (Habing 1968). Exploring all combinations of the parameters listed in Table 2 leaves us with 129,600 one-zone Cloudy models.

Table 2. Parameters for the One-zone Cloudy Models

Model Parameter	Min. Value	Max. Value	Step Size
${\rm{log}}({n}_{{\rm{H}}}/{{\rm{cm}}}^{-3})$	−4	7	1
${\rm{log}}(Z/{Z}_{\odot })$	−2	0.5	0.5
${\rm{log}}({N}_{{\rm{H}}}/{{\rm{cm}}}^{-2})$	17	22	0.5
${\rm{log}}({F}_{{\rm{FUV}}}/{G}_{0})$	−7	4	1
${\rm{log}}$(DTM ratio)	−2	−0.2	0.2

Download table as:  ASCIITypeset image

The one-zone Cloudy models must now be sampled according to the gas density PDFs found in Section 2.3, such that each combination of 〈n_SFR〉_V and 〈n_H〉_V in the galaxy simulations corresponds to a density PDF-weighted sum of the 12 different one-zone Cloudy models along the $\mathrm{log}{n}_{{\rm{H}}}$ axis that correspond to that region's other parameters ($\mathrm{log}Z$, $\mathrm{log}{N}_{{\rm{H}}}$, $\mathrm{log}{F}_{\mathrm{FUV}}$, and $\mathrm{log}\,$ DTM). For this step, six grid points in 〈n_H〉_V (from 10⁻⁴ to 10² cm⁻³) and four grid points in 〈n_SFR〉_V (from 10^−2.5 to 10^0.5 M_⊙ yr⁻¹) are used. For each combination of 〈n_H〉_V and 〈n_SFR〉_V we take the gas density PDF that comes closest in terms of both values, and shift the center of the lognormal to exactly match 〈n_H〉_V . This shift is performed to ensure that the PDFs generated are not only centered on the six chosen 〈n_H〉_V grid points, but can fill out the density space and generate a smooth total PDF as shown in Figure 7. The sampling of one-zone Cloudy models is carried out for all combinations of N_H, F_FUV, Z, and DTM ratio, leading to a look-up table of 259,200 different combinations of one-zone Cloudy models, with one such table for each spectral line considered.

2.4.1. User-defined Sub-grid Attenuation Function

We apply sígame to a simulated galaxy sample from the simba cosmological galaxy formation simulation (Davé et al. 2019). The simba simulations are run using the meshless finite-mass hydrodynamics technique in the gizmo N-body plus hydrodynamics code (Hopkins 2015, 2017). The main simba run evolves 1024³ gas elements and 1024³ dark matter particles within a 100h⁻¹ Mpc volume (simba-100) from z = 249 → 0. To improve the dynamic range and test resolution convergence, we also use a higher-resolution run with 512³ gas elements and 512³ dark matter particles within a 25h⁻¹ Mpc volume (simba-25).

simba includes a range of sub-grid models for galaxy formation physics, including H₂-based star formation, two-phase kinetic galactic winds, torque-limited and Bondi black hole accretion, three types of AGN feedback, and a sub-grid model to form and destroy dust during the simulation run (Li et al. 2018, 2019). The galaxy properties in simba have been compared to various observations across cosmic time and shown to provide reasonable agreement (e.g., Rodríguez Montero et al. 2019; Thomas et al. 2019; Appleby et al. 2020; Lower et al. 2020). For more details we refer to Davé et al. (2019).

We extracted samples of galaxies from the z = 0 snapshots of simba-100 and simba-25 using the caesar galaxy and halo catalog generator, which identifies galaxies using a six-dimensional friends-of-friends algorithm applied to dense gas (n > 0.13 cm⁻³) and stars (Thompson 2015). Before making any selection cuts, a total of 49,215 and 2463 galaxies were found by caesar for simba-100 and simba-25, respectively. To ensure we have a reasonably well-resolved gas distribution within the galaxy, we only consider galaxies with >300 gas elements in both simulation boxes, corresponding to gas masses of at least 5.6 × 10⁹ M_⊙ and 0.7 × 10⁹ M_⊙ in simba-100 and simba-25, respectively. For comparison the mass resolution of the arepo-M51 simulation, which will be used to fragment the gas as described in Section 2.3, is a few solar masses. From these we select a test sample of 400 galaxies from simba-100 and another sample of 240 galaxies from simba-25. The samples were selected to span a wide range in stellar mass, SFR, and gas mass. Only galaxies found to have nonzero mean SFR over the past 100 Myr were included.

Figure 5 shows the positions of the simba-100 and simba-25 galaxy samples in the SFR–M_⋆ space. The distribution of all simba-100 and simba-25 z = 0 galaxies identified with caesar and within the axis limits is shown underneath with logarithmic hexbin contours. The agreement with observations by Salim et al. (2018) is generally good, though simba overestimates the SFR over the range M_⋆ ∼ 10^9.5–10¹⁰ M_⊙ (see discussion of this in Davé et al. 2019). Overall, simba reproduces the observed SFR–M_* distribution fairly well, in terms of both the main sequence and the quenched fractions (Davé et al. 2019), as well as the bimodality in specific SFR (Katsianis et al. 2021).

Zoom In Zoom Out Reset image size

For our simba-100 sample, hydrogen densities range from 0.5 to 66 cm⁻³ and SFR-weighted mean metallicities 〈Z〉_SFR span 0.1–1.9 Z_⊙. By weighing the metallicity by SFR, we are approximating the metallicity that would be observed using nebular emission lines since these lines primarily trace star-forming regions. Global galaxy properties such as stellar mass, gas mass, SFR, 〈Z〉_SFR, and radius can be found in Table 4 together with their line luminosities in Table 5 in Appendix A and online.

Figure 6 illustrates the outcome of modeling and propagating the stellar light with Skirt in step 1 of sígame (see Section 2.1). Skirt outputs the radiation in cells with sizes ranging from 19 pc to 6.3 kpc, and a mass-weighted distribution in cell sizes that peaks at ∼200 pc (for the simba-100 galaxies). For the purposes of this paper, we found that a total of 10⁸ photons per galaxy was enough to give stable results, corresponding to more than ∼2500 photons per star particle. See Appendix B for a test increasing the number of photon packets to 10⁹ that resulted in a negligible change in total FUV luminosity and FUV flux distribution. For the dust component ("MediumSystem" in Skirt), we directly use the dust masses that are calculated on-the-fly in simba (Li et al. 2019). The resulting total infrared (3–1100 μm) luminosities span from 2.41 × 10⁷ L_⊙ to 2.73 × 10¹¹ L_⊙ and FUV (6–13.6 eV) luminosities from 5.78 × 10⁷ L_⊙ to 1.09 × 10¹⁰ L_⊙ using the simba-100 dust masses directly. The mass-weighted metallicities of our model galaxies range from 0.02 to 1.9 Z_⊙ for the simba-100 sample, so we set the metallicity of intervening gas to 1 Z_⊙ in step 2 of sígame (see Section 2.2).

Zoom In Zoom Out Reset image size

After step 3 of sígame (see Section 2.3), we can construct galaxy-wide density PDFs, and examples for the same galaxies as in Figure 6 are shown in Figure 7. The original density PDF in the cell data from Skirt (solid line) is compared to the derived PDF adopting a single lognormal with Mach number 10 and forcing parameter b = 1/3 versus the approach described in Section 2.3 in dotted and dashed lines, respectively.

Zoom In Zoom Out Reset image size

Figure 8 shows moment0 maps of CO(1–0), [C ii], and [O iii]88 emission constructed for the same three galaxies from simba-100 shown in Figure 6, using the output from completing step 4 of sígame (see Section 2.4). The output data cube from sígame contains spatial information, cell sizes, and cell luminosities that were combined to derive the surface brightness of each pixel in the final image, weighting each cell by its volume filling factor within the column covered by each pixel. ²⁰ We note that the calculation of the moment0 maps assumes an ISM fully transparent to the FIR line emission, an assumption that might not hold for shorter wavelengths. Some expected differences in how the emission lines trace the ISM can be observed, such as the preference for [C ii] to arrive from PDRs near star-forming regions compared to the broader CO(1–0) emission, and the relatively concentrated [O iii]88 emission restricted to H ii regions experiencing harder radiation fields from young stars. The more concentrated [C ii] and [O iii]88 emission relative to the CO(1–0) is also reflected in slightly smaller half-light radii compared to the CO(1–0) maps as illustrated with green dashed circles.

Zoom In Zoom Out Reset image size

In order to compare different assumptions in the sub-grid process, we perform four different runs on the same z ∼ 0 sample from the simba-100 box as well as one run on galaxies selected in the simba-25 box. The names of the different runs are listed in Table 3. In the first run (100Mpc_M10 in Table 3), we report the line luminosities that sígame returns when the density fragmentation is carried out the "traditional way" with a single lognormal density PDF corresponding to a Mach number 10 and a forcing parameter of 1/3. A value of 10 for the Mach number is supported by observations of clouds in the solar neighborhood with typical sound speeds of 0.2–0.3 km s⁻¹ and velocity dispersions of several km s⁻¹ (see, e.g., Goldreich & Kwan 1974; Kainulainen & Tan 2013; Hennebelle & Inutsuka 2019). Next, we applied sígame v2 to the same set of galaxies, using the modifications to the code described in Leung et al. (2020) together with a new Cloudy grid at z = 0 (100Mpc_SIGAMEv2). The next simulation run adopts the new gas fragmentation scheme described in this paper, and will be referred to as our "default run" (100Mpc_arepoPDF). The default settings are also applied to a smaller sample of 240 galaxies in the simba-25 box as a convergence test (25Mpc_arepoPDF). Finally, we also investigate the resulting line emission when not including the additional sub-grid attenuation function as described in Section 2.4.1 and otherwise adopted in the default run (100Mpc_arepoPDF_no_ext).

We will be comparing the different simulation runs with a broad selection of nearby observed galaxies. Since the samples of simulated galaxies were selected to span a wide range in stellar mass, SFR, and gas mass, they were not tailored to match any specific observed sample. We will be comparing the range in simulated and observed line luminosity as a function of SFR, looking for agreements on the order of magnitude and deferring a more careful comparison with observations to a future study. The sample of observed galaxies comprises all nearby galaxies observed with Herschel or the Infrared Space Observatory that we could find in the literature for which SFR estimates could also be made. In order to best compare with observations, we estimate the SFR of the model galaxies as a mean over the past 100 Myr, and not the instantaneous SFR of gas particles. Furthermore, all SFRs derived using a Salpeter (1955) IMF have been converted to the equivalent SFR of a Chabrier (2003) IMF as in the simba simulation, by adopting a −0.24 dex correction (e.g., Mitchell et al. 2013). The observed SFRs were obtained with various methods summarized here. For the sample of mixed-type galaxies in Kamenetzky et al. (2016) and the luminous infrared galaxies (LIRGs) in Díaz-Santos et al. (2013), the FIR luminosities, L_FIR (40–120 μm and 42.5–122.5 μm in the two respective cases), of those papers are converted into SFRs according to Kennicutt (1998). For the sample of dwarf galaxies of Cormier et al. (2015), SFRs come from Rémy-Ruyer et al. (2015), who derived SFR from FUV and Hα luminosities, corrected for dust attenuation. The galaxy sample of Schruba et al. (2012) contains CO(2–1) luminosities for 16 dwarf galaxies (with one galaxy in common with the Cormier et al. 2015 sample) as well as CO(1–0) luminosities for nearby galaxies compiled from the literature, of which we show 22 CO(1–0) observations from the HERA CO-Line Extragalactic Survey (HERACLES; Leroy et al. 2009) and 20 from Calzetti et al. (2010). For the 16 dwarf galaxies, SFRs in Schruba et al. (2012) were calculated as a combination of FUV and 24 μm emission following the approach in Bigiel et al. (2008) and Leroy et al. (2008). For the SFRs of the 50 nearby galaxies we use the values compiled from literature in Schruba et al. (2012), and refer to the references therein for specific methods. Finally, for the sample of main-sequence and starburst galaxies of Brauher et al. (2008), we convert the reported 63 μm flux to an SFR using the 70 μm monochromatic SFR conversion relation of Calzetti et al. (2010).

The comparison between simulation runs relative to observed line luminosity–SFR relations is made in Figure 9. The figure shows the offset in line luminosity from a power-law fit made to the observed line luminosity as a function of SFR for each of the eight emission lines considered here. The standard deviation of the offsets for the observed galaxies is indicated with gray bars in the background. The different sígame runs listed in Table 3 can now be compared for each of the eight emission lines considered here, with the fine-structure lines sorted in order of increasing critical density (Table 1). Due to the mixed sample of observed galaxies that we compare to, we define a "good agreement" here as when the quartile distribution of simulated line luminosities overlaps with the 2σ spread around the observed line luminosity–SFR relation.

Zoom In Zoom Out Reset image size

For run 100Mpc_M10, the resulting deviations from observations are shown with orange bars. Surprisingly, this simple assumption does a very good job for the ionized species of [N ii]122, [N ii]205, [O iii]88, and [C ii], for which good agreement is reached, but it completely underestimates CO emission. The CO(1–0), (2–1), and (3–2) lines fall below the observed line–SFR relation by median deviations of 3.4, 3.9, and 3.5 dex, respectively. The underestimation of the molecular lines is not surprising, given that the density PDF does not extend significantly beyond densities >10² cm⁻³, where molecules form (see Figure 7).

With run 100Mpc_SIGAMEv2 (purple bars), we test how well the previous sígame framework, which has so far only been applied at z ∼ 2 and above, works for the simba-100 z = 0 sample. While this version of the code gives a good agreement with observations for [C ii], [O i]63, CO(1–0), and CO(3–2), it significantly underestimates the [N ii] and [O iii]88 lines with median deviations larger than 1 dex. For all lines except [O i]63, this version also produces outliers with extremely low line luminosities, as shown with open circles (and indeed, the plot does not extend to show all outliers). We attribute this disagreement to the simplified gas fragmentation scheme used in sígame v2, in which gas particles were split into either GMCs following a locally observed mass spectrum or diffuse gas, with no intermediate ISM phase and no sanity checks on the density PDF.

The default run of the current version of sígame, 100Mpc_arepoPDF (green bars) comes closest to all line luminosities simultaneously, with the fewest outliers. In total, six out of the eight lines considered here are in good agreement with the observed relations, while the [O i]63 and CO(3–2) line luminosities are overestimated by median deviations of 1.3 and 1.0 dex, respectively, compared to the observations. We do not consider the overestimation of CO(3–2) an outstanding problem due to the lack of observations compiled here (our observed comparison sample comprises 21 galaxies).

The test run without additional extinction, 100Mpc_arepoPDF_no_ext (brown bars), looks very similar to the default run in Figure 9, with the most noticeable difference being an increased [O i]63 line prediction that falls above the observed luminosities by a median deviation of 1.9 dex. In our default run, with the user-defined attenuation function described in Section 2.4.1, the overestimation of [O i]63 is brought down by about 0.5 dex to 1.3 dex above the observed [O i]63–SFR relation, and the tail of low CO luminosities is reduced in length by several dex. This suggests that the inclusion of additional attenuation is necessary in order for enough [O i] and CO to form. An effect on the [N ii] and [O iii]88 lines could also be expected since these lines originate in dense H ii regions. Although the current version of sígame does not explicitly model H ii regions (see the discussion in Section 4.3), the additional attenuation at densities above 10² cm⁻³ should also affect these lines, but this is not clear from the distributions in Figure 9 alone. In order to investigate the effect of additional attenuation, Figure 10 compares the default run (100Mpc_arepoPDF) to that without additional attenuation (100Mpc_arepoPDF_no_ext) as a function of galaxy properties M_⋆, M_gas, SFR, and 〈Z〉_SFR. Only the CO and [O i]63 line luminosities are consistently higher and lower, respectively, in the default run, in agreement with Figure 9. The deviations in the [N ii], [O iii]88, and [C ii] line luminosities are negligible at the lower end of M_⋆, M_gas, SFR, and 〈Z〉_SFR, but tend toward negative values at the higher end. In particular, the [N ii] and [O iii]88 lines are more strongly affected by the additional attenuation for high values of SFR and 〈Z〉_SFR (rightmost two panels). This behavior can be explained by considering what happens in the one-zone Cloudy models of which our look-up table consists. At higher metallicities, increasing self-shielding of nitrogen and oxygen results in more singly ionized nitrogen and doubly ionized oxygen relative to higher ionization states since the one-zone Cloudy models include self-shielding. This in turn results in higher [N ii] and [O iii]88 luminosities with increasing metallicity for a given one-zone Cloudy model and hence a strong dependence on metallicity in run 100Mpc_arepoPDF_no_ext. However, once additional attenuation is added at densities above 100 cm⁻³, the luminosity of those one-zone models is negligible as there is no ionizing radiation to ionize nitrogen and oxygen. The difference between the two runs is therefore more pronounced at high metallicity, where models with densities above 100 cm⁻³ have had their [N ii] and [O iii]88 luminosities dramatically reduced. Since there is a weak correlation between SFR and 〈Z〉_SFR in our sample, the same can be said for the evolution with SFR seen in Figure 10.

Zoom In Zoom Out Reset image size

Finally, the test run with default settings at higher mass resolution, 25Mpc_arepoPDF (blue bars), shows quartile ranges that overlap with those of the simba-100 sample, indicating that the sígame results do not depend significantly on mass resolution in the underlying simulation. However, this analysis alone does not reveal any bias due to the different galaxy properties of the two samples. In Figure 14 in Appendix D we investigate in more depth how well the two simulation runs agree with one another for a range of physical galaxy properties, finding that for similar galaxies, the luminosities of [N ii]122, [N ii]205, [O iii]88, [C ii], and [O i]63 tend to be overestimated by up to ∼0.5–1.5 dex in simba-25 compared to simba-100, in particular at high SFR and 〈Z〉_SFR. It is interesting to note that 100Mpc_arepoPDF_no_ext behaves more similarly to 25Mpc_arepoPDF than to 100Mpc_arepoPDF, and it is hard to say whether this is a resolution issue or due to the attenuation function present in the latter two runs.

We can compare our [C ii]–SFR relation with that found through similar techniques and hence better understand the decisive differences between the models. Figure 11 shows the [C ii]–SFR relations found by sígame for the simba-100 galaxies with the 100Mpc_arepoPDF run (green contour lines) together with the relations derived by Popping et al. (2019) and Ramos Padilla et al. (2021). From the 100Mpc_arepoPDF run, we also show the distribution of only main-sequence (MS) galaxies, chosen to be between the 16th and 84th percentiles of the Salim et al. (2018) relation (purple dashed contour lines). The distribution of sígame simulated galaxies generally matches that of the observed galaxies well, but fails to reproduce some of the galaxies of lower [C ii] luminosity and the nondetections of Brauher et al. (2008). Selecting only MS galaxies does not change this picture significantly. Our simulated galaxies lie above the simulated sample of Popping et al. (2019), but in apparent extension of the samples of Ramos Padilla et al. (2021), made with cosmological simulations from the Evolution and Assembly of Galaxies and their Environments (EAGLE) project (Crain et al. 2015; Schaye et al. 2015). In the latter study, GMCs are created in postprocessing following Olsen et al. (2015), and Cloudy is used to derive line intensities as in newer versions of sígame. One interesting avenue of further study would be to apply sígame to EAGLE, to fully understand how intrinsic differences between EAGLE and simba affect this comparison.

Zoom In Zoom Out Reset image size

Another interesting comparison can be made to the CO line modeling of Vallini et al. (2018, hereafter V18) in which the sub-grid density PDF was determined using the Mach number of the underlying simulation together with a time-dependent evolution of the high-density power-law tail. V18 find good agreement with the observed relation between GMC virial mass and CO luminosity in nearby galaxies, the CO excitation in nearby starburst galaxies and the CO excitation of a low-metallicity GMC in the LMC. In contrast, our 100Mpc_M10 run, which comes closest to the method of V18 does not reproduce the CO line emission well. Two key differences that likely cause this difference are: (1) in V18, the Mach number is allowed to change and is in fact inherited from the underlying simulation, and (2) the underlying simulation used in V18 is very different from simba in the sense that they reach much higher densities before any postprocessing is performed. The density PDF of Althæa (the ramses zoom-in simulation used by V18) peaks at 300 cm⁻³ and reaches densities above 1000 cm⁻³. In our case, simba hardly reaches 100 cm⁻³, and hence using a lognormal of any Mach number (even one higher than 10) does not result in densities high enough to excite CO. It could be said that the approach in V18 works well if the underlying simulation reaches certain minimum densities, but for cosmological simulations, applying a lognormal (+power law) centered at the mean volume-averaged cell density will not be enough on its own.

We have only shown results for the present framework using one type of simulation, namely the arepo suite of M51 realizations. We also examined an arepo simulation of the central molecular zone (CMZ), the arepo-CMZ simulation (Sormani et al. 2020). The arepo-CMZ simulation contains relatively compact clouds with extreme velocity dispersions and a threshold density for the formation of sink particles an order of magnitude higher than that for the M51 simulation, meaning that the ISM is allowed to fragment further before forming sink particles. By binning the density PDFs of the arepo-CMZ simulation in the same way as done for the M51 simulation (Figure 3), we can compare the PDF shapes and the resulting line luminosity for two different simulations. At high volume-averaged gas densities, the mean PDF shapes are similar in the two simulations, but at lower densities ($\mathrm{log}{n}_{{\rm{H}}}\lt 0$), the arepo-CMZ simulation returns PDFs with relatively more mass at higher densities. The resulting line luminosities using the arepo-CMZ PDF table are in general higher (up to 0.27 dex on average for [O iii]88) except for the CO lines, which are lower by up to 0.19 dex on average. In Figure 15 in Appendix E we show the line luminosity for bins in mass-weighted hydrogen density for all cells in the simba-25 galaxy shown in the top panels of Figure 6 with different assumptions for the sub-grid density PDF. We conclude that the galaxy region sampled can be important, although a better test would be to compare with an entirely different galaxy simulation (not just a special region like the CMZ).

Another caveat associated with the approximations used for the arepo-M51 simulations is a tendency for the clouds to be too long-lived, which might in turn yield too little SFR for a given surface density. Of relevance to sígame, this could mean that we are overestimating the number of clouds (i.e., dense gas mass fraction in the PDF) for a given SFR volume density. However, choosing a completely different galaxy simulation technique (e.g., models run with adaptive mesh refinement codes such as ramses, enzo, or athena—Teyssier 2002; Bryan et al. 2014; Kim & Ostriker 2017) might return different PDFs, the effect of which we have not studied. Instead, we have focused on making the current version of sígame flexible and publicly available so that the user can include any high-resolution simulation for the gas fragmentation if they so desire.

As discussed in the previous section, this version of sígame fragments the gas density on scales smaller than the Skirt cell sizes (typically larger than ∼15 pc), thereby allowing for clumps of higher densities on scales not resolved by the parent simulation. In a similar manner, we could also expect the radiation field to be fragmented such that the mass-weighted FUV flux distribution on parsec scales can differ substantially from the overall cell mean flux value and allow for parsec-size features such as H ii regions or shielded, dense, molecular cores. Judging from the overall agreement between simulated and observed lines, we speculate that most of the FIR lines considered here arise from structures of larger extent where parsec resolution in radiation is not necessary, with the exception of [O i]63, which is overestimated. The continued overestimation of the [O i]63 and CO(3–2) line emission by our default run (see Figure 9) suggests that a more careful treatment of the sub-grid attenuation and/or other features remains missing, however. A proper treatment of the flux distribution on these scales would require considerations about the properties of natal clouds around young stars, as demonstrated with the one-dimensional stellar feedback model warpfield (Rahner et al. 2017). In addition, [C ii] and [O iii]88 remain slightly overestimated, both of which originate mainly in the atomic ISM phase in our model, hinting at an overestimation of atomic gas mass versus ionized gas mass in our gas fragmentation scheme.

4.3.1. Missing Treatment of Active Galactic Nuclei

In starburst galaxies and mergers, shocks are expected to act as an additional heating source in the ISM. Although the next version of Cloudy may include a treatment for shock-heated gas, the current version of Cloudy iterates to find a thermal and ionization equilibrium, thereby ignoring any shocked state of the gas that might be present in the simulation. One way to treat shock-heated gas would be to set up a separate grid of models, following the technique used for MAPPINGS III (Allen et al. 2008). Turbulence (2–10 km s⁻¹ in velocity dispersion) and/or density-dependent magnetic fields are also believed to play a role in setting the [O i]63 emission, through their effect on line width, shielding and pumping, and the chemical and thermal state of the gas (Canning et al. 2016). However, testing the effect of these would require Cloudy models of more than one zone where optical depth is taken into account.