# The Streams of the Gaping Abyss: A Population of Entangled Stellar Streams Surrounding the Inner Galaxy

Rodrigo Ibata, Khyati Malhan, Nicolas Martin

► To cite this version:

Rodrigo Ibata, Khyati Malhan, Nicolas Martin. The Streams of the Gaping Abyss: A Population of Entangled Stellar Streams Surrounding the Inner Galaxy. The Astrophysical Journal, 2019, 872 (2), pp.152. 10.3847/1538-4357/ab0080. hal-02393372

HAL Id: hal-02393372

https://hal.science/hal-02393372v1

Submitted on 12 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The Streams of the Gaping Abyss: A Population of Entangled Stellar Streams Surrounding the Inner Galaxy

Rodrigo A. Ibata<fnref n="1" />, ∗, Khyati Malhan<fnref n="1" />,<fnref n="2" />, and Nicolas F. Martin<fnref n="1" />,<fnref n="3" /> †

Received 2018 September 25; revised 2018 December 24; accepted 2019 January 19; published 2019 February 20

## Abstract

We present the discovery of a large population of stellar streams that surround the inner Galaxy, found in the Gaia DR2 catalog using the new STREAMFINDER algorithm. Here we focus on the properties of eight new high-significance structures found at heliocentric distances between 1 and 10 kpc and at Galactic latitudes $|b| > 20^{\circ}$, named Slidr, Sylgr, Ylgr, Fimbulthul, Svöl, Fjörm, Gjöll, and Leiptr. Spectroscopic measurements of seven of the streams confirm the detections, which are based on Gaia astrometry and photometry alone, and show that these streams are predominantly metal-poor. The sample possesses diverse orbital properties, although most of the streams appear to be debris of inner-halo globular clusters. Many more candidate streams are visible in our maps but require follow-up spectroscopy to confirm their nature. We also explain in detail the workings of the algorithm and gauge the incidence of false detections by running the algorithm on a smooth model of the Gaia catalog.

**Key words:** galaxies: formation – Galaxy: halo – Galaxy: stellar content – Galaxy: structure – surveys

## 1. Introduction

The advent of the second data release (DR2) of the Gaia mission (Gaia Collaboration et al. 2018a; Lindegren et al. 2018) has opened up a new panoramic window onto our Galaxy. The detailed astrometric measurements now show the motions of over a billion stars in our immediate environment, allowing us to deduce the kinematical structure of our host galaxy (Gaia Collaboration et al. 2018b). One of the main aims of this endeavor is to probe the formation history of the Milky Way, by identifying the structures that formed in situ and those that were incorporated into it during the many galactic merging events that are expected to have contributed to its present ensemble.

One class of remnant structures that are “fossils” of the Galactic formation process are the so-called stellar streams, which are the remains of satellites that have been disrupted by the tidal forces of their host galaxy. They provide a very exciting new avenue to assess the distribution of the dark matter, by allowing us to directly measure accelerations. The basic reason for this is that a low-mass stream follows closely its progenitor’s orbit, so a stream locus approximates an orbit in the Galaxy. Thus, mapping the line-of-sight velocity and proper-motion gradients along a stream allows us to uncover the acceleration that the stream is subject to. With several streams on different orbits it should be possible to build up a reliable 3D map of the acceleration field throughout the region of the Galaxy traversed by such streams. These ideas can be extended to higher-mass streams, by accounting for self-gravity and the consequent offset between the stream and the orbit of the progenitor (Varghese et al. 2011).

Streams can also reveal the small-scale granularity of the dark matter (Ibata et al. 2002; Johnston et al. 2002), which is an issue of fundamental importance since one of the central predictions of ΛCDM cosmology is the existence of thousands of dark matter subhalos in Milky-Way-sized galaxies (Klypin et al. 1999; Moore et al. 1999). The dark matter substructures heat the tidal tail, broadening the stream and also increasing the

dispersion of the stream in the space of the integrals of motion. Gaps in the stream may also be induced by the close passage of these dark perturbers (Carlberg 2012; Erkal et al. 2016; Price-Whelan & Bonaca 2018).

The ideal tidal streams to search for dynamical heating and gaps are those that are old (and hence long), since these have the greatest cross section for impacts from ΛCDM substructures and have existed for long enough in the potential to feel many close passages. Indeed, finding a single unheated long (i.e., old) low-mass halo stream would provide a severe challenge to ΛCDM, as it would provide incontrovertible proof that dark matter is smooth on small (subgalactic) scales. Alternatively, obtaining positive proof of the heating action of dark matter lumps would require one either to find streams that do not stray into the Galactic disk, where they may be heated by interactions with baryonic substructures (Amorisco et al. 2016), or to perform a very careful modeling of the stream population as an ensemble and include the effects of nonaxisymmetric perturbations such as that due to the Galactic bar (Hattori et al. 2016; Price-Whelan et al. 2016; Pearson et al. 2017).

The necessary first step toward achieving these important goals of studying the dark matter, as well as the structure, dynamics, composition, and assembly history of our Galaxy, is to find samples of dynamically cold stellar streams. While large photometric sky surveys like the Sloan Digital Sky Survey, Pan-STARRS1, and the Dark Energy Survey have been successfully trawled for these structures (Grillmair & DIONATOS 2006; Grillmair 2009; Bernard et al. 2016; Shipp et al. 2018), there is a clear advantage to now try to search for stellar streams in the Gaia survey given its exquisite kinematic measurements. However, the streams of interest will be mixed in a complex fashion with the Milky Way’s field populations over vast swathes of sky, so disentangling them from each other and from the Galaxy requires careful analysis. It was in view of these opportunities with Gaia (and other new surveys) that we developed an algorithm, which we called the STREAMFINDER, to search for stream-like structures (Malhan & Ibata 2018, hereafter Paper I).

<footnote n="1">¹Observatoire Astronomique, Université de Strasbourg, CNRS, 11, rue de l'Université, F-67000 Strasbourg, France; rodrigo.ibata@astro.unistra.fr</footnote>

<footnote n="2">²The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, AlbaNova, 10691 Stockholm, Sweden</footnote>

<footnote n="3">³Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>
In a previous contribution (Malhan et al. 2018, hereafter Paper II), we applied an updated version of the algorithm of Paper I, which will be presented below, to undertake a first search for tidal streams in the Gaia DR2 catalog. Searching the high-latitude sky at $|b| > 30^\circ$ and at heliocentric distances $d > 5$ kpc, we reported the detection of five new streams (named Gaia 1–5).

The aim of the present paper is to continue this search, extending the sky to $|b| > 20^\circ$ and to distances $d > 1.0$ kpc. An early result of the present effort was the detection of the 75°-long “Phlegethon” stream, which was presented in Ibata et al. (2018, hereafter Paper III).

The Gaia DR2 survey is the principal data set used here. We extinction-corrected this survey using the Schlafly & Finkbeiner (2011) corrections to the Schlegel et al. (1998) extinction maps, assuming the extinction ratios $A_G/A_V = 0.85926$, $A_{GBP}/A_V = 1.06794$, and $A_{GRP}/A_V = 0.65199$, as listed on the web interface to the PARSEC isochrones (Bressan et al. 2012). Henceforth all magnitudes and colors refer to these extinction-corrected values. We also excised all stars in the Gaia DR2 catalog that lie close to certain compact structures. These included objects within two tidal radii of the globular clusters listed in Harris (2010), objects within 7 half-light radii of the Galactic satellites listed in McConnachie (2012), and objects within a $3^\circ$ radius of M31 and $1^\circ$ radius of M33. We also remove stars in appropriately sized regions (typically $\sim 0^\circ 5$) around the open clusters NGC 188, Berkeley 8, NGC 2204, NGC 2243, NGC 2266, Melotte 66, NGC 2420, NGC 2682, and NGC 6939. This pruned catalog is used both in the search for streams below and in the construction of the contamination model.

The layout of the paper is as follows. In Section 2 we explain the new algorithm, paying particular attention to how it differs from the procedure described in Paper I. Section 3 describes the artificial data that were processed to establish the incidence of false-positive detections. The maps of the detected streams are presented in Section 4, and then we describe each new high-significance stream individually in Section 5. In Section 6, we discuss these results and draw the conclusions of our study.

## 2. STREAMFINDER Improvements

Although the original STREAMFINDER concept of Paper I worked well at high Galactic latitude, in regions of relatively uniform contamination, its shortcomings became evident when probing closer to the Galactic plane. There, the strong gradient in the background produced an unacceptably large population of false positives. We decided therefore to redesign the algorithm to undertake a full likelihood analysis, as described below.

Conceptually, we proceed in a similar manner to Paper I, by considering, in turn, each star *i* in the survey down to a chosen magnitude limit. The question we aim to answer is, if star *i* forms part of a stream, does the maximum likelihood stream solution, centered on *i*, imply that the putative stream is significant? To this end, we search for the most likely stream model $P_{stream}(\theta)$ and the corresponding stream fraction $\eta$ that maximize

$$ \ln \mathcal{L} = \sum_{\text{data}} \ln [\eta P_{\text{stream}}(\theta) + (1 - \eta) P_{\text{cont}}], \quad (1) $$

where $\theta$ are the stream-fitting parameters and the probability density function $P_{cont}$ is a model of the “contamination”

(sometimes called a “background”) from nonstream stars. By comparing to the likelihood of the no-stream case (when $\eta = 0$), Equation (1) easily allows us to measure the stream detection significance. In the analysis below, the vast majority of the Gaia survey stars are found to have $\eta = 0$ (i.e., the most likely stream solution is to have no stream at all).

### 2.1. Stream Model

