WHERE ARE THE FOSSILS OF THE FIRST GALAXIES? I. LOCAL VOLUME MAPS AND PROPERTIES OF THE UNDETECTED DWARFS

Mia S. Bovill; Massimo Ricotti

doi:10.1088/0004-637X/741/1/17

1. INTRODUCTION

Simulations of the formation of the first galaxies have matured; however, there are few observational constraints on the models. The overabundance of dark matter satellites near the Milky Way when compared to observations of luminous satellites (missing galactic satellites) and lack of dwarfs in the voids (void phenomenon) suggest that these early galaxies may be too faint to be detected directly at high redshift even for the James Webb Space Telescope (Ricotti et al. 2008; Johnson et al. 2009). However, we can detect their fossil remnants in the local universe. Hierarchical formation scenarios predict these first galaxies formed before reionization in dark matter minihalos with masses ≲ 10⁸ M_☉. Those that survive to the present constitute, in part, a sub-population of satellites around larger halos. Recent observational and theoretical advances allow us to compare simulated primordial galaxies to the observations of the faintest galaxies in the Local Group (Belokurov et al. 2006, 2007; Irwin et al. 2007; Walsh et al. 2007; Willman et al. 2005a, 2005b; Zucker et al. 2006a, 2006b; Geha et al. 2009) and constrain models of star formation in the early universe.

The formation of the first galaxies before reionization is regulated by complex feedback effects acting on cosmological distance scales. These self-regulation mechanisms have dramatic effects on the number and luminosity of the first galaxies. The gas in minihalos with a circular velocity v_max < 20 km s⁻¹ is heated to T ≲ 10,000 K during virialization. At this temperature, a gas of primordial composition is unable to cool and initiate star formation unless it can form and retain a sufficient density of H₂. Although H₂ is destroyed by dissociating UV radiation in the Lyman–Warner bands, its formation can be catalyzed by hydrogen ionizing radiation, via the formation of H⁻ (Haiman et al. 1996). Thus, radiative transfer is necessary to simulate H₂ formation and destruction in the optically thick early universe. Along with other relevant physics, three-dimensional radiative transfer was included in Ricotti et al. (2002a, 2002b, 2008), hereafter referred to as pre-reionization simulations.

Ricotti et al.'s pre-reionization simulations produce a population of primordial galaxies with stellar properties consistent with a subset of the classical dSphs (Ricotti & Gnedin 2005, hereafter RG05) and a majority of the recently discovered ultra-faint dwarfs (Bovill & Ricotti 2009; Ricotti 2010, hereafter BR09). Over the last five years, the ultra-faints doubled the census of Milky Way and M31 satellites (Willman 2010 and references therein). The new Milky Way dwarfs have V-band luminosities <10⁵ L_☉, half-light radii r_hl ∼ 20–300 kpc, metallicities [Fe/H] ≲ −2, mass-to-light ratios ≳ 100 M_☉/L_☉, stellar velocity dispersions ∼2–10 km s⁻¹, and are 30–400 kpc from the Milky Way. The new M31 dwarfs have L_V ≳ 10⁵ L_☉ with r_hl systematically greater than their Milky Way counterparts. Except for Leo T, both groups are dominated by an old metal-poor population and are devoid of gas. In addition, a subset of the M31 dwarfs are falling into their host halo for the first time (Majewski et al. 2007). These new dwarfs provide the observational laboratory we need to test our model of the properties and distribution of the fossils of the first galaxies.

The large scatter in the properties of Ricotti et al. galaxies with the same mass follows from the physics governing star formation in minihalos. Figure 1 shows that in minihalos the mean star formation efficiency, f_*, is much smaller than the star formation efficiency per free-fall time ( epsilon _* = 10%) that appears as a free parameter in the adopted sub-grid recipe for star formation: ${\dot{\rho }_*}=\epsilon _* \rho _g/t_{{\rm dyn}}$ . In contrast, in halos with mass M > 10⁸–10⁹ M_☉ we find f_* ≃ epsilon _*. The value of f_* is smaller than epsilon _* because in small halos the gas available for star formation is reduced with respect to the mean cosmic value. Two effects are dominant in reducing the amount of gas available for star formation in small mass galactic halos: (1) photo-heating of the intergalactic medium (IGM) reduces the amount of gas that falls into the potential wells with respect to the mean cosmic value; and (2) massive stars inside luminous galaxies ionize and heat their interstellar medium (ISM) expelling most of the gas available for further star formation (see Ricotti et al. 2008). Supernova explosions also contribute to produce galactic winds but UV radiation from massive stars operates on a shorter timescale and is most effective in clearing out the majority of the gas in the halo. Thus, the number of massive stars that can be produced in such small mass galaxies is regulated by feedback loops. As a result, 〈f_*〉 is much smaller than epsilon _* and its value is rather insensitive to the assumed value of epsilon _*.

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** Mean star formation efficiency 〈f_*〉 as a function of v_max for all galaxies from RG05 at redshifts z = 14.6, 12.5, 10.2, and 8.3. Here, the circular velocity of the galaxies is v_max = 17km s⁻¹(M/10⁸ M_☉)^1/3[(1 + z)/10]^1/2, calculated from their total mass M and redshift z. The star formation efficiency is the fraction of baryons within the virial radius of the halo that is converted into stars, f_* ≡ (Ω_m/Ω_b)M_*/M.
Download figure:
Standard image High-resolution image

Although, as shown in Figure 1, the mean star formation efficiency in the first galaxies is roughly a power law 〈f_*〉∝M²∝v⁶_cir, the scatter about this relationship is large, especially for the smallest mass halos in the plot (Ricotti et al. 2008). Halos with masses M < 10⁷ M_☉ can be either luminous or dark, depending on the environment in which they reside. For a given mass, halos that are relatively isolated and far from luminous galaxies tend to be the least luminous or completely dark. The proximity to a source of ionizing radiation and metals has a positive effect on triggering star formation in dark halos. Ionizing radiation stimulates molecular hydrogen formation and enhances the cooling rate of the gas (Ricotti et al. 2001, 2002a). Similarly, metal pollution stimulates gas cooling. Hence, the first galaxies are highly biased. The most luminous of them are most likely to merge into more massive objects and do not survive to the present. This explains why in the simulations shown in the present study, contrary to our naive initial expectation, we did not find fossil dwarfs with L_V > 10⁶ L_☉ around the Milky Way.

Gnedin & Kravtsov (2006, hereafter GK06) use this approximation in conjunction with high-resolution N-body simulations of the Local Group to evolve a population of dwarf galaxies around a Milky Way mass halo from z = 70 to z = 0. For details of the simulations, see Section 2 in GK06. GK06 defines a fossil as a simulated halo which survives to z = 0 and remains below the critical circular velocity of v_filt = 30 km s⁻¹ with no appreciable tidal stripping. They calculate the probability, P_S(v_max, r), of a luminous halo with a given maximum circular velocity v_max to survive from z = 8 (the final redshift of the RG05 simulation) to z = 0. For a given v_max, the number of surviving dwarfs at z = 0 is N(v_max, z = 8)P_S(v_max, r), where P_S is the survival probability for a satellite at a distance r from the host halo. The surviving halos are assigned a luminosity based on the L_V versus v_max relationship from RG05. At z = 0, GK06 has a population of dwarf galaxies with a resolution limit of v_max = 13 km s⁻¹. The halos are statistically assigned luminosities from the L_V–v_max relation given in Figure 3 of GK06. Unfortunately, with this relation, the 13 km s⁻¹ corresponds to a lower luminosity limit of L_V ∼ 10⁵ L_☉, which includes Leo T and Canes Venatici I, but excludes all the other new ultra-faint Milky Way satellites.

In this paper, we describe and test a novel method of generating N-body initial conditions which allows us to follow the evolution, merger rates, and tidal destruction of pre-reionization halos to present day and to overcome some of the limitations of the GK06 method. The initial distribution of particles in the N-body simulations represents the position, velocity, and mass distribution of the dark and luminous halos extracted from pre-reionization simulations. Our simulations have a sufficiently large volume and dynamical range to explore the distribution of fossil galaxies outside the Local Group, and in nearby filaments and voids using limited computational resources.

Our method improves on the GK06 work by removing the constraints that preclude a comparison of the GK06 simulations with the observed distributions of the ultra-faints. (1) Due to the resolution of their N-body simulations, GK06 cannot resolve dwarfs with circular velocity, v_max, <13 km s⁻¹, which roughly corresponds to a simulated dwarf with L_V < 10⁵ L_☉. With two exceptions, no ultra-faint dwarfs have L_V > 10⁵ L_☉ (CVn I, Zucker et al. 2006b and Leo T, de Jong et al. 2008). (2) The statistical matching of the baryonic properties of the pre-reionization halos to equivalent z = 0 halos in their N-body simulation does not allow GK06 to account for mergers of pre-reionization halos after reionization. While the majority of mergers would not involve two luminous pre-reionization halos, the effect cannot be ruled out a priori. (3) The GK06 statistical matching also does not account for the clustering bias of the most luminous pre-reionization halos. The formation efficiency of H₂ is dependent on stochastic effects, so the most luminous pre-reionization halos form in the highest density regions of the Ricotti et al. (2002b) simulations and are more likely to have undergone a merger with another massive, luminous pre-reionization halo. (4) Extracting the baryonic properties of their fossils at z = 0 from the final output of the pre-reionization simulation does not allow GK06 to account for cosmic variance. By z = 0, the faster evolution of structure in overdense regions (i.e., Local Group) and the slower structural evolution of underdense regions (i.e., Local Void) have produced significant variance in the numbers and types of objects seen in both.

The paper is organized as follows. In Section 2, we describe our initial conditions in detail before comparing our results to traditional CDM N-body initial conditions in Section 3. We show that with this method we can easily achieve the resolution necessary to study the distribution of ultra-faint dwarfs at z = 0 in a volume similar to the Local Volume. Initial results are in Section 4 and the remainder are presented in Bovill & Ricotti (2011), Paper II of this series. Our simulations are compared to the GK06 results in Section 4.2 in which we propose a maximum luminosity threshold for primordial fossils. Finally, we compare the properties of our simulated primordial dwarfs to observations of Local Group dwarfs in Section 4.3 and present observational tests for the primordial formation model in Section 5. A summary and conclusions are presented in Section 6.

