COSMOLOGICAL ZOOM SIMULATIONS OF z = 2 GALAXIES: THE IMPACT OF GALACTIC OUTFLOWS

Our simulations were run with an extended version of the N-body + smoothed particle hydrodynamics (SPH) cosmological galaxy formation code Gadget-2 (Springel 2005). Gadget-2 combines an entropy-conserving formulation of SPH (Springel & Hernquist 2002) along with a tree-particle-mesh algorithm for computing gravitational forces. Additions to the public version of the code have been described in Oppenheimer & Davé (2008) and include models for gas cooling, star formation, chemical enrichment, and galactic outflows.

We include photoionization heating starting at z = 9 via a spatially uniform, optically thin UV background (Haardt & Madau 2001). We account for radiative cooling from primordial gas assuming ionization equilibrium as in Katz et al. (1996) and metal-line cooling using the collisional ionization equilibrium tables of Sutherland & Dopita (1993). We track the production of four metal species (C, O, Si, and Fe) from Type II SNe, Type Ia SNe, and asymptotic giant branch (AGB) stars using metallicity-dependent yields as described in Oppenheimer & Davé (2008). We also account for energy feedback from Type II and Type Ia SNe, and mass loss from AGB stars.

Star formation is modeled using the sub-grid prescription of Springel & Hernquist (2003a). Gas particles that are sufficiently dense to become Jeans unstable (n > 0.13 cm⁻³) are treated as a two-phase interstellar medium (ISM) consisting of hot gas that condenses into cold star-forming clouds via a thermal instability (McKee & Ostriker 1977). Stars form from the cold phase and thermal feedback from Type II SNe causes the evaporation of cold clumps into the hot medium. Star formation is implemented probabilistically, such that at any given step, a sufficiently dense gas particle can spawn a star based on a Schmidt (1959) law. The star particle's mass is half the original gas particle mass. The resulting SFRs are tuned to be in accord with the observed Kennicutt (1998a) relation.

Our simulations include a galactic outflow mechanism that imparts kinetic energy to gas particles which is identical to that used in larger-scale simulations of, e.g., Davé et al. (2011b). The outflow is described by two parameters, the wind speed v_w and the mass loading factor η, which is the mass outflow rate in units of the SFR. If a particle is eligible to form stars, it is likewise eligible to be kicked into an outflow, with a probability given by η times the star formation probability. If selected, the gas particle is kicked with a velocity v_w in the direction v × a, where v and a are the particle's instantaneous velocity and acceleration. Hydrodynamic forces are turned off until the particle reaches a density of 0.013 cm⁻³, or a time of 20 kpc/v_w has passed (Oppenheimer & Davé 2008). Decoupling from hydrodynamics allows outflowing gas to escape from the galactic ISM and travel to large distances as observed in z ≈ 2 galaxies (Steidel et al. 2010), a process not properly captured at the resolution of our simulations, and also yields results that are less sensitive to numerical resolution (Springel & Hernquist 2003b). Similar implementations of galactic outflows have been used in recent zoom simulations (Genel et al. 2012b) as well as full cosmological simulations (Barai et al. 2013; Puchwein & Springel 2013).

Choices for v_w and η define our "wind model." Throughout this paper we compare the following wind models:

• 1.  
  No wind model (nw): We turn-off galactic winds (η = 0).
• 2.  
  Constant wind model (cw): v_w = 680 km s⁻¹ and η = 2 for all galaxies.
• 3.  
  Momentum-driven wind model (vzw): The kick velocity scales with galactic velocity dispersion, σ, and the mass loading factor scales as 1/σ, as in the momentum-conserving case (Murray et al. 2005). Following Oppenheimer & Davé (2008) we take $v_{\rm w} = 3\sigma \, \sqrt{f_{\rm L} -1}$ and η = σ₀/σ, where f_L is the (metallicity dependent) critical luminosity necessary to expel gas from the galaxy and σ₀ = 150 km s⁻¹. We run an on-the-fly galaxy finder to calculate galaxy masses that are then converted to σ using standard relations (Mo et al. 1998).

Throughout this paper we assume a Chabrier (2003) initial mass function (IMF) and a ΛCDM concordance cosmology with parameters Ω_Λ = 0.72, Ω_M = 0.28, Ω_b = 0.046, h = 0.7, σ₈ = 0.82, and n = 0.96, consistent with five-year Wilkinson Microwave Anisotropy Probe (WMAP) data combined with baryon acoustic oscillations and SNe constraints (Komatsu et al. 2009) as well as the final nine-year WMAP data (Hinshaw et al. 2013).

Current high-resolution near-infrared observations being carried out by adaptive optics-assisted 10 m class telescopes as well as the Hubble Space Telescope are able to resolve ∼1 kpc scales at z ∼ 2 (e.g., Förster Schreiber et al. 2011; Kartaltepe et al. 2012). In order to get full advantage of these observations, we used the "zoom-in" technique (e.g., Navarro & White 1994) to carry out cosmological simulations that follow the evolution of galaxies down to z ∼ 2 with sub-kiloparsec spatial resolution. The simulations presented here have been used in Anglés-Alcázar et al. (2013) in a different context, where we constrain the growth of massive black holes at the centers of galaxies on cosmological timescales and discuss the implications of the observed black hole–galaxy correlations.

