Cosmological Galaxy Formation Simulations Using SPH


We present the McMaster Unbiased Galaxy Simulations (MUGS), the first 9 galaxies of an unbiased selection ranging in total mass from 5 M to 2 M simulated using n-body smoothed particle hydrodynamics (SPH) at high resolution. The simulations include a treatment of low temperature metal cooling, UV background radiation, star formation, and physically motivated stellar feedback. Mock images of the simulations show that the simulations lie within the observed range of relations such as that between color and magnitude and that between brightness and circular velocity (Tully-Fisher). The greatest discrepancy between the simulated galaxies and observed galaxies is the high concentration of material at the center of the galaxies as represented by the centrally peaked rotation curves and the high bulge-to-total ratios of the simulations determined both kinematically and photometrically. This central concentration represents the excess of low angular momentum material that long has plagued morphological studies of simulated galaxies and suggests that higher resolutions and a more accurate description of feedback will be required to simulate more realistic galaxies. Even with the excess central mass concentrations, the simulations suggest the important role merger history and halo spin play in the formation of disks.

galaxies: evolution — galaxies: formation — methods: N-Body simulations

1 Introduction

Forming a galaxy like our own Milky Way remains a challenge for the currently accepted Cold Dark Matter (CDM) cosmogony. The Milky Way is comprised of three distinct, stellar components: a flattened, rotating disk; a compact, central and spheroidal bulge; and a diffuse, spherical halo of stars. Any consistent cosmogony needs to explain the origin and evolution of each of these components. CDM posits that the energy budget of the Universe is currently dominated by vacuum energy (), and the majority of the mass is invisible (dark) and only interacts with baryons via gravity. Thus, in the early Universe, thermal baryonic pressure did not support the dark matter, and because it is non-relativistic (cold) the dark matter first collapsed into small structures. Subsequently, the small structures merged hierarchically to form larger structures like the Milky Way.

The CDM paradigm provides a simple explanation for the formation of the stellar halo: stars formed early in small satellite galaxies, which got tidally stripped as their orbits brought them inside the tidal radius of the main galaxy. While early observations indicated that stars in the Milky Way’s halo might have condensed out of a monolithically collapsing gas cloud (Eggen et al., 1962), later observations found instead that formation through mergers like those proposed in the CDM paradigm are more likely (Searle & Zinn, 1978). Today, digital surveys of the sky reveal structures in the Milky Way’s stellar halo such as streams and the remnant cores of dwarf galaxies (Majewski, 1993; Belokurov et al., 2006). These are the exact signatures left by tidally disrupted satellites in simulations (Bullock & Johnston, 2005; Abadi et al., 2006).

Unlike halos, disks are not so neatly explained by CDM. Although conservation of angular momentum naturally creates rotating disks, the hierarchical buildup of structure impairs their formation, since disks form most efficiently in a quiet environment where gas cools and collapses smoothly. In CDM-inspired simulations of substructure mergers, satellites orbiting disks tidally heat stars turning thin disks into thick disks (Toth & Ostriker, 1992; Quinn et al., 1993; Velazquez & White, 1999; Kazantzidis et al., 2008). Simulations of larger galactic mergers transform disks into centrally concentrated, spheroidal systems as the disks experience significant angular momentum loss (Barnes & Hernquist, 1996; Cox et al., 2006). However, the observed distribution of galaxy morphologies can be reproduced with simple, analytic models. In these models uniformly rotating spheres collapse into centrifugally supported disks (Fall & Efstathiou, 1980; Dalcanton et al., 1997; Mo et al., 1998; van den Bosch, 2001). The spheres start with angular momentum and mass profiles predicted using simulations of CDM structure formation. The contradiction between mergers and disk formation may indicate that halo spin, not merger history, plays the dominant role in determining the morphology of galaxies.

The bulge is a spheroid of stars at the center of a galaxy. The spheroidal shape indicates that they formed as the result of mergers. However, other evidence reveals that some bulges may have a secular origin. Kormendy (1993) suggested that observations of rapidly rotating bulges indicates the existence of “pseudo-bulges”, and recent simulations show that the central regions of isolated disks can buckle and cause stars to evolve through “peanut” shaped orbits into spheroidal distributions (Debattista et al., 2004).

Whether disks or spheroids form affects many galactic properties like their kinematics, color, light distribution, star formation history and metallicity in addition to morphology. Disk kinematics are dominated by ordered circular velocity, while in spheroidal components random velocity dominates over circular velocities. Because star formation usually happens in galaxy disks where gas densities are sufficiently high, galaxies with more prominent disks have more recent star formation and display bluer colors, while spheroids tend to be more metal-rich and red. The prominence of morphological components also has an impact on the radial distribution of galactic light profiles. Galaxies with more prominent spherical components exhibit more centrally concentrated light profiles.
Models of galaxy formation require high resolution, hydrodynamic numerical simulations. Analytic modeling can evolve a CDM-motivated Gaussian density field into a spectrum of mass structures and populate those structures with stars such that they match the observered luminosity function of galaxies (Cole et al., 2000; Benson et al., 2003; Somerville et al., 2008). However, because of the non-linear interactions of processes such as gas cooling, merging, tidal stripping, star formation, stellar feedback and active galactic nuclei, it is difficult for these models to predict galaxy morphologies, though Benson & Devereux (2009) and Dutton (2009) represent recent attempts.

As an alternative, numerical simulations allow us to study how halos, disks, and bulges were created and evolve. Simulations that include gas are able to follow more physical processes than simulations that only track the gravitational interaction of dark matter. Gas can be modeled using smoothed particle hydrodynamics (hereafter SPH), which partitions the gas in the Universe into particles and with a Lagrangian approach follows the motions of those particles. SPH is effective because it concentrates computational resources on the high density regions of a simulation, where galaxies form.

Several studies of galaxies in a cosmological context have generated individual objects that are similar to the observed local galaxies (Governato et al., 2004; Robertson et al., 2004; Okamoto et al., 2005; Brook et al., 2006; Governato et al., 2007; Scannapieco et al., 2008, 2009; Ceverino & Klypin, 2009; Sánchez-Blázquez et al., 2009; Governato et al., 2009; Piontek & Steinmetz, 2009; Martínez-Serrano et al., 2009). These require high resolution and large computational resources in order for several important properties of the simulations to converge, such as the galactic structure, motions and star formation history. As a consequence of the high computational cost, the initial conditions have been carefully chosen to maximize the chance that the simulation will produce the desired type of object (usually a disk galaxy).

The previous cosmological simulations have shown that the collapse of gas is more complicated than smooth collapse into centrifugally supported disks. They show that a significant fraction of gas flows into galaxies along cold filaments (Kereš et al., 2005; Brooks et al., 2009). This raises the question: are centrifugal forces all that should support disks? Observations suggest that disks are supported by an equipartition of energy between thermal, magnetic, and cosmic ray pressure support (Cox, 2005).

1.1 Stellar Feedback

One way to introduce pressure support into simulations is by harnessing the energy massive stars release in stellar winds and supernovae explosions. Recent simulations continue to suffer from the overcooling that has long plagued morphological studies of simulated galaxies. Navarro & Benz (1991) described how angular momentum transferred from gas in the disk to the dark matter causes excess central mass concentrations, which leads to massive bulges that do not compare well to observations of disk galaxies whose bulges are fainter and less massive than their disks (Allen et al., 2006). Navarro & Benz (1991) instead proposed that stellar feedback could eliminate a significant amount of low angular momentum gas. They implemented a method to kinematically excite gas particles around recently formed stars, but this showed little improvement in eliminating low angular momentum gas (Navarro & Steinmetz, 2000).

Stellar feedback plays a larger role in the development of satellite galaxies that merge with the parent galaxy. The effects of feedback in satellites can contribute to the final morphology of the main galaxy. If a satellite brings in more gas, it will contribute to make a larger disk; more stars will produce a larger spheroid. In the main galaxies, stellar feedback may also determine the rate at which gas loses angular momentum and migrates through the disk into the dense, star forming center.

