Revealing the Complicated Story of the Cetus Stream with StarGO

The Milky Way is expected to go through frequent merger events with nearby dwarf satellite galaxies during its assembly history. Each merger event will leave streams in the stellar halo, the kinematic information of which encodes the disruption history of their progenitor systems. According to the hierarchical structure formation paradigm, in general streams that stay in the outer stellar halo come from systems accreted later compared to those in the inner region (Bullock & Johnston 2005; Cooper et al. 2010; Amorisco 2017). Moreover, the effects of both tidal disruption and dynamical friction are milder as the Galactocentric distance increases. Thus, the streams in the outer stellar halo are less phase-mixed and the clustering features should be better preserved in phase space. However, the identification of these distant streams has previously been hindered by the lack of suitable proper motions. This situation is now changing thanks to the Gaia survey (Gaia Collaboration et al. 2016, 2018), which is delivering precise proper motions across the entire sky down to ∼21st magnitude. Accompanied by large spectroscopic surveys such as SDSS (York et al. 2000) and LAMOST (Cui et al. 2012; Deng et al. 2012; Zhao et al. 2012), we are now able to search for halo substructures using precise 6D phase space for the first time.

The most significant halo substructure found in the modern era is the disrupted Sagittarius dwarf galaxy (Ibata et al. 1994, 1995; Yanny et al. 2000; hereafter Sgr), and its associated streams (Mateo et al. 1996; Ibata et al. 2001; Majewski et al. 2003). Near the Sgr trailing arm there is another stream, named the Cetus Polar Stream (CS), which was first discovered by Newberg et al. (2009) using photometry and spectra from SDSS/SEGUE data release 7 (Abazajian et al. 2009). They also found that a distant halo globular cluster (GC), NGC 5824, is likely to be associated with CS, using an orbit deduced from stream's radial velocity. Yam et al. (2013), hereafter Yam13, updated the results of the distance and velocity measurement of CS using SDSS data release 8, including an N-body simulation to show that CS can be reproduced by a disrupted dwarf galaxy with a mass of 10⁸M_⊙.

Prompted by the association between NGC 5824 and the CS, a number of studies have investigated the nature of this GC to test whether it may be the disrupted core of the progenitor. Da Costa et al. (2014) found there is an intrinsic metallicity dispersion around 0.3 dex from low-resolution spectroscopy of Ca ii triplet absorption lines of 118 red giant members. However, a separate study using a restricted sample of 26 bright stars does not show a metallicity spread (Roederer et al. 2016). Instead, they found large internal dispersion of [Mg/Fe] of 0.28 dex, which is akin to massive metal-poor GCs such as ω Cen. The most recent study, using high-resolution spectra of 117 giant stars of NGC 5824 (Mucciarelli et al. 2018), confirms the lack of spread in [Fe/H] and the unusually large spread in [Mg/Fe]. It is also worth noting that deep photometric data from MegaCam and the Dark Energy Camera (Walker et al. 2017; Kuzma et al. 2018) find no evidence of tidal tails around NGC 5824, although the latter study does find that it has a very extended envelope with a radius of 230 pc. Such a feature may also be generated from the tidal stripping of a progenitor dwarf galaxy.

In this study we identify substructures in the outer stellar halo by crossing matching spectroscopic surveys with Gaia DR2. More specifically, we focus on the identified CS and its possible association with NGC 5824. The details of our catalog are described in Section 2. The methodology of the group identification and the map of the MW GCs is presented in Section 3. We show the identified streams and their associated GCs in Section 4. Conclusions and discussions are presented in Section 5.

For our analysis we combine two spectroscopic catalogs to create our parent sample. First, we use a catalog of K giants identified using a support vector machine classifier from LAMOST DR5 (Liu et al. 2014). We augment this using a catalog of BHBs from SDSS/SEGUE (Xue et al. 2008, 2011). We cross-match these stars with Gaia DR2 to obtain proper motions, resulting in a final parent sample of around 100,000 stars.

