Rotation rates, sizes, and star formation efficiencies of a representative population of simulated disc galaxies
We examine the rotation rates, sizes, and star formation (SF) efficiencies of a representative population of simulated disc galaxies extracted from the Galaxies-Intergalactic Medium Interaction Calculation (gimic) suite of cosmological hydrodynamic simulations. These simulations include efficient, but energetically feasible supernova feedback, but have not been tuned in any way to produce ‘realistic’ disc galaxies. Yet, they generate a large number of discs, without requiring extremely high resolution. Over the wide galaxy stellar mass range, , the simulations reproduce the observed Tully-Fisher relation, the rotation curves of disc galaxies in bins of stellar mass, the mass-size relation of disc galaxies, the optical rotation to virial circular velocity ratio (’’), and the SF efficiencies of disc galaxies as inferred from stacked weak lensing and stacked satellite kinematics observations. They also reproduce the specific star formation rates of galaxies but predict too low levels of star formation for low-mass galaxies, which is plausibly due to the finite resolution of the simulations. At higher stellar masses, , the simulated galaxies are too concentrated and have too high SF efficiencies. We conjecture that this shortcoming reflects the neglect of feedback from accreting supermassive black holes in these simulations. We conclude that it is possible to generate a representative population of disc galaxies that reproduces many of the observed trends of local disc galaxies using standard numerical hydrodynamic techniques and a plausible implementation of the “subgrid” astrophysical processes thought to be relevant to galaxy formation.
keywords:galaxies: evolution — galaxies: formation — galaxies: general — galaxies: haloes — galaxies: stellar content — galaxies: structure
Understanding the formation of normal disc galaxies remains one of the key challenges in astrophysics today. A general physical picture has been in place for some time (e.g., Fall & Efstathiou 1980; White & Frenk 1991; Mo et al. 1998), but attempts to form realistic disc galaxies, and to test the general framework in detail, using ab initio cosmological hydrodynamic simulations, have had only limited success thus far. Specifically, the simulations have historically yielded galaxies with too large of a fraction of their baryons locked up in stars (the “overcooling problem”, e.g., Katz et al. 1996), sizes that are too small (the “angular momentum problem”, e.g., Navarro & Steinmetz 2000), and rotation curves that decline with radius (rather than being approximately flat as observed) due to the presence of an excessively large spheroidal component (e.g., Abadi et al. 2003). It is plausible that all three of these problems are related.
It has become increasingly clear from many different lines of evidence that feedback plays a very important role in regulating the star formation (SF) of galaxies and also in setting their internal structure. Modelling feedback from supernovae (SNe) in simulations of galaxy formation is challenging, as the injected thermal energy tends to be radiated away before it can have any hydrodynamical effect (e.g., Katz et al. 1996; Balogh et al. 2001). This is likely due to the fact that simulations still do not have sufficient resolution (nor all the physics necessary) to model the interstellar medium correctly. Cosmological simulations must therefore resort to subgrid recipes to try to overcome the overcooling problem (see the discussions in Dalla Vecchia & Schaye 2008, 2012), which is plausibly the fundamental underlying cause of the other problems mentioned above (e.g., Sales et al. 2010; Scannapieco et al. 2011). Typically, the SN energy is injected kinetically rather than thermally (e.g., Navarro & White 1993) and/or radiative cooling is suppressed by hand for a short time in the immediate aftermath of energy injection (e.g., Gerritsen 1997), although other solutions exist (e.g., Dalla Vecchia & Schaye 2012).
In the past few years there has been a good deal of progress on the
formation of more realistic discs galaxies in cosmological
simulations. A number of groups have been carrying out very high
resolution “zoomed” simulations of individual systems primarily
using a thermal SN feedback subgrid prescription (e.g.,
Okamoto et al. 2005, 2010; Governato et al. 2007, 2010, 2012; Stinson et al. 2010; Agertz et al. 2011; Guedes et al. 2011; Avila-Reese et al. 2011; Macciò et al. 2012; Brook et al. 2012; see also Piontek & Steinmetz 2011 for zoomed simulations with a kinetic feedback model).
Generally speaking, these studies have found that higher
numerical resolution, together with increased density thresholds for
SF, result in much more efficient feedback which, in turn, allows the
simulations to reproduce many of the properties of normal disc
galaxies, including reasonable SF efficiencies, disc sizes, and
approximately flat rotation curves. Along with these successes, however, there remain important shortcomings, such as the apparent inability of current simulations (and semi-analytic models) to match the observed specific star formation rates
While the ability of recent high resolution simulatons to reproduce a number of important observables is certainly encouraging, caution is warranted in applying the understanding derived from them to the galaxy population in general. Firstly, the feedback implemented in many of these studies has been tuned to produce ‘realistic’ disc galaxies. However, given the uncertain (and possibly degenerate) nature of subgrid processes, it is dangerous to claim success on the basis of a handful of simulated galaxies. The degree of realism should be judged on the basis of examining many different observables over a wide range of mass and redshift. Secondly, the systems that are chosen to be re-simulated at very high resolution have typically, although not always, been selected on the basis of their merger history - the common lore being that quiescent merger histories give rise to disc galaxies. However, the link between merger history and galaxy morphology is complicated (e.g., Okamoto et al. 2005; Zavala et al. 2008; Scannapieco et al. 2009; Sales et al. 2012) and, in any case, the fact that the majority of isolated normal galaxies in the local Universe are spirals implies that no such selection based on merger history should be necessary.
An alternative approach, therefore, is to try to simulate large, representative populations of galaxies without tuning the feedback or explicitly selecting particular merger histories in order to produce ‘realistic’ disc galaxies at the present day. The downside of this approach is that the resolution must be reduced, and so too are the kinds of questions one can ask of the simulations. Also, there is obviously no guarantee that reasonable disc (or elliptical) galaxies will be produced in the simulation, as there has been no explicit tuning to get the ‘right’ answer. At present, it is computationally feasible to simulate large populations () of normal galaxies with kpc scale resolution. This is the approach we adopt in the present study.
In particular, we use the Galaxies-Intergalactic Medium Interaction Calculation (gimic) suite of cosmological hydrodynamical simulations, which were carried out by the Virgo Consortium and are described in detail by Crain et al. (2009) (see also Schaye et al. 2010). The suite consists of a number of large, nearly spherical regions ( Mpc in radius) extracted from the Millennium Simulation and re-simulated at higher resolution with gas dynamics, including subgrid prescriptions for metal-dependent radiative cooling, SF, chemodynamics, and a kinetic implementation of SN feedback.
The resolution of gimic is somewhat lower than in the studies noted above that have looked at the formation of individual disc galaxies. Our work is therefore complementary to these very high-resolution studies in that, instead of aiming for maximum resolution on a single system (or a handful of systems), we elect to run lower resolution simulations of a large population of galaxies, so that we may compare our results to large observational samples of galaxies and place them in a wider context. We note that a number of recent zoom simulation-based studies have implied that very high resolution is required in order to form reasonable disc galaxies at all, but the accuracy of this statement depends crucially on the nature of the implementation of important subgrid physics (particularly SN feedback) in the simulations. As we have shown elsewhere (Crain et al. 2010; Font et al. 2011; McCarthy et al. 2012; Sales et al. 2012), the gimic simulations have no difficulty in forming disc-dominated (or even bulgeless) galaxies even though they do not have extremely high resolution (in any case, we show in the Appendix that our results are robust to changes in the resolution).
In the present study, we examine the relation between stellar mass and rotation velocity of disc galaxies (i.e., the ‘Tully-Fisher’ relation), rotation curves in bins of stellar mass of the disc galaxies, the relation between stellar mass and size of disc galaxies, the SF efficiency of both disc and elliptical galaxies, the relation between stellar mass and optical rotation to virial circular velocity ratio (), and the relation between stellar mass and specific star formation rate (sSFR). We compare the simulations to a variety of recent local observations over a wide galaxy stellar mass range of . Throughout we adjust the observational stellar mass measurements so that they are appropriate for a Chabrier (2003) stellar initial mass function (IMF) (e.g., we add 0.05 dex [subtract 0.1 dex] in the case a Kroupa [‘diet Salpeter’] IMF), which is what is adopted in the gimic simulations. For stellar masses derived using the -colour scalings of Bell et al. (2003) or Bell & de Jong (2001), we also apply a small adjustment to the mass using the expression given in the appendix of Li & White (2009) (see also Dutton et al. 2011) so that the stellar masses are consistent with those derived from more accurate ‘SED fitting’ of five-band SDSS data (Blanton & Roweis, 2007).
The present study is organised as follows. In Section 2 we present a description of the gimic simulations. In Section 3 we present our main results on the Tully-Fisher relation, rotation curves, mass-size relation, and SF efficiencies. In Section 4 we present some scaling arguments that link the faint-end slope of the galaxy stellar mass function to the slope of the Tully-Fisher relation. We summarize and discuss our findings in Section 5.
A detailed description of the gimic simulations can be found in Crain et al. (2009). We therefore present only a brief summary here.
The suite consists of a set of hydrodynamical re-simulations of five nearly spherical regions ( Mpc in radius) extracted from the (dark matter) Millennium Simulation (Springel et al., 2005). The regions were selected to have overdensities at that represent , where is the root-mean-square deviation from the mean on this spatial scale. The 5 spheres are therefore diverse in terms of the large-scale structure that is present within them. For the purposes of the present study we select ‘central’ galaxies only, which are defined as the most massive self-gravitating collection of stars, gas and dark matter (i.e., the most massive subhalo) in a friends-of-friends group. Systems and their substructures are identified in the simulations using the Subfind algorithm of Dolag et al. (2009), that extends the standard implementation of Springel et al. (2001) by including baryonic particles in the identification of self-bound substructures.
The cosmological parameters adopted for gimic are the same as those assumed in the Millennium Simulation and correspond to a CDM model with , , , (where is the rms amplitude of linearly evolved mass fluctuations on a scale of Mpc at ), km s Mpc, , (where is the spectral index of the primordial power spectrum). The value of is approximately 2-sigma higher than has been inferred from the most recent CMB data (Komatsu et al., 2011), which will affect the abundance of the simulated galaxies somewhat, but should not significantly affect their internal properties.
The simulations were evolved from to using the TreePM-SPH code gadget-3 (last described in Springel 2005), which has been modified to incorporate new baryonic physics. Radiative cooling rates for the gas are computed on an element-by-element basis using pre-computed tables generated with CLOUDY (Ferland et al., 1998) that contain cooling rates as a function of density, temperature, and redshift and that account for the presence of the cosmic microwave background and photoionisation from a Haardt & Madau (2001) ionising UV/X-Ray background (see Wiersma et al. 2009a). This background is switched on at in the simulation, where it ‘reionises’ the simulation volume.
Star formation is tracked in the simulations following the prescription of Schaye & Dalla Vecchia (2008). Gas with densities exceeding the critical density for the onset of the thermo-gravitational instability is expected to be multiphase and to form stars ( cm; Schaye 2004). Because the simulations lack the physics to model the cold interstellar gas phase, an effective equation of state (EOS) is imposed with pressure for densities
The timed release of individual elements by both massive (Type II SNe and stellar winds) and intermediate mass stars (Type Ia SNe and asymptotic giant branch stars) is included following the prescription of Wiersma et al. (2009b). A set of 11 individual elements are followed in these simulations (H, He, C, Ca, N, O, Ne, Mg, S, Si, Fe), which represent all the important species for computing radiative cooling rates, in the presence of a UV/X-ray background.
Feedback from SNe is incorporated using the kinetic wind model of Dalla Vecchia & Schaye (2008) with the initial wind velocity, , set to km/s and the mass-loading parameter (i.e., the ratio of the mass of gas given a velocity kick to that turned into newly formed star particles), , set to . This corresponds to using approximately 80% of the total energy available from SNe for a Chabrier IMF, which is assumed in the simulation. Note that unlike many previous works that use similar feedback schemes, the “kicked” particles are not hydrodynamically decoupled at any stage. This leads to efficient entrainment of the ISM by outflows which can significantly affect the properties of galaxies (see Dalla Vecchia & Schaye 2008). The above choice of feedback parameters yields a good match to the peak of the SF rate history of the universe. The simulation feedback parameters were chosen without any regard to observables other than the peak of the cosmic SF history.
To test for numerical convergence, the gimic simulations have been run at three levels of numerical resolution: ‘low’, ‘intermediate’, and ‘high’. The low-resolution simulations have the same mass resolution as the Millennium Simulation, while the intermediate- and high-resolution simulations have and times better mass resolution, respectively. Only the and volumes were run to at high resolution, owing to the computational expense of these calculations. We adopt the high-resolution simulations in our main analysis and make comparisons to the intermediate-resolution simulations to test for numerical convergence. The high-resolution runs have a dark matter particle mass M and an initial gas particle mass
We present a convergence study in the Appendix which shows that the simulation results pertaining to star formation efficiency, kinematics, and sizes are converged for systems with at least star particles. In the analysis below, we adopt a conservative lower stellar mass cut of , which corresponds to approximately star particles. We note, however, that a larger number of star particles appears to be required before converged sSFRs are obtained (see the Appendix).
Lastly, the gimic simulations are based on the SPH formalism. A relevant issue is whether differences in simulation techniques (e.g., SPH, fixed mesh, moving mesh) are important. It has been shown using idealised tests that standard SPH (which was adopted for the gimic simulations) has a number of problems relating to resolving fluid instabilities and the treatment of weak shocks (e.g., Agertz et al. 2007; Mitchell et al. 2009; Vogelsberger et al. 2012; Sijacki et al. 2012). While it would be prudent to address these issues in the future (perhaps using modifications to the SPH formalism that appear to mostly resolve these issues; e.g., Price 2012; Read & Hayfield 2012), we point out that both the OWLS (Schaye et al., 2010; Sales et al., 2010) and Aquila (Scannapieco et al., 2011) projects have clearly shown that the differences obtained using different hydrodynamic techniques are small compared to the differences in the results that arise from the use of different subgrid prescriptions.
Below we examine the stellar mass-rotation velocity (‘Tully-Fisher’) relation (Section 3.1), rotation curves in bins of stellar mass (Section 3.2), the mass-size relation (Section 3.3), the SF efficiency of the simulated galaxies (Section 3.4), the stellar mass- relation (Section 3.5), and the stellar mass-sSFR relation (Section 3.6). In each case we compare our results to a variety of recent local observations. We focus on the wide stellar mass range, (which corresponds roughly to a halo mass range of in the simulations, where is the total mass within the radius, , that encloses a mean density of 200 times the critical density). The lower stellar mass limit is motivated by the numerical convergence study in the Appendix, while the upper mass limit corresponds to the most massive galaxies in the simulation volumes. Within this mass range, there are 714 simulated galaxies in total in the GIMIC high-resolution volumes. This includes both disc- and spheroid-dominated galaxies.
When making comparisons to the observed Tully-Fisher relation, the mass-size relation, the observed rotation curves of disc galaxies, the , and the sSFRs, it is appropriate to select only disc-dominated systems from the simulations. One way to do this is kinematically; e.g., to use a cut in the fractional contribution of ordered rotation to the total kinetic energy and/or angular momentum (e.g., Abadi et al. 2003; Sales et al. 2010). However, in general it is not possible to do this kind of selection for observed galaxies. Therefore, we instead adopt a simple procedure for assigning galaxies to either the disc- or spheroid-dominated categories based on the distribution of their stellar light. For each galaxy we produce a projected azimuthally-averaged surface brightness profile
where is the projected distance from the galaxy center, is the effective radius, is the surface brightness at , is the Sersic index, and is a function of only that must be computed numerically (see, e.g., Graham et al. 2005). Finally, to approximately mimic the depth of SDSS observations, we restrict the fit to radii where mags. arcsec. The results are not sensitive to the exact surface brightness cut adopted.
We use the best-fit Sersic index, , to assign galaxies to the disc-dominated and spheroid-dominated categories. The case of corresponds to a purely exponential distribution, which is known to describe the disc component of galaxies very well, while corresponds to the classical de Vaucouleurs distribution, which has historically been used to describe the surface brightness profiles of elliptical galaxies. It has been shown by Blanton et al. (2003b) using a large sample of galaxies in the SDSS that there is an approximately bimodal distribution in at fixed galaxy luminosity (or stellar mass), in analogy to the well-known bimodal distribution of galaxy colours at fixed luminosity. Red galaxies are typically those with large Sersic indices (), while blue galaxies correspond to those with relatively small values of on average. Following Shen et al. (2003), we use as our cut for distinguishing between disc- and spheroid-dominated systems. With this cut, there are 432 disc-dominated and 282 spheroid-dominated galaxies over the range (see Fig. 4 below, which shows the distribution of with galaxy stellar mass).
3.1 Stellar mass - rotation velocity relation
In Fig. 1 we examine the relation between galaxy stellar mass and circular velocity () for disc galaxies, which is commonly referred to as the ‘Tully-Fisher’ (hereafter TF) relation. The top panel shows the relation for the case where the circular velocity is measured at the 3D radius that contains 80% of the total i-band flux
We focus first on the top panel of Fig. 1. The cyan points with error bars correspond to the recent measurements of Reyes et al. (2011). These authors selected a broad sample of galaxies from the SDSS and derived rotation curves from follow-up H observations (some of which were already presented in Pizagno et al. 2007). Although Reyes et al. (2011) did not focus specifically on disc galaxies, the requirement that the galaxies have usable H rotation curves effectively weeds ellipticals out of their sample. We use their estimates of the rotation velocity, , and their published -band magnitudes and colours to convert to stellar mass using the colour scalings of Bell et al. (2003). We select only galaxies from Reyes et al. (2011) with axis ratios in order to minimise inclination uncertainties and extinction corrections (as in Dutton et al. 2010).
Over the wide range the simulated galaxies fall approximately on top of the observed relation. To our knowledge, this is the first time a cosmological hydrodynamic simulation has produced a population of galaxies that is consistent with the observed TF relation. The agreement over this mass range is not perfect (the observed relation appears to be slightly steeper than the simulated one), but it is clearly a very large improvement over previous simulation studies that looked at representative populations (e.g., Navarro & Steinmetz 2000). We note that de Rossi et al. (2010) have also recently examined the TF relation of a population of galaxies simulated with slightly higher resolution but in smaller boxes ( Mpc on a side). Their simulated TF relation has approximately the correct slope but no direct comparison was made with observations so it is unclear whether the simulation reproduces the observed zero point as well.
At higher masses (), the simulated galaxies have circular velocities that are significantly larger than the observed rotation velocities. This is likely the result of overcooling in massive halos (indeed, see Fig. 4 below). This may signal that an increase in the efficiency of SN feedback is required at high masses and/or that another form of feedback (e.g., from Active Galactic Nuclei, hereafter AGN) is energetically important at these mass scales (e.g., Benson et al. 2003). On the other hand, de Rossi et al. (2010), who also implemented only SN feedback, did not find such a sharp change in the slope of their TF relation at high masses, although due to their smaller simulation box size they had a much smaller number of high mass galaxies in their sample. In practice, the mass scale at which overcooling begins to be an issue will depend on the details of the subgrid prescriptions and the adopted parameters, as well as on the resolution. However, as many studies have demonstrated, it is not obvious that the problems at the massive end can be solved with only feedback from star formation.
It also appears that the simulated TF relation has smaller scatter than the observed one. However, it is not clear whether the difference in scatter is a genuine discrepancy: since we have measured circular velocities, rather than gas rotation velocities (which will not be perfectly circular), the scatter in the simulated TF relation is likely to be underestimated. Particle noise (due to our finite resolution) inhibits the accuracy to which the gas rotation velocities can be determined in the simulations, which is why we have used circular velocities instead. In addition, scatter in the observed TF relation can be introduced by errors in corrections for disc inclination and dust attenuation (although we have attempted to mitigate this by selecting observed galaxies within a particular axis ratio range), as well as errors in distance measurements. A much more careful analysis, taking into account galaxy selection criteria and observational uncertainties, is required to make quantitative statements about the degree of similarity (or lack thereof) in the scatter in the observed and simulated galaxy populations.
In the bottom panel of Fig. 1 we compare the observed and
simulated TF relations using rather than . The
red points with error bars represent the 21 cm measurements of
Verheijen (2001), who measured the TF relation for a volume-limited
complete sample of spiral galaxies in the nearby Ursa Major Cluster.
We derive stellar masses for their galaxies using their -band
magnitudes and colours and the colour scalings of
Bell & de Jong (2001). We also show the best-fit power-law relation
In accordance with the results shown in the top panel of Fig. 1, we find that the simulated TF relation (using ) agrees well with the observed one over the wide range while at higher masses the simulated galaxies have too high circular velocities with respect to observed disc galaxies of the same mass.
In terms of the degree of scatter in the relation, we note that the role of distance uncertainties in the observed relation should be minimal for the Verheijen (2001) results, since their galaxies belong solely to the Ursa Major cluster. Interestingly, the scatter in the relation of Verheijen (2001) and Avila-Reese et al. (2008) appears to be similar to that of the simulations for galaxies with .
A caveat to bear in mind when comparing the simulated and observed relations, is that for the simulated galaxies is the maximum circular velocity within , whereas for the observations it is the maximum rotation velocity within the region out to which it is possible to measure gas rotation speeds (typically kpc). We have investigated what happens to the simulated - relation when we limit the maximum radius to 2 effective radii (to crudely mimic the extent of the HI gas). We find a nearly identical relationship for (the measured is the true one). However, at lower masses, the radius where the circular velocity curve peaks is beyond 2 effective radii and, consequently, the measured is lower than the true one, but only by dex on average. This has the effect of shifting the low mass simulated galaxies down so that they tend to lie in between the results of Avila-Reese et al. (2008) (and Verheijen 2001) and Bell & de Jong (2001), yielding even better agreement between the simulations and the observations.
As discussed at the beginning of Section 3, we selected the disc-dominated systems from the simulations for comparison to the observed TF relation (of disc galaxies). However, it is worthwhile to comment briefly on the behaviour of the spheroid-dominated galaxies, which are not shown in Fig. 1. Over the mass range where the disc-dominated simulated galaxies agree well with the observations (), the spheroid-dominated galaxies follow the disc-dominated trend closely. In other words, our results are insensitive to the cut in Sersic index for this mass range. It is only at the highest masses where we see a noticeable difference, in the sense that the spheroid-dominated galaxies have lower circular velocities ( and at a fixed stellar mass. However, as we will show in Section 3.3, it is clear that both disc- and spheroid-dominated galaxies at high stellar masses suffer from strong overcooling with respect to observed galaxies.
The TF relation encapsulates the rotation of a disc galaxy with a single number. However, since neither simulated nor observed galaxies have perfectly flat rotation curves, it is possible for simulations to lie on the observed TF relation without reproducing the rotation curves of real disc galaxies. Furthermore, it has been argued that the shapes of the rotation curves of normal disc galaxies are inconsistent with the predictions of the CDM model (e.g., Salucci et al. 2012), as deduced from dark matter only simulations. With this motivation, we next examine the rotation curves of observed and simulated galaxies as a function of galaxy stellar mass.
3.2 Rotation curves in bins of stellar mass
In Fig. 2 we compare observed galaxy rotation curves with predicted galaxy circular velocity curves in 9 equal-width bins of spanning the range . The continuous black curves (connecting filled black circles) represent the median circular velocity profiles of the simulated galaxies and the error bars represent the 5 and 95 percentiles. The dot-dashed red, dashed purple, and dotted orange curves represent, respectively, the median circular velocity profiles of the dark matter, stars, and gas separately. The shaded vertical region delineates the range softening lengths ( kpc), where the circular velocities of the simulated galaxies are expected to be unreliable. The cyan curves represent fits to the H rotation curves of Reyes et al. (2011), which we reproduce using the best-fit arctan parameters provided by Reyes et al. (2011). The best-fit curves are plotted out to the radius that encloses of the i-band flux.
The agreement with the amplitude and shape of the observed rotation curves over the mass range is encouraging. The observed and simulated velocity curves both show a gentle rise from kpc out to about kpc, on average, before levelling off to an approximately constant value at larger radii. We note that over this range in galaxy stellar mass, dark matter is the most important contributor to the circular velocity curves in the simulated galaxies (except in the innermost regions of the last mass bin), as can easily be seen by comparing the dot-dashed red curve (dark matter) to the dashed purple (stars) and dotted orange (gas) curves.
Strong deviations of the simulated circular velocity curves from the observed rotation curves are evident at higher stellar masses, however, in the sense that the simulated galaxies have mass distributions which are too concentrated, particularly in the very innermost regions ( kpc). This is clearly due to the stellar distribution being much too compact and dominant in these systems. It is worth pointing out that many observed massive galaxies do show evidence for declining rotation velocity curves (e.g., Persic et al. 1996; Courteau 1997; Verheijen 2001) and that such behaviour cannot be captured by the arctan functional form adopted by Reyes et al. (2011) in fitting their rotation curves. However, it is very unlikely that this accounts for most of the difference at high stellar masses. Indeed, the fact that for these high-mass galaxies there is a normalisation offset at relatively large radii, strongly suggests that these galaxies have too high SF efficiencies.
A comparison of the rotation curve results to the TF results in Section 3.1 reveals a nice consistency, in that the rotation curves look sensible over the range of stellar mass for which the simulated galaxies lie along the TF relation, and vice versa.
In terms of the spheroid-dominated galaxies in our overall sample
(which are not used in Fig. 2), there are no significant
differences between disc-dominated and spheroid-dominated galaxies
in the circular velocity
curves until high stellar masses () are
reached. At high stellar masses the circular velocity curves of the
spheroid-dominated galaxies are systematically lower in amplitude and
somewhat flatter than for the disc-dominated galaxies of the same stellar mass
3.3 Stellar mass - size relation of disc galaxies
In Fig. 3 we plot the relation between galaxy stellar mass and size. We use the radii which enclose 80% (top panel) and 50% (bottom panel) of the i-band flux in projection. The filled black circles with error bars correspond to the median and 16 and 84 percentiles in bins of galaxy stellar mass for the simulated gimic galaxies. The radii and are derived by fitting a Sersic law to the azimuthally-averaged i-band surface brightness profiles of the simulated galaxies, as described above. For comparison, we plot the median disc relations (and 1-sigma scatter) of Shen et al. (2003) (derived from SDSS data) and Baldry et al. (2012) (derived from GAMA data) in the bottom panel, and individual measurements of and from Reyes et al. (2011) in the top panel. Note that the comparison to both Shen et al. (2003) and Baldry et al. (2012) is fully consistent, in that both studies also derive by fitting Sersic laws to the azimuthally-average surface brightness profiles of the observed galaxies. Reyes et al. (2011), on the other hand, perform somewhat more complex surface brightness modelling (they fit a 2D double exponential), which we do not attempt for the simulated galaxies.
Over the range the simulated galaxies lie approximately between the median relations of Shen et al. (2003) and Baldry et al. (2012) and have significant overlap with both (bottom panel). The simulated galaxies are also consistent with data of Reyes et al. (2011) over this stellar mass range (taking into account error bars and scatter), although the trend with mass starts to deviate from the observations at masses of . The general agreement at low-moderate masses is encouraging and, together with the agreement with the rotation curves shown above, indicates that the simulated galaxies do not suffer from the well-known ‘angular momentum problem’. This is not due to resolution, as results do not change significantly if we adopt the intermediate-resolution gimic simulations instead (see the Appendix).
For the simulated galaxies are obviously too compact compared to observed galaxies of the same stellar mass. This is consistent with our findings from the Tully-Fisher relation and the rotation curves in bins of stellar mass presented above.
3.4 Stellar mass - star formation efficiency relation
We have shown that over the range of stellar mass, , the gimic high-resolution simulations approximately reproduce the observed TF relation as well as the rotation curves and galaxy sizes as a function of stellar mass. Given that for this mass range the rotation velocity in the simulations is dominated by dark matter, the agreement with the observed TF relation and rotation curves would imply that the models should have approximately the correct SF efficiency, . Here is the total current stellar mass of the galaxy (which excludes surrounding satellites), is the total mass (gas+stars+dark matter) within , and is the universal ratio of baryonic to total matter densities (Komatsu et al., 2011). The SF efficiency is therefore the fraction of the total mass accounted for by stars, relative to the maximal case where all the baryons are converted into stars. Cosmological hydrodynamic simulations have had a notoriously difficult time in forming normal galaxies with the correct , generally generating values of that are a factor of a few larger than observed (see, e.g., Guo et al. 2010, Dutton et al. 2010, and Avila-Reese et al. 2011, who compare the results of a number of recent hydrodynamic simulations with observations). This is often referred to as the ‘overcooling problem’ of cosmological simulations. How do the gimic high-resolution simulations fair in this regard?
In Fig. 4 we plot the trend between and stellar mass for the simulated galaxies at .
Each of the dots in Fig. 4 represents a simulated galaxy. The dots are colour-coded according to the fitted Sersic index. Note that we include both disc- and spheroid-dominated galaxies in this plot.
The dotted orange curve is the relation derived by Guo et al. (2010) by applying the abundance matching technique to the SDSS galaxy stellar mass function
Consistent with the results derived in Sections 3.1 and 3.2, the simulated galaxies have SF efficiencies that are quite comparable to what is derived from stacked WL/SK analyses of disc galaxies over the range in galaxy mass, . The correct fraction is achieved through efficient (but energetically feasible) SN feedback which ejects gas from the system (see, e.g., Fig. 7 of Crain et al. 2010). This is ultimately what also allows the simulations to reproduce the TF relation and the shape and normalization of the rotation curves over this mass range. We point out that unlike some recent studies of individual galaxies (e.g., Stinson et al. 2010; Macciò et al. 2012; Brook et al. 2012), we have not had to invoke multiple sources of feedback in order to reproduce the observed SF efficiency. All that has been done is to adjust the mass-loading of the SN feedback so that the simulations roughly reproduce the peak of the cosmic SF rate density of the Universe.
Interestingly, most galaxies in the simulations fall below the relation derived from applying the abundance matching technique to the galaxy stellar mass function (with or without scatter included) for . Dutton et al. (2010) previously pointed out the factor of inconsistency between the abundance matching and direct measure studies. It is presently unclear what the origin of the discrepancy is. The discrepancy cannot be explained by differences in the cosmological parameters adopted for the Millennium Simulation (from which Guo et al. 2010 derive their halo mass function) and current best estimates, as Dutton et al. (2010) have shown that the discrepancy exists for other abundance matching-based results (e.g., Moster et al. 2010; Behroozi et al. 2010) which adopted parameter values more in line with those derived from recent CMB analyses. Furthermore, the discrepancy is not simply due to the fact that the abundance matching technique is applied to the stellar mass function without regard to galaxy morphology, since disc galaxies dominate the stellar mass function at stellar masses where the discrepancy exists (thus, the abundance matching results should reflect those of disc galaxies at these masses).
However, the fact that the simulations approximately match the TF relation and velocity rotation curves over the same range of stellar mass, suggests that the systematic offset may lie mostly with the abundance matching results. Consistent with this hypothesis is the normalisation offset between the TF relation deduced from the abundance matching technique
Indeed, the agreement with the TF relation in the top panel of Fig. 1 implies that the simulated galaxies have the correct stellar mass fractions within over the range , since is a direct measure of the total mass (stars+gas+DM) within that radius (under the assumption that non-circular motions are small). Assuming that the observed and simulated galaxies have similar (dark) matter distributions beyond this radius (which is implied by the WL/SK data since these shows that an NFW profile (Navarro et al., 1996, 1997) works well at these radii), there is therefore consistency between the match to the WL/SK data in Fig. 4 and the TF relation in the top panel of Fig. 1.
A possible systematic error in the abundance matching formalism is the implicit assumption that satellite and central galaxies of fixed halo mass have the same
There is a small number of apparent outliers from the main trend at low to moderate stellar masses. We have investigated the nature of these outliers and found that they always correspond to galaxies in close proximity to a massive galaxy group. The dark matter associated with the galaxy is clearly truncated (which is why is higher than typical). This may be because these galaxies have passed through the more massive neighbour (and hence have been tidally stripped) in the past, but have since come back out to larger radii (see the population of “backsplash” subhaloes on extreme orbits identified in cosmological simulations by Gill et al. 2005; Ludlow et al. 2009).
Returning to the general trend, for stellar masses the simulations clearly have too high stellar mass fractions with respect to both the abundance matching results and the stacked measurements of disc and especially elliptical galaxies. Furthermore, the simulations do not produce a bi-modal distribution of stellar mass fractions (corresponding to discs and ellipticals) over the range as implied by stacked WL/SK measurements. It is conceivable that feedback from AGN, which was not included in our simulations and which is expected to become relatively more important at high masses (and for galaxies with relatively large spheroidal components), may resolve these issues. Indeed, McCarthy et al. (2010) have used simulations from the OverWhelmingly Large Simulations project (Schaye et al., 2010), which uses the same simulation code as gimic, to show that the inclusion of AGN feedback in cosmological simulations can reduce the stellar mass of the most massive galaxies by up to an order of magnitude (see Fig. 7 of that study), yielding agreement with observations on the mass scales corresponding to groups. Finally, we point out that current state-of-the-art very high resolution zoomed simulations (that neglect AGN feedback) also suffer from overcooling at these high masses
3.5 Stellar mass - relation
Related to the star formation efficiency, , is the ratio of galaxy rotation velocity to its virial circular velocity, which we define here (following Reyes et al. 2012; see also Dutton et al. 2010) as , where . If baryonic processes do not significantly affect the total matter distribution at small radii (i.e., so that it approximately traces an unmodified NFW profile) we would expect
In Fig. 5 we examine the relation between galaxy stellar mass and for our simulated disc galaxies. We make comparisons to the observations of Reyes et al. (2012), who measured this relation for a disc galaxy sample that is representative of the Tully-Fisher sample of Reyes et al. (2011) (see Fig. 1). Reyes et al. (2012) employed stacked galaxy-galaxy lensing to measure a mean (and thus a mean ) in bins of stellar mass. Consistent with our previous results, the simulated galaxies lie approximately on top of the observed relation for , while at higher masses the simulated galaxies have too high galaxy circular velocities with respect to the virial circular velocity of the halo in which they live.
The red dashed line in Fig. 5 represents the ratio expected in the absence of significant contraction or expansion of the inner potential well. If we restrict the comparison to the mass range , both the gimic simulations and the observations are consistent with a mild degree of expansion at low stellar masses () and a similarly mild amount of contraction at higher stellar masses (). At very high masses (), the simulated potential wells have clearly undergone too much contraction as a result of overcooling.
3.6 Stellar mass - sSFR relation
As discussed in the Introduction, there has been considerable progress in recent years in the ability of simulations to reproduce the structural and dynamical properties of normal disc galaxies at . This is due in large part to the implementation of strong feedback mechanisms in the simulations. However, it has been pointed out that in spite of this progress, there remain a number of important problems, including the apparent inability of simulations with strong feedback to reproduce the sSFRs of galaxies locally and at high redshift, particularly for low-mass galaxies (e.g., Davé et al. 2011; Avila-Reese et al. 2011; but see Brook et al. 2012). Here we examine the trend between galaxy stellar mass and sSFR at for disc galaxies in the gimic simulations.
In Fig. 6 we compare the relation between galaxy stellar mass and sSFR for disc galaxies in the gimic simulations with the observations of Salim et al. (2007). Salim et al. (2007) estimated star formation rates using UV data from GALEX for a large sample of 50,000 optically-selected galaxies from the SDSS. The black shaded region in Fig. 6 represents the galaxy mass range for which the sSFRs are sensitive to numerical resolution and are likely not converged.
At stellar masses of , the simulated galaxies roughly reproduce the observed trend, which is encouraging. Note that this (rough) agreement was not a given, as the feedback parameters in gimic have been tuned to match only the peak of the cosmic star formation rate density (at ), not the star formation rates of galaxies at the present day. Note, however, that the agreement at the highest masses () must be fortuitous, as we have already demonstrated that these galaxies suffer from overcooling (therefore the star formation rates of these galaxies are too high given their halo masses).
At lower stellar masses () the simulated galaxies fall below the observed relation of Salim et al. (2007). However, as we demonstrate in the Appendix, the sSFRs of galaxies with such low masses are rather poorly converged at in our simulations. It is possible that a further increase in resolution would result in yet higher sSFRs, possibly resolving the discrepancy at low stellar masses. Consistent with this hypothesis is the fact that Brook et al. (2012) have used very high resolution zoomed simulations of low-mass galaxies and are able to reproduce the relation of Salim et al. (2007) down to low stellar masses. Thus, it remains to be seen whether the problem of some simulations in not reproducing the sSFRs of low-mass galaxies is a fundamental one.
4 Connecting star formation efficiency, the stellar mass function, and the TF relation
In this section we briefly comment on the link between the star formation efficiency (i.e., the relation), the galaxy stellar mass function, and the TF relation.
The virial mass of a halo of mass can be expressed in terms of its virial circular velocity, , as
assuming a Hubble constant of km s Mpc (e.g., Mo et al. 1998). Thus, if we knew the relation between galaxy stellar mass, , and we could convert eqn. (2) into a TF relation.
For halos with masses below , the halo mass function is well represented by a power-law of the form
then matching abundances of haloes and galaxies
Inserting the above into eqn. (2) yields a TF relation of the form
Over the range (which avoids resolution issues at lower masses and the knee in the stellar mass function at higher masses, which would violate the power-law assumption in eqn. 4), the galaxy stellar mass function in gimic is characterised relatively well by a power-law with . From eqns. (5) and (6) this yields and over the same range in stellar mass. These expectations match what we actually measure over this range in gimic well.
It is therefore possible, at least in principle, to constrain the faint-end slope of the galaxy stellar mass function from the slope of the TF relation (or vice versa). On this basis, the observed TF relations in Fig. 1 suggest that must be relatively steep, , which is close the value of found by Baldry et al. (2012) for blue galaxies.
A rigorous statement to this effect would necessitate relaxing a number of assumptions. Firstly, we have assumed power-law dependencies for the halo and stellar mass functions over this mass range. This holds to a high level of accuracy for the halo mass function but should be checked for the stellar mass function. Secondly, we have assumed that there is no scatter in the - relation. We have explored the effects of adding log-normal scatter to an idealised powerlaw relation between and and find that it does not bias the recovered slope of the resulting stellar mass function. However, if the amplitude of the scatter were a strong function of and/or if the scatter were not log-normal, then this can result in a bias in . Finally, we have derived the TF relation in eqn. 6 in terms of the virial circular velocity, whereas what is measured is closer to . If dark matter dominates the dynamics, then to a good approximation . This breaks down when applied to a very wide range of halo masses owing to the mass-concentration relation of dark matter haloes, an effect, however, which is unimportant for the mass range under consideration here. Potentially more relevant here is that baryonic physics can either contract or expand the dark matter distribution (e.g., Duffy et al. 2010; Dutton et al. 2011) and can therefore introduce a halo mass dependence in the relationship between and . In the gimic simulations, we find these effects to be small over the mass range considered here (i.e., ).
5 Summary and Discussion
We have used the gimic simulations to study the mass distribution and SF efficiency of a large, cosmologically representative population of simulated disc galaxies over the wide range of stellar mass at . These simulations include efficient, but energetically feasible supernova feedback, but have not been tuned in any way to produce “realistic” disc galaxies. The main results of our study may be summarised as follows:
Over the stellar mass range , the simulations approximately reproduce: the observed Tully-Fisher relation (this is true for both the velocity measured at the radius that encloses 80% of the galaxy’s i-band light, , and the peak velocity, ), the observed rotation velocity curves of disc galaxies, the observed size-mass relation of disc galaxies, the relation stellar mass () and optical rotation to virial circular velocity ratio (’’), and the trend between SF efficiency and galaxy stellar mass, as inferred from stacked weak lensing and stacked satellite kinematics studies. Over this mass range, dark matter dominates the rotation curves of galaxies.
A comparison of both the observed and simulated values of ’’ with that predicted by dark matter only simulations indicates that the potential wells of galaxies with relatively low stellar masses () have undergone a relatively mild amount of expansion, while the potential wells of galaxies with higher stellar masses () have undergone a similarly mild degree of contraction (see Sec. 3.6).
For stellar masses , the simulated galaxies rotate faster than observed galaxies of the same stellar mass, are too compact, and have too large a fraction of their baryons locked up in stars (i.e., the simulations suffer from overcooling).
The simulations also reproduce the specific star formation rates of galaxies (with ) but have too low levels of SF in lower mass galaxies, which is plausibly due to the finite resolution of the simulations.
There is an intriguing discrepancy between the SF efficiency of disc galaxies derived from abundance matching and ‘direct’ measurements derived from stacked weak lensing/satellite kinematics for , which was previously pointed out by Dutton et al. (2010). The fact that our simulations simultaneously reproduce the TF relation, rotation curves in bins of stellar mass, and the SF efficiency derived from direct measurements over this mass range, suggests that the unknown systematic may lie mostly with the abundance matching results. One possibility is that the implicit assumption that satellites (at accretion) and central galaxies of fixed halo mass have the same star formation efficiency is incorrect (e.g., Rodríguez-Puebla et al. 2011.
Simple scaling arguments presented in Section 4 (see also Mo et al. 2010) show that it is possible to constrain the faint-end slope of the galaxy stellar mass function from the slope of the TF relation (or vice versa). Current measurements of the TF relation suggest that the faint slope should be relatively steep (), which is consistent with some recent measurements (e.g., Baldry et al. 2012).
The success of the gimic simulations in reproducing the observed trends at low to moderate stellar masses is directly attributable to efficient SN feedback. As shown previously by Crain et al. (2010), the SN feedback in gimic drives vigorous galactic winds, particularly at high redshift, which eject a significant fraction of the baryons from haloes with masses and pollute the intergalactic medium with metals. Thus, comparisons with observations of galactic outflows and the intergalactic medium represent promising further ways to test the model. We point out that efficient feedback was achieved without the need for extremely high resolution and is a result of the way in which the feedback was implemented in the simulations. Recent zoomed simulations in which the feedback was injected thermally and the cooling turned off temporarily (e.g., Governato et al. 2010; Stinson et al. 2010; Guedes et al. 2011; Brook et al. 2012; Macciò et al. 2012), appear to require much higher resolution before the feedback becomes similarly efficient.
For , the galaxies in the gimic simulations suffer from strong overcooling, as do all other current cosmological simulations of disc galaxy formation that invoke only SN feedback. We conjecture that another energetically important form of feedback is required at high masses, the likely candidate being that from accreting supermassive black holes (Benson et al., 2003; Di Matteo et al., 2005; Bower et al., 2006; Croton et al., 2006; Booth & Schaye, 2009). McCarthy et al. (2010) have shown, using simulations that form part of the OverWhelmingly Large Simulations project (Schaye et al., 2010), that AGN feedback is crucial in reproducing the hot gas and stellar component of systems with halo masses only slightly higher than explored here (corresponding to galaxy groups). It is therefore certainly plausible that AGN are relevant over the mass scale where gimic (and other current simulations) fails.
While the results of the present study certainly give reason for optimism, any successful model must not only reproduce these (and other) trends at but also their evolution to high redshift, which is by no means automatic. In a future study we intend to explore a wider range of observables (for both discs and ellipticals) over a wider range of redshifts and with a wider range of models, some of which include feedback from AGN.
The authors thank the referee, Vladimir Avila-Reese, for a helpful and constructive report. They also thank Qi Guo for providing her abundance matching results, Ivan Baldry for providing the GAMA stellar mass-size relation, and Aaron Dutton, Frank van den Bosch, Reina Reyes, Chris Brook, and Fabio Governato for helpful comments. IGM is supported by an STFC Advanced Fellowship at the University of Birmingham. The simulations presented here were carried out by the Virgo Consortium for Cosmological Supercomputer Simulations using the HPCx facility at the Edinburgh Parallel Computing Centre (EPCC) as part of the EC’s DEISA ‘Extreme Computing Initiative’, the Cosmology Machine at the Institute for Computational Cosmology of Durham University, and Darwin at the University of Cambridge. CSF acknowledges a Royal Society Wolfson Research Merit Award. This work was supported by Marie Curie Initial training Network CosmoComp (PITN-GA-2009-238536), ERC Advanced Investigator grant COSMIWAY and an STFC rolling grant to the Institute for Computational Cosmology. This research was also supported in part by the National Science Foundation under Grant No. NSF PHY11-2595.
Appendix: Numerical convergence
Here we investigate the sensitivity of our results to the numerical resolution. In Fig. 7 we compare the relations of the gimic intermediate- and high-resolution simulations. In general there is good agreement when there is a minimum of star particles (modulo some sensitivity to the exact mass scale where feedback stalls). We conservatively impose a cut of star particles on the high-resolution simulation (corresponding to ) in the paper.
In Fig. 8 we compare the relations of the gimic intermediate- and high-resolution simulations. Again, there is reasonably good agreement when there is a minimum of star particles, but with some differences occurring at the transition from the effective to ineffective feedback regimes.
In Fig. 9 we compare the sSFR relations of the gimic intermediate- and high-resolution simulations. There is reasonably good agreement for galaxies with . At lower masses, the sSFR rate decreases significantly with decreasing resolution.
- pagerange: Rotation rates, sizes, and star formation efficiencies of a representative population of simulated disc galaxies–LABEL:lastpage
- pubyear: 2008
- The specific star formation rate is defined as the ratio of the star formation rate to stellar mass of a galaxy.
- Gas particles are only placed on the EOS if their temperature is below K when they cross the density threshold and if their density exceeds 57.7 times the cosmic mean. These criteria prevent SF in intracluster gas and in intergalactic gas at very high redshift, respectively (Schaye & Dalla Vecchia, 2008).
- When a star particle forms, it has the same mass as the gas particle from which it formed (the gas particle is removed). Thus the initial stellar mass is typically M. However, star particles lose mass over cosmic time due to stellar evolution.
- The i-band luminosities are computed by treating each star particle as a simple stellar population (SSP). The simulations adopt a Chabrier IMF and store the age and metallicity of the particles. We use these to compute a spectral energy distribution for each star particle using the GALAXEV model of Bruzual & Charlot (2003). The i-band luminosity is obtained by integrating the product of the SED with the SDSS i-band transmission filter function. Our calculation neglects the effects of dust attenuation.
- We use the best-fit Sersic model to determine the projected radius that encloses 80% of the flux.
- Strictly speaking, this refers to the relation of Bell & de Jong (2001), where is the best-fit asymptotic velocity at large radii. For galaxies with flat rotation curves is obviously equivalent to . For galaxies with rising rotation curves, e.g., low-mass galaxies, as is determined from the observed rotation curve whereas represents an extrapolated quantity.
- Naively, one might expect spheroid-dominated, low-angular momentum galaxies to have more peaked circular velocity curves than disc-dominated galaxies. This would indeed be true for galaxies of fixed stellar mass fraction. However, at fixed stellar mass the situation is reversed, since at the high mass end disc-dominated galaxies tend to have a higher stellar mass fraction than spheroid-dominated galaxies (see Fig. 4), which is found to be the case observationally (e.g., Dutton et al. 2010). The higher stellar mass fraction leads to more peaked circular velocity curves for disc-dominated galaxies.
- The stellar masses of Li & White (2009) were derived using the five-band SDSS photometric data (and redshift) as described in Blanton & Roweis (2007) and assuming a Chabrier IMF, consistent with what is adopted in our cosmological simulations.
- We note the ‘direct’ measurements of halo mass are model-dependent, in that they assume a NFW halo profile with a particular mass-concentration relation that does not account for possible correlations between the scatter in concentration and galaxy morphology/colour.
- Guo et al. (2010) use abundance matching to derive a relation between and and use the relation between and from the Millennium Simulation to derive the relation between and . If dark matter dominates the rotation curves, then this is equivalent to the TF relation.
- Note, however, that usually abundance matching is performed for satellite galaxies at the time of accretion.
- The code used for the gimic simulations was also used in the Aquila project, which simulated a single halo. The results of that individual galaxy are completely in line with the trends present in Figs. 1 and 3 of the present study.
- Here we used the Millennium Simulation halo catalog to measure the median value of (1.16) for haloes in the mass range . Note that the ratio varies with halo mass, owing to the fact that lower mass haloes are more concentrated on average, but the dependence is very weak.
- Note that in Section 3.4 we argued there may be a systematic problem with galaxy stellar masses inferred through the abundance matching technique. However, this appears to be a problem in the normalisation, rather than in the slope (which is very similar to that derived from stacked WL/SK). Here we are interested in how the slopes of the relation, the TF relation, and the stellar mass function are related.
- Abadi, M. G., Navarro, J. F., Steinmetz, M., & Eke, V. R. 2003, ApJ, 591, 499
- Agertz, O., Moore, B., Stadel, J., et al. 2007, MNRAS, 380, 963
- Agertz, O., Teyssier, R., & Moore, B. 2011, MNRAS, 410, 1391
- Avila-Reese, V., Zavala, J., Firmani, C., & Hernández-Toledo, H. M. 2008, AJ, 136, 1340
- Avila-Reese, V., Colín, P., González-Samaniego, A., et al. 2011, ApJ, 736, 134
- Baldry, I. K., Driver, S. P., Loveday, J., et al. 2012, MNRAS, 421, 621
- Balogh, M. L., Pearce, F. R., Bower, R. G., & Kay, S. T. 2001, MNRAS, 326, 1228
- Behroozi, P. S., Conroy, C., & Wechsler, R. H. 2010, ApJ, 717, 379
- Bell, E. F., & de Jong, R. S. 2001, ApJ, 550, 212
- Bell, E. F., McIntosh, D. H., Katz, N., & Weinberg, M. D. 2003, ApJS, 149, 289
- Benson, A. J., Bower, R. G., Frenk, C. S., et al. 2003, ApJ, 599, 38
- Blanton, M. R., Hogg, D. W., Bahcall, N. A., et al. 2003a, ApJ, 592, 819
- Blanton, M. R., Hogg, D. W., Bahcall, N. A., et al. 2003b, ApJ, 594, 186
- Blanton, M. R., & Roweis, S. 2007, AJ, 133, 734
- Booth, C. M., & Schaye, J. 2009, MNRAS, 398, 53
- Bower, R. G., Benson, A. J., Malbon, R., et al. 2006, MNRAS, 370, 645
- Brook, C. B., Stinson, G., Gibson, B. K., Wadsley, J., & Quinn, T. 2012, MNRAS, 424, 1275
- Bruzual, G. Charlot, S. 2003, MNRAS, 344, 1000
- Chabrier, G. 2003, PASP, 115, 763
- Conroy, C. et al. 2007, ApJ, 654, 153
- Courteau, S. 1997, AJ, 114, 2402
- Crain, R. A. et al. 2009, MNRAS, 399, 1773 .
- Crain, R. A., McCarthy, I. G., Frenk, C. S., Theuns, T., Schaye, J. 2010, MNRAS, 407, 1403
- Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11
- Dalla Vecchia, C. Schaye, J. 2008, MNRAS, 387, 1431
- Dalla Vecchia, C., & Schaye, J. 2012, MNRAS, submitted (arXiv:1203.5667)
- Davé, R., Oppenheimer, B. D., & Finlator, K. 2011, MNRAS, 415, 11
- de Rossi, M. E., Tissera, P. B., & Pedrosa, S. E. 2010, A&A, 519, A89
- Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604
- Dolag, K., Borgani, S., Murante, G., Springel, V. 2009, MNRAS, 399, 497
- Duffy, A. R., Schaye, J., Kay, S. T., et al. 2010, MNRAS, 405, 2161
- Dutton, A. A., Conroy, C., van den Bosch, F. C., Prada, F., & More, S. 2010, MNRAS, 407, 2
- Dutton, A. A., et al. 2011, MNRAS, 416, 322
- Fall, S. M., & Efstathiou, G. 1980, MNRAS, 193, 189
- Ferland, G. J., Korista, K. T., Verner, D. A., Ferguson, J. W., Kingdon, J. B., Verner, E. M. 1998, PASP, 110, 761
- Font, A. S., McCarthy, I. G., Crain, R. A., Theuns, T., Schaye, J., Wiersma, R. P. C., Dalla Vecchia, C. 2011, MNRAS, 416, 2802
- Gerritsen, J. P. E. 1997, Ph.D. Thesis
- Gill, S. P. D., Knebe, A., & Gibson, B. K. 2005, MNRAS, 356, 1327
- Governato, F., Willman, B., Mayer, L., et al. 2007, MNRAS, 374, 1479
- Governato, F., Brook, C., Mayer, L., et al. 2010, Nature, 463, 203
- Governato, F., Zolotov, A., Pontzen, A., et al. 2012, MNRAS, in press (arXiv:1202.0554)
- Graham, A. W., Driver, S. P., Petrosian, V., et al. 2005, AJ, 130, 1535
- Guedes, J., Callegari, S., Madau, P., & Mayer, L. 2011, ApJ, 742, 76
- Guo, Q., White, S., Li, C., & Boylan-Kolchin, M. 2010, MNRAS, 404, 1111
- Haardt, F. & Madau, P. 2001, in Clusters of galaxies and the high redshift universe observed in X-rays, ed. D.M. Neumann & J.T.T. Van.
- Jenkins, A., Frenk, C. S., White, S. D. M., et al. 2001, MNRAS, 321, 372
- Katz, N., Weinberg, D. H., & Hernquist, L. 1996, ApJS, 105, 19
- Kennicutt, R. C., Jr. 1998, ApJ, 498, 541
- Klypin, A., & Prada, F. 2009, ApJ, 690, 1488
- Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18
- Li, C., & White, S. D. M. 2009, MNRAS, 398, 2177
- Ludlow, A. D., Navarro, J. F., Springel, V., et al. 2009, ApJ, 692, 931
- Macciò, A. V., Stinson, G., Brook, C. B., et al. 2012, ApJ, 744, L9
- Mandelbaum, R., Seljak, U., Kauffmann, G., Hirata, C. M., & Brinkmann, J. 2006, MNRAS, 368, 715
- Mandelbaum, R., Seljak, U., & Hirata, C. M. 2008, JCAP, 8, 6
- McCarthy, I. G., Schaye, J., Ponman, T. J., et al. 2010, MNRAS, 406, 822
- McCarthy, I. G., Font, A. S., Crain, R. A., et al. 2012, MNRAS, 420, 2245
- Mitchell, N. L., McCarthy, I. G., Bower, R. G., Theuns, T., & Crain, R. A. 2009, MNRAS, 395, 180
- Mo, H. J., Mao, S., & White, S. D. M. 1998, MNRAS, 295, 319
- Mo, H., van den Bosch, F. C., & White, S. 2010, Galaxy Formation and Evolution. Cambridge University Press, 2010. ISBN: 9780521857932
- More, S., van den Bosch, F. C., Cacciato, M., et al. 2011, MNRAS, 410, 210
- Moster, B. P., Somerville, R. S., Maulbetsch, C., et al. 2010, ApJ, 710, 903
- Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563
- Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
- Navarro, J. F., & Steinmetz, M. 2000, ApJ, 538, 477
- Navarro, J. F., & White, S. D. M. 1993, MNRAS, 265, 271
- Okamoto, T., Eke, V. R., Frenk, C. S., & Jenkins, A. 2005, MNRAS, 363, 1299
- Okamoto, T., Frenk, C. S., Jenkins, A., & Theuns, T. 2010, MNRAS, 406, 208
- Persic, M., Salucci, P., & Stel, F. 1996, MNRAS, 281, 27
- Piontek, F., & Steinmetz, M. 2011, MNRAS, 410, 2625
- Pizagno, J., et al. 2007, AJ, 134, 945
- Price, D. J. 2012, MNRAS, 420, L33
- Read, J. I., & Hayfield, T. 2012, MNRAS, submitted (arXiv:1111.6985)
- Reed, D. S., Bower, R., Frenk, C. S., Jenkins, A., & Theuns, T. 2007, MNRAS, 374, 2
- Reyes, R., Mandelbaum, R., Gunn, J. E., Pizagno, J., & Lackner, C. N. 2011, MNRAS, 417, 2347
- Reyes, R., Mandelbaum, R., Gunn, J. E., et al. 2012, MNRAS, in press (arXiv:1110.4107)
- Rodríguez-Puebla, A., Avila-Reese, V., Firmani, C., & Colín, P. 2011, RMxAA, 47, 235
- Sales, L. V., Navarro, J. F., Schaye, J., Vecchia, C. D., Springel, V., & Booth, C. M. 2010, MNRAS, 409, 1541
- Sales, L. V., et al. 2012, MNRAS, submitted (arXiv:1112.2220)
- Salim, S., Rich, R. M., Charlot, S., et al. 2007, ApJS, 173, 267
- Salucci, P., Wilkinson, M. I., Walker, M. G., et al. 2012, MNRAS, 420, 2034
- Scannapieco, C., White, S. D. M., Springel, V., & Tissera, P. B. 2009, MNRAS, 396, 696
- Scannapieco, C., et al. 2011, MNRAS, submitted (arXiv:1112.0315)
- Schaye, J. 2004, ApJ, 609, 667
- Schaye, J., Dalla Vecchia, C. 2008, MNRAS, 383, 1210
- Schaye, J., et al. 2010, MNRAS, 402, 1536
- Schulz, A. E., Mandelbaum, R., & Padmanabhan, N. 2010, MNRAS, 408, 1463
- Shen, S., Mo, H. J., White, S. D. M., et al. 2003, MNRAS, 343, 978
- Sijacki, D., Vogelsberger, M., Keres, D., Springel, V., & Hernquist, L. 2012, MNRAS, submitted (arXiv:1109.3468)
- Spergel, D. N., et al. 2007, ApJS, 170, 377
- Springel, V., White, S. D. M., Tormen, G., Kauffman, G. 2001, MNRAS, 328, 726
- Springel, V. 2005, MNRAS, 364, 1105
- Springel, V. et al. 2005, Nature, 435, 629
- Stinson, G. S., Bailin, J., Couchman, H., et al. 2010, MNRAS, 408, 812
- Verheijen, M. A. W. 2001, ApJ, 563, 694
- Vogelsberger, M., Sijacki, D., Keres, D., Springel, V., & Hernquist, L. 2012, MNRAS, submitted (arXiv:1109.1281)
- White, S. D. M., & Frenk, C. S. 1991, ApJ, 379, 52
- Wiersma, R. P. C., Schaye, J., Smith, B. D. 2009a, MNRAS, 393, 99
- Wiersma, R. P. C., Schaye, J., Theuns, T., Dalla Vecchia, C., Tornatore, L. 2009b, MNRAS, 399, 574
- Zavala, J., Okamoto, T., & Frenk, C. S. 2008, MNRAS, 387, 364