Due to insufficient resolution, there is currently no satisfactory treatment for stellar feedback. Multiphase gas particles attempt to capture the phenomenology of the ISM inside individual particles (Springel & Hernquist, 2003a), but it is difficult to determine the appropriate pressure for each particle to exert on the others and sometimes a stiff equation of state needs to be enforced (Springel & Hernquist, 2005). If the gas dynamics are separated into hot and cold gas (Ritchie & Thomas, 2001), it is ambiguous when and how much gas should move from one phase to the other. Disabling the gas cooling reproduces stellar feedback (Gerritsen, 1997; Thacker & Couchman, 2000; Stinson et al., 2006), but resolution is often insufficient to turn off cooling in the proper amount of gas for the right length of time, and SPH does not allow single particles to create their own outflow. Driven winds reproduce galactic outflows (Navarro & Steinmetz, 2000; Springel & Hernquist, 2003b; Oppenheimer & Davé, 2006; Okamoto & Frenk, 2009), but it is yet to be determined whether wind particles should be allowed to interact with surrounding gas and whether they provide sufficient pressure support to forming disks. Many potential avenues need to be followed to see how each feedback recipe affects galaxy formation differently. To date, no stellar feedback model has been effective at reducing the central mass concentration in high mass galaxies, though strong adiabatic feedback combined with altering the star formation threshold to a higher density in very high resolution simulations of low mass galaxies has had the most success (Mashchenko et al., 2008; Governato et al., 2010).

Computational resources can now support surveys that include a range of galaxies simulated at moderate (1 million dark matter particles inside ) resolution (Scannapieco et al., 2009; Okamoto & Frenk, 2009; Piontek & Steinmetz, 2009) with different stellar feedback recipes. Each survey uses SPH, and each of the previous surveys have used gadget. Scannapieco et al. (2009) separated hot and cold gas and used supernova feedback as a conduit from the cold to the hot phase. Okamoto & Frenk (2009) used multiphase particles from Springel & Hernquist (2005) combined with driven winds of different strengths. Scannapieco et al. (2009) found lower disk fractions than late type spirals in their simulations. Piontek & Steinmetz (2009) tested several stellar feedback recipes and did not find striking success with any of them.

Observational samples of galaxies from galaxy redshift surveys, like the 2dFGRS (Colless et al. 2001) and SDSS (York et al. 2000), now contain millions of objects, allowing a much more complete view of not only the properties of typical galaxies, but also a quantification of how galaxies are distributed within the multivariate parameter space of galaxy properties. In contrast, while many researchers are performing sophisticated galaxy formation simulations, each can only produce a handful of galaxies. Evaluating the success of these simulations requires a larger sample of simulated galaxies. When only a small number of simulations exist, it is easy to find a good observational match for any one simulation; however, when a sample of simulated galaxies exist that predict a mean and spread for any galactic property, these predictions can be directly confronted with the large observational samples.

In order to address these problems, we have begun the McMaster Unbiased Galaxy Simulations (MUGS) project. The goal of MUGS is to generate a large sample of sophisticated galaxy formation simulations that sample potential sites of L* galaxy formation in an unbiased manner for direct comparison to the large observational samples now available. In this paper, we describe the methodology used for MUGS and present an overview of the first 9 simulations, particularly focusing on the relative formation of disks and bulges.

MUGS provides an extended look at galaxies simulated using similar physics to Governato et al. (2007) and Governato et al. (2009). Namely, supernovae are modeled with adiabatic “blastwave” feedback described in Stinson et al. (2006). In §2, we describe how we created the initial conditions for MUGS. §3 details the algorithm that evolves the simulations. §4 examines the properties of the galaxies including their brightness, color, mass-to-light ratios, density profiles, bulge-to-total ratio, star formation history, and metallicity.

2 Simulations

The simulations that comprise the MUGS sample each use the volume renormalization techniques from many previous simulations. The technique allows high resolution in a cosmological context at reasonable computational cost. It focuses resolution in one specific region of a cosmological volume while simulating the rest of the volume at lower resolution. The surrounding volume provides large scale density waves and impart tidal torques on the region of interest (Quinn & Binney, 1992).

We selected our galaxies from a cosmological cube 50 Mpc on a side containing dark matter particles that was evolved to . The simulation is based on a 4096 realization of the CMBFAST (Seljak & Zaldarriaga, 1996) power spectrum initially subsampled to 256 to create a uniform resolution, dark matter-only volume. It uses a WMAP3 CDM cosmology with = 73 km s Mpc, =0.24, =0.76, =0.04, and =0.76 (Spergel et al., 2007).

The uniform volume was evolved to =0 at which point the friends-of-friends algorithm was used to find the virialized halos with a linking length = inter-particle separation (Davis et al., 1985). Every group with a mass between M (705 particles in ) and M (2820 particles in ) was examined to ensure that it did not evolve closer than 2.7 Mpc from any halo more massive than M. Massive structures contain hot gas that would significantly alter the evolution of a galaxy of interest. Simulating structures larger than M at the resolution of these galaxies is currently computationally unfeasible. Out of the 36,193 halos found with friends-of-friends, 761 were in the right mass range, and 276 of those were sufficiently isolated. From that sample, 9 halos were selected using a random number generator for more detailed simulation without regard for spin parameter or merger history.

Particles within 5 of each group’s center at z=0 were traced to their positions in the initial conditions to specify the region of interest. In an effort to efficiently use computational resources without loss of physical reality, the regions of interest have a non-spherical shape. Minimizing the number of particles is critical since the n-body tree code must reconstruct the entire tree at each minor timestep, an operation that depends directly on the number of particles. The central region was filled with a regular grid of particles to achieve an effective resolution of 2048 at the center. Surrounding the non-spherical central region is a spherical region with a radius 1.2 times the maximum radius of the central region. This immediately surrounding region is populated with particles for an effective resolution of . Outside this are three spherical regions equally spread in radius with effective resolutions of , , and . The outskirts of the 50 Mpc cube are filled at an effective resolution of . Each of these regions contains progressively more massive particles corresponding to the reduced resolution. Using this scheme, there are between 4 and 10 times more particles in the high resolution region than in the surrounding low resolution regions combined. The regular grids of particles in each region were perturbed using the Zel’dovich approximation with subsampled force resolutions matching the particle resolutions. This dark matter-only configuration was evolved to =0. The particles that ended up inside in the dark matter-only resimulation had a fraction of their mass converted into a separate, neighboring gas particle with a mass corresponding to the cosmic baryon fraction, . We then reapplied the Zel’dovich perturbations to this new regular grid. Figure 1 shows the initial configuration for a sample galaxy. In the high resolution region, dark matter particles have a mass of M and gas particles have an initial mass of M. Stars form with a mass of M. Each particle uses a gravitational softening length of 312.5 pc.

Figure 1: Projection of initial particle configuration for g15784. The red particles represent gas and the light blue represents dark matter. Gas is only included in the high resolution region at the center. The dots become less dense further away from the center because of the decreased resolution and higher particle masses at larger radii.

2.1 Galaxy Sample

One of the main goals of this study is to develop a better understanding of what sorts of galaxy morphologies are created using a wide range of initial conditions. The randomly selected MUGS sample spans a wide range of merger histories and halo angular momenta. The merger history is characterized as the redshift at which the galaxy obtained half of its final mass. The angular momentum of each halo is measured using (Bullock et al., 2001). Figure 2 shows how the and half mass redshifts of the halos in our 9 halo sample compare to the distribution of halos that passed our selection criteria. The values shown in Figure 2 are the values from the low resolution, uniform dark matter only simulation. The values found when the galaxies are run at high resolution including baryonic physics vary modestly from these values. Figure 2 shows that in the initial sample of 9 galaxies does not include any from the high tail. As the sample grows in the future, it will more fully reproduce the multivariate distribution of halo properties. We note that only one halo with half mass redshift greater than 1.5 has a values greater than 0.05, but that some of the halos that accrete half of their mass later have values that surpass 0.1. Ryden (1988) point out that late accretion contributes more angular momentum to halos because of their larger turn around radius. It is expected that galaxies with the quietest merger histories will form galaxies with the most prominent disks.

Figure 2: Three projections of the galaxy sample in mass, spin, and half mass redshift space. The light dots represent all the galaxies in our sample mass range from the uniform volume. The plus signs indicate the galaxies presented here.

3 Code

