GROUP FINDING IN THE STELLAR HALO USING PHOTOMETRIC SURVEYS: CURRENT SENSITIVITY AND FUTURE PROSPECTS

We use a catalog of 59392 M-giants identified from the 2MASS survey in a companion paper (Sharma et al. 2010; Paper I). A brief description of the catalog and the selection procedure is repeated below. The 2MASS all sky point source catalog contains about 471 million objects (the majority of which are stars) with precise astrometric positions on the sky and photometry in three bands J, H, K_s. The survey catalog is 99% complete for K_s < 14.3. An initial sample of candidate M-giants was generated by applying the selection criteria:

$$\begin{eqnarray} 10.0 & < & K_s < 14.0 \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} &&J-K_s > 0.97 \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} J-H & < & 0.561(J-K_s)+0.36 \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} &&J-H > 0.561(J-K_s)+0.19 \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} &&K_s{\rm sin}(b) > 14.0\, {\rm sin}(15^{\circ }). \end{eqnarray} \tag{ 5 }$$

All magnitudes in the above equations are in the intrinsic, dereddened 2MASS system (labeled with subscript 0 hereafter), with dereddening applied using the Schlegel et al. (1998) extinction maps. These selection criteria and the dereddening method are similar to those used by Majewski et al. (2003) to identify the tidal tails of Sagittarius dwarf galaxy.

Distance estimates for these stars were made by assuming a uniform metallicity of [Fe/H] = −1.0 and an age of 13.1 Gyr and calculating a linear fit to the color–magnitude relation of the giant stars using a theoretical isochrone from the Padova database (Bertelli et al. 1994; Marigo et al. 2008; Bonatto et al. 2004). The best-fit relationship is M_K = 3.26 − 9.42(J − K_s)₀ where M_K is the absolute magnitude of the star and (J − K_s)₀ its color. Using this relation the distance modulus is then given by

$$\begin{eqnarray} &&\mu =(K_s)_0-3.26+9.42(J-K_s)_0. \end{eqnarray} \tag{ 6 }$$

The specific values of age and metallicity were adopted because (1) they roughly correspond to the expected values for the stellar halo and (2) they lead to the identification of the largest known stellar halo structure in the data, the tails of the Sagittarius dwarf galaxy, with maximum clarity.

To make theoretical predictions of structures in the stellar halo, we use the 11 stellar halo models of Bullock & Johnston (2005), which were simulated within the context of the ΛCDM cosmological paradigm. These simulations follow the accretion of individual satellites modeled as N-body particle systems onto a galaxy whose disk, bulge, and halo potential is represented by time dependent analytical functions. Semi-analytical prescriptions are used to assign a star formation history to each satellite and a leaky accreting box chemical enrichment model is used to calculate the metallicity as a function of age for the stellar populations (Robertson et al. 2005; Font et al. 2006). The three main model parameters of an accreting satellite are the time since accretion, t_acc, its luminosity, L_sat, and the circularity of its orbit, defined as = J/J_circ (J being the angular momentum of the orbit and J_circ the angular momentum of a circular orbit having same energy). The distribution of these three parameters describes the accretion history of a halo. To study the sensitivity of the properties of structures in the stellar halo to accretion history, we additionally use a set of six artificial stellar halo models (referred to as non-ΛCDM halos) from Johnston et al. (2008) that have accretion events that are predominantly (1) radial ( < 0.2), (2) circular (>0.7), (3) old (t_acc > 11 Gyr), (4) young (t_acc < 8 Gyr), (5) high luminosity (L_sat > 10⁷ L_☉), and (6) low luminosity (L_sat < 10⁷ L_☉).

2.2.1. Generation of Synthetic Surveys from Simulations

Each N-body particle in the simulations represents a population of stars having a certain stellar mass, a distribution of ages and a monotonic age–metallicity relation. The real data from surveys on the other hand are in the form of color–magnitude combinations of individual stars, shaped by observational or other selection criteria. The procedure to convert a simulated stellar halo model to a synthetic survey consists of the following four steps⁶.

• 1.  
  Spawning the stars. We use the Chabrier (2001) exponential initial mass function (IMF) to generate stars corresponding to the given stellar mass of an N-body particle. If M is the stellar mass associated with the N-body particle and m_mean the mean mass of the IMF distribution then the number of generated stars is given by approximately M/m_mean. We use stochastic rounding to convert M/m_mean to an integer—i.e., if the fractional part of M/m_mean is less than a Poisson distributed random number with a range between 0 and 1, we increment the integral part by 1. Each generated star is assigned a mass that is randomly drawn from the IMF, an age that is randomly drawn from the given age distribution and metallicity computed from the age–metallicity relation.
• 2.  
  Assigning color and magnitude. Next, a finely spaced grid of isochrones obtained from the Padova database (Bertelli et al. 1994; Marigo et al. 2008; Bonatto et al. 2004; Girardi et al. 2004) is interpolated to calculate the color and magnitude of each generated star. The absolute magnitude is converted to an apparent magnitude based on the distance of the N-body particle from a given point of observation. In this paper we assume the point of observation to be at 8.5 kpc from the Galactic Center in the plane of the Galactic disk. Finally, depending upon the color–magnitude limits of the survey being simulated, the stars are either accepted or rejected.
