Statistics of Milky Way-mass halos

There’s no place like home? Statistics of Milky Way-mass dark matter halos

Michael Boylan-Kolchin, Volker Springel, Simon D. M. White and Adrian Jenkins
Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany
Institute for Computational Cosmology, Department of Physics, University of Durham, South Road, Durham DH1 3LE, UK
e-mail: mrbk@mpa-garching.mpg.de
Abstract

We present an analysis of the distribution of structural properties for Milky Way-mass halos in the Millennium-II Simulation (MS-II). This simulation of structure formation within the standard CDM cosmology contains thousands of Milky Way-mass halos and has sufficient resolution to properly resolve many subhalos per host. It thus provides a major improvement in the statistical power available to explore the distribution of internal structure for halos of this mass. In addition, the MS-II contains lower resolution versions of the Aquarius Project halos, allowing us to compare our results to simulations of six halos at a much higher resolution. We study the distributions of mass assembly histories, of subhalo mass functions and accretion times, and of merger and stripping histories for subhalos capable of impacting disks at the centers of halos. We show that subhalo abundances are not well-described by Poisson statistics at low mass, but rather are dominated by intrinsic scatter. Using the masses of subhalos at infall and the abundance-matching assumption, there is less than a 10% chance that a Milky Way halo with will host two galaxies as bright as the Magellanic Clouds. This probability rises to for a halo with . The statistics relevant for disk heating are very sensitive to the mass range that is considered relevant. Mergers with infall mass : redshift zero virial mass greater than 1:30 could well impact a central galactic disk and are a near inevitability since , whereas only half of all halos have had a merger with infall mass : redshift zero virial mass greater than 1:10 over this same period.

keywords:
methods: -body simulations – cosmology: theory – galaxies: halos
pagerange: There’s no place like home? Statistics of Milky Way-mass dark matter halosApubyear: 2010

1 Introduction

Our understanding of galaxy formation is strongly influenced by our own Galaxy. Several of the apparent points of tension between observations and the predictions of the now-standard Cold Dark Matter (CDM) structure formation model – for example, the abundance of dark matter subhalos in comparison to the relative paucity of observed Local Group satellites and the existence of the thin Galactic disk in light of the ubiquity of dark matter halo mergers – are most clearly demonstrated from Milky Way data. While it is clearly an excellent testing ground for a wide array of phenomena related to galaxy formation, the Milky Way (MW) is nonetheless a single object and may not be representative in its properties. This is particularly true with respect to infrequent, stochastic phenomena such as major mergers: while CDM robustly predicts the mean number of mergers onto halos of a given mass, the halo-to-halo scatter is large.

The situation is similar for cosmological re-simulations of individual MW-mass dark matter halos. In the quest to understand the Milky Way’s formation, numericists have been simulating individual MW-mass dark matter halos at increasingly high resolution. By investing substantial computational resources, the most recent generation has reached one billion particles per halo, with particle masses (Diemand et al., 2008; Springel et al., 2008; Stadel et al., 2009). Only a small number of such simulations can be performed and there is no guarantee that the simulated halos will be “typical.”

In order to explore whether various characteristics of the Milky Way (e.g., its satellite galaxy population, thin disk, or merger history) are consistent with CDM expectations, we must understand the statistical properties predicted for MW-mass halos. The advent of modern galaxy surveys has enabled such an undertaking observationally: the Sloan Digital Sky Survey (SDSS; York et al. 2000) has provided a huge sample of galaxies with stellar masses similar to that of the MW, although only mean halo properties averaged over large numbers of similar galaxies can readily be measured (Mandelbaum et al., 2006; More et al., 2009). Predicting the formation histories and redshift zero properties of the halos of a representative population of MW-mass galaxies requires numerical simulations that combine large volume with high mass resolution. Resolving the internal structure of ten subhalos per host halo requires approximately one hundred thousand particles per host, which is not trivial when combined with the goal of obtaining a statistical sample of MW-mass halos. This would, for example, require particles in the volume probed by the Millennium Simulation (MS; Springel et al. 2005), a 100-fold increase over the particle number actually used.

The only computationally feasible approach at present is to use a smaller simulation volume combined with higher mass resolution. In this paper, we use the Millennium-II Simulation (MS-II; Boylan-Kolchin et al. 2009), which was designed to meet these requirements. It allows us to investigate statistical properties of MW-mass halos, including their assembly histories and their subhalo content. Throughout this work, we also compare results for our representative sample to the six halos from the Aquarius simulation series of Springel et al. (2008). These cosmological resimulations of Milky Way-mass dark matter halos have approximately one thousand times better mass resolution than the MS-II (see Section 2.2).

The basic properties of the MS-II and of the sample of MW-mass dark matter halos we investigate are given in Section 2, together with a convergence study for subhalos. Section 3 discusses the assembly history of MW-mass halos in terms both of their virial mass and of their central gravitational potential, and relates these to their concentrations and angular momenta. Section 4 investigates the subhalo population of MW-mass halos at , including the subhalo mass function and whether its scatter is Poissonian, the distribution of masses for the most massive subhalo in each host, and typical accretion times for subhalos as a function of their current mass. In Section 5, we explore the merger histories of MW-mass halos, focusing on mergers that are likely to result in galaxy-galaxy mergers. Section 6 contains a discussion of our results and their implications, while Section 7 gives our conclusions. Some properties of satellite subhalos are quantified in the Appendix. Throughout this paper, all logarithms without a specified base are natural logarithms.

2 Numerical Simulations

2.1 The Millennium-II Simulation

The Millennium-II Simulation is a very large cosmological -body simulation of structure formation in the CDM cosmology. The simulation is described fully in Boylan-Kolchin et al. (2009). For completeness, we review some of its salient features here.

The MS-II follows billion particles, each of mass , from redshift to in a periodic cube with side length . The Plummer-equivalent force softening of the simulation is and was kept constant in comoving units for the duration of the simulation.

The cosmological parameters used in the MS-II are identical to those adopted for the Millennium Simulation and the Aquarius simulations:

(1)

where is the Hubble constant at redshift zero in units of , is the rms amplitude of linear mass fluctuations in spheres at , and is the spectral index of the primordial power spectrum. While the values of these parameters were originally chosen (for the MS) for consistency with results from the combination of WMAP 1-year and 2dF Galaxy Redshift Survey data (Spergel et al., 2003; Colless et al., 2001), some are now only marginally consistent with more recent analyses. In particular, and are now the preferred values based on WMAP 5-year results combined with baryon acoustic oscillation and Type Ia supernova data (Komatsu et al., 2009). For the most part, these differences are unimportant for the results presented in this paper: the abundance of halos in the WMAP5 and Millennium-II cosmologies differs by less than 7% at , for example. A discussion about possible effects of varying the assumed cosmology is contained in Sec. 6.3.

In addition to the raw simulation output of six phase-space coordinates and one identification number per particle, stored at 68 outputs spaced according to equation 2 of Boylan-Kolchin et al. (2009), additional information about dark matter structures and merging histories in the MS-II was computed later and saved.111Dark matter halo catalogs (both FOF and SUBFIND) and merger trees from the MS-II are publicly available at
http://www.mpa-garching.mpg.de/galform/millennium-II

Friends-of-Friends (FOF) groups were identified on-the-fly during the simulation using a linking length of . All groups at each time with at least 20 particles were stored in halo catalogs. FOF groups were searched for bound substructure using the SUBFIND algorithm (Springel et al., 2001). All substructures containing at least 20 bound particles were deemed to be physical subhalos and were stored in subhalo catalogs at each snapshot. Each subhalo in the MS-II was assigned a mass equal to the sum of its constituent particles; this is the mass for subhalos we use throughout this paper. All subhalos also have a maximum circular velocity , defined via

(2)

