Simulating the dust content of galaxies: successes and failures
Ryan McKinnon, Paul Torrey, Mark Vogelsberger, Christopher C. Hayward, Federico Marinacci
MMon. Not. R. Astron. Soc. , 1–17 (???) Printed 25 May 2017 (MN L A TEX style file v2.2)
Simulating the dust content of galaxies: successes and failures
Ryan McKinnon (cid:63) , Paul Torrey , , Mark Vogelsberger , Christopher C. Hayward , † ,and Federico Marinacci Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA TAPIR, California Institute of Technology, Pasadena, CA 91125, USA Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138, USA
Accepted ???. Received ???; in original form ???
ABSTRACT
We present full volume cosmological simulations using the moving-mesh code
AREPO tostudy the coevolution of dust and galaxies. We extend the dust model in
AREPO to includethermal sputtering of grains and investigate the evolution of the dust mass function, the cosmicdistribution of dust beyond the interstellar medium, and the dependence of dust-to-stellarmass ratio on galactic properties. The simulated dust mass function is well-described by aSchechter fit and lies closest to observations at z = 0 . The radial scaling of projected dustsurface density out to distances of Mpc around galaxies with magnitudes < i < issimilar to that seen in Sloan Digital Sky Survey data, albeit with a lower normalisation. At z = 0 , the predicted dust density of Ω dust ≈ . × − lies in the range of Ω dust values seenin low-redshift observations. We find that dust-to-stellar mass ratio anti-correlates with stellarmass for galaxies living along the star formation main sequence. Moreover, we estimate the µ m number density functions for simulated galaxies and analyse the relation betweendust-to-stellar flux and mass ratios at z = 0 . At high redshift, our model fails to produceenough dust-rich galaxies, and this tension is not alleviated by adopting a top-heavy initialmass function. We do not capture a decline in Ω dust from z = 2 to z = 0 , which suggests thatdust production mechanisms more strongly dependent on star formation may help to producethe observed number of dusty galaxies near the peak of cosmic star formation. Key words: dust, extinction – galaxies: evolution – galaxies: ISM – methods: numerical
The dust content of high-redshift galaxies provides insight into starformation and metal enrichment at early times, and the abundanceof dusty, starburst galaxies at submillimetre wavelengths (Smail,Ivison & Blain 1997; Barger et al. 1998; Hughes et al. 1998; Blainet al. 1999b; Eales et al. 1999; Scott et al. 2002) has implicationsfor theories of galaxy formation and evolution (Blain et al. 1999a;Chary & Elbaz 2001; Dunne, Eales & Edmunds 2003; Haywardet al. 2013; Casey, Narayanan & Cooray 2014). Models are chal-lenged to explain the presence of such galaxies and the key en-vironmental factors that contribute to their growth. Highlightingthis difficulty, there are recent observations of dusty galaxies atextremely high redshift, including HFLS3 at z = 6 . with dustmass M dust = 1 . × M (cid:12) (Riechers et al. 2013), A1689-zD1 at z = 7 . with M dust = 4 × M (cid:12) (Watson et al. 2015), and twogravitationally-lensed dusty sources at z = 5 . (Vieira et al. 2013;Hezaveh et al. 2013; Weiß et al. 2013).While (ultra)luminous infrared galaxies are roughly a thou- (cid:63) E-mail: [email protected] † Moore Prize Postdoctoral Scholar in Theoretical Astrophysics sand times more abundant at high redshift ( z ∼ − ) than atlow redshift (Chapman et al. 2005; Lagache, Puget & Dole 2005),not all high-redshift star-forming galaxies are dust-rich. A promi-nent example is Himiko, a z = 6 . galaxy with star-formationrate (SFR) roughly M (cid:12) yr − but very weak dust emission(Ouchi et al. 2013). The fact that some actively star-forming galax-ies are dust-rich while others are dust poor motivates a closer studyof high-redshift galaxies to better understand their formation.One important statistic is the dust mass function (DMF),whose evolution in time tracks dust growth across large popula-tions of galaxies. The DMF was first measured at low redshift aspart of the SCUBA Local Universe Galaxy Survey (Dunne et al.2000; Dunne & Eales 2001; Vlahakis, Dunne & Eales 2005). Evo-lution in the DMF has been studied over < z < (Eales et al.2009; Dunne et al. 2011; Clemens et al. 2013), and observationsfrom the Herschel ATLAS (Eales et al. 2010) find that the largestgalaxies at z = 0 . contained roughly five times more dust thanthose in the local universe (Dunne et al. 2011). The DMF has beenestimated for < z < using observations and number counts ofsubmillimetre galaxies, with dust-rich galaxies showing the mostchange compared to the present day (Dunne, Eales & Edmunds2003). Given the correlation between bolometric luminosity or SFR c (cid:13) ??? RAS a r X i v : . [ a s t r o - ph . GA ] M a y R. McKinnon et al. and dust obscuration (Wang & Heckman 1996; Adelberger & Stei-del 2000; Reddy et al. 2006), the evolution of the DMF is connectedto changes in the luminosity function. Luminosities of star-forming,dust-obscured galaxies at high redshift have been analysed in sur-vey data (Reddy et al. 2006; Dey et al. 2008; Magdis et al. 2012;Magnelli et al. 2012; Lo Faro et al. 2013; Sklias et al. 2014), andgalaxies at z ∼ are noticeably more luminous than their localcounterparts for fixed dust obscuration (Reddy et al. 2010). How-ever, the DMF remains less studied than statistics like the galaxystellar mass function, particularly in the high-redshift regime whereobservations are challenging.To approach this problem from a theoretical perspective, anumber of models have been developed to study the populationof submillimetre galaxies. Many of these models employ radiativetransfer to self-consistently track absorption and reradiation of stel-lar light by dust (e.g. using the GRASIL (Silva et al. 1998),
SUNRISE (Jonsson 2006),
RADISHE (Chakrabarti & Whitney 2009), or
ART (Yajima et al. 2012) codes) and to estimate submillimetre flux den-sities and number counts. Radiative transfer can be combined withsemi-analytic models (Baugh et al. 2005; Swinbank et al. 2008) orhydrodynamical simulations of galaxies (Chakrabarti et al. 2008;Narayanan et al. 2009, 2010; Hayward et al. 2011, 2012, 2013)to investigate how various galactic properties impact submillimetregalaxies.Such simulations have shown that flux densities in the SCUBA µ m (Holland et al. 1999) and AzTEC . mm (Wilson et al.2008) bands can be well estimated from a galaxy’s SFR and dustmass (Hayward et al. 2011). This agrees with findings that dustobscuration correlates with SFR (Adelberger & Steidel 2000). Fur-thermore, a top-heavy initial mass function (IMF), at least in star-bursts, may help to explain number counts of submillimetre galax-ies and their dust content (Baugh et al. 2005; Swinbank et al. 2008;Michałowski, Watson & Hjorth 2010).The predictive capability of such semi-analytic and radiativetransfer models motivates the inclusion of dust physics directly intogalaxy formation simulations where more diverse samples of galax-ies can be studied and the evolution of quantities like the DMFcan be traced. The direct treatment of dust in cosmological sim-ulations of uniform volumes provides the opportunity to investi-gate which environmental factors most contribute to the forma-tion of dusty, submillimetre galaxies. It also enables comparisonwith a variety of observations that cannot be fully tested in simu-lations of individual galaxies. These include the radial scaling ofprojected dust surface density around galaxies to distances of sev-eral Mpc (M´enard et al. 2010), the relation between SFR, stellarmass, and dust mass at low redshift and out to z = 2 . (da Cunhaet al. 2010; Dunne et al. 2011; Skibba et al. 2011; Bourne et al.2012; Cortese et al. 2012; Davies et al. 2012; Rowlands et al. 2012;Smith et al. 2012; Clemens et al. 2013; Santini et al. 2014; R´emy-Ruyer et al. 2015), and estimates of the cosmic dust density pa-rameter Ω dust and its evolution (Fukugita & Peebles 2004; Driveret al. 2007; M´enard et al. 2010; Dunne et al. 2011; Fukugita 2011;De Bernardis & Cooray 2012; M´enard & Fukugita 2012; Clemenset al. 2013; Thacker et al. 2013).In previous work (McKinnon, Torrey & Vogelsberger 2016,hereafter M16), we introduced a dust model accounting for theproduction of dust through stellar evolution, accretion in the in-terstellar medium (ISM) via collisions with gas-phase metals, andnon-thermal sputtering in supernova (SN) shocks that returned dustto the gas phase, and performed zoom-in simulations of a suiteof eight Milky Way-sized haloes. Here, we extend the dust model from M16 and perform the first cosmological simulations of galaxypopulations in which dust is directly treated.This paper is organized as follows. In Section 2, we describethe inclusion of new physics into our existing galaxy formationmodel and detail the initial conditions used for our simulations. InSection 3, we present our results and compare with existing data,and, in Section 4, we discuss the implications of our findings ina broader context. Finally, Section 5 summarizes our results andoffers an outlook on future work. We perform cosmological simulations using the moving-mesh code
AREPO (Springel 2010). The simulations incorporate the galaxyformation physics described in Vogelsberger et al. (2013). Briefly,this galaxy formation model includes gravity, hydrodynamics, pri-mordial and metal-line cooling (Wiersma, Schaye & Smith 2009),black hole growth (Sijacki et al. 2007), star formation (Springel& Hernquist 2003), stellar evolution, chemical enrichment trackingnine elements (H, He, C, N, O, Ne, Mg, Si, and Fe), and stellar andactive galactic nuclei feedback. It has been used in previous cosmo-logical simulations, including the Illustris simulation, that trace theevolution of the galaxy stellar mass function, luminosity function,mass-metallicity relation, and other quantities and shows broadagreement with observations (Vogelsberger et al. 2012; Torrey et al.2014; Vogelsberger et al. 2014a,b; Genel et al. 2014). In addition tothis galaxy formation model, we employ a modified version of thedust model from M16, which is described below and changes ourtreatment of dust in the circumgalactic medium (CGM).
The dust model in M16 accounts for dust production from agingstellar populations, grain growth, destruction in SN shocks, and theadvection and transport of dust in galactic winds. Dust is injectedinto the ISM as stars evolve off the main sequence, with dust massescalculated using stellar nucleosynthetic yields and estimated graincondensation efficiencies. The time-scale for grain growth throughcollisions between gas-phase atoms and grains depends on localgas density and temperature, while the time-scale for dust destruc-tion through SN sputtering scales inversely with the local SNe rate.Here, we also model the evolution of dust in galactic haloes. Thephysics of dust grains in hot gas has been studied in detail and in-cludes sputtering, cooling, and grain-grain collisions (Ostriker &Silk 1973; Burke & Silk 1974; Salpeter 1977; Barlow 1978; Draine& Salpeter 1979b; Itoh 1989; Tielens et al. 1994; Dwek, Foster &Vancura 1996; Smith et al. 1996). Thermal sputtering allows forthe erosion of dust grains by energetic atoms, and it can limit thedepletion of gas-phase metals onto grains in hot parts of a galactichalo (Burke & Silk 1974; Barlow 1978; Draine & Salpeter 1979b)and possibly enrich the intergalactic medium with metals (Bianchi& Ferrara 2005). Thermal sputtering affects grain lifetimes and gascooling in the intracluster medium (Yahil & Ostriker 1973; Dwek& Arendt 1992; McGee & Balogh 2010). Hydrogen and helium arethe main sputtering agents, and predictions of thermal sputteringrates indicate that sputtering overwhelms dust growth via accretionof gas-phase atoms for K < T < K (Draine & Salpeter1979b). The strength of thermal sputtering is expected to declinesharply below T ∼ K (Ostriker & Silk 1973; Barlow 1978;Draine & Salpeter 1979b; Tielens et al. 1994; Nozawa, Kozasa &Habe 2006). c (cid:13) ??? RAS, MNRAS000
The dust model in M16 accounts for dust production from agingstellar populations, grain growth, destruction in SN shocks, and theadvection and transport of dust in galactic winds. Dust is injectedinto the ISM as stars evolve off the main sequence, with dust massescalculated using stellar nucleosynthetic yields and estimated graincondensation efficiencies. The time-scale for grain growth throughcollisions between gas-phase atoms and grains depends on localgas density and temperature, while the time-scale for dust destruc-tion through SN sputtering scales inversely with the local SNe rate.Here, we also model the evolution of dust in galactic haloes. Thephysics of dust grains in hot gas has been studied in detail and in-cludes sputtering, cooling, and grain-grain collisions (Ostriker &Silk 1973; Burke & Silk 1974; Salpeter 1977; Barlow 1978; Draine& Salpeter 1979b; Itoh 1989; Tielens et al. 1994; Dwek, Foster &Vancura 1996; Smith et al. 1996). Thermal sputtering allows forthe erosion of dust grains by energetic atoms, and it can limit thedepletion of gas-phase metals onto grains in hot parts of a galactichalo (Burke & Silk 1974; Barlow 1978; Draine & Salpeter 1979b)and possibly enrich the intergalactic medium with metals (Bianchi& Ferrara 2005). Thermal sputtering affects grain lifetimes and gascooling in the intracluster medium (Yahil & Ostriker 1973; Dwek& Arendt 1992; McGee & Balogh 2010). Hydrogen and helium arethe main sputtering agents, and predictions of thermal sputteringrates indicate that sputtering overwhelms dust growth via accretionof gas-phase atoms for K < T < K (Draine & Salpeter1979b). The strength of thermal sputtering is expected to declinesharply below T ∼ K (Ostriker & Silk 1973; Barlow 1978;Draine & Salpeter 1979b; Tielens et al. 1994; Nozawa, Kozasa &Habe 2006). c (cid:13) ??? RAS, MNRAS000 , 1–17 imulating the dust content of galaxies We outline below the inclusion of thermal sputtering into thedust model used in M16. It is expected that other grain destruc-tion mechanisms, like grain-grain collisions and cosmic ray-drivensputtering, are subdominant compared to non-thermal SN shocksand thermal sputtering (Barlow 1978; Draine & Salpeter 1979a,b;Jones et al. 1994). We follow thermal sputtering prescriptions asused in previous galaxy modelling (Tsai & Mathews 1995; Hi-rashita et al. 2015) for simplicity of implementation.Following Equation 14 in Tsai & Mathews (1995), we esti-mate the sputtering rate for a grain of radius a in gas of density ρ and temperature T as d a d t = − (3 . × − cm s − ) (cid:18) ρm p (cid:19) (cid:20)(cid:18) T T (cid:19) ω + 1 (cid:21) − , (1)where m p is the proton mass, ω = 2 . controls the low-temperaturescaling of the sputtering rate, and T = 2 × K is the tempera-ture above which the sputtering rate is approximately constant. Thisempirical fitting formula approximately captures the temperaturedependence of sputtering rates derived in theoretical calculationsof collisions between spherical grains and impinging gas particles,which we outline in Appendix A (Barlow 1978; Draine & Salpeter1979b; Tielens et al. 1994). The associated sputtering time-scalefor the grain is given by Equation 15 in Tsai & Mathews (1995), τ sp = a (cid:12)(cid:12)(cid:12)(cid:12) d a d t (cid:12)(cid:12)(cid:12)(cid:12) − ≈ (0 . Gyr ) (cid:18) a − ρ − (cid:19) (cid:20)(cid:18) T T (cid:19) ω + 1 (cid:21) , (2)where a − is the grain size in units of . µ m and ρ − is the gasdensity in units of − g cm − , which corresponds to an effectivenumber density of n ≈ × − cm − . This time-scale is similarto the approximate sputtering time-scale given in Equation 44 ofDraine & Salpeter (1979b), where detailed projectile calculationswere performed. Given a grain of constant internal density ρ g andmass m g = 4 πa ρ g / , Equation (2) implies that mass changes ac-cording to the time-scale | m/ ˙ m | = τ sp / . In our model, we trackthe total dust mass for five chemical species (C, O, Mg, Si, and Fe)within each gas cell, but we do not track the grain size distribu-tion. To account for the effect of thermal sputtering on M i ,dust , thespecies i dust mass within each gas cell, during every time-step wecalculate the dust mass loss rate (cid:18) d M i ,dust d t (cid:19) sp = − M i ,dust τ sp / , (3)where τ sp is computed using Equation (2) and the local gas den-sity and temperature. We fix the grain radius a = 0 . µ m as ourmodel is not equipped to sample across a grain size distribution.This choice for a is motivated by the facts that the grain size dis-tribution for dust produced by AGB stars is thought to peak near . µ m (Groenewegen 1997; Winters et al. 1997; Yasuda & Kozasa2012; Asano et al. 2013b) and SNe are expected to form grains with a (cid:38) . µ m (Bianchi & Schneider 2007; Nozawa et al. 2007). Weshow in Appendix B that our results do not strongly depend on thisgrain size assumption.Combining with the dust accretion and SN-based destructionrates from Equations 4 and 5 of M16, the net rate of dust masschange is given by d M i ,dust d t = (cid:18) − M i ,dust M i ,metal (cid:19) (cid:18) M i ,dust τ g (cid:19) − M i ,dust τ d − M i ,dust τ sp / , (4)where M i ,metal is the metal mass of species i in the cell and thegrowth and destruction time-scales τ g and τ d depend on the localdensity, temperature, and Type II SN rate, as indicated by Equa-tions 5 and 7 in M16. This dust mass rate is computed on a cell-by- cell basis for every species and used to update dust masses in everytime-step.We summarize the set of parameters and quantities that char-acterize our fiducial dust model in Table 1. The dust model usedin this work differs from that of M16 in two respects: (i) it in-cludes thermal sputtering, and (ii) the dust growth parameters havebeen changed slightly to follow from Hirashita (2000), which of-fers a more detailed analysis of dust growth time-scales in molec-ular clouds. As shown in Section 3, this latter change was adoptedto lessen depletion at low redshift compared to the M16 model andbetter match the observed DMF and cosmic dust density parameter. We simulate uniformly-sampled cosmological volumes of comov-ing side length L = 25 h − Mpc with combined gas and dark mat-ter particle numbers of × , × , and × . The grav-itational softening length is held constant in comoving units until z = 1 , after which point it is fixed to the same physical value. Themaximum physical gravitational softening length is h − pc forthe run with × particles.We use Λ CDM cosmological parameters from the reanaly-sis of Planck data (Planck Collaboration et al. 2014) by Spergel,Flauger & Hloˇzek (2015) of Ω m = 0 . , Ω b = 0 . , Ω Λ = 0 . , σ = 0 . , n s = 0 . , and H =100 h km s − Mpc − = 68 km s − Mpc − . Initial conditions withthese parameters are generated at z = 127 using MUSIC (Hahn &Abel 2011) and iterated forward using
AREPO . We use the
SUB - FIND algorithm (Springel et al. 2001; Dolag et al. 2009) for iden-tifying gravitationally-bound structure and calculating gas, stellar,and dust mass components within galaxies. Galactic quantities arecomputed within twice the stellar half-mass radius.Table 2 provides details about the simulations, including soft-ening lengths and particle resolutions. Our fiducial simulations areperformed at three resolution levels and use the fiducial galaxy for-mation and feedback parameters from Vogelsberger et al. (2013),including the fiducial Chabrier (2003) IMF. We perform two ad-ditional simulations to explore how sensitive our results are to thedust model we use and to the choice of IMF.First, the “M16 model” run uses the dust model from M16,which lacked thermal sputtering and had dust growth parameterstuned to Milky Way-sized galaxies and resulting in fairly high de-pletion. The dust growth parameters used in this work and outlinedin Table 1 lead to weaker dust growth than in M16.Second, the “top-heavy IMF” run uses fiducial dust physicsbut an IMF of the form Φ( m ) ∝ m − . , with the same lower masslimit of . M (cid:12) and upper mass limit of M (cid:12) as the fiducialIMF. The m > M (cid:12) portion of the Chabrier (2003) IMF adoptsthe power law Φ( m ) ∝ m − . , and so the top-heavy IMF we ex-periment with increases the exponent by one and extends the powerlaw to the full mass range. While our top-heavy IMF is independentof galaxy properties, we note that previous works have used evenmore top-heavy IMFs in starbursts, like Φ( m ) ∝ m − (Baughet al. 2005; Swinbank et al. 2008). It is thought that a top-heavyIMF may help form large amounts of dust at high redshift. Becausethe stellar feedback prescription in Vogelsberger et al. (2013) istuned to a Chabrier (2003) IMF, for our top-heavy simulation we fixthe IMF-dependent quantities in the model of Springel & Hernquist(2003) to their fiducial values. To be precise, the parameters β , themass fraction of stars with m > M (cid:12) , and (cid:15) SN , the IMF-averagedenergy returned by supernovae per solar mass of stars formed, arekept at their values for a Chabrier (2003) IMF to ensure the stellar c (cid:13) ??? RAS, MNRAS , 1–17 R. McKinnon et al.
Table 1.
Summary of parameters used in various components of the full fiducial dust model. The dust condensation efficiencies that we use to computedust produced from stellar evolution are unchanged from those given in Table 2 of M16. The dust accretion parameters, which affect the growth time-scalecalculated in Equation 5 of M16, differ slightly from those used in M16 and are based on Equation 12 in Hirashita (2000).Parameter Value Descriptionthermal sputtering a . grain radius, in units of [ µ m ] dust accretion ρ ref . × − reference density roughly corresponding to n H = 100 cm − , in units of [ g cm − ] T ref reference temperature, in units of [ K ] τ refg . dust growth time-scale when T = T ref and ρ = ρ ref , in units of [ Gyr ] SN-based destruction E SNII, [10 erg ] Table 2.
Summary of simulation parameters and resolutions used in this work. Here, N is the total particle number, including equal numbers of dark matterand gas cells to start; (cid:15) is the maximum physical gravitational softening length, attained at z = 1 ; m dm is the dark matter resolution; and m gas is the targetgas mass for each cell in the (de-)refinement scheme. The last column describes the physics for each simulation. The M16 model refers to dust model used inM16, which lacked thermal sputtering and adopted dust growth parameters τ ref = 0 . Gyr, ρ ref = 2 . × − g cm − , and T ref = 20 K, leading to strongerISM dust growth. The top-heavy IMF run uses a pure power-law IMF of the form Φ( m ) ∝ m − . over the mass range from . M (cid:12) to M (cid:12) . It adoptsfiducial dust physics.Name Volume N (cid:15) m dm m gas physics [( h − Mpc ) ] [ h − kpc ] [ h − M (cid:12) ] [ h − M (cid:12) ] L25n128 × . . × . × fiducialL25n256 × .
25 6 . × . × fiducialL25n512 × .
625 8 . × . × fiducialM16 model × .
25 6 . × . × model used in M16 (see caption)top-heavy IMF × .
25 6 . × . × IMF has the form Φ( m ) ∝ m − . feedback model is not strongly affected by the top-heavy IMF. Intheory, to keep the Springel & Hernquist (2003) model consistentwith the top-heavy IMF, we would need to increase β , which wouldin turn affect the SN II rate and SN-driven sputtering of dust. How-ever, such changes would affect stellar feedback and its ability toreproduce the galaxy stellar mass function, which would compli-cate the interpretation of our results. To summarize, the top-heavyrun adopts fiducial dust physics (including thermal sputtering andthe Hirashita (2000) growth time-scale parameterisation) and formass return uses a power-law IMF of the form Φ( m ) ∝ m − . ,but it keeps the stellar feedback routines calibrated to the fiducialIMF. We first use our highest resolution fiducial simulation to visual-ize the distribution of dust and its redshift evolution. Figure 1shows projections of gas density, gas temperature, and dust den-sity through a slice of the full simulation volume at z = 2 , ,and . The dust surface density peaks in gas-rich halo centres,where the efficient production of dust by stars and short time-scalesfor grain-atom collisions overcome the presence of SN sputtering.Comparing different redshifts shows that the distribution of dust isrearranged through mergers, as demonstrated by the largest halo at z = 0 . It is also clear that large filaments of cold, diffuse gas farfrom potential minima – and thus sources of dust formation – haveessentially no dust. Recalling the temperature dependence of the sputtering time-scale given in Equation (2), we can see in Figure 1 that several ofthe largest haloes at z = 0 have temperatures above T ≈ K.At these temperatures, the thermal velocity is high enough to erodegrains. Lower mass haloes witness lower temperatures where thethermal sputtering rate falls off sharply. Regardless of halo size ortemperature, dust in the cool ISM is largely unaffected by thermalsputtering, and it is interesting to consider the DMF correspondingto this diverse sample of simulated galaxies.
We plot simulated DMFs at z = 2 . , . , and . in Figure 2 for ourfiducial runs at three different resolutions. Figure 2 also shows theDMFs for our two model variations, one using the dust model fromM16 and the other using a top-heavy IMF. Dust masses are com-puted within twice the stellar half-mass radius. We compare witha variety of observational data (Dunne, Eales & Edmunds 2003;Vlahakis, Dunne & Eales 2005; Eales et al. 2009; Dunne et al.2011; Clemens et al. 2013; Clark et al. 2015), although we notethat high-redshift observations are limited to the massive end of theDMF. These data have been corrected to the cosmology describedin Section 2.2. We also standardise the data to the dust mass ab-sorption coefficient κ (850 µ m ) = 0 . cm g − adopted in Dunne c (cid:13) ??? RAS, MNRAS000
We plot simulated DMFs at z = 2 . , . , and . in Figure 2 for ourfiducial runs at three different resolutions. Figure 2 also shows theDMFs for our two model variations, one using the dust model fromM16 and the other using a top-heavy IMF. Dust masses are com-puted within twice the stellar half-mass radius. We compare witha variety of observational data (Dunne, Eales & Edmunds 2003;Vlahakis, Dunne & Eales 2005; Eales et al. 2009; Dunne et al.2011; Clemens et al. 2013; Clark et al. 2015), although we notethat high-redshift observations are limited to the massive end of theDMF. These data have been corrected to the cosmology describedin Section 2.2. We also standardise the data to the dust mass ab-sorption coefficient κ (850 µ m ) = 0 . cm g − adopted in Dunne c (cid:13) ??? RAS, MNRAS000 , 1–17 imulating the dust content of galaxies z = 2 gas density temperature dust density z = 1 z = 0 − − Σ [M (cid:12) pc − ] 10 T [K] 10 − − − − Σ dust [M (cid:12) pc − ] Figure 1.
Projections of gas density, temperature, and dust density (left, middle, and right columns) at z = 2 , , and (top, middle, and bottom rows) for thehighest resolution simulation. Densities are given in physical units, and the scale bar for each redshift indicates a physical distance of Mpc. Projections wereperformed about the centre of the simulated volume, with a height and width of h − Mpc and a depth of . h − Mpc in comoving units. The distributionof dust largely traces that of gas. et al. (2011). We compare data from Dunne, Eales & Edmunds(2003) with simulated galaxies at z = 2 . , which is the value usedin that work to compute dust masses for galaxies without spectro- This choice of dust mass absorption coefficient means that dust mass datareported in Clemens et al. (2013), which used a value of κ (850 µ m ) smallerby roughly a factor of two, have been halved for this work. scopic redshifts. However, as noted in Section 3 of Dunne, Eales& Edmunds (2003), the estimated dust masses are largely insensi-tive to this choice of redshift. For the z = 1 . panel, we plot datafrom Eales et al. (2009) over the range . < z < . . FromVlahakis, Dunne & Eales (2005), we include both the directly-measured DMF and the DMF extrapolated over a larger dust massrange using IRAS PSCz data, the latter of which is given the suffix“ex” in the legend for Figure 2. c (cid:13) ??? RAS, MNRAS , 1–17 R. McKinnon et al. M dust [M (cid:12) ]10 − − − − − − Φ [ M p c − d e x − ] z = 2 . M dust [M (cid:12) ] z = 1 . M dust [M (cid:12) ] z = 0 . M dust [M (cid:12) ]10 − − − − − − Φ [ M p c − d e x − ] z = 2 . M dust [M (cid:12) ] z = 1 . M dust [M (cid:12) ] z = 0 . Figure 2.
Simulated DMFs (coloured lines) for three resolution levels (top row) and model variations (bottom row) as compared with observations (blackpoints) for z = 2 . , . , and . (left, middle, and right panels, respectively). For Eales et al. (2009), we plot data from . < z < . . From Vlahakis,Dunne & Eales (2005) we include both the directly measured DMF and the IRAS PSCz-extrapolated DMF, the latter of which has the suffix “ex” in the legend.For Dunne et al. (2011), we use the . < z < . set of data and cap the uncertainty at dex for two data points to improve readability. Observations havebeen corrected to conform to the cosmology detailed in Section 2.2. While the fiducial simulated DMFs offer a decent fit to observations at z = 0 . , they failto produce an abundance of dust-rich galaxies at higher redshift. While our fiducial model offers a reasonable fit to observeddata at z = 0 . down to the resolution limit, it does not repro-duce the abundance of high dust mass galaxies near z = 2 . and z = 1 . . At z = 2 . , the number density of simulatedgalaxies with M dust ≈ M (cid:12) is similar to that for observedgalaxies with M dust ≈ M (cid:12) . Although our simulated value of Φ( M dust ≈ × M (cid:12) ) increases by over dex from z = 2 . to . , we still have difficulty producing enough dust-rich galaxies at z = 1 . . The nature of dust processes makes the DMF behave ina much more dynamic way than the galaxy stellar mass function,since there is a diversity of ways for dust to grow (e.g. stellar injec-tion of grains and collisions with gas in the ISM) and be destroyed(e.g. SN shocks and thermal sputtering). This same core galaxyformation model without dust tracking had success in matching thegalaxy stellar mass function’s gradual flattening at the low massend from high to low redshift as more galaxies gain stellar mass(Vogelsberger et al. 2013; Torrey et al. 2014; Vogelsberger et al.2014a,b; Genel et al. 2014). However, the DMF does not evolve in such a monotonic fashion: galaxies at z = 2 . tend to be more dust-rich than at z = 0 . (Dunne, Eales & Edmunds 2003), and, evenfrom z = 0 . to z = 0 . , galactic dust masses decline by about afactor of five (Dunne et al. 2011). It is worth noting that for massbins of width . dex, a DMF value of Φ = 10 − Mpc − dex − corresponds to roughly two galaxies in our fiducial volume. Thus,the massive ends of our DMFs are sensitive to Poissonian statistics.In Appendix C, we simulate a volume eight times as large downto z = 2 . and investigate its DMF. The greater sample size pro-vided by a larger volume does not alleviate the absence of verydusty galaxies. Furthermore, the fiducial runs do not display robustconvergence as the resolution is increased. At z = 0 . , the DMFfalls off at the high-mass end more quickly with increasing resolu-tion. As a result, number densities associated with the L25n512 runlie above those for the L25n128 run in some mass bins, while thetrend is reversed in other bins. We note in Section 3.3 that volume-averaged quantities like comoving dust density display better con-vergence properties. c (cid:13) ??? RAS, MNRAS000
Simulated DMFs (coloured lines) for three resolution levels (top row) and model variations (bottom row) as compared with observations (blackpoints) for z = 2 . , . , and . (left, middle, and right panels, respectively). For Eales et al. (2009), we plot data from . < z < . . From Vlahakis,Dunne & Eales (2005) we include both the directly measured DMF and the IRAS PSCz-extrapolated DMF, the latter of which has the suffix “ex” in the legend.For Dunne et al. (2011), we use the . < z < . set of data and cap the uncertainty at dex for two data points to improve readability. Observations havebeen corrected to conform to the cosmology detailed in Section 2.2. While the fiducial simulated DMFs offer a decent fit to observations at z = 0 . , they failto produce an abundance of dust-rich galaxies at higher redshift. While our fiducial model offers a reasonable fit to observeddata at z = 0 . down to the resolution limit, it does not repro-duce the abundance of high dust mass galaxies near z = 2 . and z = 1 . . At z = 2 . , the number density of simulatedgalaxies with M dust ≈ M (cid:12) is similar to that for observedgalaxies with M dust ≈ M (cid:12) . Although our simulated value of Φ( M dust ≈ × M (cid:12) ) increases by over dex from z = 2 . to . , we still have difficulty producing enough dust-rich galaxies at z = 1 . . The nature of dust processes makes the DMF behave ina much more dynamic way than the galaxy stellar mass function,since there is a diversity of ways for dust to grow (e.g. stellar injec-tion of grains and collisions with gas in the ISM) and be destroyed(e.g. SN shocks and thermal sputtering). This same core galaxyformation model without dust tracking had success in matching thegalaxy stellar mass function’s gradual flattening at the low massend from high to low redshift as more galaxies gain stellar mass(Vogelsberger et al. 2013; Torrey et al. 2014; Vogelsberger et al.2014a,b; Genel et al. 2014). However, the DMF does not evolve in such a monotonic fashion: galaxies at z = 2 . tend to be more dust-rich than at z = 0 . (Dunne, Eales & Edmunds 2003), and, evenfrom z = 0 . to z = 0 . , galactic dust masses decline by about afactor of five (Dunne et al. 2011). It is worth noting that for massbins of width . dex, a DMF value of Φ = 10 − Mpc − dex − corresponds to roughly two galaxies in our fiducial volume. Thus,the massive ends of our DMFs are sensitive to Poissonian statistics.In Appendix C, we simulate a volume eight times as large downto z = 2 . and investigate its DMF. The greater sample size pro-vided by a larger volume does not alleviate the absence of verydusty galaxies. Furthermore, the fiducial runs do not display robustconvergence as the resolution is increased. At z = 0 . , the DMFfalls off at the high-mass end more quickly with increasing resolu-tion. As a result, number densities associated with the L25n512 runlie above those for the L25n128 run in some mass bins, while thetrend is reversed in other bins. We note in Section 3.3 that volume-averaged quantities like comoving dust density display better con-vergence properties. c (cid:13) ??? RAS, MNRAS000 , 1–17 imulating the dust content of galaxies In Figure 2, the dust model used in M16, which had a strongerdust growth mechanism and lacked thermal sputtering, differs themost from the fiducial model at z = 0 . and overproduces dust-rich galaxies. This is consistent with the finding in M16 that strongdust growth can overdeplete gas-phase metals at late times. Theseresults are largely unaffected by the inclusion of thermal sputteringsince the DMFs in Figure 2 isolate dust in the fairly cool ISM.However, the M16 model does predict more dust-rich galaxies at z = 2 . than the fiducial model and lies close to the Dunne, Eales& Edmunds (2003) data. A dust growth mechanism that allows formore variation among galaxies of different masses and SFRs maybe needed to form dust-rich galaxies at high redshift but also avoidoverproducing dust at low redshift.Similarly, the run with a top-heavy IMF, Φ( m ) ∝ m − . ,produces more dust than the fiducial L25n256 run at all redshifts.This is consistent with a top-heavy IMF shifting the galaxy stellarmass function towards lower masses due to shorter average stellarlifetimes. However, even this top-heavy IMF is unable to produceenough dust-rich galaxies at z = 2 . . This suggests that the ten-sion between our fiducial model and high-redshift observations ofmassive, dusty galaxies cannot be remedied by a variation in IMF.For the fiducial z = 0 . results, we fit data from the L25n512run with a Schechter function (Schechter 1976) of the form Φ( M dust ) ∆ M dust = Φ ∗ (cid:18) M dust M ∗ dust (cid:19) α exp (cid:18) − M dust M ∗ dust (cid:19) ∆ (cid:18) M dust M ∗ dust (cid:19) (5)to determine the best-fitting slope parameter α , characteristic dustmass M ∗ dust , and normalisation factor Φ ∗ . We obtain α = − . , M ∗ dust = 3 . × M (cid:12) , and Φ ∗ = 2 . × − Mpc − dex − .For comparison, the best-fitting Schechter function in Dunne et al.(2011) for . < z < . produces α = − . , M ∗ dust =3 . × M (cid:12) , and Φ ∗ = 5 . × − Mpc − dex − . Relativeto this observational data, the L25n512 run yields a similar slopeparameter, and though it predicts a lower turnover mass and highernormalisation factor, Figure 16 in Dunne et al. (2011) demonstrateshow these parameters are degenerate and anti-correlated. The visualizations in Figure 1 suggest that lines of sight far fromgalaxies suffer little dust extinction. We can directly quantify thisby considering the dust surface density in galactic haloes. Oneobservational technique to detect dust in haloes involves cross-correlating the brightness of quasars with the position of galaxiesto infer reddening from dust (M´enard et al. 2010). This correlationis used to estimate galactic reddening and infer dust surface densityprofiles, with the mean dust surface density following the scaling Σ dust ∝ r − . . This relation has been reproduced by analytic halomodels (Masaki & Yoshida 2012), and a similar technique has beenused to study the distribution of dust on larger scales in galaxy clus-ters (McGee & Balogh 2010).In Figure 3, we show the dust surface density profile as a func-tion of projected radial distance in physical units around galacticcentres at z = 0 . , averaging over all galaxies with < i < to match the magnitude cut used in M´enard et al. (2010). We cal-culate apparent magnitudes for simulated galaxies using the proce-dure outlined in Section 3.1 of Torrey et al. (2014). Briefly, we usethe stellar population synthesis model of Bruzual & Charlot (2003)and assign a luminosity to each star particle as a function of its age,initial stellar mass, and metallicity, and then we set each galaxy’s r [kpc]10 − − − − − Σ du s t [ M (cid:12) p c − ] L25n128L25n256L25n512Menard+ (2010)
Figure 3.
Dust surface density ( Σ dust ) as a function of projected radiusabout galactic centres at z = 0 . out to distances of Mpc in physicalunits. For each simulation, we show the mean dust surface density pro-file averaged across galaxies with < i < , using projections alongthe z axis of our box to ensure random orientations. This enables compar-ison with observational data from M´enard et al. (2010), shown in black,where galaxy position and quasar brightness correlations are used to inferreddening and SMC-type dust is assumed. The simulated dust surface den-sity scaling has lower normalisation than the observed result, particularly atlarge radii. luminosity to be the sum of the luminosities for constituent starparticles. To determine apparent magnitudes, we use the luminos-ity distance D L = 1598 Mpc for z = 0 . in our cosmology. Weperform projections for individual galaxies along the z axis of oursimulated box, resulting in random orientations with respect to theprojection axis. Every projection is carried out in a cylindrical vol-ume centered on the galactic potential minimum, using a radius of Mpc and a half-height of Mpc. Reducing this cylinder heightby a factor of two leaves the profiles in Figure 3 virtually unchangedfor radii less than Mpc and lowers them by by about . dex at themaximum radius of Mpc. The mean dust surface density profilefor all galaxies in this magnitude range is the result shown in Fig-ure 3. As noted in M´enard et al. (2010), on scales larger than thevirial radius, the dust surface density profile may be influenced bydust from surrounding or overlapping galaxies.The simulated dust surface density profiles appear well-converged out to r ≈ kpc, with the two highest-resolution runsshowing slightly greater dust surface density out to Mpc scales.Compared to the observed Σ dust ∝ r − . scaling, the simulatedprofiles are steeper for r < kpc and flatter for r > kpc.The dust surface density from our even highest resolution run stilllies below the observed data, with the tension largest for large radialdistances. In Figure 4, we show the comoving cosmic dust density ρ dust as afunction of redshift for our simulations and compare with observa-tional data at low redshift. The observational data have been cor- c (cid:13) ??? RAS, MNRAS , 1–17 R. McKinnon et al. z ρ du s t [ M (cid:12) M p c − ] L25n128L25n256L25n512Fukugita+ (2004)Driver+ (2007)Menard+ (2010)Dunne+ (2011)Fukugita (2011)De Bernardis+ (2012)Menard+ (2012)Clemens+ (2013)Thacker+ (2013) − − − − Ω du s t Figure 4.
Evolution of the comoving cosmic dust density ( ρ dust ; left axis)and associated dust density parameter ( Ω dust = ρ dust /ρ c ; right axis) as afunction of redshift for our three resolution simulations (coloured lines). Re-cent observations are shown in black (Fukugita & Peebles 2004; Driver et al.2007; M´enard et al. 2010; Dunne et al. 2011; Fukugita 2011; De Bernardis& Cooray 2012; M´enard & Fukugita 2012; Clemens et al. 2013; Thackeret al. 2013). Filled points include the contribution of dust from haloes, whileopen points and shaded regions track only galactic dust. The z = 0 value of Ω dust ≈ . × − is similar to observational estimates, though the simu-lated cosmic dust density doesn’t display the observed decline from z = 1 to z = 0 . rected as in Figure 2 to conform to the cosmology that we adopt forour simulations and, where appropriate, the dust mass absorptioncoefficient used in Dunne et al. (2011). In our simulations, the cos-mic dust density is computed by summing the dust masses of all gascells and dividing by the total comoving volume of (25 h − Mpc ) .We also show the cosmic dust density parameter Ω dust ≡ ρ dust /ρ c .In our fiducial simulations, the cosmic dust density increases byover dex from z = 5 to z = 0 , although there is very little evo-lution for z < . as the cosmic star formation rate density de-clines. Compared to the DMFs presented in Figure 2, the cosmicdust density results presented Figure 4 show stronger convergence.In particular, the L25n256 and L25n512 runs produce ρ dust valuesthat differ by less than . dex for z < . Even the low-resolutionL25n128 run displays the same qualitative behaviour with a lowernormalisation. This suggests that convergence is less of an issuewhen looking at volume-integrated dust quantities.At z = 0 , the L25n512 run reaches the values ρ dust ≈ × M (cid:12) Mpc − and Ω dust ≈ . × − , in rough agreementwith various low-redshift observations. In comparison, integratingthe best-fitting Schechter function for the L25n512 DMF yields Ω dust ≈ × − for the dust content of the ISM at z = 0 . However,we do not reproduce the observed decline in ρ dust by about a factorof three from z ≈ . to z = 0 seen in Herschel ATLAS data(Dunne et al. 2011). This decline has have a similar effect on theDMF for this redshift range, causing a drop in the Schechter func-tion parameter M ∗ dust and shifting the DMF to lower dust masses.The redshift behaviour of the cosmic dust density is not as well-studied as those of the cosmic star formation rate density and stel-lar mass density (e.g. see Madau & Dickinson 2014, and referencestherein) and would benefit from additional observations. We notethat observations of the stellar mass density ρ ∗ show an increase of more than . dex from z = 5 to z = 0 and a flattening for z < ,results similar to our simulated ρ dust evolution.Observational estimates of ρ dust at low redshift have beenobtained in a number of ways. One method includes fitting aSchechter function to DMF data and integrating it against dust massto find ρ dust = Γ(2 + α ) Φ ∗ M ∗ dust , where α , Φ ∗ , and M ∗ dust are thebest-fitting Schechter parameters (Dunne et al. 2011). Others as-sume a constant ratio between dust mass and B -band luminosityand calculate ρ dust by scaling the observed cosmic luminosity den-sity (Driver et al. 2007). The integrated dust density can also beestimated by transforming the luminosity function obtained fromphotometric surveys (Clemens et al. 2013) or derived from far-infrared power spectrum measurements (De Bernardis & Cooray2012; Thacker et al. 2013), or by combining a constant dust-to-metal ratio, mean ISM metallicity, and cool gas density parame-ter (Fukugita & Peebles 2004). We note, however, that these cal-culations tend to underestimate or neglect dust in galactic haloes,which is thought to contribute almost as much to Ω dust as ISM dust(M´enard et al. 2010; Fukugita 2011). Dust surface density profileslike in Figure 3, observationally obtained through quasar-galaxyreddening measurements (M´enard et al. 2010; M´enard & Fukugita2012), can be integrated out to the virial radius to estimate the halocomponent of dust mass and in turn a value of ρ dust that includescontributions from the ISM and CGM.The measurement of Ω dust by M´enard et al. (2010) in Fig-ure 4 accounts for dust in galactic haloes and lies above other low-redshift observations. This suggests that calculations of ρ dust and Ω dust using galactic DMFs tracing ISM luminosity or metallicitydata may be underestimating the true cosmic dust density, espe-cially in cases where galactic outflows can drive dust away fromthe ISM. The results in Figures 3 and 4 also show that Σ dust and ρ dust are interconnected: the dust content of the ISM cannot be var-ied independently of the dust content in galactic haloes, and ρ dust is influenced by dust in both of these regions. To a large degree, ρ dust determines the normalisation of quantities like Σ dust and canbe used to put constraints on the typical dust surface density seenfor individual galaxies. Figure 5 shows two-dimensional histograms indicating the averagedust mass of galaxies on and around the star formation main se-quence at z = 2 . and . . The average dust mass tends to increasewith both stellar mass and SFR as seen in both starburst (Mag-nelli et al. 2012; Santini et al. 2014) and local galaxies (Draineet al. 2007; Leroy et al. 2007; Kennicutt et al. 2009; Galametz et al.2011; Skibba et al. 2011; Fisher et al. 2013). For fixed stellar massor SFR, average dust mass increases from z = 2 . to z = 0 . , evenas the global SFR density and thus the stellar injection rate of dustdrops.We analyse the dependence of dust mass on stellar mass andSFR using a least-squares fit to the functional form log (cid:18) M dust M dust , (cid:19) = α log (cid:18) M ∗ M (cid:12) (cid:19) + β log (cid:18) SFRM (cid:12) yr − (cid:19) , (6)where M dust , , α , and β are free parameters. We apply this fit toall galaxies within σ of the star formation main sequence, us-ing the best-fitting relations from Figure 7 of Torrey et al. (2014).(The z = 2 . main sequence relation in that work is used forour z = 2 . panel.) We also impose the cut M ∗ > M (cid:12) toavoid galaxies that lie at the poorly-resolved end of the galaxystellar mass function. At z = 2 . , the best-fitting parameters are c (cid:13) ??? RAS, MNRAS000
Evolution of the comoving cosmic dust density ( ρ dust ; left axis)and associated dust density parameter ( Ω dust = ρ dust /ρ c ; right axis) as afunction of redshift for our three resolution simulations (coloured lines). Re-cent observations are shown in black (Fukugita & Peebles 2004; Driver et al.2007; M´enard et al. 2010; Dunne et al. 2011; Fukugita 2011; De Bernardis& Cooray 2012; M´enard & Fukugita 2012; Clemens et al. 2013; Thackeret al. 2013). Filled points include the contribution of dust from haloes, whileopen points and shaded regions track only galactic dust. The z = 0 value of Ω dust ≈ . × − is similar to observational estimates, though the simu-lated cosmic dust density doesn’t display the observed decline from z = 1 to z = 0 . rected as in Figure 2 to conform to the cosmology that we adopt forour simulations and, where appropriate, the dust mass absorptioncoefficient used in Dunne et al. (2011). In our simulations, the cos-mic dust density is computed by summing the dust masses of all gascells and dividing by the total comoving volume of (25 h − Mpc ) .We also show the cosmic dust density parameter Ω dust ≡ ρ dust /ρ c .In our fiducial simulations, the cosmic dust density increases byover dex from z = 5 to z = 0 , although there is very little evo-lution for z < . as the cosmic star formation rate density de-clines. Compared to the DMFs presented in Figure 2, the cosmicdust density results presented Figure 4 show stronger convergence.In particular, the L25n256 and L25n512 runs produce ρ dust valuesthat differ by less than . dex for z < . Even the low-resolutionL25n128 run displays the same qualitative behaviour with a lowernormalisation. This suggests that convergence is less of an issuewhen looking at volume-integrated dust quantities.At z = 0 , the L25n512 run reaches the values ρ dust ≈ × M (cid:12) Mpc − and Ω dust ≈ . × − , in rough agreementwith various low-redshift observations. In comparison, integratingthe best-fitting Schechter function for the L25n512 DMF yields Ω dust ≈ × − for the dust content of the ISM at z = 0 . However,we do not reproduce the observed decline in ρ dust by about a factorof three from z ≈ . to z = 0 seen in Herschel ATLAS data(Dunne et al. 2011). This decline has have a similar effect on theDMF for this redshift range, causing a drop in the Schechter func-tion parameter M ∗ dust and shifting the DMF to lower dust masses.The redshift behaviour of the cosmic dust density is not as well-studied as those of the cosmic star formation rate density and stel-lar mass density (e.g. see Madau & Dickinson 2014, and referencestherein) and would benefit from additional observations. We notethat observations of the stellar mass density ρ ∗ show an increase of more than . dex from z = 5 to z = 0 and a flattening for z < ,results similar to our simulated ρ dust evolution.Observational estimates of ρ dust at low redshift have beenobtained in a number of ways. One method includes fitting aSchechter function to DMF data and integrating it against dust massto find ρ dust = Γ(2 + α ) Φ ∗ M ∗ dust , where α , Φ ∗ , and M ∗ dust are thebest-fitting Schechter parameters (Dunne et al. 2011). Others as-sume a constant ratio between dust mass and B -band luminosityand calculate ρ dust by scaling the observed cosmic luminosity den-sity (Driver et al. 2007). The integrated dust density can also beestimated by transforming the luminosity function obtained fromphotometric surveys (Clemens et al. 2013) or derived from far-infrared power spectrum measurements (De Bernardis & Cooray2012; Thacker et al. 2013), or by combining a constant dust-to-metal ratio, mean ISM metallicity, and cool gas density parame-ter (Fukugita & Peebles 2004). We note, however, that these cal-culations tend to underestimate or neglect dust in galactic haloes,which is thought to contribute almost as much to Ω dust as ISM dust(M´enard et al. 2010; Fukugita 2011). Dust surface density profileslike in Figure 3, observationally obtained through quasar-galaxyreddening measurements (M´enard et al. 2010; M´enard & Fukugita2012), can be integrated out to the virial radius to estimate the halocomponent of dust mass and in turn a value of ρ dust that includescontributions from the ISM and CGM.The measurement of Ω dust by M´enard et al. (2010) in Fig-ure 4 accounts for dust in galactic haloes and lies above other low-redshift observations. This suggests that calculations of ρ dust and Ω dust using galactic DMFs tracing ISM luminosity or metallicitydata may be underestimating the true cosmic dust density, espe-cially in cases where galactic outflows can drive dust away fromthe ISM. The results in Figures 3 and 4 also show that Σ dust and ρ dust are interconnected: the dust content of the ISM cannot be var-ied independently of the dust content in galactic haloes, and ρ dust is influenced by dust in both of these regions. To a large degree, ρ dust determines the normalisation of quantities like Σ dust and canbe used to put constraints on the typical dust surface density seenfor individual galaxies. Figure 5 shows two-dimensional histograms indicating the averagedust mass of galaxies on and around the star formation main se-quence at z = 2 . and . . The average dust mass tends to increasewith both stellar mass and SFR as seen in both starburst (Mag-nelli et al. 2012; Santini et al. 2014) and local galaxies (Draineet al. 2007; Leroy et al. 2007; Kennicutt et al. 2009; Galametz et al.2011; Skibba et al. 2011; Fisher et al. 2013). For fixed stellar massor SFR, average dust mass increases from z = 2 . to z = 0 . , evenas the global SFR density and thus the stellar injection rate of dustdrops.We analyse the dependence of dust mass on stellar mass andSFR using a least-squares fit to the functional form log (cid:18) M dust M dust , (cid:19) = α log (cid:18) M ∗ M (cid:12) (cid:19) + β log (cid:18) SFRM (cid:12) yr − (cid:19) , (6)where M dust , , α , and β are free parameters. We apply this fit toall galaxies within σ of the star formation main sequence, us-ing the best-fitting relations from Figure 7 of Torrey et al. (2014).(The z = 2 . main sequence relation in that work is used forour z = 2 . panel.) We also impose the cut M ∗ > M (cid:12) toavoid galaxies that lie at the poorly-resolved end of the galaxystellar mass function. At z = 2 . , the best-fitting parameters are c (cid:13) ??? RAS, MNRAS000 , 1–17 imulating the dust content of galaxies M ∗ [M (cid:12) ]10 − − − S F R [ M (cid:12) y r − ] z = 2 . M ∗ [M (cid:12) ] z = 0 . M du s t [ M (cid:12) ] Figure 5.
Star formation main sequence at z = 2 . (left) and z = 0 . (right) for our L25n512 run, where the colour of each bin denotes the average dust massof galaxies whose stellar mass and SFR fall in those intervals. At both redshifts, dust mass tends to increase with stellar mass and SFR, though the correlationis stronger for stellar mass. M dust , = 2 . × M (cid:12) , α = 0 . , and β = 0 . , while at z = 0 . they are M dust , = 3 . × M (cid:12) , α = 0 . , β = 0 . .Dust mass is largely predicted by stellar mass, although there isalso a weak scaling with SFR.Figure 5 suggests that it is reasonable to associate the dusti-est galaxies with those that are the most star-forming. Such a pro-cedure is used in some hydrodynamical simulations where dust isnot directly treated in order to study the highly-luminous, dust-richsubmillimetre population. For example, Dav´e et al. (2010) assumesa z ∼ submillimetre galaxy number density of . × − Mpc − based on observations (Chapman et al. 2005; Tacconi et al. 2008)and finds that an SFR cut of around M (cid:12) yr − produces a galaxypopulation with a similar number density. This highly star-formingpopulation is taken to be the submillimetre set.However, we caution that the relation between dust mass, stel-lar mass, and SFR is complex, especially for submillimetre galax-ies. Our boxes of side length h − Mpc have difficulty captur-ing the nuclear starbursts and main sequence outliers (Sparre et al.2015; Sparre & Springel 2016) that tend to simultaneously increasethe SFR and decrease the dust mass (Hayward et al. 2011). Radia-tive transfer predicts that submillimetre flux scales more stronglywith dust mass than SFR, and Figure 5 indicates dust mass is moststrongly predicted by stellar mass. Comparing a low-stellar massstarburst with a higher-mass, lower-SFR main sequence galaxy, thelatter may have higher submillimetre flux because its increased dustmass more than compensates for its smaller SFR.
A galaxy’s dust-to-stellar mass ratio can increase not only throughdust injected into the ISM during stellar evolution, but also throughsubsequent dust growth in collisions with interstellar gas. Chem-ical evolution models (e.g. based on the work in Edmunds 2001)suggest that stellar injection of dust by itself – with no ISM dustgrowth – only produces dust-to-stellar mass ratios around − or less. Even in the extreme scenario where SNe produce more dustthan observed and condense nearly all ejected metals, this only re-sults in dust-to-stellar mass ratios near − (Dunne et al. 2011;Bourne et al. 2012). We predict some galaxies with dust-to-stellarmass ratios around − , and such dust-to-stellar mass ratios havebeen observed (Bourne et al. 2012; Cortese et al. 2012). Unless dustcondensation efficiencies are much higher than expected, this sug-gests that ISM dust growth is an important contributor to high dust-to-stellar mass ratios. We note that previous works have analyzedthe relative strengths of interstellar dust growth and stellar injec-tion of grains for increasing a galaxy’s dust-to-gas ratio (Mattsson,Andersen & Munkhammar 2012; Mattsson & Andersen 2012). Bystudying a population of galaxies, we can investigate both galaxieswhose dust-to-stellar mass ratios are driven by stellar injection ofgrains and by ISM dust growth.In Figure 6, we show the distribution of our simulated galaxiesas a function of dust-to-stellar mass ratio and stellar mass as well asdust-to-stellar mass ratio and specific star-formation rate (sSFR) at z = 2 . , . , and . . We compare with multiple sources of obser-vational data, detailed in Table 3 and meant to capture a variety ofmorphological types, colours, and metallicities, and note that sSFRanti-correlates with stellar mass (Brinchmann et al. 2004; Salimet al. 2007; Karim et al. 2011; Whitaker et al. 2012; Abramsonet al. 2014; Knebe et al. 2015). Several of the data sets in Table 3were already binned across stellar mass or sSFR or provided a meanresult which we show, together with quoted uncertainties, withoutmodification in Figure 6. For those unbinned data sets comprisedof numerous individual galaxy observations, we manually bin thedata to improve plot readability and compute sample standard de-viations in log-space to obtain symmetric error bars for Figure 6.The median sSFR drops by over dex from z = 2 . to . as the median dust-to-stellar mass ratio is largely unchanged. Thescatter in dust-to-stellar mass ratio at fixed sSFR slightly increasestowards low redshift. However, the slope of the dust-to-stellar massratio versus stellar mass relation does not change appreciably from c (cid:13) ??? RAS, MNRAS , 1–17 R. McKinnon et al.
Table 3.
Observational references with dust-to-stellar mass ratio data shown in Figure 6. We provide an approximate redshift range corresponding to ourcosmology for each sample and list the redshift at which data is plotted in Figure 6. In the last column, we briefly characterize each sample of galaxies andclarify which data we use. Several references provided already-binned data with uncertainties capturing scatter about the mean. For those references thatprovided quantities on a galaxy-by-galaxy basis, we binned data ourselves, calculating uncertainties in log-space to provide symmetric error bars for Figure 6.Reference Abbreviation Redshift Range Redshift Panel Notesda Cunha et al. (2010) D10 z < . z = 0 . we bin the sample of galaxies observed in all four IRAS bandsDunne et al. (2011) D11 z < . z = 0 . we use the mean result for these late-type galaxiesSkibba et al. (2011) S11 z < . z = 0 . we use the mean result for galaxies of all morphological typesBourne et al. (2012) B12-B z < . z = 0 . already-binned sample of galaxies with blue g − r colourB12-G z < . z = 0 . already-binned sample of galaxies with green g − r colourB12-R z < . z = 0 . already-binned sample of galaxies with red g − r colourCortese et al. (2012) C12-N z < . z = 0 . already-binned sample of H I -normal galaxiesC12-D z < . z = 0 . already-binned sample of H I -deficient galaxiesDavies et al. (2012) D12 z < . z = 0 . we bin this sample of bright galaxiesRowlands et al. (2012) R12 z < . z = 0 . we bin the sample of early-type galaxiesSmith et al. (2012) S12 z < . z = 0 . we bin the sample of early-type galaxies, excluding non-detectionsSantini et al. (2014) S14 . < z < . z = 1 . we bin this sample of galaxiesS14 . < z < . z = 2 . we bin this sample of high-redshift galaxiesR´emy-Ruyer et al. (2015) R15 z < . z = 0 . we bin the sample covering a roughly dex metallicity range M ∗ [M (cid:12) ]10 − − − − − M du s t / M ∗ z = 2 . M ∗ [M (cid:12) ] z = 1 . M ∗ [M (cid:12) ] z = 0 . − − − − − sSFR [yr − ]10 − − − − − M du s t / M ∗ z = 2 . − − − − − sSFR [yr − ] z = 1 . − − − − − − sSFR [yr − ] z = 0 . Figure 6.
Dust-to-stellar mass ratio ( M dust /M ∗ ) as a function of stellar mass (top row) and specific star-formation rate (sSFR; bottom row) at z = 2 . , . ,and . (left, middle, and right panels). For each redshift, the logarithmic distribution of simulated galaxies in the L25n512 run is given by a two-dimensionalhistogram, with bluer colours indicating greater density. Dotted black lines mark the median value in the distribution for each axis. References for the binnedobservational data (red points) are given in Table 3. Dust-to-stellar mass ratio anti-correlates with stellar mass at both high and low redshift.c (cid:13) ??? RAS, MNRAS000
Dust-to-stellar mass ratio ( M dust /M ∗ ) as a function of stellar mass (top row) and specific star-formation rate (sSFR; bottom row) at z = 2 . , . ,and . (left, middle, and right panels). For each redshift, the logarithmic distribution of simulated galaxies in the L25n512 run is given by a two-dimensionalhistogram, with bluer colours indicating greater density. Dotted black lines mark the median value in the distribution for each axis. References for the binnedobservational data (red points) are given in Table 3. Dust-to-stellar mass ratio anti-correlates with stellar mass at both high and low redshift.c (cid:13) ??? RAS, MNRAS000 , 1–17 imulating the dust content of galaxies z − − M du s t / M ∗ all galaxies10 M (cid:12) ≤ M ∗ < M (cid:12) M (cid:12) ≤ M ∗ < M (cid:12) M (cid:12) ≤ M ∗ < M (cid:12) M (cid:12) ≤ M ∗ < M (cid:12) Dunne+ (2011)
Figure 7.
Evolution of the galactic dust-to-stellar mass ratio ( M dust /M ∗ )as a function of redshift for the L25n512 run. We compute the median dust-to-stellar mass ratio for all galaxies (red line) and for different stellar massbins (blue lines, with deeper shades indicating greater stellar mass), exclud-ing galaxies with no stellar component. Observational data using galaxieswith spectroscopic redshifts from the Herschel ATLAS (Dunne et al. 2011)are shown in black, though we note these observations were not binned bymass. The median dust-to-stellar mass ratio increases by about a factor of1.7 from z = 5 to z = 0 , with the most massive galaxies displaying lessevolution in dust-to-stellar mass ratio. z = 2 . to z = 0 . . At high redshift, the observed dust-to-stellarmass ratios for M ∗ (cid:38) M (cid:12) from Santini et al. (2014) are largerthan what we predict. However, we note that most of the high-redshift galaxies in Santini et al. (2014) have SFR (cid:38) M (cid:12) yr − ,making them more star-forming than nearly all galaxies in our sim-ulation.To better understand the connection between dust-to-stellarmass ratio and stellar mass, in Figure 7 we plot the median dust-to-stellar mass ratio as a function of redshift for dex stellar massbins ranging from M (cid:12) to M (cid:12) , along with the median ra-tio across all simulated galaxies. The results from Figure 6 – thatthe dust-to-stellar mass ratio decreases with increasing stellar mass,while the overall median dust-to-stellar mass ratio does not sub-stantially increase with time – are confirmed in Figure 7. The na-ture of the galaxy stellar mass function implies that at nearly everyredshift the median dust-to-stellar mass ratio lies nearer to the ra-tio for M (cid:12) < M ∗ < M (cid:12) galaxies than for more massivegalaxies.The overall median dust-to-stellar mass ratio increases by justunder a factor of two from z = 5 to , and it is largely flat for z < . Low stellar mass galaxies display similar behaviour, while thedust-to-stellar mass ratio for large galaxies with M ∗ > M (cid:12) isroughly dex below the value for galaxies with M (cid:12) < M ∗ < M (cid:12) . However, we do not capture the decrease in dust-to-stellarmass ratio by a factor of two observed in Herschel ATLAS data for z < . (Dunne et al. 2011), a result similarly shown in Figure 4.While the dust-to-stellar mass scaling at z = 0 is in decentagreement with observations, the relation between dust-to-gas ra-tio and gas-phase metallicity displays more tension and a greatersensitivity to the parameters of our model. In Figure 8, we plot the dust-metallicity relation at z = 0 for the L25n256 and M16model simulations and offer a comparison to observations (Leroyet al. 2011; R´emy-Ruyer et al. 2014) and modelling (Asano et al.2013a; Zhukovska 2014; Popping, Somerville & Galametz 2016).It is clear that the L25n256 run fails to match slope of the ex-pected dust-metallicity relation: despite offering a reasonable fit tothe z = 0 DMF in Figure 2, the dust-metallicity relation is far tooflat. The L25n256 run is similar to the semi-analytic model of Pop-ping, Somerville & Galametz (2016) for
12 + log( O / H ) < , butdeviates strongly at high metallicity. The dust growth timescale inlarge galaxies in the L25n256 run seems to be too long, preventingthem from rapidly growing their dust mass. On the other hand, theM16 model – which employs a stronger ISM dust growth mecha-nism and lacks thermal sputtering – displays a dust-metallicity re-lation whose slope better matches observations but whose normal-isation is too high. However, the M16 model significantly overpro-duces dust-rich galaxies in its z = 0 DMF. Figure 8 offers anotherlook at the DMF tension in Figure 2 and highlights the difficulty inproducing enough dust to match the dust-metallicity relation whileavoiding a DMF with too many dust-rich galaxies.
We can combine the direct dust mass tracking in our work with stel-lar population synthesis postprocessing to make predictions aboutthe observational properties of simulated galaxies. One such prop-erty is the dust-to-stellar flux ratio ( f dust /f ∗ ), which measures theflux reradiated by dust grains as a fraction of unextincted stellarflux and has been observed to correlate with dust mass and infraredluminosity (Skibba et al. 2011). Below, we use the FSPS (Conroy,Gunn & White 2009; Conroy & Gunn 2010) stellar population syn-thesis code to estimate each galaxy’s bolometric luminosity and inturn its dust-to-stellar flux ratio.The dependence of optical depth on host galaxy properties waspreviously modelled in Jonsson et al. (2006). Following Table 2 inJonsson et al. (2006), we estimate the bolometric attenuation to be τ = 0 . (cid:18) Z . (cid:19) . (cid:18) SFRM (cid:12) yr − (cid:19) . (cid:18) M b M (cid:12) (cid:19) − . , (7)where M b is the galaxy’s total baryon mass. We calculate τ directlyon a galaxy-by-galaxy basis using our SUBFIND output. We com-pute an unattenuated synthetic spectral energy distribution (SED)for every star particle as a function of its initial mass, age, andmetallicity, assuming a Chabrier (2003) IMF, and define a galaxy’sSED to be the sum of those from constituent star particles. For agalaxy with bolometric luminosity L and optical depth τ computedin Equation (7), we calculate the dust luminosity L dust using Equa-tion 6 in Jonsson et al. (2006), L dust /L = 1 − (1 /τ )(1 − e − τ ) . (8)The stellar luminosity is then L ∗ = L − L dust . In essence, stel-lar flux is computed by integrating the attenuated stellar SED, anddust flux is obtained by integrating over the difference betweenthe unattenuated and attenuated stellar SEDs. This calculation ofthe dust-to-stellar flux ratio is simpler than estimating and remov-ing stellar emission in the mid-infrared from a dust-extincted SED(Draine et al. 2007; Mu˜noz-Mateos et al. 2009; Skibba et al. 2011),although only possible when postprocessing simulated galaxies.Figure 9 shows the distribution of galaxies in the L25n512 runas a function of dust-to-stellar flux and mass ratios at z = 0 , withobservational data from the Herschel KINGFISH Survey overlaid(Skibba et al. 2011). While the range of dust-to-stellar mass ratios c (cid:13) ??? RAS, MNRAS , 1–17 R. McKinnon et al. . . . . .
012 + log(O / H)10 − − − − − − D L25n256R´emy-Ruyer+ (2014)Zhukovska+ (2014)Popping+ (2016) . . . . .
012 + log(O / H)10 − − − − − − D M16 modelLeroy+ (2011)Asano+ (2013)
Figure 8.
Simulated dust-metallicity relations at z = 0 for the fiducial L25n256 run (left) and M16 model (right). The two-dimensional histograms usea logarithmic colourscale to indicate number density, with bluer colours denoting a greater number of galaxies at a given dust-to-gas ratio and metallicity.Red points and lines denote observational data (Leroy et al. 2011; R´emy-Ruyer et al. 2014) and results from analytic and semi-analytic modelling (Asanoet al. 2013a; Zhukovska 2014; Popping, Somerville & Galametz 2016). To improve readability, we omit error bars. From Asano et al. (2013a) we use the τ SF = 5 Gyr model and from Zhukovska (2014) we use the “6x 500 Myr bursts τ SF = 2 Gyr” model, which were compiled by R´emy-Ruyer et al. (2014). − − f dust /f ∗ − − − − − M du s t / M ∗ Skibba+ (2011)
Figure 9.
Dust-to-stellar mass ratio ( M dust /M ∗ ) as a function of dust-to-stellar flux ratio ( f dust /f ∗ ) at z = 0 . The two-dimensional histogram showsthe number distribution of galaxies from the L25n512 run on a logarithmicscale, with bluer colours denoting greater counts. We compare with obser-vational data (red circles) from the Herschel KINGFISH Survey (Skibbaet al. 2011). does tend to match these observations, the simulated dust-to-stellarflux ratios are biased to smaller values than in the Herschel KING-FISH set. The two-dimensional distribution of simulated galaxiespeaks near M dust /M ∗ ≈ × − and f dust /f ∗ ≈ × − ,though there is significant scatter in the dust-to-stellar flux ratios,with a number of galaxies recording f dust /f ∗ > . The largest ten-sion with observations comes at low dust-to-stellar flux ratio, where we predict numerous galaxies with large dust-to-stellar mass ra-tio. From Figure 6, we know that galaxies with high dust-to-stellarmass ratios tend to have low stellar masses and thus low metallic-ities and SFRs. This may drive down the optical depths calculatedin Equation (7) and thus the dust-to-stellar flux ratios. Several ofthe galaxies in Skibba et al. (2011) with low dust-to-stellar fluxratios are early-types, which are known to have smaller dust-to-stellar mass ratios on average than spirals (Rowlands et al. 2012;Smith et al. 2012). In any case, the scatter in our simulated resultsdoes confirm the observation that the dust-to-stellar flux ratio cov-ers roughly three orders of magnitude and does not effectively con-strain the mass ratio (Skibba et al. 2011).Previous works have demonstrated how hydrodynamical sim-ulations without direct dust tracking can be coupled with radia-tive transfer to study galactic flux densities at submillimetre wave-lengths (Chakrabarti et al. 2008; Narayanan et al. 2009, 2010; Hay-ward et al. 2011, 2012, 2013). Performing dust radiative transfer onsimulations of isolated and merging disc galaxies, Hayward et al.(2011) developed fitting functions to estimate submillimetre fluxdensities in the SCUBA µ m and AzTEC . mm bands as afunction of SFR and dust mass as well as dust luminosity and dustmass. While these relations were derived from simulations inves-tigating number counts of bright submillimetre galaxies, here weapply them to our full sample of galaxies to demonstrate how cos-mological simulations can benefit from results obtained through ra-diative transfer calculations.To construct submillimetre number densities, we define theHubble parameter H ( z ) = H (cid:112) Ω m (1 + z ) + Ω Λ and comovingdistance l c ( z ) = (cid:90) z c d z (cid:48) H ( z (cid:48) ) . (9)If Φ( S, z ) denotes the comoving number density of galaxies withsubmillimetre flux S at redshift z per unit logarithmic flux, then φ ( S ) = (cid:16) π deg − (cid:17) (cid:90) ∞ πl c ( z ) Φ( S, z ) d l c ( z ) (10) c (cid:13) ??? RAS, MNRAS000
Dust-to-stellar mass ratio ( M dust /M ∗ ) as a function of dust-to-stellar flux ratio ( f dust /f ∗ ) at z = 0 . The two-dimensional histogram showsthe number distribution of galaxies from the L25n512 run on a logarithmicscale, with bluer colours denoting greater counts. We compare with obser-vational data (red circles) from the Herschel KINGFISH Survey (Skibbaet al. 2011). does tend to match these observations, the simulated dust-to-stellarflux ratios are biased to smaller values than in the Herschel KING-FISH set. The two-dimensional distribution of simulated galaxiespeaks near M dust /M ∗ ≈ × − and f dust /f ∗ ≈ × − ,though there is significant scatter in the dust-to-stellar flux ratios,with a number of galaxies recording f dust /f ∗ > . The largest ten-sion with observations comes at low dust-to-stellar flux ratio, where we predict numerous galaxies with large dust-to-stellar mass ra-tio. From Figure 6, we know that galaxies with high dust-to-stellarmass ratios tend to have low stellar masses and thus low metallic-ities and SFRs. This may drive down the optical depths calculatedin Equation (7) and thus the dust-to-stellar flux ratios. Several ofthe galaxies in Skibba et al. (2011) with low dust-to-stellar fluxratios are early-types, which are known to have smaller dust-to-stellar mass ratios on average than spirals (Rowlands et al. 2012;Smith et al. 2012). In any case, the scatter in our simulated resultsdoes confirm the observation that the dust-to-stellar flux ratio cov-ers roughly three orders of magnitude and does not effectively con-strain the mass ratio (Skibba et al. 2011).Previous works have demonstrated how hydrodynamical sim-ulations without direct dust tracking can be coupled with radia-tive transfer to study galactic flux densities at submillimetre wave-lengths (Chakrabarti et al. 2008; Narayanan et al. 2009, 2010; Hay-ward et al. 2011, 2012, 2013). Performing dust radiative transfer onsimulations of isolated and merging disc galaxies, Hayward et al.(2011) developed fitting functions to estimate submillimetre fluxdensities in the SCUBA µ m and AzTEC . mm bands as afunction of SFR and dust mass as well as dust luminosity and dustmass. While these relations were derived from simulations inves-tigating number counts of bright submillimetre galaxies, here weapply them to our full sample of galaxies to demonstrate how cos-mological simulations can benefit from results obtained through ra-diative transfer calculations.To construct submillimetre number densities, we define theHubble parameter H ( z ) = H (cid:112) Ω m (1 + z ) + Ω Λ and comovingdistance l c ( z ) = (cid:90) z c d z (cid:48) H ( z (cid:48) ) . (9)If Φ( S, z ) denotes the comoving number density of galaxies withsubmillimetre flux S at redshift z per unit logarithmic flux, then φ ( S ) = (cid:16) π deg − (cid:17) (cid:90) ∞ πl c ( z ) Φ( S, z ) d l c ( z ) (10) c (cid:13) ??? RAS, MNRAS000 , 1–17 imulating the dust content of galaxies − − − − S [mJy]10 φ [ d e g − d e x − ] L25n256M16 modelBlain+ (1999)Cowie+ (2002)Smail+ (2002)Knudsen+ (2008)Chen+ (2013)
Figure 10.
Simulated number density functions in the SCUBA µ mband for the fiducial L25n256 run (red) and the M16 model (green). Fluxesare computed using the luminosity- and dust mass-dependent fitting func-tions provided in Hayward et al. (2011). Black points mark observationsfor µ m (Blain et al. 1999a; Smail et al. 2002; Cowie, Barger & Kneib2002; Knudsen, van der Werf & Kneib 2008; Chen et al. 2013) as com-piled by Casey, Narayanan & Cooray (2014). Because the fits in Haywardet al. (2011) were not designed for z (cid:46) , we show two versions of thenumber density functions: one integrated all the way down to z = 0 (solidlines) and the other only down to z = 1 (dashed lines). The versions showsimilar behaviour. While the M16 model offers a better fit to the observedsubmillimetre number densities and the high-redshift DMF in Figure 2, itsignificantly overproduces dust in the z = 0 DMF. is the number of galaxies with flux S per square degree per unitlogarithmic flux. Simplifying, we calculate submillimetre numbercounts using φ ( S ) = (cid:16) π deg − (cid:17) (cid:90) ∞ πl c ( z ) Φ( S, z ) cH ( z ) d z. (11)In practice, we numerically integrate Equation (11) using Φ( S, z ) values constructed from simulation output at discrete redshifts.Number density functions for simulated galaxies at µ mare shown in Figure 10 for our fiducial L25n256 run and the M16model. We compare with various observational data (Blain et al.1999a; Smail et al. 2002; Cowie, Barger & Kneib 2002; Knud-sen, van der Werf & Kneib 2008; Chen et al. 2013). Fluxes inthis band are computed using dust luminosities and dust massesfollowing Equation 2 and Appendix A in Hayward et al. (2011).Dust luminosities are calculated as in the construction of Figure 9.Because the fits provided in Hayward et al. (2011) are not de-signed to apply to z (cid:46) , we compute number density functionsin two ways: one using z = 0 as the lower limit of the integra-tion in Equation (11) and the other using z = 1 . This variationchanges results only slightly. Results for this submillimetre bandshow that the number density of galaxies declines over the flux in-terval − mJy < S < − mJy accessible in both simulations.The M16 model offers a much better fit to the high-flux observa-tions than the L25n256 run, in large part because the high-redshiftDMF for M16 model in Figure 2 contains many more dust-richgalaxies. However, despite the more realistic submillimetre num- ber counts, the M16 model has tension of its own: its z = 0 DMFcontains far too many galaxies with M dust > M (cid:12) . This tensionhighlights the need to form enough dust at high redshift to generaterealistic submillimetre number counts while preventing an excessof dust at low redshift. Futhermore, as noted in our discussion ofthe star formation main sequence in Figure 5, the box length of L = 25 h − Mpc used in this work makes it difficult to truly probethe submillimetre regime and uncover possible exponential cutoffsin the submillimetre number density functions. This may be worthpursuing in cosmological simulations of larger volumes or in semi-analytic models with dust tracking.
We have presented full volume cosmological simulations with amodel for dust production and destruction to study the coevolutionof dust and galaxies. Our model offers rough agreement with low-redshift observations of the DMF, cosmic dust density, and relationbetween dust-to-stellar mass ratio and stellar mass, but it also high-lights limitations that appear more fundamental. Despite offeringa reasonable match to the z = 0 DMF, the fiducial model fails tocapture the abundance of dusty galaxies at high redshift and insteadproduces galaxies whose dust masses grow roughly in a monotonicfashion. It has been suggested that perhaps the dust-rich galaxies athigh redshift have extra-high star-formation efficiencies, are moreefficient at forming dust from stars, or feature more top-heavy IMFsthan low-redshift galaxies (Dunne et al. 2011). In this scenario, themost dusty galaxies at z = 2 . could evolve to the present withmuch lower dust masses after consuming their gas and dust in starformation. The ability of the fiducial model to roughly match the z = 0 DMF but not capture the decline in dusty galaxies fromhigh to low redshift suggests that dust evolution processes may bemore dependent on host galaxy properties like SFR or gas frac-tion than assumed in this work. For example, dust yields in stel-lar ejecta may be a function of local ISM density or temperature,which evolve with redshift. To account for the observed shift to-wards lower masses in the DMF from z = 2 . to z = 0 , we needefficient sputtering of grains in dust-rich galaxies. This would bet-ter enable different galaxies to have diverse dust mass histories andperhaps lead to DMF behaviour that more closely follows the cos-mic SFR evolution (e.g. see Madau & Dickinson 2014, and refer-ences therein).Our analysis of the DMF, dust surface density profiles, andcosmic dust density evolution also highlights the observational un-certainties that make it challenging to obtain reliable dust mass es-timates. For example, the dust surface density comparison in Fig-ure 3 shows that our simulated dust radial profiles are closest toobservations from SDSS out to r ≈ kpc but lie below observa-tions by up to dex at larger radii. The cosmic dust density in oursimulations could easily absorb a factor of two or three increase andstill be consistent with observations in Figure 4, and such a normal-isation change would help boost the dust surface density profiles ona global scale. This change is plausible since dust condensation ef-ficiencies in stellar ejecta and ISM growth time-scales are not well-constrained. Independent of such a normalisation shift, it is alsopossible that thermal sputtering is slightly too strong or galacticoutflows too weak, limiting the amount of dust in galactic haloes.The sample of SDSS galaxies used to construct surface den-sity profiles in M´enard et al. (2010) has a redshift distribution thatpeaks at z ≈ . but with a full-width at half-maximum of . .While we do not predict much evolution in the cosmic dust den- c (cid:13) ??? RAS, MNRAS , 1–17 R. McKinnon et al. sity for z (cid:46) , the fact that Figure 3 is constructed at a fixed red-shift of z ≈ . may introduce some deviation in surface densityprofiles from the observed result. The observed reddening signalwas also tested for possible systematic effects (e.g. by subsamplingquasars according to magnitude bins or using sky regions with dif-ferent Galactic reddening) and found to be robust. The calculationsin M´enard et al. (2010) assume SMC-type dust, motivated in partby the observation that few high-redshift galaxies share the . µ mextinction curve bump characteristic of the Milky Way, but adopt-ing Milky Way-type dust would change dust masses about a fac-tor of two. Given the difficulties in estimating dust masses fromreddening signals, including weak constraints on dust mass absorp-tion coefficients, the discrepancies between simulated and observeddust surface density profiles could be influenced by inaccuracies inmodelled physics like thermal sputtering or galactic outflows or byuncertain observational assumptions. The M´enard et al. (2010) re-lation can possibly be used as a constraint on outflow physics, asthe dust surface density profiles at large radii are likely sensitive tothe outflow model and the coupling of dust and gas in winds. In oursimulations, we assume winds have the same depletion as the ISMfrom which they are launched (e.g. if a wind particle is created in acell where 10% of metals are locked in dust, then 10% of the metalcontent of the wind particle is assumed to be dust), but alternativemodels that couple dust more strongly to winds may help drive dustto Mpc distances. Simulations by Zu et al. (2011) suggest that re-producing the M´enard et al. (2010) relation without galactic windsis difficult, and the strength of outflows could be used to constrainenrichment in the intergalactic medium.While dust masses can be hard to estimate, the findings inSection 3.5 demonstrate the connection between a galaxy’s stel-lar mass and its dust content. For example, Figure 7 provides amethod to calculate a galactic dust mass when only the stellar massand redshift are known, using the dust-to-stellar mass ratio in theappropriate stellar mass bin. Both Figures 6 and 7 highlight thedependence of dust-to-stellar mass ratio on stellar mass and whyassuming a uniform dust-to-stellar mass ratio is not ideal. Previousobservational studies have shown that gas fraction and moleculargas fraction decrease with stellar mass (Leroy et al. 2008; Daddiet al. 2010; Geach et al. 2011; Saintonge et al. 2011; Poppinget al. 2012; Bauermeister et al. 2013; Tacconi et al. 2013; Boselli,Cortese & Boquien 2014; Bothwell et al. 2014; Morokuma-Matsui& Baba 2015), and this result has been reproduced in galaxy for-mation simulations and semi-analytic models (Hopkins et al. 2009;Obreschkow et al. 2009; Dav´e et al. 2010; Lagos et al. 2011a,b;Duffy et al. 2012; Fu et al. 2012; Genel et al. 2014; Popping,Somerville & Trager 2014; Lagos et al. 2015; Narayanan et al.2015). Additionally, simple dust and chemical evolution modelssuggest that the dust-to-stellar mass ratio increases with gas frac-tion (Dunne et al. 2011), a result that has been seen in HerschelReference Survey data (Cortese et al. 2012).Together, these findings indicate the dust-to-stellar mass ratioshould be largest in low stellar mass systems, a result that agreeswith the negative slope of dust-to-stellar mass ratio versus stellarmass shown in Figures 6 and 7. These less massive galaxies havehigh sSFRs that allow the injection of dust from stellar sourcesmore quickly than in larger systems, and since their sSFRs peaklater than those of more massive galaxies (Cowie et al. 1996),smaller galaxies can see more growth in the dust-to-stellar massratio. Knowing a galaxy’s stellar mass, and its star-formation his-tory, allows us to better estimate its dust content. In this work, we extended the dust model in the moving-mesh code
AREPO to account for thermal sputtering of grains and performedcosmological hydrodynamical simulations to analyse the evolutionof dust in a diverse sample of galaxies. We studied the evolutionof the DMF, the radial distribution of dust in galactic haloes and onMpc scales, and the contribution of dust to the cosmic mass budget.Also, we explored how a galaxy’s SFR and stellar mass impacts itsdust content. Our main conclusions are as follows:(i) Our model broadly reproduces the observed z = 0 DMFover the range of masses accessible in a (25 h − Mpc ) volume.The DMF is presented for simulations at three resolutions, with thehighest-resolution simulation softening z = 0 gravitational forceson scales of h − pc.(ii) The mean dust surface density profile for simulated galaxieswith < i < at z = 0 . declines with radial distance, similarto the Σ dust ∝ r − . scaling seen in SDSS data out to projecteddistances of Mpc, although the normalisation of the simulateddust surface density lies up to dex below observations for r (cid:38) kpc.(iii) The cosmic dust density parameter at z = 0 is estimated tobe Ω dust = 1 . × − , close to values obtained from low-redshiftobservations. We see little evolution in Ω dust for z (cid:46) . , in tensionwith power spectrum-derived measurements that show a decline ofroughly . dex. This conflict is consistent with our model’s under-production of dusty galaxies for the high-redshift DMF.(iv) At both high and low redshift, dust mass increases with stel-lar mass along the star formation main sequence. This suggests thatsemi-analytic or galaxy formation models without dust trackingcan estimate dust content using the star formation main sequence.Semi-analytic models may also benefit from fitting functions forsubmillimetre number densities.(v) The dust-to-stellar mass ratio is predicted to anti-correlatewith stellar mass at high and low redshift, and this relation paral-lels observations at z = 0 . Less massive systems witness growthin dust-to-stellar mass ratio over < z < . Our simulated galax-ies also agree well with the observed distribution of dust-to-stellarmass ratio versus sSFR at z = 0 .(vi) By combining direct dust mass tracking with stellar popu-lation synthesis postprocessing, we predict dust-to-stellar mass andflux ratios for our simulated galaxies at z = 0 and compare to ob-servations. Coupling with empirical relations from radiative trans-fer simulations, we estimate the high-redshift submillimetre num-ber density functions for our sample of galaxies at µ m.(vii) While our model reproduces the observed z = 0 DMFfairly well, it is unable to capture the abundance of dust-rich galax-ies at high redshift. Instead, the simulated DMF evolves in a fairlymonotonic fashion. Adopting a top-heavy IMF does increase theabundance of high-redshift dusty galaxies but not to the extent seenin observations.(viii) To better match the observed DMF evolution, we mayneed physical prescriptions that are more closely connected to thebehaviour of the cosmic SFR density and produce more dustygalaxies near the peak of star formation. For example, adoptingnon-constant dust condensation efficiencies that vary with ISMdensity and temperature may allow the largest galaxies to more ef-ficiently produce dust at high redshift but limit dust formation atlower redshifts where star formation is less efficient and the DMFshifts towards lower masses.The dust model presented in this work yields low-redshift re- c (cid:13) ??? RAS, MNRAS000
AREPO to account for thermal sputtering of grains and performedcosmological hydrodynamical simulations to analyse the evolutionof dust in a diverse sample of galaxies. We studied the evolutionof the DMF, the radial distribution of dust in galactic haloes and onMpc scales, and the contribution of dust to the cosmic mass budget.Also, we explored how a galaxy’s SFR and stellar mass impacts itsdust content. Our main conclusions are as follows:(i) Our model broadly reproduces the observed z = 0 DMFover the range of masses accessible in a (25 h − Mpc ) volume.The DMF is presented for simulations at three resolutions, with thehighest-resolution simulation softening z = 0 gravitational forceson scales of h − pc.(ii) The mean dust surface density profile for simulated galaxieswith < i < at z = 0 . declines with radial distance, similarto the Σ dust ∝ r − . scaling seen in SDSS data out to projecteddistances of Mpc, although the normalisation of the simulateddust surface density lies up to dex below observations for r (cid:38) kpc.(iii) The cosmic dust density parameter at z = 0 is estimated tobe Ω dust = 1 . × − , close to values obtained from low-redshiftobservations. We see little evolution in Ω dust for z (cid:46) . , in tensionwith power spectrum-derived measurements that show a decline ofroughly . dex. This conflict is consistent with our model’s under-production of dusty galaxies for the high-redshift DMF.(iv) At both high and low redshift, dust mass increases with stel-lar mass along the star formation main sequence. This suggests thatsemi-analytic or galaxy formation models without dust trackingcan estimate dust content using the star formation main sequence.Semi-analytic models may also benefit from fitting functions forsubmillimetre number densities.(v) The dust-to-stellar mass ratio is predicted to anti-correlatewith stellar mass at high and low redshift, and this relation paral-lels observations at z = 0 . Less massive systems witness growthin dust-to-stellar mass ratio over < z < . Our simulated galax-ies also agree well with the observed distribution of dust-to-stellarmass ratio versus sSFR at z = 0 .(vi) By combining direct dust mass tracking with stellar popu-lation synthesis postprocessing, we predict dust-to-stellar mass andflux ratios for our simulated galaxies at z = 0 and compare to ob-servations. Coupling with empirical relations from radiative trans-fer simulations, we estimate the high-redshift submillimetre num-ber density functions for our sample of galaxies at µ m.(vii) While our model reproduces the observed z = 0 DMFfairly well, it is unable to capture the abundance of dust-rich galax-ies at high redshift. Instead, the simulated DMF evolves in a fairlymonotonic fashion. Adopting a top-heavy IMF does increase theabundance of high-redshift dusty galaxies but not to the extent seenin observations.(viii) To better match the observed DMF evolution, we mayneed physical prescriptions that are more closely connected to thebehaviour of the cosmic SFR density and produce more dustygalaxies near the peak of star formation. For example, adoptingnon-constant dust condensation efficiencies that vary with ISMdensity and temperature may allow the largest galaxies to more ef-ficiently produce dust at high redshift but limit dust formation atlower redshifts where star formation is less efficient and the DMFshifts towards lower masses.The dust model presented in this work yields low-redshift re- c (cid:13) ??? RAS, MNRAS000 , 1–17 imulating the dust content of galaxies sults in rough agreement with a number of observables across adiverse sample of galaxies, but it also highlights areas of tension.In particular, this model fails to predict the abundance of dust-richgalaxies at high redshift and the slight decline in Ω dust as galaxiesevolve towards low redshift. Furthermore, to truly probe the high-redshift submillimetre regime and the massive end of the DMFwill require larger cosmological volumes. Nonetheless, this workdemonstrates how simulations of large galaxy populations can beused to study the evolution of dust across diverse environments andthe distribution of dust on cosmological scales. ACKNOWLEDGEMENTS
We thank Volker Springel for providing us with access to
AREPO .We also thank the referee for their constructive feedback.The simulations were performed on the joint MIT-Harvardcomputing cluster supported by MKI and FAS. RM acknowl-edges support from the DOE CSGF under grant number DE-FG02-97ER25308. MV acknowledges support through an MIT RSCaward. CCH is grateful to the Gordon and Betty Moore Founda-tion for financial support.
REFERENCES
Abramson L. E., Kelson D. D., Dressler A., Poggianti B., Glad-ders M. D., Oemler, Jr. A., Vulcani B., 2014, ApJL, 785, L36Adelberger K. L., Steidel C. C., 2000, ApJ, 544, 218Asano R. S., Takeuchi T. T., Hirashita H., Inoue A. K., 2013a,Earth, Planets, and Space, 65, 213Asano R. S., Takeuchi T. T., Hirashita H., Nozawa T., 2013b, MN-RAS, 432, 637Barger A. J., Cowie L. L., Sanders D. B., Fulton E., Taniguchi Y.,Sato Y., Kawara K., Okuda H., 1998, Nature, 394, 248Barlow M. J., 1978, MNRAS, 183, 367Bauermeister A., Blitz L., Bolatto A., Bureau M., Teuben P.,Wong T., Wright M., 2013, ApJ, 763, 64Baugh C. M., Lacey C. G., Frenk C. S., Granato G. L., Silva L.,Bressan A., Benson A. J., Cole S., 2005, MNRAS, 356, 1191Bianchi S., Ferrara A., 2005, MNRAS, 358, 379Bianchi S., Schneider R., 2007, MNRAS, 378, 973Blain A. W., Jameson A., Smail I., Longair M. S., Kneib J.-P.,Ivison R. J., 1999a, MNRAS, 309, 715Blain A. W., Kneib J.-P., Ivison R. J., Smail I., 1999b, ApJL, 512,L87Boselli A., Cortese L., Boquien M., 2014, A&A, 564, A65Bothwell M. S. et al., 2014, MNRAS, 445, 2599Bourne N. et al., 2012, MNRAS, 421, 3027Brinchmann J., Charlot S., White S. D. M., Tremonti C., Kauff-mann G., Heckman T., Brinkmann J., 2004, MNRAS, 351, 1151Bruzual G., Charlot S., 2003, MNRAS, 344, 1000Burke J. R., Silk J., 1974, ApJ, 190, 1Casey C. M., Narayanan D., Cooray A., 2014, Phys. Rep., 541, 45Chabrier G., 2003, PASP, 115, 763Chakrabarti S., Fenner Y., Cox T. J., Hernquist L., Whitney B. A.,2008, ApJ, 688, 972Chakrabarti S., Whitney B. A., 2009, ApJ, 690, 1432Chapman S. C., Blain A. W., Smail I., Ivison R. J., 2005, ApJ,622, 772Chary R., Elbaz D., 2001, ApJ, 556, 562 Chen C.-C., Cowie L. L., Barger A. J., Casey C. M., Lee N.,Sanders D. B., Wang W.-H., Williams J. P., 2013, ApJ, 762, 81Clark C. J. R. et al., 2015, MNRAS, 452, 397Clemens M. S. et al., 2013, MNRAS, 433, 695Conroy C., Gunn J. E., 2010, ApJ, 712, 833Conroy C., Gunn J. E., White M., 2009, ApJ, 699, 486Cortese L. et al., 2012, A&A, 540, A52Cowie L. L., Barger A. J., Kneib J.-P., 2002, AJ, 123, 2197Cowie L. L., Songaila A., Hu E. M., Cohen J. G., 1996, AJ, 112,839da Cunha E., Eminian C., Charlot S., Blaizot J., 2010, MNRAS,403, 1894Daddi E. et al., 2010, ApJ, 713, 686Dav´e R., Finlator K., Oppenheimer B. D., Fardal M., Katz N.,Kereˇs D., Weinberg D. H., 2010, MNRAS, 404, 1355Davies J. I. et al., 2012, MNRAS, 419, 3505De Bernardis F., Cooray A., 2012, ApJ, 760, 14Dey A. et al., 2008, ApJ, 677, 943Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS,399, 497Draine B. T. et al., 2007, ApJ, 663, 866Draine B. T., Salpeter E. E., 1979a, ApJ, 231, 438Draine B. T., Salpeter E. E., 1979b, ApJ, 231, 77Driver S. P., Popescu C. C., Tuffs R. J., Liske J., Graham A. W.,Allen P. D., de Propris R., 2007, MNRAS, 379, 1022Duffy A. R., Kay S. T., Battye R. A., Booth C. M., Dalla VecchiaC., Schaye J., 2012, MNRAS, 420, 2799Dunne L., Eales S., Edmunds M., Ivison R., Alexander P.,Clements D. L., 2000, MNRAS, 315, 115Dunne L., Eales S. A., 2001, MNRAS, 327, 697Dunne L., Eales S. A., Edmunds M. G., 2003, MNRAS, 341, 589Dunne L. et al., 2011, MNRAS, 417, 1510Dwek E., Arendt R. G., 1992, ARA&A, 30, 11Dwek E., Foster S. M., Vancura O., 1996, ApJ, 457, 244Eales S. et al., 2009, ApJ, 707, 1779Eales S. et al., 2010, PASP, 122, 499Eales S., Lilly S., Gear W., Dunne L., Bond J. R., Hammer F., LeF`evre O., Crampton D., 1999, ApJ, 515, 518Edmunds M. G., 2001, MNRAS, 328, 223Fisher D. B., Bolatto A., Drory N., Combes F., Blitz L., Wong T.,2013, ApJ, 764, 174Fu J., Kauffmann G., Li C., Guo Q., 2012, MNRAS, 424, 2701Fukugita M., 2011, arXiv:1103.4191Fukugita M., Peebles P. J. E., 2004, ApJ, 616, 643Galametz M., Madden S. C., Galliano F., Hony S., Bendo G. J.,Sauvage M., 2011, A&A, 532, A56Geach J. E., Smail I., Moran S. M., MacArthur L. A., LagosC. d. P., Edge A. C., 2011, ApJL, 730, L19Genel S. et al., 2014, MNRAS, 445, 175Groenewegen M. A. T., 1997, A&A, 317, 503Hahn O., Abel T., 2011, MNRAS, 415, 2101Hayward C. C., Jonsson P., Kereˇs D., Magnelli B., Hernquist L.,Cox T. J., 2012, MNRAS, 424, 951Hayward C. C., Kereˇs D., Jonsson P., Narayanan D., Cox T. J.,Hernquist L., 2011, ApJ, 743, 159Hayward C. C., Narayanan D., Kereˇs D., Jonsson P., Hopkins P. F.,Cox T. J., Hernquist L., 2013, MNRAS, 428, 2529Hezaveh Y. D. et al., 2013, ApJ, 767, 132Hirashita H., 2000, PASJ, 52, 585Hirashita H., Nozawa T., Villaume A., Srinivasan S., 2015, MN-RAS, 454, 1620Holland W. S. et al., 1999, MNRAS, 303, 659 c (cid:13) ??? RAS, MNRAS , 1–17 R. McKinnon et al.
Hopkins P. F. et al., 2009, MNRAS, 397, 802Hughes D. H. et al., 1998, Nature, 394, 241Itoh H., 1989, PASJ, 41, 853Jones A. P., Tielens A. G. G. M., Hollenbach D. J., McKee C. F.,1994, ApJ, 433, 797Jonsson P., 2006, MNRAS, 372, 2Jonsson P., Cox T. J., Primack J. R., Somerville R. S., 2006, ApJ,637, 255Karim A. et al., 2011, ApJ, 730, 61Kennicutt, Jr. R. C. et al., 2009, ApJ, 703, 1672Knebe A. et al., 2015, MNRAS, 451, 4029Knudsen K. K., van der Werf P. P., Kneib J.-P., 2008, MNRAS,384, 1611Lagache G., Puget J.-L., Dole H., 2005, ARA&A, 43, 727Lagos C. D. P., Baugh C. M., Lacey C. G., Benson A. J., KimH.-S., Power C., 2011a, MNRAS, 418, 1649Lagos C. d. P. et al., 2015, MNRAS, 452, 3815Lagos C. D. P., Lacey C. G., Baugh C. M., Bower R. G., BensonA. J., 2011b, MNRAS, 416, 1566Leroy A., Bolatto A., Stanimirovic S., Mizuno N., Israel F., BotC., 2007, ApJ, 658, 1027Leroy A. K. et al., 2011, ApJ, 737, 12Leroy A. K., Walter F., Brinks E., Bigiel F., de Blok W. J. G.,Madore B., Thornley M. D., 2008, AJ, 136, 2782Lo Faro B. et al., 2013, ApJ, 762, 108Madau P., Dickinson M., 2014, ARA&A, 52, 415Magdis G. E. et al., 2012, ApJ, 760, 6Magnelli B. et al., 2012, A&A, 539, A155Masaki S., Yoshida N., 2012, MNRAS, 423, L117Mattsson L., Andersen A. C., 2012, MNRAS, 423, 38Mattsson L., Andersen A. C., Munkhammar J. D., 2012, MNRAS,423, 26McGee S. L., Balogh M. L., 2010, MNRAS, 405, 2069McKinnon R., Torrey P., Vogelsberger M., 2016, MNRAS, 457,3775M´enard B., Fukugita M., 2012, ApJ, 754, 116M´enard B., Scranton R., Fukugita M., Richards G., 2010, MN-RAS, 405, 1025Michałowski M. J., Watson D., Hjorth J., 2010, ApJ, 712, 942Morokuma-Matsui K., Baba J., 2015, MNRAS, 454, 3792Mu˜noz-Mateos J. C. et al., 2009, ApJ, 701, 1965Narayanan D., Cox T. J., Hayward C. C., Younger J. D., HernquistL., 2009, MNRAS, 400, 1919Narayanan D. et al., 2010, MNRAS, 407, 1701Narayanan D. et al., 2015, Nature, 525, 496Nozawa T., Kozasa T., Habe A., 2006, ApJ, 648, 435Nozawa T., Kozasa T., Habe A., Dwek E., Umeda H., TominagaN., Maeda K., Nomoto K., 2007, ApJ, 666, 955Obreschkow D., Croton D., De Lucia G., Khochfar S., RawlingsS., 2009, ApJ, 698, 1467Ostriker J., Silk J., 1973, ApJL, 184, L113Ouchi M. et al., 2013, ApJ, 778, 102Planck Collaboration et al., 2014, A&A, 571, A16Popping G., Caputi K. I., Somerville R. S., Trager S. C., 2012,MNRAS, 425, 2386Popping G., Somerville R. S., Galametz M., 2016,arXiv:1609.08622Popping G., Somerville R. S., Trager S. C., 2014, MNRAS, 442,2398Reddy N. A., Erb D. K., Pettini M., Steidel C. C., Shapley A. E.,2010, ApJ, 712, 1070 Reddy N. A., Steidel C. C., Fadda D., Yan L., Pettini M., ShapleyA. E., Erb D. K., Adelberger K. L., 2006, ApJ, 644, 792R´emy-Ruyer A. et al., 2014, A&A, 563, A31R´emy-Ruyer A. et al., 2015, A&A, 582, A121Riechers D. A. et al., 2013, Nature, 496, 329Rowlands K. et al., 2012, MNRAS, 419, 2545Saintonge A. et al., 2011, MNRAS, 415, 32Salim S. et al., 2007, ApJS, 173, 267Salpeter E. E., 1977, ARA&A, 15, 267Santini P. et al., 2014, A&A, 562, A30Schechter P., 1976, ApJ, 203, 297Scott S. E. et al., 2002, MNRAS, 331, 817Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MN-RAS, 380, 877Silva L., Granato G. L., Bressan A., Danese L., 1998, ApJ, 509,103Skibba R. A. et al., 2011, ApJ, 738, 89Sklias P. et al., 2014, A&A, 561, A149Smail I., Ivison R. J., Blain A. W., 1997, ApJL, 490, L5Smail I., Ivison R. J., Blain A. W., Kneib J.-P., 2002, MNRAS,331, 495Smith M. W. L. et al., 2012, ApJ, 748, 123Smith R. K., Krzewina L. G., Cox D. P., Edgar R. J., Miller, IIIW. W., 1996, ApJ, 473, 864Sparre M. et al., 2015, MNRAS, 447, 3548Sparre M., Springel V., 2016, MNRAS, 462, 2418Spergel D. N., Flauger R., Hloˇzek R., 2015, Phys. Rev. D, 91,023518Springel V., 2010, MNRAS, 401, 791Springel V., Hernquist L., 2003, MNRAS, 339, 289Springel V., White S. D. M., Tormen G., Kauffmann G., 2001,MNRAS, 328, 726Swinbank A. M. et al., 2008, MNRAS, 391, 420Tacconi L. J. et al., 2008, ApJ, 680, 246Tacconi L. J. et al., 2013, ApJ, 768, 74Thacker C. et al., 2013, ApJ, 768, 58Tielens A. G. G. M., McKee C. F., Seab C. G., Hollenbach D. J.,1994, ApJ, 431, 321Torrey P., Vogelsberger M., Genel S., Sijacki D., Springel V.,Hernquist L., 2014, MNRAS, 438, 1985Tsai J. C., Mathews W. G., 1995, ApJ, 448, 84Vieira J. D. et al., 2013, Nature, 495, 344Vlahakis C., Dunne L., Eales S., 2005, MNRAS, 364, 1253Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V.,Hernquist L., 2013, MNRAS, 436, 3031Vogelsberger M. et al., 2014a, Nature, 509, 177Vogelsberger M. et al., 2014b, MNRAS, 444, 1518Vogelsberger M., Sijacki D., Kereˇs D., Springel V., Hernquist L.,2012, MNRAS, 425, 3024Wang B., Heckman T. M., 1996, ApJ, 457, 645Watson D., Christensen L., Knudsen K. K., Richard J., GallazziA., Michałowski M. J., 2015, Nature, 519, 327Weiß A. et al., 2013, ApJ, 767, 88Whitaker K. E., van Dokkum P. G., Brammer G., Franx M., 2012,ApJL, 754, L29Wiersma R. P. C., Schaye J., Smith B. D., 2009, MNRAS, 393, 99Wilson G. W. et al., 2008, MNRAS, 386, 807Winters J. M., Fleischer A. J., Le Bertre T., Sedlmayr E., 1997,A&A, 326, 305Yahil A., Ostriker J. P., 1973, ApJ, 185, 787Yajima H., Li Y., Zhu Q., Abel T., 2012, MNRAS, 424, 884Yasuda Y., Kozasa T., 2012, ApJ, 745, 159 c (cid:13) ??? RAS, MNRAS000
Hopkins P. F. et al., 2009, MNRAS, 397, 802Hughes D. H. et al., 1998, Nature, 394, 241Itoh H., 1989, PASJ, 41, 853Jones A. P., Tielens A. G. G. M., Hollenbach D. J., McKee C. F.,1994, ApJ, 433, 797Jonsson P., 2006, MNRAS, 372, 2Jonsson P., Cox T. J., Primack J. R., Somerville R. S., 2006, ApJ,637, 255Karim A. et al., 2011, ApJ, 730, 61Kennicutt, Jr. R. C. et al., 2009, ApJ, 703, 1672Knebe A. et al., 2015, MNRAS, 451, 4029Knudsen K. K., van der Werf P. P., Kneib J.-P., 2008, MNRAS,384, 1611Lagache G., Puget J.-L., Dole H., 2005, ARA&A, 43, 727Lagos C. D. P., Baugh C. M., Lacey C. G., Benson A. J., KimH.-S., Power C., 2011a, MNRAS, 418, 1649Lagos C. d. P. et al., 2015, MNRAS, 452, 3815Lagos C. D. P., Lacey C. G., Baugh C. M., Bower R. G., BensonA. J., 2011b, MNRAS, 416, 1566Leroy A., Bolatto A., Stanimirovic S., Mizuno N., Israel F., BotC., 2007, ApJ, 658, 1027Leroy A. K. et al., 2011, ApJ, 737, 12Leroy A. K., Walter F., Brinks E., Bigiel F., de Blok W. J. G.,Madore B., Thornley M. D., 2008, AJ, 136, 2782Lo Faro B. et al., 2013, ApJ, 762, 108Madau P., Dickinson M., 2014, ARA&A, 52, 415Magdis G. E. et al., 2012, ApJ, 760, 6Magnelli B. et al., 2012, A&A, 539, A155Masaki S., Yoshida N., 2012, MNRAS, 423, L117Mattsson L., Andersen A. C., 2012, MNRAS, 423, 38Mattsson L., Andersen A. C., Munkhammar J. D., 2012, MNRAS,423, 26McGee S. L., Balogh M. L., 2010, MNRAS, 405, 2069McKinnon R., Torrey P., Vogelsberger M., 2016, MNRAS, 457,3775M´enard B., Fukugita M., 2012, ApJ, 754, 116M´enard B., Scranton R., Fukugita M., Richards G., 2010, MN-RAS, 405, 1025Michałowski M. J., Watson D., Hjorth J., 2010, ApJ, 712, 942Morokuma-Matsui K., Baba J., 2015, MNRAS, 454, 3792Mu˜noz-Mateos J. C. et al., 2009, ApJ, 701, 1965Narayanan D., Cox T. J., Hayward C. C., Younger J. D., HernquistL., 2009, MNRAS, 400, 1919Narayanan D. et al., 2010, MNRAS, 407, 1701Narayanan D. et al., 2015, Nature, 525, 496Nozawa T., Kozasa T., Habe A., 2006, ApJ, 648, 435Nozawa T., Kozasa T., Habe A., Dwek E., Umeda H., TominagaN., Maeda K., Nomoto K., 2007, ApJ, 666, 955Obreschkow D., Croton D., De Lucia G., Khochfar S., RawlingsS., 2009, ApJ, 698, 1467Ostriker J., Silk J., 1973, ApJL, 184, L113Ouchi M. et al., 2013, ApJ, 778, 102Planck Collaboration et al., 2014, A&A, 571, A16Popping G., Caputi K. I., Somerville R. S., Trager S. C., 2012,MNRAS, 425, 2386Popping G., Somerville R. S., Galametz M., 2016,arXiv:1609.08622Popping G., Somerville R. S., Trager S. C., 2014, MNRAS, 442,2398Reddy N. A., Erb D. K., Pettini M., Steidel C. C., Shapley A. E.,2010, ApJ, 712, 1070 Reddy N. A., Steidel C. C., Fadda D., Yan L., Pettini M., ShapleyA. E., Erb D. K., Adelberger K. L., 2006, ApJ, 644, 792R´emy-Ruyer A. et al., 2014, A&A, 563, A31R´emy-Ruyer A. et al., 2015, A&A, 582, A121Riechers D. A. et al., 2013, Nature, 496, 329Rowlands K. et al., 2012, MNRAS, 419, 2545Saintonge A. et al., 2011, MNRAS, 415, 32Salim S. et al., 2007, ApJS, 173, 267Salpeter E. E., 1977, ARA&A, 15, 267Santini P. et al., 2014, A&A, 562, A30Schechter P., 1976, ApJ, 203, 297Scott S. E. et al., 2002, MNRAS, 331, 817Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MN-RAS, 380, 877Silva L., Granato G. L., Bressan A., Danese L., 1998, ApJ, 509,103Skibba R. A. et al., 2011, ApJ, 738, 89Sklias P. et al., 2014, A&A, 561, A149Smail I., Ivison R. J., Blain A. W., 1997, ApJL, 490, L5Smail I., Ivison R. J., Blain A. W., Kneib J.-P., 2002, MNRAS,331, 495Smith M. W. L. et al., 2012, ApJ, 748, 123Smith R. K., Krzewina L. G., Cox D. P., Edgar R. J., Miller, IIIW. W., 1996, ApJ, 473, 864Sparre M. et al., 2015, MNRAS, 447, 3548Sparre M., Springel V., 2016, MNRAS, 462, 2418Spergel D. N., Flauger R., Hloˇzek R., 2015, Phys. Rev. D, 91,023518Springel V., 2010, MNRAS, 401, 791Springel V., Hernquist L., 2003, MNRAS, 339, 289Springel V., White S. D. M., Tormen G., Kauffmann G., 2001,MNRAS, 328, 726Swinbank A. M. et al., 2008, MNRAS, 391, 420Tacconi L. J. et al., 2008, ApJ, 680, 246Tacconi L. J. et al., 2013, ApJ, 768, 74Thacker C. et al., 2013, ApJ, 768, 58Tielens A. G. G. M., McKee C. F., Seab C. G., Hollenbach D. J.,1994, ApJ, 431, 321Torrey P., Vogelsberger M., Genel S., Sijacki D., Springel V.,Hernquist L., 2014, MNRAS, 438, 1985Tsai J. C., Mathews W. G., 1995, ApJ, 448, 84Vieira J. D. et al., 2013, Nature, 495, 344Vlahakis C., Dunne L., Eales S., 2005, MNRAS, 364, 1253Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V.,Hernquist L., 2013, MNRAS, 436, 3031Vogelsberger M. et al., 2014a, Nature, 509, 177Vogelsberger M. et al., 2014b, MNRAS, 444, 1518Vogelsberger M., Sijacki D., Kereˇs D., Springel V., Hernquist L.,2012, MNRAS, 425, 3024Wang B., Heckman T. M., 1996, ApJ, 457, 645Watson D., Christensen L., Knudsen K. K., Richard J., GallazziA., Michałowski M. J., 2015, Nature, 519, 327Weiß A. et al., 2013, ApJ, 767, 88Whitaker K. E., van Dokkum P. G., Brammer G., Franx M., 2012,ApJL, 754, L29Wiersma R. P. C., Schaye J., Smith B. D., 2009, MNRAS, 393, 99Wilson G. W. et al., 2008, MNRAS, 386, 807Winters J. M., Fleischer A. J., Le Bertre T., Sedlmayr E., 1997,A&A, 326, 305Yahil A., Ostriker J. P., 1973, ApJ, 185, 787Yajima H., Li Y., Zhu Q., Abel T., 2012, MNRAS, 424, 884Yasuda Y., Kozasa T., 2012, ApJ, 745, 159 c (cid:13) ??? RAS, MNRAS000 , 1–17 imulating the dust content of galaxies Zhukovska S., 2014, A&A, 562, A76Zu Y., Weinberg D. H., Dav´e R., Fardal M., Katz N., Kereˇs D.,Oppenheimer B. D., 2011, MNRAS, 412, 1059
APPENDIX A: ANALYTICAL THERMAL SPUTTERINGCALCULATIONS
The empirical thermal sputtering rate given in Equation (1) fallsoff quickly for T (cid:46) K. In this section, we detail how such atemperature dependence arises from analytical calculations of col-lisions between gas atoms and grains and the subsequent erosion ofgrains from sputtering. We refer the reader to the existing literature(Barlow 1978; Draine & Salpeter 1979b; Tielens et al. 1994) formore thorough analysis.Following Equations 4.19 and 4.20 in Tielens et al. (1994),consider a grain of radius a in a medium of temperature T . Then,the number of particles sputtered off of the grain surface per unittime is given by d N sp d t = 2 πa (cid:88) i n i (cid:104) Y i v (cid:105) , (A1)where the leading factor of two accounts for collisions at non-normal angles, πa is the grain cross-section, and n i and (cid:104) Y i v (cid:105) arethe number density and Maxwell-Boltzmann distribution-averagedproduct of sputtering yield and velocity for a gas ion of species i .Here, Y i measures the number of particles sputtered from the grainper gas ion collision (e.g. studied in detail in Section 4.1 of Tielenset al. 1994), and the Maxwell-Boltzmann distribution correspondsto temperature T .Suppose the grain has mass m and uniform internal density ρ g and that particles sputtered from the grain surface have mass m sp . For example, a carbonaceous grain might have m sp = m C ,the mass of a carbon atom. The grain mass loss rate d m d t = 4 πa ρ g d a d t (A2)implies that the change in grain radius per unit time due to thermalsputtering is given by d a d t = n H m sp ρ g (cid:88) i A i (cid:104) Y i v (cid:105) , (A3)where A i is the abundance of gas ions of species i . This sputter-ing rate is a function of temperature due to its averaging over aMaxwell-Boltzmann distribution. Combined with analytic modelsof sputtering yields, the thermal sputtering rate shows a sharp drop-off for T (cid:46) K. Thus, the empirical formula given by Equa-tion (1) captures the essential temperature dependence and normal-isation of the thermal sputtering rate and avoids the need to calcu-late Maxwell-Boltzmann distribution-averaged sputtering integralsin our simulation code.
APPENDIX B: VARIATION OF GRAIN SIZEPARAMETER
The fiducial L25n256 run assumed a grain size of a = 0 . µ mto estimate thermal sputtering rates in Equation (2). To investi-gate the sensitivity of our results to this choice, we performed twoadditional runs at the same resolution level: one with a smallergrain size ( a = 0 . µ m) and another with a larger grain size( a = 1 µ m). Figure B1 shows the cosmic dust density and dust z ρ du s t [ M (cid:12) M p c − ] L25n256small grainlarge grainFukugita+ (2004)Driver+ (2007)Menard+ (2010)Dunne+ (2011)Fukugita (2011)De Bernardis+ (2012)Menard+ (2012)Clemens+ (2013)Thacker+ (2013) − − − − Ω du s t Figure B1.
Same as Figure 4, except using the medium-resolution sim-ulation with three grain size parameters. The fiducial L25n256 run uses a = 0 . µ m, while the small and large grain runs use a = 0 . µ m and a = 1 µ m, respectively. Evolution in the cosmic dust density is not sensi-tive to the choice of grain size parameter. density parameter for < z < , which was previously studied inFigure 4.First, the results yield the correct qualitative behaviour: thesputtering time-scale estimated in Equation (2) is longer for largergrains, and we see that by z = 0 the cosmic dust density is largestfor the a = 1 µ m run and smallest for the a = 0 . µ m run. How-ever, the dust densities predicted by these three runs differ by lessthan a factor of two. Figure 4 demonstrated that change in cos-mic dust density when improving resolution from the L25n128 toL25n256 run was just as large as varying the grain size parameterby a factor of . Also, the observational data shown for compar-ison indicate that there are larger uncertainties when estimating thecosmic dust density through a variety of means (e.g. quasar-galaxyreddening correlations, power spectrum measurements, DMF inte-gration, etc.).Thus, the results presented in Section 3 are not sensitive to ourchoice of a , especially when considering the combined uncertain-ties in dust condensation efficiencies, dust mass absorption coeffi-cients, and the amount of dust in galactic haloes. APPENDIX C: VARIATION OF SIMULATION VOLUME
In addition to the fiducial runs, we also simulate a (50 h − Mpc ) volume down to z = 2 . with × dark matter and gas par-ticles to start. This run, labelled L50n512, uses the same fiducialparameters from Table 1 and offers the same spatial and mass res-olution as the L25n256 run, but in a volume eight times as large.Figure C1 shows the DMFs for the L25n256 and L50n512 runs at z = 2 . , which are nearly identical. The number of galaxies in themass bin covering . M (cid:12) (cid:54) M dust < . M (cid:12) has increasedfrom in the L25n256 run to in the L50n512 run, a changesimilar to the factor of eight increase in volume between these runs.The L50n512 run also forms galaxies in the next highest massbin, which had no galaxies in the L25n256 run. Figure C1 suggeststhat the normalisation offset between the simulated DMFs and ob-servations is not the result of limited statistics. c (cid:13) ??? RAS, MNRAS , 1–17 R. McKinnon et al. M dust [M (cid:12) ]10 − − − − − − Φ [ M p c − d e x − ] L25n256L50n512Dunne+ (2003)
Figure C1.
Comparison of the z = 2 . DMFs for the L25n256 run(red) and the L50n512 run (green), the latter of which simulates a (50 h − Mpc ) volume with the same resolution as the L25n256 run. Thisvolume is eight times larger than the fiducial volume and enables us tosample more galaxies. The normalisation of the DMF is not substantiallychanged by an increase in simulated volume. c (cid:13) ??? RAS, MNRAS000