FORMING REALISTIC LATE-TYPE SPIRALS IN A ΛCDM UNIVERSE: THE ERIS SIMULATION

The formation of realistic late-type spirals has been a longstanding problem of galaxy formation in a ΛCDM universe. Within this framework, baryons condense at the center of dark matter halos and acquire angular momentum through tidal torques from nearby structures (Fall & Efstathiou 1980). A centrifugally supported baryonic disk forms, with a size that depends on the fraction of the original angular momentum that is retained during the contraction. In numerical simulations of this process, however, a fundamental "angular momentum problem" arises, as galaxies are produced with a baryonic component that is quite deficient in angular momentum compared to real spirals (Navarro & Benz 1991; Navarro & Steinmetz 2000).

Aside from artificial losses of angular momentum caused by insufficient resolution and other numerical effects (Okamoto et al. 2003; Governato et al. 2004; Kaufmann et al. 2007), this failure has traditionally been traced back to the very nature of the hierarchical buildup of structures: dynamical friction transfers the orbital angular momentum of merging substructures to the outer halo, and causes the associated cold baryons to sink to the center of the proto-galaxy and form a spheroid rather than a disk (e.g., Maller & Dekel 2002).

A popular solution to the angular momentum problem envisions energy injection from supernovae (SNe) and evolving stars as a mechanism to prevent efficient gas cooling and condensation and to remove low-angular-momentum material from the central part of galaxies. Modern simulations with improved resolution and more effective recipes for SN feedback (Robertson et al. 2004; Governato et al. 2007; Scannapieco et al. 2009; Stinson et al. 2010; Piontek & Steinmetz 2011; Brooks et al. 2011) have yielded rotationally supported disks with realistic exponential scale lengths, not only in galaxies formed in relative isolation but also in those that are accreted by massive groups with a dominant central elliptical (Feldmann et al. 2010, 2011). They have also modified the standard picture of gas accretion and cooling onto galaxy disks: for galaxies up to Milky Way masses, gas acquired through filamentary "cold flows" that was never shock-heated to the halo virial temperature is largely responsible for star formation in the disk at all times (Brooks et al. 2009; Kereš et al. 2009; Ceverino et al. 2010). Yet, these simulations typically continue to produce centrally concentrated systems, with rotation curves that rise steeply toward the center: simulated disk galaxies fall exclusively in the S0 or Sa category, leaving late-type spirals with negligible bulges and large disks with flat rotation curves—such as our own Milky Way—as an unsolved puzzle. Two recent exceptions are the simulations of Agertz et al. (2011) and Governato et al. (2010). In the first, replicas of Sb/Sc galaxies with moderate bulges were obtained with a low efficiency of star formation that may implicitly mimic the bottleneck of the conversion of atomic gas into molecular, at the expense of producing stellar disks that are much more massive than expected at a given halo mass (e.g., Guo et al. 2010). In the second, a realistic bulgeless dwarf galaxy with a shallow central dark matter profile was generated by resolving the inhomogeneous interstellar medium (ISM) and the process of energy injection from multiple SNe in clustered star-forming regions. In this paper we extend the latter approach to massive galaxy scales, and present initial results from a new smooth particle hydrodynamic (SPH) cosmological simulation of high dynamic range that includes radiative cooling, heating from a cosmic UV field, SN feedback, and a star formation recipe based on a high gas density threshold as in Governato et al. (2010). It is this last feature, we argue, that is key to the formation of a more realistic massive late-type spiral in ΛCDM.

Dubbed "Eris," the simulation described in this paper is part of a campaign of extreme resolution simulations of the formation of Milky Way sized galaxies (Diemand et al. 2007, 2008). It was performed in a Wilkinson Microwave Anisotropy Probe three-year cosmology, Ω_M = 0.24, $\Omega _\Lambda =1-\Omega _M$, Ω_b = 0.042, H₀ = 73  km s⁻¹ Mpc⁻¹, n = 0.96, σ₈ = 0.76, running the parallel, spatially and temporally adaptive, treeSPH-code GASOLINE (Wadsley et al. 2004) for 1.5 million cpu hours. The target halo was identified at z = 0 in a low-resolution, dark-matter-only, periodic box of 90 Mpc on a side. It was chosen to have a similar mass as the Milky Way and a rather quiet late merging history, i.e., to have had no major mergers (defined as mass ratio ⩾1/10) after z = 3. New initial conditions were then generated with improved mass resolution, centered around a Lagrangian subregion of 1 Mpc on a side, using the standard "zoom-in" technique to add small-scale perturbations. High-resolution particles were further split into 13 million dark matter particles and an equal number of gas particles, for a final dark and gas particle mass of m_DM = 9.8 × 10⁴  M_☉ and m_SPH = 2 × 10⁴  M_☉, respectively. The gravitational softening length, _G, was fixed to 120 physical pc for all particle species from z = 9 to the present, and evolved as 1/(1 + z) from z = 9 to the starting redshift of z = 90.