We focus on the distant halo in this work, where Gaia parallaxes are not suitable. Instead we obtain distance estimates from multi-band photometry (Xue et al. 2008, 2014) and, in the case of the K-giant stars, spectroscopic metallicity from LAMOST. The mean uncertainties in distances are about 13% and 5% for K-giant stars and BHB stars, respectively. Using the spectra of K giants, the LAMOST pipeline gives line-of-sight velocities with uncertainties of around 5–20 km s⁻¹ (Wu et al. 2011, 2014). There is a known systematic offset in the LAMOST radial velocity zero-point. This has been quantified for LAMOST DR5 by Tian et al. (2018), who compared to Gaia DR2 RVS velocities to obtain an offset of 5.4 km s⁻¹. We add this to our LAMOST radial velocities. Note that unlike previous versions of the LAMOST catalog, the DR5 pipeline has an improved algorithm for estimating uncertainties. This means that prior issues with overestimated uncertainties (e.g., Schönrich & Aumer 2017) have been alleviated and the uncertainties are now believed to be robust. For BHBs, line-of-sight velocities are derived from the SEGUE Stellar Parameter Pipeline (SSPP), which delivers uncertainties of 5–15 km s⁻¹ (Xue et al. 2008). Both pipelines provide measurements of [Fe/H] to an accuracy of around 0.1–0.2 dex.⁶ When we study the metallicity distributions of our identified halo groups, unless otherwise stated, we only include those members with S/N > 20.

Because we are interested in streams and overdensities that are not yet phase-mixed, we focus our study on the outer halo, retaining only stars with heliocentric distances greater than 15 kpc. In order to minimize the contribution from the Sagittarius stream we further reject stars with [Fe/H] > −1.5 dex. This cut should not remove a significant number of halo groups because we believe most accreted halo structures will come from dwarf galaxies, which typically have metallicities below −1.5 dex. The final sample has 5256 stars with 6D phase-space information and [Fe/H]. For our selected K giants, G-band magnitudes from Gaia DR2 (Lindegren et al. 2018) are in the range 14–18 mag, which have corresponding proper motion uncertainties in the range 0.04–0.28 mas yr⁻¹. Therefore, taking into account our distance uncertainties of around 15%, this means that typical transverse velocity uncertainties are around 10 km s⁻¹. The proper motion uncertainties of the final selected BHBs are slightly larger, at around 0.13–0.6 mas yr⁻¹, owing to their fainter magnitudes (G ∼17–19 mag). This yields typical transverse velocity uncertainties of 20 km s⁻¹. After applying the above cuts, the radial velocity uncertainties are around 5–10 km s⁻¹ for both populations. The uncertainty in each of the six phase-space dimensions is taken into account when we apply our group identification procedure (see Section 3.2 for more details).

Before going into the details of our approach, we briefly introduce our new cluster identification method StarGO (Yuan et al. 2018), which combines a self-organizing map (SOM) and a hierarchical group identification algorithm. The SOM trains the neural network to learn the data structure from an n-D input space and project it onto a 2D map. The clustering signatures within the n-D data set can be extracted from the visualization of the 2D neuron map. This is realized by introducing a weight vector for each neuron of the network, which has the same dimensions as the input vector. The unsupervised learning is performed by iteratively updating the weight vectors so that they become closer to the input vectors from the data set. This process is complete when all of the weight vectors reach convergence. The trained neural network is visualized by the difference in weight vectors among adjacent neurons, which is denoted by u. Neighboring neurons with u < u_thr are identified as groups, where u_thr is determined based on the distribution of u, as well as the topological structure exhibited by the neuron map (see Section 3.2 for a demonstration of how this works in practice). After the group identification procedure, we associate each star to its best matching unit (BMU) on the 2D map, which is defined to be the neuron with the weight vector closest to the input vector of that star. Stars associated with the identified neuron groups form the corresponding star groups. The performance of StarGO in relation to the commonly used Friends-of-Friends algorithm has been quantified in Yuan et al. (2018).

3.1. Input Space

3.2. Group Identification

We feed the input data into a 150 × 150 neural network, using 400 iterations of the learning process. Because we use two angles to represent the direction of the angular momentum, we use a great-circle distance to calculate the difference between neurons. Therefore, in our final trained map the difference between neighboring neurons i and j (hereafter u_i,j) is defined by the following equations:

