THE EFFECTS OF PATCHY REIONIZATION ON SATELLITE GALAXIES OF THE MILKY WAY

Our modeling is based on the Aquarius⁴ simulations (Springel et al. 2008), a suite of six highly resolved Milky Way-sized dark matter halos. These are re-simulations of six halos from a lower resolution version (900³ particles) of the Millennium-II simulation, which is a periodic box of size 100 h⁻¹ Mpc containing 2160³ particles (Boylan-Kolchin et al. 2009). The six halos were selected on the basis of final mass, and by not having a massive close neighbor at z = 0 (no late-time major merger); otherwise they were selected randomly. We note that one of these halos (the F halo) is somewhat of an outlier in terms of its merger history, in that it experienced its last major merger at z = 0.7.

The same cosmological parameters were adopted as for the original Millennium simulation (Springel et al. 2005): Ω_m = 0.25, Ω_Λ = 0.75, σ₈ = 0.9, n_s = 1 and Hubble constant H₀ = 100 h km s⁻¹ Mpc⁻¹ = 73 km s⁻¹ Mpc⁻¹. While these are consistent within the uncertainties with the parameters estimated from the three-year Wilkinson Microwave Anisotropy Probe (WMAP) data (Spergel et al. 2007), the value of σ₈ = 0.9 is slightly higher than that inferred from the most recent measurements (Komatsu et al. 2011). In principle, the adopted σ₈ could lead to an overestimate of the number of halos collapsing around the reionization epoch, but studies like Boylan-Kolchin et al. (2011a) suggest this would have a relatively minor effect.

All six halos (labeled "A" to "F") were simulated with at least two different resolutions; we use the higher of the two ("level 2;" softening length of 65 pc) in our analysis. Table 1 lists the particle mass, final halo mass, and virial radius of each of the halos. Mass and radius are listed as M₂₀₀ and r₂₀₀, defined as the mass enclosed in a sphere with mean density 200 times the critical density, and the corresponding virial radius.

Snapshots are stored at 128 output times, from redshift 127 to 0, equally spaced in log a where a = 1/(1 + z). For each snapshot, subhalos are identified using the SUBFIND algorithm (Springel et al. 2001), and linked between snapshots in a merger tree. Since SUBFIND requires a minimum of 20 gravitationally bound particles to define a subhalo, the smallest halo masses we resolve are of order a few times 10⁵ solar masses. The merger trees are our starting point for identifying subhalos in the simulation that may host present-day luminous satellite galaxies.

The Aquarius halos are drawn from a larger parent simulation, so we can use the parent box to calculate how regions around these halos reionize depending on their environment. We use the semi-analytic approximation of Zahn et al. (2007), based on the excursion set treatment of Furlanetto et al. (2004a, 2004b). This semi-analytic model assumes that galaxies reside in halos above some minimum mass, M_min, and that the number of ionizing photons produced by a galaxy is directly proportional to its host halo mass. A given region is considered ionized when the collapse fraction—i.e., the fraction of matter that lies in halos above the minimum mass—in the region exceeds some threshold value. This threshold depends on the efficiency of the ionizing sources and is smaller for models with more efficient sources. For the purpose of modeling reionization, we fix the minimum mass at M_min = 10⁸ M_☉, and consider several models for the source efficiencies, spanning a range of possibilities for the timing and duration of reionization.

We take the initial conditions of the Millennium-II simulation as input to generate reionization maps for the cube of 100 h⁻¹ Mpc per side. Reionization redshifts are tagged spatially, at a resolution of 512 cells per side, or about 195 h⁻¹ kpc. As an example, a slice through one of the reionization cubes is shown in Figure 1, illustrating the patchiness of the process.⁵

Zoom In Zoom Out Reset image size

