CHARACTERIZING THE SHARDS OF DISRUPTED MILKY WAY SATELLITES WITH LAMOST

Though the stellar halo of the Milky Way (MW) contains only $\sim 1\%$ of the Galaxy's total stellar mass, long dynamical times in the halo mean that much of the fossil record of the formation and evolution of our Galaxy is preserved in the dynamical and chemical signatures of halo stellar populations. In the prevailing Λ-Cold Dark Matter cosmological paradigm, massive galaxies grow by the agglomeration of many smaller sub-galactic fragments, or subhalos, that contribute their constituent dark matter, stars, and gas to the larger host galaxy (e.g., White & Rees 1978; White & Springel 2000). Through detailed analysis of the kinematics and chemistry of stellar populations in the MW halo, we can thus assess the relative importance of the monolithic collapse suggested by Eggen et al. (1962) as a possible origin for the stellar halo, and the stochastic process of accretion as discussed by, e.g., Searle & Zinn (1978).

Deep, large-area photometric surveys such as the Two Micron All-Sky Survey (2MASS; Skrutskie et al. 2006) and the Sloan Digital Sky Survey (SDSS; York et al. 2000) have laid bare the ubiquitous substructure lurking in the low surface brightness stellar halo. The numerous substructures that have been discovered as stellar overdensities include the Sagittarius (Sgr) tidal stream (e.g., Ibata et al. 1994; Majewski et al. 2003; Belokurov et al. 2006; Koposov et al. 2012), the large, complex overdensity in Virgo (e.g., Vivas et al. 2001; Newberg et al. 2007; Carlin et al. 2012; Duffau et al. 2014), the Triangulum–Andromeda (TriAnd) cloud (e.g., Majewski et al. 2004; Rocha-Pinto et al. 2004; Martin et al. 2007; Sheffield et al. 2014; Price-Whelan et al. 2015), and the Orphan Stream (e.g., Grillmair 2006; Belokurov et al. 2007; Newberg et al. 2010), among others (for a review of currently known structures, see Grillmair & Carlin 2016).

It is thus becoming clear that not only is much of the outer Galactic halo made up of accreted debris, but that this is an ongoing process. Tidal remnants that are visible as number-density enhancements are only the most recent infall events, and do not constitute a representative sample of satellite accretion at all epochs of MW evolution. Fortunately, accretion remnants remain relatively coherent in phase space over much longer timescales (e.g., Helmi & White 1999), making it possible to identify relics of infall events dating back many Gyr ago (note, however, that debris from the most ancient accretion events will have phase-mixed, and become indistinguishable from a smooth halo at present day; see, e.g., Johnston et al. 2008). In addition, the long-lived low-mass stars from accreted dwarf galaxies and star clusters retain the chemical signatures of their parent satellite at the time of their formation. It is thus possible to recreate a portion of the hierarchical formation history of the Galaxy by detailed analysis of the stellar halo, by which one can determine the fraction of the halo that resides in substructure at present, and use this to assess the relative contribution of in situ star formation (i.e., stars formed from gas residing in the deep potential well of the MW) and accretion (stars formed in dwarf galaxy or stellar cluster potentials).

Full cosmological models of MW-mass galaxies with sufficient resolution to include all the relevant physics at all physical scales are challenging.⁷ However, with the combination of cosmologically motivated models and semi-analytic techniques, robust predictions have been made to guide our intuition about the formation of the Galactic halo. Bullock & Johnston (2005, hereafter "BJ05") created model stellar halos from a suite of 11 N-body models of satellites from cosmologically motivated merger trees, to which stellar populations were assigned in a manner that reasonably matches the properties of MW satellites. Because these halos are generated purely from accreted satellites, they offer an extreme case to compare with our MW halo data. Though they are created entirely from satellites, these halos contain a smooth population of older debris in the inner halo (inside ${R}_{{\rm{GC}}}\sim 20$ kpc), with more recent, unmixed accretion events predominating in the outer halo. Analysis of these halos by Johnston et al. (2008) showed that different types of structures (great circles, clouds, and mixed remnants) result from different families of satellite orbits. Substructure dominates in the mock halos at large radii, but is also biased toward cloud-like remnants, which result from disruptions on radial orbits. Johnston et al. (2008) found that the most recently accreted substructures should be more metal-rich than the smooth halo (which consists predominantly of metal-poor, phase-mixed early debris). Finally, these authors concluded that on average $\sim 10\%$ (ranging from 1% to 50%) of stars should be currently visible in substructure, with the fraction depending on the epoch of accretion.

Abadi et al. (2006) used a suite of eight high-resolution cosmological (N-body+hydrodynamic) simulations of MW-like galaxies to explore the outer stellar halo. They found that outside ∼20 kpc, $\sim 95\%$ of the stars are accreted; the in situ stars found outside this radius were kicked up in merger events. This leads to a "break" in the radial surface brightness profile (SBP), where there is excess luminosity beyond the break radius compared to an extrapolation of the inner SBP. This steepening of the mass profile with radius was shown to be consistent with the spatial distribution of MW and M31 outer-halo globular clusters, suggesting that they also originated in accretion events. Zolotov et al. (2009) used high-resolution cosmological N-body+SPH simulations of 4 MW-like galaxies to assess the fraction of in situ halo stars. In the highest-resolution of their simulations, $\sim 20\%$ of the halo stars were formed in situ, and are currently located within ∼20–30 kpc of the galaxy center. These in situ stars form deep in the potential well of the host, and are kicked out via merger events. Zolotov et al. (2009) found that the fraction of in situ halo stars is higher (20%–50%) in the inner halos of hosts with relatively quiescent recent merger histories compared to those with more recent merger activity.

