Identifying Stellar Streams in Gaia DR2 with Data Mining Techniques

TL;DR

This paper uses data mining on Gaia DR2 data to identify stellar streams by analyzing integrals of motion, discovering seven streams including four new ones, and linking them to the Gaia-Enceladus progenitor.

Contribution

It introduces a novel combination of data mining, numerical, and statistical techniques to detect stellar streams in large astronomical datasets, including previously unknown streams.

Findings

Seven stellar streams identified, four new discoveries.

Streams are old, metal-poor, and tightly clustered in integral of motion space.

New streams are likely fragments of Gaia-Enceladus progenitor.

Abstract

Streams of stars from captured dwarf galaxies and dissolved globular clusters are identifiable through the similarity of their orbital parameters, a fact that remains true long after the streams have dispersed spatially. We calculate the integrals of motion for 44855 stars, to a distance of 4 kpc from the Sun, which have full and accurate 6D phase space positions in the Gaia DR2 catalogue. We then apply a novel combination of data mining, numerical and statistical techniques to search for stellar streams. This process returns seven high-confidence streams (including four that were not previously known), all of which display tight clustering in the integral of motion space. Colour-magnitude diagrams indicate that these streams are relatively simple, old, metal-poor populations. A combined evaluation of the kinematics and colour-magnitude properties suggests that the previously undiscovered streams are fragments of the Gaia-Enceladus progenitor. The success of this project demonstrates the usefulness of data mining techniques in exploring large datasets.

Tables

Table 1. Table 1: The ranges of the integrals of motion when they have become standardised.

	$L_z$	$L_{⟂}$	$Energy$
Max	10.3	7.05	19.03
Min	-5.48	-1.32	-1.15

Table 2. Table 2: Summary data of all streams we identify in the replica data set. Integrals of motion values are means.

Stream	H99	KGS-1	KGS-2	BGS-1	BGS-2	BGS-3	BGS-4
Population	64	52	27	38	31	13	9
Range (kpc)	1.5	2	4	2	1.5	1.5	1.5
$L_Z$ ($\frac{kmkpc}{s}$)	1200	-1800	-3100	430	-520	540	560
${\sigma }_{L_Z}$	250	420	200	180	240	190	80
$L_{⟂}$ ($\frac{kmkpc}{s}$)	2000	1600	460	400	180	1100	630
${\sigma }_{L_{⟂}}$	300	350	160	180	100	100	60
$E$ ($\frac{k{m}^2}{s^2}$)	-11000	-97000	-65000	-92000	-96000	-110000	-110000
${\sigma }_E$	8000	14000	13000	7000	8500	5700	2900

Table 3. Table 3: Summary data of all streams we identify in the second, more restricted, data set. Integrals of motion values are means.

Stream	H99	KGS-1	KGS-2	BGS-5	Gaia-Enc
Population	41	10	16	53	319
Range (kpc)	1.5	2	4	2	1.5
$L_Z$ ($\frac{kmkpc}{s}$)	1200	-1900	-3000	-700	50
${\sigma }_{L_Z}$	480	210	100	240	210
$L_{⟂}$ ($\frac{kmkpc}{s}$)	2100	1400	380	290	660
${\sigma }_{L_{⟂}}$	390	250	130	140	190
$E$ ($\frac{k{m}^2}{s^2}$)	-11000	-95000	-74000	-115000	-110000
${\sigma }_E$	11000	5900	6800	8000	8200

Table 4. Table 4: Gaia source ID numbers for our detected stream stars. The full machine-readable table is available online; these first few rows are shown here as a guide to form and content. \say Enc is the label for Gaia-Enceladus.

Gaia source ID	Stream
29380806120144384	Enc
30130265029073536	Enc
52624073910834176	Enc
53152251811120128	Enc
133539707282858880	Enc

Equations

$$f(\hat{y})=\frac{1}{2\pi^{n/2}}|\Sigma|^{\frac{1}{2}}\;\rm{exp}\Big{(}\frac{-(\hat{y}-m)^{\prime}|\Sigma|^{-1}(\hat{y}-m)}{2}\Big{)}$$

$$m=[\alpha ,\delta ,\pi ,{\mu }_{\alpha },{\mu }_{\delta },v_{\mathrm{rad}}]$$

$$\Sigma =\left(\begin{array}{ccc}{\sigma }_{RA}^2& \dots & \rho (\alpha ,v_{\mathrm{rad}}){\sigma }_{RA}{\sigma }_{v_{\mathrm{rad}}}\\ ⋮& \ddots & ⋮\\ \rho (\alpha ,v_{\mathrm{rad}}){\sigma }_{RA}{\sigma }_{v_{\mathrm{rad}}}& \dots & {\sigma }_{v_{\mathrm{rad}}}^2\end{array}\right)$$

$${\sigma }_i=\sqrt{N_i^{Model}}$$

$$N_i^{Gaia}-N_i^{Model}>4{\sigma }_i$$

Full text

††pubyear: 2019††pagerange: Identifying Stellar Streams in Gaia DR2 with Data Mining Techniques–A

1 Introduction

Stellar streams are the footprints of our Galaxy’s evolution. Their creation is a by-product of two possible Galactic processes: the end-result of the tidal disruption of globular clusters, or the accretion of dwarf galaxies (Price-Whelan & Bonaca, 2018; Helmi et al., 1999b). Locating and classifying both varieties of streams is, therefore, a means to uncover our Galaxy’s past (Johnson & Soderblom, 1987). Furthermore, the discovery of streams offers several additional benefits to the understanding of Galactic processes. The total number of streams can in principle place a lower limit on past accretion events in the Galactic halo (Bullock & Johnston, 2005), while those created as a result of gravitational disruption can offer insight in the formation and evolution of globular clusters (Balbinot & Gieles, 2018; Bose et al., 2018). Streams can also act as probes into the Milky Way dark matter halo, providing the means to map its mass distribution and shape (Johnston et al., 1996; Ibata et al., 2001; Koposov et al., 2010; Law & Majewski, 2010; Bowden et al., 2015; Bovy et al., 2016; Malhan et al., 2018a; Malhan & Ibata, 2019).

