Draftversion April15,2015 PreprinttypesetusingLATEXstyleemulateapjv.11/10/09 SIMULATING DEEP HUBBLE IMAGES WITH SEMI-EMPIRICAL MODELS OF GALAXY FORMATION Manuchehr Taghizadeh-Popp1,2, S. Michael Fall1, Richard L. White1, and Alexander S. Szalay2 Draft version April 15, 2015 ABSTRACT We simulate deep images from the Hubble Space Telescope (HST) using semi-empirical models of 5 galaxy formation with only a few basic assumptions and parameters. We project our simulations all 1 thewaytothe observationaldomain,addingcosmologicalandinstrumentaleffects totheimages,and 0 analyze them in the same way as real HST images (“forward modeling”). This is a powerful tool 2 for testing and comparing galaxy evolution models, since it allows us to make unbiased comparisons r betweenthepredictedandobserveddistributionsofgalaxyproperties,whileautomaticallytakinginto p account all relevant selection effects. A Oursemi-empiricalmodelspopulateeachdarkmatterhalowithagalaxyofdeterminedstellarmass 3 andscaleradius. Wecomputetheluminosityandspectrumofeachsimulatedgalaxyfromitsevolving 1 stellarmass using stellar populationsynthesismodels. We calculate the intrinsic scatterin the stellar mass-halo mass relation that naturally results from enforcing a monotonically increasing stellar mass ] along the merger history of each halo. The simulated galaxy images are drawn from cutouts of real A galaxiesfromtheSloanDigitalSkySurvey,withsizesandfluxesrescaledtomatchthoseofthe model G galaxies. The distributions of galaxy luminosities, sizes, and surface brightnesses depend on the adjustable . h parameters in the models, and they agree well with observations for reasonable values of those pa- p rameters. Measured galaxy magnitudes and sizes have significant magnitude-dependent biases, with - both being underestimatednear the magnitude detection limit. The fractionofgalaxiesdetected and o fraction of light detected also depend sensitively on the details of the model. r t Subject headings: galaxies: evolution — galaxies: formation — galaxies: fundamental parameters — s a galaxies: general — galaxies: statistics — large-scale structure of universe [ 2 1. INTRODUCTION analyze these images in the same way as real HST im- v ages, extracting catalogs from the simulated images to The Hubble Space Telescope (HST) has spent much 6 detect and measure the fluxes, sizes, and other proper- of its operational lifetime staring into deep space, sur- 5 ties of the galaxies. Our simulated images include both 5 veying galaxies in their infancy and youth. The moti- cosmologicaleffects (projectionalong pencil beams, red- 1 vation for these deep surveys (the Hubble Deep Field shiftingofpassbands,dimmingoffluxandsurfacebright- 0 and its successors)was to obtain time-lapse images that ness) and instrumental effects (point-spread function . would reveal how galaxies formed and evolved. While 1 [PSF], pixelation, noise, sky backgrounds) for any given we have made great progress in interpreting the deep 0 HST camera,filter, andexposuretime. We then extract HST surveys, this program remains challenging and far 5 catalogs of objects in the images with the widely used from complete (Galametz et al. 2013; Guo et al. 2013b; 1 SExtractor software (Bertin & Arnouts 1996). Thus, Bouwens et al. 2014, and references therein). There are : our procedure automatically takes into account all rele- v two major reasons for this. First, the samples of high- i redshift galaxies have been severely edited by selection vantselectioneffects,allowingustomakeunbiasedcom- X parisons between the predicted and observed distribu- effects, primarily limits on flux and surface brightness, r effectivelybiasingtheobservableuniversetowardbright, tionsofluminosities,sizes,andotherpropertiesofgalax- a ies. This is the approach recommended by textbooks compact galaxies. Second, on the theoretical side, there on statistical inference: map the predictions all the way are still many significant gaps in our understanding of intotheobservationaldomainandmakethecomparisons the physical processes that affect the baryonic compo- there, often called “forward modeling.” nents ofgalaxies(stars,gas,anddust) andthe radiation Theforwardmodelingmethodinnotsensitivetosmall they emit. These uncertainties are reflected in the many errors in the luminosities, sizes, and other properties of free parameters ofthe semi-analyticalmodels and in the galaxies, or even to the exact definitions of these quan- analogous sub-grid physics of the hydrodynamical mod- tities, because such errors affect measurements of both els. In this paper, we present a new approach to the the simulated and realimages in the same way. In other analysis and interpretation of deep galaxy surveys that words, these errors cancel out of the comparisons of the addresses both of these issues. simulated and real distributions of luminosity, size, and To account for selection effects, we create simulated so forth. This is one of the main advantages of the for- HST images of model galaxy populations, and we then ward modeling approach. 1Space Telescope Science Institute, 3700 San Martin Drive, In contrast, nearly all work in this field is based on BaltimoreMD21218,USA;[email protected] theopposite,butsimpler,approachofcomparingpredic- 2DepartmentofPhysicsandAstronomy,JohnsHopkinsUni- tionswithobservationsinthetheoreticaldomain(“back- versity,3400NorthCharlesStreet,Baltimore,MD21218, USA 2 Taghizadeh-Poppet al. ward modeling”). Some exceptions are the reconstruc- tions of galaxy properties (such as luminosity or size). tionofmockimagesordatastartingfromsemi-analytical Our approach can be used to inform the design of fu- models, (e.g. Blaizot et al. 2005; Overzier et al. 2013) ture surveys (e.g., choice of filters and exposure times) or hydrodynamical simulations (e.g Lotz et al. 2008; by addressing directly the question of which data are Devriendt2010;Mozena2013),althoughthelattersuffer mostuseful to discriminate between different theoretical from unrealistic star formation histories and too rapid models. This will be especially valuable in planning the growth of stellar masses at early times (Bouche et al. deepest surveys with the James Webb Space Telescope 2010). In this paper, we go beyond the creation of (JWST). mock galaxy images by deriving simulated distributions We emphasize that the main purpose of this paper of galaxy properties and comparing them with the cor- is to demonstrate the validity and utility of a general responding observated distributions. method for analyzing and interpreting deep galaxy sur- To limit the number of assumptions and parame- veys. This is an initial, exploratory study. We regard ters in our models, we adopt a semi-empirical descrip- our specific implementation of the method and the first tion of galaxy evolution. This description is based on results obtained from it and presented here as being il- the evolution of dark matter halos in cosmological N- lustrative rather than definitive. There is scope for fur- body simulations, which is now well understood, in con- ther development of the method and refinement of the trast to the evolution of the baryonic parts of galax- results. Nevertheless, the overall agreement we find be- ies. The main assumption of the semi-empirical descrip- tween our simulations and observations—with no fine- tion is that most of the information needed for simu- tuning of parameters—is remarkable and encouraging. lating a population of galaxies is already encoded in This paper is arranged as follows. §2 explains the the merger trees of their dark halos. Each halo is as- method for modeling and building simulated universes, sumed to host one model galaxy, and the properties of fromtheselectionofsemi-empiricalmodelsandthedark that galaxy are then uniquely determined by those of matter simulation to the building of simulated HST im- its halo, including its mass and size. The advantage of ages. The main results, presented in §3, include the im- this method is that it sidesteps much of the complex plementation of stellar mass-halo mass relations (§3.1), and uncertain baryonic processes in galaxy formation; present-day mass and luminosity distributions derived the disadvantage is that it likely oversimplifies some as- from the semi-empirical models of a simulated universe pects ofthese processes. This semi-empiricaldescription (§3.2), and a comparisonof the luminosity, size and sur- has been developed in numerous studies over the past face brightness distributions extracted from the simu- decade (e.g., Vale & Ostriker 2004; Conroy et al. 2006; lated images with measurements from real HST images Conroy & Wechsler 2009; Guo et al. 2010; Moster et al. (§3.3). We devote §4 to an analysis of the cosmic star 2010; Behroozi et al. 2010, 2013; Guo & White 2014; formation rate density. Lastly, §5 discusses the results Moster et al. 2013; Kravtsov2013; Reddick et al. 2013). and summarizes our conclusions. We derive the radiative spectrum of each simulated galaxy from its star formation history using stellar pop- 2. BUILDINGASIMULATEDUNIVERSE ulation synthesis models. The star formation history in In this section we describe in detail the steps required turn follows from the growth of its halo mass, including to build a self-consistent simulated universe and asso- both smooth accretion and discrete mergers with other ciated HST images. We also describe the parameters halos. Theradiativespectrumalsodependsonthemetal- chosen for our reference model for the universe, as well licity ofthe stellar populationandthe absorptionbygas as variations on those parameters explored in the other and dust in the galaxy and by gas in the intergalactic models. medium. In our implementation of this method, we use cutout 2.1. The Dark Matter Simulation images of real galaxies in the Sloan Digital Sky Survey Following standard practice, we use a ΛCDM simula- (SDSS) as templates for the visual appearance of our tion as the three-dimensional (3-D) skeleton of our sim- modelgalaxies,withtheirfluxesandsizesrescaledtore- ulated universes, placing a model galaxy at the center flectgalaxyevolutionaccordingtoourmodels. Thisgen- of each dark (sub)halo1. That defines the spatial distri- eratesmuchmorerealisticmorphologiesthanthesmooth bution and number density of galaxies as a function of S´ersic bulge+disk light profiles commonly used in pre- redshift. Most of the information needed in our method vious semi-empirical or semi-analytical models and im- is in fact provided by the halo mass and size evolution proves modeling of the detection incompleteness effects along halo merger trees, which does not involve any fit- that arise from the internal clumpiness of real galaxies. ting of free parameters (see §2.2 and §2.4). Our approachgivesus a valuable science toolfor com- For the merger trees, we use the milli-Millennium cos- paringmodelsofgalaxyevolution. Inthispaperwebuild mologicaldarkmattersimulation(mMS)(Springel et al. a reference model using plausible choices of parameters, 2005; Lemson & Springel 2006). Its smaller size com- and explore other models by changing one parameter at pared to the full Millennium or Millennium II simula- a time. Thus, we can test the sensitivity of the simu- tions (Springel et al. 2005; Boylan-Kolchinet al. 2009) lated universe to each of these parameters by compar- makesiteasiertouseforthis exploratorywork;weshow ing the statistics derivedfrom their respective simulated below that the coarser mass resolution in the mMS is images to those from reference model and real HST im- adequate for our simulations (§3.3.1). Halos are defined ages. Moreover,sincetheinput propertiesofeachmodel usingfriends-of-friendsgroupsasexplainedinGuo et al. galaxy on the simulated images are known, we are able toquantifythemodel-dependentoutput galaxydetection 1 We use “halo” and “sub-halo” interchangeably, following efficiency by comparing the input and output distribu- Guoetal.(2010). Simulating Hubble Images With Galaxy Formation Models 3 (2010). Bound dark matter structures or (sub)halos are composedofthemostmassivemain(orcentral)sub-halo surrounded by satellite sub-halos. 1 1 The mMS has the same cosmology and particle mass (1.18×109M⊙) as the much larger Millennium simula- tion, but with both a smaller box size (85.62 Mpc) and )0 ]1 a reduced number of particles (2703). The cosmologi- un s cal parameters used are the ones obtained by WMAP1 M (Spergel et al. 2003), i.e., ΩM = 0.25, ΩΛ = 0.75, h0 = M [s9 Guo 0tWh.7eM3dAaiPnffde7reσ(nK8co=emb0ae.t9tsw.ueAeetnsaetlx.hp2el0aW1in1M)edcAobPsym1Goaluonogdieestthdaeol.es(st2an0no1dt3aaar)fd-, log(8 BM es h = r o’0’ o. 0 z 2 i (5(zz =M=02h))alo fect significantly the relevant aspects of dark matter ’’ (z=4) structure. In fact, a smaller WMAP7 σ8 = 0.807 is 7 ’’ (z=6) counterbalancedbyagreaterΩM =0.272,whichresults, ’’ (z=8) for example, in the WMAP1 and WMAP7-derived halo ’’ (z=10) mass functions being very similar at z =0. 6 9 10 11 12 13 14 15 2.2. Constructing Stellar Mass-Halo Mass Relations log(M [M ]) with an Intrinsic Scatter halo sun Weobtainthestellarmassofmodelgalaxiesinthesim- Fig. 1.— Stellar mass Ms versus halo mass Mhalo relations ulation using semi-empirical modeling. This approach testedinthis paper. Included arethemodelsofGuoetal.(2010) and Behroozietal. (2013), as well as a linear SMHM relation defines the stellar mass M of the galaxy as a function s Ms = 0.025Mhalo. The low and high mass tails are necessarily of the dark matter mass of the halo, Ms = Ms(Mhalo), extrapolations of the fits due to the limitedresolution of the DM with the function Ms(Mhalo) defined to be the stellar simulations(logMhalo<10.8fortheGuoetal.(2010)modeland mass-halo mass (SMHM) relation. Although this is a logMhalo<10.3 forthe Behroozietal. (2013) model). Note that thislastmodelevolves withredshift,withM∗ (characterizing the simple one-to-one relation, it can be readily modified to kneeintheSMHMrelation)becomingsmallerathigherredshift. include statistical scatter or redshift dependence (e.g., Behroozi et al. 2013). In this paper we adopt several SMHM relations from the literature, and use them for ingModel3,weadoptthefullredshiftdependenceofthe building our simulated images. We introduce a novel, Behroozi et al. SMHM relation. Our Model 4 does not self-consistent approach that naturally adds scatter to involve halo abundance matching, but is a simple linear the SMHM mass relation. relation given by M = 0.025M , where the slope has s halo As a measure of the halo mass, we use the virial been chosen by eye to coincide with the other relations mass M (mass enclosed inside the maximum radius in the vicinity of M∗. We include this last model not vir within which the mean density is 200 times the critical because it is realistic but to study the sensitivity of our value), which is obtained from the value-added catalog results to a SMHM relation that is very different from ofGuo et al.(2010)basedonthe milli-Millenniumsimu- those of our first three models (based on the results of lation. Inthiscatalog,themassofacentralhaloisgiven Guo et al. 2010 and Behroozi et al. 2013) by M , whereas for a satellite halo it is the maximum Inthispaper,wedonotexplicitlyimposerandomscat- vir M ever attained before becoming a satellite. In the terintheSMHMrelation. Insteadweexplorethescatter vir semi-empirical approach, this preserves the stellar mass that emerges naturally from the dark matter simulation of satellite galaxies even while the outer parts of their as a consequence of one basic assumption: we assume darkhalosarebeingtidallystrippedastheyorbitwithin that the stellar mass along merger trees is a monotoni- a central halo. cally increasing function of time. Thisisphysicallyplau- Figure 1 shows some of the SMHM relations found sible because the stellarmass is concentratedinthe cen- in the literature. Our reference model (or Model 1) ter of halos due to dissipation in the baryons, and as a forsimulatingtheuniverseusestheredshift-independent resultittends tobe retainedduringmergers. Thisis the SMHM relation from Guo et al. (2010). They computed simplestphysicallymotivatedrulewehavefoundforcre- this relationusinga haloabundancematchingtechnique ating consistent galaxy stellar masses from dark matter that equates quantiles of the low redshift stellar mass simulations. The alternative is to assume that galaxies function from SDSS (Li & White 2009) to those of the loseandgainmasswilly-nillyashalomassesdecreaseand present-day halo virial mass function from the Millen- increase and as sub-halos merge; that appears much less niumsimulations. The quantile matchingis made overa plausible based on our current understanding of galaxy range of masses around the typical value M∗ (i.e., the formation and dynamics. Another alternative is to im- “knee” in the SMHM relation as well as in the mass pose scatter on the SMHM relation in a predetermined function), with extrapolations to the very low and high manner (as in the approach adopted by Behroozi et al. mass regimes, where both the stellar and halo mass dis- (2013)) tributions are not well constrained. Another option is Our assumption of monotonically increasing stellar the redshift-evolving SMHM relation of Behroozi et al. mass naturally leads to scatter in the SMHM relation (2013). This relation is based on stellar mass functions as a consequence of the following three related effects: andcosmicstarformationratesuptoz =8. Forournon- evolving Model 2, we adopt the Behroozi et al. (2013) 1. In dark matter simulations, individual halos can SMHM relation evaluated at z = 0, while for our evolv- decrease in mass from one time step to the next, 4 Taghizadeh-Poppet al. for example in events of tidal stripping. In that 4 Halo R case, we follow the approach of Guo et al. (2010) vir Galaxy R anddonot reducethestellarmassaccordingly,but 50 retain the stellar mass present before the decrease in M . 3 halo 2. In halo mergers, the dark matter mass of a de- scendant M can be smaller than the sum halo,desc ) of the progenitor dark matter masses ΣMhalo,prog (x2 P as some particles become unbound during the col- lision. This is not a rare occurrence in the simula- tions. Inordertoconservethestellarmasscontent, wemustbreakwiththeone-to-onefixedSMHMre- 1 lation. 3. Observed SMHM relations are intrinsically non- linear (see Figure 1). That leads to a conflict 0 with the assumed monotonic growth of M in halo s −1.0 −0.5 0.0 0.5 1.0 mergers. Consider the case when two halos with x = log(R) − <log(R)> M ∼ M∗ merge to create a new halo halo,prog halo with M > M∗ . The SMHM relation halo,desc halo Fig.2.— Probability density distributions of the (median sub- increases more slowly than linearly above Mh∗alo, tracted)logarithmsofRvir(darkmatterhalovirialradius)andR50 which means that M for the new halo is supposed (SDSS galaxy r-band half-Petrosian-flux radius). Note that the s to be smaller than the sum of the stellar masses radiusdistributionsresembleGaussiancurves,asexpected forap- proximately log-normal random variables. The median values are for the merging halos. That conflicts with the as- hlogRvir[Mpc]i=-1.159and hlogR50[Mpc]i=-2.824, withstandard sumption that the stellar mass along the merger deviationsσ(logRvir[Mpc])=0.176andσ(logR50[Mpc])=0.285. tree must increase monotonically. To eliminate these conflicts we adjust stellar masses have forced the stellar mass to increase monotonically retroactivelyas follows. If a halo is found to have a stel- with time, negative star formation rates are automati- larmass(accordingtotheimposedSMHMrelation)that cally excluded. decreaseswithtime, wedecreasethe stellarmassesofits Wethussimply reconstructthe spectrumofthe model immediateprogenitorstoenforceamonotonicincreasein galaxy (and derived photometry) as the sum of a series thestellarmasscontentofdarkmatterhalosacrosstheir of uniform starbursts between each of the time steps in merger trees. This adjustment is applied recursively to the simulation. When a galaxy is placed in a simulated all the progenitorsof halos with modified stellar masses. image at a particular redshift z (implying an age t ), g g Note that we are not assuming reverse causality with its rest-frame luminosity as a function of wavelength is this scheme! Our assumptionis that lowermasshalos in computed using the star-formation history up to time denserenvironments(whicharegoingtomergeinthefu- t using the evolutionary stellar synthesis code GALAXEV g ture)havetheirstarformationratessuppressedbythese by Bruzual & Charlot (2003). Note that the redshift z g environments. Thisprocedureisdescribedinmoredetail is not required to be at one of the discrete simulation in Appendix A. timesteps,sincewecanintegratethestarformationrate One natural consequence of this procedure is that, at from the previous time step to the actual time t . (The g agivenhalomass,there isa scatterinthe stellarmasses need of similar interpolation schemes has also been no- that tends to fall slightly below the one-to-one SMHM ticed by Yip 2010). In our reference model, we adopt relationatsomepointsinthemergerhistory. Theresults the Chabrier initial mass function, with a fixed solar withthemodifiedSMHMrelationswillbeshownin§3.1. metallicity(Z =0.02)andthestandarddustmodelfrom Charlot & Fall (2000) (τ = 1 and µ = 0.3). Our mod- ν 2.3. Illuminating Galaxies in Dark Matter Halos eling of galaxy spectra is flexible, as we can in principle elaboratethis model with additional variables,suchas a Theluminosityandspectrumofamodelgalaxyatany redshift-dependent metallicity. time are determined by its star formation history com- putedalongthepastmergerhistoryofitshostdarkmat- 2.4. Deriving Galaxy Sizes from Dark Matter Halo Sizes terhalo. Starformationisimpliedwhenthe stellarmass ofadescendanthaloisgreaterthanthetotalstellarmass We adjust the size of a model galaxy so that it is al- ofitsprogenitorsinconsecutivesimulationtimesteps,or waysafixedfractionofthe evolvingsize ofits darkmat- when a single halo increases its dark matter mass (and ter halo (e.g., Fall & Efstathiou 1980; Kravtsov 2013). hence its stellar mass) due to accretion of surrounding We determine the constant of proportionality between darkmatterparticles. Thestellarmassincreasebetween the galaxy size and the halo size by comparing their re- simulationtimestepsimpliesastarformationrate,which spective z = 0 distributions, as in the halo abundance isusedinstellarpopulationsynthesismodelstocompute matching method. However, we match only their medi- the emittedspectrumofthe galaxy. We assumethatthe ans instead of the full distributions, which suffer from starformationratebetweentimestepsisuniform,sothat incompleteness in the tails. It is well known that size the star formation history is completely determined by distributions for halos and galaxies are close to being the stellar mass history of a halo. Note that since we log-normal (e.g., Shen et al. 2003). Since we choose to Simulating Hubble Images With Galaxy Formation Models 5 make galaxy size linearly proportional to halo size, the improvementovermodelsthatusesmoothanalyticalpro- scale parameter and the median subtracted logarithmic files. distributions of both distributions should be the same. Figure 2 shows the distributions of halo virialradii from 2.6. Assembling the Simulated Image Millennium and r-band half Petrosian flux radii R50 for Once model galaxy properties are calculated, we gen- our main sample of SDSS galaxies from §2.5. The latter erate pencil-beams through our simulated volume and has been corrected for incompleteness using the VMAX project them onto the plane of the sky. We then simu- method (Schmidt 1968), as shown with a similar galaxy late ACS/WFC camera images, with their visible filters, sample by Taghizadeh-Popp et al. (2012). We find the as well as corresponding WFC3/IR images (sampled to relation R50 = 0.022Rvir. The dispersions are different theACS/WFCpixelsize),withtheirinfraredfilters. We (indicating that our assumption is not completely accu- include realistic PSFs for both cameras. rate) but are similar enough for our purposes. As in the To determine the 3-D structure within these pencil case of the spectra, we interpolate the galaxy size be- beams, we use a Monte Carlo method. This approach tween discrete time steps to the redshift zg of the model is much simpler and faster than other alternatives, such galaxy. asthe replicationof a simulationbox muchsmaller than thedepthspannedbythesimulatedimage. Wefirstsam- 2.5. Galaxy Image Cutouts from SDSS ple a random redshifts z from a distribution that gives g Ourmethod placescutouts ofrealgalaxyobservations constantcomovingvolumeperredshiftinterval. Thenwe ontooursimulatedimage,whichhasadvantagesoverus- randomly select a dark matter halo at z = 0, choose all ing smooth analytic galaxy light profiles. Galaxies often ofitsprogenitorsfoundatz ,andplacetheminthesim- g have a clumpy internal structure, which affects source ulated image (at redshift z ) according to their relative g detections. A real galaxy could be detected as two or 3-D spatial positions in the simulation box (interpolat- moreseparateobjects,especiallyifits surfacebrightness ing between time steps). Although the large-scale cor- approaches the image noise level. Real galaxy cutouts relations are discarded, we still preserve the short-range recreate this effect and are more realistic than smooth correlations between progenitors, while reducing consid- bulge and disk light profiles adopted in previous semi- erably the computation time. empirical and semi-analytical models. Of course, this Thevisualappearanceofthemodelgalaxyonthesim- effectisnotasimportantwhentheapparentgalaxysizes ulatedimage is givenby the best matchingSDSS galaxy arecomparabletothewidthofthepointspreadfunction, cutout (Appendix B). We select the SDSS band whose as may happen in the high redshift limit. redshifted central wavelength λ (to redshift z ) is the c g A database of SDSS galaxy images is used for the closest to the band of the simulated image. Due to the cutouts. We select from the SDSS image database the redshift of wavelengths, we end up using the bluest vis- real galaxy that is the closest neighbor to the model ible SDSS bands to represent most high-redshift simu- galaxy in a multi-dimensional space of galaxy proper- lated galaxies in the visible and infrared HST filters, ties. Once the closest matching SDSS galaxy is found, since SDSS lacks the UV and far UV bands that would werescaleitsfluxandsizeinordertomakethemequalto be more suitable for this redshift regime3. Since the u those of the model galaxy. No free parameters are fitted band is noisy in SDSS, we use the g band as the bluest orrequiredinthisstep. ThefulldetailsoftheSDSSdata bandpass. Infact,allmodelgalaxiesplacedonthesimu- and the selection procedure is described in Appendix B. lated image are represented by a g-band cutout at z >2 One might wonder whether SDSS galaxies are clumpy for ACS/WFC images or z >3 for WFC3/IR images. enough to be accurate models of high-redshift galaxies, We rescale the flux of the image cutout to match that since galaxies at higher redshifts tend to be more irreg- of the model galaxy. We apply to the model spectrum ular than local galaxies. The changes in morphologyare theeffectsofredshift,cosmologicaldistancedimmingand due both to the shift ofopticalbandfilters intothe rest- intergalactic absorption (as in Madau 1995), before ap- frame ultravioletfor distant objects and also to a higher plying the photometric filter response. merger rate and dynamically less-relaxed structures in Werescalethesizesofthegalaxycutoutsplacedonthe the early universe. simulated image according to one of the following rules: While these effects are worthy of further exploration inthe future, forthis paper the SDSS imagesare a good 1. No size scaling: The proper size of the original basis for simulations. Our analysis of galaxy counts and SDSS physical galaxy is left intact, irrespective of detection efficiencies relies on observations and simula- the model galaxy size and hence the size of its tions in the HST WFC3/IR F160W filter (λ= 1.6µm); halo. The physical size is scaled to the apparent this implies that objects with redshifts less than 3 have size on the image using the angular diameter dis- anSDSSfilter(fromgriz)thatisattheappropriaterest- tance DA(zg). frame wavelength, and the different morphologies in the 2. Size scaling: The proper size of the SDSS galaxy ultravioletareirrelevant. Moreover,galaxiesatredshifts is scaled to match the physical size of the model beyond 3 tend to be sufficiently compact that their in- galaxy,whichin turn is a fixed fractionof the halo ternalclumpinesshaslittleimpactontheirdetectability. size(asdescribedabove). Theapparentsizeisthen The median FWHM size of detected galaxies with z >3 computed from D (z ). in our simulations is 0.3 arcsec, which is only twice the A g FWHM of the 1.6µm WFC3 PSF. While there is room 3 A cross-match between SDSS and GALEX might be useful for improvement in modeling the internal structure of in the future for adding ultraviolet bands to our suite of galaxy galaxies, particularly for bluer filters in the 1 < z < 3 cutouts,butthatisoutofthescopeofthisexploratorystudyand redshift range, the SDSS cutout images are certainly an analysis. 6 Taghizadeh-Poppet al. TABLE 1 GalaxyEvolution Models Model(shorttitle) Details Model1 • Stellarmass-halomassrelationfromGuoetal.(2010) (Reference Model) • DustmodelfromCharlot&Fall(2000) • Fixedsolarmetallicity(Z =0.02)atallredshifts • Apparent sizes of SDSS galaxy cutouts on image are scaled to the theoretical size pre- dictedbythedarkmatter halosize • SDSS galaxy cutouts are chosen to be the closest match to the theoretically predicted modelgalaxyu−r colorandstellarmass. Model2 SameasModel1,butusingtheBehroozietal.(2013)SMHMrelationevaluatedatz=0. Model3 SameasModel1,butusingtheBehroozietal.(2013)redshiftdependentSMHMrelation. Model4 SameasModel1,butusingthealinearSMHMrelationMs=0.025Mhalo. Model5 SameasModel1,butusingnodustmodelforgalaxies. Model6 SameasModel1,butusingafixedverylowmetallicity(Z =0.0001)atallredshift. Model7 Same as Model 1, but without rescalingthe intrinsicsize of SDSS galaxy cutouts (except fortheangulardiameterdistancescaling,alsoappliedinthereferencemodel). Model8 SameasModel1,butthePetroR50radiusaswellasu−r andMs areusedformatching SDSStomodelgalaxies. Note. —Descriptionofeightdifferentmodelsusedforbuildingasimulateduniverse. Ourreferencemodel(Model1)contains themostplausibleparametersandsub-models. Othermodelsaredefinedbychangingoneoftheseparametersatatime. Inthe firstmodel, the z=0size-massrelationofgalax- galaxy surveys we use AUTO magnitudes and Petrosian ies is preserved at all redshifts, whereas in the sec- radii R , as calculated by SExtractor. The Petrosian P ond model, the size-mass relation evolves with redshift, radii of our input galaxy cutouts, as given by the SDSS drivenbythegrowthofthedarkhalos. Galaxysizestend pipeline, differ from those returned by SExtractor on tobesmallerathigherredshiftsinthesecondmodeldue the already simulated images, probably due to different to evolution in the halo sizes. definitions for the radius in the SDSS pipeline and in Toaddinstrumentaleffects,weconvolvewiththeHST SExtractor. Since we observe a linear offset between point spread function, apply the HST instrument ef- the distributions of these radii, we change the input ficiency and pixelation, and add noise. The HST in- SDSS values to match those from SExtractor by using strument modeling uses the same software (pysynphot, logR (SExtractor)=logR (SDSS)+0.41. P P Tiny Tim),instrumentalparametersandskybackground asusedinthestandardHST ExposureTimeCalculators, providing a high fidelity model of the real HST perfor- mance. 2.7. Source Detection and Photometry of Simulated Images The detection, extraction, and photometry of galaxies are performed by running SExtractor (Bertin & Arnouts 1996) on the simulated images. Here we directly follow the method and parameters described in Galametz et al. (2013) and designed for the analysis of real HST images. The complete output catalog merges SExtractor runs using two different detection modes. The Cold mode is used for detecting 3. RESULTS bright and extended sources, while the Hot mode The following subsections present our results for eight is optimized for extracting faint and small objects. different models of simulated universes. These models After extraction, detected or output galaxies on the are detailed in Table 1, where the reference model in- final simulated image are matched (using the detected volves the most common choices of parameters as ex- position and luminosities) to the input galaxies on the plainedintheprevioussections,andthesevenremaining same image before adding the instrumental effects. This models vary one parameter of the reference model at a provides a direct measure of the detection completeness time. In models 2–4 we change the SMHM relation; in by comparing what SExtractor detects to what was models 5–8 we explore extreme values of other models originally placed on the image. parameter (often deliberately unrealistic values) to test To compare our simulated galaxies to those from real the sensitivity of the observations to these parameters. Simulating Hubble Images With Galaxy Formation Models 7 Model 1 Model 1 Model 2 Model 2 3.0 4 2.5 4 2.5 0 0 0 0 1 3 1 2.0 1 3 1 2.0 8 8 1.5 8 8 1.5 2 2 1.0 1.0 6 6 6 6 1 1 ) 0.5 ) 0.5 ]sun4 All redshift 0 4 z = 1.63 0.0 ]sun4 All redshift 0 4 z = 1.63 0.0 M 10 12 14 10 12 14 M 10 12 14 10 12 14 [s [s M Model 1 Model 1 M Model 2 3.0 Model 2 2.5 ( ( log 0 2.50 2.0 log 0 2.50 2.0 1 2.01 1 1 1.5 2.0 1.5 8 1.58 8 1.58 1.0 1.0 1.0 1.0 6 6 6 6 0.5 0.5 0.5 0.5 4 z = 3.87 4 z = 7.27 4 z = 3.87 4 z = 7.27 0.0 0.0 0.0 0.0 10 12 14 10 12 14 10 12 14 10 12 14 log(M [M ]) log(M [M ]) halo sun halo sun Model 3 Model 3 3.0 Model 4 Model 4 3.0 4 0 0 2.5 0 4 0 2.5 1 3 1 2.0 1 3 1 2.0 8 2 8 1.5 8 8 1.5 2 6 6 1.0 6 6 1.0 1 1 ) 0.5 ) 0.5 ]sun4 All redshift 0 4 z = 1.63 0.0 ]sun4 All redshift 0 4 z = 1.63 0.0 M 10 12 14 10 12 14 M 10 12 14 10 12 14 [s [s M Model 3 Model 3 M Model 4 3.0 Model 4 2.5 ( ( g 2.5 2.0 g lo 0 0 lo 0 2.50 2.0 1 2.01 1.5 1 2.01 1.5 8 1.58 8 1.58 1.0 1.0 1.0 1.0 6 6 6 6 0.5 0.5 0.5 0.5 4 z = 3.87 4 z = 7.27 4 z = 3.87 4 z = 7.27 0.0 0.0 0.0 0.0 10 12 14 10 12 14 10 12 14 10 12 14 log(M [M ]) log(M [M ]) halo sun halo sun Fig.3.— Modifiedstellarmass-halomassrelationsforthefourmodelstestedinthispaper,obtainedafterretroactivelyreducingtheMs valuespredictedbytheone-to-oneSMHMrelationsinFig.1toenforceamonotonicallyincreasingMs asafunctionoftimealongmerger trees. Datafromthreeredshifttimestepsaswellasthecombinationfromalltimestepsisshown. Thecolorsshowthelog-scalednumber counts inbinsofsize∆logMhalo=0.06by∆logMs=0.01. 0 0 0 0 1. 1. 1. 1. x) x) x) x) < < < < P(X0.5 P(X0.5 P(X0.5 P(X0.5 Model 1 Model 2 Model 3 Model 4 0 0 0 0 0.0 1 2 0.0 1 2 0.0 1 2 0.0 1 2 x = D log(Ms[Msun]) x = D log(Ms[Msun]) x = D log(Ms[Msun]) x = D log(Ms[Msun]) Fig. 4.— Cumulative distributions of the difference between the modified SMHM relation shown inFig. 3 and the original one-to-one relationfromFig. 1(values fromall redshifttime steps areincluded). The50%, 67% and 95% quantiles areindicated by dashed vertical lines. Abouthalfofthehalosineachmodelhavetheirmassesunmodified(x=0). 8 Taghizadeh-Poppet al. 1 1 -0 -0 1 1 2 2 3c-10 3c-10 p p M M M/s-30 L/r-30 g 1 g 1 o o l Model 1 l DN/4 Model 2 DN/4 Model 1 D-10 MDMod heal l4o Millenium II D-10 MMooddeell 23 Model 4 5 DM halo Milli−Millenium 5 Model 5 -0 SDSS, Li (2009) -0 Model 6 1 1 SDSS at z=0.1 SDSS, Baldry (2008) 5 6 7 8 9 10 11 12 5 6 7 8 9 10 11 12 log(M [M ]) log L /L s sun r sun Fig.5.— Presentdaystellarmassfunctionforasimulateduni- Fig. 6.— Global luminosity functions (r-band) calculated af- verse according to the different models of stellar mass-halo mass ter evolving our simulated universe to z = 0.1 under six different relations tested in this paper: Model 1 (Guoetal. 2010), Model models. AlsoshownistheluminosityfunctionoflocalSDSSgalax- 2 (Behroozi etal. 2013) and Model 4 (a linear SMHM relation). ies measured by Blantonetal. (2005). Note the similarity with The open symbols show the stellar mass functions at z = 0 from Fig.5. Thedownturninthecomputedluminositydistributionsat themilli-MillenniumandMillenniumIIsimulations,computedby lowluminositiesisduetothefiniteparticlemassresolutioninthe converting dark matter halo masses into stellar masses using the milli-Millenniumsimulation. Guoetal. (2010) relation. Observed stellar mass functions mea- sured from local SDSS galaxies (Baldryetal. 2008; Li&White (2010) SMHM relation, agrees with the observed stel- 2009)arealsoshown. lar mass function from Li & White (2009) in the range 3.1. Modified SMHM Relations with Natural Intrinsic Ms > 109M⊙. That is expected since this SMHM re- Scatter lation was derived from halo abundance matching on mostly the same data. We also confirm that the stellar Asdescribedin§2.2andAppendix A,wemodifiedthe mass function predicted by the reference model fits per- fixed SMHM relation, selectively reducing the values of fectlythestellarmassfunctionderivedbycombiningthe M assignedto halos to enforce a monotonically increas- s Guo et al. SMHM relation with the halo mass function ing M along merger trees. This leads to the natural s fromthemilli-Millenniumsimulation,asrequiredbyour dispersion of M values shown in Figures 3 and 4 for s semi-empiricalmodeling. However,thestellarmassfunc- the various SMHM relations explored in this paper. For tion of the reference model has an artificial downturn at all models, about half of halos do not have their stellar Ms < 107M⊙, which corresponds to the mass resolution masses modified (meaning that they lie on the imposed of the milli-Millennium simulation in the identification SMHM relation). The scatter in the M -M distri- s halo of friends-of-friends groups (composed of a minimum of butions reaches about 0.07–0.12 dex in our simulations, 20darkmatterparticles). Similardownturnsarepresent similartothe scatterinferredindirectly fromtheoryand in the stellar mass functions of the other SMHM rela- observations (e.g., Reddick et al. 2013; Behroozi et al. tions tested. In contrast, the Millennium II simulation 2013). The linear SMHM relation (Model 4) is the least (Boylan-Kolchinet al. 2009), which has a dark matter affectedbyourmodificationsandshowsthesmallestdis- particle mass 125 times smaller than the mMS, shows a persionbecauseitlackstheSMHMnon-linearitypresent power-law tail at the low-mass end. in the other models. For the non-linear Models 1–3, the The halo abundance matching in the Guo model cov- scatterissmallestnearMhalo ∼1012M⊙becausethestel- ered the range 108.3 < Ms/M⊙ < 1011.8, which misses lar mass corrections are minimized near the peak of the SMHM relation. the upturn at Ms < 109M⊙ that is seen in the more complete SDSS stellar mass function from Baldry et al. (2008)(shownbytheblacklineinFig.5). TheSMHMre- 3.2. Simulated Universe Models at Low Redshift lation from Behroozi et al. (2013) (Model 2, yellow line) We evolve all models of simulated universes to the fits this upturn much better since it was built using the presentday. In this section we comparethe model prop- Baldry et al. (2008) stellar mass function. erties with low-redshift observations to test the overall The linear SMHM relation (Model 4) tracks the ap- accuracy of our method. proximately M−2 power-law mass distribution of dark matter halos. The proportionality constant for Model 4 3.2.1. Stellar Mass Functions was chosen to roughly match the observed stellar mass Figure 5 compares the stellar mass function of model function around its knee at M ∼ M∗. Again, a clear s galaxiesatz =0tothosederivedfromSDSS.Thestellar mass resolution cutoff is present at lower stellar masses, massfunction ofthe referencemodel, with its Guo et al. alsoshowninthenon-linearrelations. ThislinearSMHM Simulating Hubble Images With Galaxy Formation Models 9 model is obviously in strong conflict with the observa- ulatedvisible(ACS/WFC)andinfrared(WFC3/IR)im- tions at both higher and lower masses. agesusing filters and exposure times appropriatefor the GOODS, CANDELS and HUDF surveys. Image sizes 3.2.2. Luminosity Functions andpixelscalesarethoseofasingleACS/WFCexposure (a 200′′×200′′ field with 0.0495′′ pixels). Note that this Figure 6 compares the simulated and observed lumi- image area is small compared with the areas surveyed nosityfunctionsatz =0.1. Tofirstorder,theshapesare by many HST projects (e.g., GOODS and CANDELS), Schechter functions. The present-day luminosity func- which encompass many ACS/WFC fields; we considered tionsfromourmodelsaresimilarinshapetotheirrespec- it appropriate to begin with a more modest sky area for tive stellar mass functions in Figure 5. This is expected this exploratory project. since the r-band luminosity roughly traces the old stel- lar population that constitutes most of the stellar mass in present-day galaxies (e.g., Bell et al. (2003)). Mod- 3.3.1. Results from the Reference Model els 5 (no dust) and 6 (low-metallicity) show good agree- Figure 7 shows a comparisonbetween reference-model ment with the Blanton et al. (2005) observed luminos- simulated and real HST ACS/WFC images. At first ity function in the high-luminosity tail. Since our refer- glance, the spatial distribution of galaxies and associ- encemodelcontainsdustandmetals,itshigh-luminosity ated sizes seem to be very similar. As expected, the tail is shifted towardfainter luminosities with respect to HUDF-depth imagesshowthatmanyobjects arehidden those of Models 5 and 6. (Note that the luminosity in by noise in the GOODS-depth simulated images. the r-band is reduced by dust absorption even though Figure 8 shows scatter plots of the input (true model) the bolometric luminosity is conserved.) On the other values of the apparent sizes and surface brightnesses of hand,theflatterslopeatlowluminosities forthesethree bothdetected(usingSExtractor)andundetectedgalax- models mimics the flat tail present in their correspond- ies in the HUDF-depth image. Most of the detected ingstellarmassfunctions,whichistheresultofapplying galaxiescanbeselectedviaacutatm .29,which F160W halo abundance matching to the stellar mass function of doesnotdependstronglyonredshift. Inthelowerpanels Li & White(2009)asexplainedearlier. Abetterfitisob- of Figure 8, which plot surface brightness versus mag- tainedfor Models 2and3 with the SMHM relationfrom nitude, it can be seen that there is a small, redshift- Behroozi et al.(2013),whichincludesthelow-luminosity dependent population of low surface brightness galax- upturn. The luminosity function derivedfrom the linear ies that are brighter than the magnitude threshold but SMHM relation for model galaxies tracks the approxi- nonetheless are not detected. mately power-law mass distribution of dark matter ha- Figure9plotsthestellarmassesofgalaxiesversustheir los. input model apparent magnitudes. The stellar mass de- The perceptive reader may have noticed that even creases strongly with increasing redshift. However, an though semi-empirical modeling promises a perfect importantconclusionfromthisplotisthateventheleast match of the predicted to observed universe, the global massive detected galaxies are still generally above the luminosity function predicted by our reference model mass thresholds set by both the SDSS survey (M & s does not agree perfectly with the observed luminosity 106M⊙,Kauffmann et al.2003)andthemilli-Millennium function from SDSS at z ∼ 0. The reason for this is that the z ∼0 SMHM relation from Guo et al. (2010) is simulation (Ms & 107.1M⊙, Guo et al. 2010). Figure 10 shows the typical apparent magnitude values found at not fully consistent with our method for convertingstel- themasscompletenesslimitofthemilli-Millenniumsim- lar mass into light. They use the stellar mass function ulation as a function of redshift. Although the median derived from Li & White (2009), which is not directly magnitude initially dims as redshift increases, the curve observed but is inferred from the observed SDSS galaxy flattens out at m ∼ 32 for z > 3; although the luminosity function at z ∼ 0 and an assumed mass-to- F160W galaxiesatthoseredshiftshavesmallermasses,theyalso light ratio. On the other hand, our mass-to-light ratio have higher star formation rates and younger popula- is computed from the dark matter halo masses and the tions due to the requirementthat their stellar masses be SMHMrelationusingBruzual & Charlot(2003)spectral assembledina shorttime. Thatapproximatelycompen- models. We can in principle address this issue by com- satesforthe cosmologicaldimming dueto theincreasing putingourownSMHMrelationusinganiterativeprocess luminosity distance. that compares our measuredluminosity functions to the Biases — As we noted above, Figure 8 shows that a observations. This level ofprecisionis not needed in our relatively simple cut on the input (model) magnitude of present exploratorystudy, but it will be a useful longer- galaxiespredictswithreasonableaccuracywhichgalaxies term goal for our approach. will be detectable in the images. However, detection is not the whole story. It is also necessary to measure the 3.3. Simulated HST Images and Derived Statistics galaxyproperties,andthosemeasurementscanbebiased In this section, we show simulated HST images from through complex effects related to the morphology and our models and test the sensitivity of statistics derived internal structure of the galaxy. Detection of an object from the images to the model parameters. We focus on doesnotensurethatitsmagnitude,colors,orsizecanbe theluminosityandsizedistributions,andweassessboth measured correctly. This is an area where our forward- biases in measured parameters and source detection in- modeling approach provides more reliable results than completeness. For most tests we compare the perfectly previous approaches. known input values of size and luminosity (as given by Figures 11 and 12 show strong biases in the themodels)tothecorrespondingoutput valuesmeasured SExtractor-derived measurements of galaxy sizes and bySExtractorfromtheimages. Ourdatacomprisesim- apparent magnitudes. From galaxies detected in the 10 Taghizadeh-Poppet al. SimulatedGOODSImage RealGOODSimage ReferenceModel SimulatedHUDFImage RealHUDFimage ReferenceModel Fig.7.—SimulatedACS/WFCF850LP+F606W+F435Wimagesbuiltfromourreferencemodel(GOODSandHUDFdepths),compared totheequivalentrealHST images. Thesameexposuretimesanddisplaycontrastareusedforthecomparisonateachdepth. Thefieldof viewis2400×1200ACSpixels(∼1/6ofthefullACS/WFC fieldofview). HUDF-depth simulated image, we compare the dif- extended objects may be detectable only if noise ferences between the true (input) and SExtractor- fluctuations make their surface brightness appear measured (output) values of m and logR . Our brighter than the detection limit. F160W P major conclusions are: Figure 11 also reveals a small bias toward larger mea- 1. There is a significant magnitude-dependent lumi- suredoutputsizes forsources∼1 magbrighterthanthe nositybias,withthemeasuredm magnitudes detection limit. This bias is largest in the z > 6 red- F160W on average fainter than the true galaxy magni- shift bin, where the difference reaches around 0.2 dex, tudes (Fig.11). Forfainter galaxies,onlythe com- andmightbe relatedto intrinsic peculiarities inthe way pact, high surface brightness nuclei are detectable SExtractormeasures sizes. above the image noise level, while the lower sur- These effects are certainly detection biases and not face brightness extended envelopes are hidden in a trend resulting from supposed evolutionary effects in the noise. This bias is small for bright galaxies, the models. Figure 12 shows a scatter plot of the mea- but reaches median values of ∼0.4 magnitude and sured(output)sizesversusmeasuredmagnitudeforboth extreme values of & 1 magnitude at the detection GOODSandHUDF-depthsimulatedimages. Thestrong limit. This effect does not depend stronglyon red- bias in the SExtractor-measured sizes at the magni- shift. tude detection limit are easily visible at both depths, and the bias begins at a magnitude that is determined 2. There is also a strong bias in the measured sizes by the detection limit rather than at any physically de- near the magnitude detection limit (Fig. 11). termined luminosity. Note that the underlying images The SExtractor-measured sizes are smaller than (before adding noise) are exactly the same in these two the original input sizes of model galaxies around cases: they are equivalent to observing the same sky re- m ≃ 29. The explanation appears to be the gion twice with different exposure times. This rules out F160W sameasfortheluminositybias: thedetectablepart anyartifactintroducedbythemodel. Notehowthemea- ofafaintgalaxyissignificantlymorecompactthan sured R extends down to values close to the PSF size P the true Petrosian radius due to masking of the (FWHM=0.151′′fortheF160Wbandimage)atthemag- more extended component by noise. The differ- nitude detection limit. ence between output and input values of logR We notethatthe truebiasescouldbe evenlargerthan P has a median of 0.15 dex, with maximum values those measured in our simulations. Real galaxy light reaching even to 0.8 dex. Note that there is also a profiles have extended wings, whereas our galaxy image Malmquist bias affecting these distributions: faint cutouts only include the light present within two Pet-