We calculate four different models, including a very efficient reionization that completes around redshift 12, and three gradually more extended models. The distributions of reionization redshifts in each model, over the cells in the box, are summarized in Table 2 and Figure 2. In the table, 〈z_reion〉 and σ_z refer to the mean and standard deviation of the reionization redshifts across the cells in the box, and z_end is the redshift where the process is complete. The quantity τ_e is the Thomson optical depth; the seven-year WMAP value for this is τ_e = 0.088 ± 0.015 (Komatsu et al. 2011).

Zoom In Zoom Out Reset image size

With these maps, we can then find the reionization redshift of any given subhalo as follows. At each snapshot, we locate the position of the subhalo in the parent box, and look up the redshift at which this region of the parent box reionizes. The first snapshot where this redshift is earlier than the redshift of the snapshot, flags the subhalo as having entered a reionized region, and we define this as the reionization redshift of the subhalo. Note that this allows not only for different reionization redshifts for the six main Aquarius halos, but also that their subhalos ending up as satellites are not required to reionize at the same time as the host.

We use the level-2 merger trees of the six Aquarius halos to determine which dark matter halos can host luminous satellite galaxies based on several criteria. These again reflect our assumptions about when a subhalo is able to accrete gas from the IGM, cool its gas, and form stars, following methods developed in, e.g., Madau et al. (2008), Koposov et al. (2009), and Busha et al. (2010). We only consider subhalos that end up within the virial radius of their host halos (see Table 1).

First, we require that subhalos must grow larger than a minimum v_max of 15 km s⁻¹ before star formation can take place, and also that this happens before the subhalo is reionized. This corresponds roughly to the virial temperature where atomic hydrogen is no longer able to cool (∼10⁴ K). We do not explicitly include H₂-cooling "minihalos" as sites for low-mass star formation before reionization, but recognize that such halos may have been the sites of the very first star formation, whose metals seeded the next generation of lower-mass, long-lived stars (Bromm et al. 1999, 2003, 2009; Yoshida et al. 2004).

Once a subhalo enters a reionized region (as described in Section 2.2), the minimum v_max needed to sustain star formation increases, reflecting both the increased temperature of the IGM to a few times 10⁴ K due to photoionization, and the possibility of photoevaporation of gas out of small halos (Barkana & Loeb 1999; Iliev et al. 2005). Gnedin (2000) found the filtering mass to correspond to the scale at which the gas fraction in halos is significantly reduced due to reionization. Gnedin's expression, however, does not take into account that only the cold gas remains available for star formation; Muñoz et al. (2009) and Busha et al. (2010) argue that only halos that grow to a viral temperature ∼10⁵ K will be able to sustain star formation after reionization. Accordingly, we only let star formation continue post-reionization in subhalos whose v_max is greater than 50 km s⁻¹. This represents a rather abrupt reionization effect; we also run a model where this suppression is weaker, with the threshold set at 30 km s⁻¹.

In addition, we assume that once a satellite starts interacting with the main halo, its gas is stripped and no further star formation takes place. Thus, there are two important redshifts associated with each satellite: z_re and z_infall, the times when it is reionized and falls into the main halo, respectively.

In order to account for tidal stripping of satellites, which has been proposed as another important mechanism to resolve the missing satellite problem (Kravtsov et al. 2004), we tag the 1% most bound particles of the subhalos as star-carrying when they are accreted onto the main halo. We then track all these tagged particles to redshift zero, and determine the final luminosity of a satellite based on its remaining tagged particles. In this way, we are able to take into account tidal stripping of satellites and other dynamical effects, under the assumption that the stars are found in the bottom of the potential well.

Of course, a number of infalling systems will be completely disrupted, and contribute to the accreted stellar halo. In Figure 3, we show the projected distribution of all the dark matter particles that were tagged as stars on infall, in reionization scenario 1. Surviving satellites here appear as bright, concentrated dots, while the particles from the shredded systems make up the diffuse halo. Since the accreted stellar halos and streams forming in the Aquarius simulations have been studied extensively in Cooper et al. (2010) and Helmi et al. (2011), we will not discuss them further here. We note, however, that although our model is very simple, it still reproduces most of the features seen in Figure 6 of Cooper et al. (2010)—we take this as justification that our simple model still reasonably captures the relevant processes.

