SIGAME v3: Gas Fragmentation in Post-processing of Cosmological Simulations for More Accurate Infrared Line Emission Modeling
Karen Pardos Olsen, Blakesley Burkhart, Mordecai-Mark Mac Low, Robin G. Tre?, Thomas R. Greve, David Vizgan, Jay Motka, Josh Borrow, Gergö Popping, Romeel Davé, Rowan J. Smith, Desika Narayanan
DDraft version February 8, 2021
Typeset using L A TEX twocolumn style in AASTeX63 s´ıgame
V3: GAS FRAGMENTATION IN POST-PROCESSING OF COSMOLOGICALSIMULATIONS FOR MORE ACCURATE INFRARED LINE EMISSION MODELING
Karen Pardos Olsen, Blakesley Burkhart,
2, 3
Mordecai-Mark Mac Low,
4, 3
Robin G. Treß, Thomas R. Greve,
6, 7, 8
David Vizgan, Jay Motka, Josh Borrow, Gerg¨o Popping, Romeel Dav´e,
12, 13, 14
Rowan J. Smith, and Desika Narayanan
16, 17, 6 Department of Astronomy and Steward Observatory, University of Arizona, Tucson, AZ 85721, USA Department of Physics and Astronomy, Rutgers, The State University of New Jersey, 136 Frelinghuysen Rd, Piscataway, NJ 08854, USA Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA Department of Astrophysics, American Museum of Natural History, New York, NY 10024, USA Universit¨at Heidelberg, Zentrum f¨ur Astronomie, Institut f¨ur theoretische Astrophysik, Albert-Ueberle-Str. 2, 69120 Heidelberg, Germany Cosmic Dawn Center (DAWN) DTU-Space, Technical University of Denmark, Elektrovej 327, DK-2800 Kgs. Lyngby, Denmark Dept. of Physics & Astronomy, University College London, Gower Place, WC1E 6BT London Department of Astronomy, Wesleyan University, Middletown, CT 06457, USA Institute for Computational Cosmology, Department of Physics, University of Durham, South Road, Durham DH1 3LE, UK European Southern Observatory, 85478 Garching bei M¨unchen, Germany Institute for Astronomy, Royal Observatory, Univ. of Edinburgh, Edinburgh EH9 3HJ, UK Department of Physics and Astronomy, University of the Western Cape, Robert Sobukwe Road, Bellville 7535, South Africa South African Astronomical Observatories, Observatory, Cape Town 7925, South Africa Jodrell Bank Centre for Astrophysics, University of Manchester, Oxford Road, Manchester M13 9PL, UK Department of Astronomy, University of Florida, 211 Bryant Space Sciences Center, Gainesville, FL 32611 USA University of Florida Informatics Institute, 432 Newell Drive, CISE Bldg E251, Gainesville, FL 32611
Submitted to ApJABSTRACTWe present an update to the framework called SImulator of GAlaxy Millimeter/submillimeter Emis-sion ( s´ıgame ). s´ıgame derives line emission in the far-infrared (FIR) for galaxies in particle-basedcosmological hydrodynamics simulations by applying radiative transfer and physics recipes via a post-processing step after completion of the simulation. In this version, a new technique is developed tomodel higher gas densities by parametrizing the gas density probability distribution function (PDF)in higher resolution simulations for use as a look-up table, allowing for more adaptive PDFs than inprevious work. s´ıgame v3 is tested on redshift z = 0 galaxies drawn from the simba cosmologicalsimulation for eight FIR emission lines tracing vastly different interstellar medium phases. Includingdust radiative transfer with Skirt and high resolution photo-ionization models with
Cloudy , thisnew method is able to self-consistently reproduce observed relations between line luminosity and starformation rate in all cases, except for [N II ]122, [N II ]205 and [O I ]63, the luminosities of which areoverestimated by median factors of 1 .
6, 1 . . (cid:46)
200 pc).
Keywords: galaxies, ISM, simulation, modeling, radiative transfer
Corresponding author: Karen Pardos [email protected] INTRODUCTIONThe physical and chemical state of the interstellarmedium (ISM) plays a key role in determining the starformation rate (SFR) of a galaxy and is therefore of im-mense importance for understanding galaxy evolution. a r X i v : . [ a s t r o - ph . GA ] F e b Olsen et al.
Macroscopic galactic events like starbursts and mergerscreate pressure and density waves that determine whereand how the giant molecular clouds (GMCs) form inwhich stars and planets grow (e.g. Colombo et al. 2014;Sun et al. 2018; Chevance et al. 2020; Alves et al. 2020).Resolved observations of nearby clouds and detailedsimulations that track individual cores in GMCs frominitial perturbations to final collapse, can probe theactual density distribution above the highest densityachievable in cosmological simulations. Both observa-tions and numerical work to this effect have determinedthat the dense gas in molecular clouds (with extinctions A v > n > cm − ) typically has a den-sity PDF with a power-law tail at the high end set bygravitational collapse (Klessen 2000; Kainulainen et al.2009; Collins et al. 2012; Girichidis et al. 2014; Burkhartet al. 2015; Schneider et al. 2015; Lombardi et al. 2015;Burkhart et al. 2017; Mocz et al. 2017; Padoan et al.2017; Chen et al. 2018; Alves et al. 2017). The lowdensity portion of the density probability density func-tion (PDF) might follow a log-normal shape but this ishighly uncertain due to selection effects, completeness,and the influence of stellar feedback (Alves et al. 2017;Chen et al. 2018).On galactic scales, it becomes too computationally ex-pensive for simulations to model the fragmentation andcollapse of each cloud, not to mention the subsequentstellar feedback processes. Instead, parametrizationsare used as for example in the cosmological simulations simba (Dav´e et al. 2019) and IllustrisTNG (Weinbergeret al. 2017; Pillepich et al. 2018). In simba , an H -basedSFR recipe is used where the amount of H is estimatedfrom the metallicity and local column density followingthe sub-grid model of Krumholz & Gnedin (2011). InIllustrisTNG, gas with density n (cid:38) . − is allowedto form stars in accordance with the empirically definedKennicutt-Schmidt relation. Both simulations lack theresolution to track individual stars, and instead assumea Chabrier (2003) initial mass function for stellar popu-lations formed.Synthetic observations of simulated galaxies are ourmost promising way of directly comparing models withobservations of the ISM. Observations of the ISM typi-cally target its most effective cooling channels, namelyfar-infrared (FIR) emission lines. These were discov-ered by the Kuiper Airborne Observatory (Watson et al.1984; Stacey et al. 1991; Lord et al. 1996). The subse-quent Infrared Space Telescope (Kessler et al. 1996) and
Herschel Space Observatory ( Herschel hereafter; Pil-bratt et al. 2010) showed that the strongest emissionlines come from a handful of species, namely neutral andionized oxygen, singly ionized carbon, ionized nitrogen as well as molecular lines such as the CO rotational lines.Together, these species probe all phases of the ISM fromGMCs to photon-dominated regions (or photodissocia-tion regions; PDRs), warm neutral medium and hot ion-ized gas. Existing telescopes like ALMA, NOEMA andSOFIA as well as upcoming space and balloon missions(e.g. GUSTO and ASTHROS) provide an ever-growingdatabase of emission line observations of galaxies at allredshifts, requiring detailed modeling of the same linesfor their correct interpretation. By taking a snapshotof a hydrodynamic simulation containing gas, stars anddark matter, line emission can be calculated as a post-processing step for each cell or particle in the simulationusing physically motivated recipes.All efforts to simulate emission lines have had to de-termine small-scale values for: 1) the structure and localspectral shape of the radiation field, in particular in thefar-ultraviolet (FUV), 2) the structure of gas density,and 3) the local chemistry and level populations. Dueto computational constraints, most of these are not cal-culated with enough precision in the simulation itself,and must hence be derived in post-processing based onsub-grid models. Below, we summarize the current ap-proaches:1. For the FUV, the current practice typically in-volves adopting the solar neighborhood FUV inter-stellar radiation field strength G (Habing 1968)and scaling it by the global SFR of the simulatedgalaxy compared to that of the Milky Way (e.g.Vallini et al. 2019) to get a galaxy-averaged FUVflux; or by using the local SFR density to get a lo-cal kpc-scale FUV flux (e.g. Olsen et al. 2017; Pop-ping et al. 2019). Pallottini et al. (2019) demon-strated the use of on-the-fly radiative transfer for ahigh-resolution adaptive mesh refinement simula-tion of a z ∼ In this work, we introduce a post-processingradiative transfer calculation for large-box simula-tions at z = 02. Retrieving the high density ( n > cm − ) gasduring post-processing of a cosmological simula-tion is crucial to model several FIR lines, as canbe seen from the critical densities of typical FIRlines listed in Table 1. This step is typically doneby adopting a clumping factor that artificially in-creases the H production (e.g. Narayanan et al.2008; Gnedin et al. 2009; Dav´e et al. 2016; Lupi & IGAME v3
Inthis work, we employ simulations of higher spatialresolution as look-up tables and use analytic re-lations of density PDFs from Burkhart (2018) toinfer the PDF on sub-grid scales in a cosmologicalsimulation.
3. Once the radiation field and density are deter-mined, there are many tools publicly available thatsimultaneously solve for ionization state, temper-ature and line excitation, with some of them beingtailored to a specific ISM phase. For an overviewof widely used tools to solve for the chemistryand generate line emission in the ISM, see Olsenet al. (2018), which also discusses the limitationsof each approach.
In this work, we use the photo-ionization code
Cloudy (Ferland et al. 2017) toset the shape of the spectrum locally as well as cal-culate chemistry and line emission.
As previously mentioned, in addition to inferring thePDF in higher resolution simulations as described inpoint 2 above, we also use a theoretical framework toestimate the amount of self-gravitating gas not resolvedin the galaxy-sized simulation. In order to alleviate ten-sions with observations and dense gas models that usea log-normal PDF, Burkhart & Mocz (2019) derive adense gas and star formation model based primarily onthe presence of a power law tail. The work provides ananalytic expression connecting the transitional columndensity value between log-normal and power-law to thewidth and of the log-normal and the slope of the power-law which we will use here (see also Collins et al. 2012;Burkhart et al. 2017; Burkhart 2018).With this paper, we present a new version of ourpost-processing framework Simulator of Galaxy Millime-ter/submillimeter Emission ( s´ıgame ; Olsen et al. 2017)to model FIR/mm line emission in particle-based hy-drodynamic simulations, and test it on a large sample
Table 1.
Critical densities of important ISM cooling linesin the FIR. CO critical densities were calculated usingthe spontaneous emission coefficients from the LAMDAdatabase as accessed on 21 Aug 2020 (Sch¨oier et al. 2005),using a temperature of 100 K and typical collision crosssection of 10 − cm − . The [C II ] critical densities arefrom Goldsmith et al. (2012), calculated at a temperatureof 500 K. The remaining critical densities are from Maddenet al. (2013) at 300 K.Line n crit (cm − ) Origin[C II ] 16 a , 2 . × b ,4 . × c All ISM, but mainlyPDRs[O I ]63 4 . × b PDRs[O
III ]88 510 a H II regions, radi-ation dominated byOB stars[N II ]122 310 a Ionized ISM[N II ]205 48 a Ionized ISMCO(1-0) 650 c Molecular ISMCO(2-1) 6 . × c Molecular ISMCO(3-2) 2 . × c Molecular ISM a For collisions with electrons. b For collisions with hydrogen atoms. c For collisions with H2 molecules. of emission lines for comparison with observations. Sec-tion 2 gives a quick overview of the new structure of s´ıgame and lists the updates made since the previousversion. The test sample of galaxies used for the testsare described in Section 3. Section 4 provides more de-tails on the steps that s´ıgame goes through. Section 5compares and validates the code by comparison with ob-servations and the previous version of s´ıgame . Finally,Section 6 discusses the results and provides caveats forusing this method and we conclude in Section 7.Throughout this paper, we adopt a flat cold darkmatter cosmology with cosmological parameters Ω Λ =0 . m = 0 . b = 0 . h = 0 .
678 (Planck Collaboration et al.2016). S´IGAME VERSION 3 OVERVIEW s´ıgame derives FIR/mm line emission of simulatedgalaxies by post-processing a particle-based cosmologi-cal hydrodynamics simulation. Figure 1 provides a quicklook at the overall structure of s´ıgame version 3. Thissection will summarize the updates made since the lastversion of the code, with more details below (Sect. 4).
Olsen et al.
Stellar and dust info
Simulation data
Gas: [x,y,z,m,n H ,Z,m dust ]Stars: [x,y,z,m,Z,age] NIR/FUV flux ratio F
FUV
SKIRT
Radiative transfer with dustabsorption
Analysis of datacubesN H Cloudy models of N H (NIR/FUV) Gas fragmentation with AREPO n H n H , Z Z F FUV N H Z F
FUV N H Cloudy models of line emission (2)(1)(3) (4)
Figure 1.
Flow chart illustrating the structure of SIGAME v3.
Main updates
Leung et al. (2020) adapted version 2 of s´ıgame toapply it to the simba simulation for the first time tostudy the [C II ] luminosity function at z = 6. With thisversion 3, the main updates are: • Radiative transfer of stellar FUV emission throughdust: s´ıgame v2 used a library of starburst99 (Leitherer et al. 2014) synthesis models for stellarpopulations to derive the local FUV field, with-out taking into account the potentially importanteffect of dust absorption. In s´ıgame v3 a fulldust radiative transfer calculation is performed us-ing the three-dimensional (3D) radiative transfercode
Skirt (Camps & Baes 2020) and the actualdust-to-metal ratios from the simba simulation (Liet al. 2019). • Fragmentation of the gas on sub-grid scales: Inprevious versions of s´ıgame , gas particles fromcosmological simulations were divided into a dif-fuse and a dense gas phase, the latter being fur-ther divided into spherical GMC-like structureswith masses > M (cid:12) . With the new version of s´ıgame , gas particles are first distributed on an adaptive grid with Skirt and then fragmented ac-cording to the local gas and SFR volume densitiesusing the output from a high resolution simulation,and thereby avoiding the assumption of sphericalclouds. • High resolution thermo-chemistry of the ISM: Asin previous versions of s´ıgame , we use the photo-ionization code
Cloudy (this time version 17.02)to calculate the chemical and thermal state of thegas and the resulting line emission. However, be-cause of the decision to fragment the gas to smallermass and size scales, we have moved to using the grid command feature of
Cloudy , allowing fora look-up table of higher resolution. THE z = 0 GALAXY SAMPLEWe apply s´ıgame to a simulated galaxy sample fromthe simba cosmological galaxy formation simulation(Dav´e et al. 2019). The simba simulations are run usingthe meshless finite mass hydrodynamics technique in the gizmo N-body plus hydrodynamics code (Hopkins 2015,2017). The main simba run evolves 1024 gas elementsand 1024 dark matter particles within a 100 h − Mpcvolume ( simba -100) from z = 249 →
0. To improve thedynamic range and test resolution convergence, we alsouse a higher-resolution run with 512 gas elements and IGAME v3 Figure 2.
The position on the SFR-stellar mass plane of theselected test sample at z = 0 of 400 simba -100 galaxies (cir-cles) and 240 simba -25 galaxies (triangles), with both sam-ples color-coded by SFR-weighted metallicity. For compari-son, the grey hexbin contours show the overall distributionof z = 0 star-forming galaxies in simba -100 after discardinggalaxies with no star formation according to caesar or fewerthan 300 gas particles. The observed z = 0 main sequencefrom two studies are also shown (Salim et al. 2018; Speagleet al. 2014) with shaded regions referring to 16th and 84thpercentiles in the case of Salim et al. (2018) and 0.2 dex forSpeagle et al. (2014). dark matter particles within a 25 h − Mpc volume( simba -25). simba includes a range of sub-grid models for galaxyformation physics, including H -based star formation,two-phase kinetic galactic winds, torque-limited andBondi black hole accretion, three types of AGN feed-back, and a sub-grid model to form and destroy dustduring the simulation run (Li et al. 2018, 2019). Thegalaxy properties in simba have been compared to var-ious observations across cosmic time and shown to pro-vide reasonable agreement (e.g. Thomas et al. 2019; Ap-pleby et al. 2020; Lower et al. 2020; Rodr´ıguez Monteroet al. 2019). For more details we refer to Dav´e et al.(2019).We extracted samples of galaxies from the z = 0snapshots of simba -100 and simba -25 using the cae-sar galaxy and halo catalog generator, which identifiesgalaxies using a six-dimensional friends-of-friends algo-rithm applied to dense gas ( n > .
13 cm − ) and stars(Thompson 2015). Before making any selection cuts, atotal of 20,428 and 2,463 galaxies were found by caesar for simba -100 and simba -25, respectively. To ensure wehave a reasonably well-resolved gas distribution withinthe galaxy, we only consider galaxies with >
300 gas el- ements in both simulation boxes, corresponding to a gasmass of at least 5 . × M (cid:12) and 0 . × M (cid:12) in simba -100 and simba -25, respectively. From these we select atest sample of 400 galaxies from simba -100 and anothersample of 240 galaxies from simba -25. The samples wereselected to span a wide range in stellar mass, SFR andgas mass. Only galaxies found to have non-zero meanSFR over the past 100 Myr were included.Figure 2 shows the positions of the simba -100 galaxysample in the SFR- M (cid:63) space. The distribution of all16,617 simba -100 z = 0 galaxies identified with caesar and having non-zero mean SFRs is shown underneathwith hexbin contours. The agreement with observationsby Salim et al. (2018) is generally good, though simba overestimates the SFR over the range M (cid:63) ∼ . –10 M (cid:12) (see discussion of this in Dav´e et al. 2019).Overall, simba reproduces the observed SFR- M ∗ dis-tribution fairly well, both in terms of the main sequenceand the quenched fractions (Dav´e et al. 2019), as well asthe bimodality in specific SFR (Katsianis et al. 2020).For our sample of 400 galaxies, hydrogen densitiesrange from 0.5 to 66 cm − and SFR-weighted meanmetallicities (cid:104) Z (cid:105) SFR span 0.1 to 2.4 Z (cid:12) . By weighingthe metallicity by SFR, we are approximating the metal-licity that would be observed using nebular emissionlines since these lines primarily trace star-forming re-gions. Global galaxy properties such as stellar mass,gas mass, SFR, (cid:104) Z (cid:105) SFR , and radius can be found in Ta-ble 4 together with their line luminosities in Table 5 inAppendix A and online.We note that the version of s´ıgame as presented inthis paper is set up to take the snapshot informationfrom simba . In the near future, we plan to adapt thecode so that it can also accept output from higher reso-lution, non-cosmological, zoom particle-based, as well asmesh-based simulations, but for the purpose of testingthe new version, we use the cosmological particle-basedformat to which s´ıgame has previously been applied. S´IGAME V3The new s´ıgame framework consists of four majorsteps, each indicated on Figure 1 with numbered boxes.The following four subsections describe the algorithmsthat are applied within each step.4.1.
Dust attenuation with
Skirt (step 1)
Figure 3 illustrates the outcome of modeling and prop-agating the stellar light with the 3D radiative trans-fer code
Skirt (Camps & Baes 2020).
Skirt requiresinformation on the dust and stellar distribution, mostof which can be extracted from the simba snapshotsthemselves.
Skirt then employs the Monte Carlo tech-
Olsen et al.
Figure 3.
Examples of three simba -25 galaxies going through step 1 of s´ıgame , where stellar light is propagated, attenuatedand re-emitted by dust grains using the radiative transfer code
Skirt . Left column : Surface density maps of the total gas massin simba . Center column : Positions and ages of all stellar population particles with ages below 1 Gyr. The sizes of the stellarcircles scale with stellar population mass.
Right column : Maps of observed FUV (6–13.6 eV) flux derived with
Skirt as observedfrom a distance of 10 Mpc, showing increased FUV flux in areas close to young stellar populations and attenuation by regionsof dense gas.
IGAME v3 pho-tons per galaxy was enough to give stable results. SeeAppendix B for a test increasing the number of photonpackets to 10 that resulted in a negligible change in to-tal FUV luminosity and FUV flux distribution. In thefollowing, we provide more detail on the specific inputparameters used for and the operation of Skirt .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). For the dust types andcomposition, we use the built-in THEMIS model (Joneset al. 2017) for dust in the diffuse interstellar mediumthat can be invoked in
Skirt . This dust mixture con-tains two families of dust particles; amorphous silicateand amorphous hydrocarbon (we refer to Jones et al.2017 for more detail).For the stellar component (“SourceSystem” in
Skirt ),we adopt a Chabrier (2003) IMF as also used in simba .A smoothing length for each stellar particle must be pro-vided to
Skirt , but as simba lacks this parameter, weassumed a smoothing length that scales linearly withstellar mass and spans 100 to 300 pc. We note that theexact smoothing length prescription for the stars has lit-tle effect on the final flux distribution. Finally,
Skirt requires the initial mass of each stellar particle and notthe current stellar mass as reported in the output from simba , since the spectral energy distributions used by
Skirt already take into account the mass evolution ofthe population based on its age and metallicity. Theinitial stellar mass, before mass loss due to stellar evo-lution, is calculated using the publicly available Python-FSPS code, which itself is a python translation of theFlexible Stellar Population Synthesis code, in order toconvert current stellar mass and formation time into ini-tial stellar masses (Conroy et al. 2009; Conroy & Gunn2010; Foreman-Mackey et al. 2014).The resulting total infrared (3–1100 µ m) luminositiesspan from 2 . × L (cid:12) to 2 . × L (cid:12) and FUV (6–13.6 eV) luminosities from 5 . × L (cid:12) to 1 . × L (cid:12) using the simba -100 dust masses directly. The post-processing with Skirt returns a spectrum for each cellof each galaxy on an oct-tree adaptive grid. The restof s´ıgame works on this grid structure, for which gasdensities are calculated with the
SWIFTsimIO pack-age (Borrow & Borrisov 2020), and metallicity, dust-to-mass (DTM) ratio, element abundances, and SFR vol-ume density are inherited from the nearest fluid elementin the simba simulation. 4.2.
Attenuation by intervening gas with
Cloudy (step 2)
In addition to being attenuated by dust grains, in-terstellar light from stars is also absorbed by gas, inparticular in the hydrogen-ionizing regime at photon en-ergies above 13 . Skirt a spectrum is gener-ated in each cell that has been attenuated by dust only,which tends to absorb radiation in the FUV at 6–13 . Cloudy models over a range of total hydrogen column densi-ties, keeping the metallicity at 1 Z (cid:12) and maintainingthe dust size distribution and abundance to that appro-priate for the ISM of the Milky Way as given by the
Cloudy “grains ISM” command. The latter includesboth a graphitic and silicate component and generallyreproduces the observed overall extinction properties fora ratio of extinction per reddening of RV ≡ Av/E(B-V)= 3.1 (we refer to the
Cloudy manual for details). Weuse the local ISM abundance table of
Cloudy , which isa mean of the warm and cold phases of the ISM fromCowie & Songaila (1986) and Savage & Sembach (1996).In these
Cloudy models, and all other
Cloudy mod-els in this work, the CMB at z = 0 is included as anadditional blackbody radiation field. For the incidentcontinuum, the local ISRF in Cloudy is used. The re-sulting transmitted spectra are shown in Figure 4. Thespectrum is not attenuated significantly for column den-sities below around 10 cm − , and above this columndensity the incident spectra are attenuated both aboveand below 13 . N H > cm − . Assuming that the attenuation in theNIR-to-FUV flux ratio is primarily due to dust, we cannow match the derived NIR/FUV flux ratio with thatof Skirt in each simulation cell to translate the dustextinction into a column density of gas that the stellarlight has passed through. This method comes with twoapproximations: 1) it assumes a constant DTM ratio,which is inconsistent with simba where dust abundancesare tracked separately, 2) it requires the metallicity ofthe intervening gas, but since we do not know the exactpath of each photon, we fix the metallicity to solar. Themass-weighted metallicities of our model galaxies rangefrom 0.03 to 0.74 Z (cid:12) (SFR-weighted from 0.1 to 1.9 Z (cid:12) ),so we believe this to be a reasonable approximation.
Olsen et al.
Figure 4.
Transmitted spectra for a set of
Cloudy models with hydrogen density 1 cm − , solar dust and metallicity abundancesand a range in total hydrogen column density shown in different colors. The NIR and FUV photometric bands that we useare indicated with shaded regions and a vertical dashed line highlights the 13.6 eV ionization potential of hydrogen. Significantattenuation only occurs at column densities N H > . cm − . Below that value the spectra become indistinguishable from theinput spectrum. Gas fragmentation (step 3)
Due to the limited resolution of cosmological simula-tions, the individual cells do not contain densities higherthan ∼
10 cm − (see Appendix C). However all the FIRemission lines considered here can or will only be excitedat much higher densities, as can be seen from the criticaldensities listed in Table 1. All high density clumps andfilaments on scales below the typical minimum smooth-ing scale of ∼ simba galaxies are unresolved.In order to perform the dust radiative transfer, Skirt introduces cells of smaller sizes, with cell sizes rangingfrom 19 pc to 6 . simba -100 galaxies) anda mass-weighted distribution in cell sizes that peaks at ∼
200 pc, but this procedure does not increase the den-sity of the gas.To mitigate the lack of resolution in density, this ver-sion of s´ıgame turns towards simulations made at muchhigher spatial resolution in order to use them as look-up tables that can help describe the fragmentation ofthe gas on sub-grid scales. The high-resolution dataset chosen comes from a simulation performed with themoving-mesh code arepo (Springel 2010) of an inter-acting M51-like galaxy model (Tress et al. 2020). Thesimulation (hereafter described as arepo -M51) reachedsub-parsec resolution in dense gas and allowed for ananalysis of the formation and destruction of GMCs andhow the evolution of GMCs depends only weakly ongalaxy-galaxy interactions.4.3.1.
The arepo -M51 simulation
For details on the simulation setup we refer to Tresset al. (2020), summarizing here only the components ofrelevance to this project. The arepo -M51 simulationwas designed to be able to resolve and study GMCsin the context of a galaxy interaction. The simula- tion setup includes a time-dependent, non-equilibrium,chemistry network tracking hydrogen and CO chemistry,and follows star formation and feedback processes reach-ing sub-parsec resolutions at densities n (cid:38)
100 cm − (see Fig. 3 of Tress et al. 2020). The calculations wereexecuted using arepo , a moving-mesh magnetohydro-dynamic (MHD) code coupled to an oct-tree gravitysolver, from which only the hydrodynamic capabilitieswere used (Weinberger et al. 2020). At each time-stepthe code constructs the Voronoi grid given a set of meshgenerating points that are constrained to move followingthe local velocity of the fluid. This pseudo-Lagrangianmoving-mesh technique is naturally adaptive, allowingfor high dynamic range in spatial scales down to thesubstructure of GMCs.The galaxy model and the interaction with a compan-ion were adjusted to resemble the M51 system and thechemical network followed the cooling and self-shieldingof molecular gas from foreground interstellar radiation.This allowed the formation of GMCs where runaway col-lapse leads to star formation (SF). Dense, gravitation-ally bound, and collapsing gas is accreted onto collision-less sink particles that abstract the last stages of the SFprocess to a sub-grid model due to the limited resolutionof even this higher resolution simulation.At gas densities ρ c > − g cm − ( n H (cid:38)
600 cm − )the gas is tested to determine if it is bound and collaps-ing and if so a sink particle is created. Dense gas with ρ > ρ c that comes within an accretion radius of 2.5 pcwill then be accreted by the sink particle if bound toit. Given these densities and sizes, the sink particlesrepresent small stellar (sub)-clusters. At these scalesSF is still fairly inefficient, so only 5 % of the accretedgas mass is converted into stars. In this sense everysink particle consists of a stellar component and a gas IGAME v3 arepo -M51 simulation include the lack of early stellarfeedback such as winds, jets, and ionising radiation andthe absence of magnetic fields.4.3.2.
Parametrizing the density PDF
The arepo -M51 simulation contains star-bursting re-gions as well as regions with hardly any ongoing starformation. The question we pose is; how does the gasdensity PDF change from one region to another and canwe parametrize it for use in simba ? To investigate this,we divide the arepo -M51 simulation volume into cubesof 200 pc on a side and calculate the gas density PDFwithin each cube. The region size of 200 pc was cho-sen in order to represent typical
Skirt cell sizes, andat the same time be large enough to contain enough ele-ments in arepo -M51 to properly sample the gas densityPDF. We expect the PDF to move towards higher den-sities as the volume-averaged gas density of a region in-creases, but also as SFR density increases (Kainulainenet al. 2009). Figure 14 in Appendix D shows the result-ing mean PDFs of cells within 200 pc regions using tealdashed lines. In grey we show the 1 σ spread for 12bins of volume-averaged density, (cid:104) n H (cid:105) V , and SFR den-sity, (cid:104) n SFR (cid:105) V . A red vertical dashed line in each panelindicates the critical density, above which the cell datastarts to be incomplete as some gas cells are convertedinto sink particles. Due to the spread in PDFs, as indicated with greyareas in Figure 14, we cannot just interpolate betweenPDFs. Instead, we make a parametric fit, such thatthe fit parameters depend on (cid:104) n H (cid:105) V and (cid:104) n SFR (cid:105) V . Forthe parametric fit, we build on previous work showingthat the density distribution of a collapsing cloud canbe approximated by a log-normal with a power-law tail.Simulations of gravo-magneto-hydrodynamic turbulentfragmentation have shown that the density field of anisothermal star-forming cloud is well approximated by alog-normal distribution at low density and a power-lawdistribution at high density (Klessen 2000; Collins et al.2012; Burkhart et al. 2015; Mocz et al. 2017; Padoanet al. 2017; Burkhart 2018): p LN + P L ( s ) = N √ πσ s e ( s − s σ s , s < s t N e − αs , s > s t , (1)where N is the normalization constant, s is the meanlogarithmic overdensity, and s is the logarithmic over-density: s ≡ ln( ρ/ρ ) (2)where ρ is volume-averaged density of the entire cloud.The transition density, s t , above which the distribu-tion approximates a power-law, was shown by Burkhart(2018) to relate to the log-normal width, σ s , and power-law slope, α , as: s t = ( α − / σ s (3)In Figure 14 in Appendix D, a fit is made to all meanPDFs by combining a log-normal function with a power-law at the high-density end, keeping σ s and α as freeparameters. However, this fit does not take into ac-count gas mass locked in sink particles. The gas in sinkparticles where no SN has exploded will be dense anddistributed in relatively undisturbed GMCs, while gasin other sink particles with supernova remnants (SNR)has already been at least partly dispersed in the sur-rounding ISM due to SN feedback. In order to accountfor the dense gas mass in sink particles, we divide thesink particles in two categories: • Category I sink particles have stellar populationsthat produce no SN explosions. Here, we expectthe gas to be relatively undisturbed and we countall gas mass in the particle as dense gas followinga power-law density PDF regardless of the sinkparticle age. The percentage of sink particles inthis category is only 0 .
13 %. • Category II sink particles have formed or willproduce at least one SN. Due to feedback from0
Olsen et al. the SN explosion, a GMC will not survive forlong in this environment. In a study looking forOH(1720 MHz) masers in the inner Galaxy, He-witt & Yusef-Zadeh (2009) found that dense gascan survive accompanying 15% of SNRs. Consid-ering that SNRs are visible only during the first ∼ . . yror ∼ . ∼ × M (cid:12) and beginto reduce in gas mass after a period of roughly 10 . yrwhen SNe start to disperse the gas. The sharp cutoff ininitial gas mass at 2 × M (cid:12) is a result of retainingnominal resolution in the collisionless particles by allow-ing new sink particles to form instead of having singleparticles with anomalously high mass. Another fractionof sink particles starts out with much lower masses andsome particles never accrete enough mass to form a mas-sive star and produce a SN (teal circles).The additional sink gas mass is expected to changeonly the high density portion of each PDF, as it hasalready passed the density threshold of ρ c = 600 cm − .We first attempt to accommodate the sink gas mass byiterating on the power-law slope until enough additionalmass has been added. This method results in the orangedashed curves in Figure 12 in Appendix D, while the or-ange solid curves are combined log-normal and power-law fits made to the new distributions. This methodworks in all but the most dense and star-forming cells(bottom right panel in the same figure), where the sinkgas mass fraction of 93 % results in a nonphysical posi-tive power-law slope. In the latter case we fix α to -1.5as expected for evolved clouds undergoing rapid gravi-tational collapse (Girichidis et al. 2014) while allowingthe width, σ s , to change, moving the center such thatthe new log-normal matches the slope of the originallog-normal on the low-density side.The resulting fits, shown with orange curves in Fig-ure 14, are used as look-up tables for each cell in the Figure 5.
Distribution of sink particle gas mass and age.The teal circles and teal shaded area illustrate the selectioncriteria of category I and II sink particles, respectively, asdescribed in Section 4.3.
Skirt output, by interpolating in (cid:104) n H (cid:105) V and (cid:104) n SFR (cid:105) V .For cell densities below the minimum grid point in the arepo -M51 look-up tables, we adopt a single, narrowlog-normal PDF corresponding to a Mach number of 1and a turbulent forcing parameter b = 1 / n H = 10 − cm − to n H = 10 cm − . Figure 6.
Examples of density PDFs for the same three simba -25 galaxies as shown in Figure 3, here shown with dif-ferent colors.
Examples of the resulting galaxy-wide density PDFscan be seen in Figure 6 for the same galaxies as shownin Figure 3. The original density PDF in the cell data
IGAME v3 Table 2.
Parameters for the one-zone
Cloudy models.Model parameter Min value Max value step sizelog n H [cm − ] -4 7 1log Z [Z (cid:12) ] -2 0.5 0.5log N H [cm − ] 17 22 0.5log F FUV [ G ] -7 4 1log DTM ratio -2 -0.2 0.2 from Skirt (solid line) is compared to the derived PDFadopting a single log-normal with Mach number 10 andforcing parameter b = 1 / FIR line emission with
Cloudy (step 4)
The final step of s´ıgame v3 is to derive the line emis-sion, having gathered all the necessary information bypost-processing the simulation. To this end, we use alibrary of one-zone
Cloudy models that span the pa-rameter ranges listed in Table 2. For each value of the F FUV flux, the cosmic ray ionization rate ξ is scaled as ξ = ( F F UV /G ) ξ , (4)where ξ = 10 − s − is taken as the average Milky Wayvalue (Indriolo et al. 2007) and the Solar neighborhoodFUV flux G = 1 . × − erg cm − s − (Habing 1968).The one-zone Cloudy models must now be sampledaccording to the gas density PDFs found in Section 4.3,such that each combination of (cid:104) n SFR (cid:105) V and (cid:104) n H (cid:105) V inthe galaxy simulations corresponds to a density PDF-weighted sum of the 12 different one-zone Cloudy mod-els along the log n H axis that correspond to that re-gion’s other parameters (log Z , log N H , log F F UV andlog DTM). For this step, six grid points in (cid:104) n H (cid:105) V (from10 − to 10 cm − ) and four grid points in (cid:104) n SFR (cid:105) V (from10 − . to 10 . M (cid:12) yr − ) are used. For each combina-tion of (cid:104) n H (cid:105) V , (cid:104) n SFR (cid:105) V we take the gas density PDFthat comes closest in terms of both values, and shift thecenter of the lognormal to exactly match (cid:104) n H (cid:105) V . Thisshift is performed to ensure that the PDFs generated arenot only centered on the six chosen (cid:104) n H (cid:105) V grid points,but can fill out the density space and generate a smoothtotal PDF as shown in Fig. 6. The sampling of one-zone Cloudy models is carried out for all combinationsof N H , F FUV , Z , and DTM ratio, leading to a look-up table of 190,080 different combinations of one-zone Cloudy models, with one such table for each spectralline considered. 4.4.1.
User-defined sub-grid attenuation function
The
Skirt calculation yields the average interstellarradiation field in cells with sizes ranging from 19 pc to6 . simba -100 galaxies), but due to the res-olution limit of the underlying simba simulation, addi-tional attenuation by parsec-size sub-structures is nottaken into account by Skirt . As a way of compensatingfor this lack of resolution, we have included in the frame-work an optional user-defined function to add extinctionon sub-grid scales. In practice, this is done by generat-ing a one-zone
Cloudy grid as described in Section 4.4but for which the incident spectra were attenuated usingthe “extinguish” command in
Cloudy . This commanddiminishes the incident spectrum with a simple power-law corresponding to an intervening slab of gas of fixedcolumn density, which is set to 10 cm − here. For the simba simulations at hand, we test the framework with avery simple extinction function, that only adds extinc-tion to Cloudy grid cells of density above 10 cm − .When constructing the look-up table by sampling thegas density PDFs (Section 4.3), these extinguished mod-els were sampled for densities n H > cm − ( ξ wasnot modified for these models). This procedure roughlymimics a scenario in which high density gas regions arethe most shielded from ionizing photons. For a differentbase simulation, the attenuation function might need tobe revised or completely omitted in the very high reso-lution case. In the following section we will compare theresults with and without this additional function. VALIDATION OF THE CODEThe measure by which we evaluate the success of agiven sub-grid model is how well it matches observedFIR emission line luminosities of galaxies. In this paper,we restrict the comparison to nearby galaxies observedwith
Herschel or the
Infrared Space Observatory as foundin the literature. 5.1.
Simulation runs
In order to compare different assumptions in the sub-grid process, we perform four different runs on thesame z ∼ simba -100 box as wellas 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 ( in Table 3), we reportthe line luminosities that s´ıgame returns when thedensity fragmentation is carried out the “traditionalway” with a single log-normal density PDF correspond-ing to a Mach number 10 and a forcing parameter of1/3. Next, we applied s´ıgame v2 to the same setof galaxies, using the modifications to the code de-scribed in Leung et al. (2020) together with a new2 Olsen et al.
Table 3.
The different s´ıgame runs compared in this study.Run name Description100Mpc M10 Adopting one log-normal withMach number 10 with a forcingparameter 1/3 for the sub-griddensity profile (cf. Section 4.3).100Mpc SIGAMEv2 Comparison run with s´ıgame v2applied to the simba -100 sam-ple.100Mpc arepoPDF Adopting density PDFs with apowerlaw tail as parametrizedfor a higher resolution arepo -M51 simulations (cf. Sec-tion 4.3). This is our “defaultrun”.25Mpc arepoPDF Same as “arepoPDF” but forthe simba -25 simulation volume(cf. Section 3).100Mpc arepoPDF no ext Same as “arepoPDF” butwithout the additional ex-tinction function described inSection 4.4.1.
Cloudy grid at z = 0 ( ). Thenext simulation run adopts the new gas fragmentationscheme described in this paper, and will be deferredto as our “default run” ( ). Thedefault settings are also applied to a smaller sampleof 240 galaxies in the simba -25 box as a convergencetest ( ). Finally, we also investigatethe resulting line emission when not including the ad-ditional sub-grid attenuation function as described inSection 4.4.1 and otherwise adopted in the default run( ).5.2. Observations used
In the comparison with observations, we include allgalaxies found with total line luminosities and SFR es-timates. This is a rough validation check, since bothobserved and simulated galaxies span a wide range ofgalaxy types, but even so serves to show whether themodel can at all reproduce the observed range in values.Almost all of the FIR emission lines considered herehave been tested as tracers of SFR with varying de-grees of success. As a result, the literature holds avast number of nearby galaxies with both infrared emis-sion line measurements and derived SFRs. In order tobest compare with observations, we estimate the SFRof the model galaxies as a mean over the past 100 Myr,and not the instantaneous SFR of gas particles. Fur- thermore, all SFRs derived using a Salpeter (1955) IMFhave been converted to the equivalent SFR of a Chabrier(2003) IMF as in the simba simulation, by adopting a − .
24 dex correction (e.g. Mitchell et al. 2013). The ob-served SFRs were obtained with various methods sum-marized here; For the sample of mixed type galaxies inKamenetzky et al. (2016) and the LIRGS in D´ıaz-Santoset al. (2013), the far-infrared luminosities, L FIR (40–120 µ m and 42 . . µ m in the two respective cases),of those papers are converted into SFRs according toKennicutt (1998). For the sample of dwarf galaxies ofCormier et al. (2015), SFRs come from R´emy-Ruyeret al. (2015) who derived SFR from FUV and H α lu-minosities, corrected for dust attenuation. The galaxysample of Schruba et al. (2012) contains CO(2-1) lumi-nosities for 16 dwarf galaxies (with one galaxy in com-mon with the Cormier et al. (2015) sample) as well asCO(1-0) luminosities for nearby galaxies compiled fromthe literature, of which we show 22 CO(1-0) observationsfrom the HERA CO-Line Extragalactic Survey (HER-ACLES; 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 and24 µ m emission following the approach in Bigiel et al.(2008) and Leroy et al. (2008). For the SFRs of the 50nearby galaxies we use the values compiled from litera-ture in Schruba et al. (2012), and refer to the referencestherein for specific methods. Finally, for the sample ofmain sequence and starburst galaxies of Brauher et al.(2008), we convert the reported 63 µ m flux to an SFRusing the 70 µ m monochromatic SFR conversion relationof Calzetti et al. (2010).5.3. Line luminosity vs. SFR correlations
Figure 7 shows line luminosities as a function ofSFR for all 8 emission lines listed in Table 1 withdefault settings and the simba -100 galaxy sample (see Table 3). Galaxy symbols arecolor-coded by their molecular gas mass surface den-sity, calculated using the amount of H predicted by Cloudy within the gas half-mass radius from caesar .The panels are organized such that the left column con-tains line emission tracing primarily PDRs and ionizedregions of the ISM, while the right shows tracers of pri-marily neutral and molecular regions in the ISM. In allpanels, a black dashed line indicates a log-linear fit tothe observed line luminosities binned in terms of SFR.For panels with [C II ], [O III ]88, and [O I ]63 line emis-sion, the relations from De Looze et al. (2014) for localstarburst and dwarf galaxies are also shown.For all emission lines except the [N II ] and [O I ]63lines, there is a good overall agreement with observa- IGAME v3 Figure 7.
Correlations between line emission and SFR for the s´ıgame results compared to observations of nearby galaxies:LIRGs (D´ıaz-Santos et al. 2017), dwarf galaxies (Madden et al. 2013; Cormier et al. 2015), main-sequence and starburst galaxies(Brauher et al. 2008), (U)LIRGs (Farrah et al. 2013; Zhao et al. 2016), the randomly selected 0 . < z < .
02 COLD GASSgalaxy sample of Accurso et al. (2017), and the mixed sample of Kamenetzky et al. (2016) and Schruba et al. (2012). For thedwarf galaxy sample, we only include seven galaxies with galaxy-integrated luminosities, following the criteria of Accurso et al.(2017). For the sample of Kamenetzky et al. (2016), we include only galaxies with optical sizes smaller than 47” as a veryconservative measure to insure that the
Herschel /PACS field of view included all line emission. See the text for how SFR wascalculated. Olsen et al. tions. The simulated line luminosities span the range inobserved line luminosities for the same range in SFRs.The largest discrepancy between observations and sim-ulations can be found for the [N II ] and [O I ]63 lines,originating in primarily H II regions and dense, neutralgas, respectively (Table 1). Here, the median deviationof our models above the log-linear fit to all observedsamples, is 1 .
6, 1 . . II ]122, [N II ]205and [O I ]63, respectively. Both [NII] lines are excited bycollisions with electrons, with [N II ]122 having a highercritical density than [N II ]205 (Table 1). The slightlyhigher offset for [N II ]122 compared to [N II ]205, trans-lates into a higher [N II ]122/205 line ratio than observed,indicating that our simulated galaxies have on averagehigher electron densities, since this line ratio is expectedto increase with electron density (e.g. Herrera-Camuset al. 2016). [O I ]63 on the other hand arises in PDRsof much higher densities, and rather than depending ondensity, the excitation of [O I ]63 depends on the tem-perature of the gas as set by the FUV radiation thatheats the gas primarily through the photoelectric effect.In the current s´ıgame framework, the local FUV fluxis derived with Skirt radiative transfer together withan optional additional attenuation as described in Sec-tions B and 4.4.1, respectively. Without the last stepof additional attenuation, we see an even larger mediandeviation of 2 . I ]63 luminosity.In all cases of CO line emission, a large fraction ofour simulated galaxies with SFR (cid:46) (cid:12) /yr displaya reduced CO luminosity, forming a small ”tail” oflow luminosity galaxies reaching down below the ob-served relation by ∼ (cid:46) (cid:12) pc − ) andlow metallicities ( (cid:104) Z (cid:105) SFR (cid:46) . (cid:12) ). In the panelwith L CO(1 − , we see that this is in agreement withthe mixed and dwarf galaxy sample of Schruba et al.(2012) in CO(2-1) and CO(2-1), respectively (stars).This agreement is presumably due to the overall lowmetallicity of the dwarf sample, extending down to ∼ .
03 Z (cid:12) (M81 Dw A; Schruba et al. 2012), althoughour galaxy sample only goes down to 0.1 Z (cid:12) in SFR-weighted metallicity. While the observed galaxies wouldseem to agree with our simulations in terms of COluminosity, it is worth pointing out that these galax-ies have much lower [C II ] luminosities than our simu-lated galaxies as shown in Figure 8 which shows CO(1-0)against [C II ] luminosity for the , and runs. The galaxy samples of Kamenetzky et al. (2016), Cormier et al. (2015), and Accurso et al. (2017) are alsoshown with the latter two sample galaxies color-codedby metallicity for comparison with the simulated galax-ies. To guide the eye, a log-linear fit is made to the ob-servations and shown in dashed black line. As dictatedby the brown contour lines, the omission of additionalattenuation at high densities results in larger deviationsfrom the observed [C II ]-CO relation. Figure 8. [C II ] versus CO(1-0) luminosity for the s´ıgame v3 results ( run; green contours, run; brown dashed contoursand run; blue dotted contours) com-pared to the observed galaxy samples of Kamenetzky et al.(2016), Cormier et al. (2015) and Accurso et al. (2017). Thesimulated galaxies and (where possible) the observed galax-ies have been color-coded by metallicity. The dashed lineindicates a log-linear fit to the observations alone. Comparison of different s´ıgame runs
Figure 9 condenses the data in Figure 7 into mediandeviations between the simulated galaxy line luminosi-ties and a linear fit to the observed ones (black dashedlines in Fig. 7) as a function of SFR. The standard de-viation of the offsets for the observed galaxies is indi-cated with a dark grey envelope. The different s´ıgame runs listed in Table 3 can now be compared for each ofthe eight emission lines considered here, with the fine-structure lines sorted in order of increasing critical den-sity (Table 1).For run , the resulting deviations fromobservations are shown with orange bars. Surprisingly,this simple assumption does a very good job for the ion-ized species of [N II ]122, [N II ]205, [O III ]88 and [C II ],but completely underestimates CO emission and to alesser degree [O I ]63 as well. The [O I ]63, CO(1-0), (2-1) and (3-2) lines fall below the observed line-SFR re-lation by median deviations of 0 .
6, 3 .
1, 4 . . IGAME v3 Figure 9.
Box plot performance comparison of the different s´ıgame runs listed in Table 3 for all lines considered here, with thefine-structure lines sorted in order of increasing critical density. Light and dark grey shaded areas represent 1 and 3 σ spread inthe observations. The boxes extend from the lower to the upper quartile values, whiskers to a factor 1 . ρ c rit of atomic hydrogen for [C II ]. respectively. The underestimation of the molecular linesis not surprising, given that the density PDF does notextend significantly beyond densities > cm − wheremolecules form (see Fig. 6).With run (purple bars), we testhow well the previous s´ıgame framework, which has sofar only been applied at z ∼ simba -100 z = 0 sample. While this version of the codegives good results for [C II ] and [O I ]63, and comes closerto the observed CO luminosities than the run, it significantly underestimates the [N II ] lines. Forall lines except [O I ]63, this version also produces out-liers with extremely low line luminosities, as shown withopen circles (and indeed, the plot does not extend toshow all outliers).The test run without additional extinction, (brown bars), displaysa better agreement with the observed CO and [O III ]88luminosities than the run, but our PDFframework fares worse for all other lines. In particularthe [O I ]63 line predictions fall above the observed lu-minosities by a median deviation of 2.0 dex, and thereare still some very low CO luminosity predictions.The default run of the current version of s´ıgame , (green bars) comes closest to allline luminosities simultaneously, with the least numberof outliers. The [N II ]122 and [O I ]63 line predictionsshow the poorest match with observations, remainingoverestimated by median deviations of 1 . . I ]63 line over-predictions are reducedby about 1 dex from 2.0 dex to 1.2 dex above the obser-vational [O I ]63–SFR relation, and the tail of low COluminosities is reduced in length by several dex. Thecontinued overestimation of the [O I ]63 and [N II ] lineemission, however, suggests that a more careful treat-ment of the sub-grid attenuation and/or other featuresremains missing. A proper treatment of the flux dis-tribution on these scales would require considerationsabout the properties of natal clouds around young starsas demonstrated with the one-dimensional stellar feed-back model warpfield (Rahner et al. 2017). [C II ] and[O III ]88 remain slightly overestimated, both of whichoriginate mainly in the atomic ISM phase in our model,hinting at an overestimation of atomic gas mass vs ion-ized gas mass in our gas fragmentation scheme. We donot consider the overestimation of CO(3-2) a problemdue to the lack of observations compiled here and thepossibility that the relation with SFR is not linear asthe large deviations suggest.Finally, the test run with default settings at highermass resolution, (blue bars), showsdeviations consistent with the simba -100 sample, indi-cating that the s´ıgame results converge well with in-creasing mass resolution. In Fig. 15 in Appendix E weinvestigate in more depth how well the two simulationruns agree with one another for a range of physicalgalaxy properties, finding that for similar galaxies, theluminosities are within 1 dex of each other with slightlyhigher deviations towards low values of M (cid:63) , M gas , SFR6 Olsen et al. and (cid:104) Z (cid:105) SFR where the gas particle number in simba -100drops below 500.5.5.
Moment 0 maps of line emission
Fig 10 shows moment0 maps of CO(1-0), [C II ] and[O III ]88 emission constructed for the same galaxies from simba -100 shown in Figure 3. The output datacubefrom s´ıgame contains spatial information, cell sizes andcell luminosities that were combined to derive the sur-face brightness of each pixel in the final image, weighingeach cell by its volume filling factor within the columncovered by each pixel . We note that the calculation ofthe moment0 maps assumes an ISM fully transparentto the FIR line emission, an assumption that might nothold for shorter wavelengths. Some expected differencesin how the emission lines trace the ISM can be observed,such as the preference for [C II ] to arrive from PDRs nearstar forming regions compared to the broader CO(1-0)emission, and the relatively concentrated [O III ]88 emis-sion restricted to H II regions experiencing harder radi-ation fields from young stars. DISCUSSION6.1.
Comparison with previous models
We can compare our [C II ]-SFR relation withthat found through similar techniques and hence bet-ter understand the decisive differences between themodels. Figure 11 shows the [C II ]-SFR relationsfound by s´ıgame for the simba -100 galaxies with the and run (green contour lines) to-gether with the relations derived by Popping et al.(2019) and Ramos Padilla et al. (2020). From the run, we also show the distribu-tion of only MS galaxies, chosen to be within the 16 and84 percentile of the Salim et al. (2018) relation (purpledashed contour lines). The distribution of s´ıgame simu-lated galaxies are generally too high compared with ob-servations and exhibits a much shallower slope towardslower SFRs, where our model seems to overpredict the[C II ] luminosities by about 1 dex. When selecting onlyMS galaxies, the agreement with observations improvesonly at the high SFR-end of (cid:38) (cid:12) yr − .6.2. Choice of simulation for the fragmentation task
We have only shown results for the present frameworkusing one type of simulation, namely the arepo suite ofM51 realizations. We also examined an arepo simula-tion of the central molecular zone (Sormani et al. 2020),but found no significant change from the PDFs extracted https://github.com/jaymotka/moment0 maps and shown in Figure 14. However, choosing a completelydifferent galaxy simulation technique (e.g. models runwith adaptive mesh refinement codes such as ramses , enzo , or athena Teyssier 2002; Bryan et al. 2014; Kim& Ostriker 2017) might return different PDFs, the effectof which we have not studied. Instead, we have focusedon making the current version of s´ıgame flexible andpublicly available so that the user can include any highresolution simulation for the gas fragmentation if theyso desire.6.3.
Small-scale radiation field structure
As discussed in the previous section, this version of s´ıgame fragments the gas density on scales smaller thanthe
Skirt cell sizes (typically larger than ∼
15 pc),thereby allowing for clumps of higher densities on scalesnot resolved by the parent simulation. In a similar man-ner, we could also expect the radiation field to be frag-mented such that the mass-weighted FUV flux distri-bution on parsec scales can differ substantially from theoverall cell mean flux value and allow for parsec-size fea-tures such as H II regions or shielded, dense, molecularcores. Judging from the overall agreement between sim-ulated and observed lines, we speculate that most of theFIR lines considered here arise from structures of largerextent where parsec resolution in radiation is not nec-essary, with the exception of [O I ]63, which is heavilyoverestimated.6.3.1. Missing treatment of active galactic nuclei
While the co-evolution of active galactic nuclei (AGN)and their host galaxies is simulated in simba , the cur-rent version of s´ıgame does not include any effects ofAGN presence, such as additional heating and radiation.The additional X-ray heating has been shown to increaseexcitation of high-J CO lines at high redshifts (Valliniet al. 2019).6.3.2.
Missing shock-heating of the gas
In starburst galaxies and mergers, shocks are ex-pected to act as an additional heating source in theISM. Although the next version of
Cloudy may in-clude a treatment for shock-heated gas, the current ver-sion of
Cloudy iterates to find a thermal and ionizationequilibrium, thereby ignoring any shocked state of thegas that might be present in the simulation. One wayto treat shock-heated gas, would be to set up a sep-arate grid of models, following the technique used for
MAPPINGS III (Allen et al. 2008). Turbulence (2—10km s − in velocity dispersion) and/or density-dependentmagnetic fields are also believed to play a role in settingthe [O I ]63 emission, through their effect on line width,shielding and pumping, chemical and thermal state of IGAME v3 Figure 10.
Moment0 maps of the same galaxies shown in Fig 3, for three of the emission lines investigated here. the gas (Canning et al. 2016). However, testing the ef-fect of these would require
Cloudy models of more thanone zone where optical depth is taken into account. CONCLUSIONSThis paper introduces an improved algorithm for esti-mating FIR line emission from large-volume cosmologi-cal galaxy simulations using coarse numerical resolution.We post-process such simulations using a sub-grid modelthat estimates the distribution of dense gas up to densi-ties of ∼ cm − . This fragmentation scheme results ingas at lower densities being compressed to higher den-sities, and hence affects emission lines tracing all ISMphases. We test the scheme on eight low-to-mediumdensity ISM tracers.The density distribution on sub-grid scales is mod-eled by sorting and parametrizing resolved regions in higher-resolution simulations and interpolating on theparametrized functions to set the density distributionfor the cosmological simulations. This statistical ap-proach avoids any assumption about the size and shapeof molecular clouds, the turbulence within those clouds,or the existence of pressure equilibrium between thedense and more diffuse ISM phases. As a demonstra-tion of this scheme, we use data with sub-parsec reso-lution from a model of M51 with arepo (Tress et al.2020). Density PDFs are sampled in 200 pc regions andbinned in terms of the volume-averaged hydrogen andSFR density of each region. The resulting mean PDFsare parametrized and used to generate PDFs by inter-polation for each resolution element of the cosmologi-cal simulation to generate gas densities ∼ –10 cm − otherwise not present in the cosmological simulation re-8 Olsen et al.
Figure 11.
Comparison of s´ıgame v3 [C II ]-SFR rela-tion ( ; green, MS-only galaxies: pur-ple dashed contour lines) with other recent modeling efforts:a multi-phase model of the ISM used as a post-processingstep for the EAGLE cosmological hydrodynamical simula-tions (cyan, orange, and dark blue contour lines; RamosPadilla et al. 2020), and the z = 0 [C II ]-SFR relation foundby Popping et al. (2019) using SAMs (black line), with a 1- σ shaded region. Observations of local dwarf and starburstgalaxies by De Looze et al. (2014) are shown as grey symbolsfor comparison with the models. sults. As a future extension, this approach could also beused to develop a new SFR prescription for cosmologicalsimulations.The local radiation field strength is determined withthe radiative transfer code Skirt , which calculates thelocal dust-attenuated stellar spectrum. The
SWIFT-simIO package is used to map gas properties from thecosmological simulation output to the cell structure from
Skirt . Additional attenuation by gas is implementedthrough a set of
Cloudy models, by matching the trans-mitted spectrum in the NIR-to-FUV of a certain col-umn density of gas to the dust-attenuated spectrum of
Skirt for each gas cell. Line emission and chemical in-formation are derived from an extensive grid of 95,040
Cloudy n H , (cid:104) n SFR (cid:105) V , N H , F FUV , Z and DTM ratio. Finally, to compensate for lack ofresolution in the cosmological simulations used here, weadd attenuation on sub-grid scales with a simple func-tion that can be modified by the user to work on othersimulation types. In this case, we add attenuation for Cloudy grid models with n H > cm − .As validation of the method, results are compared toobserved relations between line luminosity and SFR fora diverse sample of nearby galaxies and eight different emission lines. The method is also checked against thedefault option of adopting a log-normal density profiledrawn from Mach 10 isothermal turbulence for the sub-grid density profile, although we expect this method tofare poorly in the case of cosmological simulations. Thevalidation leaves us with the following conclusions: • The standard method in the literature of adoptingthe density PDF of Mach 10 isothermal turbulencefor fragmenting gas on sub-grid scales results inFIR fine-structure line emission that agrees wellwith observations for lines originating mainly inthe ionized ISM, but drastically underestimateslines from the neutral and molecular regions. Forexample, the CO(1-0), (2-1) and (3-2) rotationallines fall below the observed line-SFR relation bymedian deviations of 3 .
1, 4 . . • The novel method presented here of using a highresolution simulation and a new analytic PDF ap-proach for the gas fragmentation results in emis-sion lines that match observations well regardlessof their origin in the ISM, with the noticeable ex-ception of [N II ]122, [N II ]205 and [O I ]63 that areoverestimated by on average 1 .
6, 1 .
2, and 1 . • We find that the resulting line emission of [O I ]63in particular is highly dependent on the user-defined attenuation function, without which theoverestimation of the [O I ]63 luminosities is about1 dex higher still. This underlines the need fora more careful treatment of the radiation on sub-grid scales ( <
200 pc), where denser regions shouldproduce additional attenuation of the interstellarradiation field than what is included in this frame-work by default.We have presented s´ıgame v3 which is a flexibleframework that can be adapted to any cosmologicalor galaxy simulation, and now with the option to fur-ther fragment the gas on sub-grid scales using a high-resolution simulation of choice. This tool may be use-ful for the interpretation of current data from e.g.ALMA, NOEMA and SOFIA as well as from futurespace and balloon missions such as JWST, GUSTO andASTHROS.
IGAME v3
Skirt , Robert Thompson for developing cae-sar , and the yt team for development and support of yt (Turk et al. 2011). The authors are also grateful tothe entire simba collaboration that helped in the anal-ysis and understanding of the simba output. We ac-knowledge use of the python SWIFTsimIO (Borrow &Borrisov 2020). This research has made use of theNASA/IPAC Extragalactic Database (NED), which isfunded by the National Aeronautics and Space Admin-istration and operated by the California Institute ofTechnology. KPO is funded by NASA under award No80NSSC19K1651. BB is partly funded by the PackardFellowship for Science and Engineering. JB is supportedby STFC studentship ST/R504725/1. RJS acknowl-edges an STFC ERF fellowship (grant ST/N00485X/1)and HPC from the DiRAC facility (ST/P002293/1). M-MML was partly funded by NSF grant AST18-15461.RT was supported by DFG SFB 881 “The Milky WaySystem” and the Excellence Cluster STRUCTURES(EXC 2181-390900948), as well as ERC Synergy GrantECOGAL (project ID 855130). The Cosmic Dawn Cen-ter of Excellence is funded by the Danish National Re-search Foundation under grant No. 140.
Software: python arepo (Springel 2010),
Cloudy version 17.2 (Ferland et al. 2017),
Skirt ver-sion 9.0 (Camps & Baes 2020),
SWIFTsimIO (https://github.com/SWIFTSIM/swiftsimio; Borrow & Bor-risov 2020))0
Olsen et al.
REFERENCES
Accurso, G., Saintonge, A., Catinella, B., et al. 2017,MNRAS, 470, 4750, doi: 10.1093/mnras/stx1556Allen, M. G., Groves, B. A., Dopita, M. A., Sutherland,R. S., & Kewley, L. J. 2008, ApJS, 178, 20,doi: 10.1086/589652Alves, F. O., Cleeves, L. I., Girart, J. M., et al. 2020, ApJL,904, L6, doi: 10.3847/2041-8213/abc550Alves, J., Lombardi, M., & Lada, C. J. 2017, A&A, 606, L2,doi: 10.1051/0004-6361/201731436Appleby, S., Dav´e, R., Kraljic, K., Angl´es-Alc´azar, D., &Narayanan, D. 2020, MNRAS, 494, 6053,doi: 10.1093/mnras/staa1169Astropy Collaboration, Robitaille, T. P., Tollerud, E. J.,et al. 2013, A&A, 558, A33,doi: 10.1051/0004-6361/201322068Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846,doi: 10.1088/0004-6256/136/6/2846Borrow, J., & Borrisov, A. 2020, The Journal of OpenSource Software, 5, 2430, doi: 10.21105/joss.02430Brauher, J. R., Dale, D. A., & Helou, G. 2008, ApJS, 178,280, doi: 10.1086/590249Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014,ApJS, 211, 19, doi: 10.1088/0067-0049/211/2/19Burkhart, B. 2018, ApJ, 863, 118,doi: 10.3847/1538-4357/aad002Burkhart, B., Collins, D. C., & Lazarian, A. 2015, ApJ,808, 48, doi: 10.1088/0004-637X/808/1/48Burkhart, B., & Mocz, P. 2019, ApJ, 879, 129,doi: 10.3847/1538-4357/ab25edBurkhart, B., Stalpes, K., & Collins, D. C. 2017, ApJL,834, L1, doi: 10.3847/2041-8213/834/1/L1Calzetti, D., Wu, S. Y., Hong, S., et al. 2010, ApJ, 714,1256, doi: 10.1088/0004-637X/714/2/1256Camps, P., & Baes, M. 2020, Astronomy and Computing,31, 100381, doi: 10.1016/j.ascom.2020.100381Canning, R. E. A., Ferland, G. J., Fabian, A. C., et al.2016, MNRAS, 455, 3042, doi: 10.1093/mnras/stv2390Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392Chen, H. H.-H., Burkhart, B., Goodman, A., & Collins,D. C. 2018, The Astrophysical Journal, 859, 162,doi: 10.3847/1538-4357/aabaf6Chevance, M., Kruijssen, J. M. D., Vazquez-Semadeni, E.,et al. 2020, SSRv, 216, 50,doi: 10.1007/s11214-020-00674-xCollins, D. C., Kritsuk, A. G., Padoan, P., et al. 2012, ApJ,750, 13, doi: 10.1088/0004-637X/750/1/13Colombo, D., Hughes, A., Schinnerer, E., et al. 2014, ApJ,784, 3, doi: 10.1088/0004-637X/784/1/3 Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833,doi: 10.1088/0004-637X/712/2/833Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486,doi: 10.1088/0004-637X/699/1/486Cormier, D., Madden, S. C., Lebouteiller, V., et al. 2015,A&A, 578, A53, doi: 10.1051/0004-6361/201425207Cowie, L. L., & Songaila, A. 1986, ARA&A, 24, 499,doi: 10.1146/annurev.aa.24.090186.002435Dav´e, R., Angl´es-Alc´azar, D., Narayanan, D., et al. 2019,MNRAS, 486, 2827, doi: 10.1093/mnras/stz937Dav´e, R., Thompson, R., & Hopkins, P. F. 2016, MNRAS,462, 3265, doi: 10.1093/mnras/stw1862De Looze, I., Cormier, D., Lebouteiller, V., et al. 2014,A&A, 568, A62, doi: 10.1051/0004-6361/201322489D´ıaz-Santos, T., Armus, L., Charmandaris, V., et al. 2013,ApJ, 774, 68, doi: 10.1088/0004-637X/774/1/68—. 2017, ApJ, 846, 32, doi: 10.3847/1538-4357/aa81d7Farrah, D., Lebouteiller, V., Spoon, H. W. W., et al. 2013,ApJ, 776, 38, doi: 10.1088/0004-637X/776/1/38Federrath, C., Klessen, R. S., & Schmidt, W. 2008, ApJL,688, L79, doi: 10.1086/595280Ferland, G. J., Chatzikos, M., Guzm´an, F., et al. 2017,RMxAA, 53, 385. https://arxiv.org/abs/1705.10877Foreman-Mackey, D., Sick, J., & Johnson, B. 2014,python-fsps: Python bindings to FSPS (v0.1.1), v0.1.1,Zenodo, doi: 10.5281/zenodo.12157Girichidis, P., Konstandin, L., Whitworth, A. P., &Klessen, R. S. 2014, ApJ, 781, 91,doi: 10.1088/0004-637X/781/2/91Gnedin, N. Y., Kravtsov, A. V., & Chen, H.-W. 2008, ApJ,672, 765, doi: 10.1086/524007Gnedin, N. Y., Tassis, K., & Kravtsov, A. V. 2009, ApJ,697, 55, doi: 10.1088/0004-637X/697/1/55Goldsmith, P. F., Langer, W. D., Pineda, J. L., &Velusamy, T. 2012, ApJS, 203, 13,doi: 10.1088/0067-0049/203/1/13Habing, H. J. 1968, BAN, 19, 421Harris, C. R., Millman, J., van der Walt, S. J., et al. 2020,Nature, doi: 10.1038/s41586-020-2649-2Herrera-Camus, R., Bolatto, A., Smith, J. D., et al. 2016,ApJ, 826, 175, doi: 10.3847/0004-637X/826/2/175Hewitt, J. W., & Yusef-Zadeh, F. 2009, ApJL, 694, L16,doi: 10.1088/0004-637X/694/1/L16Hopkins, P. F. 2015, MNRAS, 450, 53,doi: 10.1093/mnras/stv195—. 2017, arXiv e-prints, arXiv:1712.01294.https://arxiv.org/abs/1712.01294Hunter, J. D. 2007, Computing in Science & Engineering, 9,90, doi: 10.1109/MCSE.2007.55
IGAME v3 Indriolo, N., Geballe, T. R., Oka, T., & McCall, B. J. 2007,ApJ, 671, 1736, doi: 10.1086/523036Inoue, S., Yoshida, N., & Yajima, H. 2020, arXiv e-prints,arXiv:2008.12484. https://arxiv.org/abs/2008.12484Jones, A. P., K¨ohler, M., Ysard, N., Bocchio, M., &Verstraete, L. 2017, A&A, 602, A46,doi: 10.1051/0004-6361/201630225Kainulainen, J., Beuther, H., Henning, T., & Plume, R.2009, A&A, 508, L35, doi: 10.1051/0004-6361/200913605Kamenetzky, J., Rangwala, N., Glenn, J., Maloney, P. R.,& Conley, A. 2016, ApJ, 829, 93,doi: 10.3847/0004-637X/829/2/93Katsianis, A., Xu, H., Yang, X., et al. 2020, MNRAS,doi: 10.1093/mnras/staa3236Kennicutt, Robert C., J. 1998, ARA&A, 36, 189,doi: 10.1146/annurev.astro.36.1.189Kessler, M. F., Steinz, J. A., Anderegg, M. E., et al. 1996,A&A, 500, 493Kim, C.-G., & Ostriker, E. C. 2017, ApJ, 846, 133,doi: 10.3847/1538-4357/aa8599Klessen, R. S. 2000, ApJ, 535, 869, doi: 10.1086/308854Kroupa, P. 2001, MNRAS, 322, 231,doi: 10.1046/j.1365-8711.2001.04022.xKrumholz, M. R., & Gnedin, N. Y. 2011, ApJ, 729, 36,doi: 10.1088/0004-637X/729/1/36Leitherer, C., Ekstr¨om, S., Meynet, G., et al. 2014, TheAstrophysical Journal Supplement Series, 212, 14,doi: 10.1088/0067-0049/212/1/14Leroy, A. K., Walter, F., Brinks, E., et al. 2008, AJ, 136,2782, doi: 10.1088/0004-6256/136/6/2782Leroy, A. K., Walter, F., Bigiel, F., et al. 2009, AJ, 137,4670, doi: 10.1088/0004-6256/137/6/4670Leroy, A. K., Usero, A., Schruba, A., et al. 2017, ApJ, 835,217, doi: 10.3847/1538-4357/835/2/217Leung, T. K. D., Olsen, K. P., Somerville, R. S., et al. 2020,arXiv e-prints, arXiv:2004.11912.https://arxiv.org/abs/2004.11912Li, Q., Narayanan, D., & Dav´e, R. 2019, MNRAS, 490,1425, doi: 10.1093/mnras/stz2684Li, Q., Narayanan, D., Dav`e, R., & Krumholz, M. R. 2018,ApJ, 869, 73, doi: 10.3847/1538-4357/aaec77Lombardi, M., Alves, J., & Lada, C. J. 2015, A&A, 576, L1,doi: 10.1051/0004-6361/201525650Lord, S. D., Hollenbach, D. J., Haas, M. R., et al. 1996,ApJ, 465, 703, doi: 10.1086/177455Lower, S., Narayanan, D., Leja, J., et al. 2020, ApJ, 904,33, doi: 10.3847/1538-4357/abbfa7Lupi, A., & Bovino, S. 2020, MNRAS, 492, 2818,doi: 10.1093/mnras/staa048 Madden, S. C., R´emy-Ruyer, A., Galametz, M., et al. 2013,PASP, 125, 600, doi: 10.1086/671138Maeder, A. 2009, Physics, Formation and Evolution ofRotating Stars (Berlin Heidelberg: Springer),doi: 10.1007/978-3-540-76949-1Mitchell, P. D., Lacey, C. G., Baugh, C. M., & Cole, S.2013, MNRAS, 435, 87, doi: 10.1093/mnras/stt1280Mocz, P., Burkhart, B., Hernquist, L., McKee, C. F., &Springel, V. 2017, ApJ, 838, 40,doi: 10.3847/1538-4357/aa6475Narayanan, D., & Krumholz, M. R. 2014, MNRAS, 442,1411, doi: 10.1093/mnras/stu834Narayanan, D., Li, Y., Cox, T. J., et al. 2008, ApJS, 174,13, doi: 10.1086/521776Ocvirk, P., Gillet, N., Shapiro, P. R., et al. 2016, MNRAS,463, 1462, doi: 10.1093/mnras/stw2036Ocvirk, P., Aubert, D., Sorce, J. G., et al. 2020, MNRAS,496, 4087, doi: 10.1093/mnras/staa1266Olsen, K., Greve, T. R., Narayanan, D., et al. 2017, ApJ,846, 105, doi: 10.3847/1538-4357/aa86b4Olsen, K., Pallottini, A., Wofford, A., et al. 2018, Galaxies,6, 100, doi: 10.3390/galaxies6040100Olsen, K. P., Greve, T. R., Brinch, C., et al. 2016, MNRAS,457, 3306, doi: 10.1093/mnras/stw162Olsen, K. P., Greve, T. R., Narayanan, D., et al. 2015, ApJ,814, 76, doi: 10.1088/0004-637X/814/1/76Padoan, P., Haugbølle, T., Nordlund, ˚A., & Frimann, S.2017, ApJ, 840, 48, doi: 10.3847/1538-4357/aa6afaPadoan, P., & Nordlund, ˚A. 2011, ApJ, 730, 40,doi: 10.1088/0004-637X/730/1/40Pallottini, A., Ferrara, A., Decataldo, D., et al. 2019,MNRAS, 487, 1689, doi: 10.1093/mnras/stz1383Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010,A&A, 518, L1, doi: 10.1051/0004-6361/201014759Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS,473, 4077, doi: 10.1093/mnras/stx2656Planck Collaboration, Ade, P. A. R., Aghanim, N., et al.2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830Popping, G., Narayanan, D., Somerville, R. S., Faisst,A. L., & Krumholz, M. R. 2019, MNRAS, 482, 4906,doi: 10.1093/mnras/sty2969Price-Whelan, A. M., Sip˝ocz, B. M., G¨unther, H. M., et al.2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4fRahner, D., Pellegrini, E. W., Glover, S. C. O., & Klessen,R. S. 2017, MNRAS, 470, 4453,doi: 10.1093/mnras/stx1532Ramos Padilla, A. F., Wang, L., Ploeckinger, S., van derTak, F. F. S., & Trager, S. 2020, arXiv e-prints,arXiv:2011.13441. https://arxiv.org/abs/2011.13441 Olsen et al.
R´emy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2015,A&A, 582, A121, doi: 10.1051/0004-6361/201526067Rodr´ıguez Montero, F., Dav´e, R., Wild, V., Angl´es-Alc´azar,D., & Narayanan, D. 2019, MNRAS, 490, 2139,doi: 10.1093/mnras/stz2580Salim, S., Boquien, M., & Lee, J. C. 2018, ApJ, 859, 11,doi: 10.3847/1538-4357/aabf3cSalpeter, E. E. 1955, ApJ, 121, 161, doi: 10.1086/145971Savage, B. D., & Sembach, K. R. 1996, ARA&A, 34, 279,doi: 10.1146/annurev.astro.34.1.279Schneider, N., Bontemps, S., Girichidis, P., et al. 2015,MNRAS, 453, L41, doi: 10.1093/mnrasl/slv101Sch¨oier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., &Black, J. H. 2005, A&A, 432, 369,doi: 10.1051/0004-6361:20041729Schruba, A., Leroy, A. K., Walter, F., et al. 2012, AJ, 143,138, doi: 10.1088/0004-6256/143/6/138Sormani, M. C., Tress, R. G., Glover, S. C. O., et al. 2020,MNRAS, 497, 5024, doi: 10.1093/mnras/staa1999Speagle, J. S., Steinhardt, C. L., Capak, P. L., &Silverman, J. D. 2014, ApJS, 214, 15,doi: 10.1088/0067-0049/214/2/15Springel, V. 2010, MNRAS, 401, 791,doi: 10.1111/j.1365-2966.2009.15715.xStacey, G. J., Geis, N., Genzel, R., et al. 1991, ApJ, 373,423, doi: 10.1086/170062Sun, J., Leroy, A. K., Schruba, A., et al. 2018, ApJ, 860,172, doi: 10.3847/1538-4357/aac326Teyssier, R. 2002, A&A, 385, 337,doi: 10.1051/0004-6361:20011817Thomas, N., Dav´e, R., Angl´es-Alc´azar, D., & Jarvis, M.2019, MNRAS, 487, 5764, doi: 10.1093/mnras/stz1703 Thompson, R. 2015, SPHGR: Smoothed-ParticleHydrodynamics Galaxy Reduction.http://ascl.net/1502.012Tress, R. G., Smith, R. J., Sormani, M. C., et al. 2020,MNRAS, 492, 2973, doi: 10.1093/mnras/stz3600Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, TheAstrophysical Journal Supplement Series, 192, 9,doi: 10.1088/0067-0049/192/1/9Vallini, L., Tielens, A. G. G. M., Pallottini, A., et al. 2019,MNRAS, 490, 4502, doi: 10.1093/mnras/stz2837Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020,Nature Methods, 17, 261,doi: https://doi.org/10.1038/s41592-019-0686-2Watson, D. M., Genzel, R., Townes, C. H., Werner, M. W.,& Storey, J. W. V. 1984, ApJL, 279, L1,doi: 10.1086/184242Weinberger, R., Springel, V., & Pakmor, R. 2020, ApJS,248, 32, doi: 10.3847/1538-4365/ab908cWeinberger, R., Springel, V., Hernquist, L., et al. 2017,MNRAS, 465, 3291, doi: 10.1093/mnras/stw2944Wes McKinney. 2010, in Proceedings of the 9th Python inScience Conference, ed. St´efan van der Walt & JarrodMillman, 56 – 61, doi: 10.25080/Majora-92bf1922-00aWu, X., Kannan, R., Marinacci, F., Vogelsberger, M., &Hernquist, L. 2019, MNRAS, 488, 419,doi: 10.1093/mnras/stz1726Yungelson, L. R., van den Heuvel, E. P. J., Vink, J. S.,Portegies Zwart, S. F., & de Koter, A. 2008, A&A, 477,223, doi: 10.1051/0004-6361:20078345Zhao, Y., Lu, N., Xu, C. K., et al. 2016, ApJ, 819, 69,doi: 10.3847/0004-637X/819/1/69
IGAME v3 A. PROPERTIES OF simba -100 MODEL GALAXY SAMPLETable 4 lists the physical properties of the simba -100 galaxy sample while Table 5 lists their line luminosities, includinga few ( L [OI]145 , L [CI]610 and L [CI]371 ) not used in the validation of the method. Table 4.
Physical properties of the sample of 400 simba -100 galaxiesused in this work. A full version of this table is available online.Galaxy M (cid:63) M g as SFR (cid:104) Z (cid:105) SFR (cid:104) Z (cid:105) mw index [10 M (cid:12) ] [10 M (cid:12) ] [M (cid:12) yr − ] [Z (cid:12) ] [Z (cid:12) ]0 0.0415 0.6653 0.35 0.16 0.021 0.0459 0.7686 0.37 0.10 0.032 0.0492 0.6191 0.18 0.12 0.053 0.0500 0.9535 0.18 0.13 0.024 0.0532 0.8681 0.38 0.16 0.045 0.0534 0.6343 0.35 0.14 0.056 0.0541 1.0024 0.54 0.15 0.047 0.0556 0.7015 0.17 0.15 0.038 0.0565 0.6324 0.52 0.11 0.029 0.0595 0.5862 0.69 0.15 0.0610 0.0618 0.9894 0.55 0.17 0.0411 0.0621 0.6317 0.55 0.14 0.0712 0.0643 0.7573 0.53 0.14 0.0313 0.0650 0.7082 0.18 0.14 0.0614 0.0660 1.0680 0.36 0.19 0.0515 0.0698 0.9691 0.54 0.16 0.0516 0.0703 0.6174 0.37 0.16 0.0517 0.0708 0.9082 0.35 0.18 0.0418 0.0716 0.9648 0.17 0.20 0.0419 0.0721 0.6205 0.72 0.10 0.0820 0.0727 1.1279 0.89 0.19 0.05 Olsen et al. T a b l e . L i n e l u m i n o s i t i e s f o r t h e s a m p l e o f s i m b a - l a x i e s u s e d i n t h i s w o r k . A f u ll v e r s i o n o f t h i s t a b l e i s a v a il a b l e o n li n e . G a l a xy L [ C II ] L [ C I ] L [ C I ] L [ N II ] L [ N II ] L [ O I ] L [ O I ] L [ O III ] L C O ( − ) L C O ( − ) L C O ( − ) i nd e x [ L (cid:12) ][ L (cid:12) ][ L (cid:12) ][ L (cid:12) ][ L (cid:12) ][ L (cid:12) ][ L (cid:12) ][ L (cid:12) ][ L (cid:12) ][ L (cid:12) ][ L (cid:12) ] . e + . e + . e + . e + . e + . e + . e + . e + . e + . e + . e +
02 17 . e + . e + . e + . e + . e + . e + . e + . e + . e - . e + . e +
01 27 . e + . e + . e + . e + . e + . e + . e + . e + . e + . e + . e +
03 34 . e + . e + . e + . e + . e + . e + . e + . e + . e + . e + . e +
03 41 . e + . e + . e + . e + . e + . e + . e + . e + . e - . e + . e +
01 51 . e + . e + . e + . e + . e + . e + . e + . e + . e - . e + . e +
01 61 . e + . e + . e + . e + . e + . e + . e + . e + . e + . e + . e +
02 73 . e + . e + . e + . e + . e + . e + . e + . e + . e + . e + . e +
03 85 . e + . e + . e + . e + . e + . e + . e + . e + . e - . e + . e +
01 91 . e + . e + . e + . e + . e + . e + . e + . e + . e - . e + . e +
01 109 . e + . e + . e + . e + . e + . e + . e + . e + . e + . e + . e +
02 111 . e + . e + . e + . e + . e + . e + . e + . e + . e - . e - . e -
01 126 . e + . e + . e + . e + . e + . e + . e + . e + . e - . e + . e +
01 131 . e + . e + . e + . e + . e + . e + . e + . e + . e - . e + . e +
02 141 . e + . e + . e + . e + . e + . e + . e + . e + . e + . e + . e +
02 159 . e + . e + . e + . e + . e + . e + . e + . e + . e + . e + . e +
02 166 . e + . e + . e + . e + . e + . e + . e + . e + . e - . e + . e +
02 179 . e + . e + . e + . e + . e + . e + . e + . e + . e + . e + . e +
02 188 . e + . e + . e + . e + . e + . e + . e + . e + . e + . e + . e +
03 192 . e + . e + . e + . e + . e + . e + . e + . e + . e - . e - . e - IGAME v3 B. CONVERGENCE TESTS WITH
SkirtSkirt uses a fixed set of photon packets when iterating for a solution to the radiative transfer problem. We tested
Skirt with different photon packet numbers for one large simba -100 galaxy in particular, namely the one shown inthe top row of Figure 3. The resulting total FUV luminosity and FUV flux distribution can be seen in Figure 12 forphoton packet sizes of 10 , 10 , 10 and 10 . There is only a negligible change in FUV luminosity of 0 .
24% comparedto using the lowest photon packet number, but when looking at the flux distribution it becomes clear that at least 10 packets are necessary for a stable result, whereas increasing the number to 10 has little effect. We therefore settledon 10 photon packets. Figure 12.
Convergence tests of the
Skirt
FUV output for different total numbers of photon packets used in the calculation. Olsen et al. C. DISTRIBUTION OF CELL PHYSICAL PARAMETERSFigure 13 shows the volume-averaged gas and SFR densities for all cells in all galaxies of our sample. The gridpoint values in (cid:104) n H (cid:105) V and (cid:104) n SFR (cid:105) V used to sample the Cloudy grid are also shown in Figure 13. The cells in the simba galaxies have a larger spread in (cid:104) n SFR (cid:105) V and go to lower (cid:104) n H (cid:105) V than the chosen grid points, reflecting a largerparameter space than what is found in arepo -M51. However, we do not expect this to be a problem since the effectof (cid:104) n SFR (cid:105) V on the density PDF is lesser than that of (cid:104) n H (cid:105) V as seen in Figure 14 and the regions of low (cid:104) n H (cid:105) V aretreated separately. Figure 13.
Contour plot showing the cell distribution in density and SFR volume density for all 400 sample galaxies. Thehorizontal dashed lines indicate the center of the (cid:104) n SFR (cid:105) V bins shown in Figure 14, while the vertical dashed lines correspond tothe densities used to sample the simba gas densities in Section 4.4. Cells with (cid:104) n SFR (cid:105) V = 0 M (cid:12) yr − kpc − (cid:104) n SFR (cid:105) V = 10 − M (cid:12) yr − kpc − IGAME v3 D. COMPILATION OF PDFS FROM arepo -M51
Figure 14.
Volume-weighted gas density PDFs from the selected arepo -M51 run and fitted functions constructed as describedin Section 4.3. Each panel indicates a specific bin in gas density, (cid:104) n H (cid:105) V , and SFR volume density, (cid:104) n SFR (cid:105) V . The grey contourscorrespond to the 1- σ spread around the mean volume-weighted PDF (dashed teal curve) made from AREPO gas cells withinregions of (200 pc) . The solid teal curves are lognormal+power-law fits to the mean PDF. A vertical dashed red line indicatesthe critical density above which sink particles may form at which point the density PDF from gas cells alone is no longercomprehensive. The orange dashed curves account for the dense gas in sink particles by adding this mass at the high densitiesin the form of a modified power-law slope, and the orange solid curves are the resulting log-normal+power-law fit to the newdistribution. In cases where the modified power-law slope turns positive (last panel), the algorithm is instead allowed to searchfor a new log-normal center, keeping the power-law slope fixed at − . Olsen et al. E. COMPARISON OF SIMILAR GALAXIES IN simba -100 AND simba -25In order to make a fair comparison of the line emission calculated by s´ıgame for the different simba volumes usedhere, we have selected a handful of galaxies in both simba -100 and simba -25 that have similar global properties interms of M (cid:63) , M gas , SFR and (cid:104) Z (cid:105) SFR . We identify “galaxy pairs”, by searching for the smallest distance in the 4Dparameter space spanned by the M (cid:63) , M gas , SFR and (cid:104) Z (cid:105) SFR values, all in log units and normalized to lie in therange from 0 to 1. The result can be seen in Figure 15 in which the luminosity in simba -100 is subtracted from theluminosity in simba -25 for each galaxy pair to give a deviation in luminosity, as function of M (cid:63) , M gas , SFR and (cid:104) Z (cid:105) SFR . The colored lines show the closest pair for the parameter of that panel while the shaded regions show therange in luminosities from the 25th to 75th percentile for all galaxy pairs in that bin. A dashed horizontal line showsa deviation of 0, signalling that both galaxies have equal line luminosity. The vast majority of deviations are within1 dex, and we attribute the larger deviations at low values of M (cid:63) and M gas to the lower mass resolution in simba -100compared to simba -25. To illustrate the lack in resolution, a vertical dashed line indicates the minimum gas massin our simba -100 sample with at least 500 gas particles which is 5 . × M (cid:12) . The mass resolution of simba -25 isa factor of 8 higher, with the initial gas particle mass in simba -100 1 . × M (cid:12) , compared to 2 . × M (cid:12) in simba -25. Our simba -25 selection contains galaxies with gas masses down to 7 . × M (cid:12) , with the constraint thatthey contain at least 500 gas particles. Figure 15.
Deviation in line luminosity between galaxies of similar M (cid:63) , M gas , SFR and (cid:104) Z (cid:105) SFR in simba -100 vs simbasimba