Stellar streams are found using a variety of different search methods. A common theme amongst them is to find clusters of measured stellar parameters in astronomical surveys. Clusters often indicate stars of the same origin. Options range from abundance-, kinematic-, and location-based techniques, or combinations of all three.

Streams formed from disrupted globular clusters maintain their chemical similarity, offering a means to detect them if their populations can be resolved from the stellar background with sufficient contrast. Matched Filter techniques incorporate the colour-magnitude weighting of stars to find structures that belong to the same origin (Grillmair et al., 1995; Rockosi et al., 2002; Balbinot et al., 2011). An abundance-based technique, however, is reliant on spectroscopic data to separate stream members from the background stars. If data quality is poor, it can lead to missed detections (Malhan et al., 2018b). Additionally, detection can fail if there is significant dispersion in the abundance patterns of the stream. Notable streams that have been identified using the matched filter technique include the Palomar 5 tidal tails (Odenkirchen et al., 2001), GD-1 (Grillmair & Dionatos, 2006), Orphan (Belokurov et al., 2006), Eridanus, Palomar 15 (Myeong et al., 2017), and the streams detected in the Dark Energy Survey (Shipp et al., 2018).

Finding clusters of stars within kinematic spaces can also lead to the detection of streams. Conventional techniques rely on either space-velocities to locate co-moving groups of stars or identifying clusters in a data set’s \sayintegrals of motion values. Integrals of motion refer to the angular momentum actions defined by Binney & Tremaine (2008) and total orbital energy of a star. These quantities are ideal for stream searches as they are argued to be conserved, or quasi conserved (Helmi et al., 1999b; Helmi et al., 1999a; Myeong et al., 2018b; Koppelman et al., 2018; Meingast et al., 2019). Both techniques depend on locating clustered sets of kinematic parameters and proving these clusters to be statistically significant.

Successful stream searches using velocity parameters have resulted in the detection of the Aquarius, the VelHel-1 to -9 streams, and the S1-4, C1-2, Rg1-7 and Cand8-17 sub-structures (Williams et al., 2011; Helmi et al., 2017; Myeong et al., 2018a). While difficulties can arise if there is a strong velocity gradient along the stream, the technique is viable for detecting streams which have not undergone a significant degree of phase mixing. Phase mixing is a collective name for kinetic energy exchanges between the stream stars and other constituents of the Galaxy (Klement, 2010).

If phase mixing is an issue, integrals of motion can be relied upon instead of velocities to search for streams kinematically. Relying on conserved quantities has led to the discoveries of the H99, $S_{1},S_{2}$ and $S_{3}$ streams (Helmi et al., 1999b; Klement et al., 2009) and the rediscovery of S1-4, C1-2, Rg1-7 and Cand8-17 in a follow-up study provided by Myeong et al. (2018b). Integrals of motion provide robust information for stream searches because they are conserved over long periods. However, they require knowledge of the six-dimensional phase space positions of each star (Malhan et al., 2018b). This information is usually difficult to obtain for a large population of stars and is very challenging to acquire for very distant stars. This impediment limits the scope of integrals of motion-based studies to the local Galaxy. Nevertheless, many studies use the integrals of motion to search for streams. For example, in Gaia DR2 Koppelman et al. (2018) searched halo stars within 1 kpc of the Sun, and found 5 potential stream candidates with clearly identifiable clustering in their integrals of motion.

Attempts to mitigate the weaknesses in each of these approaches usually involve combining methods. The Cetus Polar Stream (Newberg et al., 2009) was detected using both metallicity and velocity information to differentiate the stream from the tidal tail of the Sagittarius Dwarf galaxy (Law et al., 2005). The recently developed STREAMFINDER algorithm (Malhan & Ibata, 2018) combines aspects of spatial, chemical and kinematic techniques in conjunction with probabilistic arguments to search for stellar streams and characterise their shapes and their distribution in the Galaxy. STREAMFINDER has had a great deal of success in locating previously unknown streams within the Gaia database, naming these streams after great lakes in Greek and Norse mythology (Ibata et al., 2018; Ibata et al., 2019). Combined search attempts are ideal if the required information is available since an agreement between multiple independent techniques adds to the legitimacy of any detection made.

With the continual increase in the size of wide-field astrometric, imaging and spectroscopic surveys (Luo et al., 2015; Blanton et al., 2017; Kunder et al., 2017; Casey et al., 2017; Buder et al., 2018; Gaia Collaboration et al., 2018a), the necessity of using algorithmic search methods is becoming more and more apparent (Lloyd, 1982; Davis et al., 1985). Data mining algorithms that search for clusters in data sets are well suited for this task as they can handle more challenging and higher-dimensional problems objectively, allowing one to focus primarily on verifying and characterising the properties of the discovered streams. In this study, we use the clustering algorithm DBSCAN (Ester et al., 1996) to search for clustering in the integral of motion space using stars with full 6D phase space information in the Gaia DR2 data set (Gaia Collaboration et al., 2018a), and analyse the colour-magnitude diagrams of these candidate streams. Our goal is to replicate the results of Koppelman et al. (2018, 2019) while locating additional stream candidates, using a more stringent data quality selection on the Gaia data set.