Zoom In Zoom Out Reset image size

We determine satellite candidates for all four of our reionization calculations. In addition, we do the same for four "instantaneous" redshift models, where we set the reionization redshift equal to the mean reionization redshift in the entire box. Alvarez et al. (2009) showed that the distribution of reionization redshifts for Milky Way-sized halos is nearly similar to the global spatial distribution, making the global mean an appropriate choice. Comparing the resulting satellite populations from the patchy and instantaneous reionization models then allows us to isolate any effects of reionization not being spatially uniform.

In order to compare our satellite candidates to the observed Milky Way population, we need to assign luminosities to the subhalos. There are many different procedures for doing this; see, e.g., Koposov et al. (2009) for a discussion. Here, we adopt a simple model where the star formation rate (SFR) is proportional to the mass, that is M_star = ∫ × f_gas × M_DM dt, where = 1/τ is a constant that sets the star formation efficiency (i.e., τ can be considered a timescale for which a subhalo would convert all of its gas into stars). We take f_gas × = 0.08 × 10⁻¹⁰ yr⁻¹, tuned both so that we get a reasonable satellite population, but also so that the mass of the accreted stellar halo (i.e., the total mass of satellites that fall in and are completely shredded) falls in a mass range around ∼10⁹ M_☉ (Bell et al. 2008; Tumlinson 2010). This second criterion is insensitive to the reionization prescription, because the mass of the stellar halos is dominated by a few significant progenitors, while the smaller systems affected by reionization contribute a small fraction of the total mass. While in principle we could increase this product to get a better fit for the more extreme reionization scenarios, this would then lead to very massive stellar halos.

Our model is deliberately simple because we wish to assess and compare the impact of different reionization assumptions without getting addressing uncertainties in the baryonic physics. While we do not explicitly model effects like the impact of supernova feedback on the satellites, we assume that the typical effect is captured by setting (i.e., the star formation efficiency) appropriately. (For an example of a more extensive approach in modeling the baryonic content of dark halos in cosmological simulations, we refer the reader to Font et al. 2011.)

Having determined the stellar mass formed, we then assign each stellar population an associated present-day luminosity, using the stellar population synthesis (SPS) models of Bruzual & Charlot (2003), assuming a Salpeter initial mass function, and a metallicity of Z = 0.0001 (the lowest considered in their model, and which captures the typical range of the ultra-faint satellites we are most interested in, e.g., Kirby et al. 2008; Simon et al. 2010; Norris et al. 2010). We then track the particles tagged with stellar populations to redshift zero, and so determine the present-day luminosities of the satellite candidates picked out as described earlier.

As an alternative to the simple parametric method described above, we also try the abundance matching technique developed by Busha et al. (2010) to assign luminosities to dark matter halos. This is an extrapolation of the fit of Blanton et al. (2005) based on galaxies in the SDSS (down to M_r = −12.375) to fainter magnitudes, by a power law:

$$\begin{equation} M_V - 5 \log h = 18.2 - 2.5 \log \left[ \left(\frac{v_{{\rm max}}}{1 \,\mbox{km \,s}^{-1}} \right) ^{7.1} \right]. \end{equation} \tag{ 1 }$$

The appeal of the abundance matching procedure is that there are no parameters to tune, using a model derived from statistics of observed galaxies. While it is an extrapolation at UFD luminosities, Busha et al. (2010) showed that when applied to the Via Lactea subhalos, the abundance matching model did well in reproducing the observed luminosity function (but see Boylan-Kolchin et al. 2011b for a discussion of potential problems with abundance matching for the Milky Way dwarfs). For purposes of comparison, then, we also separately assign magnitudes to our subhalos using Equation (1). In order to select an appropriate v_max for the calculation, we use the peak v_max of the subhalos during the time they were able to form stars, as defined in Section 3.1. For most subhalos, this either corresponds to the time of reionization or time of accretion.