• 3.  
  Distributing the stars in space. The position coordinates of first star that is spawned by an N-body particle, is assumed to be same as that of the spawning particle. Each subsequent star spawned by the N-body particle is distributed in a spherical region around the particle. The radius R of the spherical region is determined by the N-body particle's 32nd nearest neighbor in (X, Y, Z) position space. For this the multi-dimensional density estimation code EnBiD (Sharma & Steinmetz 2006) is used. The stars are radially distributed according to the Epanechnikov kernel function (1 − r'²)r'² where r' = r/R and r' varies from 0 to 1. A surface density map in the X–Z plane of one of the ΛCDM halos sampled using the above scheme is shown in Figure 1. It can be seen that the structures are quite adequately reproduced.
• 4.  
  Observing our stars. In simulations the distance of the stars is exactly known but in real surveys the distance is estimated by using observed properties of the selected stars. We use the same procedure to estimate the observed distances for the stars in our synthetic surveys.

Zoom In Zoom Out Reset image size

2.2.2. Synthetic Surveys and Data Sets

In this paper we use the hierarchical group finder EnLink that can cluster a set of data points defined over an arbitrary space (described in detail in Sharma & Johnston 2009). Paper I contains a complete discussion of the optimum parameters for the group finder, as well as the choice of coordinate system for the data. We include a brief outline of each of these below.

EnLink identifies the peaks in the density distribution in the data space as group centers. The region around each peak, which is bounded by an isodensity contour corresponding to the density at the valley, is associated with the group. The number of neighbors employed for this density estimation was chosen to be k = 30: a smaller value makes the results of the clustering algorithm sensitive to noise in the data, while a larger value means that small structures go undetected.

EnLink is built around a metric that can automatically adapt to a given data set, and hence can be applied to spaces of arbitrary numbers and types of dimension. However, it was shown in Paper I that the most effective way to deal with photometric data sets where the angular position of stars is well known, but the distance can be very uncertain is to transform to a Cartesian coordinate system and use a simple Euclidean metric. Prior to the transformation, a new distance estimate r' = 5(log(r/(10 pc))) − μ₀ is adopted, where μ₀ is a constant that determines the degree to which the radial dimension is ignored or used. The optimum value of μ₀ was found to be 8. Note that the clustering results were not overly sensitive to the exact choice of μ₀.

A group finding scheme in general also generates false groups due to Poisson noise in the data. In order to make meaningful comparison between different data sets with different number of data points it is imperative that one uses a consistent and objective scheme to screen out spurious groups. In EnLink spurious groups that can arise due to Poisson noise in the data are screened out by looking at the significance

$$\begin{eqnarray} &&S=\frac{\ln (\rho _{\rm max})-\ln (\rho _{\rm min})}{\sigma _{\ln \rho }}. \end{eqnarray} \tag{ 10 }$$

The contrast between the density at the peak (ρ_max) and valley (ρ_min) where it overlaps with another group can be thought of as the signal of a group and the noise in this signal is given by the variance σ_ln(ρ) associated with the density estimator. An advantage of such a definition is that the resulting distribution of significance of spurious groups is almost a Gaussian and hence the cumulative number of groups as function of significance can be written as

$$\begin{eqnarray} &&G(>S)=\left(1-{\rm erf}\left(\frac{S f_1}{\sqrt{2}}\right)\right)\frac{f_2 N}{k}, \end{eqnarray} \tag{ 11 }$$

where N is the number of data points, k is the number of neighbors used for density estimation, and f₁ and f₂ are constants for a given k and the dimensionality of data d. We adopt a threshold S = S_Th for each of our synthetic surveys (depending on the number N of stars in the survey) that satisfies the condition G(>S_Th) = 0.5. This means that the expected number of spurious groups in our analysis is 0.5. In EnLink for a chosen S_Th, all groups with S < S_Th are denied the status of a group and are merged with their respective parent groups.

A value of f₁ = 1.0 and f₂ = 15.5/(d^2.1k^0.2) were empirically derived by Sharma & Johnston (2009) and shown to be valid for a wide range of cases. An improved version of the formula with $f_1=\sqrt{(d/4.0)(1-2.3/k)}$ and f₂ = 0.4 was presented in Bland-Hawthorn et al. (2010). But these formulas were derived for the case where the data is generated by a homogeneous Poisson process and are reasonably accurate for most applications. For non-homogeneous data the presence of density gradients on large scales makes the formula slightly inaccurate. Hence, for greater precision one should derive the values of f₁ and f₂ for a model data whose large-scale distribution is similar to the data being studied but otherwise does not have any substructures in it. We consider two models here: (1) a Gaussian sphere and (2) a Hernquist sphere with a scale radius of 15 kpc, which was shown to be an appropriate description of the stellar halo on large scales by Bullock & Johnston (2005). The sample size of models was set to 10⁷. A value of f₁ = 0.93 and f₂ = 0.4 was found to fit the distribution of significance of spurious groups and this is shown in Figure 2. The distribution of significance averaged over 11 ΛCDM halos for the rest of the data sets is also shown. At small S the number of groups are dominated by the spurious groups and the curves lie on the predicted relationship. At large S due to presence of real structures the curves flatten out. The symbols in the figure denote the mean adopted value of S_Th for each data set. The distribution of stars in data set S5 is very similar to data set S3 except for the fact that it is not over the whole sky. Hence, we revise the value of S_Th for the data set S5 slightly from 4.55 as shown in the plot to 4.75. Note, for data sets with sample size less than 10⁵ the adopted relation is found to slightly underestimate the number of spurious groups but we ignore this. Even if this fact is taken into account this would only lead to a very minor revision in values of S_Th for data sets S2 and S4 and since the corresponding curves are quite flat in this region the number of detected groups would also not change much.