To build the stream model $P_{stream}(\theta)$, we use a simple leapfrog scheme to integrate an orbit both forward and backward from the star under consideration, under the influence of an acceleration field provided by the realistic Galactic density distribution of Dehnen & Binney (1998, their model 1, which is axisymmetric and contains a bulge, disk, thick disk, interstellar medium, and halo components). The orbit defines the central line for a proposed stream. This line is projected into the space of the following six observed quantities in Gaia DR2: sky position $\alpha, \delta$, proper motions $\dot{\mu}_\alpha(\equiv\mu_\alpha \cos(\delta))$, $\mu_\delta$, and color-magnitude position $G, G_{BP}-G_{RP}$, as described in Paper I. To make this projection, we adopt the recent accurate determinations of the galactocentric distance of $R_\odot = 8.122 \pm 0.031$ kpc (Gravity Collaboration et al. 2018) and of the circular velocity at the solar neighborhood of $v_c(R_\odot) = 229.0 \pm 0.2$ km s$^{-1}$ (Eilers et al. 2019). Given that $v_c(R_\odot) + V_\odot = 255.2 \pm 5.1$ km s$^{-1}$ (Reid et al. 2014), we take the V component of the peculiar velocity of the Sun to be $V_\odot = 26.2$ km s$^{-1}$ (which we keep fixed). The U and W components of the Sun’s peculiar velocity are taken from Schönrich et al. (2010). Henceforth, for convenience, we will drop the asterisk superscript from $\mu_\alpha$.

With Gaia DR2 data we have access to very precise sky positions and superb, but still noisy, proper motions. The line-of-sight velocities $v_{los}$ of the survey stars are effectively all unknown. For distant halo stars, of which the majority are faint in Gaia (see Ibata et al. 2017, their Figure 1), parallaxes are effectively uninformative, yet the distances can be constrained via the excellent Gaia photometry, given that plausible populations of stream stars should follow color-magnitude tracks and luminosity functions of plausible progenitors. For this purpose we choose the PARSEC single stellar population (SSP) models (Bressan et al. 2012) of different metallicities and ages to serve as templates. Nevertheless, a given model may yield multiple distance solutions (as we discuss in Paper I), while SSP models of different metallicity or age will also yield different distance solutions.

The stream model is simply a smeared-out version of the orbit, blurred with Gaussian widths in each observed parameter in such a way as to display similar sky width, depth, velocity distribution, and color-magnitude dispersion as known streams (or as the type of stream one wishes to detect). For the present work, we adopted a stream length of $L = 20^\circ$. For each star, we find the minimum angular distance from the stream and calculate the corresponding perpendicular distance $s$, assuming that the star has the same line-of-sight distance $d$ as the orbit. We adopt a Gaussian stream model width of $w_s = 100$ pc (we use $w$ to denote model Gaussian dispersions and $\sigma$ to denote measurement uncertainties) and a velocity dispersion of $w_v = 3$ km s$^{-1}$ in each kinematic dimension (which is converted into a model dispersion in proper motion $w_\mu$ given $d$).

The probability density model of the stream is then

$$ P_{\text{stream}}(\theta) = P_{\text{length}} \times P_{\text{width}} \times P_{\text{LF}} \times P_{\text{color}} \times P_{\mu}, \quad (2) $$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Example of a slice through the “fine” spatial probability density function $P_{\text{fine}}(\alpha, \delta, G, G_{\text{BP}}-G_{\text{RP}})$. The color–magnitude slice is centered at $G_{\text{BP}}-G_{\text{RP}} = 1.1$, $G = 18.1$, with bin width of 0.05 mag × 0.10 mag, respectively. The panels show a zenithal equal-area projection of the counts (in linear scale) centered on the north Galactic pole (left) and south Galactic pole (right). These probability density distributions are constructed to follow the Galactic populations in a smooth manner.</figcaption></figure>

where $P_{\text{length}} = 1/L$ is the (uniform) probability of lying at the measured position along the stream and $P_{\text{width}} = N(\Delta s, w_\delta)$ is the Gaussian probability of lying at the measured distance $\Delta s$ perpendicular to the stream. $P_{\text{LF}}$ is the probability that the observed star has been drawn from the luminosity function of the adopted SSP model (placed at the stream distance $d$). $P_{\text{color}} = N(\Delta(G_{\text{BP}}-G_{\text{RP}}), \sigma_{G_{\text{BP}}-G_{\text{RP}}})$ is the Gaussian probability that the measured color of the star is consistent with that expected of the SSP model (again at distance $d$) given the measured magnitude $G_0$ and the color uncertainty $\sigma_{G_{\text{BP}}-G_{\text{RP}}}$. Finally,

$$ P_{\mu} = \frac{1}{2\pi\sigma_{\mu_{\alpha}}\sigma_{\mu_{\delta}}\sqrt{1 - \rho^2}} \times \exp\left(-\frac{1}{2(1 - \rho^2)}\left[\frac{\Delta_{\mu_{\alpha}}^2}{\sigma_{\mu_{\alpha}}^2} + \frac{\Delta_{\mu_{\delta}}^2}{\sigma_{\mu_{\delta}}^2} - \frac{2\rho\Delta_{\mu_{\alpha}}\Delta_{\mu_{\delta}}}{\sigma_{\mu_{\alpha}}\sigma_{\mu_{\delta}}}\right]\right) \quad (3) $$

is the probability that the measured proper motions agree with those of the closest point on the orbit model (with measured offsets from the model of $\Delta_{\mu_{\alpha}}$ and $\Delta_{\mu_{\delta}}$). The Gaia proper-motion uncertainties $\sigma_{\mu_{\alpha}}$ and $\sigma_{\mu_{\delta}}$ are thus taken into account along with their measured correlation $C \equiv \text{pmra}_{\text{pmdec\_corr}}$ (discussed in Lindegren et al. 2018), which is subsumed into the parameter

$$ \rho = \frac{C \sigma_{\mu_{\alpha}} \sigma_{\mu_{\delta}}}{\sqrt{(\sigma_{\mu_{\alpha}}^2 + w_{\mu}^2)(\sigma_{\mu_{\delta}}^2 + w_{\mu}^2)}} \quad (4) $$

This last equation can be derived by considering the convolution of the 2D covariance matrix with an isotropic 2D Gaussian of dispersion $w_\mu$. This allows us to incorporate the effect of the model dispersion. While we have included the measured correlation between $\mu_\alpha$ and $\mu_\delta$, we ignore the other elements of the (5 × 5) astrometric covariance matrix, as they are much less important: correlations with sky position will be negligible given the size of $w_s$, and we can completely ignore

the correlations of proper motions with parallax, as the latter is not used in the likelihood calculation.

## 2.2. Contamination Model

To be usable in Equation (1), the contamination model $P_{\text{cont}}$ also needs to be a probability density function of the same six observables $(\alpha, \delta, \mu_\alpha, \mu_\delta, G, G_{\text{BP}}-G_{\text{RP}})$ as $P_{\text{stream}}$. By “contamination” what we mean here, of course, are the components of the Milky Way that are not stream-like. We could have used a Galactic synthesis model, such as the Gaia Universe Model Snapshot (GUMS; Robin et al. 2012), for this purpose, but we decided instead to build an empirical contamination model that we hoped would better trace the complexities of our Galaxy.

We smoothed out spatially the Gaia catalog by redrawing each survey star 1000 times from a 2D Gaussian in tangent-point coordinates using a dispersion of 2°. The resulting catalog was then binned into a 4D grid in position and color-magnitude. We adopted a spatial binning scheme in zenithal equal-area Galactic polar projection (the projection used in the maps in Paper II) with a bin size of 1°4 × 1°4. The color-magnitude intervals were chosen to be 0.05 × 0.10 mag. This binned array allows us to construct a probability density function in position and color-magnitude; since this will provide the fine spatial sampling of the contaminating population, let us call this $P_{\text{fine}}(\alpha, \delta, G, G_{\text{BP}}-G_{\text{RP}})$. We present an example of a slice through this probability density function in Figure 1.

However, what we need to model are the correlations between color-magnitude and proper motion as a function of position on the sky. To this end we employed the Gaussian Mixture Model (GMM) decomposition implemented in the Armadillo software package (Sanderson & Curtin 2016). The GMM representation needs to have sufficient flexibility to reproduce the complex color-magnitude behavior of Galactic populations, while also following their proper-motion behavior. By experimenting with the GMM software, we found that we could obtain a suitably smooth fit in the four dimensions of color-magnitude and proper motion with 100 GMM components, including cross-terms. To mitigate against varying depth due mainly to extinction, we limited the GMM fits to $G = 20$ mag.

In Figure 2 we show (for a representative high-latitude field) the observed proper-motion distribution extracted from a small color-magnitude box compared to the smooth GMM representation fit to the full color-magnitude sample. The model contours can be seen to provide a reasonable representation of the proper-motion distribution. In Figure 3 we have marginalized over proper motion to show the full color-magnitude distribution in the same spatial region as Figure 2. Although it is obviously challenging to reproduce the color-magnitude behavior of stars in a typical Galactic field with Gaussian components, the software again achieves a visually acceptable representation.

For the GMM modeling we found it necessary to employ 4 × 4 larger spatial pixels of size 5°6 × 5°6, so that the spatial regions would contain sufficient stars (>50,000 stars at the Galactic poles) for the GMM algorithm to yield smooth results. We label the probability density function derived in this way $P_{\text{GMM}}(\mu_\alpha, \mu_\delta, G, G_{\text{BP}}-G_{\text{RP}}|\alpha, \delta)$, which is continuous and smooth in color-magnitude and proper motion by construction.

For each star, we cut through the GMM model at the value of its photometric $G, G_{\text{BP}}-G_{\text{RP}}$ measurements, to obtain the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Example of a slice through the GMM derived for a large 5°6 × 5°6 pixel centered at $ℓ = 105°2, b = -60°8$. The contours show the GMM model summed over the color–magnitude range $0.7 < G_{BP}-G_{RP} < 0.9, 16 < G < 17$, and the circles show the data over the same range.</figcaption></figure>