The simulations were evolved using the parallel SPH code gasoline (Wadsley et al., 2004). gasoline solves the equations of hydrodynamics, and includes radiative cooling. Gravity is calculated for each particle using the Barnes & Hut (1986) tree algorithm with tree elements that span at most = 0.7 of the size of the tree element’s distance from the particle. gasoline is multistepping so that each particle calculates its forces once every gravitational timestep , where = 0.175, is the gravitational softening length (312.5 pc), and is the acceleration. For gas particles, the timestep must also be less than , where = 0.4, is the gas smoothing length and is the sound speed.

The cooling is calculated from the contributions of both primordial gas and metals as . The primordial cooling follows the non-equilibrium evolution of internal energy along with three primordial gas species (HI, HeI, and HeII). H cooling is not included. The scheme uses the collisional ionisation rates reported in Abel et al. (1997), the radiative recombination rates from Black (1981) and Verner & Ferland (1996), and bremsstrahlung and line cooling from Cen (1992). The metal cooling grid is constructed using CLOUDY (version 07.02, last described by Ferland et al. (1998)), assuming ionisation equilibrium. A uniform ultraviolet ionizing background, adopted from Haardt & Madau (in preperation; see Haardt & Madau (1996)), is used in order to calculate the metal cooling rates self-consistently. The cooling lookup table is linearly interpolated in three dimensions (i.e., , z, T) and scaled linearly with metallicity. The energy integration independently uses a semi-implicit stiff integrator for each particle with the compressive heating and density (i.e. terms dependent on other particles) assumed to be constant over the timestep.

The star formation and feedback recipes are the “blastwave model” described in detail in Stinson et al. (2006) with additional improvements as described in § 3.1. They are summarized as follows. Gas particles must be dense () and cool ( = 15,000 K) to form stars. A subset of the particles that pass these criteria are randomly selected to form stars based on the commonly used star formation equation,


where is mass of stars created, is a constant star formation efficiency factor, is the mass of gas creating the star, is how often star formation is calculated (1 Myr in all of the simulations described in this paper) and is the gas dynamical time. The constant parameter, , is tuned to 0.05 so that the simulated Isolated Model Milky Way used in Stinson et al. (2006) matches the Kennicutt (1998) Schmidt Law, and then is left fixed for all subsequent applications of the code. This star formation and feedback treatment was one of the keys to the success of Governato et al. (2007) in producing realistic spiral galaxies in a cosmological simulation and the success of Brooks et al. (2007) in matching the observed mass-metallicity relationship.

At the resolution of these simulations, each star particle represents a large group of stars (6.32 M). Thus, each particle has its stars partitioned into mass bins based on the initial mass function presented in Kroupa et al. (1993). These masses are correlated to stellar lifetimes as described in Raiteri et al. (1996). Stars larger than 8 explode as supernovae during the timestep that overlaps their stellar lifetime after their birth time. The explosion of these stars is treated using the analytic model for blastwaves presented in McKee & Ostriker (1977) as described in detail in Stinson et al. (2006). While the blast radius is calculated using the full energy output of the supernova, less than half of that energy is transferred to the surrounding ISM, ergs. The rest of the supernova energy is assumed to be radiated away. Iron and Oxygen are produced in SNII according to the analytic fits used in Raiteri et al. (1996):


The iron, oxygen, and the supernova energy ejected from SNII are distributed to the same gas within the blast radius. Each SNIa produces 0.63 Iron and 0.13 Oxygen (Thielemann et al., 1986) and ejects it into the nearest gas particle for SNIa.

3.1 Quantized Stellar Feedback

One of the aspects of stellar feedback not given detailed consideration in Stinson et al. (2006) is the clustered nature of star formation. In Stinson et al. (2006), supernova feedback is chronologically distributed according to the stellar initial mass function and lifetimes. The combination of Padua group stellar lifetimes with the Kroupa et al. (1993) IMF results in a constant, small energy release each feedback timestep (1 Myr, concurrent with star formation) for the 35 Myr until the lowest mass stars that explode as supernovae (8 M) explode. Since the blastwave radius and cooling shutoff time are calculated from the energy released during a timestep, the Stinson et al. (2006) method results in small blastwaves. McCray & Kafatos (1987) describe how the accumulated energy of stellar winds and supernova feedback create large superbubbles around star clusters. With a stochastic treatment of the energy release timing, we use the accumulated energy of all the stellar feedback that should result from a star particle to produce a larger, more realistic blastwave. The larger blastwaves should provide more pressure support to make more extended disks. Using a Kroupa et al. (1993) initial mass function, one supernova mass ( 8 M) star is created for every 200 M of stars that form. Combining globular cluster (Harris, 1996), open cluster (Piskunov et al., 2008), and embedded cluster (Lada & Lada, 2003) catalogues, we estimate that a typical star forms in a cluster of 4000 M. Thus, we require a minimum energy release of 20 supernovae ( ergs) for the MUGS simulations.

We stochastically determine when a star particle releases feedback energy. The probability that energy is released at each feedback timestep is


where is the number of supernovae calculated to explode during that star formation timestep and is the “supernova quantum”, the number of supernova required per explosion. If the probability is greater than a random number selected between 0 and 1, supernovae’s worth of energy is released. This causes SN energy to be released sporadically over the 35 Myr until the largest star remaining is M.

4 Results

In this section, we examine the galaxies that form from our set of cosmological initial conditions. First, we compare the properties of mock images of the simulated galaxies with observed galaxies. Second, we examine which parameters have the greatest impact on how much disk or spheroid forms. Third, we compare how much light is produced in our galaxies with observations of the Tully-Fisher relationship and weak lensing mass measurements. Finally, we present the star formation histories and metallicity trends for the simulated galaxies.

4.1 Tabulated Results

In order to identify the main halo and its satellites, we used the Amiga Halo Finder (AHF) (Knollmann & Knebe, 2009). AHF overcomes the difficulties Friends-of-Friends has in separating neighboring groups, by instead identifying density peaks using an adaptive mesh algorithm. The adaptive mesh algorithm naturally leads to identification of substructure, which is critical for analyzing cosmological simulations. Table 1 summarizes key properties for each of the galaxies found with AHF in the sample. Galaxies are identified using the group number from the original friends-of-friends galaxy catalog. The columns are described below.

Galaxy Mass Gasless
( M) (all masses in M) DM
g1536 7.0 0.025 2.9 1.8 0.159 5.1 6.0 1.8 3.9 2.4 1.4 5.3 1640
g5664 5.2 0.025 3.4 1.3 0.164 3.8 4.8 1.3 3.3 1.8 1.1 4.0 0
g7124 4.5 0.039 0 1.5 0.163 2.5 4.8 0.12 3.2 1.1 1.1 3.4 147
g15784 14 0.0345 3.42 1.4 0.150 10 11 3.3 5.5 4.8 2.6 11 2
g21647 7.7 0.048 0.18 0.53 0.168 5.8 7.1 1.1 3.3 2.7 1.6 5.8 6
g22437 8.8 0.013 1.56 1.2 0.170 7.6 7.3 0.56 5.6 3.5 1.6 6.6 1
g22795 8.7 0.011 0.032 1.3 0.144 6.0 6.4 0.26 5.4 2.7 1.5 6.6 76
g24334 7.7 0.041 0.085 1.2 0.162 7.1 11 1.1 5.4 3.3 2.4 8.2 60
g25271 13 0.014 2.9 1.3 0.147 8.9 11 2.2 7.2 4.0 2.4 10 3
Table 1: Simulation data

The columns are defined as follows. Mass is the total mass in M located inside at of the final simulation including baryonic physics. is the spin parameter of that matter defined as . is the redshift of the last major merger, which is defined as the redshift when the AMIGA halo finder no longer distinguished a satellite from the main halo where the satellite was greater than one-third the mass of the main halo at some point in its history. is the redshift when the main halo was one-half its mass at . is the mass of gas inside at . is the baryon fraction of the final halo (). is the mass of stars inside at . and are the mass of stars kinematically classified as part of the disk and bulge, respectively, as described in § 4.4. All the gas and stellar masses are reported in M. The numbers of gas, star, and dark matter particles inside at are , , and , respectively.

