A Low-Mass Stellar-Debris Stream Associated with a Globular Cluster Pair in the Halo

According to the $\Lambda$CDM cosmological model, the Milky Way (MW) has grown to its current size through mergers with numerous neighboring dwarf galaxies. Thanks to the advent of the $Gaia$ mission (Gaia Collaboration et al., 2018), the stellar debris from relatively massive merger events such as the $Gaia$-Enceladus-Sausage (GES; Belokurov et al., 2018; Haywood et al., 2018; Helmi et al., 2018; Myeong et al., 2018a) and the Sequoia (Seq; Myeong et al., 2019) have been identified in the inner stellar halo. In the outer stellar halo, the full 6D panoramic portrait of the Sagittarius (Sgr) stream has been obtained for the first time (Antoja et al., 2020; Ibata et al., 2020; Ramos et al., 2020). These well-studied substructures and streams are the fossils from dwarf galaxies with dark matter halo masses of 10¹⁰ – 10¹¹M_⊙. The relics from less-massive dwarf galaxies engulfed by the MW are far more difficult to identify. However, based on the hierarchical paradigm of galaxy formation, the majority of the building blocks of the MW are expected to be small ($\lesssim$ 10⁹M_⊙). The identification of numerous minor mergers is thus essential for unraveling the complete assembly history of the MW.

Also importantly, low-mass dwarf galaxies have relatively short star-formation histories, and thus can provide direct records of the high-redshift ($z\sim$ 5) Universe; the epoch when globular clusters (GCs) were also formed. These clusters later fell into the MW as their host galaxies were disrupted, thus we would expect MW GCs to be connected to halo substructures. Indeed, such associations have been seen in M31 from recent photometric studies (see, e.g., Mackey et al., 2010; Huxor et al., 2011; Mackey et al., 2019). Similar relationships are less obvious in our Galaxy, as their detection relies on spectroscopic data for numerous individual stars. The few substructures in the MW known to be associated with GCs, the Sgr stream, the GES, and Sequoia, are all expected to have originated from massive accreted dwarf progenitors. Although GCs likely also populated less-massive progenitors, as has been found in nearby dwarf galaxies (see, e.g., Georgiev et al., 2010), their dwarf progenitors are expected to have been fully disrupted in the outer halo before sinking deep into the Galactic potential.

In order to identify the substructures associated with GCs stripped from lower-mass progenitors, we employ a sample of halo stars comprising two types of old stars, blue horizontal-branch (BHB) and RR Lyrae (RRL) stars. Such stars are not only representatives of the ancient halo, but are also excellent tracers of structure, owing to their precise distance estimates. Previous studies of substructure identification in dynamical space are limited to the inner-halo region, due to the lack of good distance estimates for more distant stars. The large range of distances of this halo sample, with an uncertainty as low as 5$\%$ (based on photometry only), allows us to identify dynamically tagged groups (DTGs) in the outer halo. In Sec. 2, we combine the SDSS BHB and SDSS+LAMOST RRL catalogs, and cross-match with $Gaia$ DR2 (Lindegren et al., 2018), yielding $\sim$ 7600 stars with full 6D phase-space information. The group-identification approach is discussed in detail in Sec. 3. The DTGs with GC associations are presented in Sec. 4, including both existing and newly identified substructures. A summary and brief conclusions are provided in Sec. 5.

We combine a previous SDSS BHB catalog (Xue et al., 2008) with the recently released SDSS+LAMOST RRL catalog (Liu et al., 2020) to create a halo sample of 7640 stars that have full 6D kinematic information available. For BHB stars, line-of-sight velocities are derived from the SEGUE Stellar Parameter Pipeline (SSPP; Lee et al., 2008a, b), with uncertainties of 5 km s^-1 to 15 km s^-1 (Xue et al., 2008). The velocities of RRLs are taken from Liu et al. (2020), which utilizes empirical templates to fit velocity curves of multiple measurements from the LAMOST (Deng et al., 2012; Zhao et al., 2012) and the SDSS/SEGUE (Yanny et al., 2009) surveys. Depending on the number of measurements, the velocity precision varies from 5 km s^-1 to 15 km s^-1. The distance estimates for both types of stars are obtained from multi-band photometry with mean uncertainties of about 5 $\%$ (Xue et al., 2008, 2014; Liu et al., 2020). We then cross-match the halo sample with $Gaia$ DR2 (Lindegren et al., 2018). Given the magnitude range of the sample ($G\sim$ 17 – 19), the errors of the proper motion measurements range from 0.13 to 0.60 mas yr^-1. Combining with the distance uncertainty of 5$\%$, the typical transverse velocity uncertainty is about 15 km s^-1 for the majority of stars in the sample, located 10 to 20 kpc from the Sun. This is equivalent to the uncertainty of the line-of-sight velocities. The resulting errors in the orbital parameters and other dynamical properties are sufficiently small to enable detection of significant groups in dynamical space. For the MW GCs, we employ the catalog from Harris (2010), with proper motions determined by Vasiliev (2019a).