conditional probability $P_{GMM}(\mu_{\alpha}, \mu_{\delta}|\alpha, \delta, G, G_{BP}-G_{RP})$. Clearly, this means that we have to normalize the function so that $\iint P_{GMM}d\mu_{\alpha}d\mu_{\delta} = 1$.

The final contamination model is then

$$ \begin{aligned} P_{\text{cont}}(\alpha, \delta, \mu_{\alpha}, \mu_{\delta}, G, G_{\text{BP}}-G_{\text{RP}}) \\ = P_{\text{fine}}(\alpha, \delta, G, G_{\text{BP}}-G_{\text{RP}}) \\ \times P_{\text{GMM}}(\mu_{\alpha}, \mu_{\delta}|\alpha, \delta, G, G_{\text{BP}}-G_{\text{RP}}). \end{aligned} \quad (5) $$

This procedure may seem convoluted to the reader, but it has allowed us to build a Galaxy model that follows the density of the smooth populations with a relatively fine spatial and color–magnitude mesh, while simultaneously following the behavior of the populations in proper motion. In contrast, our experiments of using the GMM method in six dimensions were not at all satisfactory: we could not reproduce the complex structure in this 6D space and also obtain a smooth model.

This contamination model passed hundreds of visual spot-checks, in which we compared the observations to the model at different locations on the sky in different 2D slices through the G, $G_{BP}-G_{RP}$, $\mu_{\alpha}$, $\mu_{\delta}$ parameter space. In these visual checks we ensured that the model was not underfitting the data; we found occasionally that with ~50 GMM components or less the model would only reproduce the disk dwarf sequence (the feature at $G_{BP}-G_{RP} \gtrsim 1.5$ in Figure 3). We also ensured that the model was not overfitting the data, which would be visible as a lumpy model in color–magnitude space or in proper-motion space. Despite these checks, we cannot of course be certain that our very complex contamination model provides a suitable representation of the Galaxy in this multidimensional parameter space. A particular worry is that artifacts in the contamination model could give rise to false-positive detections. To address this issue, in Section 3 below we will take recourse in a simulation of the Galaxy, and later we will show

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Example of the color–magnitude behavior of the GMM. We show the GMM model (contours) for the same sky region as shown previously in Figure 2, marginalized over proper motion, along with the data (blue circles) that were used to build it. The model is a 4D function built out of 100 Gaussian components, yet it manages to represent the data in an apparently smooth way. In all of our hundreds of visual inspections of the GMM model over the sky a similar smooth behavior was seen.</figcaption></figure>

that the real Gaia data show stream-like structures unlike those found in the smooth simulation.

## 2.3. Parameter Fitting

Within a given run of the STREAMFINDER algorithm the free parameters are thus $\theta = [v_{\text{los}}^{\text{model}}, \mu_{\alpha}^{\text{model}}, \mu_{\delta}^{\text{model}}]$, and the stream fraction parameter $\eta$. The distance degeneracy is explored by rerunning the algorithm with different plausible SSP models. We find the maximum value of $\ln L$ by first sampling over $v_{\text{los}}^{\text{model}}, \mu_{\alpha}^{\text{model}},$ and $\mu_{\delta}^{\text{model}}$. Here, just as in Paper I, we sample $v_{\text{los}}^{\text{model}}$ over the escape velocity ($\pm 600 \text{ km s}^{-1}$) uniformly in steps of $20 \text{ km s}^{-1}$, and we also sample $\pm 1\sigma$ in each of the proper-motion directions in three uniform steps. For a given proposed $\theta$, we hunt for the value of $\eta$ that yields the peak likelihood, using Brent's optimization method (Press et al. 1992). The best trial solutions thus obtained are provided as starting input for a “downhill simplex” routine (Press et al. 1992) that refines $\theta$ and $\eta$ to maximize $\ln L$, iterating until convergence is found to an absolute tolerance of $10^{-4}$ in $\ln L$.

It can be appreciated from the above discussion that this is a computationally very intensive task, since at the location of each star we calculate many thousands of orbits and calculate the stream likelihood by summing over all stars in the survey; the process being repeated until every star in the survey has been examined in this way.

We decided not to use the Gaia parallaxes in the likelihood calculation since the algorithm was initially intended to detect distant halo streams for which we expected the parallax measurements to be too poor to be useful. In retrospect this was probably a mistake, which we may fix in future versions of the software. However, it will likely be challenging to implement
<figure><img src="image_6.png" /><figcaption>Figure 4. Ratio of the number of stars in the Gaia catalog compared to the GUMS simulation.</figcaption></figure>

this upgrade, as it will require us to include another dimension in the contamination model, but one that is correlated with color and magnitude. In the meantime, although the parallaxes are not directly included in the likelihood calculation, we still use the measurements by simply filtering the STREAMFINDER output to those solutions where the model distances are coherent with the measured Gaia parallaxes to within 2σ. Although the heliocentric velocities of the streams are calculated internally within the algorithm, we do not yet include this additional observable dimension in the stream model; we intend to rectify this for the Gaia DR3 release, which is expected to include the radial velocities of stars down to a limiting magnitude of G ~ 16.5.

### 3. Comparison Data

In Section 4 below we will present the maps of stream-like structures found by the STREAMFINDER algorithm, in which many unexpected features are detected. It is natural to wonder whether these apparent streams are real or false positives created by our method of hunting for linear overdensities along orbits. In order to assess these concerns, we decided to also apply the algorithm to the GUMS (Robin et al. 2012), which provides a realistic simulation of the stellar density and kinematics of the smooth components of the Milky Way.

The GUMS simulation does not include an estimate of the observational uncertainties, however. So as to model the Gaia DR2 observations as closely as possible, we chose to directly reuse the uncertainties in the real data. This was implemented as follows. For every artificial star in GUMS, we searched for the spatially closest real counterpart in Gaia within a 1° radius with a G-band magnitude within 0.1 mag of the GUMS star. If no such real star was found, the search was repeated doubling the magnitude difference criterion (repeated, if necessary, until a match was found). The proper-motion and parallax uncertainties of the matched Gaia star were then assigned to the simulated GUMS star, together with the observed proper-motion correlation, and these Gaussian uncertainties were used to resample the GUMS proper motions and parallaxes, thus making the properties of these parameters closer to the observations than in the original GUMS simulation. We apply exactly the same mask to remove stars around the locations of the Galactic satellites and globular clusters (so as to generate the same gaps in the survey, even though these satellites are not present in the simulation).

Down to $G_0 = 20$, and at the Galactic latitudes of interest to the present study ($|b| \gtrsim 20°$), the number density of stars in the GUMS simulation is significantly higher than in the Gaia DR2 catalog. The Gaia/GUMS ratio is shown in Figure 4 and

averages 0.6 for the sky at $|b| > 20°$. The reason for the overestimate in the GUMS model is not clear to us, but we believe that the discrepancy is not problematic for the present purposes, as long as we correct for it. To that end, we take the same $1.4 \times 1.4$ spatial bins as before (used in Figure 4) and simply randomly discard GUMS stars until we obtain the number density observed in Gaia DR2. This procedure has the added advantage of imprinting the effects of the actual Gaia scanning law on the simulated data, rendering it more realistic.

We then regenerate the contamination model $P_{cont}$ directly from the simulated data.

### 4. Results

The algorithm was run over metallicity values of $[Fe/H] = \{-2.0, -1.6, -1.4\}$ dex, and we integrated orbits starting from Gaia stars with magnitudes $G_0 < 19.5$ and $|b| > 20°$. The magnitude cut at $G_0 = 19.5$ was chosen so as to avoid edge effects stemming from the faint limit of our GMM model at $G_0 = 20.0$ (see Figure 3). A single SSP model age of 12.5 Gyr was adopted. Since the allowed half-length of the orbits was set to $10°$, this meant that we actually examined data at all latitudes $|b| > 10°$. The extension to lower latitude will be presented in future contributions, but this is computationally extremely expensive owing to the high source density (the analysis presented here already took over a million CPU hours).

Finally, following Lindegren et al. (2018), we filtered the output stars to retain only those having flux excess $E(\equiv phot_bp RP_excess_factor)$ in the range

$$1 + 0.015(G_{BP} - G_{RP})^2 < E < 1.3 + 0.06(G_{BP} - G_{RP})^2. \quad (6)$$

This removes stars that have suspicious photometry, which amount to ≈10% of the sources in streams. We note, however, that the maps are qualitatively identical with or without this filter.

Inspection of the STREAMFINDER output shows a clear dichotomy between the stream-like features with low proper motion ($\sqrt{\mu_e^2 + \mu_b^2} < 2$ mas yr$^{-1}$) and those of higher proper motion. The solutions with low proper motion form a relatively coherent pattern on the sky, as we show in Figure 5, which displays the stream-like features with $\sqrt{\mu_e^2 + \mu_b^2} < 2$ mas yr$^{-1}$ that are identified by the algorithm at distances between 1 and 10 kpc when assuming a model of metallicity $[Fe/H] = -2.0$ (other choices of the template metallicity result in similar sky distributions). The distance limits set the range that the algorithm is allowed to search over, and having only a narrow distance interval such as this helps speed up the computation. The top two rows of panels show the proper-motion distributions in $\mu_e$ and $\mu_b$ (deliberately using the same color scale as will be used later in Figures 6–8), the third row shows the heliocentric distance distribution, while the bottom row shows the z-component of angular momentum $L_z$ of the stream solutions. We have set the detection threshold here, and in all subsequent maps, to $8\sigma$. These features are not present in the equivalent maps based on the GUMS comparison data, as we show in Figure 20 in the Appendix, and so they do not appear to be an obvious artifact of the method. Nevertheless, these structures are not streams of low-mass accretions that we aim to identify here and constitute a contaminating population for the present work (we defer their analysis to a subsequent follow-up study). The bottom row of panels shows that they possess
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Spatial distribution of stars that possess low proper motion ($\sqrt{\mu_{\ell}^2 + \mu_b^2} < 2$ mas yr$^{-1}$) and that are identified as stream-like with $>8\sigma$ confidence in the distance interval 1 kpc $< d_{\odot} < 10$ kpc, assuming a metal-poor template with age = 12.5 Gyr and [Fe/H] = -2.0. The first and second rows of panels show the distribution of proper motion in $\mu_{\ell}$ and $\mu_b$, respectively; the third row of panels shows the distance solution calculated by the algorithm with this stellar population template, and the bottom row of panels shows the z-component of angular momentum ($L_z$) of the best-fit orbit. The first and second columns of panels show zenithal equal-area projections centered on the north and south Galactic poles, respectively.</figcaption></figure>