$$\begin{eqnarray}&&{I}_{i,j}=\sin {\left(\delta {\theta }_{i,j}/2\right)}^{2}+\cos {\theta }_{i}\cos {\theta }_{j}\sin {\left(\delta {\phi }_{i,j}/2\right)}^{2},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}\begin{array}{rcl}{u}_{i,j} & = & {\mathrm{log}}_{10}\left[{\left(\delta {E}_{i,j}/{E}_{\mathrm{norm}}\right)}^{2}+{\left(\delta {L}_{i,j}/{L}_{\mathrm{norm}}\right)}^{2}\right.\\ & & +\left.{\left(2/\pi \cdot \arcsin \sqrt{{I}_{i,j}}\right)}^{2}\right].\end{array}\end{eqnarray} \tag{ 2 }$$

The trained neural network is visualized by the grayscale map of u (see Figure 1(b)), where u represents the difference in weight vectors between adjacent neurons. The clouds in grayscale correspond to distinctive structures within the underlying data set, with darker shades denoting higher values of u, i.e., steeper gradients in the weight vectors.

Zoom In Zoom Out Reset image size

We first mark neurons with u < u_5% as light blue pixels and neurons with u < u_16% as gold pixels in Figure 1(c), where u_5% and u_16% denote the 5th and 16th percentiles of the distribution of u, respectively (see Figure 1(a)).

The light blue patches of neurons are very similar (among the top 5% of all neurons) in the weight vector space. They are connected by the gold neurons, which are the top 16% of all the neurons. These highlighted neurons reveal the clustering structures from the input space. The inner dense cores of clusters are associated with the light blue neurons, and the less dense outskirts are associated with the gold neurons. In the trained map we can easily see three regions that are populated by significant light blue patches connected by gold neurons: one in the middle, one in the lower left corner, and one in the lower right corner. We consider these groups as candidate groups, and try to identify their group members systematically.

We then mark the neurons with the lowest clustering significance as purple pixels, where we use a threshold u > u_thr. This threshold parameter u_thr is adjusted in order to maximize the size of each group, as described below. For Figure 1(c) we have chosen u_thr = u_60%, i.e., the purple pixels correspond to neurons with clustering significance in the bottom 40%. These pixels form the boundaries separating candidate neuron groups into isolated islands, such as the islands corresponding to the salmon and cyan candidate groups in Figures 3(Ia) and (IIa). Neurons within the islands have similar weight vectors and these weights are significantly different from the neurons in the purple boundary regions. These islands represent the strongest clustering signals in the underlying 4D data set.

In order to systematically identify clustered star groups, we start with u_thr0 = u_100%. At this stage our neuron map consists of one large group. As we gradually decrease the value of u_thr, boundary (purple) pixels begin to populate the map and separate islands begin to emerge (analogous to the rising tide cutting off Mont-Saint-Michel from the mainland). The value of u for which the first significant island emerges is denoted as u_thr1, which in this case corresponds to u_66%. This produces the salmon island in Figure 3(Ia). We then continue decreasing u_thr until the next large island forms, which is at u_thr2 = u_60% (corresponding to the cyan island in Figure 3(IIa)). As we keep decreasing u_thr the island in the lower right corner appears when u_thr3 = u_30% (the magenta island in Figure 3(IIIa)). At this threshold, the region populated by light blue patches in the middle of the map is still connected with the gold patches in the upper right corner (Figure 3(IIIa)). We could continue the above prescription to isolate this candidate group and other, smaller groups. We defer the analysis of these new groups to a subsequent paper.

Observational uncertainties are taken into account by utilizing the 2D neuron map. For each star we generate 1000 Monte Carlo realizations, given each star's uncertainty in distance, radial velocity, and proper motion. We assume all uncertainties follow Gaussian distributions. After the creation of these 1000 realizations, for each one we identify its BMU in the trained map. Note that the BMU for each realization is not necessarily the same as the BMU of the fiducial star.