The main component of a FOF halo is generally a dominant subhalo222This is sometimes referred to as the “central” subhalo and corresponds closely to the standard picture of a smooth, virialized dark matter halo whose particles make up most of the bound component of the FOF group. FOF halos can also contain non-dominant subhalos (“satellites”); these typically contain in total only a small fraction of the mass of the dominant subhalo, even in the highest resolution simulations currently possible (Gao et al., 2004; Springel et al., 2008; Diemand et al., 2008; Stadel et al., 2009).

Each FOF halo containing at least one subhalo has an associated virial mass and virial radius, defined as the mass and radius of a sphere containing an average density times the critical value Throughout this paper, we use , which is derived from the spherical top-hat collapse model (e.g., Gunn & Gott 1972; Bryan & Norman 1998); we denote the virial mass and the virial radius . Other common choices include (corresponding to and ) and [corresponding to and ]. At , for the cosmology in equation (2.1) and for any given halo. For dominant subhalos, is comparable to, but usually slightly less than, .

Merger trees were constructed at the subhalo level by requiring each subhalo to have at most one unique descendant. Descendants were determined by considering all subhalos at snapshot containing at least one particle from a given subhalo at snapshot and computing a weighted score based on the particles’ binding energies from . The descendant-finding algorithm also checks possible descendants at to account for cases where a subhalo temporarily disappears from the catalogs, usually because it passes very near another, more massive subhalo. Note that the merger trees are not constructed to explicitly follow mergers of FOF groups but rather to follow mergers of individual subhalos; the merger of two subhalos is often delayed relative to the merger of their original host FOF groups by several dynamical times. Such subhalo mergers correspond much more closely to the intuitive notion of disjoint objects merging to form a monolithic descendant than do mergers of FOF groups.

2.2 The Aquarius Simulations

The Aquarius Project (Springel et al., 2008) is a suite of cosmological resimulations of Milky Way-mass dark matter halos. Six halos, denoted Aq-A through Aq-F, were selected at from a cosmological simulation known as hMS (Gao et al., 2008). The hMS is a lower resolution version of the MS-II: it follows particles in the same volume, with the same cosmological parameters, and with the same amplitudes and phases333The amplitudes and phases are identical for all perturbations with wavenumber smaller than the Nyquist frequency of the hMS; higher wavenumber modes in the MS-II were randomly sampled from the same underlying power spectrum for the initial power spectrum as the MS-II. The halos were selected randomly from the full set of Milky Way-mass halos in the hMS simulation subject only to a weak isolation criterion: the nearest halo with mass greater than half that of the candidate halo was required to be at least from the candidate halo at . This isolation criterion is not particularly restrictive: only 22% of the halos from the mass-selected sample in the MS-II described below fail to meet this requirement. Once selected, the halos were re-simulated at up to five different levels of resolution. The highest level, which used a particle mass of and a force softening of parsecs, led to over 1.5 billion particles within at . This resolution was used for only one of the six halos but all six were resimulated at the next resolution level, for which is at most . These are the “level 2” simulations that we compare to the MS-II halos in this paper. Subhalo catalogs and merger trees analogous to those from the MS-II were also created for the Aquarius simulations. Further details of the Aquarius Project are presented in Springel et al. (2008).

2.3 Milky Way-mass halos

Recent estimates of the mass of the Galaxy’s dark matter halo range from (Wilkinson & Evans, 1999; Sakamoto et al., 2003; Battaglia et al., 2005; Dehnen et al., 2006; Xue et al., 2008; Li & White, 2008). To bracket these estimated values of and to assess any possible trends with halo mass, we select all halos with from the redshift zero output of the MS-II. The resulting sample of 7642 total “Milky Way-mass” halos forms the basis of the analysis in this paper. With the mass resolution of the MS-II – – this selection criterion means that the MW-mass halos considered here are represented by 50,000-500,000 particles. The typical virial radius of these halos is , which is 200 times larger than the Plummer-equivalent force softening of the MS-II.

An alternative to selecting halos by mass is instead to select by , which is more observationally accessible. Such a selection is not without complications, however, as the presence of a disk and bulge is likely to significantly alter in real galaxies. In the models of Klypin et al. (2002), for example, the adiabatic response of the dark matter halo to the formation of the MW causes an increase in of between 19% and 100%. The virial mass of the halo should not be modified by these effects.

2.4 Subhalos

Figure 1: A comparison of various relations for satellite subhalos with at redshift zero. Top: the median and relation for subhalos of Milky Way-mass halos (black) and for subhalos of all halos (red) from the MS-II. The best power law fit to the median relation is plotted in blue. Bottom: A direct comparison of for the six level 2 Aquarius halos (filled black circles) and for the corresponding six halos in the MS-II (open magenta squares); the best-fit relation from the upper panel is also plotted as a blue dashed line. The two data sets agree well for , while the MS-II points lie systematically above the Aquarius points at lower due to mass resolution effects. The six central halos are also included in the bottom panel (upper right) and are offset from the satellite subhalo relation because they are not affected by tidal stripping.

Halo substructure is a direct result of the hierarchical nature of halo formation in the CDM cosmology. Properties of subhalos such as their abundance, spatial distribution, and internal structure can provide an important link between simulations and observations, since (at least some) subhalos are expected to be the hosts of satellite galaxies. One should consider only subhalos that are reliably resolved, however, as subhalo structure can be modified by numerical effects at low resolution. The upper panel of Fig. 1 shows the relationship between at (denoted ) and at () for all satellite subhalos in our MW-mass halo sample that lie within distance from the center of their host halos (black lines); the three lines show the , , and percentiles of the distribution. Also plotted is the same relation for all host halos in the MS-II (red dashed lines). The Milky Way sample agrees very well with the full simulation sample.

The upper panel of Fig. 1 also shows the best power law relation fitted over the range (blue dashed curve) of . The relation follows the power law fit closely for , while there is a strong break at low . This is due to mass resolution, as is confirmed by the lower panel of Fig. 1. This panel shows the same relation but now for all six of the level 2 Aquarius halos (black circles) and for the corresponding six halos in the MS-II (open magenta squares). While the level 2 Aquarius data follow the power-law relation down to the smallest plotted [], the MS-II data start to differ systematically at . Note that the host halos follow a relationship with a similar slope but with larger mass at fixed . This is in part because central subhalos are not affected by tidal stripping. In addition, the way SUBFIND assigns particles to bound structures, starting with all particles bound to the dominant subhalo, contributes to the offset. Given the results of Fig. 1, we adopt 150 particles (), corresponding to , as the minimum number of particles for which we consider a subhalo well enough resolved to estimate . This corresponds to a cut where the relation deviates by about 10% from the power law fit in the upper panel of Fig. 1.

Subhalos in high resolution simulations can lose substantial mass as they are tidally stripped and gravitationally shocked within the potential of their hosts (e.g., Tormen et al. 1998; Klypin et al. 1999; Ghigna et al. 2000; Hayashi et al. 2003; Kravtsov et al. 2004b; Gao et al. 2004; Boylan-Kolchin et al. 2009). The putative galaxies within these subhalos should be much more resilient to stripping and disruption both because they sit in the densest regions of dark matter halos and because the condensation of baryons to the centers of dark matter halos deepens the gravitational potential (White & Rees, 1978).

This differential stripping, affecting dark matter more strongly than stellar components, means that the relationship between the stellar content and the instantaneous dark matter properties of a subhalo is generally very complicated. It is much more likely that the mass (or ) a subhalo had when it was still an independent halo (i.e., a dominant subhalo in a FOF group) can be related to its stellar content in a straightforward manner. Indeed, several studies have shown that many features of the observed galaxy distribution can be reproduced in -body simulations simply by assuming that the stellar mass of a galaxy is a monotonically increasing function of the maximum dark matter mass ever attained by the subhalo that surrounds it. This “abundance matching” assumption (Vale & Ostriker, 2004) can reproduce observed galaxy clustering and its evolution with redshift (Conroy et al., 2006; Brown et al., 2008; Conroy & Wechsler, 2009; Wetzel & White, 2010; Guo et al., 2010).