2. NUMERICAL METHOD

To achieve the resolution necessary to study the ultra-faint dwarfs in a z = 0 volume equivalent to the Local Volume, we developed a method for generating initial conditions for N-body simulations which provides the required mass resolution, while using only limited computational resources. Our simulations allow us to trace the merger rate and tidal stripping of the first galaxies from reionization to the modern epoch. Traditional initial conditions for CDM simulations begin with an evenly distributed grid of uniform mass particles before the positions and velocities are perturbed according to a given power spectrum. Our method follows the same concept, except the initial distribution of the particles is not a uniform grid, but represents the distribution of halos in the final outputs of a 1 Mpc³ high-resolution cosmological hydrodynamical simulation run to z = 8.3 (Ricotti et al. 2002b; hereafter we refer to these as the pre-reionization simulations and the halos found in their 1 Mpc³ outputs as pre-reionization halos). Thus, each particle represents a dark or luminous halo with a different mass and given stellar properties.

All of the initial conditions described in this section were runs from their initial redshift, z_init, to z = 0 using Gadget 2 (Springel 2005) on the Maryland HPCC Deepthought and analyzed with the Amiga halo finder (AHF; Knollmann & Knebe 2009; Gill et al. 2004).

To construct our high-resolution region, we produce a lattice of the pre-reionization simulation z = 8.3 output. This gives us a grid similar to that used in traditional CDM, except power on scales below 1 Mpc is already present through the positions of the pre-reionzation halos. To add the larger scale power, we perturb the particle positions and velocities of the pre-reionization halos according to a power spectrum with no power for modes l < 1 Mpc. This method is similar to the one described in Tormen & Bertschinger (1996; see Appendix A for details). The high-resolution region ∼10 Mpc on a side, with a mass resolution of ∼3.2 × 10⁵ M_☉, is embedded in a coarse resolution volume 50 Mpc on a side containing 250³ particles at z = 8.3.

When compared to traditional zoom simulations, our high-resolution region has several key differences. Primarily, each of our particles represents a resolved halo from the pre-reionziation simulations. Each of these pre-reionization halos has a set of dark matter and stellar properties derived at z = 8.3. This technique allows us to push our simulations to higher mass resolutions over a "Local Volume" sized region without a prohibitive increase in the number of particles. However, this technique precludes us from determining detailed density profiles of pre-reionization halos at z = 0. The stellar properties of the pre-reionzation halos are preserved through the unique IDs of each particle in our simulation. If, in the modern epoch, a given pre-reionization halo is in a halo whose maximum circular velocity has never exceeded the filtering velocity, it has not accreted gas from the IGM after reionization. This filtering velocity, v_filter represents the critical value for the maximum circular velocity of the halo below which star formation is suppressed by the reheating of the IGM via reionization feedback (Gnedin 2000; Babul & Rees 1992; Efstathiou 1992; Shapiro et al. 1994, 2004; Thoul & Weinberg 1996; Quinn et al. 1996; Navarro & Steinmetz 1997; Susa & Umemura 2004; Benson et al. 2006; Hoeft et al. 2006). The subsequent lack of star formation in these low-mass halos allows us to approximate its present-day observable properties from those at reionization. We consider the initial conditions built using the method described above as our first-order simulations, specifically, runs A, B, and C (Table 1). The initial conditions for run D, which are significantly different than those described for runs A–C, are described in the next section. A detailed methodology for both sets of initial conditions is given in Appendices A and B.

Table 1. Simulation Runs

Name	IC Method	Volume	HR Volume	Mass Resolution		z_init
		(Mpc³)	(Mpc³)	(10⁶ M_☉)	(kpc)
A	First order	50³	∼9³	3.16	1	8.3
B	First order	50³	∼9³	1.0	1	8.3
C	First order	50³	∼9³	0.316	1	8.3
D	Second order	50³	∼9³	0.316	1	10.2

Note. From left to right, the columns are (1) run identifier, (2) type of initial conditions, and approximate reionization model, (3) size of low-resolution volume in Mpc³, (4) approximate size of the high-resolution volume in Mpc³, (5) mass of dark tracer particles in 10⁶ M_☉, (6) softening length in kpc, and (7) initial redshift of the zoom simulation.

Download table as: ASCII Typeset image

For the remainder of this work, we focus on run C since runs A and B do not have the resolution necessary to study the dwarf populations inside the Milky Way halo.

2.1. Approximating Cosmic Variance

For our first-order simulations (see Table 1), we assume that every part of our "Local Volume" evolves at the rate associated with the mean density of the universe, before and after reionization. However, there are deviations from this mean due to linear perturbations on large (>1 Mpc) scales. The evolution of a given region depends on its mean density with regions of higher density evolving faster than their lower density counterparts (Cole 1997; Reed et al. 2007; Crain et al. 2009). As a result, halos will collapse, and form stars, at later times in the voids compared to the filaments. To account for this effect, we relate the overdensity or underdensity of each region of our high-resolution region to the speed of its evolution. We express the evolution of a region as a function of its densities as z_eff = z_init + Δz, where z_eff is the effective redshift, z_init is the redshift of the simulations and the effective redshift of a region whose local density is the average density of the universe, ρ_o(z_o), and

$\begin{equation} \Delta z = (1+z_{{\rm init}})[(1+\delta)^{-0.6}-1)] \end{equation} \tag{ 1 }$

is the correction to z_init due to δ, the local overdensity or underdensity of a given region.

To approximate this variance effect, we produce a set of second-order initial conditions as in Cole (1997) (run D). The primary difference between runs D and C lies in the construction of the high-resolution volume. Instead of using a single pre-reionization simulation output at z = 8.3 (runs A–C), we use outputs at multiple redshifts (z = 8.3–14) to approximate the differential evolution of the universe up to z_init = 10.2. Before constructing our high-resolution region, we calculate the effective redshift of each 1 Mpc³ subvolume. Each subvolume is then assigned a pre-reionization output based on its effective redshift, with the lowest density voids at z_eff = 14 and highest density regions at z_eff = 8.3. For additional details, see Appendix B.

In addition to accounting for cosmic variance, comparisons of runs A–C and run D allow us to probe two different reionization scenarios. Since we cannot account for baryonic evolution after "reionization" when the pre-reionization outputs are transformed into our N-body simulation, we assume that the photoevaportation and reheating during reionization precludes any further baryonic evolution in the minihalos. During reionization by UV, we assume that the IGM throughout our volume was reheated to ∼10⁴ K (Ricotti & Ostriker 2004). We also assume that the entire volume was reionized at ∼z_init. For runs A–C, this approximates reheating at z_init ∼ 8.3 by UV photons generated by stars in the first galaxies (Sokasian et al. 2004; Wise & Cen 2009). Since the voids evolve at a slower rate than the filaments, using the same pre-reionization output for our entire simulation is effectively allowing the low density regions to evolve for a longer time before their IGM is reheated to 10⁴ K. This is consistent with reionization and reheating beginning in the filaments before spreading into the voids. Since low-mass halos in the voids would have had more time to accrete gas and form stars before the reheating cut off their gas supply, we expect the voids in runs A–C to be significantly brighter than those in run D (Figure 2)

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** Large-scale structure of the same region of our simulations for the low-resolution simulations (top), run C (middle), and run D (bottom). For runs C and D, white shows the halos with no luminous component while the larger, red points show the luminous pre-reionization halos. The color of the latter does not depend on luminosity. The bar across the top of each panel shows the scale in Mpc.
Download figure:
Standard image High-resolution image

Since each 1 Mpc³ subvolume in run D used a pre-reionization output consistent with its effective redshift, the entire high-resolution region has been given the same amount of time to evolve. When we transition from the pre-reionization output to our N-body simulations, run D does not allow low-mass halos in the low density regions to continue to evolve as the denser filaments are reheated. Instead, run D approximates a universe in which all of space is reionized and reheated at approximately the same time. Uniform reheating of the filaments and voids is a characteristic of reionization and reheating from X-rays produced by remnants of the first stars accreting from the ISM at higher redshift (Venkatesan et al. 2001; Ricotti & Ostriker 2004; Ricotti et al. 2005; Ripamonti et al. 2008; Shull & Venkatesan 2008). X-rays could also be produced by primordial black hole binaries (Mirabel et al. 2011; Saigo et al. 2004). When compared to reheating by UV radiation, X-ray reheating produces noticeably darker voids.

2.2. A Note on the Halo Occupation Distribution

The luminosities of the z = 0 halos are determined by the luminosities of their component pre-reionization halos. These pre-reionization luminosities are taken directly from the Ricotti et al. (2002a, 2002b) pre-reionization simulations and are determined by the feedback prescriptions and star formation efficiencies used in that work. Predictions made based on the resulting luminosity function and galactocentric radial distribution are a result of the primordial formation model we assume for the smallest dSphs. Note that the match of luminosity and dark matter mass in this simulation is not statistical, but rather a direct result of the distribution of the remnants of the first galaxies in the modern epoch.

3. TESTS OF THE METHOD

In this section, we present consistency checks of our method to confirm that it reproduces known results from previous CDM simulations. We also discuss numerical effects introduced by our use of a spectrum of masses in our high-resolution region. First, we confirm that the large-scale structure and clustering of matter are consistent with traditional CDM simulations run with constant mass per particle. Then, we see that the halo mass function is consistent with the mass function of halos derived from the Press–Schechter formalism (Press & Schechter 1974). Finally, we confirm that the number of subhalos and their Galactocentric distribution agree with the published results of the Via Lactea (Diemand et al. 2008) and Aquarius (Springel et al. 2008) simulations and that mass loss due to tidal stripping of the z = 0 halos is also in agreement with Kravtsov et al. (2004).

Figure 2 shows a region of our Local Volume 5 Mpc across at z = 0. From top to bottom are the coarse simulation; runs C and D are seen from the same viewing angle. For runs C and D, the luminous pre-reionization halos are shown as large red dots plotted over the white distribution of dark tracer particles. We find both runs C and D reproduce the large-scale structure seen in the coarse resolution simulation.

3.1. Mass Resolution