In principle, any star can be mapped to the trained neural network. However, the difference in the input vector and the weight vector of its BMU can be very large if the input vector is very different from the training sample. Mapping such stars has little meaning because they are not closely associated with any neurons. Thus, we need to define the limit of the trained neural network, within which the mapping can be considered appropriate. We first calculate the 4D distance between stars in the training sample and their BMUs, which is referred to as u_vw. Then we can get the largest distance denoted by u_vw,max out of the entire sample. This maximum offset denotes the limits of the neuron map, i.e., stars that lie beyond this limit cannot be closely attached to any neuron of the map. We remove such data points from our set of realizations if the 4D distance to its BMU is larger than u_vw,max.

We also wish to test whether the star's association to its respective group is robust. After rejecting realizations according to the preceding paragraph, we then analyze the BMU for each remaining realization. For illustrative purposes, we plot the 2D density map of all the remaining realizations for the magenta group in Figure 2(a). It can be seen that some BMUs lie outside the boundary of the identified group and we consider these invalid realizations.

Zoom In Zoom Out Reset image size

After we have carried out these two steps and obtained the total number of valid realizations for each of the candidate members, we reject the worst 16%, where "worst" corresponds to the candidate members with the fewest valid realizations. This sets the confidence level of member identification for each group individually. For OS the candidate members all have at least 200 valid realizations, which means more than 200 out of 1000 realizations for each member in OS are associated with the corresponding neuron group. For CS and Sgr the candidate members all have at least 400 valid realizations out of 1000.

3.3. Mapping the GCs of the Milky Way

The first iteration of StarGO generates the 2D neuron map shown in Figure 1, which was introduced in Section 3.2. As discussed in Section 3.2, our analysis of this map results in three identified groups, which we will now discuss. These are illustrated in Figure 3, where each row corresponds to one of the identified groups.

Zoom In Zoom Out Reset image size

4.1. The Orphan and Sagittarius Streams

In this subsection we demonstrate the efficacy of our method by presenting two groups that correspond to previously identified halo streams.

The first group identified in Section 3.2 is shown in the first row of Figure 3. It is the largest identified indivisible group in our halo catalog, consisting of 390 K giants and 132 BHBs. From its sky position (Figure 4) we can easily see that this group corresponds to the Sagittarius stream (Sgr; e.g., Majewski et al. 2003; Belokurov et al. 2006). As discussed in Section 3.3, we have mapped a sample of GCs onto the same neural network. We find that five of these GCs are closely associated with our Sagittarius group (the red stars in Figure 3). These are Whiting 1, NGC 6715 (M54), Arp 2, Terzan 8, and Pal 12, all of which have been previously associated with the Sagittarius stream (Law & Majewski 2010a, 2010b; Sohn et al. 2018).

Zoom In Zoom Out Reset image size

Next, we discuss the group in the lower left corner of our neuron map, shown in the second row of Figure 3. Although there are only 22 stars belonging to this group, it is compact in both energy and angular momentum spaces. The sky positions of these stars are shown in Figure 4, from which it is clear that this group corresponds to the Orphan Stream (OS; e.g., Belokurov et al. 2006). The distances to these stars range from 18 to 47 kpc, with a clear gradient along the stream, e.g., similar to the results from SEGUE (19–47 kpc; Newberg et al. 2010). The OS is also detected in the southern sky to decl. = −38° (Grillmair et al. 2015). A recent study by Koposov et al. (2019) shows an all-sky view of the OS, which has a ∼210° length. However, we can only identify its stream members with decl. >−20°, as our halo catalog consists of data from northern hemisphere surveys (SDSS/SEGUE and LAMOST).

4.2. The Cetus Stream and NGC 5824

The last of our three groups is shown in the third row of Figure 3. This contains 151 members from our halo catalog, and as with the other groups, it is highly clustered in the input space (see Figures 3(IIb) and (IIc)). The sky position of these members shows two separated regions, one in the southern Galactic hemisphere (at 0° ≳ R.A. ≳ 50°, corresponding to l ∼ 150°, b ∼ −40°) and one in the north (at −150 ≳ R.A. ≳ −100, corresponding to l ∼ 50°, b ∼ 50°). The former region is the most dense, with 106 stars. The distance (20–35 kpc) and sky position tells us that this corresponds to the Cetus stream (CS; e.g., Newberg et al. 2009; Yam et al. 2013). Although it overlaps with part of the Sgr stream in sky position (see Figure 4) and has a similar distance, the directions of the angular momentum of Sgr and CS are clearly different (Figure 3(Ic) and (IIIc)).