Since it is useful to know when a given subhalo first became a satellite and what its mass was at that time, we trace each subhalo back in time and find the point at which its bound mass (determined by SUBFIND) reached a maximum. We denote this redshift by and the corresponding subhalo mass and maximum circular velocity as and . The subscript “acc” is used because subhalos usually reach their maximum mass immediately before being accreted into a more massive object; in this sense, represents an accretion redshift.444Even though the criterion defining is that the subhalo’s bound mass is maximized over its history, we avoid the subscript “max” so as not to cause any confusion with the peak in the circular velocity curve , which is defined at all redshifts for which a subhalo can be identified This allows us to study quantities such as the distribution of and the average mass loss of satellite subhalos. Relationships among , , and for subhalos are presented in the Appendix.

3 Main Halos

In this section, we study some properties of the main halos themselves. Sec. 3.1 explores the assembly histories of Milky Way mass halos, both in terms of and . In Sec. 3.2, we investigate the distribution of halo concentrations. The focus of Sec. 3.3 is on halo spin parameters. Throughout this section, we compare properties of the six Aquarius halos to properties of the full MS-II sample of “Milky Way” halos.

3.1 Halo Assembly

Figure 2: Top: Mass assembly histories of Milky Way-mass halos. The median MAH is shown as a dashed curve, while the gray shaded region contains 68% of the distribution. Also plotted are the MAHs for the MS-II versions of the six Aquarius halos. The Aquarius halos provide a fairly representative sample of MAHs; halos A and C have quiescent merger histories, while halo F has a recent major merger. Bottom: accretion histories from the MS-II. As in the upper panel, the median relation is plotted as a dashed curve while the region is shown as the gray shaded region. The earliest-forming Aquarius halos in terms of MAHs – halos A and C – also assemble earliest as judged by , with changing by less than 25% from to . Halo F, on the other hand, undergoes a major merger at that changes by 25%.

The mass assembly histories (MAHs) of dark matter halos have been studied extensively via -body simulations and extended Press-Schechter (1974) theory (Lacey & Cole, 1993, 1994; Wechsler et al., 2002; van den Bosch, 2002; Zhao et al., 2003b, a; Tasitsiomi et al., 2004; Cohn & White, 2005; Neistein et al., 2006; Li et al., 2008; Zhao et al., 2009; McBride et al., 2009). As a direct result of this work, many basic properties of MAHs have been established. For example, Wechsler et al. (2002) found that MAHs of individual halos follow a one-parameter family . The parameter is directly related to the typical formation time for halos of mass and is an increasing function of mass. Tasitsiomi et al. (2004) showed that this exponential form was unable to fit individual cluster-mass MAHs in some cases and that a more general function, , provided a better match on a case-by-case basis. McBride et al. (2009) confirmed and extended this result over a wide range in halo masses.

The median MAH from the full MW sample is plotted as the dashed line in Fig. 2; the gray shaded region marks the region for the MAHs. A purely exponential form is a poor fit to the median MAH, while the modified form of Tasitsiomi et al. with and provides a fit to within 10% for . The median MAH can be fitted extremely well by

(3)

with and : this fit matches the measured median MAH to within 1.6% for . The solid lines in Fig. 2 show the MAHs for the Aquarius halos measured from the MS-II.555Figure 13 of Boylan-Kolchin et al. (2009) shows that there is excellent agreement between the MAHs of the level 2 Aquarius simulations and the corresponding halos in the MS-II. The Aquarius halos sample the full distribution of MAHs: halos A and C form early, halo F forms late, and halos D and E track the median MAH. Halo B trails the mean MAH at early times but catches up at .

While the MAHs are normalized to unity at , the scatter in grows fairly quickly with redshift, reaching a factor of 2.5 at and 4 at . The scatter in at fixed is nearly constant for . As a result, formation times defined with respect to fixed fractions of the final mass – i.e., – will all have similar scatter in , approximately 0.25, for .

While characterizing the growth of for a halo that reaches MW-like masses at is useful, it is far from a complete description of that halo’s assembly history. Many studies have shown that halo growth occurs inside-out, with the central regions built up first and the outskirts assembled later (e.g., Fukushige & Makino 2001; Loeb & Peebles 2003; Zhao et al. 2003b; Diemand et al. 2007; Cuesta et al. 2008). Furthermore, defining halo boundaries with respect to means that a halo’s virial mass will increase with time, even in the absence of any physical accretion, simply due to the re-definition of the boundary. A halo can therefore grow substantially in without any appreciable change in the mass at small radii, where the main baryonic component should lie. Accordingly, it is useful to find a way to characterize the build-up of the central regions of halos in addition to studying the growth of .

To this end, we consider the growth of in the bottom panel of Fig. 2. For the sample of Milky Way-mass halos considered here, is on average a factor of 6.1 smaller than , so probes the mass distribution on a scales times smaller than . Additionally, the gravitational potential energy per unit mass of a halo is proportional to , so studying the evolution in probes the growth of the dark matter halo’s central potential. The dashed curve in the bottom panel of Fig. 2 show the median relation while the shaded region shows the scatter. For redshifts , the median value of can be approximated to within 1.3% as

(4)

Qualitatively, the trend in is the same as in : both grow rapidly at early times and slowly at late times. At a more detailed level, however, the two show important differences in growth. – and consequently, the central potential – approaches its redshift zero value much earlier than : if we define a formation redshift as the time when reaches half of its present value, then the median is approximately 4, in comparison with . This is a reminder that the central potential of a halo is set much earlier than its virial mass. The scatter in at fixed is noticeably larger than at fixed , showing that the spread in formation times is larger if these are defined using rather than .

As in the upper panel of Fig. 2, the build-up of for the Aquarius halos (the colored lines) spans the range given by the full sample of MW-mass halos. Furthermore, the behavior of mass growth for individual halos is matched by that of . Halos A and C also assemble their central potentials at an early time, while F forms much later.

A stable central potential is likely to be conducive to disk galaxy formation and evolution, as potential fluctuations may drive “secular” evolution and help transform disks into bulges. Halos A-E have nearly constant central potentials from to (in the case of C, from ), indicating they are possible disk galaxy hosts. Halo F, which experiences a late major merger, would likely host a spheroid-dominated galaxy. These expectations are generally borne out by SPH simulations of the Aquarius halos by Scannapieco et al. (2009). Using a dark matter mass resolution 2 to 3 times better than that of the MS-II, they find that halos C, D, and E all have well-defined disks and that halos A and B also contain (smaller) disks, while halo F has a spheroid with no disk at all.

3.2 Halo Concentrations

Figure 3: The probability distribution of halo concentrations (solid black line), along with the best log-normal fit (dashed line) with mean and standard deviation . Values for the individual Aquarius halos in the MS-II are shown as colored vertical lines. See the text for a discussion of how the concentrations are computed.

While numerical simulations have shown that dark matter halos have a nearly universal density profile (Navarro et al., 1996, 1997, hereafter NFW), the radial scale of the profile does not follow the virial scaling of . An additional parameter is therefore required to specify the profile of the typical halo of a given ; a common choice is halo concentration , defined as the ratio of the virial radius to the radius at which the logarithmic slope of the density profile reaches . The scatter in at fixed halo mass is related to the diversity of halo formation histories. Previous studies have shown that average halo concentrations decrease weakly with halo mass – – with a fairly large scatter at fixed mass, (Navarro et al., 1997; Jing, 2000; Bullock et al., 2001b; Neto et al., 2007; Gao et al., 2008; Macciò et al., 2008; Zhao et al., 2009).