Zoom In Zoom Out Reset image size

An example application of the group finder to one of the halos from data set S1 is shown in Figure 1. The upper panel shows the surface density of stars in the X–Z plane while the lower panel shows the stars which have been identified as groups. It can be seen that our scheme is very successful in detecting the structures. Even faint structures that are not easily visible in the surface density maps are detected by the group finder.

EnLink is a hierarchical group finder and the reported groups obey a parent–child relationship, with one master group being at the top. In our case this master group represents the smooth component of the halo and the rest of the groups are considered as substructures lying within it. Note that since we analyze magnitude limited surveys we expect two density peak even in a smooth halo: one at the Galactic center and the other at the location of the Sun due to the presence of a large number of low-luminosity stars. Hence, we ignore all groups whose density peaks lie within 5 kpc of the Sun or the Galactic center.

When we apply EnLink to our synthetic data sets, because we know the progenitor satellites from which the stars originally came from, we can make some quantitative assessment of how well the group finder has performed. Following Paper I, we define purity as the maximum fraction of stars in a group that came from a single satellite. If I is the set of all satellites and J is the set of all groups, then

$$\begin{eqnarray} {\rm Purity}(j) & = & \frac{{\rm max}_{i \in I}\lbrace n_{ij}\rbrace }{n_{.j}}, \end{eqnarray} \tag{ 12 }$$

where n_ij is the number of stars from satellite i in group j and n_.j = ∑_in_ij is the total number of number of stars in that group. The mean value of purity P = ∑Purity(j)/|J|is a good indicator of to what extent the recovered groups correspond to real physical associations.