One of the problems inherent in running simulations where only a localized region is populated with high resolution particles is that it is possible for low resolution particles or particles from outside the gas region to pollute the region of interest and cause unphysical results. No low resolution particles lie within at of any of the MUGS simulations. Some dark matter particles (“gasless DM”) that originated outside the gas region ended up inside the virial radius. The presence of these particles without corresponding gas indicates that some of the gas in the simulation experienced slightly less pressure than in reality. The initial conditions for g1536 were the first ones that were generated, before the final criteria were firmly in place, and only contained gas particle pairs for dark matter particles that ended up inside ; as a consequence, it contains many more Bad DM particles than the other simulations.

4.2 Simulated Images

The most intuitive way to compare simulations with observations is through mock images of the simulations. Such images can be created by assigning stellar population models like Starburst99 (Leitherer et al., 1999) to each star particle to determine the color and luminosity each star particle should contribute to an image. Additionally, dust can modulate the image with extinction and scattering. Sunrise (Jonsson, 2006) is a Monte Carlo ray tracing program that produces simulated images by assuming dust exists in metal rich gas. Figure 3 shows images 50 kpc on a side that include scattering and absorption and are produced using sunrise. The image brightness and contrast are scaled using as described in Lupton et al. (2004) since disks have an exponential surface brightness profiles meaning the images span a wide range in surface brightness.

Figure 3: Edge-on mock images of the simulated galaxies using the , , and filters found using the Monte Carlo radiative transfer program sunrise. Each image is a 2D projection of a box 50 kpc on a side.

Each galaxy is aligned so that the total angular momentum of the gas inside 1 kpc is pointed upwards in Figure 3. This presents the simulated galaxies edge-on to demonstrate the relative size of the disk and spherical component. In several of the images, a thin disk of young, blue stars is surrounded by a halo of old, red stars. In other images, little disk component is evident and the spherical component dominates. We note that several of the galaxies (g5664, for example) are not perfectly aligned indicating that the disks are warped or that the stars are orbitting around a slightly different axis than the gas. g7124 appears to be elongated vertically rather than horizontally because it is dominated by its spheroid which is strongly elongated perpendicular to its angular momentum of the inner gas, i.e. it is a spheroid rotating about the minor axis.

In contrast to the images presented in Figure 3, the magnitudes and surface brightnesses used throughout the rest of the paper are derived from the face-on projection of the galaxies generated by sunrise, for which the extinction effects are minimized.

4.3 Color-Magnitude Relationship

One simple quantitative evaluation of the simulations is how the color and brightness of galaxies compare to observations. Figure 4 shows a versus absolute color-magnitude diagram (CMD) of galaxies from the Sloan Digital Sky Survey (SDSS) as the shaded two-dimensional histogram (Bailin & Harris, 2008). The colors and magnitudes have been inclination-corrected to their expected face-on values. The two well-known features of the observed CMD are the relatively narrow red sequence, which extends to very bright galaxies, and the more broad blue cloud, which is abruptly truncated at . Relatively few galaxies are observed to lie in the intermediate “green valley”. The simulated galaxies are overplotted as the labelled points.

Figure 4: Simulated galaxy color as a function of absolute magnitude in the Sloan Digital Sky Survey (SDSS) filter bandpasses overlaid on an inclination corrected color magnitude diagram of SDSS galaxies from from Bailin & Harris (2008).

The first conclusion that can be drawn from Figure 4 is that the simulated galaxies lie in observationally-populated regions of the CMD: they are representative of the colors and magnitudes of observed galaxies. Three simulated galaxies, g7124, g22795, and g25271, lie on the red sequence while the remainder are members of the blue cloud. This is a slightly smaller proportion than the of SDSS galaxies that lie on the red sequence (within 0.08 of the mode of the red sequence, i.e. the Bailin & Harris (2008) CMD parameter, ). Given the small size of the simulated sample, this difference should not be overinterpreted. However, a physical reason to expect different red sequence fractions is that the two samples are found in different environments: many observed galaxies are cluster members lying in massive halos, while the simulated galaxies all lie at the centers of galaxy-mass halos.

This is demonstrated in Figure 5, where we have used the group catalogue of Yang et al. (2007) to restrict the observational sample to central galaxies of halos with masses , and that do not have a neighbouring group with halo mass within and a projected radius of  Mpc. This is essentially identical to the selection criteria for our simulated galaxies. The relative fractions of blue versus red galaxies in Figure 5 are indistinguishable from those of our simulated sample: for the observed sample and for the simulated sample. It is also apparent from this Figure that not only do the simulated galaxies lie within the observed ranges of colour and magnitude for the global SDSS population, but they are also a particularly good match to observed galaxies that lie in halos of the correct mass and environment (although perhaps somewhat too luminous, a point we will return to in Section 4.7).

Figure 5: As in Figure 4, but with the SDSS sample restricted to isolated central galaxies with the same halo mass range as the simulated sample.

While the MUGS galaxies are representative of some observed galaxies, they are unusually “green”: the red sequence galaxies lie on the blue edge of the sequence, while the blue cloud galaxies lie in the red half of the cloud. The simulated red sequence galaxies therefore have more star formation activity than usually observed, while the blue cloud galaxies have either less recent or more ancient star formation than observed galaxies. One explanation for the sustained star formation in the red sequence galaxies is that it is a consequence of the environmental effects noted above. However, the simulated galaxies are also too blue in Figure 5, which takes these environmental effects into account. Their blue color could also highlight some numerical failure of our simulations. MUGS does not include AGN feedback, which might be able to drive significantly more star forming gas out of the central regions of galaxies. The blue color may also result from overcooling, where lack of resolution causes excess gas to be driven to the galaxy center where it cools rapidly and forms stars due to the high gas density.

The redness of the simulated blue cloud galaxies could also be a natural result expected for moderate mass galaxies evolving in an isolated environment, though Figure 5 again suggests that this is not the case. Alternatively, it could be the result of excess ancient star formation due to our neglect of the increased feedback from metal-free Population III stars, or overcooling of gas that resulted in too many stars formed in the dense centers of galaxies at early times. Overcooling can in fact simultaneously make the blue galaxies too red by building too large of a bulge (as seen in §4.4), and make the red galaxies too blue by continuing to provide cold gas to the centers of bulge-dominated galaxies at late times.

Regardless of these details, the colors and magnitudes of the simulated galaxies agree well with those of observed galaxies, giving us confidence that analyzing them will provide a template of how galaxies form in the real universe.

4.4 Bulge vs. Disk

Another quantitative comparison between simulated galaxies and observations is how the light and matter are distributed. Figure 6 shows the face-on band surface brightness profiles of the simulated galaxies each fit with the sum of a de Vaucouleurs law (with effective radius, ) and exponential disk (with scale lengths, ) profiles. We calculate a ratio of the light from the bulge to the total light of the galaxy using


(Binney & Merrifield, 1998). The calculated ratios shown in Figure 6 are all higher than 0.5, indicating that bulges dominate the emission from all of our galaxies. In contrast, observed ratios for galaxies with classical bulge/disk profiles are often , (see, e.g., figure 8 of Allen2006 base on the observed Millenium Galaxy Catalog). Again, the high B/T ratios in the simulations indicate that too many stars form in the central region.

Figure 6: Surface brightness profiles of face-on mock -band images of the simulated galaxies, created using sunrise. The profiles are fit with the sum of an law and an outer exponential disk. The bulge-to-total ratios are calculated using equation 5.

Galaxies can also be separated into bulge and disk components based on their kinematics. For this analysis, disks are aligned such that the angular momentum of the gas within 30 kpc of the center of the halo points along the z-axis. Figure 7 shows the ratio of the -component of the specific angular momentum vector to the total specific angular momentum for every star particle, as a function of its three dimensional radius from the center. Thus, Figure 7 represents a comparison between the angular momentum of the gas and the angular momentum of the stars. Each star within 40 kpc of the galactic center is classified as belonging to the disk if greater than 0.75 of its total angular momentum is aligned with the gas, or to the bulge otherwise. The ratio of the masses of these components is given as the bulge to disk ratio. Note, though this is not the same bulge-to-total ratio derived photometrically as in Figure 6, the results are similar as the mass of stars in the bulge is always greater than the mass in the disk. Since the 40 kpc radius is much larger than what is generally considered part of the bulge, hereafter we refer to the kinematically defined component as the spheroid.