We can estimate a concentration parameter for each halo in our MW-mass sample by assuming that each halo’s density structure can be fitted with an NFW profile; we then have and therefore, . The measured values for the Aquarius halos in the MS-II differ from those in the level 2 Aquarius simulations by less than 10%, showing this radius to be accurately determined at MS-II resolution. The probability distribution function for from our MW-mass halo sample is plotted as a solid black line in Fig. 3. The dashed line shows the best-fitting normal distribution with and . The dark (light) shaded region marks the () range of this distribution. There is a slight but noticeable excess of halos at low concentrations compared to the best-fitting log-normal distribution; Neto et al. (2007) have shown that this is due to a population of unrelaxed halos.

Vertical colored lines mark the concentrations of the Aquarius halos in the MS-II, computed as described above. Halos B, D, E, and F lie near the peak of the distribution while halos A and C lie in the high concentration tail (though see the caveat about halo A below). This behavior is expected given the results of Fig. 2: the concentration of a halo is (loosely) related to the density of the universe when the halo formed. Early-forming halos such as A and C are more concentrated because the universe was denser when they formed.

It is important to note that we did not fit the density profiles of these halos to a specific parametric form such as NFW when estimating the concentrations presented here. Instead, and were measured directly from the simulation data using SUBFIND. Navarro et al. (2010) have made a thorough study of the structure of the Aquarius halos at , concluding that their density profiles are not perfectly universal and therefore cannot be fully described by a model with two scale parameters but no shape parameter (such as NFW). In agreement with Neto et al. (2007) and Gao et al. (2008), they found significantly better fits with the three parameter model of Einasto (1965), although the mass range spanned by the Aquarius halos is too small to see the systematic redshift dependence of the shape parameter that was measured in the earlier work. We emphasize that profile fitting typically attempts to represent the widest possible range in radius but tends to focus on the inner regions, usually starting at the minimum resolvable radius. This is in part because inner regions of halos are more relaxed and also more regular than the outskirts. For example, Navarro et al. show that the Aquarius A-1 halo is fitted extremely well by an Einasto profile over the range but that at larger radii, its density profile rises above this model due to substructure and recently accreted material. The latter is a substantial contribution to the mass and leads to an increase in the virial radius of the halo relative to an extrapolation of the profile fitted to the inner region. This is not the case for any of the other halos (see fig. 3 of Navarro et al.) In this sense, the measured concentration of the A halo is larger than it “should be.”

3.3 Halo Spin

Figure 4: The probability distribution for the spin parameter for Milky Way-mass halos in the MS-II (solid black curve), along with regions (shaded dark and light gray, respectively). Dashed black line: best-fitting log-normal distribution, with mean and standard deviation . Dashed cyan line: fit to the distribution given in equation (6) with and .

The angular momentum content of a galaxy is often assumed to be connected to the spin of its dark matter halo. In the commonly-adopted model of Mo, Mao, & White (1998), the specific angular momentum of a galactic disk is a fixed fraction of that of its host halo. Many more recent versions of the Mo, Mao, and White model (e.g., Dutton et al. 2007) modify this assumption only slightly. Quantifying the angular momentum content of dark matter halos may therefore help to connect dissipationless -body simulations to the properties of disk galaxies.

A convenient parametrization of a dark matter halo’s angular momentum is the modified spin parameter (Bullock et al., 2001a), defined as

(5)

where is the specific angular momentum of the halo. The distribution of spin parameters for dark matter halos as a function of their mass has been studied by several authors (Barnes & Efstathiou, 1987; Warren et al., 1992; Cole & Lacey, 1996; Lemson & Kauffmann, 1999; Bullock et al., 2001a; Bailin & Steinmetz, 2005; Macciò et al., 2007; Bett et al., 2007; Macciò et al., 2008). One of the most recent, and extensive, studies is that of Bett et al. (2007), who showed that the distribution of halo spins in the MS is well-fitted by the function

(6)

In Figure 4, we show the probability distribution function of for the MS-II MW-mass halos (solid curve), with and regions given by the dark and light shaded regions, respectively. The best-fitting log-normal distribution has and and is shown as a dashed black line. The dashed cyan line shows a fit to equation (6) with and ; this provides a significantly better description of the data than the log-normal fit.

The values for all of the Aquarius halos in the MS-II are shown as vertical tick marks in Fig. 4. While all of the halos lie within the region, five of the six lie below the median . The lone exception, halo F, had a recent major merger (see Fig. 2). Recent major merger remnants tend to have spin parameters that are higher than typical (Vitvitska et al., 2002; Maller et al., 2002). The spin parameter of these halos usually decreases after the major merger due to further accretion, which brings in additional non-aligned angular momentum. It therefore seems that at least five, and possibly all six, of the Aquarius halos lie in the lower half of the distribution.

4 subhalos

Figure 5: The differential subhalo mass function. Data from the 2039 MW-mass halos in the MS-II with are shown as the black data points, while the best-fitting relation using equation (8) is shown as a dotted black line. Results are plotted both for redshift zero subhalo masses (, lower points) and for masses at accretion (, upper points; offset 0.5 dex vertically for clarity). Extrapolation of the MS-II results to low agrees well with results from the individual level 2 Aquarius simulations (colored lines).

The abundance of subhalos in MW-mass halos has been the topic of fierce debate since late in the 1990’s, when cosmological simulations were finally able to reliably resolve dark matter substructure within galaxy-scale halos. These studies predicted that MW-mass halos should have hundreds of dark matter subhalos with , while only of order ten such satellite galaxies were seen in the Local Group (Klypin et al., 1999; Moore et al., 1999). Since then, the discovery of several faint Local Group satellites in SDSS data (e.g., Belokurov et al. 2007) has extended the range of observed satellite luminosities by almost four orders of magnitude, to , roughly doubling the observational sample. On the other hand, the dynamic range of dark matter simulations has increased by an additional three decades in mass and the number of resolved dark matter subhalos has increased by a factor of approximately one thousand: Springel et al. (2008) find subhalos with inside of the Aq-A-1 simulation at (see also Diemand et al. 2008 and Stadel et al. 2009). It is unambiguous that the observed number of luminous satellites with falls far short of the expected number of dark matter subhalos above the same threshold for standard CDM models. The resolution of this discrepancy likely lies in the strong dependence of galaxy formation physics on (see Kravtsov 2010 for a recent review) combined with the limited sky coverage and the surface brightness detection threshold of the SDSS (Tollerud et al., 2008; Walsh et al., 2009). Substructure in galaxy-mass halos is therefore still of great interest for clarifying the properties both of dark matter and of galaxy formation at low mass (and low star formation efficiency). In this section, we explore the subhalo content of the MW-mass halo sample from the MS-II at .

In the analysis that follows, the dominant subhalo is excluded in all cases. Furthermore, we consider only subhalos lying within ; from their host; this provides consistency when comparing halos of different host mass. It is very important to note that using a different limiting radius such as or will strongly affect subhalo abundance statistics because the radial distribution of subhalos is heavily weighted toward large radii (De Lucia et al., 2004; Gao et al., 2004; Springel et al., 2008). This choice also affects statistics measuring the dynamical evolution of subhalos, as subhalos falling into their host for the first time experience tidal stripping that depends on the local density and therefore, on the pericenters of their orbits about the halo.

4.1 Subhalo mass function

Figure 6: The cumulative subhalo mass function for MW-mass halos having . Results are shown both in terms of (lower solid line) and (upper solid line). The dotted curves show equation (7) using the parameters determined by fitting the differential mass functions (Fig. 5).

