The Streams of the Gaping Abyss: A population of entangled stellar streams surrounding the Inner Galaxy

TL;DR

This paper reports the discovery of numerous stellar streams around the inner Galaxy using Gaia data and a new detection algorithm, revealing their properties and origins, primarily as debris from globular clusters.

Contribution

It introduces the STREAMFINDER algorithm for detecting stellar streams and presents eight new high-significance streams with detailed properties and spectroscopic confirmation.

Findings

Discovered eight new stellar streams in Gaia DR2 data.

Streams are mostly metal-poor and likely originate from globular clusters.

The STREAMFINDER algorithm effectively identifies streams with low false positive rates.

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\deg$, named Slidr, Sylgr, Ylgr, Fimbulthul, Sv\"ol, Fj\"orm, Gj\"oll 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.

Tables

Table 1. Table 1 Gjöll stream stars observed with CFHT/ESPaDOnS.

RA J2000	Dec J2000	$G_0$	${\left(G_{\mathrm{BP}}-G_{\mathrm{RP}}\right)}_0$	$v_{helio}$	CFHT	Observation	Exposure	New
$\left[\mathrm{deg}\right]$	$\left[\mathrm{deg}\right]$	$\left[\mathrm{mag}\right]$	$\left[\mathrm{mag}\right]$	$\left[\mathrm{km}{s}^{-1}\right]$	odometer	date	[s]	sample?
65.581723	-0.863916	14.026	0.974	$-75.54\pm 0.32$	2329638	2018-10-24	600	N
74.627326	-6.423376	13.157	1.041	$-6.27\pm 0.36$	2329640	2018-10-24	300	Y
69.793188	-1.536265	15.753	0.841	$-28.28\pm 0.80$	2335537	2018-11-20	1200	N
72.848250	-6.758777	16.396	0.645	$-9.79\pm 1.20$	2335539	2018-11-20	1500	N
82.104259	-13.340007	15.799	0.796	$79.04\pm 2.08$	2335541	2018-11-20	1200	Y

Table 2. Table 2 Orbital properties of the eight star streams presented here.

Name	$d$	$r_{\mathrm{Peri}}$	$r_{\mathrm{Apo}}$	$L_z$
	$\left[\mathrm{kpc}\right]$	$\left[\mathrm{kpc}\right]$	$\left[\mathrm{kpc}\right]$	$\left[\mathrm{km}{s}^{-1}\mathrm{kpc}\right]$
Slidr	$3.65\pm 0.09$	$4.45\pm 0.08$	$42.4\pm 3.2$	$-1458\pm 9$
Sylgr	$4.10\pm 0.02$	$2.49\pm 0.02$	$19.4\pm 0.5$	$-671\pm 7$
Ylgr	$9.53\pm 0.26$	$10.30\pm 0.04$	$30.6\pm 0.9$	$3087\pm 36$
Fimbulthul	$4.22\pm 0.01$	$1.50\pm 0.01$	$6.53\pm 0.01$	$490\pm 1$
Svöl	$7.76\pm 0.21$	$2.53\pm 0.13$	$19.1\pm 1.1$	$-793\pm 30$
Fjörm	$4.90\pm 0.07$	$6.93\pm 0.02$	$31.4\pm 0.6$	$-2373\pm 7$
Gjöll	$3.38\pm 0.10$	$7.96\pm 0.22$	$31.9\pm 4.4$	$2721\pm 159$
Leiptr	$7.92\pm 0.42$	$12.8\pm 0.6$	$85.2\pm 26$	$4689\pm 314$

Table 3. Table 3 Stream candidate stars with measured radial velocity in the SDSS, LAMOST and Gaia-RVS catalogues.

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

Equations

$$\ln \mathcal{L}=\sum _{\mathrm{data}}\ln [\eta {\mathcal{P}}_{\mathrm{stream}}(\theta )+(1-\eta ){\mathcal{P}}_{\mathrm{cont}}],$$