We selected two ∼[5 h⁻¹ Mpc]³ regions for re-simulation from an intermediate-resolution full cosmological simulation surrounding two central galaxies characterized by similar stellar masses (∼3 × 10¹⁰ M_☉) but different morphologies and merger histories. The parent simulation had 2 × 256³ gas+dark matter particles in a [24 h⁻¹ Mpc]³ box and was run including momentum-driven winds. As described in Anglés-Alcázar et al. (2013), zoom initial conditions were generated by identifying all dark matter particles within the virial radius of each selected z = 2 galaxy, tracing their positions back to their locations on the initial grid. Refinement regions were initially defined by the grid cells containing these particles and, subsequently, significantly enlarged by an iterative cleaning procedure incorporating an increasing number of neighboring cells. High-resolution regions were populated with a large number of lower mass particles (yielding × 64 mass resolution increase for our highest resolution simulations) and small-scale power spectrum fluctuations were applied according to the spatial resolution of the refined grid. Additionally, two nested concentric layers of progressively lower resolution were defined surrounding the high-resolution regions in order to reduce numerical artifacts due to the difference in particle masses and to ensure that the large-scale gravitational torques acting on re-simulated halos are accurately represented.

The resulting zoom simulations have (high-resolution) gas particle mass m_gas ≈ 2.3 × 10⁵ M_☉, dark matter particle mass m_DM ≈ 1.2 × 10⁶ M_☉, and softening length ≈ 0.47 h⁻¹ kpc comoving (∼224 pc physical at z = 2), equivalent to 2 × 1024³ particles homogeneously distributed in a [24 h⁻¹ Mpc]³ box. Additionally, we run similar zoom simulations with a factor of 2 lower spatial resolution and a factor of eight lower mass resolution in order to test our results for numerical convergence. A total of 215 snapshot files logarithmically spaced in the redshift range z = 2–11 were produced for each simulation, corresponding to time intervals ranging from ∼5 to 25 Myr.

Galaxies were identified by means of the Spline Kernel Interpolative Denmax algorithm (skid⁸) at all available snapshots independently, and are thus defined as bound groups of star-forming gas particles (i.e., gas particles with densities above the threshold for star formation) and star particles (see Kereš et al. 2005). We used a spherical overdensity algorithm to associate each skid-identified galaxy with a dark matter halo, where the virial radius was defined to enclose a mean density given by Kitayama & Suto (1996). Overlapping halos were grouped together so that every final halo has a central galaxy and a number of satellite galaxies by construction.

The sample of galaxies presented here was defined by selecting the eight most massive central galaxies located within the re-simulated volumes, with the additional requirement of having no contamination by low-resolution dark matter particles within their virial radius; several central galaxies were located near the boundaries of the high-resolution region and were rejected because of contamination. Even though this sample is not mass nor volume complete, our galaxies are representative of "normal" z = 2 systems, with similar gas fractions and SFRs relative to similar mass galaxies from the parent population.

The full evolution of each z = 2 galaxy was reconstructed by identifying its most massive progenitor at all available redshift snapshots. Galaxies were also identified across simulations with different wind models to allow for a detailed model comparison. Following Anglés-Alcázar et al. (2013), we calculated structural and kinematic properties of each galaxy relative to the position of its most bound gas particle, which is a more meaningful and more stable definition of the nominal center of the galaxy compared to that computed by skid, especially during close galaxy encounters and galaxy mergers.

Each galaxy is identified not only across simulations with different wind models but also back in time at all available snapshots, as the central galaxy in the resimulated halo. Figure 3 shows examples of the time evolution of galaxies in terms of their stellar mass, gas fraction, SFR, and metallicity. Since we are interested in the evolution of galaxies over cosmological timescales, the time evolution of all physical quantities has been smoothed over time intervals of ∼50 Myr. Major mergers can still be identified as abrupt changes in the stellar mass of galaxies (e.g., galaxy g222 between redshifts z = 2–3 and g54 at z ≈ 4). These events are generally followed by a temporary increase in the SFR of galaxies but their overall evolution is dominated by smooth accretion, as has been reported for simulations of high-redshift galaxies elsewhere (e.g., Dekel et al. 2009; Davé et al. 2010).

Zoom In Zoom Out Reset image size

Simulations with no winds clearly result in higher stellar masses for all galaxies at all times. Very early on, the SFR of galaxies increases very rapidly for simulations with no winds. When galactic outflows are included, galaxies regulate themselves resulting in significantly lower SFRs relative to the no wind model at z > 3. For momentum-driven winds, recycling of the outflowing gas back into galaxies occurs on a timescale less than the Hubble time, providing an additional gas supply which is not available in simulations with no winds (Oppenheimer & Davé 2010). This additional gas supply results in comparatively higher SFRs toward the end of the simulation at z = 2. Constant winds are, however, more efficient in removing cold gas from galaxies owing to the higher outflow velocities (v_w = 680 km s⁻¹) for comparatively similar mass loading factors (η = 2)—for the three most massive systems in our galaxy sample (g222, g2403, and g54), we find typical wind velocities v_w ≈ 500–600 km s⁻¹ and mass loadings η ≈ 1.4–1.9 in the redshift range z = 8 → 2 for the momentum conserving model. This results in systematically lower SFRs for the constant wind model at all times.

Interestingly, specific SFRs follow a common trend regardless of the wind model adopted: galaxies have the highest specific SFRs very early on (as high as 10 Gyr⁻¹ at z = 8) and then decrease monotonically with decreasing redshift. To first order, outflows reduce SFRs by a factor 1/(1 + η) relative to halo accretion rates (Davé et al. 2012), which leads to specific SFRs that remain the same across wind models (SFRs and stellar masses are reduced by the same amount). Wind recycling accretion, however, yields increasing specific SFRs at later times, especially in the momentum-driven wind model.