Our simulations produce maps of the present-day distribution and properties of pre-reionization fossils in a 5³ Mpc³ volume around a Milky Way type halo and in local filaments and voids. One of our goals is to map the distribution and properties of fossil galaxies outside the large hosts. This is done to quantify the number and properties of luminous dwarfs in the voids if stars are formed in minihalos before reionization. These dwarfs would have evolved in relative isolation, and, if found by observations, would represent unambiguous and unperturbed fossils of the first galaxies. However, at this time, the only observational sample to which we can compare our simulations is the classical dSphs and ultra-faint dwarfs near the Milky Way and M31. The faintest known dwarfs (L_V < 10³ L_☉) are found at less than 50 kpc from the Galactic center. Observations are likely incomplete at R > 50 kpc (approximately one-quarter the virial radius; Koposov et al. 2008; Walsh et al. 2009). To compare our simulations to observations of the faintest known dwarfs, we must resolve halos within 100 kpc of the Milky Way center.

Run A, with a minimum particle mass of 3.5 × 10⁶ M_☉, was not able to resolve subhalos within 200 kpc of the Milky Way mass halos. In run C, we increase our mass resolution to 3.5 × 10⁵ M_☉ by increasing the number of pre-reionization halos in the initial conditions. By decreasing our minimum pre-reionization halo mass to 3.5 × 10⁵ M_☉, we are able to resolve subhalos at R > 50 kpc (see Figure 3). At z = 0, a luminous pre-reionization halo, evolving in isolation, is surrounded by a cloud of lower mass pre-reionization halos and tracer particles. The number of dark particles increases with the total mass of the luminous pre-reionization halo and the mass resolution of the simulation. The detectability of the lowest mass halos by AHF is dependent on the ability of the luminous pre-reionization halos to accrete and retain their clouds of tracer particles. The larger number of low-mass pre-reionization halos and tracer particles in runs C and D will allow more pre-reionization halos to accrete large enough clouds to be detected as a present-day halo.

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** We show images of our three Milky Ways. From top to bottom, MW.1 (1.82 × 10¹² M_☉) from run C, MW.2 (0.87 × 10¹² M_☉), and MW.3 (1.32 × 10¹² M_☉) from run D. The left panels show both the dark (white) and luminous (red) pre-reionization halos. The right panels show only the luminous pre-reionization halos in gray scale with the brightest pre-reionization halos in white. In the right panels our Milky Ways have been rotated ∼180° relative to the view in the left-hand panels. The bar across the top of each panel shows the scale in Mpc
Download figure:
Standard image High-resolution image

Resolving subhalos near a large galaxy is complicated by the background density field of the host halo and the stripping of the clouds of tracer particles during tidal interactions. To resolve a subhalo in the inner 100 kpc of Milky Way mass halo, the pre-reionization halo must retain enough of its cloud to be considered a bound system. In addition, it must have a high enough central density to be seen against the background of the host halo. The effect of the larger mass of the pre-reionization halo on the central concentration of the subhalo will be discussed in Section 3.2. For AHF, the lower limit to robustly detect halos at z = 0 is a cloud of ∼50 tracer particles (Knollmann & Knebe 2009).

Our simulations cannot provide information on the z = 0 stellar properties of a halo for a pre-reionization halo which has undergone significant tidal disruption. Beyond the stripping of the accumulated dark cloud described above, our simulations do not allow for the breaking apart of the pre-reionization halos. We only consider a pre-reionization halo unaffected by tides if its cloud of dark particles remains intact and detectable. This negates comparisons within 50 kpc of a host halo where a significant number of the observed ultra-faint dwarfs have been modified by tides and our luminous pre-reionization halos are stripped of their clouds of tracer particles. In the galactocentric radial distributions and luminosity functions presented in Paper II of this series, we only include the simulated and observed sample at R > 50 kpc.

3.2. Softening Length

In this section, we discuss one of the most prevalent numerical effects of using a spectrum of particle masses in our high-resolution region instead of uniform particle masses. The effects of this spectrum of masses primarily manifest in the lower mass halos and are sensitive to our choice of the softening length, epsilon . We also show the mass functions from runs C and D for our chosen softening length.

Typically, the softening length is set at 2% of the average distance between particles in comoving coordinates. For a representative volume of the universe with particles of uniform mass, epsilon = 0.02N^−1/3 Mpc, where N is the number of particles per Mpc³. For the high-resolution region, we have a particle mass range between 3.5 × 10⁵ M_☉ and 2.5 × 10⁸ M_☉, requiring particle softening lengths from 0.5 kpc to 2 kpc. The public version of Gadget 2 does not have the capability of assigning softening lengths to each particle. Therefore, we must choose a single softening length for all the particles in the high-resolution region. To determine the optimal value of epsilon , we have run the same initial conditions with softening lengths in our high-resolution region of epsilon = 0.1 kpc, 1 kpc, and 5 kpc. We find that the best results for epsilon = 1 kpc (corresponding to a uniform particle mass of ∼10⁷ M_☉).

The right panel of Figure 4 shows the mass functions of runs C and D compared to the Press–Schechter mass function run with epsilon = 1 kpc. For C and D, we see a deficit in the number of 10⁹–10¹¹ M_☉ halos when compared to the Press–Schechter and an over abundance of M < 10⁷ M_☉ halos. The deficit for larger halos may result from the location of our high-resolution region. The Press–Schechter is the mass function of a typical volume of the universe. Our high-resolution region is underdense, containing three filaments bordering a void. The overabundance for M < 10⁷ M_☉ halos has a slope similar to the initial halo mass function from the pre-reionization simulations. At those masses, the z = 0 halos are dominated by one pre-reionization halo. This suggests that the steeper slope of the mass function at low masses is a numerical effect reflecting the behavior of the z = 8.3 mass function from the pre-reionization simulations.

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** Left: mass functions for run C evolved with three different softening lengths, 100 pc (red), 1 kpc (black), and 5 kpc (blue). Note that while there is a negligible difference at large masses, = 1 kpc gives us the least residual when compared to the expected mass function. The dotted line is the Press–Schechter for a ∼7 Mpc³ volume, equivalent to the mass of the bound halos. Right: mass function of all halos found by AHF in our high-resolution region for run C (gray line) and run D (black line). The dotted line is the Press–Schechter for a ∼7 Mpc³ volume at z = 0, equivalent to the mass of the bound halos. In both mass functions, we only include z = 0 halos which contain only high-resolution particles.
Download figure:
Standard image High-resolution image

When epsilon is set lower than 1 kpc (red curve in Figure 4), low-mass halos with one or more luminous pre-reionization halo are preferentially destroyed by numerical effects. Statistically, luminous pre-reionization halos are more massive than their dark counterparts, hence they migrate to the centers of their modern halos via dynamical friction. Any two-body interaction between a luminous pre-reionization halo and lower mass dark tracer particle will result in artificial heating. Over the entire simulation, such interactions artificially heat the cloud of tracer particles until it disperses. We find that for epsilon = 0.1 kpc only the most massive pre-reionization halos with the deepest potentials are able to retain their clouds. Isolated pre-reionization halos are surrounded by an extremely tenuous cloud of low-mass dark particles, which is not detected by AHF as a bound halo.

Using epsilon > 1 kpc also artificially decreases the number of the low-mass halos (blue curve in Figure 4). Unlike the deep potentials of the massive halos, the potentials of halos with masses M < 10⁸ M_☉ are relatively shallow. If epsilon is too large, the low-mass potentials will be flattened to the point where the pre-reionization halos are unable to accrete the tracer particles required for AHF detection. In halos with M ≳ 10⁹ M_☉, this effect is minimal. However, we are primarily interested in halos with M < 10⁹ M_☉.

3.3. Subhalo Scale Comparisons

In this section, we study the distribution of subhalos around our Milky Ways. We use runs C and D to explore the simulated distribution of z = 0 subhalos around our Milky Way mass hosts. Comparisons are made with other CDM simulations and with observations.

In each simulation, we search for Milky Way type halos, using observational and theoretical constraints. This gives us a range of halo masses for candidate Milky Ways of ∼0.6 × 10¹² M_☉–4 × 10¹² (Watkins et al. 2010; Kallivayalil et al. 2009; Zaritsky et al. 1989; Klypin et al. 2002) and upper mass estimates for the Local Group of ∼5.3 × 10¹² (Li & White 2008; van der Marel & Guhathakurta 2008). These criteria give us three Milky Ways, one in the run C and two in the run D, respectively (Table 2). All three hosts have masses on the low end of the observed Milky Way mass range.

Table 2. The Milky Way Mass Halos

Name	Run	Mass	R_vir	v_max
		(10¹² M_☉)	(kpc)	(km s⁻¹)
MW.1	C	1.82	248.1	203.4
MW.2	D	0.87	222.6	196.6
MW.3	D	1.32	194.7	177

Note. The columns from left to right are (1) the Milky Way identifier, (2) the simulation run each Milky Way is embedded in, (3) the AHF-derived virial mass of each halo in 10¹² M_☉, and (4) the AHF-derived v_max of the halo in km s⁻¹.

Download table as: ASCII Typeset image

The Milky Way halo in run C, MW.1, is in one of the highest density regions of our volume, with a companion galaxy of mass 10¹¹ M_☉ at a distance of 2 Mpc. The Milky Ways in run D have masses of 0.87 × 10¹² M_☉ for MW.2 and 1.32 × 10¹² M_☉ for MW.3. Though they are both in filaments, the nearby environments of MW.2 and MW.3 differ (see Figure 2). MW.2 sits at the intersection of three filaments, and there are ∼10¹¹ M_☉ halos within 1.5 Mpc. In contrast, MW.3 is only 1–2 Mpc away from a complex of galaxies with masses ∼10¹¹ M_☉ that appears to be in the process of merging to form another Milky Way mass system. For our comparisons with traditional simulations, and with observations, we use all three Milky Way mass halos. This allows us to explore differences between the first and second order as well as variations introduced by environmental effects.