negative $L_z$ (meaning that they are prograde), and we found that we could eliminate this population by setting a higher detection threshold for prograde stars. However, in the following we choose to simply filter out the stars for which the most likely stream solutions have $\sqrt{\mu_{\ell}^2 + \mu_b^2} < 2$ mas yr$^{-1}$, as we can then present maps with a uniform detection threshold that does not depend on $L_z$.

The complement to the sample in Figure 5 is shown in Figure 6, which displays the stream solutions with $\sqrt{\mu_{\ell}^2 + \mu_b^2} > 2$ mas yr$^{-1}$ in the distance range between 1 and 10 kpc derived using the SSP template with age 12.5 Gyr and metallicity [Fe/H] = -2.0. Figures 7 and 8 show similar information but using a template with metallicity of [Fe/H] = -1.6 and [Fe/H] = -1.4, respectively. The layout
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Same as Figure 5, but showing the stars with proper motion $\sqrt{\mu_{\ell}^2 + \mu_{b}^2} > 2$ mas yr$^{-1}$ that are identified as likely stream members using a >8σ confidence threshold. In addition to the Phlegethon stream reported in Paper III, we detect here the structures labeled “Slidr,” “Slygr,” “Ylgr,” “Fimbulthul,” “Svöl,” “Fjörm,” “Gjöll,” and “Leiptr.”</figcaption></figure>

of the panels is identical to those in Figure 5, and the color-coding is identical as well.

The well-studied GD-1 stream (Grillmair & Dionatos 2006) is detected at high confidence (up to 20σ in portions of the stream) in all three metallicity maps, although its contrast decreases when using the higher-metallicity templates. We note that the imposed upper distance limit of 10 kpc artificially limits this stream to $\ell \gtrsim 110^\circ$. The only other previously detected stream present in our maps is the Jhelum structure

(Shipp et al. 2018). The detection of GD-1 and Jhelum shows that the algorithm works as intended and is able to recover low surface brightness streams.

A large number of additional stream features are detected by the algorithm at this 8σ threshold. We have marked and named 10 other streams that are distinct and coherent in all of the observed parameters and clearly distinguishable from possible residual features of the Galactic disk or bulge. In addition to the “Phlegethon” stream, detected by our algorithm and reported in
<figure><img src="image_9.png" /><figcaption>Figure 7. Same as Figure 6, but for stars with [Fe/H] = -1.6.</figcaption></figure>

Paper III, the new detections are “Slidr,” “Sylgr,” “Ylgr,” “Fimulthul,” “Svöl,” “Fjörm,” “Gjöll,” and “Leiptr.” These names are taken from Norse mythology, representing 8 of the 11 rivers that existed in the void (or gaping abyss) at the beginning of the world, with Slidr and Gjöll also being streams of the underworld (as is the Phlegethon in Greek mythology).

In Figures 21–23 in the Appendix, we reproduce exactly the same maps as for Figures 6–8, but for stream-like detections derived from the GUMS simulation. These false positives are found toward the bulge region of the simulation in areas of very

high source density. Their distribution in proper motion is markedly different from the detections we report here.

A number of other stream candidates can also be seen in these maps, especially in Figure 8. They are, however, of somewhat lower significance than the eight new obvious structures we have named, or they have kinematic properties that make them not so easily distinguishable from the false positives in our processing of the GUMS catalog.

We now proceed to examine each of the eight streams separately, first the streams in the northern Galactic hemisphere
<figure><img src="image_10.png" /><figcaption>Figure 8. Same as Figure 6, but for stars with [Fe/H] = -1.4.</figcaption></figure>

in order of increasing longitude in Figure 6, and then the streams in the southern Galactic hemisphere, also in order of increasing longitude.

# 5. The Streams

## 5.1. Slidr

The Slidr stream can be perceived as a 32°-long prograde structure in the northern hemisphere maps. Although it is detected with all of the chosen metallicity templates, it appears strongest in the [Fe/H] = -1.6 map (Figure 7). In Figure 9 we show the parameter correlations of stars of this structure that

the algorithm identifies; these Gaia DR2 stars are displayed as red circles with their corresponding uncertainties. A total of 156 stars are detected at >8σ confidence.

In order to obtain a global fit to these data, we reused the procedure developed in Paper III to fit the “Phlegethon” stream. This refitting is necessary because the orbit fits computed by the STREAMFINDER are optimized for speed rather than accuracy and reliability (for instance, by using the above-mentioned “downhill simplex” method, its solutions can get caught in local likelihood maxima), and furthermore the search routine only has access to data over the chosen stream search length L. The dedicated orbit-fitting software implements a
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Properties of the Slidr stream. From top to bottom, the panels show proper motion $\mu_\alpha$ and $\mu_\delta$, parallax $\varpi$, heliocentric line-of-sight velocity $v_{\text{helioc}}$, and Galactic latitude position of the stream stars, as a function of Galactic longitude $\ell$. The red circles are the 156 members of the Gaia DR2 catalog with five-component astrometric solutions that are identified by the STREAMFINDER software as stream stars using an $8\sigma$ threshold. The uncertainties on these points are shown with error bars. The large circles in the fourth panel mark those stream members whose radial velocity has been measured by the SDSS, LAMOST, or Gaia RVS surveys. As explained in the text, we fitted an orbit model to the Gaia data, without using any of the radial velocity information. The corresponding best-fit model is shown with a dashed line and can be seen to give a reasonable representation of the parameter profiles. The radial velocity profile is effectively a prediction of the model and can be seen to pass close to all but one of the metal-poor stars with [Fe/H] < -1.6, confirming the reality of the structure. Additionally, the bottom panel also shows the $E(B-V)$ extinction map, extracted from Schlegel et al. (1998), from which the reader may verify that the detected stream does not follow any obvious structure in the interstellar extinction.</figcaption></figure>

Markov chain Monte Carlo (MCMC) exploration for the initial conditions of an orbit that best fits the data in position, proper motion, and parallax. The orbit integrator and Galactic potential model are identical to those used in the STREAMFINDER, and the parameter space is sampled using the affine-invariant “walkers” approach of Goodman & Weare (2010). Contrary to Paper III, however, we did not attempt to fit any radial velocity measurements, partially because we do not have this information for all of the streams presented here, but more importantly because we desired to obtain predictions for the radial velocity gradient based on the Gaia astrometry alone.

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Same as Figure 9, but for the Sylgr stream. A total of 103 stars are identified as candidate members of this structure at >8$\sigma$ confidence.</figcaption></figure>

The algorithm maximizes the likelihood of the orbit given the data, adopting the following probability density function for the stream:

$$P_{\text{fit}}(\theta) = P_{\text{width}} \times P_{\mu} \times P_{\varpi}. \quad (7)$$

The terms $P_{\text{width}}$ and $P_{\mu}$ here are identical to those defined above in Section 2.1 for the STREAMFINDER. The last term is a simple Gaussian probability density function for the observed parallax $P_{\varpi} = N(\Delta\varpi, \sigma_{\varpi})$, where we have ignored the intrinsic parallax dispersion, since the observational uncertainties completely dominate the dispersion of the data about the model. We ran the MCMC procedure for $1.1 \times 10^6$ iterations and rejected the first $10^5$ (burn-in) steps. The best-fit model is shown as a dashed line in Figure 9 and can be seen to give an excellent representation of the proper motions in $\mu_\alpha$ and $\mu_\delta$ (top two panels), the parallax (middle panel), and the position on the sky (bottom panel). The model distance to the structure (at $\alpha = 171.043, \delta = 6^\circ$) is $d = 3.65 \pm 0.09$ kpc.

In addition, the best-fit orbit gives a prediction of the heliocentric line-of-sight velocity of the structure (shown in the fourth panel). By chance, our Slidr sample contains 10 stars previously observed in spectroscopy by the Sloan Digital Sky Survey (SDSS) DR10 (Yanny et al. 2009; large blue circles), 14 stars observed by LAMOST DR4 (Cui et al. 2012, large
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Same as Figure 9, but for the Ylgr stream, which contains 349 candidate members.</figcaption></figure>

magenta circles), and two stars observed by the Gaia Radial Velocity Spectrometer (RVS; large orange circles). The candidate stream members with radial velocity measurements are listed in Table 1. The two Gaia DR2 RVS stars happen to be a subset of the LAMOST sample. There is a clear separation in kinematic behavior between those stars that have spectroscopically measured metallicities above and below [Fe/H] = −1.6. We display the position of stars that are more metal-poor than this limit with large filled circles, while the more metal-rich stars are marked by squares. Apart from one outlier, all of the metal-poor subsample follows the predicted line-of-sight velocity gradient. The mean metallicity of the six stream member stars in the SDSS sample is [Fe/H] = −1.80, essentially identical to the mean metallicity of [Fe/H] = −1.84 of the eight LAMOST stars that are stream members.

