Simulating galaxies in the reionization era with FIRE-2: morphologies and sizes
We study the morphologies and sizes of galaxies at using high-resolution cosmological zoom-in simulations from the Feedback In Realistic Environments project. The galaxies show a variety of morphologies, from compact to clumpy to irregular. The simulated galaxies have more extended morphologies and larger sizes when measured using rest-frame optical B-band light than rest-frame UV light; sizes measured from stellar mass surface density are even larger. The UV morphologies are usually dominated by several small, bright young stellar clumps that are not always associated with significant stellar mass. The B-band light traces stellar mass better than the UV, but it can also be biased by the bright clumps. At all redshifts, galaxy size correlates with stellar mass/luminosity with large scatter. The half-light radii range from 0.01 to 0.2 arcsec (0.05–1 kpc physical) at fixed magnitude. At , the size of galaxies at fixed stellar mass/luminosity evolves as , with –2. For galaxies less massive than , the ratio of the half-mass radius to the halo virial radius is and does not evolve significantly at –10; this ratio is typically 1–5% for more massive galaxies. A galaxy’s ‘observed’ size decreases dramatically at shallower surface brightness limits. This effect may account for the extremely small sizes of galaxies measured in the Hubble Frontier Fields. We provide predictions for the cumulative light distribution as a function of surface brightness for typical galaxies at .
keywords:galaxies: evolution – galaxies: formation – galaxies: high-redshift – cosmology: theory
High-redshift galaxies are thought to be the dominant source of cosmic reionization (e.g. kuhlen.faucher.2012:reion; robertson.2013:reion.hudf12; robertson.2015:reion.planck). The number density of these galaxies, as described by the ultraviolet (UV) luminosity function, is well constrained for galaxies brighter than at (e.g. mclure.2013:uvlf.z7to9.hudf12; schenker.2013:uvlf.z7to8.udf12; bouwens.2015:uvlf.z4to10; finkelstein.2015:uvlf.combined.field). Recently, the Hubble Frontier Fields (HFF) program (lotz.2017:hubble.frontier.field), which takes advantage of strong gravitational lensing by foreground galaxy clusters, has made it possible to estimate UV luminosity functions down to (e.g. bouwens.2017:lensing.uncertainty; livermore.2017:faint.galaxies). But one of the dominant outstanding uncertainties is the intrinsic size distribution of faint galaxies, which is necessary in order to determine the completeness of the observed sample due to surface brightness limits (bouwens.2017:small.galaxy.sizes).
There are only a few galaxies at these redshifts that have robust size measurements. oesch.2010:hudf.morphology measured the sizes of galaxies brighter than at –8 in the Hubble Ultra-Deep Field (HUDF) and found that the half-light radii of galaxies evolve according to , with –1.5 (see also, e.g. bouwens.2004:hudf.size.evolution; ferguson.2004:size.evolution; ono.2013:size.evolution.udf12; kawamata.2015:hff.size.z6to8). It is also expected from analytic models that galaxy size decreases with increasing redshift (mo.mao.white.1998:galactic.disk; mo.mao.white.1999:lbg.structure).
High-redshift galaxies tend to be intrinsically small. The half-light radii of bright galaxies () at –8 range from 0.5–1 kpc (e.g. oesch.2010:hudf.morphology). More recently, kawamata.2015:hff.size.z6to8 and bouwens.2017:small.galaxy.sizes measured the half-light radii for a sample of fainter galaxies () from the HFF. They found that the size–luminosity relation has large scatter, with half-light radii spanning more than an order of magnitude (0.1–1 kpc) at fixed UV magnitude. A fraction of these faint galaxies have extremely small sizes from a few pc to 100 pc, although these results are very uncertain because they are far below the resolution of the Hubble Space Telescope (HST).
Morphological studies of intermediate-redshift galaxies (–3) reveal that their images typically contain a number of star-forming clumps (e.g. guo.2015:clumpy.galaxy.candels), which only contribute a small fraction of the total stellar mass (e.g. wuyts.2012:candels.clumps). So far, the sizes of galaxies are measured using noise-weighted stacked images over all available bands (dominated by rest-frame UV), so it is likely that the extremely small galaxy sizes in the HFF are biased by such clumps (e.g. vanzella.2017:globular.formation). With the James Webb Space Telescope (JWST, scheduled to launch in 2018), one will be able to probe these faint, high-redshift galaxies with deeper imaging, higher spatial resolution, and at longer wavelengths. This makes it possible to compare galaxy morphology and sizes in different bands and to recover the stellar mass distribution using multi-band images via pixel-by-pixel spectral energy distribution (SED) modeling (e.g. smith.hayward.2015:sed.modeling).
The goal of this paper is to make predictions of morphology and sizes for galaxies, which can be used to plan and interpret future observations. We use a suite of high-resolution cosmological zoom-in simulations from the Feedback In Realistic Environments (FIRE) project111http://fire.northwestern.edu (hopkins.2014:fire.galaxy). The FIRE simulations include explicit treatments of the multi-phase interstellar medium (ISM), star formation, and stellar feedback. The simulations are evolved using the FIRE-2 code (hopkins.2017:fire2.numerics), which is an update of the original FIRE code with several improvements to the numerics. These simulations predict stellar mass functions and luminosity functions in broad agreement with observations at these redshifts. When evolved to , simulations with the same physics have been shown to also reproduce many other observed galaxy properties (hopkins.2017:fire2.numerics, and references therein).
The paper is organized as follows. In Section 2, we describe the simulations briefly and the methods we use to measure galaxy sizes. We present the results in Section LABEL:sec:results and discuss their implications to observations in Section LABEL:sec:discussion. We conclude in Section LABEL:sec:conclusion. We adopt a standard flat CDM cosmology with Planck 2015 cosmological parameters , , , , , and (planck.2016:cosmo.param). We use a kroupa.2002:imf initial mass function (IMF) from 0.1–, with IMF slopes of from 0.1– and from 0.5–. All magnitudes are in the AB system (oke.gunn.1983:ab.system).
2.1 The simulations
We use a suite of 15 high-resolution cosmological zoom-in simulations at from the FIRE project, spanning a halo mass range – at . These simulations are first presented in ma.2017:fire.hiz.smf. The mass resolution for baryonic particles (gas and stars) is – (more massive galaxies having larger particle mass). The minimum Plummer-equivalent force softening lengths for gas and star particles are –0.42 pc and –2.1 pc (see table 1 in ma.2017:fire.hiz.smf for details). The softening lengths are in comoving units above but switch to physical units thereafter. All of the simulations are run using exactly identical code gizmo222http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html (hopkins.2015:gizmo.code), in the mesh-less finite-mass (MFM) mode with the identical FIRE-2 implementation of star formation and stellar feedback.
The baryonic physics included in FIRE-2 simulations are described in hopkins.2017:fire2.numerics, but we briefly review them here. Gas follows an ionized-atomic-molecular cooling curve from – K, including metallicity-dependent fine-structure and molecular cooling at low temperatures and high-temperature metal-line cooling for 11 separately tracked species (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe). At each timestep, the ionization states and cooling rates H and He are calculated following katz.1996:sph.cooling, with the updated recombination rates from verner.1996:recombination.rates, and cooling rates from heavier elements are computed from a compilation of cloudy runs (ferland.2013:cloudy.release), applying a uniform but redshift-dependent photo-ionizing background from cafg.2009:uvb, and an approximate model for H ii regions generated by local sources. Gas self-shielding is accounted for with a local Jeans-length approximation. We do not include a primordial chemistry network nor Pop III star formation, but apply a metallicity floor of .
We follow the star formation criteria in hopkins.2013:sf.criteria and allow star formation to take place only in dense, molecular, and locally self-gravitating regions with hydrogen number density above a threshold . Stars form at 100% efficiency per local free-fall time when the gas meets these criteria, and there is no star formation elsewhere. The galactic-scale star formation efficiency is regulated by feedback to (e.g. orr.2017:fire.ks.law). The simulations include the following stellar feedback mechanisms: (1) local and long-range momentum flux from radiation pressure, (2) SNe, (3) stellar winds, and (4) photo-ionization and photo-electric heating. Every star particle is treated as a single stellar population with known mass, age, and metallicity. The energy, momentum, mass, and metal returns from each stellar feedback processes are directly calculated from starburst99 (leitherer.1999:sb99).
2.2 Post-processing and size definition
We use Amiga Halo Finder (ahf; knollmann.knebe.2009:ahf.code) to identify halos in our simulations. The halo mass () and virial radius () are computed by ahf, applying the redshift-dependent virial overdensity criterion from bryan.norman.1998:xray.cluster. Each zoom-in simulation contains one central halo, which is the most massive halo of the zoom-in region. In this work, we also consider other less massive halos in the zoom-in region. We restrict our analysis to halos with zero contamination from low-resolution particles, which also having more than 100 star particles and total particles within the virial radius. In Figure LABEL:fig:nhalo, we show the number of halos in our simulated sample at , 8, and 10. At a given mass, our sample includes both central halos and less massive halos in the zoom-in regions, and include simulations run with different mass resolutions. In Appendix LABEL:sec:append:res, we show that our results converge reasonably well with resolution. At a given redshift, we project all star particles inside the halo along a random direction onto a two-dimensional uniform grid. The pixel size of the grid is 0.0032 arcsec (3.2 mas), which equals to 1/10 of the pixel size of JWST’s Near-Infrared Camera (NIRCam) and corresponds to 10–20 pc in the redshift range of our interest. Each star particle is smoothed over a cubic spline kernel with a smoothing length , where is its distance to the nearest neighbor star particle. We adopt as our default value, but varying –10 only makes small difference for a small fraction of our galaxies. We only consider intrinsic morphologies and sizes and do not include dust extinction in this work. The majority of galaxies (over 95%) in our sample have intrinsic UV magnitude fainter than (stellar mass ). We find that dust attenuation has little effect on these low-mass, faint galaxies (also see ma.2017:fire.hiz.smf), so most results in this paper are not affected by dust extinction. We make projected images for stellar surface density and rest-frame UV (1500 Å) and rest-frame B-band (4300 Å) surface brightness. The rest-frame UV of galaxies at these redshifts falls in the short-wavelength channel of NIRCam (observed wavelengths 0.6–2.3 m, spatial resolution 0.032” per pixel), in F090W band for galaxies and F150W band for galaxies. The rest-frame B-band falls in NIRCam’s long-wavelength channel (2.4–5 m, spatial resolution 0.065” per pixel), in F277W band for galaxies and F444W band for galaxies. The SED of each star particle is computed using the synthesis spectra predicted by the Binary Population and Spectral Synthesis (BPASS) models (version 2.0; eldridge.2008:bpass.model)333http://bpass.auckland.ac.nz. We use the binary stellar population models in BPASS by default444Although we use starburst99 single-star models for stellar feedback in the simulations, BPASS binary models are not expected to strongly affect the feedback strength, since bolometric luminosities and supernovae rates are similar between these models (see section 2.2 in ma.2017:fire.hiz.smf, for more discussion). We prefer the binary models because they are able to reproduce the nebular emission line features in high-redshift galaxies (e.g. steidel.2016:stellar.nebular.spec; strom.2017:nebular.line.mosfire).. In Figure 2.1, we show example images for six galaxies labeled by A–F. Each panel is 2” (11.6 kpc) on each side. Some galaxies and structures are so small that they only occupy very few pixels, so we further smooth the images using a two-dimensional Gaussian kernel with a dispersion equal to the size of 10 pixels (0.032”) only for easier visualization here. Most galaxies in our sample show clumpy, irregular morphologies that cannot be well described by a simple profile (see also figures 2 and 3 in ma.2017:fire.hiz.smf; cf. jiang.2013:uv.morphology.z6; bowler.2017:bright.z7.galaxies). Therefore, we adopt a non-parametric approach to define galaxy sizes, in a way similar to curtis-lake.2016:size.evolution and ribeiro.2016:galaxy.size.evolution. For every galaxy, we place a circular aperture of 1” in diameter, whose center is located by iteratively finding the B-band surface brightness-weighted center of all pixels within the 1”-aperture, as illustrated by the white dotted circles in Figure 2.1. We visually inspect all galaxies to ensure the apertures are reasonably located. The same aperture is used for the size measurement in stellar mass, UV, and B-band luminosity for the same galaxy as follows. We sort the pixels enclosed in the 1”-diameter aperture in descending order of surface density or surface brightness, and find the number of ‘brightest’ pixels that contribute 50% of the total mass or luminosity within the 1” aperture. We calculate the area spanned by these pixels and define the ‘half-mass’ or ‘half-light’ radius as . We quote the galaxy stellar mass and luminosity as the total amount enclosed in the 1”-diameter aperture555We have checked that 1” aperture is sufficiently large for most galaxies in our sample, except for the few galaxies that are at a close encounter during a merger. In fact, galaxy F in Figure 2.1 is the mostly affected object, whose stellar mass and half-mass radius are underestimated by about .. One may also measure the half-mass or half-light radius alternatively by finding the radius that encloses half of the mass or light within some larger aperture. This is close to the commonly used algorithms in observations for size measurements, such as SExtractor and galfit (e.g. oesch.2010:hudf.morphology), where concentric circular or ellipsoid apertures are usually assumed. However, this approach suffers from two main issues when applying to clumpy, irregular galaxies in our simulations. First, these galaxies do not have a well-defined center: one may use the position of intensity peak on the image or intensity-weighted center and get very different results. Second, for multi-clump systems (e.g. galaxies B, D, E, and F in the rest-frame UV, see Figure 2.1), the size defined in this way in fact represents the spatial separation between the bright clumps. The non-parametric definition we use better reflects the physical size of individual clumps. For single-component objects, such as galaxy A in Figure 2.1 and well-ordered massive galaxies at intermediate and low redshifts, both definitions should give consistent results. Nonetheless, we note that our size measurement depends on how we smooth the star particles. In general, using a larger smoothing length results in slightly larger galaxy sizes, but the difference is usually small for most of the galaxies. We refer the readers to Appendix LABEL:sec:append for details.