As shown by several authors (e.g., Gao et al. 2004; Giocoli et al. 2008; Angulo et al. 2009), the abundance of subhalos with a given fraction of their host’s mass varies systematically, with more massive hosts having more subhalos. This reflects the later formation times of more massive halos, which leave their subhalos less time to be tidally stripped. In the MS-II, the normalization of the differential mass function increases by approximately 15% for each decade in host mass. In order to maximize our resolution, in this section we use only the 2039 MS-II host halos with when computing the mass function.

The differential mass function of (non-dominant) subhalos is shown for this sample in Fig. 5. Results are plotted as black points both for the actual redshift zero masses (; lower data points) and also for the accretion masses (; upper points, offset vertically by 0.5 dex for clarity). For comparison, subhalo mass functions obtained from the individual level 2 Aquarius simulations (thick colored lines) are also plotted in Fig. 5.

We assume that the cumulative mass function can be represented by a power law with an exponential cut-off at high masses:

(7)

The differential mass function then also has a simple functional form:

(8)

Equation (8) provides an excellent fit to the MS-II data for , , , and . This fit is plotted as a dotted line in Fig. 5. Although the range of overlap between the MS-II and Aquarius data is relatively small, extrapolating the fit from the MS-II data to low agrees extremely well with the Aquarius data all the way to the Aquarius level 2 resolution limit, .

Figure 7: The cumulative subhalo velocity function for MW-mass halos having (solid line), plotted in terms of . The dotted line shows the analog of equation (7), determined by fitting the differential velocity function; the best-fitting parameters are given in the figure.

We also use equation (8) to fit the subhalo mass function computed in terms of . The comparison of MS and MS-II results in Guo et al. (2010) suggests that per subhalo is required to obtain converged results for . We are therefore able to probe only for . This limited range makes the determination of the slope using the MS-II data alone nearly impossible. Giocoli et al. (2008) have suggested that this slope is close to that of the redshift zero cumulative mass function, and we also fit with equation (8) holding fixed to -0.935, the value obtained from fitting the differential mass function for . A good fit can be obtained with , , and and is plotted as the upper dotted line in Fig. 5. Again, the extrapolation of the MS-II fit to low masses agrees very well with subhalo mass functions computed directly from the level 2 Aquarius simulations (upper set of solid curves). This shows that our adopted value of is indeed appropriate.

The corresponding cumulative subhalo mass functions are plotted in Fig. 6 for redshift zero subhalo masses (, lower solid curve) and for subhalo masses at accretion (, upper solid curve). The dotted lines show equation (7) with parameters taken directly from fits to the differential mass functions (i.e., we do not fit the cumulative mass functions independently). These curves are excellent representations of the data for occupation numbers . At lower , the presence of a very small number of ongoing major mergers results in an excess cumulative abundance compared to the fit for .

The abundance of subhalos as a function of is a related quantity, and one that is often used when comparing simulation data to observations because is less affected than by the dynamical evolution of a subhalo within its host. Fig. 7 shows the cumulative abundance of subhalos, , in terms of (solid curve). It exhibits the same behavior as : a power law at low with an exponential cut-off at high . By fitting the differential velocity function (to ensure that errors in each bin are uncorrelated) to the same functional form as in equation (8) and converting the fit to the cumulative velocity function, we find that the slope of is -2.98 at low and that at ; this fit is shown as the dotted line in Fig. 7. This low- slope is in good agreement with results from the individual, high resolution Aquarius and Via Lactea simulations.

4.2 Scatter in the subhalo mass function

The halo-to-halo scatter about the mean relation in equation (7) is a matter of considerable interest for several reasons. A number of groups have invested considerable computational resources into simulating individual MW-mass halos at extremely high resolution. State-of-the-art simulations currently use approximately one billion particles to resolve the mass distribution within (Diemand et al., 2008; Springel et al., 2008; Stadel et al., 2009). Such calculations are extremely expensive, however, so only a small number can be performed. In order to interpret their results, we must know how much scatter is expected among halos of the same mass. This scatter is also a fundamental input parameter for halo occupation distribution (HOD) models (e.g., Ma & Fry 2000; Peacock & Smith 2000; Seljak 2000; Scoccimarro et al. 2001; Berlind & Weinberg 2002).

Figure 8: The ratio of the measured dispersion in relative to that expected from a Poisson distribution with mean value equal to . The dispersion is larger than Poissonian and the ratio of the two increases systematically to lower . The dashed gray line shows a model in which the variance is a linear sum of Poisson fluctuations and an 18% fractional intrinsic scatter [equation (9)], while the shaded gray region encompasses the same model with intrinsic scatter between 14 and 22%.

We find that the scatter in the cumulative mass function at fixed can be well-modeled by postulating that the variance in is the sum of two contributions, one due to Poisson fluctuations () and one intrinsic ():

(9)

This is demonstrated in Fig. 8. The solid curves show as a function of for host halos of different masses, from Milky Way-mass hosts to halos corresponding to rich galaxy clusters. We have used the MS for host halos with , as the volume of the MS-II does not provide sufficient statistics at this mass scale; for all other masses we use the MS-II. The dashed gray line shows the model of equation (9) with a fractional intrinsic scatter of 18%, while the shaded gray region corresponds to varying from 14% to 22%.

This analysis shows that the scatter in the subhalo mass function is nearly Poissonian for massive subhalos (, corresponding to occupation numbers of ) but is broader than Poissonian at low masses ( or ). This behavior agrees with the mass functions computed from the level 2 Aquarius simulations, which probe an additional three decades in mass to . The scatter in the subhalo abundance among these simulations is fairly small: the fractional scatter, defined as the standard deviation in normalized by the mean of , is constant over the range at approximately 12%. This is somewhat smaller than our estimate from MS-II data but is consistent given the small number of Aquarius halos. Poisson scatter would give 4.5% at and 1.5% at . This is clearly inconsistent with the Aquarius results.

Many previous studies have analyzed subhalo occupation statistics using reduced moments of the subhalo mass function, defined via

(10)

These studies, which have usually been done in the context of HOD modeling, have used the values as measures of how close the scatter in subhalo count is to Poissonian: for a Poisson distribution, for all . All studies have found that over a wide range in mass and have concluded that the scatter is indeed Poissonian (e.g., Kravtsov et al. 2004a; Zheng et al. 2005). While the finding that does indeed justify the replacement of with when computing correlation functions in HOD modeling, it does not show that the variation about at fixed can be well-described by Poisson statistics: distributions that differ substantially from Poisson can result in mere percent-level variations from .

Figure 9: The reduced second moment parameter versus for several ranges of . The behavior of is nearly independent of host mass: , with , for . Despite this, a comparison with Fig. 8 shows that the satellite HOD is not Poisson at these mass fractions.

Figure 9 shows as a function of for the same ranges of host halo mass as Fig. 8. The behavior is similar at all masses: is very close to, but slightly greater than, unity for , corresponding to occupation numbers . The deviations from are quite small – 2 to 3% – over this range. This is precisely the range where the scatter in deviates systematically and substantially from the Poisson expectation: Fig. 8 shows that the scatter is 30% broader than Poissonian at and nearly twice as broad at , whereas over this same mass range. The satellite HOD is certainly not Poisson for even though is within 2% of unity. While Poisson scatter is an excellent approximation for , it is increasingly inaccurate at higher occupation number.

Why do we find a dispersion that is significantly broader than Poisson even when is only 2% larger than unity? Let us re-write as

(11)

From this expression, it is clear that any distribution where increases more slowly than will lead to for large values of . Equation (11) also shows that the deviation from can be directly related to the fractional scatter in the model of equation (9):

(12)

Even slight departures from are therefore important: for example, if , then , i.e., a 46% intrinsic scatter. These results suggest that more care is needed in treating the scatter in HOD modeling. It is still true that specifies the error incurred by assuming Poisson statistics when computing the moment of the HOD, however. Using Poisson statistics will therefore result in error in the second moment and error in the third moment of the HOD.

