Revealing the Complicated Story of the Cetus Stream with StarGO

TL;DR

This study employs a novel unsupervised learning tool, StarGO, to identify and analyze stellar streams in the Milky Way's outer halo, revealing new features and associations with globular clusters, and suggesting a dwarf galaxy merger origin.

Contribution

The paper introduces StarGO, an innovative neural network-based method for identifying stellar streams and their associated clusters in the Milky Way's halo using multi-dimensional data.

Findings

Discovery of a northern counterpart to the Cetus stream.

Confirmation of the association between NGC 5824 and the Cetus stream.

Identification of a bifurcated southern Cetus stream with different velocities.

Abstract

We use a novel cluster identification tool StarGO to explore the metal poor ([Fe/H] $<$ -1.5) outer stellar halo (d $>$ 15 kpc) of the Milky Way using data from Gaia, LAMOST and SDSS. Our method is built using an unsupervised learning algorithm, a self-organizing map, which trains a 2-D neural network to learn the topological structures of a data set from an n-D input space. Using a 4-D space of angular momentum and orbital energy, we identify three distinct groups corresponding to the Sagittarius, Orphan, and Cetus Streams. For the first time we are able to discover a northern counterpart to the Cetus stream. We test the robustness of this new detection using mock data and find that the significance is more than 5-sigma. We also find that the existing southern counterpart bifurcates into two clumps with different radial velocities. By exploiting the visualization power of StarGO, we attach MW globular clusters to the same trained neural network. The Sagittarius stream is found to have five related clusters, confirming recent literature studies, and the Cetus stream has one associated cluster, NGC 5824. This latter association has previously been postulated, but can only now be truly confirmed thanks to the high-precision Gaia proper motions and large numbers of stellar spectra from LAMOST. The large metallicity dispersion of the stream indicates that the progenitor cannot be a globular cluster. Given the mean metallicity of the stream, we propose that the stream is the result of a merger of a low-mass dwarf galaxy that hosted a large nuclear star cluster (NGC 5824).

Equations

$$I_{i,j}=\sin (\delta {\theta }_{i,j}/2{)}^2+\cos {\theta }_i\cos {\theta }_j\sin (\delta {\varphi }_{i,j}/2{)}^2,$$

$$\begin{array}{rl}{u}_{i,j}& ={\log }_{10}[(\delta E_{i,j}/E_{\mathrm{norm}}{)}^2+(\delta L_{i,j}/L_{\mathrm{norm}}{)}^2\\ & +(2/\pi \cdot \arcsin \sqrt{I_{i,j}}{)}^2].\end{array}$$

Full text

1. Introduction

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 Galacto-centric 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 $\sim$21st magnitude. Accompanied by large spectroscopic surveys such as SDSS (York et al., 2000) and LAMOST (Cui et al., 2012; Zhao et al., 2012; Deng et al., 2012), we are now able to search for halo substructures using precise 6-D 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 N-body simulation to show that CS can be reproduced by a disrupted dwarf galaxy with mass of $10^{8}$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) find 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 find large internal dispersion of [Mg/Fe] of 0.28 dex, which is akin to massive metal poor GCs such as $\omega$ 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 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 the 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 Sec. 2. The methodology of the group identification and the map of the MW GCs is presented in Sec. 3. We show the identified streams and their associated GCs in Sec. 4. The conclusions and discussions are presented in Sec. 5.

2. Data

For our analysis we combine two spectroscopic catalogs to create our parent sample. Firstly 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.

Since we focus on the distant halo in this work, 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*-1* (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*-1*. 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 km s*-1* to 15 km s*-1* (Xue et al., 2008). Both pipelines provide measurements of [Fe/H] to an accuracy of around 0.1 to 0.2 dex111Note that the BHB metallicities in Xue et al. (2008) were obtained using a bespoke version of the SSPP..When we study the metallicity distributions of our identified halo groups, unless otherwise stated, we only include those members with S/N $>$ 20.

Since 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 helio-centric 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 5,256 stars with 6-D 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*-1*. Therefore, taking into account our distance uncertainties of around 15%, this means that typical transverse velocity uncertainties are around 10 km s*-1*. The proper motion uncertainties of the final selected BHBs are slightly larger, at around 0.13 – 0.6 mas yr*-1*, owing to their fainter magnitudes ($G$ $\sim$ 17 – 19 mag). This yields typical transverse velocity uncertainties of 20 km s*-1*. After applying the above cuts, the radial velocity uncertainties are around 5 to 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 details in Sec. 3.2).