While the detection in the southern Galactic hemisphere is not new, the detection of the northern counterpart has not previously been found. This new feature consists of 45 stars and is more diffuse on the sky. Because both counterparts come from the same single identified group in the neuron map, sharing the same total energy and angular momentum, we believe these two overdensities are from the same stream. We hereafter refer to these as the northern and southern counterparts of the CS.

4.2.1. Statistical Significance of the Northern Detection

Before we discuss the association of these two detections, we first test the significance of the new northern counterpart. This is done using a smooth halo drawn from the mock Gaia DR2 catalog of Rybizki et al. (2018), which is based on the popular Galaxia/Besancon models (Robin et al. 2003; Sharma et al. 2011). We create a sample of mock halo stars similar to our observed sample of metal-poor stars, taking halo stars (i.e., age ≥11 Gyr) with decl. ≳ −10°, d ≳ 15 kpc, [Fe/H] < −1.5 dex, and G band magnitude < 19. We then feed these mock stars through our neural network and select those that are tightly associated with the CS detection (i.e., the magenta region in Figure 3(IIIa)). We find that the smooth halo produces 3301/155,343 = 2% stars associated with the CS detection, which is far fewer than the observed fraction of 151/1679 = 9%. If we concern ourselves with only the new northern detection (restricting ourselves to b ≥ 20°), we find that the mock sample has 1795/107,486 = 1.6% of stars associated with this detection, whereas the observed sample contains 45/1105 = 4%. In order to put this on a more quantitative footing, we can calculate the corresponding binomial probability, i.e., the probability that a smooth halo can, by chance, provide 45 stars that lie within our CS detection region. This takes into account the Poisson fluctuation in the background population, under the assumption that the sample of 1105 halo stars is drawn from a smooth population. We find that the binomial probability of obtaining 45 stars in our detection region, from a sample of 1105 smooth stars, is less than 1e-7%, i.e., the significance is more than 5σ.

4.2.2. Metallicity

In Figure 5 we plot the metallicities of our CS members. We divide our sample into K giants (g-band S/N >20), and BHBs, for both northern and southern detections. The metallicities of the K giants confirm our claim that these two detections are part of the same stream. Both northern and southern detections have very similar metallicities (see Figure 5(a)), and crucially these metallicities are clearly distinct from the background halo population (see Figure 5(c)). If we apply a KS test to assess the probability that the metallicities of the northern and southern detections are consistent with the background halo population, we obtain p-values below 0.01, i.e., they are clearly inconsistent. Note that the K-giant metallicities have been calculated using the LAMOST metallicity pipeline, which only extends to around −2.5 dex (Luo et al. 2015), so the lack of a low-metallicity tail compared to the BHBs is not surprising.

Zoom In Zoom Out Reset image size

We measure the mean and dispersion of the metallicity for these stars using a Bayesian maximum likelihood technique (Foreman-Mackey et al. 2013), taking into account the observational uncertainties of around 0.2 dex. We find that the BHBs have a mean of −2.15 ± 0.04 dex and a dispersion of 0.23 ± 0.04 dex, while the K giants have a mean of −2.07 ± 0.02 dex and a dispersion of 0.12 ± 0.02 dex. The means are consistent, and although the dispersions of the two populations differ slightly, it is clear that these dispersions are inconsistent with zero at a high level of significance. We therefore conclude that the stream comes from a disrupted dwarf galaxy rather than a disrupted GC, for which we would expect a dispersion of ∼0 dex. We return to this issue in the discussion.

4.2.3. Association to NGC 5824