Gas fractions are systematically lower for the no wind simulations compared to either the constant wind or the momentum-driven wind models. While galaxies g222, g2403, and g54 have gas fractions below 20% for simulations with no winds at z = 2, galactic outflows are able to maintain gas fractions a factor of two higher at the end of the simulation, in better agreement with observations (Tacconi et al. 2010). Gas fractions may be conveniently expressed in terms of the depletion time t_dep ≡ M_gas/SFR and the specific SFR as f_gas ≡ [(1 + (t_depsSFR)⁻¹)]⁻¹ (Davé et al. 2012). Since specific SFRs are roughly insensitive to wind model early on and the depletion time is shorter in higher stellar mass galaxies ($t_{\rm dep} \propto M_{*}^{-0.3}$; Davé et al. 2011a), no wind simulations yield lower gas fractions.

Finally, metallicities are significantly higher for all galaxies in the no wind model, owing to a greater conversion of baryons into stars and a lack of ejection of metals into the IGM (Hirschmann et al. 2013). These results are qualitatively consistent with those obtained using lower-resolution non-zoom cosmological runs (Davé et al. 2011b, 2011a).

The radial profiles of galaxies reveal interesting differences among wind models, as shown in Figure 4 for galaxies g222, g2403, and g54 at z = 2. Face-on projected stellar and gas surface densities, circular velocity, and metallicity as a function of radial distance from the centers of galaxies have been calculated as average values within logarithmically spaced radial bins. The most massive galaxies (g222 and g2403) reach resolved stellar and gas surface densities in the range Σ_star = 10⁴–10⁵ M_☉ pc² and Σ_gas = 10³–10⁴ M_☉ pc² respectively, with differences up to an order of magnitude among wind models. At their centers, stellar surface densities of most galaxies reach values well above 10³ M_☉ pc² despite gravitational forces being softened at scales ∼ ≈ 224 pc (Figure 2).

Zoom In Zoom Out Reset image size

No wind simulations usually result in higher stellar surface densities at all radii compared to the vzw and cw wind models, and often steeper density gradients at the centers of galaxies. The effects of galactic winds are, however, more evident in the distribution of the star-forming gas: nw galaxies contain most of their gas within the inner kiloparsec, reaching gas surface densities which can be factors of a few above the cw and vzw models. Galactic winds yield more extended disks by removing gas from their centers and having it re-accreted at larger scales. Wind recycling, especially in the vzw model, results in significant amounts of gas extending up to a few kiloparsecs away from the galaxy center (see galaxies g222 and g54 in Figure 4). As we show in the next sections, this result is common to most of our simulated galaxies. Galaxy g2403 is one of the few galaxies that present a very compact gas distribution even in the presence of momentum-driven winds, similar to galaxies from simulations with no galactic outflows.

Most of our simulated galaxies are characterized by flat rotation curves extending up to a few kiloparsecs. The no wind model usually produces more massive and more compact galaxies and this results in higher circular velocities that peak at smaller radii compared to the vzw and cw models. Simulations with winds produce more smoothly rising rotation curves, with lower circular velocities at small radii that tend to approach values similar to the nw model at larger scales, where the contribution from the dark matter component becomes dominant (see galaxies g222 and g54 in Figure 4; this is explored in more detail in Section 4.4). Galaxy g2403 is again an atypical case, with very similar rotation curves for the nw and vzw wind models.

Finally, negative radial gas-phase metallicity gradients (higher metallicities at the centers of galaxies) are common in most simulated galaxies, in general agreement with observations (Yuan et al. 2011; Swinbank et al. 2012b; Jones et al. 2013). Nonetheless, inverted metallicity gradients have been observed in z ∼ 2–3 galaxies and attributed to either galaxy interactions or the infall of metal-poor gas into the centers of galaxies (Cresci et al. 2010; Queyrel et al. 2012; Troncoso et al. 2013). In our simulations, abrupt changes in azimuthally averaged metallicities or even positive metallicity gradients may eventually occur for galaxies undergoing mergers (e.g., galaxy g222 at R ≈ 3–4 kpc). Overall, we find that simulations with no winds usually result in higher metallicities at all radii and also show steeper metallicity gradients. Galactic outflows serve to redistribute metal-enriched gas from the central regions of galaxies over larger scales, resulting in less steep metallicity gradients.

In order to understand the global effects of galactic winds it is useful to look at their effects on different galaxy properties at fixed halo mass, since the latter is expected to be roughly insensitive to baryonic processes. In Figure 5, we plot the halo baryonic fraction, the central stellar mass to halo mass fraction, and the SFR as a function of halo mass for all simulated galaxies at z = 2. Here, each vzw galaxy is connected to the corresponding nw and cw galaxies to help identify both the global effects of outflows on the galaxy population as well as the effects on individual galaxies. Total halo masses are very similar for the nw and vzw models at z = 2, except occasionally due to the slightly different timing of merger events and our adopted definition of dark matter halos for each central galaxy (see Section 2). Halo masses are, however, systematically lower for simulations with constant winds, for which gas outflows are able to escape more easily from the halo potential well, all of which generally have escape velocities below the constant wind velocity of v_w = 680 km s⁻¹. Note that momentum-driven wind velocities are below that of the cw model in all but the most massive galaxies, especially at high redshift when galactic velocity dispersions are lower (Oppenheimer & Davé 2008).

Zoom In Zoom Out Reset image size

Figure 5 shows that halo baryonic fractions are only slightly higher for the no wind model (and higher on average than the cosmological baryonic fraction f_b = 0.165) compared to the vzw model: most of the outflowing gas does not escape the halo potential well for the momentum-driven wind model. This is despite the fact that the wind speed typically is comparable to or exceeding the escape velocity (and scales with it as well); but hydrodynamic (i.e., ram pressure) slowing is actually dominant in many cases (Oppenheimer & Davé 2008). In contrast, the high efficiency of constant winds in expelling gas from halos results in systematically lower halo baryonic fractions, a significant suppression of star formation in the central galaxy, and the consequent reduction of stellar masses for all galaxies in our sample.