$${\mathcal{P}}_{\mathrm{stream}}(\theta )={\mathcal{P}}_{\mathrm{length}}\times {\mathcal{P}}_{\mathrm{width}}\times {\mathcal{P}}_{\mathrm{LF}}\times {\mathcal{P}}_{\mathrm{color}}\times {\mathcal{P}}_{\mu }$$

$$\displaystyle\begin{split}\mathcal{P}_{\mu}&={{1}\over{2\pi\sigma_{\mu_{\alpha}}\sigma_{\mu_{\delta}}\sqrt{1-\rho^{2}}}}\times\\ &\exp(-{{1}\over{2(1-\rho^{2})}}\Bigg{[}{{\Delta^{2}_{\mu_{\alpha}}}\over{\sigma^{2}_{\mu_{\alpha}}}}+{{\Delta^{2}_{\mu_{\delta}}}\over{\sigma^{2}_{\mu_{\delta}}}}-{{2\rho\Delta_{\mu_{\alpha}}\Delta_{\mu_{\delta}}}\over{\sigma_{\mu_{\alpha}}\sigma_{\mu_{\delta}}}}\Bigg{]})\\ \end{split}$$

$$\rho =\frac{C{\sigma }_{{\mu }_{\alpha }}{\sigma }_{{\mu }_{\delta }}}{\sqrt{({\sigma }_{{\mu }_{\alpha }}^2+w_{\mu }^2)({\sigma }_{{\mu }_{\delta }}^2+w_{\mu }^2)}}.$$

$$\begin{array}{rl}{\mathcal{P}}_{\mathrm{cont}}(\alpha ,\delta ,{\mu }_{\alpha },{\mu }_{\delta },G,G_{\mathrm{BP}}-G_{\mathrm{RP}})=& \\ P_{\mathrm{fine}}(\alpha ,\delta ,G,G_{\mathrm{BP}}-G_{\mathrm{RP}})\times & \\ P_{\mathrm{GMM}}({\mu }_{\alpha },{\mu }_{\delta }|\alpha ,\delta ,G,G_{\mathrm{BP}}-G_{\mathrm{RP}}).& \end{array}$$

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

$${\mathcal{P}}_{\mathrm{fit}}(\theta )={\mathcal{P}}_{\mathrm{width}}\times {\mathcal{P}}_{\mu }\times {\mathcal{P}}_{\varpi }.$$

Full text

††slugcomment: Accepted by The Astrophysical Journal

1. Introduction

The advent of the second data release (DR2) of the Gaia mission (Lindegren et al., 2018; Gaia Collaboration et al., 2018a) has opened up an 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 three-dimensional 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 $\Lambda$CDM cosmology is the existence of thousands of dark matter sub-halos 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 $\Lambda$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 $\Lambda$CDM, as it would provide incontrovertible proof that dark matter is smooth on small (sub-galactic) 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 modelling of the stream population as an ensemble, and include the effects of non-axisymmetric 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 towards 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 sky 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).

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{\rm\,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>0.5{\rm\,kpc}$. An early result of the present effort was the detection of the $75^{\circ}$-long “Phlegethon” stream, which was presented in Ibata et al. (2018, hereafter Paper III).

The Gaia DR2 survey is the principal dataset 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_{G_{\rm BP}}/A_{V}=1.06794$ and $A_{G_{\rm RP}}/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), as well as within a $3^{\circ}$ radius of M31 and $1^{\circ}$ radius of M33. We also remove stars in appropriately-sized regions (typically $\sim 0.5^{\circ}$) 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, as well as 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 $\mathcal{P}_{\rm stream}(\theta)$ and the corresponding stream fraction $\eta$, that maximizes:

[TABLE]