Before looking at the distribution of satellites around individual Milky Ways, we check the distribution of the number of dark matter subhalos as a function of host mass. In Figure 5, we show a linear relation between host mass and the number of satellites for both runs C and D. There is good agreement with the Via Lactea and Aquarius runs when we adjust their results for our lower mass resolution. Our simulations can robustly resolve halos with M > 10⁷ M_☉ (v_max ≳ 5.5 km s⁻¹). To scale the number of subhalos within R_vir in the Via Lactea and Aquarius simulations, we use v_max ∼ 5 km s⁻¹ for Via Lactea and v_max ∼ 7 km s⁻¹ for Aquarius (from Figure 27 in Springel et al. 2008).

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** Number of satellites as a function of the mass of the host halo. The total number of satellites within R_vir for each halo is represented by circles, and the number of luminous satellites within R_vir for each halo by triangles. The results from runs C and D are shown as the open and filled symbols, respectively. The predictions from Via Lactea II (Diemand et al. 2008) and Aquarius (Springel et al. 2008), scaled to our mass resolution, are shown as the open stars. Aquarius is the star with a greater number of satellites within R_vir. The ranges of the Tollerud et al. (2008) and Walsh et al. (2009) predictions at 200 kpc are the green and purple barred lines, respectively.
Download figure:
Standard image High-resolution image

We use knowledge of the stellar properties of the pre-reionization halos to investigate the expected number of luminous satellites for a given host mass. We consider a subhalo luminous if it contains at least one pre-reionization halo with M_* > 10² M_☉, or has a z = 0 mass M > 10⁹ M_☉ (v_max > 20 km s⁻¹). To study the distribution of the number of luminous satellites, N_sat(L_V > 10² L_☉) versus M_host, we do not need to know the luminosity of the satellite at z = 0, only whether it is luminous. We find that all of the luminous subhalos in Figure 5 formed stars before reionization since we have no z = 0 halos above the threshold for post-reionization gas accretion (10⁹ M_☉: v_max = 20 km s⁻¹) which do not contain a primordial stellar population. For runs C and D, we find the number of luminous satellites increases linearly with host mass. For hosts with M < 10¹¹ M_☉, we see a larger scatter in the total number of satellites. Additionally, in that host mass range, we see greater scatter in the mapping of the total number of satellites to the number of luminous satellites. Since runs C and D contain only three Milky Way mass systems, the lack of scatter may also be due to small number statistics. The decrease in scatter may be a function of how dominant the halo is in its environment. In the filaments, a 10¹² M_☉ halo dominates the region around it, negating any environmental effects inside the virial radius. A lower mass host, however, is not able to dominate its environment. Therefore, the number of satellites for low-mass hosts will be more sensitive to the environment in which they are embedded.

The current observational sample of dwarfs is complete only to within 50 kpc (Simon & Geha 2007; Koposov et al. 2008; Walsh et al. 2009). The corrections for the detection limits of current surveys were done from a theoretical perspective by Tollerud et al. (2008). They used halos from Via Lactea I (Diemand et al. 2007), assuming a simple relationship between halo mass and luminosity for the subhalos. The range of Tollerud et al. (2008) is shown in Figure 5 as the shorter, thick black line. Unlike their work, our simulations do not assume a relationship of luminosity to halo mass. Instead, we draw the stellar properties of the z = 0 halos directly from the cosmologically consistent pre-reionization simulations. This accounts for the large scatter in stellar mass as a function of halo mass for the smallest galaxies (Ricotti et al. 2002b). Our results are consistent with the upper end of the Tollerud et al. (2008) range for the number of luminous satellites within ∼200 kpc. Based on these comparisons, the total number of subhalos and number of luminous satellites around MW.1, MW.2, and MW.3 are in agreement with results of other published works.

We next compare the distribution of maximum circular velocity for all subhalos around a Milky Way for our simulations with other CDM simulations. We find that the satellite mass functions for halos from runs C and D are consistent with one another, and results from Aquarius, Via Lactea, and Polisensky & Ricotti (2010, Figure 6). Based on this, we argue that our simulations can reproduce the number and distribution of subhalos around the Milky Ways as well as traditional N-body simulations. In the next section, we discuss the observational and theoretical criteria for a halo to be defined as a fossil of the first galaxies.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** Number of satellites within R_vir with greater than a given v_max for our two second-order (black lines) and one first-order (gray line) Milky Ways. The two versions of our method produce equivalent distributions and match the CDM simulations from Polisensky & Ricotti (2010).
Download figure:
Standard image High-resolution image

4. RESULTS

4.1. Definition of a Fossil Dwarf

For observed dwarfs, a fossil is defined as a dSph that underwent ≳ 70% of its star formation before reionization and today is a diffuse, spherical system devoid of gas (Ricotti & Gnedin 2005). These dim dwarfs populate dark matter halos whose circular velocities have never been above the filtering velocity, preventing them from accreting gas from the IGM after reionization. In our simulations, we define a fossil halo for which max(v_max(z)) < v_filt. Any halo with v_max(z = 0) < v_filt is referred to as a candidate fossil. However, in regimes where tidal stripping is considerable, there is a significant chance that a halo with a v_max < v_filt at z = 0 had a maximum circular velocity above the threshold for accretion from the IGM at an earlier time (Kravtsov et al. 2004).

Given these criteria, we classify our z = 0 halos into three populations as follows. (1) A non-fossil is a z = 0 halo for which v_max(z = 0) > v_filt. (2) Halos which are candidate fossils but for which max(v_max(z)) was above the IGM accretion threshold in the past are classified as polluted fossils. The non-fossils and a fraction of the polluted fossils accreted gas from the IGM and formed a significant population of stars after reionization. Therefore, our simulations cannot provide robust information on the non-fossil and polluted fossil stellar properties in the modern epoch. (3) For the true fossils we are able to generate detailed information on their stellar properties. A true fossil is defined as any z = 0 halo for which v_max never exceeded the IGM filtering mass, suppressing gas accretion and star formation after reionization.

To separate the polluted fossils from the true fossils of the first galaxies, we follow the v_max evolution for each candidate fossil back from z = 0 to z_init. We find that f(v_max), the fraction of candidate fossils which have max(v_max) > v_filt, as a function of their v_max(z = 0), is consistent with results found by Kravtsov et al. (2004; see Figure 7). In addition, we find that f(v_max) does not have a strong dependence on the environment of the fossils. When we compare the results for all the fossils (solid line) with those within 1 Mpc (dotted line) and 400 kpc (dashed line) of MW.2 and MW.3 we do not see a significant difference. These results are independent of the choice of the filtering velocity. For the remainder of this work, we use the term fossil in reference to only these true fossils.

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** Fraction of candidate fossils with max(v_max(z)) > v_filt, where v_filt = 20 km s⁻¹ for run D (lines) and Kravtsov et al. (2004, asterisks). The solid, dashed, and dotted lines show the fraction of true fossils for three different sub-populations. The solid line shows the relation for all the candidate fossils in run D, while the dashed and dotted lines show the fraction of true fossils for candidate fossils within 1 Mpc and 400 kpc of MW.3, respectively.
Download figure:
Standard image High-resolution image

In addition to maintaining v_max < v_filt for its entire evolution, a fossil must also survive to z = 0 without being tidally stripped. While objects that have undergone tidal stripping are unlikely to retain their pre-reionization stellar properties (Peñarrubia et al. 2008), our initial conditions do not allow us to simulate tidal effects beyond the stripping of a z = 0 halo's tracer particles. The use of N-body particles to represent pre-reionization halos forces the masses of those halos to be conserved. No matter how strong the tidal forces are, the stellar and dark matter properties will not change. This is in no way consistent with the current understanding of the effect of tidal stripping on a satellite's stellar population.

While the dark matter halo can be stripped away, leaving the stellar properties relatively intact (Choi et al. 2009; Peñarrubia et al. 2008), once 90% of the dark matter has been stripped and the mass loss reaches the outer stellar radii, the stripping of the stellar populations will occur at a faster rate than the denser dark matter cusp increasing the mass-to-light ratio of the system (Peñarrubia et al. 2008). We have no way of tracking the mass loss of an isolated pre-reionization halo to determine which components have been disrupted. We therefore err on the side of caution: we use the destruction of a z = 0 halo's dark particle cloud to flag halos which have undergone tidal stripping. Any present-day halo whose cloud of tracer particles has been destroyed or stripped down to N ≲ 50 particles will not be robustly detected as substructure and its mass will be added to that of the host galaxy. If N < 20, the z = 0 halo will not be detected at all (Knollmann & Knebe 2009). This adds a second, lower mass criteria for a pre-reionization halo to be identified as part of a fossil. It must be "found" in a halo at z = 0 by AHF to be considered a fossil. Any pre-reionization halo not in a z = 0 halo is assumed to be completely disrupted.

Given these criteria, we can say a few things about our fossil population. Our fossils are dimmer and less massive than the polluted fossils and non-fossils. As a population, they are less likely to have undergone mergers involving two or more luminous pre-reionization halos (Figure 8). We define a merge between two or more luminous pre-reionization halos as a galaxy merger. Using v_filt = 20 km s⁻¹, 25% of the fossils have two or more luminous pre-reionization halos compared to 40% of candidate fossils. The majority of true fossils (75%) contain only one luminous pre-reionization halo; however, the remainder do not represent a negligible fraction. We find the same result when using the v_filt = 30 km s⁻¹ adopted by GK06. As with the 20 km s⁻¹ case, 75% of true fossils have only one luminous pre-reionization halo. Therefore, while the majority of fossils have not undergone galaxy mergers, it is not an effect that can be ruled out.

Figure 8. Refer to the following caption and surrounding text. — **Figure 8.** Fraction of luminous true fossils which have undergone <6 galaxy mergers after reionization for v_filt = 20 km s⁻¹ (black) and v_filt = 30 km s⁻¹ (red). We define a galaxy merger as any merger in which two or more of the components contain a luminous population. For >4–5 galaxy mergers, the fraction of z = 0 true fossils becomes negligible.
Download figure:
Standard image High-resolution image

4.2. Luminosity Threshold for Fossils

Before making detailed comparisons between our simulations and observations, we compare our work and the N-body simulations in GK06. Unlike our method, which allows us to directly trace the pre-reionization halos to the present day, GK06 statistically matches pre-reionizaion halos to their counterparts at z = 0 based on their v_max at z = 8.3. To make a direct comparison with GK06 we must use our MW.1 from run C, since GK06 only used the z = 8.3 outputs from the pre-reionization simulations.