We apply the neural-network-based clustering method StarGO (Yuan et al., 2018) to search for substructures that are clustered in the 4D space of orbital energy and angular momentum, both of which are approximately conserved, even in non-spherical potentials (e.g., Helmi & de Zeeuw, 2000). The gravitational potential of McMillan (2017) is used to derive dynamical parameters with AGAMA (Vasiliev, 2019b). As in the previous work from Yuan et al. (2019, 2020), we use ($E$, $L$, $\theta$, $\phi$) as the input space, where the latter two angular parameters characterize directions of the orbital poles, and are defined as:

We use a 100$\times$100 neuron network, and follow a similar recipe as Yuan et al. (2020) to identify dynamical groups. Each grid point of the neuron map hosts one neuron, as shown in the left panel of Figure 1, which has an initially randomized 4D weight vector. For a given input vector, the neuron that has the weight vector closest to it is defined as its best-matching unit (BMU). Each neuron updates its weight vectors to come closer to the input vector in the 4D space; the learning effectiveness depends on its distance to the BMU on the 2D map. The weight vectors of the neurons close to the BMU will be assimilated into the input vector more efficiently compared to those neurons located farther away. The final result of the learning process is a converged map, after a sufficient number of iterations. This process preserves the structure of the input data, and projects it onto a 2D map.

The learning results can be revealed by differences in the weight vectors between adjacent neurons, which are defined as a 100$\times$100 $u$-matrix. In the left panel of Figure 1, the gray-colored neurons have the top 20$\%$ $u$ values, denoting the lowest 20$\%$ similarities between neighbors. These neurons form gray boundaries, and separate the others into isolated islands (see the white patches in the left panel). Compared to the boundary neurons, those in islands have more similar weight vectors, and thus correspond to stars clustered in dynamical space. The idea of group identification is to find the islands isolated by the gray boundaries as we scan the threshold ($u_{\rm thr}$) that defines the boundary. We check the significance and contamination of neuron groups at each threshold value when they appear as isolated islands. By this means, we are able to systematically identify all the significant groups of stars clustered in the input space. In Figure 1, we show all of the four groups identified in this work, for three different values of $u_{\rm thr}$. They are plotted by different colors on the neuron map shown in the left panel. The corresponding star groups form separate clusters in the input space shown in the right two panels.

Evaluation of the significance and contamination of each detected group is implemented as a post process. We first draw a Monte Carlo (MC) sample of 10,000 mock stars, based on the probability density function, in each dimension of the input space. Then we connect each mock star with its BMU on the trained neuron map, and obtain the probability ($p$) of a mock star being associated to a detected group $\mathcal{G}$ of $n$ members. The significance of $\mathcal{G}$ can be quantified by the binomial probability of detecting more than $n$ stars from the halo sample of $\mathcal{N}$ stars. The contamination can be derived as $p$/($n$/$\mathcal{N}$). We consider $\mathcal{G}$ as valid only if the significance is larger than 5 $\sigma$ and the contamination rate is less than 20$\%$. For each DTG, the valid members are re-verified by their probabilities ($p\geqslant$ 20$\%$) of being associated to the same group, after taking the observational uncertainties into account. The confidence of each DTG is derived as the average probability of its valid member stars being associated to it.

After validation of the detected groups, we check if any valid group is associated to known MW GCs. This is done by generating 1000 realizations for each GC, according to its observational uncertainties in 6D kinematics. As done for the mapping of mock stars, we connect each realization of a given GC with its BMU on the neuron map. The confidence level of the association between a GC and a DTG is quantified by the probability of the mock GC sample being associated with the same DTG. This value is used to compare the associations of different GCs with their DTGs.