As some basic tests of our model, we check how well it reproduces various features of the observed Milky Way dwarf population, such as their radial distribution, and observed mass-to-light ratios. Figure 4 shows the radial distribution of satellites in our simulation, as a fraction of total number. We find all six halos to have a satellite distribution consistent with what is observed around the Milky Way, with the possible exception of the B halo, which has a larger fraction of satellites somewhat closer in.

Zoom In Zoom Out Reset image size

As for mass-to-light ratios, we note that subhalo masses in a simulation cannot be directly compared to observed dwarf galaxy masses. However, Walker et al. (2009, 2010) and Wolf et al. (2010) show that the mass within the half-light radius is well constrained for the Milky Way dwarfs. We thus estimate the half-light radius of our satellites from their assigned luminosities by fitting a power law to the data tabulated in Walker et al. (2010), and measure the mass within the half-light radius by fitting an NFW profile (Navarro et al. 1997). The result is shown in Figure 5: the red points are simulated satellites, and the black crosses (UFDs) and triangles (classical dSphs) show the data from the Milky Way dwarf galaxies. Our satellite properties are in good agreement with the observed values, and recover the observed relation between M and r_0.5. The points shown are for satellites with the Model 1 reionization history; driving reionization to earlier times yield the same relation, but fewer satellites overall, and fewer in the intermediate r_0.5 (or luminosity) range. We conclude that our simple model is overall sufficient for reproducing the observed trends.

Zoom In Zoom Out Reset image size

A major novel factor in our work is that the reionization redshift of a halo is completely determined by its environment. This has two main consequences: first, that the six main halos will reionize at different times according to their locations in the Millennium-II box. But additionally, the subhalos of each of these main halos are not required to reionize at the same time as their main halo—their reionization times are also set by their particular path through the main box. Here, we explore the distribution in reionization redshifts for the subhalos, and how this influences the resulting number of potential satellites.

Table 3 summarizes the first effect—that the six Aquarius host halos do not reionize at the same times in the box. In particular, the F halo consistently reionizes later than the box mean, while the C halo reionizes significantly earlier. The other four halos are closer to the mean of the box (and thus slightly earlier than the median; also see Table 2).

As for the second effect: Figure 6 shows the distribution of satellite reionization redshifts for the six halos, here in the case of scenario 1. (The distributions for the other three models show similar characteristics, but with fewer satellites and narrowing distributions as the models become more efficient.) Several interesting trends are evident: first, there is a distribution—while some of them are strongly peaked (C, D, and F in particular), all show a range of redshifts. Second, the peak of the distribution does not necessarily correspond to the reionization redshift of their host halo (the dotted lines).

Zoom In Zoom Out Reset image size

The C halo is the most extreme in this regard—the peak of the satellite distribution is near the mean of the box, while the host reionized significantly earlier. Essentially, this happens because the location of the C halo is near the edge of one of these early-reionized regions in the large box (cyan in Figure 1), responsible for the early reionization of the C halo itself. Many of its satellites, however, fall in from a different direction than the bubble front, and end up with a reionization distribution peaking around the box mean. Note that Figure 6 shows only systems that survive as self-bound satellites. If we were to also include systems that are shredded and are part of the diffuse main halo by redshift zero, we would see a larger number of subhalos reionizing around the same time as the C main halo.

Having established that there is indeed a distribution of satellite reionization redshifts, the next question is whether this leads to a different satellite population between patchy and instantaneous models. Figure 7 considers the number of satellites produced in each model for different sets of assumptions: either patchy or instantaneous, and a suppression v_max of either 30 or 50 km s⁻¹. We see that regardless of cutoffs and patchy versus instantaneous, as reionization is pushed to earlier times, there is a strong decline in the number of satellites. This overall trend is similar to what was found in previous studies (Muñoz et al. 2009; Tumlinson 2010; Busha et al. 2010). For a given reionization model, we do see additional differences, however. For either choice of suppression v_max, the instantaneous reionization models systematically make 10%–20% more satellites than the corresponding patchy models; in the more extended models this effect dominates the difference between stronger and weaker suppression. As the reionization sources become more efficient and reionization occurs earlier, the patchiness is less important, and the variation in suppression v_max becomes the dominant effect in this scenario.