Figure 5 (middle panel) confirms our earlier expectations: for a given halo mass, central galaxy stellar masses are systematically higher in the absence of galactic outflows. Interestingly, SFRs for the nw and vzw models happen to be similar at z = 2 for most simulated galaxies. However, as we have seen in Figure 3 for galaxies g222, g2403, and g54, momentum-driven winds result in rather different star formation histories, with significant suppression of star formation early on and enhanced activity due to wind recycling at later times (Oppenheimer & Davé 2010). The wind recycling channel does not act on the cw simulations for the range of halo masses probed here, since these halos generally have escape velocities lower than the assumed wind speed.

Figure 6 illustrates the effects of galactic winds in the star formation histories of our simulated galaxies by showing the evolution of the cosmic star formation efficiency (cSFE) in the redshift range z = 2–6. The cSFE is the ratio of the SFR of the central galaxy to the total halo baryonic accretion rate given in Dekel et al. (2009):

$$\begin{equation} \dot{M}_{\rm acc} \simeq 6.6\, f_{\rm b}\, (M_{\rm halo}/10^{12}{M}_{\odot })^{1.15}\, (1+z)^{2.25}\, {M}_{\odot }\,{\rm yr}^{-1}. \end{equation} \tag{ 1 }$$

At z = 6, the average cSFE of galaxies for simulations including outflows (cw and vzw models) is ∼0.15, about a factor of two lower than simulations with no winds. From z = 6 to z = 4, the average cSFE increases rapidly up to ∼0.69 for nw simulations, and then decreases down to ∼0.54 at z = 2. Simulations with constant winds seem to follow a similar tendency, with the average cSFE peaking at z ∼ 4–5 and decreasing at lower redshifts, but with significantly lower values (roughly × 3, as expected from the (1 + η) suppression of star formation). Momentum-driven winds result, however, in an overall increase of cSFEs from z = 6 to z = 2, owing to the recycling of gas that was launched into winds at earlier times. The vzw model causes an effective delay in star formation by ejecting significant amounts of gas from small, early galaxies and having it reaccreted by z ∼ 2. This coincidentally produces average SFRs similar to simulations with no winds at z ≈ 2. Note that trends for individual galaxies mimic the trends for the overall cosmic star formation history in these wind models, with momentum-driven winds generally producing a later peak in cosmic SFR than no-wind or constant-wind cases (Oppenheimer & Davé 2006).

Zoom In Zoom Out Reset image size

Galaxy rotation curves contain substantial information about the structure and kinematics of galaxies. Figure 7 (top panel) shows the rotation curves of all simulated galaxies at z = 2 (dotted lines), color-coded by wind model. Here, we calculate the rotation velocity from the enclosed mass at a given radius, which we call v_circ. Our galaxies are characterized by flat rotation curves extending up to several kpc scales, with asymptotic circular velocities ranging from ∼100 to 260 km s⁻¹ at z = 2 for the vzw model. In Figure 4, we showed for galaxies g222, g2403, and g54 that galactic outflows can have a significant impact on the overall shape of their rotation curves.

Zoom In Zoom Out Reset image size

The thick solid lines in Figure 7 (top panel) show stacked rotation curves for each wind model, where circular velocities have been averaged over all galaxies as a function of radius. Simulations with winds produce galaxies with more gradually rising rotation curves compared to the more centrally peaked rotation curves of galaxies with no winds. The average circular velocities for the nw and vzw models tend to approach similar values at larger radii, where the contribution from the dark matter component becomes dominant. Simulations with constant winds are more efficient in removing gas from halos, resulting in systematically lower asymptotic circular velocities, but their shape is similar to the momentum-driven winds.

The middle panel of Figure 7 shows a more observationally motivated way of calculating the rotation curve. Rather than taking the enclosed mass and assuming full rotational support, here we directly calculate the SFR-weighted azimuthal velocities of gas particles with respect to the total angular momentum of the galaxy, and average them within logarithmically spaced radial bins, which we call v_ϕ. This is analogous to an Hα rotation curve. The rotation curves derived in this way show a significantly more gradual increase in circular velocity with radius in the inner parts, because of the increased dispersion support in the central regions. This is particularly clear for simulations with no winds, which produce more compact galaxies and centrally peaked rotation curves based on the mass profiles.

The bottom panel of Figure 7 shows the ratio of v_ϕ/v_circ. We see that the gas azimuthal velocity is, on average, lower than the inferred circular velocities. The gas reaches rotation velocities comparable to the circular velocity at ∼3 kpc away from the centers of galaxies. In the inner regions, the contribution of random motions to the dynamical equilibrium of galaxies becomes comparable to the support from ordered rotation. The discrepancy between v_ϕ and v_circ at scales <3 kpc occurs for all simulated galaxies and wind models, and it is well resolved in our simulations. Using v_ϕ to infer the mass profile of galaxies could, therefore, lead to significant misestimations (Valenzuela et al. 2007).

Figure 8 shows the build-up of rotation curves with time by comparing stacked rotation curves from z = 6 to z = 2 for our simulations with momentum-driven winds. Solid lines show v_ϕ (for the gas), and dashed lines show v_circ. We find that the average asymptotic circular velocities increase from ∼90 to 170 km s⁻¹ with decreasing redshift, as expected because our galaxies are becoming more massive. For v_circ, the peak of the (average) rotation curve is always at ∼1–2 kpc, perhaps moving inward to lower redshifts. In contrast, the location of the peak of v_ϕ occurs at 3–4 kpc, and does not vary much with redshift. Comparing v_ϕ to v_circ, we see that gas azimuthal velocities are, on average, lower than the inferred circular velocities at all redshifts. The ratio v_ϕ/v_circ evaluated at R ≈ 3 kpc increases from ∼0.48 at z = 6 to ∼0.86 at z = 2, indicating that galaxies become, on average, progressively more rotationally supported with time (and therefore with increasing mass) in the redshift range z = 6 → 2.