This radial velocity confirmation firmly establishes the
reality of the Slidr stream and proves that the algorithm is
able to uncover very low contrast structures that have been
missed in previous searches. The presence of kinematic outliers
here provides a reminder that the STREAMFINDER produces
samples of stream candidates. The 8σ selection threshold
indicates that, at the position of each identified star, there is a
stream present at >8σ confidence that shares the kinematic
properties of the star (and is consistent with its photometric

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Same as Figure 9, but for the Fimbulthul stream. This 18°-long structure contains 309 candidate members.</figcaption></figure>

properties and parallax). In particular, it does not mean that
each stream candidate star identified by the STREAMFINDER
is a member with >8σ confidence of the detected stream.

The orbit of Slidr is prograde but highly radial, with a
pericenter at 4.45 ± 0.08 kpc and apocenter at 42.4 ± 3.2 kpc.
Its z-component of angular momentum is L<sub>z</sub> = −1458 ±
9 km s<sup>−1</sup> kpc, according to the MCMC orbital fit (we adopt
the usual convention that negative L<sub>z</sub> means prograde motion).
Its orbit is shown in blue in Figure 18, and its position in the
L<sub>z</sub>—pericentric distance plane is displayed in Figure 19. These
orbital parameters are listed in Table 2, for easy comparison to
the other streams.

The effect of sample contaminants on the derived orbits is
hard to gauge. Here and in most cases below, the fact that the
orbit as derived without radial velocities agrees with the small
set of stars whose radial velocities have been measured
suggests that the presence of contaminants is not problematic.
Ultimately, however, the candidate stream stars identified by
the STREAMFINDER should be followed up with spectroscopy
to eliminate any doubts and fully determine the orbit.

5.2. Sylgr

The Sylgr stream is a 24°-long nearby structure that appears
to have highest contrast in the [Fe/H] = −1.6 map (Figure 7),
where it can be easily identified as a stream of 103 stars with
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Same as Figure 9, but for the Svöl stream, of which the algorithm finds 45 candidate members.</figcaption></figure>

highly negative $\mu_b$ in the northern hemisphere. The observed parameter gradients are shown in Figure 10. The result of the MCMC orbit-fitting procedure is shown with a dashed line and gives a reasonable representation of the structure. Of the members we identify with the STREAMFINDER, three have measured radial velocities in the SDSS; these are shown on the fourth panel of Figure 10. As before, we did not use the radial velocity information in the orbit-fitting procedure. All three of the SDSS stars lie close to the predicted radial velocity gradient of the stream, and they possess extremely low metallicities ($\langle[\text{Fe/H}]\rangle = -2.7$). It will be interesting to rerun the detection software in the region around this structure using a very metal-poor template.

The best-fit orbit shows that the structure is prograde ($L_z = -686 \pm 7 \text{ km s}^{-1}$ kpc), with a pericenter at $2.49 \pm 0.02$ kpc and an apocenter at $19.5 \pm 0.5$ kpc. The color-magnitude diagram (CMD) and orbital properties of the structure are shown with orange circles in Figures 17–19.

### 5.3. Ylgr

Figure 11 shows the properties of the 30°-long Ylgr stream. It is most clearly visible in the $[\text{Fe/H}] = -2.0$ map (Figure 6). The sample we identify contains 349 stars, three of which have been previously observed spectroscopically by the SDSS. All

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Same as Figure 9, but for the Fjörm stream. The 148 candidate members of this structure span an arc of 87° in the northern sky.</figcaption></figure>

three of those stars are metal-poor ($\langle[\text{Fe/H}]\rangle = -1.87$), and the sample displays very small scatter (0.04 dex). The most likely stream fit without radial velocity information (dashed line) does not pass through the observed SDSS radial velocity measurements. If we include the radial velocity information in the likelihood, and assuming that the star with $v_{\text{helio}} = -53.6 \text{ km s}^{-1}$ is a contaminant, we obtain the dotted-line orbital fit. The orbit of the stream is highly retrograde ($L_z = 3087 \pm 36 \text{ km s}^{-1}$ kpc), with a pericenter at $10.3 \pm 0.04$ kpc and an apocenter at $30.6 \pm 0.9$ kpc.

### 5.4. Fimbulthul

Fimbulthul is perhaps the most fascinating of the streams we have identified and will be discussed in greater detail elsewhere. A total of 309 stars over an 18° length of sky are identified as candidate members of this structure by the STREAMFINDER algorithm. Two stars of this sample happen to have been measured by the Gaia RVS. The orbital fit without radial velocity anchors is shown with the dashed line in Figure 12. However, since this fit fails to account for the observed radial velocities, we decided to refit the structure including the RVS values as constraints. The new fit is shown with the dotted line and provides a good description of the structure.
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Same as Figure 9, but for the Gjöll stream. The STREAMFINDER identifies 28 stars as being members of this system.</figcaption></figure>

The revised orbit model has a pericentric distance of 1.50 ± 0.01 kpc and an apocenter distance of 6.53 ± 0.01 kpc and is again highly retrograde ($L_z = 490 \pm 1$ km s⁻¹ kpc). These properties are very close to those of the prototypical massive cluster ω Centauri, which suggests that the Fimbulthul stream is indeed the long-sought-for tidal tail of that globular cluster.

### 5.5. Svöl

Svöl is a 26°-long stream seen in the northern hemisphere and that appears with highest contrast in the more metal-poor maps. The stars of this structure are displayed in Figure 13, to which we again applied our MCMC orbit-fitting procedure. Although one star of the sample has been observed by LAMOST, the relatively high metallicity of this object ($[Fe/H] = -1.09 \pm 0.8$) suggests that it may be a contaminant.

Nevertheless, the orbit based purely on Gaia DR2 data is well constrained: the pericenter lies at 2.53 ± 0.13 kpc, the apocenter is at 19.1 ± 1.1 kpc, and the stream is prograde, with $L_z = -793 \pm 30$ km s⁻¹ kpc. This orbit is shown in purple in Figure 18, and its angular momentum and pericenter are compared to other streams and globular clusters in Figure 19. As can be seen in Figure 19 and more quantitatively in Table 2,

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Same as Figure 9, but for the Leiptr stream, which contains 67 candidates.</figcaption></figure>

the properties of Svöl are extremely close to those of the Sylgr stream, which suggests the possibility that they are two fragments of a single structure. While we cannot rule out this possibility, we were unable to find trajectories that link the two features in less than one orbital period.

### 5.6. Fjörm

Fjörm is the longest of the streams presented here: it spans 87° in the northern Galactic hemisphere and passes close to the north Galactic pole. We identify 148 candidate members with the adopted $8\sigma$ threshold. Seven stars of this sample have spectroscopically measured radial velocities and metallicities in the SDSS, and a further two in the LAMOST survey.

In Figure 14 we show the orbital fit based purely on the Gaia DR2 measurements. The spectroscopically observed stars follow the expected velocity trend very well. All of the SDSS stars are very metal-poor, with an average metallicity of $\langle[Fe/H]\rangle = -2.20$ and sample dispersion of 0.30 dex. Of the two LAMOST stars, one has $[Fe/H] = -2.16 \pm 0.07$, in excellent agreement with the SDSS subsample, while the other star is likely a contaminant since it is metal-rich ($[Fe/H] = -0.70 \pm 0.18$) and is also a velocity outlier.

The best stream fit has a pericenter at 6.93 ± 0.02 kpc and an apocenter at 31.4 ± 0.6 kpc and is prograde with $L_z = -2373 \pm 7$ km s⁻¹ kpc.
<figure><img src="image_19.png" /></figure>

**Figure 17.** CMDs of the eight high-significance streams identified in this work. The Slidr system is shown without a color shift, but each stream thereafter has been shifted by 0.5 mag in color so as to render the distributions visible. The PARSEC isochrone model with metallicity $[Fe/H] = -1.6$ and age of 12.5 Gyr is superimposed, using the distance derived for the Slidr stream (corresponding to a distance modulus of 12.81). We also mark those stars with radial velocity measurements, as indicated in the legend. Panel (a) shows the original CMDs derived from *Gaia* DR2, while panel (b) shows the Pan-STARRS1 DR1 (PS1) photometry (Chambers et al. 2016) of the same stars. The consistency in the appearance of the *Gaia* and PS1 CMDs indicates that our detections are not due to spurious *Gaia* photometry.

## 5.7. Gjöll

At first sight, the Gjöll stream (Figure 15) could perhaps be mistaken to be part of the Galactic Anticenter Stellar Structure (GASS), the broadly distributed stellar population that is

detected over most of the anticenter direction (Slater et al. 2014) and that is probably the result of the warping and flaring of the outer disk. Our algorithm decomposes the GASS into numerous stream-like filaments. Even though it forms a 25°-long arc within less than 30° of the Galactic plane in the anticenter region, the Gjöll stream stands out as having very different kinematic properties to the GASS population, and at a heliocentric distance of $3.4 \pm 0.1$, it is substantially closer.

In the first submission of the present work, we identified 57 stars as probable Gjöll stream members, and the resulting best-fit orbit is shown in Figure 15. One of those stars is present in the SDSS spectroscopic survey (blue triangle) and possesses a radial velocity that agrees well with the orbit prediction based on the earlier *Gaia* DR2 sample.

We were awarded observing time at the Canada-France-Hawaii Telescope (CFHT) in semester 2018B to measure the radial velocities of our stream candidates, and we decided to follow up bright members of the Gjöll structure, which is particularly well located on the sky for observations in October and November. We selected five stars from the original sample for spectroscopic follow-up with the very high resolution ESPaDOnS instrument at the CFHT observatory. These stars are listed in Table 3, where we also provide the date and CFHT odometer identification of these spectroscopic observations.