In addition to the detection and characterization of high number density, coherent halo substructures resulting from recent tidal disruption, recreating the merger history of the MW requires systematic searches for kinematically cold structures from older accretion events. Schlaufman et al. (2009) searched for these "Elements of Cold HalO Substructure (ECHOS)" in 137 individual SDSS lines of sight, each covering ∼7 square degrees on the sky. Based on identified statistically significant radial velocity peaks among inner-halo main-sequence turnoff (MSTO) stars (at distances of 10–17.5 kpc from the Sun), these authors estimated that an upper limit of $\sim 1/3$ of the metal-poor MSTO stars in the inner halo reside in ECHOS. Starkenburg et al. (2009) developed a metric (the 4distance) for quantifying the separation of stars in four-dimensional position–velocity space. Their application of the 4distance to the 101 halo giants from the Spaghetti survey (Morrison et al. 2000) found that 20 stars ($\sim 20\%$ of the sample) are in groups, compared to an expectation of 9 from a randomized sample. Thus the authors place a lower limit that more than 10% of halo stars must be in substructure. The 4distance metric was also applied to a sample of ∼4000 BHB stars from SDSS by Xue et al. (2011), who showed that the observations have significantly more 4distance "pairs" than a smooth halo, but are deficient in substructure when compared to the pure-accretion halos of Bullock & Johnston (2005). This is, however, reconciled when comparing only the old stars from the models; thus, the (old, metal-poor) BHB stars must only be tracing a small fraction of the predominantly metal-rich, recent-infall debris. Most recently, the 4distance metric was also applied to a sample of 4568 SEGUE K giants by Janesh et al. (2016). By comparing the number of groups (created via friends-of-friends aggregation of pairs identified by the 4distance) to a smooth model halo, this work showed that $\sim 50\%$ of the halo stars in their sample are identified with groups (though with significant contamination from false positives). The fraction of stars in substructure increases with Galactocentric radius, and is also higher in more metal-rich populations than in the most metal-poor halo stars. A large fraction ($\gt 50\%$ beyond Galactocentric radius of 30 kpc) of the groups are noted by Janesh et al. (2016) to likely be associated with the Sgr stream. A correlation function statistic was applied by Cooper et al. (2011) to the SDSS BHB star sample from Xue et al. (2008) to search for spatial and kinematic correlations in the Galactic halo. The number of BHB stars was found to be deficient at ${R}_{{\rm{GC}}}\gt 30\;{\rm{kpc}}$ compared to the predictions of the Aquarius simulations of MW-like halos, though significant clustering is present in the BHB sample at these large radii that is consistent with model halo predictions.

Statistical characterization of the fraction of substructure in the halo need not be limited to spectroscopically identified samples. Bell et al. (2008) used SDSS DR5 photometry of MSTO stars to fit the density profile of the Galactic halo. Based on the rms residuals about this fit, this work concluded that $\gt 40\%$ of the halo is in substructure. This fraction increases with distance, nearly doubling from 10 to 35 kpc, though much of this rise can be attributed to Sgr and other known substructures. Bell et al. showed that the measured rms is consistent with some of the Bullock & Johnston (2005) models, which lends support to the notion that the MW halo is largely accreted debris. Helmi et al. (2011) analyzed Aquarius models (with stellar populations from the semi-analytic method of Cooper et al. 2010), and found that the rms stellar density in halos with no smooth component is much larger than that observed in SDSS by Bell et al. (2008). Addition of a 10% smooth halo component brings them roughly into agreement, though because the model halos are highly anisotropic, the vantage point of the observer affects this significantly.

In this work, we measure the fraction of substructure among Galactic halo stars observed by the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST) spectroscopic survey. Our method is similar to the 4distance technique developed by Starkenburg et al. (2009), but because the LAMOST depth is not uniform along different lines of sight, we choose not to use a global correlation statistic, but instead identify statistical excesses in separate regions of sky. LAMOST will eventually obtain a nearly complete magnitude-limited sample (at high Galactic latitudes) over a huge contiguous region of the northern sky. Targets for the LAMOST survey are selected with a smoothly varying selection function, making the sample a relatively unbiased data set. The combination of a vast contiguous sky area and the simple selection function make this a valuable data set for identifying substructures on a variety of spatial scales. Though the selection function for a given direction is simple, we nevertheless must model its effects over the entire survey; to do so, we compare all of our results in this work to mock observations in which we apply the LAMOST selection function to the Galaxia model (Sharma et al. 2011) in each region of sky. Finally, we compare our results to accretion-derived halos from Bullock & Johnston (2005), again observed with the LAMOST selection function.

This paper is outlined as follows. In Section 2, we discuss the LAMOST survey and the selection of Galactic halo stars for our substructure search. Section 3 outlines our technique for identifying substructures, which we dub Stellar Halo Accretion Related Debris Structures (SHARDS). In Section 4 we compare our results to the Galaxia model of the smooth halo, and to halos built entirely from disrupted satellites. We obtain estimates for the fraction of halo stars in substructure in Section 5, and characterize the overall properties of the substructures we have identified. We conclude in Section 6 with some context for our results and a discussion of upcoming work that builds upon these results.

The LAMOST survey is an ongoing effort being carried out with the ∼4 m effective aperture Guoshoujing Telescope in northern China (Cui et al. 2012; Zhao et al. 2012). The telescope has 4000 robotically positioned optical fibers arrayed over a 5° diameter focal plane, feeding 16 optical spectrographs that produce spectra with resolution $R\sim 2000$ covering wavelengths ∼3800 < λ < 9000 Å. LAMOST has completed 3 years of survey operations plus a Pilot Survey, and has internally released a total of ∼5.7 million spectra to the collaboration, spanning observation times from 2011 October–2015 May.⁸ Of these, ∼3.1 million are AFGK-type stars (mostly part of the LEGUE survey of Galactic structure; Deng et al. 2012) with estimated stellar parameters (${T}_{{\rm{eff}}}$, $\mathrm{log}g$, [Fe/H]; Wu et al. 2014), and all objects have available radial velocities or redshifts accurate to 5–10 km s⁻¹ (e.g., Luo et al. 2015; Xiang et al. 2015). The survey reaches a limiting magnitude of r = 17.8 (where r denotes magnitudes in the SDSS r band), with most targets brighter than $r\sim 17;$ LAMOST will achieve a nearly magnitude-limited survey to $r\sim 17$ over much of the high Galactic latitude northern sky.