Figure 10: The subhalo occupation distribution from host halos in the MS-II with (upper two panels) and with (bottom panel). Black data points show the measured distribution, green lines show the Poisson distribution with the same , and magenta lines show the Negative Binomial distribution with the same and variance given by equation (9) with fractional intrinsic scatter . Both distributions match the data for (top panel, corresponding to ), while the Poisson distribution is noticeably too narrow for (middle panel, corresponding to ) and much too narrow for (bottom panel, corresponding to ). The Negative Binomial distribution matches the data well at all .

The subhalo occupation distribution can be well-approximated by the Negative Binomial distribution, which is given by

(13)

Here is the usual Gamma function and the parameters and are determined by the mean and variance of the distribution:

(14)

Using our model for the variance in , the values of and can be expressed in terms of and the fractional intrinsic scatter alone:

(15)

Given , it is therefore straightforward to compute the full distribution . Note that in the limit – i.e., as the intrinsic scatter goes to zero – the Negative Binomial distribution approaches the Poisson distribution with mean , as desired.

The good agreement between the measured subhalo occupation distribution and the Negative Binomial distribution is illustrated in Fig. 10. Black data points show the subhalo occupation distribution for (top panel, corresponding to ), 15 (middle panel, corresponding to ), and 125 (bottom panel, corresponding to ). The upper two panels use MS-II hosts with , while the lower panel uses MS-II hosts that are a factor of ten more massive. Poisson (green lines) provides a fairly good match for (the dispersion is 9% smaller than that of the data) but is noticeably too narrow at (the Poisson dispersion is 25% smaller than that of the data). For , Poisson is a very poor match to the data, with a dispersion that is too small by a factor of 2.25. The Negative Binomial distribution (magenta curves) matches the data well for all three values of .

While the results presented in this section were derived using the cumulative mass function, we have confirmed that they hold equally well for the cumulative maximum circular velocity function. Specifically, the intrinsic scatter in is nearly identical to that in , i.e., the subhalo abundance also shows an intrinsic scatter of 18% when expressed in terms of . This is not true when considering , however. This definition introduces additional scatter in the subhalo abundance due to the scatter in concentrations at fixed mass: host halos of a given have a range of values.

Ishiyama et al. (2009) have recently studied a sample of 125 well-resolved halos from a cosmological simulation; their work corroborates this result. They show that the scatter in is of order for host halos with . They find the scatter to be markedly reduced when using a criterion more similar to ours, for subhalos within , signifying that subhalo properties correlate more tightly to host than to host (see also Springel et al. 2008). This is yet another reminder that results on subhalo mass and circular velocity functions are sensitive both to how the host halo sample is defined and to how subhalos are selected (e.g., the limiting radius chosen).

Finally, we have investigated whether the subhalo mass function shows systematic variation with the host’s environment (as measured by the overdensity on a scale of ) or host halo properties. No obvious correlations exist. Furthermore, the variation in the subhalo mass function between hosts can be modeled as a simple normalization shift, i.e., the mass function within individual halos is consistent with the low-mass slope measured from the ensemble average. This is in good agreement with the behavior of the mass Aquarius halos’ mass functions (e.g., Fig. 5).

4.3 Massive subhalos

The parameter in equation (7) is approximately the mass fraction relative to of the most massive subhalo in a “typical” host; therefore shows that the most massive subhalo in a Milky Way-mass halo typically has about 1% the mass of its host. The distribution of for the most massive subhalo in each of our hosts is shown in Fig. 11. Data are plotted for both (solid histogram and filled squares) and for (dashed histogram and open squares666Note that what is plotted as the dashed histogram in Fig. 11 is the probability distribution for the subhalo at with the largest , not of the most massive subhalo at .). Values for the individual Aquarius halos are shown by colored symbols. Fig. 11 confirms that the probability distribution function for peaks at for subhalo masses measured at , albeit with a large spread: the probability for to be larger than 0.025 is 29% and to be larger than 0.1 is 4.3%, while 4.8% of halos have .

Figure 11: Distribution over of the most massive subhalo within of each of our Milky Way-mass hosts, both in terms of (solid histogram) and (dashed histogram). In each case, is computed relative to . The values of for each Aquarius halo in the MS-II are also noted with square symbols (filled for , open for ). The distribution of is fairly broad and peaks at (0.035) for (). The cyan lines show the prediction for (solid) and (dashed) if the most massive subhalo is Poisson-sampled from our analytic fit to the cumulative mass function. In both cases, the Poisson prediction agrees very well with the actual distribution.

The model presented in Sec. 4.2 suggests the scatter in the subhalo mass function should be nearly Poissonian for masses relevant for the most massive subhalo. We can test this directly: if the most massive subhalo in each halo is Poisson sampled from the underlying , then the probability density for finding in the range can be directly computed as

(16)

The cyan lines in Fig. 11 show computed this way for (solid line) and (dashed line); the agreement with the actual distributions for (the histograms) is excellent. The distribution of is therefore consistent with the Poisson hypothesis, both for the redshift zero masses and for the accretion masses.

The distribution of computed with respect to can be used to estimate the probability that the Milky Way should host a galaxy at least as massive as the Large Magellanic Cloud (LMC), the most luminous MW satellite. Using the abundance matching results of Guo et al. (2010) and a stellar mass of for the LMC (Kim et al., 1998), we find that . For a Milky Way mass of , there is a 8% chance of having a satellite with this accretion mass or greater. For a mass of – consistent with the abundance matching results of Guo et al. and the Local Group timing argument value obtained by Li & White (2008), this probability rises to 27%.

This same line of reasoning can be used to compute the probability of a dark matter halo hosting two satellite galaxies at least as massive as the Small Magellanic Cloud (SMC). The SMC has a total stellar mass of approximately (Stanimirović et al., 2004), corresponding to an abundance matching mass at accretion of . If the Milky Way has a dark matter halo of mass , then there is a 3.3% chance of having a second-ranked satellite with at least this infall mass. For a halo mass of , the probability becomes 20%.

These probabilities ignore any scatter in the relation. This scatter is subdominant to the uncertainty in the Milky Way’s mass, however. For example, Guo et al. (2010) give 80% confidence intervals that correspond to on the infall halo masses of the Magellanic Clouds based on a dispersion of 0.2 in at fixed . Clearly, the satellite statistics quoted above give some support for the higher estimates of the Milky Way’s halo mass.

4.4 Subhalo accretion times

Figure 12: Redshift of accretion as a function of for all resolved subhalos in galaxy-mass (left) and cluster-mass (right) host halos. The median relation is shown as a solid curve while the 25 and 75 (10 and 90) percentiles are shown as dashed (dotted) curves. More massive subhalos at the present day were accreted later on average than low-mass subhalos; the spread of accretion times for massive objects is also smaller. Subhalos of more massive host halos are accreted at lower redshifts on average, reflecting the hierarchical build-up of dark matter halos.

Subhalos persisting to are fossils of previous halo merger events. The accretion redshifts of subhalos are therefore important quantities for understanding the dynamical evolution of subhalos and can help constrain other quantities, such as merging timescales. The median for subhalos as a function of is shown as a solid curve in Fig. 12. The 25th and 75th (10th and 90th) percentiles of the distribution are included as dashed (dotted) curves. Results for Milky Way-mass host halos are shown in the left panel. The median is a decreasing function of , which is mainly a dynamical effect: massive subhalos that are accreted at early times are able to lose their orbital energy and angular momentum and merge with the dominant subhalo by while low-mass subhalos are not (e.g., Boylan-Kolchin et al. 2008). As a result, the only remaining massive subhalos are those that were accreted at relatively late times. At , the typical of a surviving subhalo is somewhat larger than , while at , for surviving subhalos. It is difficult to compare the individual Aquarius halos, as there are only 60 satellites across all six halos satisfying and , but the median accretion time for these 60 satellites, , falls well within the expected range based on the distribution from the MS-II.