The version of the code used in this study includes Compton cooling, atomic cooling, and metallicity-dependent radiative cooling at low temperatures (Mashchenko et al. 2006). A uniform UV background modifies the ionization and excitation state of the gas and is implemented using a modified version of the Haardt & Madau (1996) spectrum. Three parameters characterize the star formation and feedback recipes: (1) the star formation threshold n_SF, (2) the star formation efficiency _SF, and (3) the fraction of SN energy that couples to the ISM _SN. Star formation occurs when cold (T < 3 × 10⁴ K), virialized gas reaches a threshold density n_SF = 5 atoms cm⁻³ and is part of a converging flow. It proceeds at a rate

$$\begin{equation} d\rho _*/dt=\epsilon _{\rm SF} \rho _{\rm gas}/t_{\rm dyn} \propto \rho _{\rm gas}^{1.5} \end{equation} \tag{ 1 }$$

(i.e., locally enforcing a Schmidt law), where ρ_* and ρ_gas are the stellar and gas densities, and t_dyn is the local dynamical time. We choose _SF = 0.1. An additional identical run with _SF = 0.05, the value adopted in Governato et al. (2007) and Brook et al. (2011), yielded a galaxy with nearly identical structural properties, and will be discussed in a forthcoming paper. Each star particle is created stochastically with an initial mass m_* = 6 × 10³  M_☉, and the gas particle that spawns the new star has its own mass reduced accordingly. A star particle represents a simple stellar population with its own age, metallicity, and a Kroupa (2001) initial stellar mass function (IMF). Each SN deposits metals and a net energy of _SN × 10⁵¹ erg into the nearest neighbor gas particles, with _SN = 0.8 (the same value adopted in previous simulations). The heated gas has its cooling shut off until the end of the snowplow phase of the SN blastwave, which is set by the local gas density and temperature and by the total amount of energy injected E (Stinson et al. 2006). For the typical ISM conditions at the threshold found in this study, this translates into regions of size R_E ∼ 30E^0.32₅₁ pc heated by individual SNe and having their cooling shut off for a timescale t_E ∼ 5 × 10⁵E^0.31₅₁ yr, where E₅₁ ≡ E/10⁵¹ erg. The energy injected by many SNe adds up to create larger hot bubbles and longer shutoff times. No feedback from an active galactic nucleus was included.

The adoption of a density threshold for star formation that is 50 times higher than in many previous lower-resolution studies is possible owing to the high mass and spatial resolution of this run, which resolves the giant cloud complexes where star formation actually occurs in the ISM and the true scale height of the neutral atomic ISM. In particular, the local Jeans length corresponding to our density threshold (for T = 10³ K, a lower bound on the typical temperature of the cold gas in the simulations) is resolved with more than five SPH smoothing lengths, thus preventing artificial fragmentation (Bate & Burkert 1997). While not as high as the value of n_SF = 100 atoms cm⁻³ used in the Governato et al. (2010) dwarf galaxy simulation, whose particle mass was substantially lower and allowed the resolution of the Jeans length of star-forming gas at much higher densities, our star formation threshold is still large enough to allow the development of a clumpy, inhomogeneous ISM with more localized energy injection by multiple overlapping SN explosions. This allows galactic pressure-driven outflows to develop and remove low-angular-momentum material. To demonstrate the important role of the star formation threshold on the structural properties of massive galaxies, we have run a low-threshold twin simulation (termed "ErisLT") with n_SF = 0.1 atoms cm⁻³. We have kept all the other simulation parameters fixed (same mass and spatial resolution and identical feedback scheme) except for the star formation efficiency parameter, _SF, which was lowered from 0.1 (Eris) to 0.05 (ErisLT) to match the observed normalization of the star formation density in local galaxies (see Governato et al. 2010). ErisLT was stopped at redshift 0.7 in order to limit the computational burden.

In this first paper, we focus on the final z = 0 state and properties of the galaxy and on the comparison with observational constraints. The main characteristics of the simulated galaxy are listed in Table 1.

Table 1. Properties of the Simulated Galaxy

Galaxy	Mvir	Vpeak	M*	fb	fcold	mDM	mSPH	G	SF	NDM	Ngas	N*	B/D	Rd
Eris (z = 0)	7.9	238	3.9	0.121	0.12	9.8	2	120	0.1	7.0	3.0	8.6	0.35	2.5
Eris (z = 1)	5.4	237	2.9	0.126	0.40	9.8	2	120	0.1	4.8	2.0	6.2	0.30	1.8
ErisLT (z = 1)	5.5	308	3.4	0.158	0.18	9.8	2	120	0.05	4.9	2.9	8.3	0.42	1.4

Notes. Columns 2, 3, 4, 5, and 6 list the virial mass (in units of 10¹¹  M_☉), peak circular velocity (in  km s⁻¹), total stellar mass of the halo (in units of 10¹⁰  M_☉), baryonic mass fraction, and cold (T < 3 × 10⁴ K) gas fraction. Columns 7 and 8 list the mass resolution of individual dark matter and SPH particles (in units of 10⁴  M_☉), and Columns 9 and 10 list the spline gravitational force softening (in pc) and the star formation efficiency. Columns 11, 12, and 13 list the total number (in units of 10⁶) of dark matter, gas, and star particles within the virial radius of the halo. Columns 14 and 15 list the bulge-to-disk ratio and disk scale length (in kpc) estimated from the i-band photometric decomposition.