To select a sample of Galactic halo giants, we use stellar distance estimates derived from the LAMOST stellar parameters (Carlin et al. 2015). This work showed that, given the typical uncertainties in stellar parameters given by the LAMOST pipeline, distances to most stars can be derived to $\lesssim 30\%$ accuracy. Our sample of halo giants is selected from the LAMOST database with the following constraints: $0.0\lt {(J-K)}_{0}\lt 2.0$, ${K}_{0}\lt 15.5$, spectra having ${\rm{S}}/{\rm{N}}\gt 5$ (in either the g or r band), derived parameters of $0\lt \mathrm{log}g\lt 3.5$, $3800\lt {T}_{{\rm{eff}}}\lt 6500$ K, and $| Z|$ > 5 kpc.⁹ The color and magnitude cuts exclude stars with spurious measurements or extreme stellar types (whose distances are not well-determined). The $| Z|$ criterion ensures that the stars in our sample are many scale heights beyond the disk plane (the disk scale height is ∼1 kpc; e.g., Robin et al. 1996; Chen et al. 2001; Siegel et al. 2002; Jurić et al. 2008), minimizing the contribution of disk stars to our sample. We are interested in assembly of a sample of well-understood, first-ascent red giant branch (RGB) stars, so the $\mathrm{log}g$ and ${T}_{{\rm{eff}}}$ cuts are intended to remove objects that could be evolved stars (e.g., RR Lyrae, horizontal-branch, and asymptotic giant branch stars), while including sufficiently warm stars to retain the most metal-poor "normal" giants. Finally, the signal-to-noise ratio (S/N) criterion, while generously allowing relatively low S/N spectra into the sample, is still likely removing many legitimate K-giant stars from our sample; Liu et al. (2014) showed that K-giants can be reliably identified from LAMOST spectra even at low S/N. For the current study, we choose to avoid including dwarf contamination that is inevitable in the low S/N K-giant sample, but will include these in future expansions of this work. The distribution of Galactocentric distances (${R}_{{\rm{GC}}}$) for this sample of 10,481 halo RGB stars is shown in Figure 1, and their density on the sky is shown in Figure 2. The sample contains stars extending beyond ${R}_{{\rm{GC}}}\gt 60\;{\rm{kpc}}$ in the halo, with a substantial number of stars between $40\lt {R}_{{\rm{GC}}}\lt 60\;{\rm{kpc}}$. In Figure 2, most of the stars are at high Galactic latitudes well away from the plane (as enforced by our $| Z| \gt 5$ kpc cut), with nearly the entire north Galactic cap region sampled substantially by LAMOST.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We wish to identify substructures that we will refer to as SHARDS. By this we mean structures that constitute statistical excesses in velocity–distance space (if we consider localized sky areas, this is essentially a four-dimensional phase space consisting of 3D position and line-of-sight velocity).

Our goal is to assess the fraction of Galactic halo stars observed by LAMOST that are part of substructures. To do so, we wish to consider not only structures that are clustered spatially on the sky, but also those that clump in velocity and distance. Accretion relics remain coherent in phase space for a longer time than they are visible in configuration space, so using velocity–distance metrics to identify structures should sample more of the accretion history of the MW than simply using spatial clustering and/or distance. Furthermore, the large areal coverage of the LAMOST survey allows us to probe much larger spatial scales than surveys that are limited to localized regions of sky.

We initially considered a method similar to that of Schlaufman et al. (2009), which treated each SDSS/SEGUE plate individually. However, because the LAMOST limiting magnitude is fairly bright ($r\sim 17.8$), our study is limited to giant stars, of which there are very few on each individual LAMOST line of sight (referred to as "plates"). Nevertheless, when the entire $\gt 3$ million-star data set currently available from LAMOST is considered, we have a sufficient number of giants to statistically probe the halo. Freeing ourselves of the constraints of individual plate areas helps us to probe larger spatial scales than were available in the ECHOS searches, but the varying selection function from plate to plate must then be accounted for in our work, as discussed further in Sections 3 and 4 below.

Our method of searching for substructures proceeds as follows. We divide the sky into $12^\circ \times 12^\circ$ regions in right ascension (R.A.) and declination (decl.), covering the entire northern sky visible to LAMOST ($-10^\circ \lt \mathrm{decl}.\lt 60^\circ$). This scale was chosen to provide a large statistical sample in each region, while not being so large that the Galactic stellar populations vary dramatically across the region. In order to identify substructure in velocity and distance, we need to know what the underlying ("smooth") halo populations (if there are any) look like. In other words, we wish to assess what structures exist in excess over what we expect from the smooth halo distributions. We adopt the Galaxia (Sharma et al. 2011) model for this purpose. Galaxia generates expected observational catalogs by sampling from empirically based density and velocity prescriptions for MW stellar populations. In particular, Galaxia uses an oblate power-law halo density profile (following Robin et al. 2003) with a velocity ellipsoid of $({\sigma }_{r},{\sigma }_{\theta },{\sigma }_{\phi })=(141,75,75)$ km s⁻¹. Each instance of the model produces a predicted set of observations given prescribed color/magnitude ranges and a region of sky. Sampling from the underlying distributions begins with a random number seed, which can be changed to generate distinct model observations in a given sky area. Thus, to sample the smooth distributions, we choose to run multiple instances of the Galaxia model initialized with different random number seeds, and build a set of mock catalogs with which to statistically compare the LAMOST observations.