Zoom In Zoom Out Reset image size

We begin by considering in detail the morphology of the two galaxies that were specifically chosen for resimulation, g54 and g222, and therefore lie near the center of the zoomed region. Recall that g54 was chosen to be somewhat more isolated and have a fairly quiet merger history, while g222 lived in a denser region with a more violent merger history.

Figure 9 shows two-dimensional projected views of galaxy g54 at four different inclination angles. The top three rows show the stellar surface density, gas surface density, and SFR density distributions resulting from our fiducial simulation including momentum-driven winds. The three bottom rows show mock Hα line intensity, line-of-sight velocity, and velocity dispersion maps of galaxy g54 that mimic the resolution of current near-infrared integral field spectroscopic observations with SINFONI. Simulated SFRs have been converted to Hα luminosity using L_Hα[erg s⁻¹] = 2.3 × 10⁴¹ SFR [M_☉ yr⁻¹] (from Kennicutt 1998b, corrected for a Chabrier (2003) IMF). Then, mock Hα line intensity maps have been obtained by (1) placing the simulated galaxy at z = 2, converting its physical size to its apparent angular size, (2) convolving the obtained Hα flux map with a 017 beam, and (3) matching the pixel size to those of typical observations (005 pix⁻¹). Mock Hα flux-weighted line-of-sight velocity and velocity dispersion maps have been obtained by degrading the spatial resolution in a similar way.

Zoom In Zoom Out Reset image size

With stellar mass ∼2.5 × 10¹⁰ M_☉, gas fraction f_gas ≈ 0.48, and SFR ∼13.6 M_☉ yr⁻¹, galaxy g54 has a prominent two-arm spiral structure extending up to ∼10 kpc scales that could be observable even at high inclination angles (i ⩽ 60°). Star-forming gas and stars generally trace the same structure, though the stars show a prominent bulge that is absent in the gaseous component. Mock Hα line intensity maps show a somewhat compact central nucleus with radius ∼1 kpc and a factor ∼2 increased brightness relative to the spiral arms (depending on inclination angle). Observations using star formation tracers and rest-frame optical observations should, thus, infer similar spiral morphologies. The velocity maps illustrate the effects of the inclination angle on the observed kinematic properties of galaxies. A very smooth velocity gradient along the morphological major axis can be clearly identified for this galaxy at a wide range of inclination angles (i ⩾ 30°). Interestingly, galaxy g54 shows a characteristic "spider diagram" pattern in the velocity iso-contours at intermediate inclination angles (i ≈ 30°–60°), as expected for inclined rotating disks. For nearly edge-on observations, line-of-sight projected peak-to-peak velocities reach values >400 km s⁻¹, consistent with the circular velocity v_circ ≈ 230 km s⁻¹ calculated at the outer edge of the disk (R ∼ 10 kpc). With a relatively flat velocity dispersion map, this galaxy would be identified as a large, thick, rotation-dominated disk. These properties are in general quite similar to, though somewhat scaled down from BX442, the z = 2.18 grand design spiral observed by Law et al. (2012). BX442 has a measured stellar mass ∼6 × 10¹⁰ M_☉, Hα-derived SFR ∼45 M_☉ yr⁻¹, inclination-corrected circular velocity v_circ ≈ 234 km s⁻¹, and very similar size and morphology compared to our simulated galaxy g54.

Figure 10 shows similar two-dimensional maps for galaxy g222, a marginally higher mass galaxy that has undergone more frequent interactions and mergers. Indeed, there is an infalling galaxy along the lower spiral arm, which is prominent in the stellar distribution but much less evident in the gas, possibly because its gas has been stripped during infall. The stellar surface density distribution shows a high concentration of stars at the center of g222 and possibly a weak, extended spiral structure. The gas and SFR distributions reveal high levels of turbulence on this galaxy. The morphology of g222 appears highly disturbed in Hα emission and it does not seem to trace the smooth stellar component as in the case of g54. The velocity maps show signs of ordered rotation in the underlying large-scale gas distribution at high inclination angles (i ⩾ 60°), but the inferred velocity gradients are significantly disturbed. The velocity dispersion map reveals an irregular gas clumpy structure with significant turbulent motions. The Hα shows up to several clumps, depending on the viewing angles, but these clumps are not apparent in the stellar distribution, which suggests these are short-lived gaseous clumps as found in similar simulations by Genel et al. (2012b), and as inferred in observations by Wuyts et al. (2012); we will examine clump properties in more detail in future work. Overall, this simulated galaxy shares various characteristic morphological and kinematic properties with typical clumpy disk galaxies observed at z ∼ 2 (Förster Schreiber et al. 2009, 2011; Genzel et al. 2011; Swinbank et al. 2012a, 2012b).

Zoom In Zoom Out Reset image size

Note, however, that star formation in our simulated galaxies appears to be in general more centrally concentrated than observed z ∼ 2 disk galaxies (Figure 1), for which a central Hα peak is sometimes not present, particularly for less massive galaxies. This might be partially explained by obscuration effects, which will be quantified in future work in order to make a more detailed comparison between observed Hα line intensity maps and the intrinsic SFR surface density distribution of simulated galaxies. One numerical issue is that the pressurized ISM resulting from the star formation prescription (Springel & Hernquist 2003a) may inhibit the formation of off-center clumps by gravitational instability, especially in galaxies undergoing a more quiescent evolution. Despite this, the two example galaxies presented here—with stellar masses M_star > 10¹⁰ M_☉ that are actually comparable to commonly observed z ∼ 2 systems—illustrate the wide diversity in morphologies predicted in these simulations, as well as its dependence on tracer (gas, stars, Hα). Disk structures are usually present, but can range from quiescent "grand-design" spirals to turbulent and clumpy disks.