where $\theta$ are the stream fitting parameters, and the probability density function $\mathcal{P}_{\rm cont}$ is a model of the “contamination” (sometimes called a “background”) from non-stream 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 $\mathcal{P}_{\rm stream}(\theta)$, we use a simple leapfrog scheme to integrate an orbit both forwards and backwards 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 $\mu^{*}_{\alpha}(\equiv\mu_{\alpha}\cos(\delta))$, $\mu_{\delta}$, and color-magnitude position $G$, $G_{\rm BP}-G_{\rm 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{\rm\,kpc}$ (Gravity Collaboration et al., 2018a), and of the circular velocity at the Solar neighborhood of $v_{c}(R_{\odot})=229.0\pm 0.2{\rm\,km\,s^{-1}}$ (Eilers et al., 2018). Given that $v_{c}(R_{\odot})+V_{\odot}=255.2\pm 5.1{\rm\,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{\rm\,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 asterix 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_{\rm 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{\rm\,pc}$ (we use $w$ to denote model Gaussian dispersions and $\sigma$ to denote measurement uncertainties), and a velocity dispersion of $w_{v}=3{\rm\,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:

[TABLE]

where $\mathcal{P}_{\rm length}=1/L$ is the (uniform) probability of lying at the measured position along the stream, and $\mathcal{P}_{\rm width}=\mathcal{N}(\Delta s,w_{s})$ is the Gaussian probability of lying at the measured distance $\Delta s$ perpendicular to the stream. $\mathcal{P}_{\rm 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_{\rm color}=\mathcal{N}(\Delta(G_{\rm Bp}-G_{\rm Rp}),\sigma_{G_{\rm Bp}-G_{\rm 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_{\rm Bp}-G_{\rm Rp}}$. Finally,

[TABLE]

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{\tt pmra\_pmdec\_corr}$ (discussed in Lindegren et al. 2018), which is subsumed into the parameter

[TABLE]

This last equation can be derived by considering the convolution of the two-dimensional covariance matrix with an isotropic two-dimensional 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\times 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 useable in Eqn. 1, the contamination model $\mathcal{P}_{\rm cont}$ also needs to be a probability density function of the same six observables $(\alpha,\delta,\mu_{\alpha},\mu_{\delta},G,G_{\rm BP}-G_{\rm RP})$ as $\mathcal{P}_{\rm 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 catalogue by re-drawing each survey star 1000 times from a two-dimensional Gaussian in tangent-point coordinates using a dispersion of $2^{\circ}$. The resulting catalogue was then binned into a 4-dimensional grid in position and color-magnitude. We adopted a spatial binning scheme in zenithal equal area (ZEA) Galactic polar projection (the projection used in the maps in Paper II) with a bin size of $1.4^{\circ}\times 1.4^{\circ}$. The color-magnitude intervals were chosen to be $0.05\,{\rm mag}\times 0.10\,{\rm 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_{\rm fine}(\alpha,\delta,G,G_{\rm BP}-G_{\rm 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{\rm\,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 modelling we found it necessary to employ $4\times 4$ larger spatial pixels of size $5.6^{\circ}\times 5.6^{\circ}$, so that the spatial regions would contain sufficient stars ($>50000$ stars at the Galactic poles) for the GMM algorithm to yield smooth results. We label the probability density function derived in this way $P_{\rm GMM}(\mu_{\alpha},\mu_{\delta},G,G_{\rm BP}-G_{\rm 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_{\rm BP}-G_{\rm RP}$ measurements, to obtain the conditional probability $P_{\rm GMM}(\mu_{\alpha},\mu_{\delta}|\alpha,\delta,G,G_{\rm BP}-G_{\rm RP})$. Clearly, this means that we have to normalize the function so that $\iint P_{\rm GMM}\mathrm{d}\mu_{\alpha}\,\mathrm{d}\mu_{\delta}=1$.

The final contamination model is then:

[TABLE]

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 both reproduce the complex structure in this six-dimensional 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 two-dimensional slices through the $G,G_{\rm BP}-G_{\rm RP},\mu_{\alpha},\mu_{\delta}$ parameter space. In these visual checks we ensured that the model was not under-fitting the data; we found occasionally that with $\sim 50$ GMM components or less, the model would only reproduce the disk dwarf sequence (the feature at G_{\rm BP}-G_{\rm RP}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}1.5 in Figure 3). We also ensured that the model was not over-fitting 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 multi-dimensional parameter space. A particular worry is that artefacts 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 show 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_{\rm los}^{\rm model},\mu_{\alpha}^{\rm model},\mu_{\delta}^{\rm model}]$, and the stream fraction parameter $\eta$. The distance degeneracy is explored by re-running the algorithm with different plausible SSP models. We find the maximum value of ${\ln}\mathcal{L}$ by first sampling over $v_{\rm los}^{\rm model}$, $\mu_{\alpha}^{\rm model}$, and $\mu_{\delta}^{\rm model}$. Here, just as in Paper I, we sample $v_{\rm los}^{\rm model}$ over the escape velocity ($\pm 600{\rm\,km\,s^{-1}}$) uniformly in steps of $20{\rm\,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}\mathcal{L}$, iterating until convergence is found to an absolute tolerance of $10^{-4}$ in ${\ln}\mathcal{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 this upgrade, as it will require us to include another dimension in the contamination model, but one which 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\sigma$. 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\sim 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 over-densities along orbits. In order to assess these concerns, we decided to also apply the algorithm to the Gaia Universe Model Snapshot (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^{\circ}$ radius with a $G$-band magnitude within $0.1{\rm\,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|\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}20^{\circ}), 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^{\circ}$. 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\hbox{{}^{\circ}\hskip-3.0pt.}4\times 1\hbox{{}^{\circ}\hskip-3.0pt.}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 $\mathcal{P}_{\rm cont}$ directly from the simulated data.

4. Results

The algorithm was run over metallicity values of ${\rm[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^{\circ}$. 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{\rm\,Gyr}$ was adopted. Since the allowed half-length of the orbits was set to $10^{\circ}$, this meant that we actually examined data at all latitudes $|b|>10^{\circ}$. The extension to lower latitude will be presented in future contributions, but this is computationally extremely expensive due 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{\tt phot\_bp\_rp\_excess\_factor})$ in the range:

[TABLE]

This removes stars that have suspicious photometry, which amount to $\approx 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^{2}_{\ell}+\mu^{2}_{b}}<2{\rm\,mas/yr}$), and those of higher proper motion. The solutions with low proper motion form a relatively coherent pattern on sky, as we show in Figure 5, which displays the stream-like features with $\sqrt{\mu^{2}_{\ell}+\mu^{2}_{b}}<2{\rm\,mas/yr}$ that are identified by the algorithm at distances between $1$ and $10{\rm\,kpc}$ when assuming a model of metallicity ${\rm[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_{\ell}$ 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 artefact 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 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^{2}_{\ell}+\mu^{2}_{b}}<2{\rm\,mas/yr}$, 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^{2}_{\ell}+\mu^{2}_{b}}>2{\rm\,mas/yr}$ in the distance range between $1$–$10{\rm\,kpc}$ derived using the SSP template with age $12.5{\rm\,Gyr}$ and metallicity ${\rm[Fe/H]=-2.0}$. Figures 7 and 8 show similar information but using a template with metallicity of ${\rm[Fe/H]=-1.6}$ and ${\rm[Fe/H]=-1.4}$, respectively. The layout 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\sigma$ 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{\rm\,kpc}$ artificially limits this stream to \ell\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}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\sigma$ threshold. We have marked and named ten 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 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, and are eight of the eleven 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, 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 towards 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 in order of increasing longitude in Figure 6, and then we will proceed to examine 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^{\circ}$-long prograde structure in the northern hemisphere maps. Although it is detected with all of the chosen metallicity templates, it appears strongest in the ${\rm[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 dots with their corresponding uncertainties. A total of 156 stars are detected at $>8\sigma$ 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 re-fitting 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 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. The algorithm maximises the likelihood of the orbit given the data, adopting the following probability density function for the stream:

[TABLE]

The terms $\mathcal{P}_{\rm width}$ and $\mathcal{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 $\mathcal{P}_{\varpi}=\mathcal{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^{\circ}$, $\delta=6^{\circ}$) is $d=3.65\pm 0.09{\rm\,kpc}$.

In addition, the best-fit orbit gives a prediction of the heliocentric line of sight velocity of the structure (shown on the fourth panel). By chance, our Slidr sample contains ten stars previously observed in spectroscopy by the SDSS DR10 (Yanny et al. 2009, large blue points), 14 stars observed by LAMOST DR4 (Cui et al. 2012, large magenta points), as well as two stars observed by the Gaia Radial Velocity Spectrometer (RVS; large orange points). The candidate stream members with radial velocity measurements are listed in Table 3. 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 ${\rm[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 members stars in the SDSS sample is ${\rm[Fe/H]=-1.80}$, essentially identical to the mean metallicity of ${\rm[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\sigma$ selection threshold indicates that, at the position of each identified star, there is a stream present at $>8\sigma$ confidence that shares the kinematic properties of the star (and is consistent with its photometric properties and parallax). In particular, it does not mean that each stream candidate star identified by the STREAMFINDER is a member with $>8\sigma$ confidence of the detected stream.

The orbit of Slidr is prograde but highly radial, with a pericenter at $4.45\pm 0.08{\rm\,kpc}$ and apocenter at $42.4\pm 3.2{\rm\,kpc}$. Its $z$-component of angular momentum is $L_{z}=-1458\pm 9{\rm\,km\,s^{-1}}{\rm\,kpc}$, according to the MCMC orbital fit (we adopt the usual convention that negative $L_{z}$ means prograde motion). Its orbit is shown in blue in Figure 18, and its position in the $L_{z}$ – 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^{\circ}$-long nearby structure that appears to have highest contrast in the ${\rm[Fe/H]=-1.6}$ map (Figure 7), where it can be easily identified as a stream of 103 stars with 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 (${\rm\langle[Fe/H]\rangle=-2.7}$). It will be interesting to re-run 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{\rm\,km\,s^{-1}}{\rm\,kpc}$), with a pericenter at $2.49\pm 0.02{\rm\,kpc}$ and an apocenter at $19.5\pm 0.5{\rm\,kpc}$. The CMD and orbital properties of the structure are shown with orange points in Figures 17–19.

5.3. Ylgr

Figure 11 shows the properties of the $30^{\circ}$-long Ylgr stream. It is most clearly visible in the ${\rm[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 three of those stars are metal-poor (${\rm\langle[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 into the likelihood, and assuming that the star with $v_{helio}=-53.6{\rm\,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{\rm\,km\,s^{-1}}{\rm\,kpc}$), with a pericenter at $10.3\pm 0.04{\rm\,kpc}$, and an apocenter at $30.6\pm 0.9{\rm\,kpc}$.

5.4. Fimbulthul

Fimbulthul is perhaps the most fascinating of the streams we have identified, and will be discussed in greater detail elsewhere (Ibata et al 2019, under review). A total of 309 stars over a $18^{\circ}$ 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 re-fit 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.

The revised orbit model has a pericentric distance of $1.50\pm 0.01{\rm\,kpc}$, an apocenter distance of $6.53\pm 0.01{\rm\,kpc}$ and is again highly retrograde ($L_{z}=490\pm 1{\rm\,km\,s^{-1}}{\rm\,kpc}$). These properties are very close to those of the prototypical massive cluster $\omega$ Centauri, and we show in Ibata et al (2019) that the Fimbulthul stream is indeed the long sought-for tidal tail of that globular cluster.

5.5. Svöl

Svöl is an $26^{\circ}$-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 (${\rm[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\pm 0.13{\rm\,kpc}$, the apocenter is at $19.1\pm 1.1{\rm\,kpc}$, and the stream is prograde, with $L_{z}=-793\pm 30{\rm\,km\,s^{-1}}{\rm\,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, 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^{\circ}$ 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 ${\rm\langle[Fe/H]\rangle=-2.20}$ and sample dispersion of $0.30$ dex. Of the two LAMOST stars, one has ${\rm[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 (${\rm[Fe/H]=-0.70\pm 0.18}$) and is also a velocity outlier.

The best stream fit has a pericenter at $6.93\pm 0.02{\rm\,kpc}$, an apocenter at $31.4\pm 0.6{\rm\,kpc}$, and is prograde with $L_{z}-2373\pm 7{\rm\,km\,s^{-1}}{\rm\,kpc}$.

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^{\circ}$-long arc within less than $30^{\circ}$ 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 1, 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 HD182572 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 on 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 due to the very low signal to noise ratio of the spectra. However, we note that the SDSS star in the earlier sample has ${\rm[Fe/H]=-1.78\pm 0.05}$, consistent with a metal-poor population.

This structure is strongly retrograde ($L_{z}=2721\pm 159{\rm\,km\,s^{-1}}{\rm\,kpc}$), which rules out a connection with the GASS. The orbit has pericenter $8.0\pm 0.2{\rm\,kpc}$ and apocenter $31.9\pm 4.4{\rm\,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^{\circ}$-long Leiptr stream (Figure 16) also lies towards the Galactic anticenter, but is seen at higher Heliocentric distances ($>5.5{\rm\,kpc}$). The largest sample of candidate members is found using the ${\rm[Fe/H]=-1.6}$ template. One of the stars we identify 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\pm 0.6{\rm\,kpc}$, and an apocenter of $85\pm 25{\rm\,kpc}$. The strongly retrograde nature of this structure is reflected in its very high $z$-component of angular momentum: $L_{z}=4689\pm 314{\rm\,km\,s^{-1}}{\rm\,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 5 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 Gaia Universe Model Snapshot (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 populations models of ancient metal-poor stars appropriate for tidal debris of disrupted globular clusters. Due to computational expense, a single age ($12.5{\rm\,Gyr}$), and only three choices in metallicity (${\rm[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{\rm\,kpc}$. Many of these have only a small number of members ($\sim 10$) that pass the selected $8\sigma$ detection threshold, but they nevertheless appear to be plausible stream candidates. Examples of such cases include the features seen near ($\ell\sim-60^{\circ}$, $b\sim 45^{\circ}$) 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 disk, showing many dozens of stream-like features. Comparison to the false positives map derived from our processing of the GUMS simulation (Figure 23) suggests that some, but probably not most, of these features are artefacts. Follow-up spectroscopy is now needed to establish the reality of those candidates.

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 5 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) into the current sample of streams, as they are also detected by our algorithm at $>8\sigma$ in the same $1$–$10{\rm\,kpc}$ distance range as the eight high-confidence streams reported here (Jhelum is more distant at $\sim 13{\rm\,kpc}$, Shipp et al. 2018). It is perhaps initially surprising to notice that six out of these ten 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 improve 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. $\sim 10{\rm\,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 $\omega$ 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{\rm\,kpc}$ and $4.88\pm 0.04{\rm\,kpc}$, respectively, and so are inconsistent with being trivially related to Sylgr and Svöl whose apocenters are $19.4\pm 0.5{\rm\,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 ${\rm[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 the their implications for the dark matter. Fjörm, the longest of these features, spans an $87^{\circ}$-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.

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 Scientique 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 web site 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 a funding from the state managed by the French National Research Agency as part of the investments for the future program. RAI and NFM gratefully acknowledge support from a “Programme National Cosmologie et Galaxies” grant.

Appendix – Maps of false positives in GUMS