The low abundance and insignificance of dark discs in simulated Milky Way galaxies
We investigate the presence and importance of dark matter discs in a sample of 24 simulated Milky Way galaxies in the apostle project, part of the eagle programme of hydrodynamic simulations in Lambda-CDM cosmology. It has been suggested that a dark disc in the Milky Way may boost the dark matter density and modify the velocity modulus relative to a smooth halo at the position of the Sun, with ramifications for direct detection experiments. From a kinematic decomposition of the dark matter and a real space analysis of all 24 halos, we find that only one of the simulated Milky Way analogues has a detectable dark disc component. This unique event was caused by a merger at late time with an LMC-mass satellite at very low grazing angle. Considering that even this rare scenario only enhances the dark matter density at the solar radius by 35% and affects the high energy tail of the dark matter velocity distribution by less than 1%, we conclude that the presence of a dark disc in the Milky Way is unlikely, and is very unlikely to have a significant effect on direct detection experiments.
keywords:cosmology: theory, dark matter – Galaxy: structure, disc – methods: numerical
The very successful CDM cosmological model is based on the assumption that around 25% of the energy density of the Universe is in the form of as-yet undetected, weakly-interacting particles that make up the dark matter (DM). Validating this assumption requires detecting the particles either indirectly through their decay or annihilation products, or directly through interaction with the atoms of a detector, or by finding evidence of their existence in particle accelerators (see Bertone et al. (2005); Bertone (2010) for reviews). In the case of direct detection experiments knowledge of the local (solar neighbourhood) DM density and its velocity distribution is essential. This is particularly important for experiments sensitive to low mass particles for which the energy required to interact with an atom in the detector can only be reached by energetic particles in the high-velocity tail of the distribution.
Historically, the distribution of the DM velocity modulus has been characterised by a Maxwellian with a peak value of ; the local dark matter density is normally taken to be . High resolution N-body simulations of halos of Milky Way type galaxies (Springel et al., 2008) have since confirmed this general picture except that the DM velocity distribution is anisotropic and better described by a multivariate Gaussian rather than a Maxwellian. These simulations have also revealed halo-to-halo variations in the velocity modulus typically of amplitude which are related to the assembly history of individual halos (Vogelsberger et al., 2009).
While high resolution -body simulations have fully characterized all the relevant properties for direct detection experiments, they ignore effects due to galaxy formation which could, in principle, modify this picture. For example, the contraction of the inner halo induced by the condensation of baryons towards the centre steepens the density profile, shifting the velocity distribution towards larger values. (Barnes & White, 1984; Blumenthal et al., 1986; Gnedin et al., 2004)
Another potentially important baryon effect is the formation of a “dark disc” perhaps facilitated by the formation of the baryonic disc (Read et al., 2008). Such a dark disc could form by the accretion of a massive satellite in the plane of the galactic disc at late times (Read et al., 2009). The presence of the galactic disc modifies the potential well (relative to the DM-only case) and this may help confine dark matter stripped from the satellite onto co-planar orbits, boosting the creation of a dark disc co-rotating with the stars. If a large fraction of the DM at the position of the Sun were co-rotating with the stellar disc, fewer particles would have a large velocity in the detector frame, thus modifying the outcome of DM direct detection experiments, potentially to a very significant extent. It is therefore of central importance for these detection experiments to quantify how common and how massive such discs are in CDM.
By comparing the morphological and kinematic properties of the Milky Way to idealized simulations of accretion events onto discs Purcell et al. (2009) concluded that the co-rotating dark matter fraction near the Sun is at most 30%. A hydrodynamical “zoom” simulation of a galaxy was found to have formed a dark disc contributing around of the DM at the solar radius co-rotating with the stars (Ling et al., 2010). On the other hand, using a chemodynamical template Ruchti et al. (2014) found no evidence for accreted stars near the Sun, leaving little room for a dark disc in our galaxy.
In this paper we search for dark discs in the twelve apostle (“A Project Of Simulating The Local Environment”) simulations (Sawala et al., 2016; Fattahi et al., 2016)) carried out as part of the “Evolution and Assembly of GaLaxies and their Environments” (eagle) programme (Schaye et al., 2015; Crain et al., 2015). The apostle simulations are designed to reproduce the kinematic properties of the Local Group. They thus have merger histories resembling that of our own Galaxy and should provide informative predictions for the abundance of dark discs. These simulations have been shown to reproduce the satellite galaxy luminosity functions of the MW and Andromeda and do not suffer from the “too-big-to-fail” problem (Boylan-Kolchin et al., 2011). They were performed with the same code used in the eagle simulation, which provides an excellent match to many observed properties of the galaxy population as a whole, such as the stellar mass function and the size distribution at low and high redshift.
Of particular relevance for this study are the rotation curves of the eagle galaxies, which agree remarkably well with observations of field galaxies (Schaller et al., 2015a). This indicates that the matter distribution in the simulated galaxies is realistic and suggests that baryon effects on the DM are appropriately modelled. The same set of simulations have been used to make predictions for indirect detection experiments (Schaller et al., 2016). One aim of this Letter is to help complete a consistent picture for both types of DM experiments from state-of-the-art simulations.
We assume the best-fitting flat CDM cosmology to the WMAP7 microwave background radiation data (Komatsu et al., 2011) (, , and ), and express all quantities without factors. We assume a distance from the Galactic Centre to the Sun of and a solar azimuthal velocity of .
2 Simulations and method
Details of the eagle code used to carry out the apostle simulations used in this work may be found in (Schaye et al., 2015; Crain et al., 2015). The apostle simulations are described by Sawala et al. (2016). Here we briefly describe the parts of the model that are most relevant to dark discs; in this section we also describe our procedure for identifying dark discs.
2.1 Simulation setup and subgrid model
The eagle simulation code is built upon the gadget code infrastructure (Springel, 2005). Gravitational interactions are computed using a Tree-PM scheme, and gas physics using a pressure-entropy formulation of smooth particle hydrodynamics (SPH) (Hopkins, 2013), called anarchy (Dalla Vecchia (in prep.), see also Schaller et al. (2015c)). The astrophysical subgrid model includes the following processes: a star formation prescription that reproduces the Kennicutt-Schmidt relation (Schaye & Dalla Vecchia, 2008), injection of thermal energy and metals in the ISM following Wiersma et al. (2009b), element-by-element radiative cooling (Wiersma et al., 2009a), stellar feedback in the form of thermal energy injection (Dalla Vecchia & Schaye, 2012), supermassive black hole growth and mergers and corresponding AGN feedback (Booth & Schaye, 2009; Schaye et al., 2015; Rosas-Guevara et al., 2015). Galactic winds develop naturally without imposing a preferred direction or a shutdown of cooling.
The free parameters of the model were calibrated (mainly by adjusting the efficiency of stellar feedback and the accretion rate onto black holes) so as to reproduce the observed present-day galaxy stellar mass function and observed relation between galaxy masses and sizes, as well as the correlation between stellar masses and central black hole masses (Schaye et al., 2015; Crain et al., 2015). Galaxies are identified as the stellar and gaseous components of subhalos found using the subfind algorithm (Springel et al., 2001; Dolag et al., 2009). No changes to the original model used by Schaye et al. (2015) were made to match the Local Group properties (Sawala et al., 2016).
2.2 Milky Way halo selection
The 24 galactic halos analyzed in this study come from twelve zoom resimulations of regions extracted from a parent N-body simulation of a volume containing particles. Halo pairs were selected to match the observed dynamical properties of the Local Group (Fattahi et al., 2016). In each volume a pair of halos with mass in the range that will host analogues of the MW and Andromeda galaxies are found. We use the two halos in each volume to construct our sample. These selection criteria ensure that our sample is not biased towards particular halo assembly histories and consists of galaxies with a similar environment to that of the MW and thus, plausibly, with a relevant star formation history. This selection contrasts with that by Read et al. (2009) where three simulated galaxies “were chosen to span a range of interesting merger histories”.
The high-resolution region of our simulations always encloses a sphere larger than centred on the centre of mass of the Local Group at . The primordial gas particle mass in the high resolution regions was set to and the DM particle mass to ; the Plummer-equivalent softening length for all particle types is (physical). DM-only simulations were run from the same initial conditions and are denoted as dmo in the remainder of this letter.
The galaxies that formed in our halos have a stellar mass in the range to and half-mass radii in the range to in reasonable agreement with observational estimates for the MW (e.g. Bovy & Rix, 2013). The Bulge-over-Total ratios of our galaxies, obtained from a kinematical decomposition, vary from to . The four haloes that formed an elliptical galaxies are kept in the sample for comparison purposes.
2.3 Stellar and dark matter velocity distributions
For each galaxy in the sample, we define a North pole axis to be in the direction of the angular momentum vector of all the stars within a spherical aperture of around the centre of potential of the halo. Note that, as in the lower-resolution study of eagle galaxies by Schaller et al. (2015b), there is no significant offset between the centre of the stellar and DM distributions. We then select matter in a torus, along the plane of the stellar disc, around the solar radius, , with both radial and vertical extents of . These tori contain of the order of DM and star particles, allowing the density and velocity distributions to be well sampled. In the case of the dmo simulations, we place the torus in the same plane as in the corresponding full baryonic run in order to make sure that any bulk halo rotation is accounted for. We note that the alternative choice of aligning the planes with the inertial axes of the halos leads to qualitatively similar results. Similarly, shifting the radius of the torus to smaller or larger values does not change the results of our study. We also generated randomly oriented tori within which we compute the mean DM density.
Within each galactic plane torus we calculate the velocity distributions of the DM and stars in the radial, azimuthal and vertical directions in bins of width of . The tori are split in a large number of angular segments to bootstrap-resample the distributions and construct an estimate of the statistical fluctuations induced by the finite particle sampling of our simulations. An example (halo , see below) of the velocity distributions for the DM and stars is shown in Fig. 1.
In this section we formulate a criterion for identifying dark discs in the simulations using the velocity distribution (as in previous studies) and apply it to our sample of galaxies. We then perform an analysis of the spatial distribution to confirm our findings. Finally, we analyse the cases that present tentative evidence for a dark disc.
3.1 Azimuthal velocity distribution models
In order to quantify the prominence of a DM disc in velocity space, we fit three different models to the azimuthal DM velocity distribution in the simulations.
Model 1: a single Gaussian, centred at with the root mean square value (rms) as the only free parameter;
Model 2: A single Gaussian, with both the centre and rms as free parameters;
Model 3: Two Gaussians, one centred at and the other one at location that is free to vary. The rms of both Gaussians, as well as their relative normalisation, are the other three free parameters.
The results of collsionless simulations are well described by the first model (Vogelsberger et al., 2009), and this also applies to the azimuthal velocity distributions extracted from our dmo halos. The best-fitting Gaussian from Model 2 is shown as a dashed black line on the bottom left panel of Fig. 1. Finally, the third model, shown as two dotted green lines also in the bottom left panel of the figure, is the model used by Read et al. (2009), which includes a halo component (as in Model 1) and a disc component represented by the second Gaussian.
The azimuthal velocity distribution of a halo with a significant dark disc would either have a single-peaked distribution (Model 2) shifted to a mean velocity comparable to the typical stellar azimuthal velocity at the radius of the Sun, or require a clear second Gaussian in addition to the halo component (Model 3).
For each galaxy we find the best-fitting parameters for all three models of the DM azimuthal velocity distributions. We then use the Akaike Information Criterion (AIC: Akaike, 1998), corrected for finite sample size (Burnham & Anderson, 2002), to select amongst the different models111For two models with and , the relative likelihood of the two models given the data is given by . The one with significantly lowest AIC, for a given halo, minimises information loss and should be favoured 222Using the Bayesian Information Criterion (e.g. Schwarz (1978)) leads to similar results..
3.2 Abundance of dark discs from velocity space analysis
The top panel of Fig. 2 shows the AIC of all three models for all halos in our sample. For of our halos, either the model with two Gaussians (Model 3, red squares) is disfavoured or all three models are too close for a decision to be made. Only halo 10 (discussed further below) is clearly better modelled by two Gaussians.
For all other halos, the simple model consisting of a single Gaussian either centred at (Model 1, blue triangles) or at an adjustable value (Model 2, yellow circles) are either slightly or strongly favoured by the velocity distributions. For those cases where the preferred mean velocity is non-zero, to establish whether the off-centre Gaussians are located at large values of , which would also indicate the presence of a dark disc, we show the position of the mean of the Gaussian in the central panel of Fig. 2. As can be seen, ignoring halo 10, the location of the centre of the shifted Gaussian (Model 2, yellow circles) is always below . This implies that baryon effects have displaced the peak of the Gaussian by less than the typical rms value of Model 1 (top dashed blue line, located at ), thus rendering this shift inconsequential for DM direct detection experiments. This shifted Gaussian centre is also located at less than , implying that the rotation identified by the model is not commensurate with the rotation speed of the stars, one of the criteria signalling the presence of a dark disc. Note also that if these slightly shifted Gaussians were to represent genuine dark discs, we would infer the presence of a counter-rotating component in halo 2. We also don’t find any correlation between the dark disc AIC values and galaxy properties such as mass, size or Bulge-to-Total ratio. We conclude that in of our representative MW halos, no dark discs able to affect the tail of the velocity distribution are detected.
3.3 Abundance of dark discs from real space analysis
If a dark disc can be identified in velocity space, it should also be identifiable in real space. To verify this, we constructed randomly orientated tori of the same size as the original torus aligned with the stellar disc plane. We compute the DM density in each of them and construct a distribution of densities whose medians and percentiles are displayed on the bottom panel of Fig. 2. We then determine the location in this distribution of the DM density corresponding to the torus aligned with the plane of the stellar disc. The presence of a dark disc should manifest as an enhanced DM density in the aligned torus compared to the other tori. However, as galaxies (especially centrals and late types) are preferentially aligned with their halo (e.g. van den Bosch et al., 2002; Bett et al., 2010; Velliscig et al., 2015), the dark matter density in the plane of the stellar disc is likely to be enhanced as a result of the anisotropic dark matter distribution even in the absence of a dark disc.
For most halos in our sample, the DM density is indeed enhanced in the plane of the stellar disc (see bottom panel of Fig. 2), but this enhancement is also present in the dmo simulations. This implies that the stellar disc that formed in this halo is aligned with the halo and not that a dark disc has formed. Only halo shows signs of a dark disc: its dmo counterpart displays a DM density lower than the median in the plane of the stellar disc but the apostle simulation displays a vastly enhanced DM density (more than above the median density). This high density is not the result of an intrinsic alignment of baryons with the halo but rather evidence for the presence of a dark disc.
Examining the variation of density with torus orientation shows that only one out of our MW halos (halo 10 again) shows signs of the presence of a dark disc. This result confirms our earlier conclusion from the velocity analysis.
3.4 Detailed analysis of the galaxy with a dark disc
We now carry out a more detailed analysis of halo 10, the only one333This is the second halo in volume AP-5 of Sawala et al. (2016). Detailed properties of this simulation volume and halo can be found in Table 2 of Fattahi et al. (2016). This galaxy has a stellar mass of and a half-mass radius of , well within our sample. With a Bulge-to-Total ratio of , it is one of the most disk-dominated system of our sample. for which a double Gaussian azimuthal velocity distribution (Model 3) is favoured and the only one for which the real space analysis also implies the presence of a significant overdensity in the plane of the stellar disc. The distributions of each velocity component for this halo are displayed in Fig. 1. For the azimuthal distribution, the best-fitting models 2 and 3 are plotted. The double Gaussian model is clearly favoured. The second Gaussian is centred on , a displacement comparable to the rms of the best-fitting Gaussian distribution for the corresponding dmo halo torus velocity distribution (). Baryonic effects have induced a clear second peak in the distribution of .
Apart from the position of the velocity peak, the other important property of a dark disc is the amount of dark matter that it contains. This can be characterized in different ways. Simply evaluating the integral of the two Gaussians required to fit the velocity distribution implies that of the total DM mass at the position of the Sun is in the component rotating with respect to the halo. Another way of measuring the mass of the rotating component is to apply a kinematical bulge/disc decomposition technique (e.g. Abadi et al. (2003); Scannapieco et al. (2012)) whereby the orbital energies are compared to the energy of the equivalent circular orbit of the same radius. The bulge (in our case the halo component) will be distributed symmetrically around (if the halo is not rotating, or slightly displaced from if it is) and the disc will appear as a clear second peak. Comparing the mass in the halo component to the total shows that of the mass is non-rotating, in agreement with the simpler estimate given above.
Finally, we find that the DM density in the plane of the disc in halo 10 is . Compared to the median of all the tori, this again indicates an excess of , in agreement with the kinematically derived dark disc mass fraction.
It is interesting to note that despite containing roughly a third of the mass at the location of the Sun, the dark disc has only a minimal impact on the distribution of the velocity modulus (top left panel of Fig. 1), which remains essentially unchanged from the dmo case (blue line). The velocity distribution of this halo, as well as its expected signal in direct detection experiments, lie well within the halo-to-halo scatter measured by Bozorgnia et al. (2016) in their sample of simulated MW analogues which includes some of the apostle simulations. More quantitatively, the DM mass with a velocity magnitude larger than has increased by only in the baryonic simulation compared to its dmo counterpart. We conclude that the one dark disc that has formed in our 24 simulations would have a minimal impact on direct dark matter detection experiments.
3.5 Origin of the dark disc
To investigate the origin of the dark disc that formed in one of our simulations we trace back in time the particles that have the largest azimuthal velocities today. We find that these particles belonged to a subhalo that merged with the central galaxy at . At the time of the merger, the subhalo had a mass of and the galaxy in it a stellar mass of , comparable to the Large Magellanic Cloud. The satellite impacted the galaxy at an angle of only above the plane of the stellar disc. Most of its material is then tidally stripped in the plane of the disc forming of a stream. The formation of a dark disc in this galaxy is thus consistent with the formation mechanism proposed by Read et al. (2009).
We searched for dark discs in a sample of 24 Milky Way analogues simulated as part of the apostle project of simulations of volumes selected to match the kinematical and dynamical properties of the Local Group. This environment is similar to that in which the Milky Way formed suggesting that the galaxies in our sample are representative of plausible formation paths for the Milky Way.
We find that only one out the 24 cases develops a dark disc aligned with the plane of the stellar disc. The dark disc was identified by fitting models with and without a disc to the azimuthal velocity distribution at the radial location of the Sun and selecting the best model according to the AIC. The identification was then confirmed by searching for unusual dark matter overdensities in the plane of the stellar disc (comparing to randomly oriented discs). None of the other 23 halos show any evidence for a dark disc. We conclude from our unbiased sample of MW halo analogues that dark discs are rare.
Our simulations reveal that rather unusual conditions are required for the formation of a dark disc. In our case, the disc resulted from a recent impact, at a very low grazing angle, of a satellite as large as the LMC. According to Purcell et al. (2009) and Ruchti et al. (2014), Milky Way kinematical data indicate that our galaxy could not have experienced an encounter of this kind.
For the dark disc that formed in our simulations we found that of the mass in a torus at the location of the Sun is rotating at a mean velocity of . However, this rotating dark matter component does not significantly modify the distribution of the velocity modulus measured in the counterpart of the halo in a dark matter only simulation. The azimuthal velocity distribution is still well fit by a Maxwellian, indicating that this disc would have a negligible impact on direct dark matter detection experiments.
We conclude that while evidence for recent satellite collisions with the Milky Way disc would be very interesting to find in, for example, the GAIA data, direct dark matter searches need not be concerned about the potentially confusing effects of a dark disc at the position of the Earth.
This work would have not be possible without Lydia Heck and Peter
Draper’s technical support and expertise. We thank Nassim Bozorgnia,
Francesca Calore and Gianfranco Bertone for useful discussions on DM
direct detection experiments.
This work was supported by the Science and Technology Facilities Council (grant number ST/L00075X/1) and the European Research Council (grant numbers GA 267291 “Cosmiway” ).
This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/H008519/1, and STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure.
- Abadi et al. (2003) Abadi M. G., Navarro J. F., Steinmetz M., Eke V. R., 2003, ApJ, 597, 21
- Akaike (1998) Akaike H., 1998, in Parzen E., Tanabe K., Kitagawa G., eds, Springer Series in Statistics, Selected Papers of Hirotugu Akaike. Springer New York, pp 199–213, doi:10.1007/978-1-4612-1694-0_15
- Barnes & White (1984) Barnes J., White S. D. M., 1984, MNRAS, 211, 753
- Bertone (2010) Bertone G., 2010, Particle dark matter: Observations, models and searches. Cambridge University Press
- Bertone et al. (2005) Bertone G., Hooper D., Silk J., 2005, Phys. Rep., 405, 279
- Bett et al. (2010) Bett P., Eke V., Frenk C. S., Jenkins A., Okamoto T., 2010, MNRAS, 404, 1137
- Blumenthal et al. (1986) Blumenthal G. R., Faber S. M., Flores R., Primack J. R., 1986, ApJ, 301, 27
- Booth & Schaye (2009) Booth C. M., Schaye J., 2009, MNRAS, 398, 53
- Bovy & Rix (2013) Bovy J., Rix H.-W., 2013, ApJ, 779, 115
- Boylan-Kolchin et al. (2011) Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2011, MNRAS, 415, L40
- Bozorgnia et al. (2016) Bozorgnia N., et al., 2016, preprint, (arXiv:1601.04707)
- Burnham & Anderson (2002) Burnham K. P., Anderson D. R., 2002, Model selection and multimodel inference: a practical information-theoretic approach. Springer Science & Business Media
- Crain et al. (2015) Crain R. A., et al., 2015, MNRAS, 450, 1937
- Dalla Vecchia & Schaye (2012) Dalla Vecchia C., Schaye J., 2012, MNRAS, 426, 140
- Dolag et al. (2009) Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497
- Fattahi et al. (2016) Fattahi A., et al., 2016, MNRAS, 457, 844
- Gnedin et al. (2004) Gnedin O. Y., Kravtsov A. V., Klypin A. A., Nagai D., 2004, ApJ, 616, 16
- Hopkins (2013) Hopkins P. F., 2013, MNRAS, 428, 2840
- Komatsu et al. (2011) Komatsu E., et al., 2011, ApJS, 192, 18
- Ling et al. (2010) Ling F.-S., Nezri E., Athanassoula E., Teyssier R., 2010, J. Cosmology Astropart. Phys., 2, 12
- Purcell et al. (2009) Purcell C. W., Bullock J. S., Kaplinghat M., 2009, ApJ, 703, 2275
- Read et al. (2008) Read J. I., Lake G., Agertz O., Debattista V. P., 2008, MNRAS, 389, 1041
- Read et al. (2009) Read J. I., Mayer L., Brooks A. M., Governato F., Lake G., 2009, MNRAS, 397, 44
- Rosas-Guevara et al. (2015) Rosas-Guevara Y. M., et al., 2015, MNRAS, 454, 1038
- Ruchti et al. (2014) Ruchti G. R., Read J. I., Feltzing S., Pipino A., Bensby T., 2014, MNRAS, 444, 515
- Sawala et al. (2016) Sawala T., et al., 2016, MNRAS, 457, 1931
- Scannapieco et al. (2012) Scannapieco C., et al., 2012, MNRAS, 423, 1726
- Schaller et al. (2015a) Schaller M., et al., 2015a, MNRAS, 451, 1247
- Schaller et al. (2015b) Schaller M., Robertson A., Massey R., Bower R. G., Eke V. R., 2015b, MNRAS, 453, L58
- Schaller et al. (2015c) Schaller M., Dalla Vecchia C., Schaye J., Bower R. G., Theuns T., Crain R. A., Furlong M., McCarthy I. G., 2015c, MNRAS, 454, 2277
- Schaller et al. (2016) Schaller M., et al., 2016, MNRAS, 455, 4442
- Schaye & Dalla Vecchia (2008) Schaye J., Dalla Vecchia C., 2008, MNRAS, 383, 1210
- Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
- Schwarz (1978) Schwarz G., 1978, The Annals of Statistics, 6, 461
- Springel (2005) Springel V., 2005, MNRAS, 364, 1105
- Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
- Springel et al. (2008) Springel V., et al., 2008, MNRAS, 391, 1685
- Velliscig et al. (2015) Velliscig M., et al., 2015, MNRAS, 453, 721
- Vogelsberger et al. (2009) Vogelsberger M., et al., 2009, MNRAS, 395, 797
- Wiersma et al. (2009a) Wiersma R. P. C., Schaye J., Smith B. D., 2009a, MNRAS, 393, 99
- Wiersma et al. (2009b) Wiersma R. P. C., Schaye J., Theuns T., Dalla Vecchia C., Tornatore L., 2009b, MNRAS, 399, 574
- van den Bosch et al. (2002) van den Bosch F. C., Abel T., Croft R. A. C., Hernquist L., White S. D. M., 2002, ApJ, 576, 21