For comparison, the right panel of Fig. 12 shows the distribution of for surviving subhalos in much more massive hosts, those with . The trend is very similar to that seen for Milky Way-mass hosts, namely, the median is a decreasing function of . The typical values are uniformly lower for the massive halo sample, however: at , the median surviving subhalo in the massive host sample was accreted at . This difference is a manifestation of the hierarchical growth of dark matter halos: more massive halos are dynamically younger and have accreted a larger fraction of their mass (and of their subhalo population) recently. Fig. 12 also indicates that a large fraction of the galaxy population in galaxy clusters has joined the cluster fairly recently, whereas most of the luminous satellites of galaxies similar to the Milky Way have been part of their host halos since . At all masses, the median value of is 2.4 (see Fig. 15), indicating that a typical subhalo has lost 60% of its mass.

The results in the left panel of Fig. 12 appear to differ substantially from those of Gao et al. (2004), who found that 90% of subhalos in a Milky Way-mass halo at have . These differences are due to the choice of definition for . Gao et al. define as the last time a subhalo entered the FOF group of the main progenitor of its host while we define to be the redshift when the subhalo’s mass was at its maximum. This latter definition is closer to the first time a subhalo enters its host’s FOF group. We have checked that using Gao et al.’s definition of gives results similar to theirs. Our results for more massive halos (the right panel of Fig. 12) agree with the findings of De Lucia et al. (2004, e.g., their fig. 10), who use the same definition of as we do.

5 Merger histories

Galaxy-mass dark matter halos are constantly bombarded by smaller halos in the CDM model. This has been the source of great concern with respect to the formation of MW-like galaxies: how is it possible to reconcile an active merging history with the existence of a thin stellar disk (e.g., Toth & Ostriker 1992; Velazquez & White 1999; Benson et al. 2004; Kazantzidis et al. 2008; Hopkins et al. 2009)? In order to understand the severity of the problem, we investigate the merger histories of MW-mass halos in this section, concentrating on mergers of objects with mass comparable to that of the MW disk.

When considering the merger history of a dark matter halo, two major questions are (1) what is the rate of infall of other, smaller dark matter halos across the halo’s virial radius and (2) what is the rate at which other dark matter halos lose their identity and merge with the center of the halo? Clearly, the set of objects that merge into the center of a halo will be a subset of those that crossed the halo’s virial radius at earlier times, but connecting the two is far from trivial. Several recent papers have addressed the frequency of halo-halo mergers. Fakhouri & Ma (2008), Guo & White (2008), and Genel et al. (2009) investigated the halo merger rate directly in the MS, while Neistein & Dekel (2008) and Cole et al. (2008) used extended Press-Schechter theory to derive merger rates and compared the results to the MS. The most relevant study for this paper is Stewart et al. (2008), who used an -body simulation to extract halo-halo merger rates for MW-mass halos. Stewart et al. found that while recent major mergers are very rare, mergers involving halos with masses of are quite common over the past 10 Gyr and mergers with halos having masses exceeding the mass of the galactic disk are virtually inevitable over the same time period. These results led the authors to conclude that the frequency of mergers may present a substantial challenge to our understanding of disk galaxy stability, as the majority of galaxies in MW-mass halos are disk-dominated (Weinmann et al., 2006; van den Bosch et al., 2007; Park et al., 2007)).

The mass resolution of the MS-II allows us to expand on the results of Stewart et al. by studying mergers with the central regions of MW-mass halos. These are more clearly relevant to the survival of stellar disks (which have characteristic radii that are less than 10% of ). Accordingly, in this section we consider only mergers between accreted subhalos and the dominant subhalos in MW-mass halos. We construct samples of merging objects by searching the main progenitor branches of the dominant subhalos for merger events with subhalos satisfying .

In order to ensure that the mergers we are considering are well resolved and are not affected by limited numerical resolution, we have computed the orbits and internal properties of subhalos immediately prior to merging. For the mergers considered here – those with for the satellite – the median separation between the satellite and the dominant subhalo immediately prior to merger is about ; in 80% (90%) of the cases, the separation is less than .

Furthermore, the satellites themselves are well-resolved immediately prior to merging: as we show in the Appendix, even halos at our lower limit of have on average 70 particles just prior to merging, while satellites with typically have 150 particles at the time of merging. Accreted halos that are capable of substantially impacting the Galactic disk can therefore be reliably tracked to within of , allowing us to separate true mergers with the centers of halos from accretion events that cross but remain at large halo-centric distances.

Figure 13: The fraction of central subhalos that have merged with an object having greater than 1, 3, 5, 7, 10, and 15% of (upper through lower curves) since redshift , as a function of . Virtually all halos have merged with a halo since , while only half have merged with a halo in the same redshift range. Recent merger events with are common: of halos have had such a merger since .

Figure 13 explores the probability for dominant MW-mass subhalos to have merged since redshift with another subhalo with (top to bottom). Approximately two-thirds of MW-mass halos have experienced a merger with a halo having since , and over 90% have had at least one such merger since . On the other hand, only 30% of halos have merged with a halo since and just 50% of halos have ever experienced such a merger. Approximately 50% of halos have had a merger with since z=1.

It is important to note that we have been quoting satellite masses in terms of the maximum dark matter mass the satellite has ever attained, . Subhalos usually lose a large fraction of their mass before merging: for our sample, the average subhalo loses 90% of its mass between accretion and merging (see Fig. 16), independent of . Subhalos with therefore typically have a mass ratio of 1:100 at merger, meaning the mass involved in a potential disk impact may be three times smaller than the disk mass rather than three times larger.

Figure 14: The cumulative probability of having a merger with for the 1st, 2nd, 3rd, 4th, and 5th-most massive central merger that a given halo experiences, both over the halo’s entire lifetime (solid lines) and since (dashed lines). The median value for the most massive merger is , while the 3rd-most massive merger typically has .

We can also investigate the distributions of merger histories in terms of the -most massive halo with which a MW-mass halo has merged. This is done in Fig. 14. The solid lines show the probability of having a merger with as a function of , while the dashed curves give the same probability for mergers at . Different colors correspond to different merger ranks: green shows the probability distribution for the most massive merger for each halo (in terms of ) while red shows the distribution for the 5th most massive merger. 90% of halos have had one merger event with since , 60% have had two such events, and 30% have had three. The median value for the most massive merger a halo has experienced since is , and the median redshift for such mergers777This is the redshift at which the satellite subhalo merges with the central subhalo, not the redshift at which the progenitor halos merge. is ; only 10% of halos have had two events with . The most massive merger event since for each of the Aquarius halos (A-F) has of 0.033, 0.078, 0.028, 0.030, 0.033, and 0.259, respectively. Halos A, C, D, and E therefore lie near the 10% most quiescent halos of Fig. 14 while halo B is fairly typical and halo F has had a recent major merger, which is rare.

A thorough understanding of the stability of galactic disks in the presence of satellite infall requires the inclusion of baryonic physics, as several studies have shown that factors such as the orbits and gas content of the merging galaxies can have a substantial impact on remnant properties (Springel & Hernquist, 2005; Kazantzidis et al., 2008, 2009; Scannapieco et al., 2009; Hopkins et al., 2009; Purcell et al., 2009; Stewart et al., 2009; Moster et al., 2010; Governato et al., 2009; D’Onghia et al., 2010) Nevertheless, we can still use the statistics computed above to learn about the frequency of satellite-disk interactions, as the dynamics that determine which subhalos merge are mostly set by dark matter halo properties such as the satellite’s orbit and mass at infall.