Figure 9 shows the galactocentric radial distribution for GK06 (blue band) and for MW.1 (black lines). Both curves only include the true fossils. For L_V > 10⁵ L_☉ (lower panel) we find that our simulations are consistent with GK06, if on the low end of their range. However, the brightest simulated true fossils in GK06 with L_V > 10⁶ L_☉ have no counterparts around MW.1. We ascribe this discrepancy to the difference in how our work follows the pre-reionization halos to the modern epoch.

Figure 9. Refer to the following caption and surrounding text. — **Figure 9.** Radial distribution of the true fossils around MW.1 in run C (black lines) and the results from GK06 (blue band) for halos with *L_V* > 10⁵ L_☉ and *L_V* > 10⁶ L_☉. We have used a v_filt = 30 km s⁻¹ threshold to determine whether a z = 0 halo is a true fossil.
Download figure:
Standard image High-resolution image

While both methods allow for the growth and stripping of a halo via accretion and tidal forces, our simulations also account for clustering of the pre-reionziation halos. The most luminous pre-reionization halos correspond to the most massive halos at z = 8.3. These 10⁷–10⁸ M_☉ galaxies are preferentially located in higher density regions within the 1 Mpc³ pre-reionization simulation. This increases the probability that the pre-reionization halos with L_V > 10⁵ L_☉ will have undergone a galaxy merger relative to those with L_V < 10⁵ L_☉. In Figure 10, we show the histogram of the number of luminous pre-reionziation halos for true fossils with L_V < 10⁵ L_☉ (left panel) and L_V > 10⁵ L_☉ (right panel). Only ∼0%–5% of the highest luminosity fossils have never undergone a galaxy merger compared to ∼90% of fossils with L_V < 10⁵ L_☉. This is independent of our choice of filtering velocity.

Figure 10. Refer to the following caption and surrounding text. — **Figure 10.** Left: histogram of the fraction of true fossils at z = 0 with a given number of galaxy mergers after reionization for v_filter = 20 km s⁻¹ (black line) and v_filter = 30 km s⁻¹ (red line). As in Figure 8, the number of galaxy mergers is a proxy for the number of luminous pre-reionization halos in a z = 0 halo. In this panel, we show only the true fossils with *L_V*(z = 0) < 10⁵ L_☉. Right: the fraction of true fossils with a given number of luminous pre-reionization halos for only those with *L_V*(z = 0) > 10⁵ L_☉ for v_filter = 20 km s⁻¹ (black line) and v_filter = 30 km s⁻¹ (red line). Note the shifted peak and different shape of the histogram in this panel.
Download figure:
Standard image High-resolution image

Why does this explain the discrepancy between our results and GK06 in Figure 9? The definition of a true fossil is a dwarf whose maximum circular velocity has never gone above the threshold for accretion for the IGM. In Figures 9 and 10, we set v_filt = 30 km s⁻¹. Since the brightest pre-reionization halos are also the most massive, one or two galaxy mergers at high redshift would be enough to push v_max above the filtering velocity and classify the halo as a non-fossil. In run C, there are only 11 true fossils with L_V > 10⁶ L_☉, none of which are within 1 Mpc of MW.1.

This gives us a maximum luminosity threshold, 10⁶ L_☉, above which an observed dwarf is unlikely to be a primordial fossil. Of the true fossil candidates identified in RG05, this puts seven into question: And I (4.49 × 10⁶ L_☉), And II (9.38 × 10⁶ L_☉), And III (1.13 × 10⁶ L_☉), And VI (2.73 × 10⁶ L_☉), Antila (2.4 × 10⁶ L_☉) and KKR 25 (1.2 × 10⁶ L_☉) around M31, and Sculptor (2.15 × 10⁶ L_☉) around the Milky Way. The remaining seven, And V, Cetus, Draco, Phoenix, Sextans, Tucana, and Ursa Minor, all have L_V < 10⁶ L_☉ and remain reasonable candidates for the fossils of the first galaxies. The classical Milky Way fossils, Sextans and Ursa Minor, as well as the ultra-faint Canes Venatici I, all have metallicity distributions suggesting star formation durations <1 Gyr and populations >10 Gyr old (Kirby et al. 2011a). These dwarfs, in addition to Sculptor, have star formation histories that are dominated by outflows, in contrast to their brighter counterparts ("polluted fossils" Fornax and Leo I & II; Kirby et al. 2011b). However, unlike the other outflow-dominated dwarfs which have relatively short star formation bursts, Sculptor has undergone star formation over several Gyrs (Babusiaux et al. 2005; Shetrone et al. 2003; Tolstoy et al. 2003) and a fraction of Draco's stars may be of intermediate age (Cioni & Habing 2005).

4.3. Fossil Properties

In this section, we present the stellar properties of our simulated true fossils and compare them with observed stellar properties of Milky Way satellites. These comparisons include V-band luminosity, L_V, half-light radius, r_hl, metallicity, [Fe/H], and mass inside the half-light radius M_1/2 (Walker et al. 2009; Wolf et al. 2010). In BR09, we showed strong statistical agreements between the stellar properties of the pre-reionization halos and the observed distribution of known classical dSph and ultra-faint dwarfs. Here we improve our previous results by relaxing some of the assumptions made in BR09.

As in GK06, BR09 assumed that none of the luminous pre-reionization halos had undergone a galaxy merger. Thus, the present-day distribution of stellar properties for the fossils would be identical to that of the pre-reionization halos. In addition, our previous work assumed that the voids were reheated to T ∼ 10⁴ K well after the clusters and filaments, as expected for UV reionization by stars. As seen in Figure 2, a universe reionized first in the clusters and filaments and then in the voids, run C, produces a larger number of luminous objects in the voids when compared to a universe reionized at the same time by redshifted X-rays from primordial black holes (Ricotti & Ostriker 2004; Ricotti et al. 2005; run D). As in BR09, for all observed stellar properties, we use the measurements with the lowest error bars.

From hierarchical formation models, we know that all halos have undergone merger/accretion events since their epochs of formation. For 60% of our pre-reionization halos, these mergers are with dark halos, producing a daughter halo with the same stellar properties as the parent. However, for all runs and all halos, regardless of their fossil status, ∼40% of the z = 0 halos contain more than one luminous pre-reionization halo. These galaxy mergers will change the stellar properties of the systems.

True fossil halos in the modern epoch derive their stellar properties solely from their pre-reionization populations. For the 75% of luminous true fossils that contain only one luminous pre-reionization halo, the z = 0 stellar properties are taken directly from those of the pre-reionization halo. We account for the reddening of the stellar population by using M^rei_*/L ∼ 5. Note that we use such a large stellar mass to light ratio to account for stellar mass lost since reionization. The stellar mass to light ratio of our simulated galaxies at z = 0 is

$\begin{equation} \frac{M_\ast ^{{\rm rei}}}{L} = \left(\frac{M_\ast ^{{\rm today}}}{L}\right) \Bigg(\frac{M_\ast ^{{\rm rei}}}{M_\ast ^{{\rm today}}}\Bigg), \end{equation} \tag{ 2 }$

where M^rei_* is the mass of the stellar population at reionization and M^today_* is the mass of the stellar population at z = 0. The ratio between them, M^rei_*/M_*^today, is between 2 and 20, depending on the IMF of the primordial stellar population. RG05 used a range of M/L ratios and found no dependence of the fossil properties on the choice of mass-to-light ratio.

For the one-quarter of true fossils that have undergone a galaxy merger, the stellar properties are calculated as follows. Throughout this section, the superscript f will denote the stellar and dark matter properties of the z = 0 halo, and the superscript i the properties of the component, luminous pre-reionization halos.

The final V-band luminosity, L^f_V of a fossil halo at z = 0, is the sum of the V-band luminosities, Lⁱ_V, of the component pre-reionization halos. We assume that stellar mass is conserved during all mergers, an assumption that will be addressed in future, higher resolution simulations.

We determine the half-light radii, r^f_hl, for z = 0 fossils using the three-dimensional r_hl from the pre-reionization simulations, with the following assumptions. (1) The dynamical evolution of the stars is decoupled from that of the dark matter. (2) The kinetic energy of the stars is conserved. (3) The collision of the luminous pre-reionization halos is elastic with respect to the stars. (4) Enough time has passed since the collision for the halo to return to an equilibrium state. Given the kinetic energy conservation of the stars

$\begin{equation} (\sigma _\ast ^f)^2 = \big(L_V^f\big)^{-1} \sum L_V^i \times (\sigma _{\ast }^i)^2, \end{equation} \tag{ 3 }$

where σⁱ_* and σ^f_* are the three-dimensional stellar velocity dispersions of the parent and daughter halos. For a halo in equilibrium, r_hl ∼ σ²_*; therefore,

$\begin{equation} r_{{\rm hl}}^f = \big(L_V^f\big)^{-1} \sum L_V^i \times r_{{\rm hl}}^i. \end{equation} \tag{ 4 }$

We use r^f_hl to calculate an average surface brightness, 〈Σ_V〉, for our fossils in units of L_☉ pc⁻². The Σ_V and r_hl distributions as a function of luminosity are shown in Figure 11.

In Figure 11, the black symbols are the observed Milky Way and M31 satellites overlaid on colored contours showing the equivalent distributions for the simulated true fossils. The cyan and red contours show the stellar properties of the fossils above and below the Sloan Digital Sky Survey (SDSS) detection limits, respectively. We see that, as in BR09, our simulations are able to reproduce the observed Σ_V and r_hl distributions for the ultra-faint and classical dSphs, with a few exceptions. We are unable to account for the ultra-faints with r_hl < 60 pc (Coma Berenics, Segue 1 and 2, Leo V, and Willman 1), all but one of which (Leo V) are within ∼50 kpc of the Milky Way.

In BR09 we called attention to an as yet undetected population of ultra-faints with surface brightnesses below SDSS limits. The existence of these dwarfs was independently proposed in Bullock et al. (2010), who named them "stealth galaxies." The detection of these ultra-faint dwarfs is a test for the fossil scenario. In this section, we summarize the properties expected of these extremely ultra-faint fossils.

In Figures 11–16, the simulated true fossils are shown as two sets of contours. Up until now, we have been comparing the ultra-faints and a subset of the classical dSphs to the simulated true fossils with Σ_V > 10^−1.4 L_☉ pc⁻². These true fossils, shown by the cyan contours, would be detectable by the SDSS (Koposov et al. 2008). The red contours show the true fossils which would remain undetected by SDSS. In Section 5, we present the existence and properties of the true fossils with surface brightnesses below the SDSS detection limits as a test for primordial star formation in minihalos. For the remainder of this section, we direct the reader to the red contours in Figures 11–16.