We describe the data set in Section 2, outlining the subset restrictions and quality criteria we applied. Section 3 explains the process of our stream finding method: 3.1 outlines how we actually applied the data-miner to the data set, 3.2 will explain how we chose our search parameters. The determination of statistical significance is explained in 3.3. Section 4 presents our results, with our findings broken up into two parts; Section 4.1 presents a comparison of our results to the literature, and our finalised results presented in 4.2. Discussion of the results and interpretation of our overall findings of the paper are presented in Section 5.

2 The Data Set

Our search for potential stream candidates use the astrometric and kinematic data available in the second Gaia data release (Gaia Collaboration et al., 2018a). In this data set, 7,224,631 stars have full six-dimensional phase space and photometry information. From this subset, we create two data sets. The first, the \sayreplica data set, uses the same selection as Koppelman et al. (2018), which we will use to compare the results of our method to theirs. The second data set, the \sayrestricted data set, is a subset of the replica data set with stricter selection criteria to help ensure that our results are reliable.

Creating the replica data set involved only accepting stars with a parallax uncertainty of less than 20%, leaving us with a total of 6,447,952 objects. Restricting the parallax uncertainties to this level enables inversion of the parallax measurement to provide a distance measure. Inverting parallaxes typically creates difficulties when attempting to incorporate uncertainties due to the hyperbolic nature of the conversion. A potential remedy for this issue is to use Bayesian distance priors to calculate the uncertainties; however, inverting parallaxes with percentage errors less than 20% is an acceptable alternative (Bailer-Jones, 2015). We use this replica data set to compare our results to Koppelman et al. (2018).

For the restricted data set, we apply a correction to the Gaia data set, and then apply an additional data quality selection. We correct for the global Gaia parallax zero-point offset by applying a $+0.052$ mas correction to the parallax measurements (Zinn et al., 2019; Leung & Bovy, 2019) which did not change the overall number of stars in the data set. The quality selection cut only accepted stars with a minimum number of radial velocity transits $rv\_nb\_transits$ $\geq 5$ (Marchetti et al., 2018). This cut reduced the data set to 4,533,363 objects red, and was the only additional cut that we applied.

We convert the equatorial sky coordinates, proper motions, parallax and radial velocity measurements from Gaia to Cartesian Galactocentric coordinates using astropy (Astropy Collaboration et al., 2013). In this system, the Sun lies on the X-axis of the Cartesian plane at a distance of 8.3 kpc from the origin (Gillessen et al., 2009) and 27 pc above the Galactic midplane (Chen et al., 2001). The Z-axis projects towards the Galactic north pole, and the Y-axis aligns with the direction of Galactic rotation. For velocity parameters we use the convention adopted by Johnson & Soderblom (1987) with the Cartesian representation of the motions of stars defined as (U,V,W) where U is the direction moving radially toward the Galactic centre, V is the circular tangential velocity in the direction of Galactic rotation, and W is the velocity directed towards the Galactic North Pole. We correct for solar motion using the peculiar velocity correction values of $(U,V,W)_{\odot}=(11.1^{+0.69}_{-0.75},12.24^{+0.47}_{-0.47},7.25^{+0.37}_{-0.36})\,\rm{kms^{-1}}$ supplied by Schönrich et al. (2010) and Bovy (2015) to adjust the circular velocity and Solar radius while also placing all stars in the Local Standard of Rest (McMillan, 2017) $\text{V}_{\text{LSR}}=232\;\rm{kms^{-1}}$.

The Z component of angular momentum $L_{Z}$, perpendicular angular momentum $L_{\perp}$111$L_{\perp}$ is calculated by adding the X and Y components of angular momentum together in quadrature: $L_{\perp}^{2}=L_{X}^{2}+L_{Y}^{2}$, and total energy $E$ values are calculated using the library package gala (Price-Whelan, 2017). To model the gravitational potential of the Galaxy, we use the potential class MilkyWayPotential in gala. The constituents of this potential consist of the following analytic models: a Hernquist potential to model the bulge and nucleus of the Galaxy (Hernquist, 1990), a Miyamoto-Nagai (Miyamoto & Nagai, 1975) potential to shape the Galactic disk (taken from Bovy, 2015), and a standard Navarro-Frenk-White potential (Navarro et al., 1996) to model the halo. Using the phase-space coordinates and the specified gravitational potential, we can calculate the integrals of motion for each star in the data set.

Following previous stream search publications, we focus on the halo of the Galaxy (Helmi et al., 1999a; Klement et al., 2009; Williams et al., 2011; Bernard et al., 2014; Koppelman et al., 2018; Ibata et al., 2019). We choose halo stars in Gaia DR2 by selecting stars with large total velocities relative to the Local Standard of Rest (Nissen & Schuster, 2010; Koppelman et al., 2018). Using a selection of $|\text{V}-\text{V}_{\text{LSR}}|>210\;\rm{kms^{-1}}$ results in a reduced subset of 152,865 stars for the replica data set, and 62,527 stars for the restricted data set. It is important to note that this approach excludes halo stars with disk-like velocities from the search sample (Nissen & Schuster, 2010; Bonaca et al., 2017; Koppelman et al., 2018, 2019) while also resulting in some contamination from thick disk stars which may have large vertical velocities.

3 Stream Searching With DBSCAN

Our stream search begins within the local Solar neighbourhood (out to 1 kpc) before extending the sample to include stars across a broader region of the Galaxy. We observe how far any streams detected within the local halo extend before we can no longer confidently identify them. We only search for potential stream targets at a distance $<4$ kpc, as increasing the range further leads to cumbersome computational times.