The data were reduced using the Libre-ESpRIT pipeline (Donati et al. 1997), and we measured the velocities of the stars by cross-correlation to our observation of the radial velocity template HD 182572 using the *fxcor* command in PyRAF.

With the more conservative parameter choices requested by the reviewers of the present manuscript and of the accompanying study of the Fimbulthul stream, the number of candidate members of Gjöll dropped to 28, and the previously cross-matched SDSS star is no longer a member of the new sample, and neither are three of the five stars we observed with ESPaDOnS. We tried to carefully track down the reason why these stars were lost to the final sample, but there are many causes: some have simply a lower stream likelihood than the adopted $8\sigma$ threshold, but others were removed by the updated filtering. Yet all of these stars follow the predicted trend derived from the earlier sample, as can be seen in the fourth panel of Figure 15. The confirmation of our radial velocity profile prediction proves that the algorithm works as intended and is capable of finding very low surface brightness streams, even close to the Galactic disk, where the density of contaminating sources is high.

We have not as yet used the ESPaDOnS spectra to measure the metallicity of the member stars, and we expect this to be challenging owing to the very low signal-to-noise ratio of the spectra. However, we note that the SDSS star in the earlier sample has $[Fe/H] = -1.78 \pm 0.05$, consistent with a metal-poor population.

This structure is strongly retrograde ($L_z = 2721 \pm 159$ km s$^{-1}$ kpc), which rules out a connection with the GASS. The orbit has pericenter $8.0 \pm 0.2$ kpc and apocenter $31.9 \pm 4.4$ kpc. As can be seen in Figure 18 (pink line), this stream stays close to the plane of the disk.

## 5.8. Leiptr

The 48°-long Leiptr stream (Figure 16) also lies toward the Galactic anticenter but is seen at higher heliocentric distances (>5.5 kpc). The largest sample of candidate members is found using the $[Fe/H] = -1.6$ template. One of the stars we identify
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Orbits of the eight streams, as derived by the MCMC orbit-fitting algorithm. The best-fit orbit solutions presented in Section 5 are shown in the *x*–*y* plane (left) and the *R*–*z* plane (right) and are integrated 1 Gyr into the past and 1 Gyr into the future. The position of the Sun in this coordinate system is (x, y, z) = (−8.122,0,0.017) kpc.</figcaption></figure>

was also measured by the Gaia RVS, and again, the best-fit orbit
passes close to the radial velocity of this star. The orbit has a
pericenter distance of 12.8 ± 0.6 kpc and an apocenter of
85 ± 25 kpc. The strongly retrograde nature of this structure is
reflected in its very high z-component of angular momentum:
$L_z = 4689 \pm 314 \text{ km s}^{-1} \text{ kpc.}$

6. Discussion and Conclusions

In this contribution, we have provided a detailed exposition
of the workings of the updated STREAMFINDER algorithm,
which was used previously to detect five streams in Paper II
(Gaia 1–5) and the Phlegethon stream (Paper III). Part of the
necessary validation process consisted of building a Gaia-like
mock catalog that is devoid of stream structures, based on the
GUMS simulation. This mock was processed in an identical
way to the real Gaia DR2 data. We thereby showed that we can
limit the occurrence of false positives by setting a detection
threshold based on our stream likelihood statistic.

In our stream search we used template stellar population
models of ancient metal-poor stars appropriate for tidal debris
of disrupted globular clusters. Due to computational expense, a
single age (12.5 Gyr) and only three choices in metallicity
([Fe/H] = −1.4, −1.6 and −2.0) were used. This rather coarse
grid in metallicity will be refined in future work, and we will
also explore different age choices.

We detected a large number of candidate stream structures in
the Gaia DR2 catalog between distances of 1 and 10 kpc. Many
of these have only a small number of members (~10) that pass
the selected 8σ detection threshold, but they nevertheless
appear to be plausible stream candidates. Examples of such
cases include the features seen near (ℓ ~ −60°, b ~ 45°) in
Figure 7. We also find that the most metal-rich map (Figure 8)
is highly structured around the bulge region and close to the

<figure><img src="image_21.png" /><figcaption>**Figure 19.** Properties of the orbits of the detected streams in *L*<sub>*z*</sub> vs. pericenter distance. We have included the Phlegethon stream from Paper III, shown in gray. The small circles represent the orbits of the Galactic globular clusters analyzed by Gaia Collaboration et al. (2018b).</figcaption></figure>

disk, showing many dozens of stream-like features. Compar-
ison to the false-positive map derived from our processing of
the GUMS simulation (Figure 23) suggests that some, but
probably not most, of these features are artifacts. Follow-up
spectroscopy is now needed to establish the reality of those
candidates.
<table><caption>**Table 1**