Download table as:  ASCIITypeset image

3.1. Structural Parameters

At the present epoch, Eris is a massive, barred, late-type spiral with structural properties consistent with those of Sb/Sbc galaxies, of which our own Milky Way is an example. It has a virial radius of R_vir = 239 kpc (defined as the radius enclosing a mean density of 93ρ_crit; Bryan & Norman 1998), total mass M_vir = 7.9 × 10¹¹  M_☉, spin parameter λ = 0.019 (à la Bullock et al. 2001), and 7.0, 3.0, and 8.6 million dark matter, gas, and star particles within R_vir, respectively. The minimum smoothing length for gas particles is five times smaller than the force softening. The total mass enclosed within 60 kpc is M_<60 = 3.3 × 10¹¹  M_☉. The rotation curve, shown in the left panel of Figure 1, has a peak circular velocity of V_peak = 238  km s⁻¹ (reached at 1.34 kpc) and a value at 8 kpc (the solar circle) of V_{c, ☉} = 206  km s⁻¹. Its overall shape out to 20 kpc, including the peak in the central bulge-dominated kpc, is reminiscent of the recent reconstruction of the Milky Way rotation curve by Sofue et al. (2009). The circular velocity decreases gently to distances of 60 kpc from its value at the solar radius, in agreement with observations of blue horizontal-branch halo stars in the Sloan Digital Sky Survey (SDSS; Xue et al. 2008). The measured V_{c, ☉}, M_<60, and M_vir agree within the errors with the values of V_{c, ☉} = 221 ± 18  km s⁻¹, M_<60 = (4.0 ± 0.7) × 10¹¹  M_☉, and M_vir = 1.0^+0.3_−0.2 × 10¹²  M_☉ derived recently for the Milky Way using the narrow GD-1 stream of stars (for V_{c, ☉}; Koposov et al. 2010) and halo stars as kinematic tracers (for M_<60 and M_vir; Xue et al. 2008).

Zoom In Zoom Out Reset image size

3.2. Brightness Profile