We stress that our comparisons throughout this work tell us how much substructure is present relative to the Galaxia model predictions for the underlying halo, so that some fraction of any observed discrepancies may be due to deficiencies in the Galactic components encoded in the model.

Finally, we must account for the inhomogeneous way that the LAMOST survey targets are selected. Plates in the LAMOST survey are separated by magnitude ranges into "Very Bright" (VB; $r\lt 14$), "Bright" (B; $14\lt r\lt 16.8$), "Medium" (M; $16.8\lt r\lt 17.8$), and "Faint" (F; $r\gt 17.8$) plates so that a variety of observing conditions can be accommodated, while also limiting the magnitude range included in each plate to ensure that sufficient S/N is achieved for a large fraction of stars in each observation. This means that the magnitude range probed in a given region of the sky may differ dramatically from a neighboring region. In such cases, using only the distance and velocity distributions of observed stars without accounting for the selection function of stars that were targeted may strongly bias the results. We thus need to "observe" the mock catalogs using the LAMOST selection function in order to obtain a valid comparison.

The process of searching for excess stellar substructures proceeds as follows.

• 1.  
  Run 100 simulations (each initialized with different random seeds) for that field with Galaxia (including only thick-disk and halo components, and using the same constraints that were applied to isolate halo stars from the LAMOST data; see Section 2).
• 2.  
  In each simulation, randomly select a star that is within 0.025 mag in ${K}_{{\rm{S}}}$ versus $(J-{K}_{{\rm{S}}})$ color–magnitude space (or, if the nearest star is more than 0.025 mag away, try 0.05 mag, then if still none, simply select the nearest) of each LAMOST object in the Galaxia output. In this way, we "observe" the Galaxia model in the same way that LAMOST has observed that part of the sky. This results in a sample consisting of 92% mock halo stars and 8% thick-disk stars (averaged over all of the mock observations, and over the entire sky), as tagged by Galaxia.
• 3.  
  Bin the simulated observations into 25 km s⁻¹ by 5 kpc bins in the ${V}_{{\rm{GSR}}}$ versus distance plane, then calculate the average number of stars found in each of these bins over the 100 simulated observations of that field of view. We take the standard deviation over the 100 mock catalogs as an error bar on the counts in each bin.
• 4.  
  Bin the observed LAMOST data in the same way, and calculate the residual number of stars for each bin, expressed as a deviation in the sense ${n}_{\sigma }=({N}_{{\rm{LAMOST}}}-{N}_{{\rm{model}}})/{\sigma }_{{\rm{model}}}$, where ${N}_{{\rm{LAMOST}}}$ is the number of observed stars in each bin, ${N}_{{\rm{model}}}$ is the average number predicted by the model, and ${\sigma }_{{\rm{model}}}$ is the standard deviation of the model counts over the 100 simulations.
• 5.  
  Define an "excess" as a bin that has $\gt 3\sigma$ excess of observed stars. Rather than trying to identify which stars should be removed as "background" stars, we retain all stars from each bin in which we have identified an excess. Because we expect some excesses simply due to statistical fluctuations (and we are not considering deficient bins where there are fewer stars than predicted), this method includes more stars than are actually present in SHARDs. We explore this effect further in Section 4.2, and assess the fraction of "smooth halo" stars that are identified as "excess" by our method. For the analysis in Sections 4 and 5, we subtract the expected number of background stars to derive the net number of excess stars in SHARDS.

Figure 3 shows examples of two fields from our study. The upper row is field 74, centered at (R.A., decl.) = (1620, 200), and the lower row is field 78, at (R.A., decl.) = (2100, 200). The left panels in each row show the observed line of sight velocity distribution (in a Galactocentric frame), ${V}_{{\rm{GSR}}}$, from LAMOST as a solid black histogram. The dashed line in each of these panels shows the mean velocity from the 100 Galaxia model realizations in that field, with the $2\sigma$ region given by the gray shaded region (where σ is the standard deviation of the model counts in each bin). These panels show a stark difference in the amount of velocity substructure exceeding the $\gt 2\sigma$ level. Field 74 (the upper panel) shows an obvious excess centered at ${V}_{{\rm{GSR}}}\sim -100$ km s⁻¹ that exceeds the $2\sigma$ shaded region in four consecutive bins. This field also has a high-velocity peak that exceeds the expected velocity. In contrast, field 78 (lower panels) has no velocity peaks exceeding the $2\sigma$ region in the ${V}_{{\rm{GSR}}}$ histogram.

Zoom In Zoom Out Reset image size

We now show that with the addition of distance information to the clustering search, substructures emerge even in cases such as field 78, which shows no obvious velocity structure in the lower-left panel of Figure 3. The center panels of both rows in Figure 3 display 2D histograms with bins of 25 km s⁻¹ in ${V}_{{\rm{GSR}}}$ and 5 kpc in (heliocentric) distance. The color of each bin corresponds to the average number of stars in that bin from the 100 model realizations, and the black points are the LAMOST data. Given this binning and the standard deviation, σ, for each bin, one can simply compare the expected number count to the observed count from LAMOST. The deviation from model expectation, ${n}_{\sigma }=({N}_{{\rm{LAMOST}}}-{N}_{{\rm{model}}})/{\sigma }_{{\rm{model}}}$, is mapped in the right panels of Figure 3, where the red-orange colors represent excesses relative to the model, and blue bins are deficits. Comparing the 2D maps to the velocity histograms, it is clear that the obvious velocity peak in field 74 (upper row) is also coherent in distance, as expected for a tidal debris structure. This feature is most likely related to the Sagittarius stream, and has velocity and distance consistent with those from the model of Law & Majewski (2010) and the SDSS observations of Belokurov et al. (2014). Note also that some features that appear as strong excess bins in the right-most panels contain only a single star; because the average occupancy of many bins (given the model convolved with the LAMOST selection function) is well below one star for many bins, a single star may be a statistically significant excess. However, for subsequent analysis we require bins to have two or more stars for consideration as SHARDS.