Stream Candidate Stars with Measured Radial Velocity in the SDSS, LAMOST, and Gaia RVS Catalogs</caption><thead><tr><th>R.A. J2000<br>(deg)</th><th>Decl. J2000<br>(deg)</th><th>G<sub>0</sub><br>(mag)</th><th>(G<sub>BP</sub> - G<sub>RP</sub>)<sub>0</sub><br>(mag)</th><th>v<sub>helio</sub><br>(km s<sup>-1</sup>)</th><th>[Fe/H]</th><th>Source</th><th>Stream</th><th>Probable<br>Member?</th></tr></thead><tbody><tr><td>166.127838</td><td>4.167013</td><td>17.722</td><td>0.746</td><td>-98.05 ± 4.17</td><td>-1.78</td><td>S</td><td>Slidr</td><td>Y</td></tr><tr><td>166.144402</td><td>4.184200</td><td>17.093</td><td>0.615</td><td>-93.08 ± 3.04</td><td>-1.88</td><td>S</td><td>Slidr</td><td>Y</td></tr><tr><td>166.883060</td><td>4.226261</td><td>17.413</td><td>0.758</td><td>11.81 ± 2.45</td><td>-1.46</td><td>S</td><td>Slidr</td><td>N</td></tr><tr><td>171.380846</td><td>4.785477</td><td>16.401</td><td>0.677</td><td>246.79 ± 3.27</td><td>-1.60</td><td>S</td><td>Slidr</td><td>N</td></tr><tr><td>160.048762</td><td>6.708863</td><td>15.948</td><td>0.652</td><td>-74.18 ± 2.93</td><td>-1.83</td><td>S</td><td>Slidr</td><td>Y</td></tr><tr><td>164.598178</td><td>7.083380</td><td>16.256</td><td>0.691</td><td>76.21 ± 2.91</td><td>-1.57</td><td>S</td><td>Slidr</td><td>N</td></tr><tr><td>171.452492</td><td>8.174430</td><td>16.485</td><td>0.689</td><td>272.06 ± 2.58</td><td>-1.13</td><td>S</td><td>Slidr</td><td>N</td></tr><tr><td>165.533067</td><td>9.028086</td><td>16.621</td><td>0.634</td><td>-94.34 ± 4.21</td><td>-1.86</td><td>S</td><td>Slidr</td><td>Y</td></tr><tr><td>157.866653</td><td>9.050309</td><td>17.496</td><td>0.760</td><td>-58.27 ± 3.85</td><td>-1.72</td><td>S</td><td>Slidr</td><td>Y</td></tr><tr><td>153.073899</td><td>13.038558</td><td>17.095</td><td>0.711</td><td>-41.80 ± 3.20</td><td>-1.71</td><td>S</td><td>Slidr</td><td>Y</td></tr><tr><td>176.467662</td><td>2.183338</td><td>12.846</td><td>1.083</td><td>101.12 ± 7.56</td><td>-1.20</td><td>L</td><td>Slidr</td><td>N</td></tr><tr><td>158.456446</td><td>5.152103</td><td>14.861</td><td>0.892</td><td>-64.11 ± 10.05</td><td>-1.71</td><td>L</td><td>Slidr</td><td>Y</td></tr><tr><td>171.674974</td><td>5.826564</td><td>16.783</td><td>0.640</td><td>213.89 ± 17.43</td><td>-1.68</td><td>L</td><td>Slidr</td><td>N</td></tr><tr><td>169.298167</td><td>6.306373</td><td>13.840</td><td>0.985</td><td>-116.06 ± 9.26</td><td>-1.69</td><td>L</td><td>Slidr</td><td>Y</td></tr><tr><td>160.048762</td><td>6.708863</td><td>15.948</td><td>0.652</td><td>-69.98 ± 17.13</td><td>-2.13</td><td>L</td><td>Slidr</td><td>Y</td></tr><tr><td>164.515970</td><td>6.835322</td><td>12.336</td><td>1.074</td><td>-99.75 ± 10.14</td><td>-1.77</td><td>L</td><td>Slidr</td><td>Y</td></tr><tr><td>169.580920</td><td>7.111003</td><td>17.198</td><td>0.673</td><td>-140.68 ± 15.56</td><td>-1.74</td><td>L</td><td>Slidr</td><td>Y</td></tr><tr><td>171.452492</td><td>8.174430</td><td>16.485</td><td>0.689</td><td>264.07 ± 11.84</td><td>-1.37</td><td>L</td><td>Slidr</td><td>N</td></tr><tr><td>165.351172</td><td>8.748006</td><td>17.658</td><td>0.730</td><td>233.61 ± 18.56</td><td>-1.38</td><td>L</td><td>Slidr</td><td>N</td></tr><tr><td>166.425204</td><td>9.396366</td><td>14.529</td><td>0.927</td><td>18.71 ± 4.90</td><td>-0.25</td><td>L</td><td>Slidr</td><td>N</td></tr><tr><td>165.450935</td><td>10.228027</td><td>15.254</td><td>0.857</td><td>-108.02 ± 12.99</td><td>-1.76</td><td>L</td><td>Slidr</td><td>Y</td></tr><tr><td>153.969672</td><td>13.511505</td><td>15.326</td><td>0.742</td><td>-55.47 ± 17.18</td><td>-1.95</td><td>L</td><td>Slidr</td><td>Y</td></tr><tr><td>152.074217</td><td>14.758481</td><td>15.889</td><td>0.621</td><td>-54.09 ± 16.97</td><td>-2.00</td><td>L</td><td>Slidr</td><td>Y</td></tr><tr><td>149.086722</td><td>16.938353</td><td>17.467</td><td>0.803</td><td>-43.31 ± 12.45</td><td>-1.48</td><td>L</td><td>Slidr</td><td>Y</td></tr><tr><td>176.467662</td><td>2.183338</td><td>12.833</td><td>1.076</td><td>108.73 ± 1.90</td><td></td><td>G</td><td>Slidr</td><td>N</td></tr><tr><td>164.515970</td><td>6.835322</td><td>12.321</td><td>1.067</td><td>-93.26 ± 2.16</td><td></td><td>G</td><td>Slidr</td><td>Y</td></tr><tr><td>180.587088</td><td>-0.344169</td><td>17.257</td><td>0.630</td><td>-208.56 ± 5.12</td><td>-2.38</td><td>S</td><td>Slygr</td><td>Y</td></tr><tr><td>182.320384</td><td>-0.241553</td><td>18.375</td><td>0.740</td><td>-184.84 ± 5.81</td><td>-2.58</td><td>S</td><td>Slygr</td><td>Y</td></tr><tr><td>182.105642</td><td>0.411199</td><td>17.473</td><td>0.608</td><td>-205.63 ± 3.28</td><td>-3.10</td><td>S</td><td>Slygr</td><td>Y</td></tr><tr><td>169.596593</td><td>-12.042540</td><td>15.359</td><td>0.979</td><td>327.50 ± 1.91</td><td>-1.88</td><td>S</td><td>Ylgr</td><td>Y</td></tr><tr><td>168.852187</td><td>-11.994377</td><td>17.616</td><td>0.818</td><td>322.58 ± 5.34</td><td>-1.92</td><td>S</td><td>Ylgr</td><td>Y</td></tr><tr><td>169.349924</td><td>-11.958677</td><td>16.684</td><td>0.733</td><td>-53.65 ± 3.05</td><td>-1.82</td><td>S</td><td>Ylgr</td><td>N</td></tr><tr><td>200.291401</td><td>-26.516432</td><td>12.131</td><td>1.164</td><td>207.99 ± 1.17</td><td></td><td>G</td><td>Fimbulthul</td><td>Y</td></tr><tr><td>202.470324</td><td>-24.439468</td><td>12.716</td><td>1.058</td><td>191.37 ± 2.62</td><td></td><td>G</td><td>Fimbulthul</td><td>Y</td></tr><tr><td>239.801474</td><td>26.730392</td><td>17.458</td><td>0.885</td><td>-31.06 ± 62.67</td><td>-1.08</td><td>S</td><td>Svöl</td><td>N</td></tr><tr><td>200.725566</td><td>23.984631</td><td>19.061</td><td>0.754</td><td>-70.33 ± 12.53</td><td>-2.09</td><td>S</td><td>Fjörm</td><td>Y</td></tr><tr><td>199.536720</td><td>24.934613</td><td>18.282</td><td>0.735</td><td>-104.11 ± 5.18</td><td>-2.07</td><td>S</td><td>Fjörm</td><td>Y</td></tr><tr><td>201.073852</td><td>32.632881</td><td>18.800</td><td>0.722</td><td>-94.45 ± 7.99</td><td>-1.67</td><td>S</td><td>Fjörm</td><td>Y</td></tr><tr><td>206.874238</td><td>40.968247</td><td>18.674</td><td>0.667</td><td>-68.97 ± 7.67</td><td>-2.07</td><td>S</td><td>Fjörm</td><td>Y</td></tr><tr><td>209.218215</td><td>41.462402</td><td>18.503</td><td>0.676</td><td>-69.67 ± 6.61</td><td>-2.61</td><td>S</td><td>Fjörm</td><td>Y</td></tr><tr><td>208.598025</td><td>42.246064</td><td>18.555</td><td>0.696</td><td>-84.00 ± 6.50</td><td>-2.44</td><td>S</td><td>Fjörm</td><td>Y</td></tr><tr><td>207.328617</td><td>42.609333</td><td>18.350</td><td>0.682</td><td>-78.93 ± 7.64</td><td>-2.45</td><td>S</td><td>Fjörm</td><td>Y</td></tr><tr><td>195.329241</td><td>4.878705</td><td>16.970</td><td>0.839</td><td>-9.97 ± 8.76</td><td>-0.70</td><td>L</td><td>Fjörm</td><td>N</td></tr><tr><td>198.416415</td><td>13.004091</td><td>15.508</td><td>0.913</td><td>-76.30 ± 14.25</td><td>-2.16</td><td>L</td><td>Fjörm</td><td>Y</td></tr><tr><td>72.0846775</td><td>-5.176483</td><td>18.721</td><td>0.859</td><td>-15.18 ± 4.38</td><td>-1.78</td><td>S</td><td>Gjöll</td><td>Y<sup>a</sup></td></tr><tr><td>89.1120223</td><td>-28.189130</td><td>12.241</td><td>1.259</td><td>190.98 ± 0.58</td><td></td><td>G</td><td>Leiptr</td><td>Y</td></tr></tbody></table>

**Notes.** The Source column lists the provenance of the radial velocity and metallicity information: S = SDSS; L = LAMOST; G = Gaia RVS. The positions correspond to the sky coordinates in Gaia DR2 of the sources identified as stream members according to the STREAMFINDER. Cross-matches with SDSS and LAMOST were performed with a 1" search radius.

<sup>a</sup> Although a likely stream member, changes to the algorithm made during the revision of this work caused this star to no longer appear in the 8σ stream candidate sample.
<table><caption>**Table 2**
Orbital Properties of the Eight Star Streams Presented Here</caption><thead><tr><th>Name</th><th>d<br>(kpc)</th><th>r<sub>Peri</sub><br>(kpc)</th><th>r<sub>Apo</sub><br>(kpc)</th><th>L<sub>z</sub><br>(km s<sup>-1</sup> kpc)</th></tr></thead><tbody><tr><td>Slidr</td><td>3.65 ± 0.09</td><td>4.45 ± 0.08</td><td>42.4 ± 3.2</td><td>-1458 ± 9</td></tr><tr><td>Sylgr</td><td>4.10 ± 0.02</td><td>2.49 ± 0.02</td><td>19.4 ± 0.5</td><td>-671 ± 7</td></tr><tr><td>Ylgr</td><td>9.53 ± 0.26</td><td>10.30 ± 0.04</td><td>30.6 ± 0.9</td><td>3087 ± 36</td></tr><tr><td>Fimbulthul</td><td>4.22 ± 0.01</td><td>1.50 ± 0.01</td><td>6.53 ± 0.01</td><td>490 ± 1</td></tr><tr><td>Svöl</td><td>7.76 ± 0.21</td><td>2.53 ± 0.13</td><td>19.1 ± 1.1</td><td>-793 ± 30</td></tr><tr><td>Fjörm</td><td>4.90 ± 0.07</td><td>6.93 ± 0.02</td><td>31.4 ± 0.6</td><td>-2373 ± 7</td></tr><tr><td>Gjöll</td><td>3.38 ± 0.10</td><td>7.96 ± 0.22</td><td>31.9 ± 4.4</td><td>2721 ± 159</td></tr><tr><td>Leiptr</td><td>7.92 ± 0.42</td><td>12.8 ± 0.6</td><td>85.2 ± 26</td><td>4689 ± 314</td></tr></tbody></table>

Note. The tabulated *d* is the heliocentric distance at the anchor point of the fit.

We selected eight high-significance streams (Slidr, Sylgr, Ylgr, Fimbulthul, Svöl, Fjörm, Gjöll, and Leiptr) out of our candidates, all of which contained (by chance) a small number of stars with extant radial velocity measurements. In five cases, our orbital models constructed from the sky position, proper-motion, and parallax information made predictions for the radial velocity trends of the streams that were close to the actual radial velocity measurements. In the two cases (Ylgr and Fimbulthul) where the default model failed, including the measured radial velocities in the fit produced a reasonable representation of the data. We cannot yet tell whether the radial velocity prediction for Svöl is correct or not, as the single spectroscopically observed star in that stream is very likely a contaminant.

We were also able to secure follow-up spectroscopy for the Gjöll stream using the ESPaDOnS high-resolution spectrograph at the CFHT. Of the five Gjöll stream star candidates that we observed based on a sample constructed before the first submission of the present work, all appear to be genuine stream members. This demonstrates the predictive power of the STREAMFINDER algorithm. However, with the more stringent criteria imposed during the revision of this work, only two of the Gjöll stars observed with ESPaDOnS are retained in our final $8\sigma$ sample. This suggests that the present algorithm parameter settings are conservative, and that we will be able to find many more good stream candidates at lower detection thresholds.

We now switch to a brief discussion of the properties of the stream population. It is natural to include Phlegethon and the retrograde GD-1 stream (Koposov et al. 2010) in the current sample of streams, as they are also detected by our algorithm at $>8\sigma$ in the same 1–10 kpc distance range as the eight high-confidence streams reported here (Jhelum is more distant at ~13 kpc; Shipp et al. 2018). It is perhaps initially surprising to notice that 6 out of these 10 systems have retrograde orbits. In contrast, most globular clusters, which after all are the progenitors of such streams, have prograde orbits, as can be seen in Figure 19.