As before, we check to see whether any MW GCs could be associated with this group and find that NGC 5824 (the red star in Figure 3(III)) is tightly associated. As shown in Figure 2(b) and discussed in Section 3.3, our Monte Carlo analysis demonstrates that the association is robust. We defer a discussion of the orbital properties of the cluster to the following subsection.

We also compare the metallicity of CS members with high-resolution studies of NGC 5824. This is shown in Figures 5(c) and (d). Roederer et al. (2016) found [Fe/H] = −1.94 ± 0.12 dex from 26 red giants stars (see the gold dashed line with a transparent band). In a later study, Mucciarelli et al. (2018) found [Fe/H] = −2.11 ± 0.01 dex (right lime dashed line) from 87 red giant branch stars, and [Fe/H] = −2.20 ± 0.01 dex (left lime dashed line) from 30 asymptotic giant branch stars. Overall, the typical metallicity of CS members agrees very well with that of NGC 5824 from these high-resolution spectroscopic studies. This fact supports the scenario that they are from the same progenitor system. We return to this issue later in the discussion.

4.2.4. Orbital Properties

We now investigate the orbit of the stream. Hitherto the orbit has not been well-determined, owing to the lack of precise proper motions (e.g., Yam13), but Gaia DR2 enables us to analyze the orbit in much greater detail. Figure 6 shows the properties of our stars in terms of heliocentric distance and radial velocity as a function of Galactic latitude, which is chosen as the proxy for angle along the stream because the orbit has a large extension in Galactic latitude.

Zoom In Zoom Out Reset image size

We start by analyzing the stars from the existing detection in the south (−70° ≲ b ≲ − 20°). We can see from the Galactic rest-frame radial velocity distribution (v_gsr) that our detection appears to be split into two clumps, one with −75 < v_gsr < −25 km s⁻¹ (cyan circles in Figure 6) and one with 50 < v_gsr < 100 km s⁻¹ (red circles in Figure 6). The former clump corresponds to the detection in Yam13, which has a negative radial velocity gradient with b. The clump at positive v_gsr is new, although it should be noted that Yam13 found a number of stars with velocities that are consistent with this clump. Due to small number statistics they were unable to make the association, but with our powerful StarGO tool and a much larger catalog from LAMOST we are now able to do so. What we appear to be detecting is phase-mixed material from the progenitor system.

To understand this better, we compare the properties of the stream to the orbit of NGC 5824, which has pericenter and apocenter at 19 kpc and 41 kpc, respectively. It is possible that they are stripped from the same progenitor system. In this scenario, even though tidal streams do not necessarily follow the precise orbit of the progenitor, it should still enable us to gain a clearer picture of its disruption. Figure 6 shows the orbit of the cluster integrated forward and backward for 1.8 Gyr, which is around three times its orbital period (T = 0.6 Gyr). The thin lines denote parts of the orbit for which decl. < −10°, i.e., these cannot correspond to our debris, as our stars are all located at decl. > −10°.

We first focus on the existing detection from Yam13, given by the cyan circles in Figure 6. The grayscale maps in the middle columns show the orbits of these cyan stars integrated for a short time in order to illustrate their gradients. From these, there are two portions of the NGC 5824 orbit that are consistent with this debris, highlighted by the transparent gold and pink bands. The former is around 1.2 Gyr ahead of the cluster (∼2 orbital periods), while the latter is around 1.5 Gyr behind (∼2.5 orbital periods). The new detection at positive radial velocities (red points) is shown in the right columns. Here there is only one portion of the orbit that is consistent with the debris, highlighted with the transparent blue band. This is around 0.3 Gyr ahead of the cluster (∼0.5 orbital periods). The new detection in the northern Galactic hemisphere (magenta points) is harder to discern due to the wide spread in radial velocities, but is consistent with a number of portions of the orbit. Note that even though the radial velocities cover a wide range, their angular momenta are all tightly concentrated (see Figure 3).

If the above picture is correct, it peculiar that there are gaps in the detecting material. Part of this can be explained by portions of the orbit that lie outside our field of view (either at negative decl. or behind the Galactic plane), but there are other portions that we could be expected to detect. In particular the trailing material from less than one orbital period behind the cluster should be visible at −70° ≲ b ≲ −20°. This portion of the stream corresponds to pericenter of the orbit and so perhaps the stream is sparsely populated here.