While field 78 (lower panels in Figure 3) shows no obvious velocity structure in the left-most panel, there are indeed significant excesses in velocity–distance phase space. We define SHARDS as all bins containing at least two stars and greater than $3\sigma$ excess in position–velocity phase space (i.e., bins in the right column of Figure 3 that are above $3\sigma$), and select all observed stars from these bins (note that in Figure 3, bins with 3 or more stars are outlined with a thick solid line, and those with two stars with a dashed outline). In field 74, which has Sagittarius debris obviously contributing significant substructure, 68 stars out of a total of 202 halo giants reside in SHARDS (after subtracting the background contribution in each bin). Field 78, which shows no obvious velocity structure in the histogram, nevertheless contains 17 stars in SHARDS out of a total of 263 in the field. This illustrates the value of using additional dimensions of information to seek substructures.

We divide the sky into 180 regions that are $12^\circ \times 12^\circ$ in R.A. and decl., spaced 12° apart on the sky (with centers at Declinations of $-4^\circ ,8^\circ ,20^\circ ,32^\circ ,44^\circ ,$ and 56°). Of these, there are 112 regions that have more than 10 LAMOST halo giants, the threshold we applied for inclusion in our search, with a total of 10,481 stars. Then, we perform the same search for SHARDS in each of these regions as in the fields illustrated above. This yields a total of 1140 (2065) stars in $\gt 3\sigma$ ($\gt 2\sigma$) excesses (after subtracting the model-predicted "expected" number), or $\sim 10.9\%$ ($\sim 19.7\%$) of the LAMOST halo giants identified with SHARDS.

3.1. Caveats

To place our observations into context, we first examine the results from running our code on the model halos from Bullock & Johnston (2005; BJ05). These halos are built purely from accreted satellites generated from cosmologically based merger trees and accretion histories. We thus use our mock observations of these halos as a test case to see what one might expect to observe if the Galactic halo consists entirely of satellite debris.

Another test we use to place our observations in context is to generate model halos using only the smooth halo and thick disk prescriptions in Galaxia, then observe these halos in the same way that LAMOST samples the sky. This provides a picture of how many statistical fluctuations and chance groupings are detected as SHARDS by our technique in the case of a smooth halo.

4.1. Comparison to Bullock & Johnston "Pure Accretion" Halo

Even if the Galactic halo is made up entirely of accreted satellites, the debris from ancient accretion events may be phase mixed, and thus no longer visible as substructure in phase space. Given that most of the accreted satellites were assimilated long ago, the fraction of halo stars that is in observable substructure will not be 100%. Furthermore, the fraction that is in substructure at present depends on the specific merger history of the Galaxy. To assess what we expect to see in the case that the Galactic halo is made up wholly of satellite remnants, we use the model halos from BJ05. In this work, cosmologically motivated initial conditions were used to give the properties of accreted satellites, each of which was then modeled via N-body simulations. The results of all of these N-body models for a given host galaxy were compiled together to create a simulated stellar halo created from the simulated, accreted satellites. The satellite properties in BJ05 match those of known MW dwarfs, and the resultant halos have roughly the same total luminosity as the MW. In these models, the destroyed satellites that contribute most of the stellar halo are accreted early (at least ∼9 Gyr ago), and their debris is mostly smoothly distributed and located in the inner halo (inside ${R}_{{\rm{GC}}}\sim 10$ kpc). Most visible substructure is seen in the outer (${R}_{{\rm{GC}}}\gt 20$ kpc) halo, where more recent accretion events deposit their debris.

The 11 simulated halos from BJ05 have been implemented in Galaxia, so that we can take simulated observations of these models and then search for SHARDS in them using the same algorithm we have used for the LAMOST data. For each of the 11 simulations, we generate a catalog using Galaxia, then observe this model by selecting the star that is closest to each star from LAMOST in K₀ versus ${(J-K)}_{0}$. This mock observed catalog is then passed through the routines outlined above to detect all SHARDS in the BJ05 model (given the LAMOST selection function). An example of the results for Field 74 from BJ05 models 1 and 8 (note that these two halos were selected arbitrarily) is given in Figure 4. These mock observational results can be compared directly with those from the top row of Figure 3, which shows the LAMOST results for Field 74. The velocity distribution (left panels of Figure 4) of the BJ05 halos is clearly broader than the expectations from Galaxia, with more stars in the wings of the ${V}_{{\rm{GSR}}}$ histograms and a deficit (relative to the MW halo prediction) near ${V}_{{\rm{GSR}}}\sim 0$ km s⁻¹. In the middle and right-hand panels of Figure 4, it becomes clear that most of the high-velocity excess in the mock halos is nearby, well-mixed debris, rather than distant distinct substructures. One obvious difference between these mock halos and the LAMOST observations in Field 74 is the dominance of the coherent Sagittarius debris structure in the LAMOST data; no such obvious substructures are visible in the BJ05 models. There is a feature in halo 8 (lower panels) that stretches from roughly (${V}_{{\rm{GSR}}}$, dist) $\sim (250$ km s⁻¹, 25 kpc) to (${V}_{{\rm{GSR}}}$, dist) $\sim (100$ km s⁻¹, 45 kpc), similar to the caustic-like features expected from some accretion relics (e.g., Johnston et al. 2008; Sanderson & Helmi 2013). However, in general these mock halos show low levels of substructure at all radii, as might be expected from a melange of accretion debris. Of 202 halo giants in Field 74, there are 13 and 31 stars identified in SHARDS in halos 1 and 8 from BJ05, respectively. In the LAMOST observations, this number is 68 of 202; however, ∼30–35 of these are likely related to the Sgr stream, so that the total is ∼35 of non-Sgr stars in SHARDS in this field. Thus the pure-accretion halo predictions are consistent with the number seen in LAMOST if we exclude Sgr, but predict fewer stars in SHARDS than observed when Sgr debris is included. In fact, in only one (halo 7, with 79) of the field 74 regions of the BJ05 halos do we find as many stars in SHARDS as we have observed with LAMOST.