Zoom In Zoom Out Reset image size

Note that the numbers shown in Figure 7 are averages; as shown in Figure 6 there is a substantial spread between the actual number of satellites between the different halos. The error bar indicates the typical rms, ± ∼ 15 satellites. However, halos A–E all show more satellite candidates in the instantaneous scenarios, while the F halo has about 10% fewer. For five out of the six halos, then, the main effect of patchy reionization is that the number of satellite candidates is reduced by 10%–20% compared to the instantaneous models. This can be understood in the context of Figure 6, where the F halo is the only halo whose majority of satellites reionize later than the mean.

In practice, the difference between the patchy and instantaneous case in Figure 7 is due to a combination of two effects: that the global mean is not necessarily the best instantaneous redshift to choose for a given halo, as well as not capturing the spread in reionization times. However, we do not find that using, e.g., the host halo reionization redshift, or the median satellite reionization redshift, leads to better agreement in each case. While we could recover the number of satellites in each case by tuning the instantaneous reionization redshift, the value of such an exercise is limited: the only way to determine that appropriate redshift would be to already know the result from the patchy model. As such, this illustrates the main problem with using instantaneous and/or spatially homogeneous reionization models. Depending on the application, the fact that there is a distribution of satellite reionization redshifts may be important to the results—and this distribution cannot be captured in a homogeneous model, nor can the appropriate instantaneous approximation be recovered without knowing the patchy result.

Having explored the reionization distributions in our simulations, we here discuss the resulting luminosity function of our simulated satellites. In order to compare them to the population of observed dwarfs, however, we first need to correct the observational sample of dwarf galaxies (data taken from Martin et al. 2008; Wadepuhl & Springel 2011) for incompleteness and biases to establish a luminosity function. We follow the method outlined in Koposov et al. (2008) and Tollerud et al. (2008), which employed models for the completeness of the SDSS Data Release 5 (DR5) to quantify the detection efficiency and thus adjust the observed luminosity function. In addition, we only include the simulated satellites that would be detectable according to the DR5 criteria:

$$\begin{equation} r_{{\rm max}} = \left(\frac{3}{4\pi f_{{\rm DR5}}} \right)^{1/3} 10^{\left(-0.6 M_V -5.23\right)/3} \mbox{ Mpc}, \end{equation} \tag{ 2 }$$

where f_DR5 = 0.194 is the fraction of the sky covered by DR5. That is, simulated satellites with r > r_max(M_V) are not included in the plots in this section. In each case, r_max is calculated with respect to a point 8 kpc from the center of the main halo.

Figure 8 shows the resulting cumulative luminosity function for the simulated satellite population of all six halos, for the four patchy reionization models and with v_max = 50 km s⁻¹. The luminosities shown are those obtained with the abundance matching technique (Equation (1)), but the resulting luminosity functions are very similar to those using our simple SPS-based model. (However, there is substantial scatter between the methods on a subhalo–subhalo basis, i.e., rms ∼1 mag.) The dotted line shows the resulting luminosity function when applying the abundance matching prescription to all subhalos in each simulation that lie within the virial radius by redshift z = 0, without assuming any suppression due to baryonic processes. As such, it gives a good illustration of the missing satellite problem: without some process such as reionization to suppress star formation and reduce the number of luminous satellites, the six Aquarius halos would predict between 500-1000 satellites around Milky Way-like galaxies.

Zoom In Zoom Out Reset image size