Figure 7: The ratio of the -component of the specific angular momentum vector to the total specific angular momentum for every star particle as a function of its three dimensional radius from the center. Stars are sorted into bins 0.04 wide in log(r) and 0.02 high in . The horizontal dashed line at =0.75 indicates the separation between disk and spheroid. The vertical dashed line at r=40 kpc indicates the radial extent inside which particles are classified as part of the spheroid.

What Creates the Disks?

The bulge/disk decompositions show that MUGS simulations remain crude representations of the formation of disk galaxies; the bulge fractions are much larger than what are observed. Still, young stellar disks are apparent in at least g1536 and g15784 in Figure 3, so we briefly examine which halo properties correlate with disk formation in the simulations. Figure 8 shows the photometric bulge-to-total ratio for the galaxies as a function of mass, final spin parameter (), and the redshift of the last major merger. The use of unbiased simulations allows us to draw conclusions from these plots because no selection criteria was based on those variables. One exception to our use of the photometric B/T ratio is g5664, which is the only galaxy with significantly different photometric and kinematic B/T ratios. The photometric decomposition for g5664 is unusual, and seems to display two distinct exponential components. The total bulge+disk fit is therefore poor and the decomposition based on that fit is suspect. Thus, for g5664, we have plotted the mean of the kinematic and photometric B/T ratios in Figure 8.

A correlation is apparent between and B/T. Galaxies with the most recent last major mergers have the largest bulges and typically the reddest colors in Figure 4. One galaxy that had its last major merger long ago, but still has a large bulge is g25271. We note that g25271 has one of the lowest spin parameter values. If we look in detail at the history of g25271, we find that it was involved in enough retrograde minor mergers that it never had a chance to grow a significant disk.

We note that the galaxies with the most significant bulge separate into two groups, one with higher spin and one with lower spin. The higher spin galaxies all experienced recent last major mergers (also see D’Onghia & Navarro, 2007). This indicates that it is only possible for galactic halos to gain angular momentum sufficient to increase its spin parameter above 0.04 if it has a late major merger. We further note that that these galaxies with a high spin parameter display the largest bulge fractions. Observed giant low surface brightness galaxies also exhibit large bulges (Lelli et al., 2010). This might confirm the predictions of analytic models of galaxy formation (e.g. Dalcanton et al., 1997) that halos with the highest spin parameters should host low surface brightness galaxies. Governato et al. (2009) showed how disks can form after late major mergers, but refrained from making any direct comparisons to LSBs.

The lower spin group experienced a wide range of merger histories. Based on their spin parameters, it is apparent that they must have undergone a significant number of retrograde mergers.

Figure 8: The photometric bulge-to-total ratios shown in Figure 6 as a function of galaxy mass, spin, and last major merger redshift. Galaxies at the bottom of these plots exhibit the strongest disks although they still have .

What Creates the Spheroids?

Every galaxy in our sample displays a dominant spheroidal stellar component in both photometric and kinematic decompositions. Since spheroids easily form as the result of merging, we examine the formation of the spheroid in the simulation that has the quietest merger history and largest disk (g15784) here to discover whether there are any significant secondary effects. We leave study of spheroid formation in the other galaxies for future work. Figure 9 shows the orbits of stars classified as spheroid stars in g15784 using the kinematic decomposition. The stars are colored by age so that recently formed stars are yellow and stars that formed long ago are blue. The edge-on view shows that many of the recently formed “spheroid” stars are actually orbitting in the disk plane. The face-on view shows that many new stars are being formed in the central region of the galaxy with radial velocities that lead to their classification as “spheroid” stars. While Figure 9 is colored to emphasize the young stars orbitting in a disk, the overall shape of the stars classified as part of the spheroid is spherical.

Figure 9: Velocities of stars kinematically classified as part of the spheroid in g15784 at z=0 are shown as the thin lines. The stars are colored by their age with more recently formed stars yellow and less recently formed stars blue. The image is 80 kpc across.

Figure 10 shows the formation history of the stars in the spheroid at =0. The spheroid stars are separated into three categories based on where they formed, 1) in satellites; 2) within 1 kpc of the center of the main galaxy; 3) or in the disk of the main galaxy. We emphasize that stars in the disk category are no longer part of the disk, but are spheroid stars that formed in the disk. The decomposition was based on the distance stars formed from the center of the main halo. The stars were sorted into bins one million years in length based on their formation time. The center of mass of the stars that formed during that time was calculated by combining the stars formation locations in each time bin with the center of the main halo found using the AMIGA halo finder (Knollmann & Knebe, 2009) at the nearest output (generally full outputs were only written every 100 Myr). The large number of stars that formed within the central kpc of the main halo made finding the star formation center straightforward for the MUGS galaxies. We calculated the three dimensional distance between the center and the stars that formed during each time interval, and placed them into radial bins 1 kpc wide. Because the stellar density drops strongly at the edge of the disk, the existence of an empty radial bin serves as a delineation between the main galaxy and satellites. We therefore classify all stars within this radius (but outside of 1 kpc) as “disk-formed” stars, and all stars outside of this radius as “satellite-formed” stars. The total bulge stellar mass is 1.1 M, the mass formed inside 1 kpc is 2.9 M, the mass formed in the disk is 3.3 M, and the mass formed in satellites is 4.1 M.

Figure 10: Formation history of stars kinematically classified as part of the spheroid in g15784 at . The solid line shows the total star formation history for all the stars that comprise the spheroid at . The blue dotted line shows the portion of that total that formed in the central 1 kpc of the main halo. The green dashed line shows the portion that formed in the disk and was heated to the spheroid. The red dot-dashed line shows the portion of the inner 40 kpc spheroid that formed in satellites.

The final merger that disturbs stars from the disk into the spheroid happens about 13 Gyr into the simulation. After that time, nearly all the stars that constitute the spheroid form in the central 1 kpc. These are the stars that appear yellow in Figure 9. Stars forming in the unresolved central kpc is not limited to the final Gyr, they are a constant source of spheroid stars throughout the evolution of g15784. For a significant fraction of the history of g15784, the mass of stars formed in the central kpc of the main halo is similar to the mass of stars formed in satellites. While the formation of stars in the central 1 kpc is reminiscent of observations of rapidly rotating pseudo bulges (Kormendy, 1993), we cannot draw any firm conclusions since the inner 1 kpc is unresolved in these simulations.

Stars that form in the disk also contribute a large fraction of stars to the final spheroid, though a merger 6 Gyr into the simulation marks the end of the dominant contribution of disk stars to the spheroid. We leave a detailed examination of how stars migrate from the disk to the spheroid for future work, but speculate that most of the migration is due to the tidal disruption caused by merging satellites given the reduction in disk contribution to the spheroid once the merging activity is complete. We note that Debattista et al. (2004) showed that secular evolution accounted for some migration of stars from the disk to the spheroid in high resolution simulations of isolated disks.

Satellites are expected to contribute a large mass of stars to the spheroid since stars are tidally stripped into spheroidal orbits. In the case of g15784’s quiet merger history, the mass of the spheroid formed in satellites is similar to that formed in the inner 1 kpc and the disk. Finally, we note that the star formation in satellites is also largely confined to the unresolved central kpc similar to the main galaxy.

The large number of stars that form in the central kpc and disk echo the findings of Zolotov et al. (2009) and indicate that the predominant spheroid component is due in large part to stars that form in the main galaxy. However, satellites contribute a significant mass of stars to the spheroid and disrupt a large fraction of stars out of the disk, so these simulations do not exclude traditional spheroid creation through satellite mergers and tidal stripping.

4.5 Density Profiles

From the simulations, one can find the three dimensional distribution of all the matter in the galaxy as opposed to just that which is visible. Figure 11 shows the resulting radial density profile. Particles are sorted by radius and placed into 1000 bins, each containing the same number of particles. The dark matter only profiles for every galaxy are shown as the gray solid lines. The dark matter and total density profiles are similar for each simulated galaxy. The dark matter profiles follow a power law or slightly shallower from the virial radius into 0.1 before growing shallower, similar to dark matter only runs (Navarro et al., 1997; Moore et al., 1998; Reed et al., 2005). The total profile continues a steady rise along the power law into 1 kpc (about 0.005 ) due to the stellar density before flattening out to a slope less than inside 1 kpc. Stars dominate the matter profile inside of 2 kpc.