Zoom In Zoom Out Reset image size

4.2. Comparison to a "Smooth" Model Halo

Given the sparse selection of halo giants from LAMOST, it is likely that some stars will be identified as SHARDS simply due to statistical noise from sampling the smooth halo distribution. In addition to the inhomogeneous subsampling of the halo, uncertainties in the distances (typically about 20%–30%, related mostly to stellar parameter errors; see Carlin et al. 2015) can also shift stars from their true location in the 2D plane of ${V}_{{\rm{GSR}}}$ versus distance, making some stars look like outliers compared to the expected populations. To quantify these effects, we generate a model halo using only the smooth prescriptions for the halo and thick disk that are implemented in Galaxia. By performing mock observations of these halos in the manner described in Section 4.1, we generate a catalog of what would be expected in the LAMOST halo-giant sample if the underlying distributions are smooth.

We run this observed catalog of "smooth halo" stars through our SHARDS pipeline. Plots of the results for field 74 are shown in Figure 5. The velocity histogram reproduces the underlying distribution well, with only a handful of bins making excursions larger than one sigma. However, our method does find a number of stellar excesses in this field. We find that 18 of 202 stars in Field 74 are in SHARDS (compare to 68 of 202 in LAMOST, or ∼35/202 if Sgr is excluded). We repeat this exercise for 10 different Galaxia smooth halo models generated with different random seeds. The average number of stars in SHARDS is 11/202 (ranging from 5 to 18) in Field 74, while the average number in the 11 BJ05 halos in this same field is 40/202 (range: 13–79). This gives an idea of the level of "background" contributed by sampling the smooth halo with the LAMOST selection function, as well as the simplistic method of denoting all stars in bins with $\gt 3\sigma$ excesses of two or more stars as SHARDS. The smooth halo simulation and the BJ05 pure accretion halos predict a range that helps us to interpret the number of LAMOST halo stars in SHARDS according to our algorithm.

Zoom In Zoom Out Reset image size

5.1. Fraction of Halo Stars in Substructure

One of the main questions we wish to answer is what fraction of stars in the Galactic halo reside in detectable phase-space substructures. While the determination of an absolute fraction is beyond the scope of this work, we can assess what fraction of the LAMOST halo stars is part of substructure, and compare this to expectations for what a similar survey would observe from the BJ05 mock halos.

Figure 6 shows the results of our substructure search in the LAMOST DR1-3 data. Each of the panels in this figure shows the fraction of observed stars that are identified as SHARDS in each 5 kpc bin in Galactocentric radius. To compute this fraction, we summed all stars in SHARDS, then subtracted the expected number of stars in each bin. This fraction thus consists of the number of excess stars in substructures. The upper row shows the fraction of all stars in bins containing two or more stars that are identified as excesses above three sigma, and the lower panels require groups to have 3 or more stars for inclusion. All panels show the observed fraction of substructure in LAMOST as a solid black histogram. The panels on the left show results from all 11 of the BJ05 halos as lines of different colors, with the average and standard deviation as a dashed black line and shaded gray region. On the right side, the dashed lines represent the average of 10 simulations of the smooth halo, with the gray shaded region encoding 1σ variation (where σ is the standard deviation) about this mean.

Zoom In Zoom Out Reset image size

The upper-right panel of Figure 6 compares the LAMOST observational results (solid black line) to the average of 10 smooth halo realizations (dashed line with 1σ shaded region), with bins required to have more than a single star in each of the SHARDS. It is immediately clear that the stochastic sampling of the halo by LAMOST yields some background of excess stars relative to the Galaxia model predictions, but that the likelihood of detecting multiple stars in a $\gt 3\sigma$ excess in the smooth halo is slim, yielding an average of less than 10% of stars in SHARDS over all radii. Thus, the fact that the LAMOST substructure fraction (solid line) is larger than 10% over nearly all radii less than 50 kpc signals a significant difference between the smooth halo and what we have observed. This fraction is consistent with the 10% lower limit on the fraction of Spaghetti survey data that were estimated by Starkenburg et al. (2009) to be associated with halo substructures.

The upper-left panel of Figure 6 compares SHARDS that contain multiple (i.e., $\gt 1$) LAMOST halo stars to similar detections in the 11 BJ05 model halos. From this figure we conclude that, given the LAMOST selection function and our SHARDS algorithm, the amount of substructure in the Galactic halo is consistent with model halos built wholly of accreted satellites. Between $15\lesssim {R}_{{\rm{GC}}}\lesssim 60\;{\rm{kpc}}$, the fraction of LAMOST stars in SHARDS follows closely the fraction from BJ05. Given that there is also an excess above the smooth halo in this distance range, it is clear that we are seeing a signal of substructure. The fraction of substructure in the pure-accretion BJ05 models is found by our method to be as low as ∼10%–20% at all radii (on average); this arises because of a combination of the LAMOST selection function and the fact that debris from the earliest accretion events is mostly well-mixed at present, and thus is no longer visible as four-dimensional substructure by our method.