We now conduct a more quantitative comparison to the structural and dynamical properties of z ∼ 2 galaxies inferred from near-infrared integral field spectroscopic observations obtained with SINFONI at the Very Large Telescope. In particular, we focus on the "Hα sample" of the SINS survey presented in Förster Schreiber et al. (2009), consisting of 62 rest-frame UV/optically selected star-forming galaxies at z ≈ 1.5–2.5, and the SHiZELS survey (Swinbank et al. 2012b), consisting of nine Hα-selected galaxies at z ≈ 0.84–2.23 drawn from the HiZELS near-infrared narrow-band survey (Geach et al. 2008). In Figure 11, we show some key structural and dynamical quantities versus stellar mass, for our simulated galaxies color-coded by wind model, as well as for the SINS (black crosses) and SHiZELS (brown diamonds) samples.

Zoom In Zoom Out Reset image size

The upper left panel of Figure 11 shows the specific SFR as a function of stellar mass for our sample of galaxies at z = 2. Our simulations naturally predict the existence of a correlation between galactic SFRs and stellar masses, M_⋆∝ SFR, fairly insensitive to stellar feedback models, due to the dominance of smooth and steady cold accretion (e.g., Finlator et al. 2006; Davé 2008; Davé et al. 2011b). As expected from models, a tight M_⋆–SFR relation has been observed out to z ∼ 2, with a roughly constant slope close to unity (e.g., Daddi et al. 2007; Elbaz et al. 2007; Noeske et al. 2007). However, the observed normalization of this relations increases more rapidly with redshift from z = 0 → 2 that current cosmological simulations predict, resulting in too low predicted specific SFRs at z ∼ 2 (Daddi et al. 2007; Davé 2008). Our zoom simulations follow this trend—Figure 11 shows that they lie below the M_⋆–SFR correlation of Daddi et al. (2007) by a factor of two to three in specific SFRs (for the vzw model), and by somewhat higher factors for the other wind models. The discrepancy is slightly worse when compared to the SINS galaxies, although the Hα selection could play some role in preferentially picking out high-SFR objects (Förster Schreiber et al. 2009). On the other hand, the difference between our simulations and the SHiZELS galaxies appears to be not as large, but note that their typical redshifts correspond to the low-z end of the redshift distribution of the SINS sample (z ∼ 1.5), below the final redshift reached by our simulated galaxies. This discrepancy between models and data has been much debated; for instance recent Herschel results suggest that Hα-inferred SFRs may slightly overestimate the true (i.e., far-infrared bolometric) SFR by some non-trivial factor (Nordon et al. 2010). Other ideas for explaining this discrepancy invoke variations in the stellar initial mass function (Narayanan & Davé 2012) or modifying the star formation law to account for metal-dependent H₂ formation (Krumholz & Dekel 2012). We will not pursue this here, except to note that our zoom simulations do not alleviate this discrepancy known from our larger cosmological runs.

Simulations with momentum-driven winds result in higher specific SFRs compared to simulations with no winds at z = 2, more in agreement with observations. In the vzw model, the mass loading factor scales as η ∼ 1/σ and outflows tend to suppress early star formation while providing high gas fractions to maintain comparatively high SFRs at later (z ∼ 2) times, as seen in Figures 3 and 6. This suggests that even stronger outflows might be needed in simulations in order to match the z ∼ 2 M_⋆–SFR relation (Davé 2008). This is not trivial, however, since outflows must not be too strong lest they fail to produce enough early metals to enrich the IGM (Oppenheimer & Davé 2006, 2008) and enough photons to reionize the universe (e.g., Finlator et al. 2012).

The middle left panel of Figure 11 shows the z = 2 size–mass relation for all simulated galaxies. The effective radius R_{1/2, sfr} is defined to enclose half of the total SFR and it is calculated as a two-dimensional projected radius averaged over 100 random viewing angles. Here we see that there is a clear separation of galaxies in the M_⋆–R_{1/2, sfr} plane for the different wind models. For a given stellar mass, simulations with no winds produce galaxies with a very compact distribution of gas, with most of the star formation happening within <1 kpc from their centers. Simulations with constant winds populate the low mass and small size region of the plot, and only the more massive galaxies in the sample are able to maintain gas disks with sizes R_{1/2, sfr} > 1 kpc. Momentum-driven winds produce more extended galaxies with sizes of several kiloparsecs, along with a minority of more compact galaxies.

Outflows affect sizes by preferentially ejecting the star-forming gas from the central regions and having it re-accrete over larger scales (Brook et al. 2012), resulting in large, extended star-forming disks that are more in agreement with SINS and SHiZELS data. Nonetheless, it is interesting that compact galaxies can also occur in the vzw case; galaxy g2403 is one such compact galaxy (see Figure 4). From our visualizations,⁹ it appears that these more compact systems generally arise during or shortly after a significant merger event (Wuyts et al. 2010; Bournaud et al. 2011), but we will quantify this more rigorously in the future.

The diversity of sizes is also reflected in the half-mass radii of the stars R_{1/2, *} (Figure 11, middle right). Generally, the sizes are comparable to the half-SFR radius, but in some cases R_{1/2, *} can be <1 kpc. This may have interesting consequences for the progenitors of so-called compact ellipticals, i.e., early-type galaxies at z ∼ 2 that show very high stellar densities. We show as the dashed line a nominal threshold stellar surface density of $M_*/R_{\rm 1/2,*}^{1.5}=10^{10.3}$ M_☉ kpc^−1.5 for compact galaxies from Barro et al. (2013, see their Figure 1). The radius used for these observations is actually the H-band half-light radius, but this should be fairly comparable to the stellar half-mass radius. Even with winds, a small fraction of our galaxies lie above this surface density threshold (i.e., below the line). This suggests that our simulations do, with some frequency, produce galaxies that have sufficient stellar densities to be the progenitors of compact ellipticals. These simulations have no mechanism for quenching star formation, and hence our compact galaxies are still star-forming; it is an additional constraint on models to form compact passive systems that whatever feedback mechanism is responsible for quenching, it operates when the galaxy is in a dense (likely post-merger) state.