We see that while satellites brighter than M_V ∼ −10 are insensitive to the reionization model, the total number and the number of faint satellites change dramatically with more efficient reionization. This is consistent with the results of Font et al. (2011), who found that photoheating is mainly important for suppressing star formation in satellites analogous to ultra-faint dwarfs. In our simulations, only Model 1, and in some cases Model 2, can match the faint end of the luminosity function. In fact, there is a factor of 3–4 difference between the total number of satellites predicted between the different reionization models for a given halo. This behavior may ultimately allow for constraining the Milky Way reionization epoch with the faintest satellites.

Another striking feature of Figure 8 is the scatter between the six halos—there is a factor of ∼3 difference between the B and D halos, which have the fewest and most satellites, respectively. Thus, the halo–halo scatter is of the same order as the differences in the impact of the various reionization models on the faint end. In order to break this degeneracy, it will be necessary to understand the general distribution of substructure around Milky Way-like halos, requiring a statistical sample of halos at sufficient resolution. Since the bright end of the luminosity function is not affected by the reionization history, however, it may be possible to break the degeneracy by constraining the halo–halo scatter using the brighter satellites, and use the fainter satellites to constrain the reionization epoch of the Milky Way.

We find that only the F halo hosts a satellite as bright as the Magellanic Clouds. This is consistent, however, with recent studies based on SDSS DR7 on the abundance of close, bright satellites around Milky Way-like galaxies (Liu et al. 2011; Boylan-Kolchin et al. 2011a; Busha et al. 2011; Guo et al. 2011; Tollerud et al. 2011), and the possibility that the Clouds have entered the Milky Way halo only recently (Besla et al. 2007, 2010). Hence, we do not consider the mismatch at the brightest magnitudes to be problematic, especially given the simplicity of our star formation model.

Similar to the brighter satellites, the total stellar masses in the accreted halo of shredded, infalling satellites do not change significantly with reionization model. This is consistent with the results of Cooper et al. (2010), who studied in detail the building of accreted stellar halos in the Aquarius simulations. They found that the significant progenitors are similar in mass to the brightest Milky Way dwarf spheroidals, with each halo assembled from less than five such objects.

To highlight some of the differences between the various model assumptions, Figure 9 shows the luminosity function of the D halo satellites for some parameter choices. The upper panels show the results using patchy reionization and v_max = 50 km s⁻¹, instantaneous reionization and v_max = 50 km s⁻¹, and patchy reionization and v_max = 30 km s⁻¹ respectively. The three lower panels show the same, but with luminosities determined by abundance matching (Equation (1)).

Zoom In Zoom Out Reset image size

While the models with post-reionization star formation suppressed at v_max = 50 km s⁻¹ give a good match to the observed luminosity function, when we decrease the suppression threshold to 30 km s⁻¹ none of the models fit the bias-corrected luminosity function well, despite the overall number of satellites only being slightly increased. This holds true both for the abundance matching-derived luminosities, and the best-fit SFR/SPS model. In particular, even if the star formation efficiency is decreased to give the correct number of satellites overall, the shape still does not match the observed luminosity function. It is possible, however, that this lower suppression model could still be a good match, if we explicitly consider other processes like supernova feedback (Li et al. 2010).

Compared to patchy reionization, the instantaneous models generally yield more satellites. While only the D halo is shown here as an illustration, as discussed in Section 4.1 this holds true for all of the Aquarius halos, except for the F halo. It is interesting to note, however, that the differences between different reionization models, luminosity prescriptions, etc., are in general dominated by the halo–halo scatter. That is, the differences in the satellite populations of the six halos for any given model are as large as the differences between the range of models we are exploring. As an example, Figure 10 shows the satellite luminosity functions for all six halos overplotted, here for reionization model 1 and the SFR luminosity prescription. While the B halo, which is the clear outlier, is also the lowest mass halo in this simulation suite, the uncertainty of where the Milky Way itself fits in this distribution makes direct comparisons to constrain models difficult.

Zoom In Zoom Out Reset image size