Approximately 90% of MW-mass halos have had a central merger with a halo of since . If such mergers destroy galactic disks, it is extremely difficult to understand the frequency of disk galaxies in MW-mass halos (; Weinmann et al. 2006; van den Bosch et al. 2007; Park et al. 2007). Around half of MW-mass halos have experienced a merger of since . These events are likely to strongly affect galaxy disks; however, the frequency with which they occur, coupled with the possibility that a non-negligible fraction of such mergers may have sufficient gas to reform a disk in the merger remnant, does not seem high enough to present a major obstacle to forming galactic disks in a majority of MW-mass halos. It is therefore essential to simulate a large number of mergers in the mass range of with a distribution of realistic orbits and galaxy models, including a variety of gas fractions. This will significantly advance our knowledge of disk formation and survivability in the CDM cosmology.

6 Discussion

6.1 How (a)typical are the Aquarius halos?

In general, the agreement between the properties of the Aquarius halos and those of the full sample of MW-mass halos from the MS-II is quite good: no halo is an outlier from all of the relations studied here. Halos A and C form somewhat earlier than average. Halo C also has a larger concentration than is typical. On the other hand, neither A nor C seems to be an outlier in the distribution of the spin parameter . Halo F is somewhat unusual in that it has a recent major merger; on the other hand, it is quite typical in most of its other properties. The Aquarius halos seem to reflect the diverse properties of Milky Way-mass halos, both in terms of assembly history and structure.

The Aquarius halos do not, however, seem to uniformly sample the spin parameter : five of the six halos lie below the median of the distribution, while the sixth is a merger remnant and may well relax to a lower value as well. We have explicitly checked that the isolation criterion in the Aquarius halo selection does not bias the spin parameters, as the distribution of from the isolated and non-isolated samples are statistically identical. The somewhat low spin parameter distribution of the Aquarius halos appears to be a statistical fluctuation. This may be connected to the relatively quiet recent merger history of the central subhalos for the Aquarius halos.

6.2 How (a)typical is the Galaxy’s halo?

While the Milky Way is frequently considered as the prototypical massive spiral galaxy, this has its roots as much in convenience as in evidence. While it is difficult to use dissipationless -body simulations to discern whether the MW is indeed a typical galaxy for a halo of its mass, we can investigate whether the Galaxy’s halo itself is typical among those of similar mass. Many possible constraints are predicated on a more precise determination of the mass of the MW’s halo, however. Current estimates vary by a factor of approximately 2 to 3 (with ). This uncertainty has a large effect on the probability of hosting massive satellite galaxies such as the SMC and LMC: the chances of having a subhalo capable of hosting the LMC or two subhalos capable of hosting the LMC and SMC are approximately 3-8% for a MW mass of but rise to 20-27% if the MW’s halo has a mass of .

On the other hand, a more massive Galactic halo implies that merger events of a given correspond to a larger fraction of the MW’s disk mass. If the MW’s thin disk must be explained via an unusually quiescent merger history over the past 10 Gyr (e.g., Hammer et al. 2007), we find that the MW’s halo must be even more unusual (in terms of merger history) if the halo is massive. If we assume the MW cannot have merged with a halo since , this corresponds to for and if . From Figure 13, we find a 50% chance for the former case but only a 20% chance for the latter. If halos with destroy thin disks, the probabilities drop to 12% and 4%, respectively.

A possible resolution of this tension – that the existence of the LMC and SMC argue for a more massive MW halo while the lack of recent mergers argues for a less massive Galactic halo – is to directly connect the two issues. A recent analysis of the LMC’s orbital history (Besla et al., 2007) based on the updated proper motions of the LMC (Kallivayalil et al., 2006) indicates that the LMC likely fell into the Milky Way’s halo within the last Gyr. The LMC’s infall mass is likely to be in excess of (and the SMC’s infall mass is likely to be only slightly less massive; Guo et al. 2010), meaning that the Milky Way has actually experienced a recent 1:10 merger at the halo level. Stewart et al. (2008) show that Milky Way-mass halos typically undergo one such merger at the halo level over the past 10 Gyr, perhaps indicating that the Milky Way’s accretion history is special only insofar as the LMC fell in so recently.

6.3 Cosmological Parameter Dependence

As noted in Sec. 2.1, the cosmological parameters used for the MS-II are slightly different from the most recent observational estimates: the values of and are approximately away from the best-fitting values of 0.81 and 0.96, respectively. Many of our results are insensitive to such changes. The abundance of MW-mass halos differs by only in the WMAP5 cosmology compared to the Millennium-II cosmology, for example, and the basic properties of halo assembly histories scale weakly with (van den Bosch, 2002). The subhalo abundance of MW-mass halos could potentially be affected by and , as decreasing tends to reduce halo and subhalo concentrations at a fixed mass (e.g., Macciò et al. 2008). This would tend to decrease the amplitude of in lower cosmologies, as subhalos would be more easily disrupted. On the other hand, halos tend to form later in such cosmologies, meaning that subhalos are subject to less dynamical evolution (see the discussion at the start of Sec. 4.1); this would tend to increase the amplitude of for a lower . These two effects compensate to some extent, reducing the overall shift.

In order to properly ascertain the effects of cosmology on the structural properties of MW-mass halos, it is necessary to run simulations that differ only in their cosmological parameters and to search for differences in statistical properties. Nevertheless, there is reason to think that changing from the cosmology of the Millennium, Millennium-II, and Aquarius simulations to the WMAP 5 cosmology will not have a strong impact. For example, the difference in amplitude of between and is only 14% (Reed et al. 2005; Diemand et al. 2008; Diemand & Moore 2009), so reducing from 0.9 to 0.81 will likely affect the amplitude of the subhalo mass function by a few percent only.

7 Conclusions

We have presented a statistical study of Milky Way-mass halos from the Millennium-II Simulation, which contains over 7000 halos in the mass range . Our principal results can be summarized as follows:

  • As several previous studies have shown, halo growth proceeds in an “inside-out” fashion: the central gravitational potential of a typical MW-mass halo reached half of its present-day value by , on average, whereas the virial mass reached half of its present-day value at .

  • The ratio of mass to the host halo mass for the most massive subhalo of a MW-mass halo has a broad distribution that peaks at 1%. The corresponding quantity computed using infall mass instead of peaks at 3.5%

  • The differential and cumulative subhalo mass functions, computed in terms of , are both well-fitted by a power law with an exponential cut-off at large . We find both the abundance per log decade in and the number of subhalos with mass greater than to be proportional to for small .

  • The scatter in the cumulative mass function at fixed is only well-approximated by a Poisson distribution at large (corresponding to mean occupation number ). Intrinsic scatter of approximately 18%, independent of the host halo mass, becomes increasingly important at lower (larger ) and is dominant over Poisson scatter at ().

  • The Negative Binomial distribution with variance in that is equal to the sum of a Poisson term and a fractional intrinsic scatter of 18% matches the subhalo occupation data well at all for host halo masses between and .

  • The statistic [equation (10)], which is frequently used to characterize deviations of the HOD from Poisson, is not discriminating: distributions that differ strongly from Poisson can lead to few percent differences of from unity.

  • Accretion redshifts of subhalos vary systematically with host halo mass and with the ratio of subhalo to host masses. Massive subhalos at were typically accreted more recently. This is a dynamical effect: massive subhalos accreted at early times either merge via dynamical friction or are heavily tidally truncated. At fixed , subhalos of more massive host halos were accreted more recently, reflecting the hierarchical build-up of dark matter halos. Subhalos with have a median accretion redshift of for hosts with and for halos with .

  • The frequency of central mergers – that is, the merger of a satellite with the central region of a MW-mass halo – is a strong function of the satellite’s . Since 90% of MW-mass halos have experienced a central merger with a satellite having while 50% have had such a merger with a satellite, and only 30% with a satellite. The compatibility of thin stellar disks with the frequent mergers expected in CDM thus depends strongly on exactly which mergers destroy disks. If