3. Method

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 dimension 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_{\mathrm{thr}}$ are identified as groups, where $u_{\mathrm{thr}}$ is determined based on the distribution of $u$, as well as the topological structure exhibited by the neuron map (see Sec. 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

In this work we apply StarGO to the metal poor halo catalog using an input space constructed from angular momentum and orbital energy. More specifically, we use total angular momentum $L$, orbital energy $E$, and two parameters, $\theta$ and $\phi$, to characterize the direction of the angular momentum vector. Our goal is to find star groups which are clustered in the ($E$, $L$, $\theta$, $\phi$) space. As we know, in an axisymmetric potential, $E$, $L$, and $L_{z}$ are approximately conserved and the latter two give $\theta=\mathrm{arccos}(L_{z}/L)$. Although $\phi=\mathrm{arctan}(L_{x}/L_{y})$ is not strictly conserved, it has been shown that $L_{x}$ and $L_{y}$ evolve coherently (e.g. Helmi & de Zeeuw, 2000; Knebe et al., 2005; Klement, 2010; Gómez & Helmi, 2010; Maffione et al., 2015) and so we take this angle as the fourth dimension of our input space. For the halo stars of extra-galactic origin, the exchange in angular momentum and orbital energy can be severe when their progenitors fall into the MW. However, stars from the same progenitor will have a similar disruption history, leading to similar changes in angular momentum and orbital energy. Although the conservation of angular momentum and orbital energy will be violated, stars from a given progenitor can remain clustered in this space, as has been shown in mock simulations of substructures stripped from merging satellites (Yuan et al., 2018). The final step is to make these quantities dimensionless, i.e. $(E-E_{\rm{min}})/E_{\mathrm{norm}}$, $(L-L_{\rm{min}})/L_{\mathrm{norm}}$, $\theta/\pi$, and $\phi/2\pi$, where $E_{\mathrm{norm}}$ and $L_{\mathrm{norm}}$ denote the ranges of the entire sample (e.g. $E_{\rm{norm}}=E_{\rm{max}}-E_{\rm{min}})$.

In a following paper, Xue et al. (in prep) will present a companion analysis using a Friends-of-Friends approach, based on a similar input space.

3.2. Group Identification

We feed the input data into a 150$\times$150 neural network, using 400 iterations of the learning process. Since 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 neighbouring neurons $i$ and $j$ (hereafter $u_{i,j}$) is defined by the following equations,

[TABLE]

[TABLE]

The trained neural network is visualized by the gray-scale map of $u$ (see Fig. 1b), where $u$ represents the difference in weight vectors between adjacent neurons. The clouds in the gray-scale 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.

We first mark neurons with $u<u_{5\%}$ as light blue pixels and neurons with $u<u_{16\%}$ as gold pixels in Fig. 1c, where $u_{5\%}$ and $u_{16\%}$ denote the 5th and 16th percentile of the distribution of $u$ respectively (see Fig. 1a).

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_{\mathrm{thr}}$. This threshold parameter $u_{\mathrm{thr}}$ is adjusted in order to maximize the size of each group, as described below. For Fig. 1c we have chosen $u_{\mathrm{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 Fig. 3Ia, 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 4-D data set.

In order to systematically identify clustered star groups, we start with $u_{\mathrm{thr0}}=u_{100\%}$. At this stage our neuron map consists of one large group. As we gradually decrease the value of $u_{\mathrm{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_{\mathrm{thr1}}$, which in this case corresponds to $u_{66\%}$. This produces the salmon island in Fig. 3Ia. We then continue decreasing $u_{\mathrm{thr}}$ until the next large island forms, which is at $u_{\mathrm{thr2}}=u_{60\%}$ (corresponding to the cyan island in Fig. 3IIa). As we keep decreasing $u_{\mathrm{thr}}$ the island in the lower right corner appears when $u_{\mathrm{thr3}}=u_{30\%}$ (the magenta island in Fig. 3IIIa). 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 (Fig. 3IIIa). 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 2-D 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 since they are not closely associated to any neurons. Thus we need to define the limit of the trained neural network, within which the mapping can be considered as appropriate. We first calculate the 4-D distance between stars in the training sample and their BMUs, which is referred to as uvw. Then we can get the largest distance denoted by uvw,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 4-D distance to its BMU is larger than uvw,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 2-D density map of all the remaining realizations for the magenta group in Fig. 2a. It can be seen that some BMUs lie outside the boundary of the identified group and we consider these as invalid realizations.

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 globular clusters of the Milky Way

In order to study the relationship between the MW GCs and the identified star groups, we attach each GC to their BMU on the trained neuron map This process requires the kinematics of each GC in order to calculate its angular momentum and orbital energy. These kinematics are taken from the Harris catalog (Harris, 1996, 2010) and Vasiliev (2018). We normalize the input vector of each GC by $E_{\mathrm{norm}}$ and $L_{\mathrm{norm}}$ from the training sample (see Section 3.1). We also exclude all GCs for which the 4-D distance to the BMU is larger than the largest offset between a star in our training sample and its BMU. The GCs that pass these cuts are mapped to their BMUs and denoted by red stars in Fig. 1c, where is can be seen that many lie clustered in this space and are associated to over-densities in the input stellar halo catalog. We plot the GCs attached to the magenta and salmon islands in Fig. 3I, III as red stars. Following the same recipe as the previous section, for each GC that is associated to one of the groups we generate 1000 realizations to test for robustness. We find that each cluster has at least 900 valid realizations. We illustrate this process in Fig. 2b for the GC lying in the magenta island, showing the 2-D density map of all realizations that can be closely attached to the neuron network.

4. Results

The first iteration of StarGO generates the 2D neuron map shown in Fig. 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 Fig. 3, where each row corresponds to one of the identified groups.

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 Fig. 3. It is the largest identified indivisible group in our halo catalog, consisting of 390 K giants and 132 BHBs. From its sky position (Fig. 4) we can easily see that this group corresponds to the Sagittarius stream (Sgr; e.g. Majewski et al., 2003; Belokurov et al., 2006a). 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 to our Sagittarius group (the red stars in Fig. 3). These are Whiting 1, NGC 6715 (M 54), Arp 2, Terzan 8, and Pal 12, all of which have been previously associated to the Sagittarius stream (Law & Majewski, 2010b, a; Sohn et al., 2018).

Next we discuss the group in the lower-left corner of our neuron map, shown in the second row of Fig. 3. Although there are only 22 stars belonging to this group, it is compact in both energy and angular momentum spaces. The sky position of these stars are shown in Fig. 4, from which it is clear that this group corresponds to the Orphan Stream (OS; e.g. Belokurov et al., 2006b). The distances to these stars range from 18 kpc to 47 kpc, with a clear gradient along the stream, i.e. 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 deg (Grillmair et al., 2015). Recent study by Koposov et al. (2019) shows an all-sky view of the OS, which has $\sim$ 210 deg length. However we can only identify its stream members with Decl. $>$ -20 deg as our halo catalog consists of data from northern hemisphere surveys (SDSS/SEGUE and LAMOST).

4.2. The Cetus Stream and NGC 5824

The final of our three groups is shown in the third row of Fig. 3. This contains 151 members from our halo catalog and, as with the other groups, it is highly clustered in the input space (see Fig. 3IIb,IIc). The sky position of these members shows two separated regions, one in the southern Galactic hemisphere (at 0 $\gtrsim$ R.A. $\gtrsim$ 50 deg, corresponding to $l\sim 150$ deg, $b\sim-40$ deg) and one in the north (at -150 $\gtrsim$ R.A. $\gtrsim$ -100, corresponding to $l\sim 50$ deg, $b\sim 50$ deg). 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 Fig. 4) and has a similar distance, the directions of the angular momentum of Sgr and CS are clearly different (Fig. 3Ic,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 (Sharma et al., 2011; Robin et al., 2003). We create a sample of mock halo stars similar to our observed sample of metal-poor stars, taking halo stars (i.e. age $\geq$ 11 Gyr) with $\rm Decl.\gtrsim-10$ deg, d $\gtrsim$ 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 which are tightly associated to the CS detection (i.e. the magenta region in Fig. 3IIIa). We find that the smooth halo produces 3,301/155,343 = 2$\%$ stars associated to the CS detection, which is far fewer than the observed fraction of 151/1,679 = 9$\%$. If we concern ourselves with only the new northern detection (restricting ourselves to $b\geq$20 deg), we find that the mock sample has 1,795/107,486 = 1.6% of stars associated to this detection, whereas the observed sample contains 45/1,105 = 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 which lie within our CS detection region. This takes into account the Poisson fluctuation in the background population, under the assumption that the sample of 1,105 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 1,105 smooth stars, is less than 1e-7%, i.e. the significance is more than 5-sigma.

4.2.2 Metallicity

In Fig. 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 Fig.5a), and crucially these metallicities are clearly distinct from the background halo population (see Fig.5c). 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), and so the lack of a low-metallicity tail compared to the BHBs is not surprising.

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\pm 0.04$ dex and a dispersion of $0.23\pm 0.04$ dex, while the K-giants have a mean of $-2.07\pm 0.02$ dex and a dispersion of $0.12\pm 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 globular cluster for which we would expect a dispersion of $\sim$ 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 to this group and find that NGC 5824 (the red star in Fig. 3III) is tightly associated. As shown in Fig. 2b 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 Fig. 5cd. Roederer et al. (2016) found [Fe/H] = -1.94$\pm$0.12 dex from 26 red giants stars (see gold dashed line with transparent band). In a later study, Mucciarelli et al. (2018) found [Fe/H] = -2.11$\pm$0.01 dex (right lime dashed line) from 87 red giant branch stars, and [Fe/H] = -2.20$\pm$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. Fig. 6 shows the properties of our stars in terms of helio-centric 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.

We start by analyzing the stars from the existing detection in the south ($-70\lesssim b\lesssim-20$ deg). We can see from the Galactic rest frame radial velocity distribution ($v_{\rm gsr}$) that our detection appears to be split into two clumps, one with $-75<v_{\rm gsr}<-25$ km/s (cyan circles in Fig. 6) and one with $50<v_{\rm gsr}<100$ km/s (red circles in Fig. 6). The former clump corresponds to the detection in Yam13, which has a negative radial velocity gradient with $b$. The clump at positive $v_{\rm 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. Fig. 6 shows orbit of the cluster integrated forward and backwards 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 ${\rm Dec}<-10$ deg, i.e. these cannot correspond to our debris, as our stars are all located at ${\rm Dec}>-10$ deg.

We first focus on the existing detection from Yam13, given by the cyan circles in Fig. 6. The gray-scale 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 ($\sim 2$ orbital periods), while the latter is around 1.5 Gyr behind ($\sim 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 ($\sim 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 Fig. 3).

If the above picture is correct, it is curious as to why 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 declinations or behind the Galactic plane), but there are other portions which 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\lesssim b\lesssim-20$ deg. This portion of the stream corresponds to pericentre of the orbit and so perhaps the stream is sparsely populated here.

From this orbital analysis, together with the homogeneous metallicities shown in Fig. 5, we conclude that these new detections (both at positive $v_{\rm gsr}$ and in the northern hemisphere) have orbits that are consistent with the Cetus stream 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 Cetus Stream. 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.

5. Conclusions and Discussions

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 Cetus Stream (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 $\gtrsim$ R.A. $\gtrsim$ 50 deg) corresponds to the previous CS from Newberg et al. (2009) and Yam13. The northern counterpart (at -150 $\gtrsim$ R.A. $\gtrsim$ -100 deg) is a new discovery and has the same angular momentum and orbital energy as the southern counterpart.

In order to understand more about the progenitor system, 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 deg (see Fig. 3IIIc). 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 $\sim$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. The fact 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 are 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 $\pm$ 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 $\sim$-2 dex suggests a stellar mass in the range 104M*⊙* – 106M*⊙*, 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 Newberg et al. (2009); Koposov et al. (2012); Yam et al. (2013). The high fraction of BHBs suggest the progenitor system had most of its stars formed 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 6-D 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, is 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 Cetus Stream. 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 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 $\omega$ Cen. The mass-luminosity relation for GCs suggests a stellar mass of over $10^{6}\leavevmode\nobreak\ {\rm M}_{\odot}$. This is comparable to, or even greater than, the stellar mass of the progenitor dwarf estimated above ($10^{4}$ to $10^{6}\leavevmode\nobreak\ {\rm M}_{\odot}$). 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 in a short timescale ($\sim$ 1 Gyr) due to dynamical friction (Binney & Tremaine, 2008). According to the quantitative comparisons between the MW GC population and cosmological zoom-in simulations from the E-MOSAICS project, the MW hosts 6$\pm$1 nuclear clusters from disrupted dwarf galaxies (Kruijssen et al., 2019). NGC 5824, like M54 and wCen, is likely to be the nuclei 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 thanks the support of ”Recruitment Program of Global Youth Experts” of China and NSFC under grants 11390371,11873052,11890694. This project is 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 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, co-founded by the National Astronomical Observatories, Chinese Academy of Sciences and the Alibaba Cloud.