The peak density ρ_peak of an unbound structure evolves in time from an initial maximum state set by the structural properties of the satellite before disruption to a final minimum state where it is fully mixed in the entire volume roughly outlined by the orbit of the satellite. In the intermediate state, transient structures such as streams, shells, and clouds are apparent. Below, we individually analyze each simulated satellite debris (a total of 1515 satellites' debris from our 11 stellar halo models) in order to understand the dependence of peak density on the properties of its progenitor satellite during each of these stages.

Initially, the structure is bound and its peak density ρ_peak is observed to depend primarily upon the luminosity L of the progenitor satellite. The exact dependence assumed in the simulations is given by ρ_peak ∝ L^0.4 and comes from fitting King profiles to observed dwarf galaxy data (Mateo 1998). In the final, completely phase mixed state the ensemble of stars in the original satellite occupy a large volume V_mix and have peak density of order L/V_mix. In a spherical potential, V_mix is expected to be inversely dependent on the circularity of the orbit because the debris occupies the region between the apocenter and the pericenter of the orbit and an eccentric orbit (with lower ) explores a wider range in radius. In the case of an axisymmetric potential the debris also occupies the region generated by the precession of the orbit about the axis of the system. Evolution in the intermediate state is governed by phase mixing, which occurs on a characteristic time scale T_mix that depends on both the circularity and the mass of the satellite (Johnston 1998; Helmi & White 1999), which in turn is correlated to luminosity L. In this stage, the density has both spatial and temporal dependence. Using the action angle formalism, Helmi & White (1999) analytically deduced the general dependence as ρ_peak ∝ 1/(r²t²_acc) for spherical potentials and as ρ_peak ∝ 1/(r²t³_acc) for axisymmetric potentials, where r is the distance of the peak density of a satellite debris from the center of the parent galaxy and t_acc is the time since accretion of the progenitor satellite.

Rather than model these multiple stages and effects independently, we instead examine empirically to what extent the peak densities of debris in our simulations can be connected to progenitor properties with a single formula. Overall, the physical effects described above lead us to expect that ρ_peak follows a 1/r² radial dependence and is largest for objects that are luminous, recently accreted and are disrupting on mildly eccentric orbits. We find that the relation

$$\begin{eqnarray} \rho ^{\rm emp}_{\rm peak} & \propto & e^{-t_{\rm acc}}\frac{L \epsilon }{r^2} \end{eqnarray} \tag{ 13 }$$

provides the least scatter about the actual value ρ_peak measured for the simulated satellite debris (Figure 3) and is a reasonable representation of trends discussed earlier. Note that in Figure 3, instead of associating ρ_peak simply with the global density maximum, we find the point where ρr² is a maximum, a choice that effectively identifies the point where the density contrast with respect to the background is greatest.

Zoom In Zoom Out Reset image size

The top panels of Figure 4 shows the average accretion history adopted in 11 ΛCDM model halos by plotting the probability of occurrence of accretion events P(E) (i.e., the number fraction of events of a given (L, t_acc in the left-hand panel, and of a given (, t_acc) in the right-hand panel). Note that the distribution of is non-uniform for recent events because we only included those that gave rise to unbound structures—events on circular orbits ( ∼ 1) remain bound for a longer time than those on more eccentric orbits ( ∼ 0). For t < 8 Gyr the probability also abruptly goes to zero for L < 10⁵ L_☉. This is because recent accretion events with L < 10⁵ L_☉ are not included in the simulations.

Zoom In Zoom Out Reset image size

We can use our understanding of the evolution of the peak density of satellite debris to develop some intuition for which of the events in the top panels of Figure 4 are most likely to be recovered as groups. An unbound structure in a stellar halo is expected to be detected as a group if the ratio ρ_peak/ρ_b > 1 (where ρ_b is the density of the background stars, i.e., halo stars after excluding the stream under consideration). Using Equation (10) this corresponds to a significance level of S ∼ 3. The radial distribution of halo stars was explored in Bullock & Johnston (2005) and they found that the slope, dln(ρ)/dln r, transitions between −1 for r < 10 kpc and −3.5 for r > 50 kpc. As a first order approximation we assume ρ_b ∝ 1/r², which when combined with Equation (13) gives for the ratio of peak to background density as

$$\begin{eqnarray} \rho _{\rm peak}/\rho _{b} & \propto & e^{-t_{\rm acc}} L \epsilon. \end{eqnarray} \tag{ 14 }$$

In the above equation, among the three parameters of an accretion event, the strongest dependence is on t_acc; so the oldest accretion events have the least probability of being detected. For given t_acc, events with high luminosity and on more circular orbits (high ) are more likely to be detected.

We confirmed this intuition by applying the group finder to the surveys if the 11 ΛCDM models our idealized data set S1—the second row of panels in Figure 4 shows the fraction of accretion events in the (L, t_acc) and (, t_acc) plane that were recovered as groups from these surveys. Formally, if P(E) is the probability of occurrence of all the accretion events and P(SE) the probability of an event being identified by the group finder in a stellar halo, then the panels plot the probability of detecting an event in a stellar halo, given by the conditional probability P(S|E) = P(SE)/P(E). A comparison of the second and top rows in the figure demonstrates that nearly all events on all types of orbits and of all luminosities are recovered for t_acc < 8 Gyr. Older events on more circular orbits are detectable as groups even further back, as might be anticipated from Equation (14). However, few accretion events older than 10 Gyr were recovered and none older than 11 Gyr. We conclude that the phase mixing of debris imposes this fundamental limit, and photometric surveys alone will never be able to explore this epoch of our Galaxy's accretion history.

Having confirmed how the recovery of a disrupted satellite as a group depends on its accretion characteristics, we now examine to what extent the distribution of group properties in a stellar halo reflects the distribution of properties of accreted objects, and hence what we might be able to say about the recent Galactic accretion history.

For example, simple intuition suggests that the number of groups should reflect the number of recent events. The top row of Table 2 shows that we recovered an average of 62 groups from our S1 ΛCDM data sets with a minimum of 36 and a maximum of 82. As anticipated, we recovered significantly more groups from those of our non-ΛCDM that had a larger number of recent accretion events (i.e., the ones biased toward recent accretion, or dominated by many low-luminosity events) and fewer groups from those that had a smaller number (i.e., the ones biased toward ancient accretion, or dominated by a few high-luminosity events). The non-ΛCDM stellar halo formed from predominantly circular events also gave rise to more events as these can be detected for longer (as demonstrated in Section 4.2 and by Equation (14))

Similarly, since accretion time is the dominant factor that determines the detectability of an accretion event, we expect the fraction of halo material detected as groups to be primarily related to the integrated mass accretion history of the halo and less sensitive to the total number of such accretion events. The top row of Table 3 shows that the fraction of stars in groups varies widely for our S1 ΛCDM data sets with a minimum of less than 2% and a maximum over 20%: the total fraction is heavily influenced by the most luminous object accreted on to the halo, and these are very few in number. The recently accreted halo (i.e., the young halo) shows the maximum fraction in groups as the structures have not yet dispersed. In spite of the total number of accretion events being very different, the low-luminosity and high-luminosity halos have similar fractions because for both of them a similar amount of mass was accreted within the last 10 Gyr. The circular halos, which can probe much older events, have higher fractions compared to the radial halos.

We would also like to know the luminosity function of accreted objects, which should be reflected in the distribution of group stellar masses m_group. The thin black lines in the top left panel of Figure 5 show the normalized, cumulative distribution of group masses recovered from our 11 ΛCDM—we found that our non-ΛCDM halos that were constructed from events that had a ΛCDM luminosity function (i.e., the old/recent and circular/radial halos) all lay within the region occupied by the ΛCDM lines (i.e., they had the same group mass function). In contrast, the brown line in this panel shows the group mass function for the high-luminosity non-ΛCDM halo, which is clearly distinct from the ΛCDM lines and biased toward higher mass groups. While the corresponding results for the low-luminosity halo are less striking in this panel (yellow line), Table 2 already demonstrates that the difference in luminosity function could be seen in the far larger number of groups recovered from this halo than for the ΛCDM halos.

Zoom In Zoom Out Reset image size

The last unknown in the accretion history is the orbit distribution. Since eccentric orbits with low have larger apocenters than more circular orbits of the same energy and a larger fraction of time during an orbit is spent near the apocenter rather than the pericenter, groups formed from radial accretion events are expected to be detected further away from the Galactic center. The probability distribution of radial distance of groups for data set S1 shown in top left panel of Figure 6 demonstrates this. The black lines are again for the ΛCDM halos. The radial (purple) and the circular (cyan) halo profiles clearly outline the black lines. More than tow-thirds of the groups in the radial halo are above 100 kpc while two-thirds of the circular halo groups are below 60 kpc. For the ΛCDM halos about 30% of the groups lie above 100 kpc.

Zoom In Zoom Out Reset image size

In summary, our analysis here using data set S1 suggests that, if the distance of the stars is known accurately, then number of groups, amount of material in substructures and the probability distributions of masses and radial distance can be related to the number, luminosity function and orbit distribution of recent accretion events. The remaining question is whether the accretion history of a halo can be constrained with real observational data.

Our results (top row Table 3) suggest that about 9% of material is in unbound recovered groups. As mentioned earlier, the fraction of material in groups is found to vary greatly from halo to halo and this is because it is typically dominated by a few luminous accretion events. A significant fraction of material is also contained in the bound groups. When these are taken into account the fraction is found to rise to 36% (last four columns of the top row in Table 3). Nevertheless, the fraction of material in unbound groups is still quite small. To explore this further, we plot in Figure 7 the fraction of material in recovered groups as a function of galactocentric radius (black solid line). The contributions of bound and unbound groups are also shown separately (blue and dark green solid lines). To create the plots the data points were tagged as recovered or un-recovered and bound or unbound and then binned by radial galactocentric distances. The dashed lines show the actual fraction of material in bound and unbound structures in the simulations. For bound groups the actual fraction is indistinguishable from the recovered fraction, which simply illustrates that all the bound groups are correctly recovered by our clustering scheme. The dashed lines, i.e., the actual fractions, show that outer halo (beyond 120 kpc) is predominantly bound whereas the inner halo (less than 60 kpc) is predominantly unbound. Next we look at the recovered fractions. The f^recovered_unbound(r) is quite low in the inner regions and then rises slowly and peaks at around 100 kpc. This rise coincides with a fall of fraction f^recovered_bound(r). This is due to the following two facts: (1) low (circularity of orbit) are more likely to be unbound than circular events since they pass close to the center and (2) low events are more likely to be found at the apocenter and their apocenter is further away than that of an high event of similar energy. The black line representing f^recovered_{bound+unbound}(r) shows that the outer halo (beyond 60 kpc) is highly structured whereas the inner halo (less than 40 kpc) which contains about 50% of the material (as shown by the red line) has very little material which can be recovered as structures. This implies that in the inner regions strong phase mixing greatly limits the amount of material that can be recovered as structures in three-dimensional configuration space. The analysis also suggests that for a given survey the fraction of material in groups will depend sensitively upon the contribution of the inner halo to the survey, which in turn is determined by the geometry and the depth of the survey.

Zoom In Zoom Out Reset image size

A few other issues which can affect the fraction of material in groups are as follows. We have here considered only the three-dimensional configuration space, additional information in the form of velocities should help detect more groups and increase the fraction of material in groups. Also, the choice of photometric selection function can decrease or increase the fraction of material in groups depending upon the contribution of the smooth component to the sample of stars in the survey (see Section 5.2 for further details). Issues related to our choice of the clustering scheme which can affect the fraction of material in groups are discussed in Section 6.2.

Figure 8 summarizes our results by plotting the number fraction of accretion events recovered as a function of accretion time, where each symbol represents a 1 Gyr interval. As expected, our idealized survey S1 recovers the largest fraction of events at all times, and the SDSS MSTO survey (S5), which covers less than 1% of the volume of any of the other surveys, recovers the fewest. However, there is no simple answer to the question of which survey is most sensitive. Our synthetic data sets represent surveys of different stellar tracers (i.e., M-giant, MSTO or RR Lyrae stars) selected with different observing strategies. As a result, these data sets explore the space around the Galaxy with a variety in the numbers of stars, depths and accuracies in distance estimates, as outlined in Section 2.2.1 and summarized in Table 1. In addition to these differences in spatial exploration, systematic differences in the stellar populations of objects of different luminosities and accretion times (e.g., as reflected in the stellar-mass/metallicity relation for Local Group dwarfs, see Larson 1974) means that the choice of tracers affects the relative number of stars contributing to each data set from different accretion events.

Zoom In Zoom Out Reset image size

Figure 9 illustrates this effect by plotting the sampling probability of accretion events in the L–t_acc plane relative to data set S1—the sampling probability being computed as the probability distribution of stars appearing in a survey in the L–t_acc plane (the L and t_acc of a star being that of its parent satellite). The top panel, for data set S2, clearly shows the expected bias of M-giants toward tracing the highest metallicity and hence highest luminosity events. In addition, M-giants are intermediate-age stars, which means that such surveys do not contain stars from ancient accretion events.

Zoom In Zoom Out Reset image size

The second two panels show the sampling probability for our deep MSTO surveys S3 and S3', which differ only in the color range from which stars are selected. Both surveys contain stars from lower luminosity objects and earlier accretion times than the S2 surveys. A comparison of the two shows that as the red edge of the color limit increases, the sampling probability increases for old and high-luminosity events (upper right-hand region in the plots) and decreasing the red edge of the color limit has the opposite effect, i.e., sampling probability increases for recent and low-luminosity events (lower left-hand region in the plots). This is because (1) the blue edge of the MSTO stars in an isochrone (knee-shaped feature in the color–magnitude diagram) shifts redward with the increase in age and metallicity of the stars and (2) the high-luminosity events are also metal rich. Increasing the photometric errors has an effect similar to increasing the color range.⁸ This is the reason why the sampling probability for high-luminosity events is slightly higher for S3 as compared to S1 although both surveys have the same color limits.

The fourth panel, for the RR Lyrae survey (data set S4), shows the strongest bias toward old and low-luminosity events as RR Lyraes are old, low-metallicity stars.

The bottom panel, for the MSTO samples with SDSS sky-coverage and magnitude limits (i.e., data set S5), indicates that stars from recent, low-luminosity events and ancient, high-luminosity events are poorly represented in this survey compared to the S1 sample. Recall that S1 surveyed all-sky to greater depth (r = 24.5) than S5, with the same color cut, but assuming perfect knowledge of distances. Figure 10 shows the sampling probability for surveys like S1 but with an SDSS-like latitude (top panel) and magnitude (middle panel) cut. Combining the two cuts (bottom) panel reproduces the bulk of the features in the lower panel of Figure 9 which demonstrates that accretion events are "missing" from the SDSS sample simply because of the small sample volume rather than errors in distance estimates.

Zoom In Zoom Out Reset image size

Overall, we conclude that the variety of observing strategies and tracers used in current and future surveys can lead to a range in the contrast that a given structure may have with respect to the smooth background. This suggests that different surveys will be sensitive to different types of accretion events.

The last five rows in Figure 4 show the fraction of events recovered as groups from our S2, S3, S3', S4, and S5 surveys of the ΛCDM stellar halo models, which can be broadly understood given the survey characteristics outlined in Table 1 and Section 5.1.

For example, relative to the other surveys, the S2 (M-giant) data sets contain a low number of stars, that explore only out to ∼100 kpc from the center of the Milky Way with modestly accurate distance estimates and a bias toward intermediate-age, high-metallicity populations. These effects combined explain the absence of groups in the S2 panels of Figure 4 detected from events that were old (M-giant stars are intermediate age), low luminosity (M-giant stars are relatively high metallicity), or on eccentric orbits (with apocenters beyond 100 kpc where debris spends most of its time). However, despite the relatively low number of stars overall, these surveys are particularly well suited to finding signatures of high luminosity, relatively recent accretion events: these stand out in M-giants at high surface brightness (because of the high-luminosity and recent accretion of the progenitor) and better contrast to the background (since ancient events which form a smooth background are not represented).

Data set S3 (MSTO stars) explores a similar volume to the S2 surveys, but contains a far greater number of stars, with less accurate distance estimates and a wider variety in the ages and metallicities of stars. Overall Figure 4 demonstrates that these MSTO surveys can recover additional groups from lower luminosity and more ancient accretion events (because the survey contains lower metallicity and older stars), as well as those on more eccentric orbits (because the large distance errors mean that a significant fraction of stars scatter into the sample from beyond 100 kpc).

Data set S3' is constructed the same way as S3, but with a finer color cut. This lowers the total number of stars, increases the accuracy of the distance estimates but decreases the contribution from older and higher metallicity stars. Since the last property decreases the contribution of the smooth background, the S3' surveys do even better than the S3 surveys at recovering groups corresponding to low-luminosity events.⁹

Data set S4—the RR Lyrae surveys—contains only a comparable number of stars to the M-giants surveys but explores to far greater depth and with much more accurate distances than S2, S3, or S3', and represent old, low-metallicity populations. Figure 9 shows that stars from intermediate-luminosity satellites are missing from the survey, and Figure 4 confirms that these events are not recovered as groups. However, S4 has a similar success as S3 and S3' in filling out other regions of accretion history space.

Finally, data set S5—the SDSS MSTO surveys—contains the same stellar populations as S3 but explores a much smaller volume because of the brighter magnitude limit. A comparison of the bottom panels of Figure 4 with and Figure 9 demonstrates how this shallower magnitude limit effectively eliminates sensitivity to many recent events which do not contribute significantly to the inner halo. Similarly, the bottom right-hand panel of Figure 4 shows that the group finder recovers no events on eccentric orbits that are likely to have debris beyond the edge of the SDSS survey.

Recall that in our idealized surveys we found that number of groups (Table 2), amount of material in substructures (Table 3) and the probability distributions of group masses (Figure 5) and radial distance (Figure 6) can be related to the number, luminosity function and orbit distribution of recent accretion events. If we now account for the sensitivity of the different surveys to accretion events of different properties outlined in the previous two sections we can examine how such surveys could combine to present a picture of recent accretion history. Note that, because of the low number of groups recovered, we have not included the results for the distribution of group properties in data set S5 in Figures 5 and 6.

For example, our S2 M-giant surveys are not expected to be sensitive to low-luminosity events or those on very radial orbits. This insensitivity accounts for why Table 2 shows that almost no groups were recovered from these surveys of our mock-halos built entirely from radial orbit events, and why there is not a striking difference in the number of groups recovered from our low-luminosity model and our standard ΛCDM models. However, the S2 panel in Figure 5 clearly shows that our high-luminosity model halo gave rise to a much larger fraction in massive groups than our ΛCDM models. Hence, we expect that 2MASS can already provide a reasonable census of recent accretion events from the upper end of the luminosity function.

Previous results for our S3 and S3' data sets suggest that information about the lower end of the luminosity function of recent accretion events might be filled in by these future MSTO surveys. Indeed, Table 2 indicates that these surveys detect a larger number of groups in our low-luminosity halo, which contained a larger number of low-luminosity accretors, when compared to our standard models. This indicates that MSTO surveys should recover debris from accreted objects—like the lower luminosity classical dwarf spheroidal and ultra-faint galaxies—that would be missed by 2MASS.

Since they explore to far greater depth, future RR Lyrae surveys (S4) will have unique power to probe the outer halo, beyond 100 kpc and assess the fraction of debris that lies at these large radii. This power is demonstrated by the very clear difference in the S4 panel of Figure 6 between the radial distribution of groups found in the RR Lyrae survey of our radial orbit halo (magenta line) compared to the ΛCDM cases (thin black lines). This difference implies that RR Lyrae surveys will flesh out our picture of the orbit distribution of accretors onto our Galaxy. Additionally they are also very efficient at detecting the low-luminosity accretion events.

The synthetic surveys that we analyze do not include other galactic components, e.g., the disk and the bulge. However, this will have little impact on the analysis we present here because (1) the volume of the stellar halo that we explore (100 kpc and beyond) is much larger as compared to the volume of the disc and (2) most stellar halo structures lie away from the disc and bulge regions. Note that even in the present analysis we have neglected groups that lie within 5 kpc of sun or the galactic center (Section 3.1). Another thing that we have neglected is the extinction which will make the detection of halo stars in the plane of the disc very difficult. In realistic surveys one might have to filter out the low-lying latitudes as was done for real 2MASS data in Paper I. However an analysis of S2 data after removing the low-latitude regions revealed that this has little impact on the number of reported groups or the ability to discern different accretion history. In general a few structures will be missed while a few others will be split and discovered as two groups.

For the MSTO surveys we have neglected the possibility of contamination by quasars (QSOs). They are expected to contaminate toward the blue color limits. In the color range 0.2 < g − r < 0.4 and r < 21.5, Jurić et al. (2008) estimate the contamination by QSOs to be 10%. Since QSOs would be distributed isotropically and in a smooth fashion they will not give rise to any false structures but due to their contribution to the background they might decrease the significance of some of the faint structures, making them difficult to identify.

A final limitation of our analysis is the fact that recent accretion events with L < 10⁵ L_☉ are not included in the stellar halo simulations that we use here. This is because the star formation recipes in the simulations were fine tuned to reproduce the observed distribution of satellites known at the time of the publication of Bullock & Johnston (2005), and the ultra-faint galaxies were discovered later on. Our analysis shows that the MSTO and RR Lyrae surveys are sensitive to low-luminosity events hence we anticipate such surveys to detect more structures than that reported by us.

In Paper I we identified structures in the 2MASS M-giant sample and found 16 groups, seven of which corresponded to known bound structures, one was probably a part of the disk, and two others were due to masks employed in the data. Excluding these, we found six unbound groups.

The 2MASS results in Tables 2 and 3 and in bold dashed lines in Figures 5 and 6 provide a detailed comparison between the properties of these observed groups and those recovered from the synthetic M-giant surveys (data set S2) of our 11 ΛCDM models (thin black lines). Overall, we find broad consistency between the observations and models. In particular, the total number of groups (six) sits close to the middle of the range found in the models (8.4 ± 8).

Note that the results reported for data set S2 are without any selection by latitude. We also repeated the analysis by applying the same latitude limits used to generate our 2MASS M-giant sample but did not find a significant change in the mean number of recovered groups for ΛCDM halos, or in the figures. However, the fraction of material in groups was found to be around 30%, which is about three times higher than that reported in Table 3 (unbound groups). This is because much of the stellar halo mass is concentrated toward the center and the plane of the galaxy. This mass near the center of the Galaxy does not contain any substructures and its exclusion results in an increase in the mass fraction of recovered groups. On the other hand in the real 2MASS M-giant data we find the mass fraction in groups to be 4.5%, which is quite low. The reason for this is that in the real 2MASS data nearly 75% of the stars are due to the LMC which is a bound satellite system and this is quite un-typical for a ΛCDM halo. For comparison in our simulated halos the maximum fraction of mass in a single bound satellite was 20%. Excluding the LMC, increases the fraction of material in groups to 18%, which is more typical of what we see in simulated halos.

We can also use our synthetic surveys to ask whether the 2MASS groups are likely to correspond to real accretion events and, if so, what type of events. Table 1 shows that the typical purity of groups recovered from the simulated halos in data set S2 was high, which supports the interpretation of the groups in 2MASS as real physical associations. Moreover, the results from Section 5 suggest that these groups correspond to recent (<10 Gyr ago), high-luminosity accretion (>5 × 10⁶ L_☉) events on orbits whose apocenters lie within 100 kpc.

Finally, we can test what the 2MASS groups can contribute to our picture of accretion history by asking what fraction of accreting objects we expect to detect as groups in this data set. Figure 8 suggests that the groups represent 60% of objects accreted within the last 6 Gyr, falling to as little as 10% of those accreted 7.5–8.5 Gyr ago and 0% of those with accreted earlier than 11.5 Gyr ago. While these numbers clearly portray the weakness of the M-giants as tracers of accretion history, they do not capture their strength. Figure 11 repeats Figure 8, but this time restricting attention to events that we know contain M-giants within the apparent magnitude range of the 2MASS survey—those that have high luminosity (>10⁷ L_☉) and are on mildly eccentric orbits (>0.8). 100% of these events are recovered that were accreted within the last 8 Gyr, and more than 50% with accretion times less than 10 Gyr ago. We conclude that the structures in 2MASS give us a fairly complete census of recent, massive accretion events along mildly eccentric orbits.

Zoom In Zoom Out Reset image size

While a full application of our group finder to the SDSS MSTO sample is beyond the scope of the current work, analysis of our S5 data sets suggests that we might expect to find 0–4 groups, corresponding to relatively high-luminosity, intermediate-age accretion events on mildly eccentric orbits. The S2 and S5 panels of Figure 4 imply that the majority (though not all) of these events would also be recovered from the 2MASS M-giant sample. Indeed, four groups that are plausibly associated with the ancestral siblings of the classical dwarf spheroidal satellites are apparent through visual inspection of SDSS MSTO sample (the tails from disruption of Sagittarius, the Monocerous ring and the Virgo overdensity, and the Orphan Stream; see, e.g., Belokurov et al. 2006), and two of these are also seen in the M-giants. This again indicates broad consistency between structures seen in the real stellar halo and those in our models. (While many more streams from globular clusters and lower-luminosity dwarfs have been found in SDSS using matched-filtering techniques, e.g., Grillmair 2009, these objects are missing from the stellar halo models, so we compare only the number of higher-luminosity streams.)

From Table 3 it can be seen that, for data set S5, the mean fraction of material in groups is 4% which increases to 6% when bound groups are also included. These results are consistent with what we see for the S3 survey, if we take into account the shallow depth of the S5 survey due to which the outer halo, which is also more structured (see Figure 7), is missed. Similarly, most bound structures are in the outer parts of the halos (see Figure 7), hence unlike other surveys, including them does not increase the fraction by much for the S5 survey. However, when compared to the results of Bell et al. (2008), where they find the fractional rms deviation from a smooth analytic model, σ/total, to be around 40% for SDSS MSTO stars, our mass fractions are apparently low. It should be noted that Bell et al. (2008) had also analyzed the same set of halos as used by us and reported good agreement with the SDSS data. Hence, the cause of the mismatch is due to the methodologies being different, which we explore below. First, we think that σ/total cannot be directly interpreted as the amount of mass in a structure, e.g., considering equally populated bins, if 10% of the bins differ in mass by order of the mass in the bins one gets σ/total = 0.33. Second, an analytical model as adopted by them might not necessarily be a good description of the smooth component of the halo. This misfit will contribute to σ/total but will not give rise to any structure in our scheme. It is also important to note that in our clustering scheme only those structures are detected which give rise to peak in the density distribution. Third, the group finder truncates a structure when its isodensity contour hits that of another structure. Hence, the envelope region around a structure can contribute to σ/total but is ignored by us. Note that by following points in the envelope regions along density gradients one can associate them with the groups but such a classification is not free from ambiguity and in general decreases the purity of the groups. Finally, in our analysis we only consider significant groups which stand out above the Poisson noise. However, there might be low-significance fluctuations which we have deliberately ignored and these will invariably contribute to σ/total. However, low-significance fluctuations, unless they are massive, which is rare, will not dominate the mass fraction in structures. To conclude we think that, the methodologies being very different, it is very difficult to interpret the results of one scheme in terms of the other.