Figure 11. Refer to the following caption and surrounding text. — **Figure 11.** Left: Figure 1 from BR09. Surface brightness and half-light radii are plotted against V-band luminosity. The small black squares show the properties of the pre-reionization halos at z = 8.3. The other black symbols show the dwarf populations for the Milky Way and M31. The asterisks are non-fossils (dIrr), crosses are polluted fossils (dE and some dSph), the filled circles and triangles are the fossils (dSph) known before 2005 for the Milky Way and M31, respectively, and the open circles and triangles are the ultra-faint populations those galaxies found since 2005. Right: surface brightness and half-light radii are plotted against V-band luminosity. The cyan contours show the distribution for the fossils from run D and the overlaid black symbols show the observed dwarfs. In this panel, we color the observed dwarfs whose half-light radii are inconsistent with our simulations green. The magenta contours show the undetectable fossils with Σ_V below the zeroth-order detection limit of the SDSS, ∼ − 1.4, (Koposov et al. 2008). In both panels, the solid black lines show the surface brightness limit of the Sloan (Koposov et al. 2008) and the dashed black lines show the trends from Kormendy & Freeman (2004) for luminous Sc–Im galaxies (10⁸ L_☉ < *L_B* < 10¹¹ L_☉).
Download figure:
Standard image High-resolution image

Figure 12. Refer to the following caption and surrounding text. — **Figure 12.** Left: the stellar mass to light ratios calculated from Illingworth (1976) and stellar velocity dispersions versus the V-band luminosities for run D (blue contours) and observations (red symbols). Symbols are the same as in Figure 11. Once again, the dashed lines show the Kormendy & Freeman (2004) trends for Sc–Im galaxies with 10⁸ L_☉ < *L_B* < 10¹¹ L_☉. Right: the M/1/2L_V vs. 1/2L_V using the Walker et al. (2009) mass estimator.
Download figure:
Standard image High-resolution image

Figure 13. Refer to the following caption and surrounding text. — **Figure 13.** Left: the mass function of the detected (solid) and undetected (dashed) fossils with *L_V* > 10² L_☉ within 1 Mpc of MW.2 (top) and MW.3 (bottom) from run D. Right: same as the left panels except the x-axis is the dynamical mass inside the half-light radius (Walker et al. 2009) calculated from the velocity dispersion and half-light radius of our fossils.
Download figure:
Standard image High-resolution image

Figure 14. Refer to the following caption and surrounding text. — **Figure 14.** [Fe/H] distribution for the true fossils plotted against the V-band luminosities for run D (blue contours) and observations (red symbols). Symbols are the same as in Figure 11. Our results agree with Salvadori & Ferrara (2009) for *L_V* > 10⁴.
Download figure:
Standard image High-resolution image

Figure 15. Refer to the following caption and surrounding text. — **Figure 15.** Maximum circular velocity, v_max of our simulated true fossils plotted against the V-band luminosities. The cyan and red contours are the same as in Figure 11. Here we show no observed dwarfs due to the lack of data.
Download figure:
Standard image High-resolution image

Figure 16. Refer to the following caption and surrounding text. — **Figure 16.** Stellar velocity dispersion, σ_* against the half-light radius, r_hl. The black symbols are the observed dwarfs and blue and red contours from run D have the same meanings as in Figure 11.
Download figure:
Standard image High-resolution image

The mass-to-light ratios and σ_* of the observed and simulated populations are shown as the top and bottom parts of the left panel in Figure 12. As in Figure 11, the five dwarfs that do not match the r_hl of the simulated fossils are marked with filled green circles. Excepting this subpopulation, the ultra-faints show the same distribution as the simulated true fossils for both M/L and stellar velocity dispersions. In the left panel of Figure 12, the masses of our simulated halos are calculated from the stellar properties using Illingworth (1976). The right panel shows the mass-to-light ratios inside the half-light radii, M/L_1/2 versus half the V-band luminosity using the Walker et al. (2009; Wolf et al. 2010) mass estimator. The latter mass estimator is more accurate for dispersion supported systems, but we note that the agreement between the mass-to-light ratios of our fossils and ultra-faint dwarfs is independent of the mass estimator we use to calculate M(σ_*, r_hl). As expected, the undetected dwarfs (red contours in Figure 12) would have M/L > 10³M_☉/L_☉, higher than even the most dark matter-dominated ultra-faint dwarfs. However, the range of their stellar velocity dispersion is 2–10 km s⁻¹, equivalent to the ultra-faint dwarfs and detectable fossils, and shows no evolution with decreasing luminosity.

As seen in the left panel of Figure 13, the mass functions of the detected and undetected fossils peak at 10⁸ M_☉. Note, however, that this mass function is for the total dark matter mass, not the dynamical mass calculated from the velocity dispersion and half-light radius. Our simulations provide us with the information needed to plot a mass function of the dynamical mass, referred to in the right panel of Figure 13 as the observed mass. For the derived observational mass function both the detected and undetected fossils peak at 2 × 10⁷ M_☉. This peak corresponds to the "common mass scale" for dwarf spheroidals (Strigari et al. 2008); however, no such sharp peak is seen in the total dark matter mass function. This is because local feedback effects in small halos produce a large scatter in the relationship between the total and stellar masses of halos. Two dwarfs with similar mass and extent of their stellar populations, and thus similar dynamical mass within the luminous radius, may be embedded within halos whose masses vary by an order of magnitude. While the more massive halo's stellar population is concentrated at the center of its potential, the less massive halo's stars fill a larger fraction of its dark matter halo. This produces either two halos with the same dark matter mass and different M₃₀₀(r_hl, σ_*) or vice versa, a common M₃₀₀ but very different dark matter masses.

For the metallicity distribution, we also use a luminosity-weighted average:

$\begin{equation} [{\rm Fe/H}]^f = \log \left(\sum 10^{[{\rm Fe/H}]^i} \times L_V^i\right) - \log \big(L_V^f\big). \end{equation} \tag{ 5 }$

The distribution of metallicity versus L_V is shown for our z = 0 fossils and the known ultra-faint and classical dwarfs. As in BR09, the fossil metallicities from run D are consistent with the observed distribution for the ultra-faint and classical dSph. We also find our results for L_V > 10⁴ L_☉ to be in agreement with Salvadori & Ferrara (2009) while for the dimmest fossils our work finds comparatively lower metallicities. The undetected dwarfs (red contours in Figure 14) have [Fe/H] < −2.5 and as low as −3.5 with slightly larger scatter than their detectable counterparts.

The maximum circular velocity versus L_V contours for our simulated true fossils are shown in Figure 15 to illustrate the following. While v_max does decrease by approximately a factor of two over four decades of luminosity, the scatter in v_max at a given L_V is large. Though a halo with v_max < 6 km s⁻¹ is likely to have a L_V < 10⁴ L_☉, there is, at most, a minimal trend of decreasing v_max with decreasing luminosity for the primordial fossils. This highlights a theme across all our stellar property comparisons. Because of the strong dependence of their stellar properties on stochastic feedback effects, there is no baryonic property that shows a strong trend with maximum circular velocity and the size of the dark matter halo.

We now briefly discuss the M31 satellite population. Figure 16 shows the σ_* plotted against r_hl on a scale similar to the top left panel of Figure 18 in Collins et al. (2010). The circles show the Milky Way dSphs, while the triangles show the dSphs associated with M31. We find that four of the six M31 dSphs plotted are within, albeit at the edges of, the contours of detectable true fossils. Like their Milky Way counterparts, the new M31 dSphs show reasonable agreement with our simulated primordial fossils except for r_hl, which are higher than expected from our simulations for two of the M31 dwarfs. However, this does not represent a major problem for our model since ∼65% of simulated true fossils with L_V > 10⁴ L_☉ have undergone one or more major mergers that may have puffed up their stellar populations. Our estimates do not account for extra heating of the stellar populations by the kinetic energy of the collision. A higher σ_* would result in a more extended stellar population in the same mass halo. We will discuss the comparison between the M31 dSphs and our simulated fossil dwarfs in an upcoming paper.

4.4. The Inner Ultra-faints

In this section, we discuss a possible origin scenario for the inner ultra-faint dwarfs, i.e., the ultra-faints whose half-light radii and mass-to-light ratios are lower than our true fossils. These dwarfs are Segue 1 and 2, Leo V, Pisces II, and Willman 1, and excepting Leo V and Pisces II (both at ∼180 kpc), all are within 50 kpc of the Milky Way. This classification of the ultra-faint Milky Way dwarfs in shown in Table 3.

Table 3. Milky Way Ultra-faint Dwarfs Classified by Their Distance from Our Galaxy (Columns) and Whether or Not They are Consistent with Our Predictions for the Fossils of the First Galaxies (Row)

	R < 50 kpc	R > 50 kpc
Inconsistent	Segue 1	Pisces II *
	Segue 2	Leo V *
	Willman 1

Consistent	Coma Ber.	Bootes I & II
		CVn I & II
		Hercules
		Leo IV & Leo T
		Ursa Major I

Notes. Note the correlation between distance and consistency. (*) Pisces II and Leo V are both on the lower end of radii expected for fossils, as such they are marked as "inconsistent," but are not as far from predictions as the "inconsistent" ultra-faints within 50 kpc.

Download table as: ASCII Typeset image

However, their mass-to-light ratios follow a shifted power law with a similar slope to the true fossils and more luminous dwarfs. The stellar velocity dispersions are in the range expected for primordial fossils, but the inner ultra-faints show an L_V–σ_* combination which would be expected for true fossils below the detection limits of SDSS (red contours in Figure 12). These properties are either directly affected by tidal stripping (r_hl and σ_*) or are derived from affected properties (Σ_V and M/L). However, the metallicity of the stars is not affected by tidal stripping.

Figure 14 shows that the metallicities of the inner ultra-faint dwarfs do not fall on the luminosity–metallicity relation. However, their scatter is consistent with expectations for true fossils. To place the Segues, Leo V, Pisces II, and Willman 1 on the luminosity–metallicity relation traced by the majority of the ultra-faints and our fossils, their luminosities would need to be increased by one to two orders of magnitude. We suggest these dwarfs may be a subset of bright primordial fossils that have been stripped of 90%–99% of their stars.