A fundamental property of clustered data points is the density contrast between points inside the cluster versus points outside. Density-Based Spatial Clustering of Applications with Noise (DBSCAN, Ester et al., 1996) is a data mining method which utilises this fact to locate clustered data points. The algorithm requires two input parameters: a characteristic distance $\epsilon$ (Ester et al., 1996; Chen et al., 2018) and a minimum cluster member number termed $MinPts$. DBSCAN selects an arbitrary point within the data set and retrieves all other data points reachable within a distance of $\epsilon$. If the number of other data points within $\epsilon$ exceeds $MinPts$, the selected data point is a core-point. DBSCAN then moves from the core-point to another point within the $\epsilon$ range and determines if that point too also satisfies the core-point criteria. If the test is successful, the cluster size expands to include that point. The process then repeats with another point within range. Data values are assigned a border-point label if they lie within a distance of $\epsilon$ of a core-point, but fail to classify as a core-point. A cluster is complete when border-points surround the core-points. At this stage, DBSCAN moves to a data point outside the cluster and repeats the process. Any data point that does not have the specified $MinPts$ number of data points within a distance of $\epsilon$, and no core-points in range is labelled as a noise-point. Every cluster is given a unique integer label, with noise-points taking a label of -1.

Figure 1 offers a visual representation of how DBSCAN operates. In the example, $MinPts=3$ and $\epsilon=1$. One can see that two points have three data values within a distance of $\epsilon$, and so have been classified as core-points. Two points have fewer than three data values within range and so are classified as border-points while two data points have no neighbours within a distance of $\epsilon$ and are given noise classifications.

A significant advantage to using DBSCAN is its versatility and applicability to data of arbitrary dimensions. Isolating clusters in higher dimensional space creates a stronger argument that an actual cluster is present, as it is less likely that the apparent structure is coincidental. Additionally, DBSCAN will produce the same cluster patterns every time provided the search parameters remain the same. These aspects are highly useful for stream detection.

3.1 Applying DBSCAN to the data set

DBSCAN works under the assumption that all data points are located precisely at their recorded position and draws clusters where the points stand. However, this is not realistic as measured values have uncertainties, which means that DBSCAN may potentially classify clusters that are statistical anomalies created by noise in the data rather than actual clustered sets of stars. Therefore, we perform searches on multiple Monte-Carlo iterations of the search subset to take the measured uncertainties into account.

Positional components, proper motions, parallaxes, and radial velocities were modelled as multivariate Gaussians in a framework similar to Marchetti et al. (2018) with the joint probability distribution of these variables represented by the following equation:

[TABLE]

which describes the probabilities with which the random variable $\hat{y}$ takes on values jointly in its $n$ elements. The variables $m$ and $\Sigma$ refer to the mean vector and covariance matrix of the distribution respectively. The mean vector is given by

[TABLE]

where $\alpha$ and $\delta$ are right ascension and declination, $\pi$ is the parallax, $\mu_{\alpha}$ and $\mu_{\delta}$ are the proper motions in right ascension and declination, and $v_{\rm rad}$ refers to the radial velocity. In addition to this, we constructed a $6\times 6$ covariance matrix of

[TABLE]

where $\sigma_{i}$ denotes the standard deviation of each variable in the mean vector, and $\rho(i,j)$ denotes the correlation coefficients of the parameters $i$ and $j$.

With these parameters we can repeatedly sample from $\hat{y}$, providing random samples of each star’s phase space parameters in line with what we would expect based on the known correlations and uncertainties. We note that radial velocities within the Gaia DR2 data set are uncorrelated to the other astrometric parameters. Therefore, correlation coefficients between radial velocity and the astrometric parameters are zero, and the standard deviations are assumed to follow a Gaussian distribution.

Using the covariance matrix, we sample 200 potential integrals of motion values for each star in the search subset. A sample of 200 was large enough to reproduce the results consistently while also minimising the search time. DBSCAN is then run on each of these subsets to produce an initial set of clusters within each random sample.

Clusters should, in principle, remain coherent even when subjected to the random deviations produced by uncertainties. If any cluster drawn by the algorithm in one Monte-Carlo iteration does not appear in another, the cluster is likely just an artefact produced by the random deviations of the Monte-Carlo sampling. We remove any clusters that do not share at least four stars in common with another cluster in one of the other Monte-Carlo sets, which effectively cleans the data set of statistical anomalies that fail to re-occur when the data is allowed to vary within the errors. This process effectively eliminates clusters that cannot be conclusively proven to exist based on the uncertainties of Gaia measurements.

At this point, the clusters found in each Monte-Carlo subset are combined to create an aggregate subset. In this subset, we apply DBSCAN one final time with re-optimised search parameters. This final application of DBSCAN acts to join any partially formed groups together with other partial pieces. After the clusters are merged, they are extracted from the data set and have all duplicates removed.

3.2 Specifying Input Parameters

When choosing our input parameters we set $MinPts=4$. This specification enables sparsely populated stream candidates to be accepted by the algorithm, while still imposing that some clustering must be present initially. Our method of choosing an $\epsilon$ value relies on the use of a 4-distance plot (Ester et al., 1996), which is a ranked plot of the standardised distance from each data point to its fourth nearest neighbour. As the magnitudes of the integrals of motion can vary significantly, these quantities need to be standardised. We do this through normalising each of the integrals of motion individually, taking the mean to be zero and scaling to unit variance. Normalisation is achieved using the StandardScaler class in the scikit-learn library (Pedregosa et al., 2011). Extracted $\epsilon$ values are measured in this standardised distance, and are unitless.