Figure 11: Particles are sorted by radius and placed into 1000 bins each containing the same number of particles. The radius of the bin is plotted as the mass weighted mean of the particle radii for each bin. The total density profile (various colours) for each halo is compared with the profile of the dark matter from the same simulations (thin grey). The arrow is placed at a radius of , the radius inside which Power et al. (2003) determines density profiles can no longer be trusted.

4.6 Rotation Curves

Rotation curves provide a comparison between the density profile found in these simulations and observations. Figure 12 shows the circular velocity, as a function of radius, where . Mock observations are not used because there is too much scatter in the star particle velocities at a given projected distance from the center of the galaxy.

The excessive bulge illustrated in §4.4 is again apparent in the rotation curves. Rather than exhibiting flat rotation curves as are observed, the simulated rotation curves all have a peak at the center due to a large central concentration of mass. Governato et al. (2007) showed that increased resolution without any changes to star formation can also limit the central concentration of matter as the lessened impact of two body interactions with halo dark matter particles minimizes angular momentum losses. Governato et al. (2010) show how higher resolution and consequent modifications to the star formation threshold density can limit the mass concentration in simulations of slightly less massive galaxies.

Figure 12: Rotation speed as a function of radius. The circular velocity () defined as in radial bins for each of our galaxies.

4.7 Tully-Fisher

The Tully-Fisher (TF) relationship shown in Figure 13 compares the luminosity of a galaxy with its rotation velocity. Velocity indicates the dynamical mass for the inner parts of a galaxy, so the TF relationship indicates how bright a galaxy of a given mass should be. Most galaxies in the range of masses that we simulated have flat rotation curves, so the radius at which the rotation velocity is measured is irrelevant. In the simulations, however, the rotation curve rises sharply initially and the declines gradually with radius. We follow the example of Governato et al. (2007) and measure the rotation velocity at 3.5. Governato et al. (2007) find in a resolution study that this is the radius at which velocity converges. It is unclear from Figure 12 that the rotation curves are flat at this radius, but they are flatter than they are at smaller radii. The I band magnitude for the galaxies generally matches the observations of the corresponding velocities, though the velocities remain slightly high. Governato et al. (2010) showed that higher resolution and higher star formation density threshold produce galaxies with more slowly rising rotation curves, so higher resolution simulations of these galaxies is an obvious next step.

Figure 13: I band absolute magnitude as a function of circular velocity at 3.5 disk scale lengths (). The large dots represent the simulated galaxies while the small dots are from the Giovanelli et al. (1997) sample of cluster galaxies and the triangles are from Courteau et al. (2007).

4.8 Mass-to-Light Ratio

The relationship between the luminosity of a galaxy and the mass of its halo is a key observational question that provides constraints on semi-analytic models and strongly affects the interpretation of observational samples with respect to theoretical predictions about halos. In Figure 14, we have plotted the luminosity of our simulated galaxies versus their total (dark plus baryonic) mass. The luminosities are K-corrected to the band (the band redshifted to ) using kcorrect v4.1.4 (Blanton et al., 2003) using the sunrise-generated SDSS magnitudes in order to facilitate making comparison to the observations.

By construction, our galaxies inhabit a relatively narrow band in halo mass, from to . Their luminosities cover a similar relative range in luminosity, from to , for a typical total mass-to-light ratio of within the virial radius.

We have overplotted a variety of observational estimates of the relationship between central galaxy luminosity and the total halo mass onto Figure 14. These observational techniques are: weak gravitational lensing (Mandelbaum et al., 2006, blue dot-dashed line for exponential-dominated late-type galaxies, and red dashed line for de Vaucouleurs-dominated early-type galaxies), a halo model simultaneously fit to galaxy clustering and the luminosity function Cacciato et al. (2009, solid green line), and the kinematics of satellite galaxies More et al. (2009, dotted black line for the mean and shaded gray region for the measured spread in halo mass at a given luminosity). Note that while all of these studies are based on SDSS data, the techniques use completely independent properties of the galaxies to determine the halo mass. We have adopted , as used in the simulations, to convert the -independent units used in these studies into physical units for direct comparison.

When comparing the observations to the simulations, it is important to realize that in the simulations, the systems are selected based on the halo mass, while the luminosity is an output of the simulations; we have therefore plotted the halo mass as the independent variable and the luminosity as the dependent variable in Figure 14. However, in the observational studies, the galaxies are selected based on their luminosity and the halo mass is what is measured; therefore, the luminosity is the independent variable and the halo mass is the dependent variable. This should be kept in mind when considering Figure 14, and is why the weak lensing result around late-type galaxies appears multi-valued.

The halo mass-luminosity relationship for the simulated galaxies has the same shape and similar normalization as the observed relationship, particularly when comparing to the Mandelbaum et al. (2006) and Cacciato et al. (2009) results. The simulated galaxies are, however, all more luminous than the observed galaxies at a given halo mass, or equivalently the observed galaxies lie in more massive halos at a given luminosity. This difference is no larger than the difference seen between observational techniques, so it may not be significant. One aspect of the simulated and observed samples that may be important is that the simulated galaxies are constrained not to lie near large group or cluster halos; given the mass dependence of halo bias, those halos are themselves likely to be more massive. However, it seems unlikely that this is the explanation given that the most isolated of the observational samples is that of More et al. (2009), with whom our results are the most discrepant. A more likely explanation is that this is another symptom of the overcooling problem, which turns too much of the gas in the system into stars and therefore results in an overluminous galaxy.

Figure 14: Galaxy luminosity in the band plotted as a function of the total halo mass including dark matter for the 9 simulated MUGS galaxies. Observational measurements of the mass-luminosity relationship are shown based on weak gravitational lensing (Mandelbaum et al., 2006), halo model fits to the galaxy clustering and luminosity function (Cacciato et al., 2009), and the kinematics of satellite galaxies (More et al., 2009, the measured dispersion is shown as the shaded region).

4.9 Star Formation Histories

One of the major ways to track the evolution of a galaxy is by examining when its stars formed. Figure 15 shows the star formation histories (hereafter, SFH) for each of our simulated galaxies. Peaks in the SFH correspond to merger events which drive starbursts. Mergers tidally disrupt disk gas and excite instabilities that cause gas overdensities (Mihos & Hernquist, 1996; Springel, 2000; Cox et al., 2006). Following the merger peaks, the shape of the star formation history is an exponential increase followed by an exponential decline as the reservoir of gas is exhausted. The galaxies in Figure 15 that have the most and latest mergers correspond to those with the least well defined disks.

Figure 15: The star formation rate plotted as a function of time for all the stars inside () and just stars that end up in the kinematically-defined disk (). Stars are sorted by their formation time into 50 Myr bins and the total mass of stars in that bin is divided by 50 Myr.

4.10 Metallicity Evolution

Trends in metallicity help track the star formation history of galaxies. Specifically, Figure 16 shows the distribution of stars in [O/Fe] as a function of [Fe/H]. Evolution of metallicity proceeds from low [Fe/H] as stars form and produce iron that is ejected during supernova explosions. Type II supernovae (SNII) occur on a shorter timescale than Type Ia, and they eject alpha elements such as carbon, oxygen, and silicon as well as iron. We use oxygen to track the abundance of alpha elements in our simulations. After stellar populations age for 40 Myr, Type Ia supernovae (SNIa) start to explode, ejecting few alpha elements, but large quantities of iron. The range of [Fe/H] over which the [O/Fe] ratio remains high and constant indicates how many SNII explode before SNIa start contributing significant quantities of iron. Thus, regions of active star formation like the Milky disk will generate long, high [O/Fe] plateaus before [O/Fe] decreases due to the influx of iron from SNIa (Bensby et al., 2005; Reddy et al., 2006). Places where star formation is not as efficient like dwarf galaxies exhibit lower [O/Fe] at low [Fe/H] distinct from abundances of stars measured in the Milky Way’s disk and halo (Venn et al., 2004). Models of galaxy formation that couple N-body simulations with analytic prescriptions for the stellar content of satellite galaxies are able to reproduce this dichotomy (Font et al., 2006).