5. OBSERVATIONAL TESTS

In this section, we present a set of observational tests that can provide support for a primordial formation scenario for the faintest Milky Way satellites.

Better determination of whether the ultra-faints are being tidally disrupted can help determine whether a subset of the faintest Milky Way satellites are pristine fossils. The ultra-faint dwarfs whose r_hl do not match our simulated true fossils, shown as filled green circles in Figures 11–16, display signs of being tidally disrupted by the Milky Way, including proximity to the Milky Way (R < 50 kpc). While Willman 1, Segue 2, Leo V, and Coma Berenticis show signs of tidal disruption, this does not prove the primordial scenario. However, it would place their origin as disrupted objects in line with our proposal in Section 4.4. If additional observations show these tidal ultra-faints are not tidally disrupted, then our primordial formation model cannot explain their current properties. The exception to this picture is Segue 1. The tidal status of Segue 1 has been recently debated (Niederste-Ostholt et al. 2009; Norris et al. 2010; Martinez et al. 2011), and Simon et al. (2011) show that the kinematics of the stars are consistent with a stellar population well within the tidal radius. The high density of the stars suggests that Segue 1 formed at high redshift in a rare, high σ peak which our 1 Mpc³ pre-reionization volume is not large enough to include.

The number of Milky Way satellites alone provides a test for star formation in minihalos. For a given filtering velocity, there is a number of satellites, N_nf, which has a v_max at infall above the filtering velocity. For v_filt = 20 km s⁻¹, N_nf is 90 ± 10 and for v_filt = 30 km s⁻¹, N_nf is 60 ± 8; the latter is equivalent to the number of currently known Milky Way satellites. If the number of satellites, N_sat, is greater than N_nf, some minihalos had to have formed stars before reionization. Conversely, if N_sat < N_nf, either no minihalos formed stars or none survived near the Milky Way.

We next outline the stellar properties we can expect of the undetected dwarfs around the Milky Way if they are part of a population of fossils of the first galaxies. The red contours of Figures 11–16 show the properties of the predicted population.

1.
Half-light radii. The undetected dwarfs should have the same distribution of half-light radii as the currently known, ultra-faint population, from ∼100 pc to ∼1000 pc.
2.
Mass-to-light ratio. The mass-to-light ratio of undetected dwarfs should generally be greater than 10³M_☉/L_☉ and as high as approximately 10⁵M_☉/L_☉ and follow a roughly linear relation for dwarfs with L_V < 10⁵ L_☉.
3.
Stellar velocity dispersion. There should be no decrease in the stellar velocity dispersion, σ_*, with V-band luminosity. This directly contradicts the decreasing σ_* with L_V seen for tidally stripped dwarfs in Figure 4 of Wadepuhl & Springel (2011).
4.
Metallicity. The undetected dwarfs should have typical [Fe/H] < −2.5, a significant number with [Fe/H] < −3.0. This would make these undetected dwarfs excellent candidates for the search for ultra-metal-poor stars (Frebel et al. 2010; Frebel & Bromm 2010).

6. SUMMARY AND CONCLUSIONS

In this paper, the first of a series, we present a new method for generating initial conditions for cosmological N-body simulations, which allows us to create simulated maps of the present-day distribution of fossils in a "Local Volume." In order to produce these maps, we assume pre-reionization fossils do not accrete gas and form stars after reionization. They are hosted in dark halos that maintain circular velocities below a critical threshold, v_filt ∼ 20–30 km s⁻¹. The precise value of v_filt depends on details of the reheating in the local IGM by stars and active galactic nuclei. Therefore, we explore different values for v_filt. For our purposes, we do not need to include gas dynamics. The lack of post-reionization baryonic evolution in the fossils allows us to simply simulate the evolution of the dark matter and stars using N-body techniques.

We have combined the results from previous cosmological simulations of the formation of the first galaxies (Ricotti et al. 2002b; Ricotti & Gnedin 2005; Ricotti et al. 2008) with N-body simulations in which each particle in the initial conditions represents a pre-reionization minihalo. Our N-body simulations zoom in on a Local Volume containing one to two Milky Ways. We follow the merger history and tidal stripping of pre-reionization fossils as they merge to form more massive galactic satellites of the Milky Way. We also trace the evolution of more massive non-fossil satellites, but we do not account for star formation taking place after reionization.

Our goal is to determine if a widespread population of primordial dwarfs is consistent with the observed population of Milky Way and Andromeda satellites, and, at the same time, matches observations of dwarfs in the Local Void. It is not well established whether halos with v_max < 20 km s⁻¹, too small to initiate collapse via Lyα cooling, remain dark or form luminous dwarfs. Our simulations are a first attempt to constrain the theory of self-regulated galaxy formation before reionization using "near field" observations. Observational tests based on our results can constrain models of star formation in minihalos before reionization.

In this paper, we present and validate our method by comparing our results with published cosmological N-body simulations. We then present maps of the Local Volume showing the distribution of stars formed before reionization in the present-day universe. We find that primordial fossils are present in the voids regardless of the details of reionization; however, reionization by X-rays produces darker voids. Finally, we show plots of the present-day properties of true fossils and compare them to observations of classical and ultra-faint dwarfs for the Milky Way and Andromeda. We show that the simulated properties of fossils agree with a subset of ultra-faint dwarfs discovered in the Milky Way and Andromeda, and reiterate the BR09 result of a large population of fossils with surface brightness below the SDSS detection limits. The properties of this "stealth" population are shown as red contours in Figures 11–16.

We find that most classical dSph satellites are unlikely true fossils of the first galaxies, even though they have properties expected of fossils: diffuse, old stellar populations with no gas (Ricotti & Gnedin 2005; Bovill & Ricotti 2009). The reason that true fossils in the Milky Way have luminosities <10⁶ L_☉ is that the most luminous pre-reionization fossils with v_max ∼ 20 km s⁻¹ form in overdense regions and are strongly clustered. Thus, they are likely to merge into more massive non-fossil dwarfs. The surviving fossils found today are a sub-population with lower typical luminosities, and form in less clustered regions in which feedback effects suppress rather than stimulate star formation.

Other results found in this first paper are highlighted below.

1.
Voids contain many low luminosity fossil galaxies. However, they have surface brightnesses and luminosities that make them undetectable by SDSS. One possible way to detect these void dwarfs is if they experience a late phase of gas condensation from the IGM as proposed in Ricotti (2009). Future and present 21 cm surveys such as ALFALFA and GALFA may be used to find these objects (Giovanelli et al. 2005; Begum et al. 2010).
2.
We find a linear scaling relation between the number of luminous satellites and the mass of host halos. The scaling has scatter similar to the relationship between the total number of sub-halos with M > 10⁷ M_☉ (v_max > 5 km s⁻¹) and the host mass, although the normalization is 3–4 times lower.
3.
Overall, ∼25% (for v_filt = 20 km s⁻¹) to ∼30% (for v_filt = 30 km s⁻¹) of the primordial fossils at the present day have undergone a merger with another luminous fossil. This fraction increases with the modern luminosity of the dwarf. Hence, the typical half-light radii of this population can be larger than the original distribution at reionization. These fossils are even harder to detect due to their lower surface brightness. This effect also increases the spread of the relationship between half-light radii versus luminosity and surface brightness versus luminosity of fossils at z = 0.
4.
Due to the dependence of the properties of primordial dwarfs on their formation environment (Ricotti et al. 2008), we find very few true fossils with L_V > 10⁶ L_☉, and none within 1 Mpc of our Milky Ways. This places the identification of some of the more luminous classical dSphs fossils in doubt.
5.
Leo V, Pisces II, Segue 2, and Willman 1 have half-light radii which are too small, and metallicities too large, for their luminosities. Due to their proximity to the Milky Way, we speculate that their stars and dark halos have been affected by tides. Hence, these ultra-faints may represent a population of massive primordial dwarfs which have lost ≳ 90% of their stars via tidal interactions.
6.
We reiterate the existence of a yet undetected population of fossils with luminosities L_V < 10⁴ L_☉ and surface brightness <10^−1.4 L_☉ pc⁻². We present plots showing, in detail, the expected properties of this population. We also note that some of the new ultra-faint satellites in Andromeda have half-light radii in agreement with the properties of the undetected fossil population, but luminosities ∼10⁵ L_☉.

The second paper of this series (Bovill & Ricotti 2011) will focus on studying the galactocentric radial distribution and cumulative luminosity functions of simulated fossil and non-fossil satellites and compare the results to the Milky Way dwarfs. Although we have not unmistakably demonstrated the existence of fossils around the Milky Way using existing observations, we propose new observational tests already feasible with current instruments that should tell us whether pre-reionization fossil formation was widespread in the early universe.

In a future paper, we will run simulations in which we follow each star and dark matter particle in our pre-reionization simulations to z = 0. These higher resolution simulations will allow us to study tidal stripping of stars and dark matter for halos within 50 kpc from the Milky Way and better understand the nature of the faintest ultra-faint dwarfs. This method will also allow us to better quantify the stellar properties of non-fossil dwarf satellites formed by mergers of fossils.

The simulations presented in this paper were carried out using computing clusters administered by the Center for Theory and Computation of the Department of Astronomy at the University of Maryland ("yorp") and the Office of Information Technology at the University of Maryland ("hpcc"). This research was supported by NASA grants NNX07AH10G and NNX10AH10G. The authors thank the anonymous referee for constructive comments and feedback. M.S.B. and M.R. also thank Stacy McGaugh, Susan Lamb, Rosie Wyse, and Derek Richardson for helpful conversations and comments on this work.

APPENDIX A: FIRST-ORDER INITIAL CONDITIONS

To study the distribution and stellar properties of the fossils of the first galaxies across a region equivalent to the Local Volume, we need an N-body simulation with the dynamical range to resolve M < 10⁷ M_☉ halos through a volume ∼10 Mpc on a side. We also want to directly trace the present-day distribution of the halos from our pre-reionization simulations, while keeping the required computational resources reasonable.