Deriving the logic behind a 4-distance plot comes from the inherent characteristics of a clustered data set. In a clustered data set, the spacing between data values falls into two categories: the distances between members inside a cluster, and distances between data values which do not belong to the same cluster. Distances between two points inside a cluster are small relative to points outside the cluster. Additionally, points which are not part of the same cluster can also experience a significant variance in their distance from one another. Therefore, two populations of distances exist based on the different variances of the distances for a given data set. Determining the $\epsilon$ value is achieved through creating a ranked scatter plot of the distances between the points in the data set, with the highest rank given to the largest distance. The two populations can then be separated visually by finding the point in the plot where the scatter plot effectively becomes smooth. This junction represents the point of separation between the two populations, and we take this to be the value of $\epsilon$. Figure 2 is an example of such a plot created during the process of determining the $\epsilon$ value for our stream search. The point at which the line changes from discrete to continuous is highlighted and taken as the $\epsilon$ value.

In order to take uncertainties into account, we create 4-distance plots for a random selection of 50 of the Monte-Carlo samples used in the search. We extract the $\epsilon$ from each of these 4-distance plots and take the average value. The average $\epsilon$ is a better measure of the clustering in the data set. Our choice to use 50 4-distance plots produced an accurate measure of $\epsilon$ without excessive time spent on calculations. The standardised $\epsilon$ values which we obtained for the replica data set were $\epsilon_{replica}=0.286$, and for the restricted data $\epsilon_{restricted}=0.306$. To provide context the ranges of the integrals of motion in this standardised form are provided in Table 1.

3.3 Statistical Significance

We determine statistical significance for our DBSCAN clusters in integral of motion space through the use of the Milky Way observation simulator galaxia (Sharma et al., 2011). galaxia is capable of producing realistic synthetic representations of Galactic observations, enabling us to apply the same selections that we applied to Gaia data to create a manufactured halo subset. With the simulated data, we can compare the clusters found by our method to the expected stellar density provided by galaxia’s output. If the number of stars in the space the cluster occupies in Gaia is significantly higher than the simulated equivalent, we have confidence that the stream is a real kinematic structure.

We created three synthetic Galactic surveys and combined them in order to produce a data set with the same number density as Gaia. Each synthetic survey was created by running galaxia with a different random seed integer, to create three plausible Galactic surveys. All default simulation parameters were used except for an apparent magnitude range of $M=4-13$ to mimic Gaia data (Gaia Collaboration et al., 2018a). With the simulations complete, local halo subsets were created using the same method as described in Section 2. Stars were then randomly drawn from this overabundant set to produce the same number of stars as the Gaia halo.

The statistical significance of the detected clusters was determined through the use of a multidimensional histogram, used to break the data into 3-dimensional cells of $L_{Z}$, $L_{\perp}$, and $E$. Bin sizes for the histogram were $\Delta L_{Z}=400~{}km~{}kpc~{}s^{-1}$, $\Delta L_{\perp}=200~{}km~{}kpc~{}s^{-1}$, and $\Delta E=8500~{}km^{2}~{}s^{-2}$. Cells that did not have stream stars within them were then removed from this grid, leaving a 3-dimensional object that captured the volume in integral of motion space where the streams resided. Our choice of cell size split the integrals of motion into a $20\times 20\times 20$ grid of equally spaced cells. This cell size was small enough to be effective in capturing the general shape of the stream in integral of motion space, while not being so small as to prevent calculating useful Poisson statistics for the stars in each cell. A corresponding object comprised of the same cells was then created in the synthetic set, and the number of stars in those cells were counted. The numbers of stars contained within each object were defined as $N_{i}^{Gaia}$ and $N_{i}^{Model}$ respectively, while we defined the standard deviation of the synthetic cells to be

[TABLE]

Our criterion for statistical significance is the same as is used in Williams et al. (2011)

[TABLE]

Upon passing this test, we declare the cluster in question to be a statistically significant grouping.

4 Results

As described in section 2, we are considering two different subsets of Gaia data: the \sayreplica data set, which uses the same selection as in Koppelman et al. (2018), and the \sayrestricted data set, which is the replica data set with a correction applied to the parallax zero-point and only considering stars with five or greater independent radial velocity measurements. We will discuss the results of our stream searches in these two data sets separately.

4.1 The replica data set

In the replica data set, we initially searched for streams out to 1 kpc of the Sun before incrementally increasing the search range out to 4 kpc. Our goal was to locate any streams residing in the local solar neighbourhood, and then test how far they extended beyond this distance. The purpose of this was to compare our results with Koppelman et al. (2018) who found five: the H99 stream, initially discovered by Helmi et al. (1999a), and four additional stream candidates. We refer to these four a streams as Koppelman-Gaia-Stream-1 (KGS-1), KGS-2, KGS-3 and KGS-4.

Through our efforts, we were able to resolve the H99, KGS-1, and KGS-2 streams using DBSCAN. Our method identified the majority of the previously reported stars said to belong to these streams and also indentifies new potential members. However, we were not able to resolve KGS-3 and KGS-4. Furthermore, it appears that some of the halo sub-structures found in Myeong et al. (2018a, b), overlap with the H99 and KGS-1 streams. The sub-structures found in these publications used a cross match of Gaia DR1 and the Sloan Digital Sky Survey data release 9, termed the SDSS-Gaia catalogue. In addition to confirming these known streams, we also identify four previously unknown sub-structures, which we will call Borsato-Gaia-Stream-1 (BGS-1), BGS-2, BGS-3, and BGS-4. Table 2 lists the number of stars, distance range, means and standard deviations of the integrals of motion for each of the seven streams, and Figure 3 shows the integral of motion space for all stars in the replica data set, with the streams plotted in separate colours. Halo stars not marked as members of a stream by DBSCAN are plotted in the background as smaller opaque points, to better orient the reader.

4.1.1 The Helmi stream