From this orbital analysis, together with the homogeneous metallicities shown in Figure 5, we conclude that these new detections (both at positive v_gsr and in the northern hemisphere) have orbits that are consistent with the CS from Yam13 as well as NGC 5824. It is very likely that all of this material comes from one disrupted system. The system leaves material at a wide range of orbital phases (about 4.5 orbital periods), and as a consequence the stream is very diffuse.

It is possible that this system could consist of multiple components that fell in together with very similar orbital properties. A plausible picture is that NGC 5824 may come from another dwarf galaxy, which fell into the MW along with the progenitor dwarf of the CS. In this case, the different portions of the stream may be stripped from these two galaxies. Group infall is the most common scenario for the accretion of satellite galaxies in zoom-in simulations of MW-like systems (Wetzel et al. 2015). There is also observational evidence supporting this picture. For example, the Leo-Crater group (de Jong et al. 2010; Belokurov et al. 2014; Torrealba et al. 2016) contains several low-mass dwarf galaxies (Leo V, Leo IV, and Crater 2) and Crater Globular Cluster. Recently, Kallivayalil et al. (2018) showed that four newly discovered ultra-faint dwarf galaxies may have fallen in with the Magellanic Clouds.

In this work we have applied a novel clustering method StarGO to a metal-poor halo catalog with full 6D phase-space information, constructed by cross-matching Gaia DR2 with spectroscopy from LAMOST DR5 and SDSS/SEGUE. After the first iteration of the workflow we are able to identify three significant star groups, two of which can be easily confirmed as the Sagittarius and Orphan streams. The third star group corresponds to the CS. Our detection populates two separate regions on the sky, one in the southern Galactic hemisphere and one in the north. The southern counterpart (at 0° ≳ R.A. ≳ 50°) corresponds to the previous CS from Newberg et al. (2009) and Yam13. The northern counterpart (at −150° ≳ R.A. ≳ −100°) is a new discovery and has the same angular momentum and orbital energy as the southern counterpart.

In order to understand the progenitor system more, we have investigated the orbit of the CS in detail, which is now possible thanks to the precise proper motion information from Gaia. Despite often being referred to as the Cetus Polar Stream, we have found that the angular momentum vector of our CS members is not actually aligned with the Galactic plane, as expected for a polar stream, and is offset by around 45° (see Figure 3(IIIc)). We therefore suggest future studies drop the term "polar" from its name. We also found the southern counterpart of the CS bifurcates in radial velocity. The clump with negative radial velocity is close to the Yam13 detection, which has a negative distance gradient with b. The new clump has a positive distance gradient with b and, when compared to the orbit of the existing material, aligns well albeit with a phase offset of around 0.9 Gyr (corresponding to ∼1.5 times the orbital period).

We believe that the new northern counterpart is a mixture of debris. In general the stream is very diffuse, as had been seen in previous studies. That means that it is hard to identify in configuration and/or velocity space. The high degree of phase mixing suggests the progenitor was massive and/or was stripped a long time ago (Johnston et al. 2001; Amorisco 2015). A suite of N-body simulations is needed in order to test these possibilities and draw quantitative conclusions about the progenitor and its disruption history.

We have checked the metallicity of the southern and northern parts and find that they are in very good agreement, noticeably offset from the background halo population. The intrinsic metallicity dispersion (i.e., after accounting for observational errors) of all the CS members from the BHB catalog is 0.23 ± 0.04 dex, which is significantly larger than the spread expected for a GC and more consistent with a dwarf galaxy. According to the stellar mass–stellar metallicity relation for Local Group dwarf galaxies from Kirby et al. (2013), the metallicity of ∼−2 dex suggests a stellar mass in the range 10⁴M_⊙–10⁶M_⊙, i.e., similar to low-mass classical MW dwarf galaxies such as Ursa Major I or Ursa Minor. The identified CS has a remarkable high ratio of BHB to K-giant stars (81/70), as has been noted previously by Newberg et al. (2009), Koposov et al. (2012), and Yam et al. (2013). The high fraction of BHBs suggests the progenitor system had most of its stars form in an early period of star formation.