To isolate a sample with little expected contribution from false positives, the bottom panels of Figure 6 show the fraction of stars in substructure, but including only stars that are in stellar excesses of 3 or more stars. As expected, the smooth halo (lower-right panel) has very few detections of such groups, while there are numerous SHARDS in LAMOST that contain $\geqslant 3$ stars out to ${R}_{{\rm{GC}}}\sim 40\;{\rm{kpc}}$ (note that the lack of such stars beyond 40 kpc is likely simply due to the paucity of observed stars at such distances). The $\gtrsim 10\%$ fraction of LAMOST stars in substructures with more than 3 stars exceeds the fraction derived from running our algorithm on the BJ05 models (lower left panel) over radii $20\lesssim {R}_{{\rm{GC}}}\lesssim 35\;{\rm{kpc}}$. This may be due to the presence of the Sgr stream over much of the LAMOST footprint; few of the BJ05 model halos contain such a massive, late-infalling satellite.

Figure 7 summarizes our conclusions from Figure 6, focusing on the results for SHARDS that contain 2 or more stars (left panel) and 3 or more stars (right panel), and are thus securely identified as substructures with little contamination from false positives. The solid black line in both panels of Figure 7 shows the fraction of LAMOST stars that are part of SHARDS as a function of Galactocentric radius. For comparison, the dotted–dashed line shows number of "false positives" as the average fraction of stars in SHARDS from the 10 simulated smooth halos, with the 1σ variation among models given by the cyan shaded region. The average of the eleven BJ05 accreted halos is also shown as a dashed line, with the variation among these mock halos shaded in magenta. The fraction of stars in substructure in the LAMOST data exceeds that in the smooth halo at all radii out to at least 40 kpc. In fact, the LAMOST results exceed even the purely accreted halos at radii between $20\lesssim {R}_{{\rm{GC}}}\lesssim 35\;{\rm{kpc}}$, and the level of substructure is consistent at nearly all radii. The fraction of LAMOST halo stars that reside in SHARDS is consistently about 10% above the floor set by the false positives at all radii out to ${R}_{{\rm{GC}}}\lesssim 40$ kpc. Beyond 40 kpc, there are too few stars in our LAMOST catalogs to draw strong conclusions.

Zoom In Zoom Out Reset image size

In Figure 8 we show the background-subtracted group size of SHARDS from LAMOST (black histogram), the accreted halos from BJ05 (blue dashed line), and the smooth Galaxia halos (dotted–dashed red histogram). As expected, there are more large groups ($N\gtrsim 3$ stars) in the LAMOST and accreted-halo results than in the smooth halo (i.e., the "false positives"). Between $3\lt N\lt 10$ stars per group, the fraction of stars in SHARDS detected in LAMOST exceeds even the fraction expected in the purely accreted halos. However, there are more stars in groups with $\gt 10$ stars in the BJ05 halos than in LAMOST, suggesting that there may be some intrinsic difference between the clustering scale of SHARDS and those predicted by models of accretion-derived halos.

Zoom In Zoom Out Reset image size

We conclude from analysis of Figures 6–8 that the population of halo stars in LAMOST is clearly inconsistent with being drawn from a smooth halo. Indeed, the fraction of MW halo stars in substructure is consistent with halos (from BJ05) created entirely from accreted satellites. Our conservative choice requiring at least 3 stars in each SHARD places a lower limit of $\sim 10\%$ of halo stars in substructure at radii ${R}_{{\rm{GC}}}\lt 40\;{\rm{kpc}}$.

5.2. Comparison of Global Properties of SHARDS to Non-SHARDS Halo Stars

We now compare the properties of stars identified as part of SHARDS in the LAMOST halo sample to those that are not members of SHARDS. Figure 9 shows the metallicity distribution of all stars that are not part of SHARDS as a solid line; members of SHARDS with three or more stars are shown as the dashed histogram.¹⁰ The SHARDS are shifted toward higher metallicities than the stars that are not in substructure, with a deficit of stars at [Fe/H] ≲ −1.6 and an excess in all bins with [Fe/H] > −1.1 relative to the stars that are not in SHARDS. Schlaufman et al. (2011) found a similar trend in the ECHOS from SDSS; the mean metallicities of ECHOS are more metal-rich than the average [Fe/H] of MSTO stars along the same line of sight. The metallicity of debris is correlated with satellite luminosity (Johnston et al. 2008), such that more luminous satellites are more enriched than their fainter counterparts. Furthermore, intact dwarf galaxies in the Local Group exhibit a clear luminosity-metallicity relation (Kirby et al. 2013), with satellites at [Fe/H] = −1.0 at luminosities of $L\sim {10}^{8}{L}_{\odot }$. Because the substructures that are observed as SHARDS (or ECHOS) result from relatively recent accretion events, their bias toward more metal-rich populations suggests that the late-infalling satellites contributing to the outer Galactic halo have been predominantly luminous dwarfs rather than metal-poor ultra-faint dwarf spheroidals. Of course, the most prominent substructure in the halo—the Sgr stream—is an ongoing accretion of a metal-enriched, luminous satellite, which may contribute many of the metal-rich SHARDS we have detected.¹¹ Indeed, typical metallicities in the Sgr stream range from about $-1.2\lt {\rm{[Fe/H]}}\lt -0.4$ (see summary in Law & Majewski 2016), which is similar to the metallicity range of the excess stars in SHARDS.

Zoom In Zoom Out Reset image size