The H99 stream remained statistically significant out to a distance of 1.5 kpc, with a total population of 64 stars. At distances greater than 1.5 kpc, our method could no longer resolve H99. The increased number of stars in the larger data set likely led to some overlap between the uncertainties of the stream stars and the surrounding population. This overlap would have prevented DBSCAN from creating a complete border. As a result, no cluster formed. We find that 55 of the 64 stars are common to the H99 stream found in Koppelman et al. (2019), though that study also found additional members as far away as 5 kpc. Of additional note is the fact that our method does not appear to cluster the H99 stars at higher energy values. The fixed $\epsilon$ value is likely preventing the inclusion of these stars as part of the cluster.

4.1.2 The Koppelman-Gaia-streams

The two streams from Koppelman et al. (2018) that we identify in the replica data set, KGS-1 and KGS-2, are depicted in blue and red in Figure 3. In Koppelman et al. (2018), the search for streams was conducted within a 1 kpc range of the Sun, and so we restricted our initial search of the replica data set to the same range, and then repeated the search in larger volumes. For the KGS-1 stream, we identify 23 stars within 1 kpc of the Sun, of which 12 are also found in Koppelman et al. (2018). When extending the search range, we found KGS-1 extended out to 2 kpc from the Sun, and containing a total of 52 members. A noticeable peculiarity of this stream are the three stars that are offset from the rest of the group, mainly in $L_{\perp}$ and $E$. Follow-up analysis of these stars revealed that they have substantial error bars. The size of these errors was likely enough for DBSCAN to include these stars in some of the iterations, meaning that these stars are false positive detections.

For KGS-2 we identify 12 stream stars within 1 kpc of the Sun, of which 10 are also reported in Koppelman et al. (2018). Moreover, the KGS-2 stream remains coherent out to a range of 4 kpc with a final member count of 27. We suspect that the stream would remain statistically significant beyond this point but the increasing number of stars obtained by expanding the search range led to impractical calculation times for our stream searches.

Figure 4 shows the overlap of our stars (plotted as coloured circles) with the H99, KGS-1 and KGS-2 streams (shown as coloured squares). Stars included in both data sets have both a circle and a square. H99 contains all stream stars out to a range of 1.5 kpc, which was the maximum range our method could resolve the stars. KGS-1 and KGS-2 contain all stream stars out to 1 kpc, the maximum range that Koppelman et al. (2018) reported. Restricting the ranges in this fashion allows one to fairly compare the results of each stream search, as the distances are consistent.

4.1.3 The Borsato-Gaia-streams

Similar to the H99 and KGS-1 and -2 streams, we first detected the BG-streams initially in the local solar neighbourhood, before examining how far they extended. The final results find that the BGS-1 stream extends out to 2 kpc, with a total of 38 members. BGS-2, BGS-3 and BGS-4 all stayed statistically significant out to 1.5 kpc and contained a final member count of 31, 13, and 9 members respectively.

4.2 The restricted data set

We applied the additional quality cuts, correcting for the global parallax offset and only accepting stars with $\geq 5$ radial velocity transits, and repeated the stream search method. Figure 5 presents the results of our search attempts with this new data set. Similarly to the previous search, our method was able to locate the H99, KGS-1, and KGS-2 streams again. The original Borsato-Gaia-Streams: BGS-1, -2, -3 and -4, were not detected in this new data set. The majority of the stars in these streams ($\sim 70\%$) were removed from the data set when the additional quality cut was made. As these streams were located in a much more dense region of integral of motion space, it is likely that the remaining stars did not form a dense enough cluster for the data-miner to relocate the streams in the restricted data set.

However, we do identify two additional substructures located in the centre of the integral of motion plot that were not apparent in the replica data set. The larger of the two appears to be the Gaia-Enceladus dwarf galaxy remnant (Helmi et al., 2018) as the cluster is located in the same region, and has its characteristic shape, a central bump in the $L_{Z}$ vs $E$ plot. The other cluster we have named BGS-5 because it does not correspond to any of the streams identified by Koppelman et al. (2018), or the Borsato-Gaia-Streams found in the replica data set. BGS-5 does appear to almost overlap with BGS-2, however, the two artefacts are distinctly separate when considering all three dimensions of the integrals of motion, and share no stars in common.

The substructure we suspect is Gaia-Enceladus contains 319 members in total and remains statistically significant out to 1.5 kpc, while BGS-5 contains 53 members and remains significant out to 2 kpc. Table 3 presents these new streams with their average integrals of motion and standard deviations, along with the recovered H99, KGS-1 and -2 streams and their revised member totals. The three rogue stars originally included in KGS-1 were not included in the revised stream, implying the quality cut acted to remove these stars.

Figures 6 and 7 show calculated orbits in Galactocentric $R$ vs $Z$ and $X$ vs $Y$ coordinates and absolute colour-magnitude diagrams for the stream members for the restricted data set. Equivalent figures for the replica data set can be found in Appendix A. Using each star’s phase space coordinates as its initial conditions, the MilkyWayPotential class in gala, and a time step of 1 Myr, we integrate forward over each star’s orbit for 1000 steps, which totals 1 Gyr. The result of this integration produces an Orbit object, which contains all the relevant parameters of the star’s orbit at each time step (Bovy et al., 2016). The orbit plots for each stream show each stream member’s orbit for 1 Gyr, with the current positions of the stars shown in yellow.