Galactic outflows seem to affect the bulge mass fraction of galaxies, as shown in Figure 11, upper right panel. To estimate the stellar bulge mass fraction, we perform a simple bulge-disk kinematic decomposition for the stellar content of all simulated galaxies. We calculate the azimuthal velocity v_ϕ of each star particle with respect to the direction of the total angular momentum of the galaxy as for our rotation curves, and estimate the bulge mass as double the mass of particles moving with v_ϕ < 0. Though this is not analogous to observational determinations of the bulge mass, it accurately characterizes the mass in the spheroidal component of our simulated systems. With winds, the bulge fraction decreases with increasing stellar mass, suggesting that galaxies start out small and dispersion-dominated, and move toward being more ordered disks. The opposite trend is seen for our no wind simulations, tracking more the canonical behavior that galaxies start out as disks and then merge together to form more dispersion-dominated systems (e.g., White & Frenk 1991). SINS observations suggest that smaller galaxies tend to be more dispersion-dominated, qualitatively favoring our wind simulations, although a more careful comparison that accurately mimics how bulge-to-disk ratios are measured in data are needed for a more definitive result.

The lower left panel of Figure 11 shows the gas velocity dispersion, σ, as a function of stellar mass at z = 2. Here σ is calculated as the spatially integrated (SFR-weighted) one-dimensional velocity dispersion calculated within R_{1/2, sfr} and averaged over 100 random orientations. Despite being extended and rotationally supported disks, most of our galaxies are characterized by high velocity dispersions (>30 km s⁻¹). This suggests that turbulent motions are significant even in rotationally supported z ∼ 2 galaxies, as inferred from observations (Law et al. 2009; Förster Schreiber et al. 2009; Wright et al. 2009; Swinbank et al. 2012b) and reported in previous simulations (e.g., Ceverino et al. 2010; Genel et al. 2012b). Our simulated galaxies are characterized by high σ values but still lower relative to observed z ∼ 2 galaxies with similar stellar masses, suggesting that turbulent motions are not fully resolved at scales comparable to the spatial resolution of our simulations. Indeed, the self-regulated multi-phase model for star-forming gas is meant to capture the turbulent pressure arising from the continuous formation and disruption of gas clouds at a sub-grid level (Springel & Hernquist 2003a). For our simulated galaxies, the effective sound speed of the star-forming gas may reach values ∼150 km s⁻¹ (Anglés-Alcázar et al. 2013), well above the resolved large-scale turbulent velocities.

For a given stellar mass, simulations including outflows result in galaxies with higher gas velocity dispersion, in better agreement with observations. Note that since outflowing gas is hydrodynamically decoupled as it is ejected from the ISM, the outflows themselves are not directly injecting turbulence. Instead, this increased turbulence relative to simulations without winds likely owes to the higher gas and disk fractions, implying lower Toomre Q parameters for similar mass galaxies (Toomre 1964) and hence increased gravitational fragmentation, and possibly to the injection of energy due to the recycling of the outflowing gas back into the galaxies, in qualitative agreement with the analytic equilibrium model of Genel et al. (2012a).

The lower right panel of Figure 11 shows the disk rotation velocity, v_d divided by velocity dispersion σ, as a function of stellar mass. Here v_d is taken as the peak azimuthal velocity (v_ϕ) from the gas rotation curves obtained in Section 4.4. Direct comparisons with data are not straightforward since rotation is computed from a variety of means in the data, and we have not tried to mimic this in detail; hence we have not plotted data here. Nonetheless, the momentum-driven wind model yields a typical v_d/σ ∼ 3, which is in broad agreement with the turbulent high-z disks seen in SINS and SHiZELS. Constant and no wind models result in slightly higher v_d/σ on average, though still within the range of the data. There is little trend with mass, much less than for the stellar bulge fraction, so even though the stellar components of higher mass galaxies are diskier (for simulations with outflows), their gas content is not more rotationally supported with increasing stellar mass. Note that the stellar component is less rotationally supported than the gas in all cases (cooling of shock-heated gas may dissipate turbulent energy) and, therefore, the gas component can be rotationally supported (v > σ) even for stellar bulge-dominated galaxies. Galaxy g2743 represents an extreme example for the simulation with constant winds, with its stellar bulge fraction close to unity and still a rotationally supported gas disk (see Figures 1 and 2).

Overall, the sizes and dynamical properties of simulated disks at z = 2 are in fair agreement with observations from the latest integral field unit studies of high-z star-forming galaxies, particularly in the case of momentum-driven winds. This occurs despite some overly simple assumptions in the modeling, such as decoupling of wind material escaping the disk, and the usual concerns about the ability of SPH to suppress viscosity and resolve dynamical instabilities (e.g., Agertz et al. 2007). Indeed, the properties of simulated galaxies are more strongly dependent on feedback models relative to the details of the specific hydrodynamic technique (e.g., Hopkins et al. 2013). This suggests that our simulations capture the dominant processes for establishing the structural properties of galaxies and provide a plausible model for the formation of disks during this cosmic epoch.