There is evidence for two components in the Galactic halo (the inner and outer halo; Carollo et al. 2007, 2010; Beers et al. 2012), with a transition at ${R}_{{\rm{GC}}}\sim 20\;{\rm{kpc}}$. In simulations of MW-like galaxies, Abadi et al. (2006) showed that a break in the SBP occurs naturally at ∼20 kpc from the galaxy center, where within this radius most of the stars are formed in situ, and kicked up during merger events, and outside 20 kpc, $\sim 95\%$ of stars are accreted. The accreted stars contribute excess luminosity beyond the break relative to an extrapolation of the inner SBP. In Figure 10, we compare the radial distribution of the "smooth" component in our LAMOST sample (stars that are not in SHARDS) to stars that are in SHARDS. The SHARDS (dashed line) deviate slightly from the "not-SHARDS" sample (solid histogram) at radii of ∼20–40 kpc, with the SHARDS more prominent in this range. While we have no way of distinguishing accreted versus in situ halo stars, this break in the relative numbers of SHARDS versus non-SHARDS stars suggests that accretion is a more important contributor to the halo beyond ${R}_{{\rm{GC}}}\sim 20\;{\rm{kpc}}$ than within this radius.

Zoom In Zoom Out Reset image size

We have presented a technique for the statistical identification of substructure in the Galactic halo using spectroscopically confirmed RGB stars from LAMOST. Our method relies on comparison to the Galaxia model (Sharma et al. 2011) along each line of sight, observing the model halos with the LAMOST selection function to ensure a valid comparison. We conservatively estimate that $\gtrsim 10\%$ of the MW halo stars from LAMOST are part of substructures that we refer to as SHARDS. We quantify the significance of our substructure detection by comparing to a smooth-halo model (Galaxia) and to mock halos consisting entirely of accreted satellites (from Bullock & Johnston 2005, BJ05). We find that the fraction of halo stars in substructure exceeds what is expected from statistical fluctuations due to incomplete sampling of a smooth halo over all Galactocentric radii ${R}_{{\rm{GC}}}\lt 40\;{\rm{kpc}}$. Additionally, the LAMOST substructure fraction follows closely the fraction of BJ05 stars that are identified with SHARDS, suggesting that the Galactic halo is consistent with a purely accreted origin.

The SHARDS we have detected are in general more metal-rich than the halo stars that are not part of substructure. The radial profile of SHARDS differs from that of the non-SHARDS population beyond ${R}_{{\rm{GC}}}\sim 20\;{\rm{kpc}}$, suggesting a break in stellar populations that may be due to a transition from well-mixed early accretion debris and in situ stars in the inner halo to late-infall accretion dominated populations in the outer halo.

In this work, we are interested in the global characteristics of the Galactic halo, but of course many known substructures are likely present in our catalog of SHARDS. A brief examination shows evidence of the Sgr stream, TriAnd, and the Virgo overdensity, and there are likely many others within this data set. We defer discussion of the previously known structures to later work, where we plan to examine them in more detail. Likewise, a detailed comparison of the chemical abundances of SHARDS to those of halo stars that are not in substructure will come in future work. Currently, the LAMOST stellar parameters pipeline derives only a bulk metallicity, [Fe/H], but not detailed abundances such as [α/Fe] or [C/Fe]. Lee et al. (2015) have shown that quality measurements of these abundances can be readily achieved from LAMOST spectra, so we expect in the near future to analyze the abundance signatures of SHARDS, as was done by, e.g., Schlaufman et al. (2012) for ECHOS from SDSS. Chemical abundances may also provide an additional dimension that can be used to distinguish accretion relics from in situ halo populations. As the LAMOST survey continues to fill more contiguous area to more uniform depths, it will be possible to characterize the level of substructure on different spatial scales, with which we can explore the luminosity function and orbit types of infalling satellites that have contributed to the halo (Johnston et al. 2008).

Our method to identify substructure is similar in many respects to other techniques that have been used (e.g., the 4distance—Starkenburg et al. 2009; Janesh et al. 2016; ECHOS—Schlaufman et al. 2009; two-point correlation function: Cooper et al. 2011; Xue et al. 2011). Our technique takes advantage of the large filling factor of the LAMOST survey on the sky, and we account for the complicated selection function of the survey. We conclude that, beyond ${R}_{{\rm{GC}}}\gt 20\;{\rm{kpc}}$, the fraction of halo giants from LAMOST that are in substructure is consistent with expectations from stellar halos (Bullock & Johnston 2005) built entirely from accreted satellites.

We thank the referee for valuable comments that greatly improved this work. J.L.C. is grateful for the hospitality and support from Chao Liu and Licai Deng during his summer in Beijing, and for the Chinese Academy of Sciences President's International Fellowship that supported this work. This work was supported by the U.S. National Science Foundation under grants AST 09-37523 and AST 14-09421, and by The Marvin Clan, Babette Josephs, Manit Limlamai, and the 2015 Crowd Funding Campaign to Support Milky Way Research. C.L. also acknowledges the Strategic Priority Research Program "The Emergence of Cosmological Structures" of the Chinese Academy of Sciences, grant No. XDB09000000, the National Key Basic Research Program of China, grants No. 2014CB845700, and the National Science Foundation of China, grants No. 11373032 and 11333003. T.C.B. acknowledges partial support for this work from grants PHY08-22648; Physics Frontier Center/Joint Institute or Nuclear Astrophysics (JINA), and PHY 14-30152; Physics Frontier Center/JINA Center for the Evolution of the Elements (JINA-CEE), awarded by the US National Science Foundation. Y.W. is supported by the National Science Foundation of China (NSFC) under grant No. 11403056. Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope; LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

This research made use of the publicly available software Galaxia and Topcat (Taylor 2005), the IDL Astronomy Library (Landsman 1993), and NASA's Astrophysics Data System. This work relied on the zen computing cluster at the National Astronomical Observatories of China, Chinese Academy of Sciences in Beijing. We thank Ricardo Muñoz for valuable discussion of an early draft of this work.