By exploiting the visualization power of StarGO, we are able to map MW GCs to the trained neural network. We find a strong dynamical association between the Sgr stream and five of its surrounding GCs, which has been reported in previous studies (Law & Majewski 2010a; Sohn et al. 2018) using both HST and Gaia DR2 data. For the CS, its possible relationship between NGC 5824 was initially postulated using radial velocities in its discovery paper (Newberg et al. 2009). Thanks to Gaia we can now obtain 6D kinematics of both the CS and NGC 5824. From its angular momentum and orbital energy we find that NGC 5824 is located within the same group as the CS members, confirming the association between the cluster and the stream. By integrating the orbit of NGC 5824 we find that the stream is consistent with material stripped.

The metallicities of our CS detection, both for K giants and BHBs, are in good agreement with that of NGC 5824 from spectroscopic studies of the cluster (Da Costa et al. 2014; Roederer et al. 2016; Mucciarelli et al. 2018), providing independent evidence of their common origin. Another possible scenario is that NGC 5824 comes from a different dwarf galaxy progenitor that shares a similar orbit and mass with the progenitor of the CS. This happens when a group of satellites fall into a large host galaxy such as the MW. Such a picture is supported by the group infall scenario from both observations and simulations.

NGC 5824 is also reported to have a small spread in [Fe/H] but a large spread in [Mg/Fe] (Roederer et al. 2016; Mucciarelli et al. 2018), which is typically found only in massive and/or metal-poor clusters such as ω Cen. The mass–luminosity relation for GCs suggests a stellar mass of over 10⁶ M_⊙. This is comparable to, or even greater than, the stellar mass of the progenitor dwarf estimated above (10⁴–10⁶ M_⊙). High-resolution hydrodynamical simulations have shown such systems can be robustly produced, with star clusters forming at the centers of some halos at high redshift (Kimm et al. 2016; Ricotti et al. 2016). Given the mass ratio between NGC 5824 and its host galaxy, if the star cluster is initially born off-center, it will inspiral into the center on a short timescale (∼1 Gyr) due to dynamical friction (Binney & Tremaine 2008). According to quantitative comparisons between the MW GC population and cosmological zoom-in simulations from the E-MOSAICS project, the MW hosts 6 ± 1 nuclear clusters from disrupted dwarf galaxies (Kruijssen et al. 2019). NGC 5824, like M54 and wCen, is likely to be the nucleus of an accreted dwarf galaxy.

Clearly this is a complex system. The fact that it appears to be almost fully phase-mixed means that orbital analyses (such as those we have presented) are always going to have limited power. Detailed N-body simulations are required to fully understand the history of the progenitor's accretion and its relation to NGC 5824. We have begun to undertake this task and will report our findings in due course.

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 wish to thank Hao Tian and Jundan Nie for sharing the Match Filter code when checking the tidal tails of NGC 5824, and João A. S. Amarante for providing an introduction to emcee. Z.Y. and M.C.S are partly supported by the National Key Basic Research and Development Program of China (No. 2018YFA0404501) and NSFC grant 11673083. Z.Y. and J.L. acknowledge the support from Special Funding for Advanced Users through LAMOST FELLOWSHIP. J.L. also acknowledges the NSFC grant 11703019. X.-X. Xue is thankful for the support of the "Recruitment Program of Global Youth Experts" of China and NSFC under grants 11390371, 11873052, and 11890694. This project was developed in part at the 2018 Gaia-LAMOST Sprint workshop, supported by the National Natural Science Foundation of China (NSFC) under grants 11333003 and 11390372. Z.Y. and Y.W. are partly supported by the China Postdoctoral Science Foundation. Y.W. also acknowledges the Young Researcher Grant of National Astronomical Observatories, Chinese Academy of Sciences; and the National Natural Science Foundation of China under grant 11803048. J.C. is supported by the Astronomical Big Data Joint Research Center, cofounded by the National Astronomical Observatories, Chinese Academy of Sciences and the Alibaba Cloud.