In this work, we focus on the substructures that are dynamically associated with MW GCs. Although numerous valid DTGs could be identified from the trained neuron map, only those having strong associations with GCs are analyzed in this work. In total, we identify four DTGs at three different values of $u_{\rm thr}$. The details of these groups are summarized in Table 1, where $n_{\rm BHB}$ and $n_{\rm RRL}$ denote the number of group members from the SDSS BHB and SDSS+LAMOST RRL samples, respectively. The contamination fraction, $\mathcal{F}_{\rm c}$, and confidence level for each DTG are listed, as well as its mean and dispersion of [Fe/H]. The assigned substructures and associated GCs are shown in the last two columns. Since our halo sample is mainly populated by stars with [Fe/H] $\lesssim-$1.5, all of the identified DTGs have very low mean metallicities ([Fe/H] $\approx$ $-$1.8 to $-$2.0). After taking into account the observational error of each group member star, the intrinsic dispersions of [Fe/H] of these DTGs are in the range of 0.25 – 0.65 dex, which excludes the possibility of their progenitors being GCs. Figure 1 clearly shows that all the DTGs stand out from the gray background of the halo sample, and form separate clumps in the input space of ($E$, $L$) and ($\theta$, $\phi$). We note that, except for DTG-1 (green), which has a retrograde orbit with positive $\theta$, the other three groups, DTG-2 (salmon), DTG-3 (magenta), and DTG-4 (blue) all have prograde orbits with negative $\theta$. The GCs associated with each DTG are embedded well-within its individual clump, shown as lime star symbols filled by the same colors as their DTGs. The confidence level of each association, derived in the same way as that of each member star, is provided in parentheses following the GC name in the last column of Table 1.

Figure 2 shows the location of the DTGs and their associated GCs in different dynamical-space visualizations, color-coded as in Figure 1. We differentiate the BHB members from the RRL members by filling the former with blue colors. The left panel of Figure 2 shows the projected action-space map. DTG-1 clearly occupies the corner of retrograde orbits, and DTG-4 is situated in the region representing radial orbits. DTG-2 and DTG-3 have fairly polar orbits, and significantly overlap with each other in this projection; note that DTG-2 has higher orbital energy than DTG-3, which makes them clearly separable in the ($E$, $L_{\rm z}$) space shown in the right panel. DTG-1 has slightly lower energy than DTG-2, but they have distinguishable distributions of $L_{\rm z}$. Among all the groups, DTG-4 has the lowest energy, as well as the lowest rotational motion. Utilizing these features of orbital properties, we analyze and assign an origin to each group below.

In this work, we employ two types of old stars, BHBs and RRLs, to construct a fair sample of $\sim$7600 stars in the ancient halo. Both are excellent distance indicators, thus we are able to obtain accurate 6D kinematic information for this halo sample, and derive their dynamical parameters. We apply the neural-network-based clustering method StarGO to this sample in the space of orbital energy and angular momentum, and identify four DTGs that are confidently associated with known MW GCs. The largest group, DTG-2, is confirmed to be the Sgr Stream, and the other two, DTG-1 and DTG-4, are very likely the debris of Sequoia and GES left in the outer halo.

We show that DTG-3 is a new stream, having a very polar orbit with inclination angle of 80^∘, which we refer to as LMS-1. It is associated with a pair of GCs (NGC 5024 and NGC 5053), which is the first example of a low-mass stellar-debris stream with embedded GCs from a disrupted low-mass dwarf galaxy. By tracing the orbit of the more-massive member of the pair, NGC 5024, we suggest that it is probably the core of the dwarf galaxy progenitor of LMS-1 and NGC 5053. We then use the orbital information of NGC 5024 as the initial condition for a dwarf satellite with a total mass of 2$\times$10⁹M_⊙, and run a N-body simulation in an analytic MW potential. The resulting stream nicely covers the observed LMS-1 in both hemispheres, with a few outliers possibly due to orbital precession.

The stellar system of LMS-1 and its pair of GCs belongs to the vast polar structure (VPOS; Riley & Strigari, 2020), initially suggested by Pawlowski et al. (2012). Although recent studies have shown that the orbital poles of the existing stellar streams and GCs do not cluster around the direction of the VPOS, the LMS-1 system adds a substantial accreted dwarf to it. We believe that numerous additional low-mass stellar-debris streams remain to be discovered that could be associated with GCs in the outer halo. This will help build a more complete MW assembly history, as well as open a new window to study the formation and evolution of GCs in ancient dwarf galaxies.

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

The authors would like to thank the referee for constructive comments. Z.Y. wishes to thank Eugene Vasiliev for insightful discussions about the GC pair, and Iulia Simion for sharing her expertise on VOD and HAC. Z.Y. acknowledges the support from the LAMOST Fellow project, and Shanghai Sailing Program (Y955051001). J.C. is supported by Young Scientists Foundation of JiangSu Province (BK20181110), the NSFC fundings (No. 11825303, 11861131006) and the Alibaba Cloud Fellowship. T.C.B. acknowledges partial support from grant PHY 14-30152, Physics Frontier Center/JINA Center for the Evolution of the Elements (JINA-CEE), awarded by the US National Science Foundation. T.C.B. is also grateful for support from a PIFI Distinguished Scientist award from the Chinese Academy of Sciences. Y.H. is supported by National Natural Science Foundation of China grants 11903027, 11973001, 11833006, 11811530289, U1731108, and U1531244, and National Key R&D Program of China No. 2019YFA0405503.