We generate hybrid initial conditions in which the distribution of matter is governed by different mechanisms on different scales. As in traditional N-body simulations, structure on scales l > 1 Mpc is set by a power spectrum. However, on scales l < 1 Mpc, the positions, velocities, and masses of the particles are set by the Press–Schechter, via outputs from our pre-reionization simulations (Ricotti et al. 2002a, 2002b, 2008). These hybrid initial conditions are set as follows. (1) We locate an analog to the Local Volume within a coarse resolution 50³ Mpc³ volume run from z = 40 to z = 0. (2) A high-resolution region is built out of the final outputs from the pre-reionization simulations at z_init = 8.3. (3) Finally, we insert our high-resolution, "Local Volume" into the larger coarse resolution simulation at z_init = 8.3 and run it to z = 0 using Gadget 2 (Springel 2005). In the rest of Appendix A, we will explain these steps in more detail.

We need to generate and run a volume large enough to contain at least one subvolume analogous to the Local Volume. Our coarse resolution simulation is a 50³ Mpc³ volume with 250³ particles run from z = 40 to z = 0. The power spectrum at z = 40 is generated by P-GenIC. At z = 0, we use HOP (Eisenstein & Hut 1998) to locate potential Milky Ways. In our "Local Volume," we look for a filament between two Virgo-like clusters with two to three halos with M ∼ 10¹² M_☉ within a 7³–10³ Mpc³ volume. Ideally, one of our Milky Ways should have an equal mass companion within 1 Mpc; however, we were not able to find such a pair in our volume.

From the location of our Milky Ways, we define our "Local Volume" as a region ∼5–10 Mpc across, centered on one of our Milky Ways. Once we have a "Local Volume" at z = 0, we estimate its equivalent volume at z_init = 8.3. We do this by tagging the coarse resolution particles in our present-day Local Volume, and, using their positions at z = 8.3, determine the equivalent rectangular prism containing the majority of the tagged particles. At this point, we have defined a high-resolution region with dimensions m × n × p.

We turn the 1 Mpc³, z = 8.3 output from the pre-reionization simulations to an equivalent cube of N-body particles as follows. First, any pre-reionization halo in the output becomes an N-body particle with a position, velocity, mass and, critically, unique ID. We then choose a mass resolution, m_trun, for our high-resolution simulations and truncate the mass function of the pre-reionization halo at that resolution. To account for the additional mass needed to bring each 1 Mpc³ to the average density of the universe, we add dark tracer particles. These tracer particles have mass m_trace ≲ m_trun. The positions and velocities of the tracer particles are determined from the position and velocities of the pre-reionization halos with m < m_trun. At the end of this process we have a 1 Mpc³ cube of N-body particles with m_min ∼ m_trace, where positions, velocities, and masses are determined by the z = 8.3 output of the pre-reionization simulations. The mass function produced by this method is shown in the upper left panel of Figure 17. The spike in the lowest mass bin shows the mass of the tracer particles. Critically, each particle has a unique ID allowing us to trace each pre-reionization halo to z = 0 and retrieve its baryonic properties in the modern epoch.

Figure 17. Refer to the following caption and surrounding text. — **Figure 17.** Truncated mass function of the pre-reionization outputs from z = 8.3 (top left), z = 10.2 (top right), z = 12.1 (bottom left), and z = 14 (bottom right). The spike in the lowest mass bin at all four redshifts is due to the dark tracer particles.
Download figure:
Standard image High-resolution image

We have generated a 1 Mpc³ cube for which each particle is a tracer for a pre-reionization halo. From the location of our Local Volume, we have a rectangular prism at z = 8.3 where our high-resolution region will go. We duplicate the 1 Mpc³ box to form the m × n × p prism used for the high-resolution region. We now have an m × n × p prism with power on l < Mpc scales. We add power on l > Mpc scale using the position shift, δx, of the coarse resolution particles via linear interpolation between them. Once we have δx, we use the linear relation δv = A(z)δx to calculate the velocity perturbation, δv, for our high-resolution particles, where A(z) is the ratio of the δv/δx at a given redshift. We now have a high-resolution region with m × n × p embedded inside a 50³ Mpc³ simulation at z = 8.3. In the coarse resolution region, all power comes from the power spectrum and in the high-resolution region the power comes from the power spectrum on l > Mpc scales and from the Press–Scheter via the pre-reionization simulations on l < Mpc scales.

We then evolve the 50³ Mpc³ volume from z = 8.3 to z = 0 using Gadget 2 (Springel 2005) with the outputs analyzed by AHF (Knollmann & Knebe 2009). These simulations are runs A–C in Table 1.

APPENDIX B: SECOND-ORDER INITIAL CONDITIONS

GK06 and our runs A–C assume the formation and evolution of all structures in the universe occur at the same rate. However, this is not valid. The rate at which structure forms and evolves is dependent on the local overdensity or underdensity of the universe, δ (Cole 1997). The average density of the universe is δ = 0.0. Regions with δ < 0.0 form structures at a slower rate and vice versa for δ > 0.0. This difference affects the evolution of our primordial fossils in the following manner. Since dark matter halos will form at later times in the voids, they will accrete gas and begin to form stars at a lower z compared to the filaments and clusters. Star formation is shut off in these minihalos by the reheating of the universe during reionization. If this occurs at the same time everywhere (X-ray reionization; Ricotti & Ostriker 2004; Ricotti et al. 2005) fewer of the minihalos in the voids will have formed stars, and those that did will have lower luminosities.

Reionization begins in the clusters and filaments and spreads to the voids, reheating the lower density regions at lower redshift. Runs A–C would not need to account for this effect if the delay in halo and star formation in the voids exactly matched their later reheating. The second-order initial conditions in run D work in the model where the delay in reionizing and reheating in the voids is not enough to balance their later structure formation.

We begin the quantification of these different rates of evolution as a function of density, by assuming all parts of the universe evolve in a self-similar manner. Therefore, a 10⁸ M_☉ halo in a cluster will collapse at a higher redshift than a halo of same mass near a Milky Way and that halo will form at a higher redshift than a 10⁸ M_☉ halo in a void. We tie the δ to an effective redshift, z_eff:

$\begin{equation} \Delta z = (1+z_{{\rm init}})((1+\delta)^{-0.6}-1), \end{equation} \tag{ B1 }$

where z_eff = Δz + z_init, with z_init = z_eff(δ = 0.0), the current redshift of the entire universe, and Δz is the shift due to the local underdensity or overdensity. Using this relation, we link the local density of the universe to the rate of structure formation via an effective redshift. The distribution of z_eff for our 50³ Mpc³ volume is shown as the red curve in Figure 18. We now take the quantified relation between z_eff and δ, and apply it to our initial conditions. Before going any further, note that the selection of our Milky Ways and determination of our high-resolution region are the same for run D as for runs A–C.

Figure 18. Refer to the following caption and surrounding text. — **Figure 18.** Fraction of 1 Mpc³ cubes with a given effective redshift, z_eff. The red curve shows the distribution for all the Mpc³ cubes in our entire 50³ Mpc³, low-resolution volume. The black histogram shows the fraction of 1 Mpc³ volumes within the high-resolution region which use a given pre-reionization output. Specifically, z = (8.3, 10.2, 12, 1, 14). Sub-regions with z_eff ∼ 17 use the z_eff = 14 pre-reionization output. The dashed vertical line shows the z_init = 10.2 for our second-order initial conditions.
Download figure:
Standard image High-resolution image

To account for the different rates of structure formation, we take additional steps when generating our high-resolution region. (1) First, for each 1 Mpc³ cube within our coarse resolution volume, we calculate local density at z = z_init. (2) From the density, we use Equation (B1) to find the shift in effective redshift due the overdensity or underdensity of each subvolume, and then calculate z_eff (Figure 18). (3) Since we have a discrete set of pre-reionization outputs at z = (8.3, 10.2, 12.1, 14), we divide the 1 Mpc³ cubes within our high-resolution region into four bins based on their densities and effective redshifts. The fraction of cubic Mpc within our high-resolution region in each effective redshift bin is shown as the black histogram overlaid on the z_eff distribution in Figure 18. Note, that both the histogram and the smoother curve follow the same general shape. (4) Finally, based on which bin each Mpc³ falls into, we assign it a pre-reionization output.

To account for the faster evolution in our high density regions, for run D we use a z_init = 10.2 instead of the z_init = 8.3 used in runs A–C. This allows us to assign a z_eff = 8.3 to place near Milky Ways in our high density region. Nothing else changes substantially, except for using the z = 10.2 output from the low-resolution simulation to generate the l > 1 Mpc structure in the high-resolution region.

All the pre-reionization outputs are truncated to the same mass resolution using the same method described in the previous section. As near as possible, we use tracer particles of the same mass in each pre-reionization output. The truncated mass functions of the three additional pre-reionization outputs used in our second-order initial conditions are the additional panels (bottom and top right) in Figure 17.

After the high resolution has been built and the large-scale modes added at z_init = 10.2, we embed it inside the coarse resolution volume region and run to the present with Gadget 2. Run D has two key differences when compared to runs A–C. First, the 50³ Mpc³ snapshot used to generate the l > 1 Mpc modes is z = 10.2 instead of z = 8.3. Second, the z = 8.3 pre-reionization output is not used for the entire high-resolution volume, but only in the overdense regions, with outputs from z = (10.2, 12.1, 14) used for the average and underdense regions.

WHERE ARE THE FOSSILS OF THE FIRST GALAXIES? I. LOCAL VOLUME MAPS AND PROPERTIES OF THE UNDETECTED DWARFS

Article metrics

Permissions

Share this article

Dates

ABSTRACT

1. INTRODUCTION

2. NUMERICAL METHOD

2.1. Approximating Cosmic Variance

2.2. A Note on the Halo Occupation Distribution

3. TESTS OF THE METHOD

3.1. Mass Resolution

3.2. Softening Length

3.3. Subhalo Scale Comparisons

4. RESULTS

4.1. Definition of a Fossil Dwarf

4.2. Luminosity Threshold for Fossils

4.3. Fossil Properties

4.4. The Inner Ultra-faints

5. OBSERVATIONAL TESTS

6. SUMMARY AND CONCLUSIONS

APPENDIX A: FIRST-ORDER INITIAL CONDITIONS

APPENDIX B: SECOND-ORDER INITIAL CONDITIONS