Our colour-magnitude diagrams were made using the $G$, $G_{\rm BP}$ and $G_{\rm RP}$ magnitudes available from Gaia. All stream stars showed appropriate Gaia data quality, satisfying phot_bp_rp_excess_factor $<1.3+0.06(G_{\rm BP}-G_{\rm RP})^{2}$ (Gaia Collaboration et al., 2018b), with the exception of one star in the H99 stream, which had an phot_bp_rp_excess_factor $=2.3334208$ and has been removed. The phot_bp_rp_excess_factor delimits photometry that can be considered to be well-behaved in the Gaia data set. We corrected for reddening and extinction using the publicly available extinction map from the Infrared Science Archive (IRSA) which uses reddening estimates provided by Schlegel et al. (1998), and converted the apparent magnitudes to absolute magnitudes using the Gaia parallaxes (Gaia Collaboration et al., 2018b). Dartmouth isochrones (Dotter et al., 2008) were fitted by creating plausible tracks that ran through the body of the stellar population. Upper and lower bounds for the isochrones were created by varying the ages and metallicities of the isochrone tracks such that the average distance between the fitted isochrone and the isochrone bounds was no larger than $\Delta|M_{v}|=0.35$ and $\Delta|G_{\rm BP}-G_{\rm RP}|=0.1$. Any points that lay outside of the isochrone bounds were considered points of contamination. To quantify this contamination we took the ratio of stars that lay outside the isochrone bounds and divided by the total population of stars.

We begin by analysing the H99 stream orbits and colour-magnitude diagrams. H99 appears to have lost all of its spatial coherence and is dispersed throughout the Galaxy. This stream is considered to belong to a dwarf galaxy progenitor with an approximate mass of $\sim M^{8}_{\odot}$ and is projected to contribute as much as $\sim 10-14\%$ of the total stars in the Galactic halo in multiple stream-like structures (Koppelman et al., 2019). Such multiple structures are likely the cause of the spatial dispersion in the orbits, as each stream-fragment may be on slightly different paths leading to a tangle of orbital paths that span a large portion of the sky. The colour-magnitude fit for this stream suggests an age of 12 Gyr and a metallicity of $[M/H]=-1.5$, which is in line with the findings of Koppelman et al. (2019). There are, however, stars which do not fall within the isochrone bounds. Using our dispersion approximation, we calculate that 32.5% fall outside the isochrone bounds. As the population of this stream does not form a tight isochrone track, this progenitor is likely a dwarf galaxy containing a range of stellar populations.

KGS-1 has a significant vertical component to its orbit, with the extreme members of the stream reaching almost 20 kpc both above and below the Galactic plane while tracking a well defined elliptical path in the $X$ vs $Y$ plane. Analysing the colour-magnitude plots suggests that this stream it is old and metal-poor, with an approximate age of 11 Gyr and a metallicity of $[M/H]=-1.5$. The majority of the stars lie within the isochrone bounds with only 30% of the stars lying outside.

The KGS-2 stream is also distinct from the other stream populations. The stars on the orbits in this stream fly out into the outer halo regions of our Galaxy before swinging back towards the Galactic centre. This stream satisfies a 10 Gyr fit with a metallicity of $[M/H]=-1.5$. 27% of the stars in this colour-magnitude plot do not belong to this population of stars.

The next halo sub-structure that we will analyse is what we believe to be a remnant of the Gaia-Enceladus dwarf galaxy remnant. This sub-structure contained by far the most stars out of all the stream candidates found. The orbital plots in Figure 7 show that these stars will span the entire Galaxy, with a large fraction of them on highly radial orbits. The colour-magnitude diagram illustrates a clear population of stars with an approximate stellar age and metallicity of 9 Gyr and $[M/H]=-1.3$. The isochrone characteristics imply a slightly younger age than Helmi et al. (2018); however, our metallicity measure is in agreement. Our dispersion measurements suggest that 24% of the stars found in this stream do not match the isochrone fits. The colour-magnitude diagram clearly shows a coherent population of old, metal-poor stars with a characteristic spread that suggests that it is a remnant dwarf galaxy.

We consider the other new substructure we resolved in the restricted data set, the BGS-5 stream. All the stars in this stream appear to have highly radial orbits, manifested in the low $L_{Z}$ values of the stream. BGS-5 is also predicted to spend a significant amount of time orbiting through the thick disk of the Galaxy, seen in the $Z$ vs $R$ plot. Using our isochrone fitting procedure, we estimated the age of this stream of stars to be 11 Gyr with a metallicity value of $[M/H]=-1.4$. These values differ significantly with the Gaia-Enceladus fit, suggesting that they are two distinct populations of stars, even though they border each other quite closely in integral of motion space (see Figure 5). The dispersion is also quite low, with only 12% of the stars falling outside the isochrone fits.

Finally, we compare the streams that we have found in this paper with the halo sub-structures found in Myeong et al. (2018a) and Myeong et al. (2018b), with the latter using the conserved actions of a star’s orbit to find stream candidates. Figures 8 and 9 indicate the locations of these structures with their integrals of motion, with the former plot illustrating likely sub-structures which share the same integrals of motion spaces as our stream. The integrals of motion for these streams were calculated using the phase space coordinates provided in Myeong et al. (2018a) and Myeong et al. (2018b). The points in these plots represent the average $L_{Z}$, $L_{\perp}$ and $E$ values of these streams, while the ellipses show the standard deviation in those values, calculated from the errors provided in Myeong et al. (2018a, b). All the streams we detected in this study have been plotted in grayscale, and all non-stream stars have been removed to better illustrate the positions of the Myeong streams. In Figure 8, a number of the streams found in Myeong et al. (2018a) and Myeong et al. (2018b) overlap with H99 and KGS-2. Of significant note is the Rg6 stream (where Rg stands for retrograde), which appears to almost entirely overlap with KGS-2 stream in the restricted data set, implying that there is a genuine chance that this sub-structure appears in both the Gaia data set and the SDSS-Gaia catalogue. It could very well be the case that Rg6 and KGS-2 are the same strucutre. The additional streams that lie close to H99 and KGS-2 without overlapping are likely to be fragments of a more extensive progenitor system. Dwarf galaxies are predicted to fragment during the accretion process (Koppelman et al., 2019), rather than forming one monolithic cluster in integral of motion space. While only three of the streams identified in the replica data set were found after the transition to the restricted data set, the persistence of these features, and the consistency with some of the substructures identified by Myeong et al. (2018a, b), give us confidence in the effectiveness and reliability of our method.