To correctly compare simulations with observations, we created artificial images of our simulations and from them measured photometric bulge-to-disk ratios and disk scale lengths. The mock images were created using the radiation transfer code Sunrise (Jonsson 2006), which produces spectral energy distributions using the age and metallicities of each simulated star particle, and takes into account the three-dimensional effect of dust reprocessing. The results for a Kroupa IMF are shown in Figure 2. A two-dimensional photometric decomposition was performed on the dust-reddened i-band light distribution with the Galfit program (Peng et al. 2002). At the present epoch, the total i-band magnitude is M_i = −21.7, and a stellar disk with a scale length R_d = 2.5 kpc dominates the light distribution (Figure 3). The disk scale length is comparable to the value R_d = 2.3 ± 0.6 kpc, adopted for the Milky Way in the compilation by Hammer et al. (2007), and with the scale length of the Milky Way thin disk, 2.6 kpc, as traced by M dwarfs in the solar neighborhood (Jurić et al. 2008). Its value also agrees with the scaling relations of spiral galaxies (Courteau et al. 2007). The SDSS u − g = 1.03 mag and g − r = 0.49 mag integrated colors, obtained directly from the Sunrise images, fall within 1σ of the mean optical colors of late-type galaxies as luminous as Eris (Blanton et al. 2003). The "downbending" observed in Eris' brightness exponential profile at about four disk scale lengths appears to be characteristic of late-type spirals (Pohlen & Trujillo 2006). As in the sample of truncated late-type spirals of Bakos et al. (2008), there is no break in the stellar surface mass density profile of Eris: rather, Eris' stellar age profile shows a "U shape" with a minimum of 6 Gyr at the break radius, explaining the origin of the break as a radial change in stellar population likely caused by the stochastic radial migration of young stars from the inner parts of the disk to the outskirts (Ros`kar et al. 2008).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Eris' bulge-to-disk ratio (as determined by a two-component fit to the i-band surface brightness profile), B/D = 0.35, is also typical of Sb spirals, which are characterized by a median (±68/2%) value log B/D = −0.53^+0.27_−0.30, and of many Sbc galaxies, which have log B/D = −0.86^+0.34_−0.40 (Graham & Worley 2008). A three-component decomposition (disk+bar+bulge) will lower the B/D ratio further. The bulge Sérsic index, n_s = 1.4, is indicative of a "pseudobulge" rather than a classical one: according to Weinzirl et al. (2009), ∼3/4 of all bright spirals have low n_s ⩽ 2 bulges. Eris' large final disk (disk-to-total ratio D/T = 0.74) is not typically found in lower-resolution simulations of Milky Way sized galaxies that impose no restrictions on merger history: e.g., only one of the eight galaxies simulated by Scannapieco et al. (2010) has a photometric D/T as large as 0.68 (and six have D/T < 0.5), and only one out of the six galaxies above M_vir = 10¹¹  M_☉ simulated by Brooks et al. (2011) has a disk-to-total ratio comparable to Eris' ("h239," which is offset, however, from the stellar-mass–halo-mass relation).

3.3. Stellar Content

Eris' total mass in baryons is M_b = 9.5 × 10¹⁰  M_☉, corresponding to a mass fraction f_b = 0.12 that is 30% lower than the universal value (for the adopted cosmology) of 0.175. Stars (and their remnants) comprise 41% of all baryons within R_vir: the total stellar mass, M_* = 3.9 × 10¹⁰  M_☉, is comparable to the value estimated for the Milky Way, (4.9–5.5) × 10¹⁰  M_☉, by Flynn et al. (2006).

To make a bias-free comparison with the stellar-mass–halo-mass relation derived from the abundance matching technique by Behroozi et al. (2010) we adopt the following procedure. We fit the SDSS u, g, r, i, z broadband colors from the mock Sunrise images with the flexible stellar population synthesis code of Conroy et al. (2009): the fit assumes a Kroupa IMF and provides a photometric stellar mass estimate of ${\cal M}_*=3.2\times 10^{10}\,\,{M_\odot }$ (C. Conroy 2011, private communication), 18% lower than the value directly measured in the simulation. The photometric stellar mass of Eris can now be weighted self-consistently against the Behroozi et al. (2010) average stellar-mass–halo-mass relation (which uses a Chabrier 2003 IMF), free of IMF systematics, after offsetting all Behroozi et al. (2010) stellar masses by 0.06 dex (to correct from Chabrier to Kroupa IMF). The comparison, depicted in the right panel of Figure 4, demonstrates that Eris' implied "baryon conversion efficiency," $\eta \equiv ({\cal M}_*/M_{\rm vir}) \times (\Omega _M/\Omega _b)=23\%$, is in excellent agreement with that predicted by the abundance matching technique. This contrasts with the recent analysis of many hydrodynamic simulations of galaxy formation by Guo et al. (2010), who show that the great majority of them lock too many baryons into stars to be viable models for the bulk of the observed galaxy population. Note that the intrinsic scatter in the stellar mass at a given halo mass is estimated to be 0.17 dex, independent of halo mass (Yang et al. 2009).

Zoom In Zoom Out Reset image size

With a circular velocity at the radius, R₈₀ = 6.8 kpc, containing 80% of the i-band flux of V₈₀ = 210  km s⁻¹, our Galaxy lies close to the Tully–Fisher relation of the Pizagno et al. (2007) galaxy sample (see the left panel of Figure 4). As discussed in Pizagno et al. (2007), the Tully–Fisher relation uses V₈₀ as the primary velocity measure rather than V_2.2, the circular velocity at 2.2 disk scale lengths, since the former is less sensitive to the degeneracies of bulge–disk decomposition. The ratio V_2.2/V₂₀₀ = 214  km s⁻¹/129  km s⁻¹ = 1.66 in Eris, where V₂₀₀ is the circular velocity at the radius enclosing a mean overdensity of 200 ρ_crit (R₂₀₀ = 177 kpc), is equal to the value suggested by the dynamical model for the Milky Way of Klypin et al. (2002). It is also consistent with the recent measurements of the virial mass of the Milky Way by Smith et al. (2007) and Xue et al. (2008), implying V_2.2/V₂₀₀ = 1.48^+0.25_−0.26 and V_2.2/V₂₀₀ = 1.67^+0.31_−0.24, respectively.⁴ Note that while there is an unsettled disagreement among estimates of the Galaxy's virial mass,⁵ Eris is likely to be 20% less massive than the Milky Way and therefore the agreement between simulated and observed kinematic data should be confirmed in future simulations of slightly more massive halos and different accretion histories. Like the Milky Way, however, Eris is offset relative to determinations using various dark halo mass tracers for late-type disk galaxies by Dutton et al. (2010), who predict the ratio V_2.2/V₂₀₀ = 1.11^+0.22_−0.20(2σ) for typical dark matter halos with V_2.2 = 220  km s⁻¹.

3.4. Gas Content

Eris' H $\mathsc{i}$ mass is $M_{{\rm H}\,\mathsc{i}}=1.9\times 10^9\,\,{M_\odot }$, comparable to the H $\scriptstyle \rm I$ mass estimated for the Milky Way disk by Nakanishi & Sofue (2003), but smaller than the value of ∼5 × 10⁹  M_☉ given by Wolfire et al. (2003). The H $\scriptstyle \rm I$-to-stellar mass ratio, 1.9 × 10⁹  M_☉/3.9 × 10¹⁰  M_☉ = 0.049, is equal to the median value observed in the Galactic All-Sky Survey (GASS; Catinella et al. 2010) for galaxies of comparable stellar mass. Eris' H $\scriptstyle \rm I$ disk extends out to about 15 kpc (six stellar disk scale lengths), similar to the size of the H $\scriptstyle \rm I$ disk of the Milky Way (Nakanishi & Sofue 2003). Clustered SN explosions create a large number of holes in the face-on H $\scriptstyle \rm I$ distribution of Eris (Figure 2) due to bubbles of hot gas expanding perpendicular to the disk. These holes are mostly located within the bright optical disk and preferentially in regions of high star formation: they are kpc in size, as observed, e.g., in the nearby low-inclination spiral galaxy NGC 6946 (Boomsma et al. 2008).

About 6.7 × 10⁹  M_☉ are found in a cold phase below T = 3 × 10⁴ K. This is comparable to the total mass of the atomic and warm ionized medium inferred for the Milky Way (e.g., Ferriére 2001). The gas mass that is hot (T > 3 × 10⁵ K) and thus potentially X-ray luminous is M_X = 3.6 × 10¹⁰  M_☉, 63% of the total gas content. For comparison, 12% of the gas is in the cold phase and 25% is in a warm phase with 3 × 10⁴ K < T < 3 × 10⁵ K. The fractions of cold, warm, and hot gases within its inner 20 kpc are 83.5%, 1.5%, and 15%, respectively. Hot gas within 20 kpc contains significant amount of angular momentum and is co-rotating with the cold disk. The hot gas baryon fraction, f_X = M_X/M_vir = 0.046, is 3.8 times smaller than the cosmological baryon fraction. This implies an average density for the hot gas that is 3.8 times smaller than assumed in the standard "cooling flow" halo model à la White & Frenk (1991), and yields a factor of 14.5 smaller X-ray emission measure. Contrary to the standard assumption that hot gas follows the radial distribution of the dark matter, the density distribution of hot gas in Eris follows a "flattened" ρ_X(r)∝r^−1.13 power-law profile out to 100 kpc (see Figure 5). This gives origin to an X-ray surface brightness profile that is not as sharply peaked as expected for hot halos with Navarro–Frenk–White (Navarro et al. 1996; NFW) profiles and that satisfies the observational constraints (e.g., Anderson & Bregman 2010). The flattened density profile we find is consistent with the results of much lower-resolution simulations by Crain et al. (2010), who identified the reason for the more extended gas distribution and weaker X-ray coronae in the entropy injection by SNe at z ∼ 1–3. Only 10% of the gas in Eris is in a very hot phase above T = 10⁶ K: the mean density of million degree gas at R ⩾ 70 kpc is n ⩽ 6 × 10⁻⁵ atoms cm⁻³, which is well within the observational constraints from O $\scriptstyle \rm VI$ absorption measurements (Sembach et al. 2003) in the halo of the Milky Way, and high enough to produce significant ram pressure stripping of dwarf spheroidal satellites (Mayer et al. 2007).

Zoom In Zoom Out Reset image size

The observed dispersion measure (DM) of pulsars in the Large Magellanic Cloud (LMC) provides another constraint to the hot halo of the Milky Way (Anderson & Bregman 2010). Of the 11 pulsars discovered by Manchester et al. (2006) in the direction of the LMC, 3 have DMs below 45  cm⁻³ pc and may be located within the Galaxy at a random position along the line of sight to the Cloud. The other eight have DMs in the range between 65 and 130  cm⁻³ pc and are thought to be associated with the LMC, with the lowest dispersion values belonging to pulsars located on the near side of the Cloud, some 50 kpc away. This leads to an estimate for the DM introduced by Galactic free electrons toward the LMC of DM ≈ 70  cm⁻³ pc (Anderson & Bregman 2010). To compare these values with prediction from the Eris simulation, we have calculated the mean integrated column density of free electrons between randomly positioned "observers" on a circle in the stellar disk at a galactocentric distance of 8 kpc, and a "pulsar" 50 kpc away. The mock pulsar was located at galactic coordinates l = 280° and b = −30° like the LMC. The predicted DM,

$$\begin{equation} {\rm DM} = \int _0^{50\,{\rm kpc}} n_e(l)dl=(62 \pm 3) \,{\rm cm^{-3} \,pc,} \end{equation} \tag{ 2 }$$

where n_e is the electron density along the line of sight, is perfectly consistent with the data. Collectively, the above arguments indicate that the reservoir of hot gas around the Eris simulated galaxy appears to satisfy pulsar DM as well as X-ray surface brightness and O $\scriptstyle \rm VI$ absorption constraints.

3.5. Star Formation and Kinematic Decomposition

The top left panel of Figure 6 shows the star formation history of all star particles identified within Eris' virial radius at z = 0 (regardless of whether they formed within the main host or in satellites), and of its kinematically decomposed present-day disk and spheroid. The decomposition technique follows Scannapieco et al. (2009), and is based on the distribution of orbital circularity parameters, j_z/j_c, of the simulated stars introduced by Abadi et al. (2003). Here, j_z is the angular momentum of each star in the z-direction (i.e., the direction defined by the total angular momentum of all gas particles within 5 kpc from the host center) and j_c is the angular momentum of a circular orbit at the same radius. Spheroidal stars are defined as those that are not part of the rotationally supported disk and therefore typically include bulge and stellar halo stars, as well as stellar bars if they are present. The distribution of circularity parameters, shown in the top right panel of Figure 6, is characterized by two peaks: one at j_z/j_c ≃ 1 that is indicative of the presence of a dominant cold disk in rotational support, and a second one at j_z/j_c ≃ 0 corresponding to a modest hot spheroidal component dominated by velocity dispersion.

Zoom In Zoom Out Reset image size

The vast majority of disk stars in the Milky Way reside in a thin disk component with exponential scale height h_z = 300 ± 60 pc (Jurić et al. 2008). A study of the vertical structure of edge-on spiral galaxies finds that the scale height of their thin disk increases systematically with circular velocity as z₀ = 610 (V_c/100  km s⁻¹)^0.9 pc (Yoachim & Dalcanton 2006), where z₀ is the scale height of a sech² profile, z₀ ≈ 2h_z at large heights above the disk plane. Eris' kinematically decomposed stellar disk agrees well with the above scaling: by fitting an exponential (sech²) profile to the simulation data, we derive a scale height of h_z = 490 pc (z₀ = 860 pc) at a galactocentric distance of 8 kpc.

Today, Eris is forming stars at a rate of SFR = 1.1  M_☉ yr⁻¹, comparable to the value of SFR = 0.68–1.45  M_☉ yr⁻¹ recently inferred for the Milky Way by Robitaille & Whitney (2010) using Spitzer data. The star formation rate declines rapidly with redshift from a plateau value of ∼10  M_☉ yr⁻¹ maintained between z = 2 and z = 5. SN feedback and photoheating by the ultraviolet radiation background efficiently quench star formation at z > 5. The rate of formation of spheroidal stars fades rapidly after redshift 3, while disk stars nearly triple their mass from z = 2 to the present.

The star formation rate surface densities Σ_SFR and H $\scriptstyle \rm I$ gas surface densities $\Sigma _{{\rm H}\,\mathsc{i}}$ (Kennicutt–Schmidt law) of Eris' disk are shown in the bottom left panel of Figure 6. The simulation data were averaged over square patches 750 pc on the side for comparison with the pixel-to-pixel observations at 750 pc resolution of seven spiral galaxies from the THINGS survey (Bigiel et al. 2008). The same THINGS data are also plotted against the total gas surface density $\Sigma _{\rm gas}=\Sigma _{{\rm H}\,\mathsc{i}}+\Sigma _{\rm H_2}$. The observed relationship between Σ_gas and Σ_SFR varies dramatically among and within spiral galaxies, and most galaxies show little or no correlation between $\Sigma _{{\rm H}\,\mathsc{i}}$ and Σ_SFR (Bigiel et al. 2008). Rather, there is a clear correlation between $\Sigma _{\rm H_2}$ and Σ_SFR (molecular Schmidt law), and gas at densities in excess of ∼10  M_☉ pc⁻² is observed to be fully molecular. At solar metallicities, this corresponds to the column of atomic hydrogen required to shield a molecular region against photodissociation (Krumholz et al. 2009). While the total mass of cold gas in Eris is comparable to the sum of the atomic and molecular gas mass of Sb galaxies such as the Milky Way, we do not model directly the formation of H₂ molecules. Indeed, as clearly shown in Figure 6, even at Eris's high resolution we do not capture kpc-sized regions with gas surface densities in excess of the characteristic shielding column. Furthermore, our gravitational softening is still large compared to the size of typical molecular clouds, which are a few tens of pc in size, and our gas density threshold for star formation is still somewhat below the density of real molecular clouds. Both effects contribute to further smooth out the high-density tail of the total gas distribution. Yet the figure shows that our simulated disk is forming stars in the same range of $\Sigma _{\rm SFR}\hbox{--}\Sigma _{{\rm H}\,\mathsc{i}}$ values observed in spiral galaxies.

The global atomic gas depletion timescale, i.e., the time needed for the present rate of star formation to consume the existing atomic gas reservoir, is $t_{{\rm H}\,\mathsc{i}}\equiv M_{{\rm H}\,\mathsc{i}}/{\rm SFR}=1.9\times 10^9\,\,{M_\odot }/1.1\,\,{M_\odot \,\rm yr^{-1}}=1.7$ Gyr in Eris today. As this is significantly smaller than the Hubble time, star formation rates and gas fractions must be set by the balance between gas accretion from the halo and stellar feedback. The depletion time for atomic gas observed in the COLD GASS survey (Saintonge et al. 2011) is on average 3 Gyr, with a large scatter from one galaxy to another.

Approximately 70% of bright spirals have B/T ⩽ 0.2 (Weinzirl et al. 2009). Recent attempts to generate such disk-dominated galaxies in cosmological simulations have failed to reproduce simultaneously their observed morphologies as well as their baryonic/stellar content (Scannapieco et al. 2010; Agertz et al. 2011). Our very high resolution Eris simulation appears to form a close analog of a Milky Way disk galaxy by capturing a realistic inhomogeneous ISM in which star formation occurs in high-density regions of mass comparable to that of giant cloud complexes. Contrary to the "inefficient star formation" prescription adopted by Agertz et al. (2011), this is achieved with a strong localized SN feedback and a high (10%) Schmidt-law efficiency. We stress that our efficiency parameter is simply phenomenological, and does not have a direct relation to the true star formation efficiency within giant molecular clouds, which is the end result of all the processes regulating star formation including feedback, and concerns scales (tens of parsecs) that are not yet resolved in cosmological simulations. While the chosen star formation and feedback prescriptions combine to expel a significant amount of baryons from the system, they preferentially remove low-angular-momentum material and leave plenty of cold gas available for disk star formation (cf. Scannapieco et al. 2009), thus allowing a good fit to the Tully–Fisher relation (cf. Piontek & Steinmetz 2011).

It is interesting to compare the properties of Eris and ErisLT at redshift 1 (see Table 1).

• 1.  
  The ErisLT control run produces a galaxy resembling an early-type Sa spiral with B/D = 0.42 (cf. Eris' B/D = 0.31), closer to (but still on the low side of) the typical outcome of previously published cosmological simulations of disk formation (e.g., Governato et al. 2010; Scannapieco et al. 2010). Its rotation curve peaks at 308  km s⁻¹ (cf. Eris' 237  km s⁻¹) and declines steeply within a few kpc from the center (Figure 1).
• 2.  
  The baryon fraction in ErisLT is higher than in Eris and close to the universal value. The difference in the baryon content of the two simulations is established at very early times (see the bottom right panel of Figure 6). The baryon content of ErisLT increases from z = 10 to z = 5, as its dark matter halo grows from M_vir = 3 × 10⁹  M_☉ to M_vir = 6 × 10¹⁰ M_☉. At higher redshifts (smaller progenitor mass), the collapse of baryons is heavily suppressed by the ultraviolet radiation background. Eris, by contrast, maintains a relatively low baryon fraction, between 64% and 73% of the cosmic value at all redshift z < 10. Note that these values are higher than the value measured at the present epoch in the bulgeless dwarf galaxy simulation of Governato et al. (2010), which is only 30% of the cosmic fraction.
• 3.  
  While ErisLT was run with a star formation efficiency that was half of that adopted for Eris, its stellar content at z = 1 is 20% higher, showing that the most important parameter in the star formation recipe is not the efficiency but rather the star formation density threshold. Indeed, by allowing the gas to reach higher densities before turning into stars, the ISM develops a more inhomogeneous structure, with important consequences on the large-scale pattern of star formation in the galaxy, and, as a byproduct, on the effect of SNe feedback (e.g., Governato et al. 2010).
• 4.  
  At z = 1, about 14% of the total gas content in ErisLT is in a cold phase, 19% is warm, and 67% is in a hot component. The same fractions are 30%, 18%, and 52% in Eris, i.e., there is two times more cold gas in Eris than in ErisLT. The difference in the density distribution and thermodynamic state of the ISM in the two simulations is clearly seen in Figure 7, which shows the local gas density versus radius of each gas particle at different epochs. The horizontal green line in each panel marks the adopted star formation threshold, and particles are color-coded according to their temperature. Stars are born only in the blue gas particles above the green lines. Because of the increased threshold, fewer regions are forming stars at any given time in Eris. According to Equation (1), however, within these cold dense clumps the number of massive stars per unit gas mass born at threshold scales as $\sqrt{n_{\rm SF}}$, i.e., it is a factor (5/0.1)^1/2 = 7 times larger in Eris versus ErisLT. As these stars explode as SNe, the energy injected at high redshift in Eris's ISM is enough to disrupt and unbind the star-forming regions, generating strong outflows that leave the galaxy and reduce its baryonic content.
• 5.  
  Figure 7 also shows that, at z = 5, star formation can only occur in Eris' inner densest regions within 10 kpc from the center. By contrast, in ErisLT, star formation is spread out over a significant fraction of the host galaxy, extending well beyond 20 kpc. At z = 3, when outflows are stronger due the high star formation rates triggered by mergers, it is the cold, low-angular-momentum gas at the center that is preferentially removed in Eris, as previously found for dwarf galaxies by Governato et al. (2010). Such outflows are weaker in the low-threshold simulation, where star formation is more diffuse and the energy injection from SNe is more evenly deposited in the ISM: cold gas accumulates in the inner regions already at very high redshifts and continues to turn into stars unabated by feedback (see also Ceverino & Klypin 2009; Robertson & Kravtsov 2008). At lower redshift ErisLT has now consumed its low-angular-momentum cold gas in star formation, while Eris has preserved a higher cold gas fraction due to the more effective regulation of star formation via feedback. By z = 1 both galaxies have formed roughly the same amount of stars, but 90% of the gas within 20 kpc from the center is still cold in Eris, compared to only 57% in ErisLT. The distributions of stellar angular momenta are also different in the two cases (see Figure 6). Gas being blown out at high redshift has a systematically lower angular momentum than the gas that gets accreted at later times (Brook et al. 2011), and this explains the flatter rotation curve of Eris relative to ErisLT.

Zoom In Zoom Out Reset image size

It is fair to point out that, while the success of our Eris simulation in matching the observations appears to be linked to the ability of correctly following star formation in an inhomogeneous ISM and regulating it with SN feedback, many avenues remain unexplored and require further investigation. Simulations of even higher resolution, approaching the true gas densities reached in star-forming giant molecular cloud complexes (about a factor of 10 higher than adopted here), are needed to test the convergence and robustness of our results, and are in the making. The density threshold should depend on metallicity and therefore on redshift (Gnedin et al. 2009), which may affect the structure of the ISM and SN feedback differently in progenitors at different epochs. The feedback model that we use is still phenomenological, and the actual mechanism of outflow generation may require the combination of more than one effect to support a large-scale blastwave (e.g., Ceverino & Klypin 2009). Lastly, the simulations presented in this paper neglect cooling by metal lines at temperatures above 10⁴ K. At the mass scale considered here, cold flows are mostly responsible for the assembly of the star-forming disk even at low redshift, as opposed to the cooling flow of the hot halo mode (Brooks et al. 2009). This suggests that the details of the cooling function for gas above 10⁴ K, namely in the hot mode, are not important for the assembly of the disk. Piontek & Steinmetz (2011) find that metal cooling gives origin to a stronger burst of star formation at high redshift in the progenitors of Milky Way sized galaxies and thus to bigger bulges. This result, however, was obtained with the canonical low star formation density threshold (n_SF = 0.1 atoms cm⁻³), and it is unclear in which temperature range the effect of metal cooling is most crucial. These authors also show that the augmented bulge can be suppressed by boosting the effect of (thermal and kinetic) SN feedback. It is conceivable that, by further increasing the density threshold for star formation toward actual giant molecular cloud densities, along with modeling the formation of molecular hydrogen, the ISM may become increasingly clumpy and dense, populating the high surface density tail of the Kennicutt–Schmidt relation that remains currently unresolved in Eris. As a result, one may expect that heating and outflows from even more localized SN feedback may become stronger in such high-density regions, and eventually offset the impact of increased cooling, as suggested by Piontek & Steinmetz (2011). We plan to explore these issues in the next generation of simulations with increased resolution, a higher star formation density threshold, and the inclusion of molecular phase physics. In addition to the effect on disk and bulge, metal cooling at T > 10⁴ K may also have an impact on the density profile and temperature of the hot halo. Piontek & Steinmetz (2011) find a significant reduction of the mass of the hot halo when metal cooling is included, under the assumption of a gas metallicity of Z = 0.5 Z_☉. This is higher than the metallicity measured in spiral galaxies with extended X-ray and Hα emission around their disk, which is in the range Z = 0.01–0.1 Z_☉ (Rasmussen et al. 2009), suggesting that the effect may have been overestimated. The mean metallicity of hot gas in Eris closely matches the latter observations, being on average Z = 0.08 Z_☉ (with Z_☉ = 0.0194; Anders & Grevesse 1989).

The last major merger in our simulations occurs at z ∼ 3, and therefore Eris is expected to show some offset from the observed structural parameters of the average spiral galaxy. The same likely applies to the Milky Way itself, which indeed closely resembles Eris. Whether or not the good match with the properties of typical spiral galaxies is related to the fact that we have selected a particularly quiet merging history, in which outflows shut off early as the rate of star formation drops after the last major merger, will have to be investigated.

At z < 3, the evolution of the B/D ratio is non-monotonic and it is seen to decrease following a major merger at z = 3 and a minor merger at z = 1, and to grow secularly at lower redshifts. The details of this evolution will be the subject of a forthcoming paper. More typical galaxies undergoing significant mergers at later times may develop smaller bulges due to the more prolonged effect of SNe outflows, which could explain why 11 out of 19 nearby massive spirals show no evidence for a classical bulge (Kormendy et al. 2010)). If true, the trend would be at odds with the standard picture in which mergers lead to earlier type objects. This would be an important new ingredient, perhaps complementary to the finding that gas-rich mergers assist the growth of larger disks (Hopkins et al. 2009; Governato et al. 2009).

This research was funded by NASA through grant NNX09AJ34G and by the NSF through grant AST-0908910 (P.M.), by the Swiss National Foundation (SNF), and by an ARCS Foundation Fellowship to J.G. Simulations were carried out on NASA's Pleiades supercomputer, the UCSC Pleiades cluster, and the Swiss National Supercomputing Center's ROSA Cray-XT5. We thank the referee for useful comments that improved this paper and acknowledge useful discussions on the topic of this paper with Oscar Agertz, Alyson Brooks, Valentino González, Fabio Governato, Dušan Kereš, Andrey Kravtsov, Mark Krumholz, Brant Robertson, Sijing Shen, and Romain Teyssier. We are indebted to Frank Bigiel for helping with THINGS data, Charlie Conroy for providing a photometric stellar mass estimate for Eris, and Elena D'Onghia for helping in generating the initial conditions and selecting the Eris halo. L.M. thanks the Aspen Center of Physics for hospitality during the early stages of the work.