Figure 16: Metallicity distribution, [O/Fe] as a function of [Fe/H], of all the stars inside of the galaxy. Stars are sorted into bins, 0.03 wide in [Fe/H] and 0.01 high in [O/Fe]. Evolution progresses from low [Fe/H] and high [O/Fe] to solar (0.0) [Fe/H] and [O/Fe]. The multiple evolutionary tracks through this diagram correspond to different stellar structures that each have their own star formation history.

Figure 16 shows that the majority of the stars in the simulations form with solar metallicity and abundances. At low [Fe/H], stars form with both high and low [O/Fe]. The peak [O/Fe] values are 0.4, which are lower than the observed abundances of the disk, which have been observed up to 0.6 Bensby et al. (2005); Reddy et al. (2006). This indicates a shortcoming of the simple power law fit used for oxygen enrichment. The fit only allows type II supernovae from stars up to 40 M to produce oxygen, so the chemical model does not capture chemical enrichment from more massive stars that produce more Oxygen compared to Iron (Woosley & Weaver, 1995).

Much like observational galactic archaeology, the distinct metallicity evolution tracks apparent in Figure 16 provide a useful alterative for discovering substructure inside galaxies as substructures exhibit stellar populations with different metallicity signatures. We leave further discussion of such methods for future work.

5 Conclusions

We presented 9 galaxies from the McMaster Unbiased Galaxy Simulations simulated using N-body gravity and SPH. The galaxies are selected from a mass range around the mass of the Milky Way and from isolated environments, but their selection was otherwise unbiased for factors such as accretion history and angular momentum.

The galaxies were examined using the radiative transfer program sunrise to enable comparisons between the simulated galaxies and real galaxies in the observed plane. The simulated galaxies have colors and magnitudes that compare well with a sample of inclination corrected isolated galaxies from SDSS, and in particular separate into the well-known red sequence and blue cloud. However, both simulated populations tend too much towards the “green valley”, indicating that they contain more old stars than blue cloud galaxies and more young stars than galaxies along the red sequence.

The surface brightness profiles of the simulated galaxies can all be fit with exponential disks combined with a central de Vaucouleurs law similar to real galaxies. However, the proportion of bulge to total light (B/T) is higher than what is typically observed. The B/T ratios are also high when the stars are decomposed into the spheroid and disk based on their kinematics. There are no galaxies with a B/T fraction less than 0.5 whereas observed samples find many galaxies with B/T 0.5. This result is similar to that found in many previous simulations. We note that many of the recently formed stars that are classified as part of the spheroid form with orbits in the disk plane in the central regions of the disk. We also note that most of the stars that comprise the spheroid formed in situ, but we leave the question of how the stars may migrate from the disk to a spherical distribution for future work.

As to the question of why galaxies form with more or less spheroid, there seems to be a modest trend with accretion history. We find that the largest disks (g1536 and g15784) form with a quiet merger history in which they had their last major merger prior to =3, while the largest bulges all resulted from recent last major mergers. There also appears to be a dependence on halo spin as g25271 has a relatively quiet merger history, but low and results in a significant bulge and red color. Since all the galaxies used in MUGS fall in a limited mass range, these conclusions do not include the impact of mass.

We compare the brightness of the final galaxies with their mass at two different radii in the final output. First, we compare our galaxies with the observed Tully-Fisher relationship that probes the amount of light a galaxy produces with the mass contained in its inner regions. Since the final mass concentration of all the galaxies is too high, we use a rotation velocity from 3.5 disk scale lengths away from the galaxy center. Previous studies have shown that this is the radius at which rotation velocities converge. The galaxies are still slightly fainter than the observed sample based on these inner velocities. Second, we compare the brightness of our galaxies with observations of the whole halo mass derived using a number of different methods. In each case, the simulated galaxies are brighter than comparable observed galaxies at the same mass. This indicates that too many stars form in the simulation. The high central velocities used in the Tully-Fisher relationship indicate that mass gets too concentrated at the centers of the halos and while the galaxies are fainter than the observed TF relationship, the lack of resolution makes it difficult to determine whether the amount of stars formed is too many or too few.

We are thus left with the challenge of creating more realistic simulations in order to obtain more accurate insights into how the important physical processes involved in galaxy formation result in the observed population of galaxies. Fortunately, there has been much recent work that guides the way forward. It has been shown repeatedly that higher resolution makes better disks (Governato et al., 2004, 2007). More recently, it has been shown that high resolution combined with clustered star formation can remove central density concentrations from dwarf galaxies (Mashchenko et al., 2008; Governato et al., 2010).

While these simulations open many possibilities for expanding our understanding of how galaxies form, they also show that there is much work left to be done before we can claim to have simulated a sample of galaxies that compares well with real ones.