5 Discussion

We have successfully demonstrated that searching for clusters of stars in the integral of motion space using data-mining techniques is a viable way to search for streams. Our method has a well-defined approach in locating and classifying streams, making the results easily reproducible. We now discuss the significance of these detections, as well as the issues and improvements we could apply in order to enhance our results.

Of primary significance is the fact that our results using the restricted data set not only verified the existence of previously discovered streams but also have provided evidence that they extend past their ranges specified initially. The extension of the search distance outside the Solar neighbourhood has led us to understand that the H99 and KG-streams likely span a large portion of the Galaxy, which is in alignment with expectations (Koppelman et al., 2019). It also may be the case that the additional stream members that we found in the \sayreplica data set but not the \sayrestricted set are also stream members. We plan to verify if this is true when the Gaia DR3222Improved photometry and astrometry are planned to be released in late 2020. See https://www.cosmos.esa.int/web/gaia/release for details. catalogue is released, and the number of RV transits for these stars has increased to satisfy the quality cut criteria.

The quality cuts are responsible for our inability to resolve the original Borsato-Gaia-Streams in the restricted data set. We analysed the populations of the Gaia-Enceladus and BGS-5 streams and found that none of the stars in these clusters were members of the original BG-streams. The removal of a substantial number of the member stars after the quality cut was applied left a cluster with a density too low for the data miner to resolve. Conversely, Gaia-Enceladus and BGS-5 were not resolved in our initial search, likely because measurement errors on a number of lower quality stars in the data set washed out any coherent shape that DBSCAN would have found. With only the quality data remaining, this was no longer the case, and these coherent bodies of stars were found.

The orbital plots of these streams depict populations with a high degree of dispersion throughout the Galaxy. While the orbit paths for these stars do overlap in some instances, it is highly unlikely their existence would be detected if it were not for locating the streams using the stars’ integrals of motion. The colour-magnitude plots show reasonably well defined populations, indicating that the stream stars are physically related and are not simply a selection of field stars.

Furthermore, our results provide evidence to suggest that dwarf galaxies fragment when accreted by their host galaxy, leading sub-structures to form within the larger cluster of the integrals of motion. We see this with the Myeong et al. (2018a) streams, as their smaller sizes overlap with the areas where our streams are known to exist. This finding may imply that these structures originated from the same progenitor, and are in effect pieces of a larger whole. Fragmentation of a stream has already been predicted to have occurred with the H99 stream (Koppelman et al., 2019). These findings suggest that while it is possible to locate stars of a similar origin using conserved quantities, it may be that other clusters of stars from the same progenitor too exist within the Galaxy, with integrals of motion which differ slightly. Further analysis of the properties of these streams are required to fully determine how and when they entered the Milky Way.

While our stream search technique has produced notable positive results, it could be improved. Our method of stream detection consistently produces results that differ slightly compared to the literature. DBSCAN failed to find two streams found in Koppelman et al. (2018) KGS-3, and KGS-4, as well as failing to include all the stars in the streams that it did find. It is likely that a systematic difference in our search procedure is driving these differences in the results. The difference may not necessarily be a negative, as our search method can find clusters in the centre of the integrals of motion distribution where the stellar density is much higher. However, there are caveats to using DBSCAN. Our method can produce false positives for cluster membership. However, this can be mitigated if higher quality data is used. Furthermore, choosing a $MinPts$ value of four and relying on a 4-distance plot to select an $\epsilon$ value works well only if clustering size remains consistent throughout the data set, which generally is not the case for streams. Ideally, a data mining method that can vary its clustering range and produce clusters of multiple sizes would be optimal. Such a method would likely be more successful in capturing the stars with higher energy values in the H99 stream in Figure 3. In addition to this, it is apparent that searching for clusters in integrals of motion is not immune to contamination of stars that are not part of the stream. More analysis is required when finding potential stream candidates in this way in order to determine whether the stars found do form a population of stars with the same origin. Caveats aside, while the algorithm may lack the required complexity to resolve all streams, and has a tendency to introduce contaminants, it is robust enough to produce consistent clustering in the data.

DBSCAN can, however, still consistently detect clusters even with some variation in the $\epsilon$ parameter. Our optimisation procedure changed the $\epsilon$ value when switching from the replica to the restricted data set, yet we were still capable of resolving the H99, KGS-1 and KGS-2 streams. The fact that results remained consistent for these streams means that DBSCAN is a robust search method, even when some leeway is permitted to vary the search parameters. The results of this study have shown that it is possible to locate stream candidates using the DBSCAN data miner to search for clusters in the measured integrals of motion of a star. While there are issues associated with the data miner, it is useful in locating clustering in integrals of motion which warrant legitimate detection of streams.

Acknowledgements

This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://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. Parts of this research were conducted by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. JDS and SLM acknowledge the support of the Australian Research Council through Discovery Project grant DP180101791. NWB and SLM acknowledge the support of the UNSW Scientia Fellowship program.

For data mining and analysis, the following software packages have been used: astropy (Astropy Collaboration et al., 2013), gala (Price-Whelan, 2017), scikit-learn (Pedregosa et al., 2011), numpy (van der Walt et al., 2011), matplotlib (Hunter, 2007), seaborn (Waskom et al., 2016), pandas (McKinney, 2010), and topcat (Taylor, 2005). This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France (DOI 10.26093/cds/vizier). The original description of the VizieR service was published in A&AS 143, 23.

This research has also made use of the NASA/ IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

Appendix A Orbits and Colour-Magnitude Diagrams for the Replica Data Set