We suspect that this overabundance of retrograde streams is due to a bias of our kinematic method, because prograde streams have lower contrast over the contaminating field star population. If this diagnosis is correct, as the astrometric accuracy of the Gaia data improves over the coming years, we should be able to pick out the prograde streams more and more cleanly. This argument would suggest that there should be many tens of streams with prograde orbits waiting to be detected out to the same distance as the retrograde streams found here (i.e., ~10 kpc). Of course, it is also possible that dynamical effects such as interactions with spiral arms and the bar have preferentially perturbed the prograde streams, phase-mixing them more rapidly and hence rendering them undetectable with our method.

Most of the streams that we have identified do not have any obvious progenitor in the sense that there is no surviving globular cluster with the same pericenter, apocenter, and $L_z$ (see Figure 19). The notable exception is Fimbulthul, which, as we will discuss at greater length in an accompanying contribution, is the trailing arm of the tidal stream of ω Centauri. Figure 19 shows that the Sylgr and Svöl streams have similar pericenter and $L_z$ to the globular clusters M10 and M12; however, the apocenters of these clusters are $5.17 \pm 0.07$ kpc and $4.88 \pm 0.04$ kpc, respectively, and so are inconsistent with being trivially related to Sylgr and Svöl, whose apocenters are $19.4 \pm 0.5$ kpc and $19.1 \pm 1.1$, respectively. What seems more likely is that Sylgr and Svöl are fragments of a now completely disrupted progenitor, and it is interesting to note in this context that the member stars of Sylgr are very metal-poor and include an SDSS star with metallicity $[Fe/H] = -3.10 \pm 0.08$.

Our analysis is beginning to reveal the accretion events that contributed to the current state of our Milky Way but that are now largely dissolved and as a consequence possess extremely low surface brightness. These streams are only detectable thanks to the exquisite astrometric information provided by Gaia DR2, which allows us to pick them out as coherently moving groups over vast swathes of sky. It is interesting to note that some metal-poor globular clusters may have formed prior to reionization (see, e.g., Brodie & Strader 2006, and references therein), and so mapping the properties of their stream remnants may allow us a new handle to probe the physics of that era.

Follow-up studies to measure the radial velocities of the member stars are now essential to be able to employ these structures to measure the Galactic potential and explore their implications for the dark matter. Fjörm, the longest of these features, spans an 87°-long arc in the northern sky and on its own should be a powerful probe of the Galactic potential. With the continued improvements in the accuracy and depth of Gaia data expected over the coming years, we foresee being able to greatly extend the detection horizon for faint stellar streams from ancient low-mass progenitors. Thus, the prospects for this field appear very bright.
<table><caption>**Table 3**
Gjöll Stream Stars Observed with CFHT/ESPaDOnS</caption><thead><tr><th>R.A. J2000<br>(deg)</th><th>Decl. J2000<br>(deg)</th><th>G<sub>0</sub><br>(mag)</th><th>(G<sub>BP</sub>–G<sub>RP</sub>)<sub>0</sub><br>(mag)</th><th>v<sub>helio</sub><br>(km s<sup>-1</sup>)</th><th>CFHT<br>odometer</th><th>Observation<br>Date</th><th>Exposure<br>(s)</th><th>New<br>Sample?</th></tr></thead><tbody><tr><td>65.581723</td><td>-0.863916</td><td>14.026</td><td>0.974</td><td>-80.44 ± 0.26</td><td>2329638</td><td>2018 Oct 24</td><td>600</td><td>N</td></tr><tr><td>74.627326</td><td>-6.423376</td><td>13.157</td><td>1.041</td><td>-10.92 ± 0.28</td><td>2329640</td><td>2018 Oct 24</td><td>300</td><td>Y</td></tr><tr><td>69.793188</td><td>-1.536265</td><td>15.753</td><td>0.841</td><td>-33.17 ± 0.65</td><td>2335537</td><td>2018 Nov 20</td><td>1200</td><td>N</td></tr><tr><td>72.848250</td><td>-6.758777</td><td>16.396</td><td>0.645</td><td>0.81 ± 3.18</td><td>2335539</td><td>2018 Nov 20</td><td>1500</td><td>N</td></tr><tr><td>82.104259</td><td>-13.340007</td><td>15.799</td><td>0.796</td><td>74.41 ± 1.51</td><td>2335541</td><td>2018 Nov 20</td><td>1200</td><td>Y</td></tr></tbody></table>

Note. The positions correspond to the sky coordinates in Gaia DR2 of the sources identified as stream members according to the STREAMFINDER. The final column lists whether the star is part of the final conservative sample of 28 stars generated during the revision of the manuscript.

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.

We thank the staff of the Canada-France-Hawaii Telescope for taking the ESPaDOnS data used here, and for their continued support throughout the project. Based on observations obtained at the Canada-France-Hawaii Telescope (CFHT), which is operated by the National Research Council of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawaii.

Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III website is http://www.sdss3.org/.

SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration, including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for

Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.

Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope, LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences.

PyRAF is a product of the Space Telescope Science Institute, which is operated by AURA for NASA.

This work has been published under the framework of the IdEx Unistra and benefits from funding from the state managed by the French National Research Agency as part of the investments for the future program. R.A.I. and N.F.M. gratefully acknowledge support from a “Programme National Cosmologie et Galaxies” grant.

## Appendix
### Maps of False Positives in GUMS

The Appendix comprises Figures 20–23.
<figure><img src="image_22.png" /><figcaption>**Figure 20.** Map of false-positive stream-like detections with low proper motion ($\sqrt{\mu_{\ell}^2 + \mu_b^2} < 2$ mas yr$^{-1}$), derived by applying the STREAMFINDER to the GUMS simulation. The layout of the panels is identical to that of Figure 5.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 21.** Map of false-positive stream-like detections in the GUMS simulation, using an SSP template with [Fe/H] = -2.0. The layout of the panels is identical to that of Figure 6, to which this figure should be compared.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Figure 22. Same as Figure 21, but for [Fe/H] = -1.6 (cf. Figure 7).</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure 23. Same as Figure 21, but for [Fe/H] = -1.4 (cf. Figure 8).</figcaption></figure>
ORCID iDs

Rodrigo A. Ibata https://orcid.org/0000-0002-3292-9709

Nicolas F. Martin https://orcid.org/0000-0002-1349-202X

References

Amorisco, N. C., Gómez, F. A., Vegetti, S., & White, S. D. M. 2016, MNRAS, 463, L17

Bernard, E. J., Ferguson, A. M. N., Schlafly, E. F., et al. 2016, MNRAS, 463, 1759

Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127

Brodie, J. P., & Strader, J. 2006, ARA&A, 44, 193

Carlberg, R. G. 2012, ApJ, 748, 20

Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv:1612.05560

Cui, X.-Q., Zhao, Y.-H., Chu, Y.-Q., et al. 2012, RAA, 12, 1197

Dehnen, W., & Binney, J. 1998, MNRAS, 294, 429

Donati, J. F., Semel, M., Carter, B. D., Rees, D. E., & Collier Cameron, A. 1997, MNRAS, 291, 658

Eilers, A.-C., Hogg, D. W., Rix, H.-W., & Ness, M. 2019, ApJ, 871, 120

Erkal, D., Belokurov, V., Bovy, J., & Sanders, J. L. 2016, MNRAS, 463, 102

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018a, A&A, 616, A1

Gaia Collaboration, Helmi, A., van Leeuwen, F., et al. 2018b, A&A, 616, A12

Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65

Gravity Collaboration, Abuter, R., Amorim, A., et al. 2018, A&A, 615, L15

Grillmair, C. J. 2009, ApJ, 693, 1118

Grillmair, C. J., & Dionatos, O. 2006, ApJL, 643, L17

Harris, W. E. 2010, arXiv:1012.3224

Hattori, K., Erkal, D., & Sanders, J. L. 2016, MNRAS, 460, 497

Ibata, R. A., Lewis, G. F., Irwin, M. J., & Quinn, T. 2002, MNRAS, 332, 915

Ibata, R. A., Malhan, K., Martin, N. F., & Starkenburg, E. 2018, ApJ, 865, 85

Ibata, R. A., McConnachie, A., Cuillardre, J.-C., et al. 2017, ApJ, 848, 129

Johnston, K. V., Spergel, D. N., & Haydn, C. 2002, ApJ, 570, 656

Klypin, A., Kravtsov, A. V., Valenzuela, O., & Prada, F. 1999, ApJ, 522, 82

Koposov, S. E., Rix, H.-W., & Hogg, D. W. 2010, ApJ, 712, 260

Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2

Malhan, K., & Ibata, R. A. 2018, MNRAS, 477, 4063

Malhan, K., Ibata, R. A., & Martin, N. F. 2018, MNRAS, 481, 3442

McConnachie, A. W. 2012, AJ, 144, 4

Moore, B., Ghigna, S., Governato, F., et al. 1999, ApJ, 524, L19

Pearson, S., Price-Whelan, A. M., & Johnston, K. V. 2017, NatAs, 1, 633

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes in FORTRAN (2nd ed.; Cambridge: Cambridge Univ. Press)

Price-Whelan, A. M., & Bonaca, A. 2018, ApJL, 863, L20

Price-Whelan, A. M., Sesar, B., Johnston, K. V., & Rix, H.-W. 2016, ApJ, 824, 104

Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130

Robin, A. C., Luri, X., Reylé, C., et al. 2012, A&A, 543, A100

Sanderson, C., & Curtin, R. 2016, JOSS, 1, 26

Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103

Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525

Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829

Shipp, N., Drlica-Wagner, A., Balbinot, E., et al. 2018, ApJ, 862, 114

Slater, C. T., Bell, E. F., Schlafly, E. F., et al. 2014, ApJ, 791, 9

Varghese, A., Ibata, R., & Lewis, G. F. 2011, MNRAS, 417, 198

Yanny, B., Rockosi, C., Newberg, H. J., et al. 2009, AJ, 137, 4377