This paper makes use of simulations performed as part of the SHARCNET Dedicated Resource project: “MUGS: The McMaster Unbiased Galaxy Simulations Project” (DR316, DR401, and DR437). We would like to thank Allison Sills, Bill Harris, and Victor Debatista for useful conversations. We would also like to thank Rok Rokar and Peter Yoachim for helpful IDL code that contributed to this project. As should be apparent from the bibliography, much of this work was inspired by Fabio Governato. GS is a Fellow of the Jeremiah Horrocks Institute and did the bulk of this work as a CITA National Fellow. JB and HC were supported by NSERC grant PHY-0205413. JW was supported NSERC grant.


  1. Abadi M. G., Navarro J. F., Steinmetz M., 2006, \mnras, 365, 747
  2. Abel T., Anninos P., Zhang Y., Norman M. L., 1997, New Astronomy, 2, 181
  3. Allen P. D., Driver S. P., Graham A. W., Cameron E., Liske J., de Propris R., 2006, \mnras, 371, 2
  4. Bailin J., Harris W. E., 2008, \apj, 681, 225
  5. Barnes J., Hut P., 1986, \nat, 324, 446
  6. Barnes J. E., Hernquist L., 1996, \apj, 471, 115
  7. Belokurov V., Zucker D. B., Evans N. W., Gilmore G., Vidrih S., Bramich D. M., Newberg H. J., Wyse R. F. G., Irwin M. J., Fellhauer M., Hewett P. C., Walton N. A., Wilkinson M. I., Cole N., Yanny B., Rockosi C. M., 2006, \apjl, 642, L137
  8. Bensby T., Feltzing S., Lundström I., Ilyin I., 2005, \aap, 433, 185
  9. Benson A. J., Bower R. G., Frenk C. S., Lacey C. G., Baugh C. M., Cole S., 2003, \apj, 599, 38
  10. Benson A. J., Devereux N., 2009, \mnras, pp 1955–+
  11. Binney J., Merrifield M., 1998, Galactic astronomy
  12. Black J. H., 1981, \mnras, 197, 553
  13. Blanton M. R., Brinkmann J., Csabai I., Doi M., Eisenstein D., Fukugita M., Gunn J. E., Hogg D. W., Schlegel D. J., 2003, \aj, 125, 2348
  14. Brook C. B., Kawata D., Martel H., Gibson B. K., Bailin J., 2006, \apj, 639, 126
  15. Brooks A. M., Governato F., Booth C. M., Willman B., Gardner J. P., Wadsley J., Stinson G., Quinn T., 2007, \apjl, 655, L17
  16. Brooks A. M., Governato F., Quinn T., Brook C. B., Wadsley J., 2009, \apj, 694, 396
  17. Bullock J. S., Dekel A., Kolatt T. S., Kravtsov A. V., Klypin A. A., Porciani C., Primack J. R., 2001, \apj, 555, 240
  18. Bullock J. S., Johnston K. V., 2005, \apj, 635, 931
  19. Cacciato M., van den Bosch F. C., More S., Li R., Mo H. J., Yang X., 2009, \mnras, 394, 929
  20. Cen R., 1992, \apjs, 78, 341
  21. Ceverino D., Klypin A., 2009, \apj, 695, 292
  22. Cole S., Lacey C. G., Baugh C. M., Frenk C. S., 2000, \mnras, 319, 168
  23. Courteau S., Dutton A. A., van den Bosch F. C., MacArthur L. A., Dekel A., McIntosh D. H., Dale D. A., 2007, \apj, 671, 203
  24. Cox D. P., 2005, \araa, 43, 337
  25. Cox T. J., Jonsson P., Primack J. R., Somerville R. S., 2006, \mnras, 373, 1013
  26. Dalcanton J. J., Spergel D. N., Summers F. J., 1997, \apj, 482, 659
  27. Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, \apj, 292, 371
  28. Debattista V. P., Carollo C. M., Mayer L., Moore B., 2004, \apjl, 604, L93
  29. D’Onghia E., Navarro J. F., 2007, \mnras, 380, L58
  30. Dutton A. A., 2009, \mnras, 396, 121
  31. Eggen O. J., Lynden-Bell D., Sandage A. R., 1962, \apj, 136, 748
  32. Fall S. M., Efstathiou G., 1980, \mnras, 193, 189
  33. Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B., Verner E. M., 1998, \pasp, 110, 761
  34. Font A. S., Johnston K. V., Bullock J. S., Robertson B. E., 2006, \apj, 638, 585
  35. Gerritsen J. P. E., 1997, Ph.D. Thesis
  36. Giovanelli R., Haynes M. P., Herter T., Vogt N. P., da Costa L. N., Freudling W., Salzer J. J., Wegner G., 1997, \aj, 113, 53
  37. Governato F., Brook C., Mayer L., Brooks A., Rhee G., Wadsley J., Jonsson P., Willman B., Stinson G., Quinn T., Madau P., 2010, \nat, 463, 203
  38. Governato F., Brook C. B., Brooks A. M., Mayer L., Willman B., Jonsson P., Stilp A. M., Pope L., Christensen C., Wadsley J., Quinn T., 2009, \mnras, 398, 312
  39. Governato F., Mayer L., Wadsley J., Gardner J. P., Willman B., Hayashi E., Quinn T., Stadel J., Lake G., 2004, \apj, 607, 688
  40. Governato F., Willman B., Mayer L., Brooks A., Stinson G., Valenzuela O., Wadsley J., Quinn T., 2007, \mnras, 374, 1479
  41. Haardt F., Madau P., 1996, \apj, 461, 20
  42. Harris W. E., 1996, \aj, 112, 1487
  43. Jonsson P., 2006, \mnras, 372, 2
  44. Kazantzidis S., Bullock J. S., Zentner A. R., Kravtsov A. V., Moustakas L. A., 2008, \apj, 688, 254
  45. Kennicutt R. C., 1998, \apj, 498, 541
  46. Kereš D., Katz N., Weinberg D. H., Davé R., 2005, \mnras, 363, 2
  47. Knollmann S. R., Knebe A., 2009, \apjs, 182, 608
  48. Kormendy J., 1993, in H. Dejonghe & H. J. Habing ed., Galactic Bulges Vol. 153 of IAU Symposium, Kinematics of extragalactic bulges: evidence that some bulges are really disks. pp 209–+
  49. Kroupa P., Tout C. A., Gilmore G., 1993, \mnras, 262, 545
  50. Lada C. J., Lada E. A., 2003, \araa, 41, 57
  51. Leitherer C., Schaerer D., Goldader J. D., González Delgado R. M., Robert C., Kune D. F., de Mello D. F., Devost D., Heckman T. M., 1999, \apjs, 123, 3
  52. Lelli F., Fraternali F., Sancisi R., 2010, ArXiv e-prints
  53. Lupton R., Blanton M. R., Fekete G., Hogg D. W., O’Mullane W., Szalay A., Wherry N., 2004, \pasp, 116, 133
  54. Majewski S. R., 1993, \araa, 31, 575
  55. Mandelbaum R., Seljak U., Kauffmann G., Hirata C. M., Brinkmann J., 2006, \mnras, 368, 715
  56. Martínez-Serrano F. J., Serna A., Doménech-Moral M., Domínguez-Tenreiro R., 2009, \apjl, 705, L133
  57. Mashchenko S., Wadsley J., Couchman H. M. P., 2008, Science, 319, 174
  58. McCray R., Kafatos M., 1987, \apj, 317, 190
  59. McKee C. F., Ostriker J. P., 1977, \apj, 218, 148
  60. Mihos J. C., Hernquist L., 1996, \apj, 464, 641
  61. Mo H. J., Mao S., White S. D. M., 1998, \mnras, 295, 319
  62. Moore B., Governato F., Quinn T., Stadel J., Lake G., 1998, \apjl, 499, L5+
  63. More S., van den Bosch F. C., Cacciato M., 2009, \mnras, 392, 917
  64. Navarro J. F., Benz W., 1991, \apj, 380, 320
  65. Navarro J. F., Frenk C. S., White S. D. M., 1997, \apj, 490, 493
  66. Navarro J. F., Steinmetz M., 2000, \apj, 538, 477
  67. Okamoto T., Eke V. R., Frenk C. S., Jenkins A., 2005, \mnras, 363, 1299
  68. Okamoto T., Frenk C. S., 2009, \mnras, 399, L174
  69. Oppenheimer B. D., Davé R., 2006, \mnras, 373, 1265
  70. Piontek F., Steinmetz M., 2009, ArXiv e-prints
  71. Piskunov A. E., Kharchenko N. V., Schilbach E., Röser S., Scholz R., Zinnecker H., 2008, \aap, 487, 557
  72. Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, \mnras, 338, 14
  73. Quinn P. J., Hernquist L., Fullagar D. P., 1993, \apj, 403, 74
  74. Quinn T., Binney J., 1992, \mnras, 255, 729
  75. Raiteri C. M., Villata M., Navarro J. F., 1996, \aap, 315, 105
  76. Reddy B. E., Lambert D. L., Allende Prieto C., 2006, \mnras, 367, 1329
  77. Reed D., Governato F., Verde L., Gardner J., Quinn T., Stadel J., Merritt D., Lake G., 2005, \mnras, 357, 82
  78. Ritchie B. W., Thomas P. A., 2001, \mnras, 323, 743
  79. Robertson B., Yoshida N., Springel V., Hernquist L., 2004, \apj, 606, 32
  80. Ryden B. S., 1988, \apj, 329, 589
  81. Sánchez-Blázquez P., Courty S., Gibson B. K., Brook C. B., 2009, \mnras, 398, 591
  82. Scannapieco C., Tissera P. B., White S. D. M., Springel V., 2008, \mnras, 389, 1137
  83. Scannapieco C., White S. D. M., Springel V., Tissera P. B., 2009, \mnras, 396, 696
  84. Searle L., Zinn R., 1978, \apj, 225, 357
  85. Seljak U., Zaldarriaga M., 1996, \apj, 469, 437
  86. Somerville R. S., Hopkins P. F., Cox T. J., Robertson B. E., Hernquist L., 2008, \mnras, 391, 481
  87. Spergel D. N., Bean R., Doré O., Nolta M. R., Bennett C. L., Dunkley J., Hinshaw G., Jarosik N., Komatsu E., Page L., Peiris H. V., Verde L., Halpern M., Hill R. S., Kogut A., 2007, \apjs, 170, 377
  88. Springel V., 2000, \mnras, 312, 859
  89. Springel V., Hernquist L., 2003a, \mnras, 339, 289
  90. Springel V., Hernquist L., 2003b, \mnras, 339, 312
  91. Springel V., Hernquist L., 2005, \apjl, 622, L9
  92. Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006, \mnras, 373, 1074
  93. Thacker R. J., Couchman H. M. P., 2000, \apj, 545, 728
  94. Thielemann F.-K., Nomoto K., Yokoi K., 1986, \aap, 158, 17
  95. Toth G., Ostriker J. P., 1992, \apj, 389, 5
  96. van den Bosch F. C., 2001, \mnras, 327, 1334
  97. Velazquez H., White S. D. M., 1999, \mnras, 304, 254
  98. Venn K. A., Irwin M., Shetrone M. D., Tout C. A., Hill V., Tolstoy E., 2004, \aj, 128, 1177
  99. Verner D. A., Ferland G. J., 1996, \apjs, 103, 467
  100. Wadsley J. W., Stadel J., Quinn T., 2004, New Astronomy, 9, 137
  101. Woosley S. E., Weaver T. A., 1995, \apjs, 101, 181
  102. Yang X., Mo H. J., van den Bosch F. C., Pasquali A., Li C., Barden M., 2007, \apj, 671, 153
  103. Zolotov A., Willman B., Brooks A. M., Governato F., Brook C. B., Hogg D. W., Quinn T., Stinson G., 2009, \apj, 702, 1058
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description