Using our plausible simulated population of z = 2 disk galaxies, we now examine the evolution of physical properties from z = 8 → 2. We focus here on the momentum-driven wind model, since from a variety of wider constraints it is our favored model of the three presented here (see, e.g., Davé et al. 2011b, 2011a). However, many of the results are broadly similar for our other wind models. We show the evolutionary tracks for our 8 galaxies in various colors in Figure 12. Since we are interested on the evolution of galaxies on cosmological timescales, all physical quantities have been averaged over time intervals of ∼150 Myr. All tracks go from left to right (i.e., lower to higher mass), with individual unit redshifts indicated by the points along the tracks starting at z = 8. The same data as shown in Figure 11 is reproduced here, but note that this is for observed z ≈ 2 galaxies and hence it is shown here for reference.

Zoom In Zoom Out Reset image size

The top left panel of Figure 12 shows the evolutionary tracks of simulated galaxies in the main sequence (M_*–sSFR) plane. All galaxies show a similar evolution of decreasing specific SFR as they increase their stellar mass, consistent with a roughly linear SFR to M_* relation with the overall normalization decreasing with redshift. Despite the increase in stellar mass, mergers cause a temporary enhancement of specific SFRs relative to the dominant decreasing trend. These variations in specific SFR are apparent in the evolutionary tracks of our simulated galaxies even after averaging over time intervals of ∼150 Myr. For the three most massive galaxies—g222 (blue), g2403 (green), and g54 (red)—major mergers can be identified in Figure 3 as abrupt changes in their stellar mass, and connected to the effects on the evolutionary tracks in Figure 12 by direct comparison to the location of the points indicating integer redshifts.

At z = 2, galaxy sizes scale with stellar mass for simulations with winds, in agreement with observations (Figure 11). Middle panels in Figure 12 show that galaxies tend to increase in size with time, as expected, but their evolutionary tracks exhibit significant variation between galaxies, resulting in a large scatter in the size–mass diagram at any given redshift. Large variations of the half-SFR radius, R_{1/2, sfr}, tend to occur in redshift intervals during which a major galaxy merger is taking place. This suggests that major mergers, despite representing only a fraction of the total mass growth in galaxies (e.g., Kereš et al. 2005; Rodighiero et al. 2011), have a significant impact on galaxy sizes, as it appears from our visualizations. The evolution of the stellar half-mass radius, R_{1/2, *}, roughly follows that of the star-forming gas and is also imprinted by significant variations occurring during galaxy mergers. Interestingly, our most compact galaxy at z = 2 (g2403; green evolutionary track) was also in the compact regime at z ∼ 4, both cases preceded by a major merger (see Figure 3, middle panel). This suggests that major mergers may drive galaxies toward the region of the size–mass diagram populated by compact ellipticals. At the earliest epochs, the radius of galaxies does not change much, and may in fact become smaller with time. This partly reflects our star formation criterion of n_H > 0.13 cm⁻³ that results in star formation occurring over a wide area when the universe is very dense; it is therefore unlikely to be a robust prediction.

The upper right panel of Figure 12 shows the stellar bulge fraction evolution. Early on, most galaxies are bulge-dominated, as a result of high merger rates and large gas reservoirs that keep the mass distribution disordered. As time proceeds, the combined effects of smooth gas accreting into galaxies with higher specific angular momentum and galactic outflows removing preferentially low angular momentum gas from their centers, cause galaxies to increase their sizes and reduce their stellar bulge fractions on average (Governato et al. 2009, 2010; Brook et al. 2011). It is interesting that this process only kicks in around z ∼ 4 to start producing disk-dominated systems, at least for the range of galaxy masses considered here. Our simulations, thus, predict that z ∼ 4 was the beginning of the epoch of disk formation for massive galaxies. By z ∼ 2, this results in a trend of decreasing bulge fraction with increasing stellar mass in this wind model. We note that this trend is opposite to simple expectations from classic hierarchical galaxy formation models (e.g., White & Frenk 1991), in which disks form first and then merge later to give rise to more dispersion-dominated systems.

The bottom left panel of Figure 12 shows how simulated galaxies evolve with redshift in the M_*–σ plane. Galaxies tend to evolve along the correlation of σ and M_* observed at z ∼ 2, extrapolated to lower masses. Disk peak rotation velocities follow a similar trend with stellar mass, giving rise to a nearly constant ratio v_d/σ ≈ 3 for all simulated galaxies independent of redshift, as shown in Figure 12, bottom right panel. Hence while the predicted v_d/σ values are similar to that observed at z = 2, it remains to be seen whether these same simulations can produce very thin disks as seen at z = 0. In order to do so, something must alter the current evolution of v_d/σ, perhaps as a result of the dropping accretion rate or else the dropping outflow rate. Note that the simulated velocity dispersions shown here have been averaged over many random directions, in analogy with the spatially integrated σ values uncorrected for average background velocity gradients inferred for the SINS galaxies (Förster Schreiber et al. 2009). Contributions from disk rotation velocities may thus overestimate σ values and underestimate the inferred rotational support of galaxies. For simulated galaxies, we can eliminate the contributions from ordered rotation simply by calculating the velocity dispersion along the line of sight perpendicular to the plane of the disk, σ_z. Interestingly, while we find σ_z < σ in most cases, as expected, σ_z is also correlated with the stellar mass of galaxies.

Overall, the redshift evolution of these properties shows a consistent buildup of size and velocity dispersion with time (and mass), although mergers can particularly impact the inferred sizes substantially over short time periods. The large scatter in observed sizes may, therefore, reflect the short-term merger history of galaxies even more reliably than instantaneous SFRs or bulge-to-disk ratios. Galaxies evolve from being bulge-dominated when small to disk-dominated when larger, with disks becoming prominent only at z ≲ 4. Despite this, galaxies are always rotationally supported even at early times when bulge-dominated, suggesting that the present-day association between small bulge-to-disk ratio and large rotational support does not necessarily apply to high-redshift galaxies.