# Charting the Galactic Acceleration Field. II. A Global Mass Model of the Milky Way from the STREAMFINDER Atlas of Stellar Streams Detected in Gaia DR3

Rodrigo Ibata, Khyati Malhan, Wassim Tenachi, Anke Ardern-Arentsen, Michele Bellazzini, Paolo Bianchini, Piercarlo Bonifacio, Elisabetta Caffau, Foivos Diakogiannis, Raphael Errani, et al.

► To cite this version:

Rodrigo Ibata, Khyati Malhan, Wassim Tenachi, Anke Ardern-Arentsen, Michele Bellazzini, et al. Charting the Galactic Acceleration Field. II. A Global Mass Model of the Milky Way from the STREAMFINDER Atlas of Stellar Streams Detected in Gaia DR3. The Astrophysical Journal, 2024, 967, 10.3847/1538-4357/ad382d. insu-04590292

HAL Id: insu-04590292

https://insu.hal.science/insu-04590292

Submitted on 28 May 2024

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
# Charting the Galactic Acceleration Field. II. A Global Mass Model of the Milky Way from the STREAMFINDER Atlas of Stellar Streams Detected in Gaia DR3

Rodrigo Ibata, Khyati Malhan, Wassim Tenachi, Anke Ardern-Arentsen, Michele Bellazzini, Paolo Bianchini, Piercarlo Bonifacio, Elisabetta Caffau, Foivos Diakogiannis, Raphael Errani, et al.

► To cite this version:

Rodrigo Ibata, Khyati Malhan, Wassim Tenachi, Anke Ardern-Arentsen, Michele Bellazzini, et al. Charting the Galactic Acceleration Field. II. A Global Mass Model of the Milky Way from the STREAMFINDER Atlas of Stellar Streams Detected in Gaia DR3. The Astrophysical Journal, 2024, 967, 10.3847/1538-4357/ad382d. insu-04590292

HAL Id: insu-04590292

https://insu.hal.science/insu-04590292v1

Submitted on 28 May 2024

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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>

# Charting the Galactic Acceleration Field. II. A Global Mass Model of the Milky Way from the STREAMFINDER Atlas of Stellar Streams Detected in Gaia DR3

Rodrigo Ibata¹, Khyati Malhan²,³, Wassim Tenachi¹, Anke Ardern-Arentsen⁴, Michele Bellazzini⁵, Paolo Bianchini¹, Piercarlo Bonifacio⁶, Elisabetta Caffau⁶, Foivos Diakogiannis⁷, Raphael Errani¹, Benoit Famaey¹, Salvatore Ferrone⁶, Nicolas F. Martin¹, Paola di Matteo⁶, Giacomo Monari¹, Florent Renaud¹,⁸, Else Starkenburg⁹, Guillaume Thomas¹⁰,¹¹, Akshara Viswanathan⁹, and Zhen Yuan¹

¹ Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France; rodrigo.ibata@astro.unistra.fr

²Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117, Heidelberg, Germany

³Myrspoven AB, Västgötagatan 1, 11827 Stockholm, Sweden

⁴Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

⁵INAF—Osservatorio di Astrofisica e Scienza dello Spazio, via Gobetti 93/3, I-40129 Bologna, Italy

⁶GEPI, Observatoire de Paris, Université PSL, CNRS, 5 Place Jules Janssen, 92190 Meudon, France

⁷Data61, CSIRO, Floreat WA, Australia

⁸University of Strasbourg Institute for Advanced Study, 5 allée du Général Rouvillois, F-67083 Strasbourg, France

⁹Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD Groningen, The Netherlands

¹⁰Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain

¹¹Universidad de La Laguna, Dpto. Astrofísica, E-38206 La Laguna, Tenerife, Spain

Received 2023 November 28; revised 2024 March 21; accepted 2024 March 21; published 2024 May 21

## Abstract

We present an atlas and follow-up spectroscopic observations of 87 thin stream-like structures detected with the STREAMFINDER algorithm in Gaia DR3, of which 28 are new discoveries. Here, we focus on using these streams to refine mass models of the Galaxy. Fits with a double-power-law halo with the outer power-law slope set to $-\beta_h = -3$ yield an inner power-law slope of $-\gamma_h = -(0.97_{-0.21}^{+0.17})$, a scale radius of $r_{0,h} = 14.7_{-1.0}^{+4.7}$ kpc, a halo density flattening $q_{m,h} = 0.75 \pm 0.03$, and a local dark matter density of $\rho_{h,\odot} = 0.0114 \pm 0.0007 M_{\odot} pc^{-3}$. Freeing $\beta$ yields $\beta = 2.53_{-0.16}^{+0.42}$, but this value is heavily influenced by our chosen virial mass limit. The stellar disks are found to have a combined mass of $4.20_{-0.53}^{+0.44} \times 10^{10} M_{\odot}$, with the thick disk contributing $12.4\% \pm 0.7\%$ to the local stellar surface density. The scale lengths of the thin and thick disks are $2.17_{-0.08}^{+0.18}$ and $1.62_{-0.13}^{+0.72}$ kpc, respectively, while their scale heights are $0.347_{-0.010}^{+0.007}$ and $0.86_{-0.02}^{+0.03}$ kpc, respectively. The virial mass of the favored model is $M_{200} = 1.09_{-0.14}^{+0.19} \times 10^{12} M_{\odot}$, while the mass inside of 50 kpc is $M_{R<50} = 0.46 \pm 0.03 \times 10^{12} M_{\odot}$. We introduce the Large Magellanic Cloud (LMC) into the derived potential models, and fit the Orphan stream therein, finding a mass for the LMC that is consistent with recent estimates. Some salient highlights include the nearby trailing arm of ω Cen, and a nearby very metal-poor stream that was once a satellite of the Sagittarius dwarf galaxy. Finally, we unambiguously detect a hot component around the GD-1 stream, consistent with it having been tidally preprocessed within its own dark matter subhalo.

Unified Astronomy Thesaurus concepts: Stellar streams (2166); Galaxy structure (622); Dark matter (353)

Supporting material: machine-readable table

## 1. Introduction

The formation of galaxies involves a complex interplay between the accretion of gas and the merging of previously bound fragments. When these fragments are of low mass compared to their host or their dominant companion, tidal forces can have time to act progressively to slowly unbind material from the lower-mass object. If the disruption is gradual, so that the injected energy is only just enough to unbind the most loosely bound stars, the ejected objects will find themselves on similar orbits to their progenitor, and diffuse along a tidal stream pattern that approximates that orbit (Combes et al. 1999). The stars that by chance acquire higher kinetic energy leave approximately through the L1 Lagrange point and race ahead of their progenitor, forming the leading tidal arm. Conversely, those stars that are perturbed into lower-

energy orbits emerge roughly through the L2 point and form the trailing arm (see, e.g., Fukushige & Heggie 2000).

Tidal disruption can take place over hundreds of millions to billions of years, with the dwindling progenitor slowly losing stars that drift away along their slightly different orbits. The requirement that the stream stars originated in their progenitor provides a physical constraint that allows us to attempt to turn back the clock and find the smooth potential and any perturbations that can permit this inversion. This is perhaps the most useful property of stellar streams, giving us a means to peer into the dynamical past of a galaxy without making assumptions about its equilibrium, and so uncover properties of the global and local dark matter distribution that are otherwise invisible. Nevertheless, it should be borne in mind that tidal disruption typically produces ejecta with a spectrum of energies, meaning that stars move along the tidal stream pattern at different speeds, so that position along the pattern is not perfectly correlated with ejection time. Furthermore, epicyclic patterns are also present in streams, which can complicate their analysis (Küpper et al. 2008).

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
<figure><img src="image_2.png" /></figure>

# Charting the Galactic Acceleration Field. II. A Global Mass Model of the Milky Way from the STREAMFINDER Atlas of Stellar Streams Detected in Gaia DR3

Rodrigo Ibata¹, Khyati Malhan²,³, Wassim Tenachi¹, Anke Ardern-Arentsen⁴, Michele Bellazzini⁵, Paolo Bianchini¹, Piercarlo Bonifacio⁶, Elisabetta Caffau⁶, Foivos Diakogiannis⁷, Raphael Errani¹, Benoit Famaey¹, Salvatore Ferrone⁶, Nicolas F. Martin¹, Paola di Matteo⁶, Giacomo Monari¹, Florent Renaud¹,⁸, Else Starkenburg⁹, Guillaume Thomas¹⁰,¹¹, Akshara Viswanathan⁹, and Zhen Yuan¹

¹ Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France; rodrigo.ibata@astro.unistra.fr

²Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117, Heidelberg, Germany

³Myrspoven AB, Västgötagatan 1, 11827 Stockholm, Sweden

⁴Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

⁵INAF—Osservatorio di Astrofisica e Scienza dello Spazio, via Gobetti 93/3, I-40129 Bologna, Italy

⁶GEPI, Observatoire de Paris, Université PSL, CNRS, 5 Place Jules Janssen, 92190 Meudon, France

⁷Data61, CSIRO, Floreat WA, Australia

⁸University of Strasbourg Institute for Advanced Study, 5 allée du Général Rouvillois, F-67083 Strasbourg, France

⁹Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD Groningen, The Netherlands

¹⁰Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain

¹¹Universidad de La Laguna, Dpto. Astrofísica, E-38206 La Laguna, Tenerife, Spain

Received 2023 November 28; revised 2024 March 21; accepted 2024 March 21; published 2024 May 21

## Abstract

We present an atlas and follow-up spectroscopic observations of 87 thin stream-like structures detected with the STREAMFINDER algorithm in Gaia DR3, of which 28 are new discoveries. Here, we focus on using these streams to refine mass models of the Galaxy. Fits with a double-power-law halo with the outer power-law slope set to $-\beta_h = -3$ yield an inner power-law slope of $-\gamma_h = -(0.97_{-0.21}^{+0.17})$, a scale radius of $r_{0,h} = 14.7_{-1.0}^{+4.7}$ kpc, a halo density flattening $q_{m,h} = 0.75 \pm 0.03$, and a local dark matter density of $\rho_{h,\odot} = 0.0114 \pm 0.0007 M_\odot\text{pc}^{-3}$. Freeing $\beta$ yields $\beta = 2.53_{-0.16}^{+0.42}$, but this value is heavily influenced by our chosen virial mass limit. The stellar disks are found to have a combined mass of $4.20_{-0.53}^{+0.44} \times 10^{10} M_\odot$, with the thick disk contributing $12.4\% \pm 0.7\%$ to the local stellar surface density. The scale lengths of the thin and thick disks are $2.17_{-0.08}^{+0.18}$ and $1.62_{-0.13}^{+0.72}$ kpc, respectively, while their scale heights are $0.347_{-0.010}^{+0.007}$ and $0.86_{-0.02}^{+0.03}$ kpc, respectively. The virial mass of the favored model is $M_{200} = 1.09_{-0.14}^{+0.19} \times 10^{12} M_\odot$, while the mass inside of 50 kpc is $M_{R<50} = 0.46 \pm 0.03 \times 10^{12} M_\odot$. We introduce the Large Magellanic Cloud (LMC) into the derived potential models, and fit the Orphan stream therein, finding a mass for the LMC that is consistent with recent estimates. Some salient highlights include the nearby trailing arm of ω Cen, and a nearby very metal-poor stream that was once a satellite of the Sagittarius dwarf galaxy. Finally, we unambiguously detect a hot component around the GD-1 stream, consistent with it having been tidally preprocessed within its own dark matter subhalo.

Unified Astronomy Thesaurus concepts: Stellar streams (2166); Galaxy structure (622); Dark matter (353)

Supporting material: machine-readable table

## 1. Introduction

The formation of galaxies involves a complex interplay between the accretion of gas and the merging of previously bound fragments. When these fragments are of low mass compared to their host or their dominant companion, tidal forces can have time to act progressively to slowly unbind material from the lower-mass object. If the disruption is gradual, so that the injected energy is only just enough to unbind the most loosely bound stars, the ejected objects will find themselves on similar orbits to their progenitor, and diffuse along a tidal stream pattern that approximates that orbit (Combes et al. 1999). The stars that by chance acquire higher kinetic energy leave approximately through the L1 Lagrange point and race ahead of their progenitor, forming the leading tidal arm. Conversely, those stars that are perturbed into lower-

energy orbits emerge roughly through the L2 point and form the trailing arm (see, e.g., Fukushige & Heggie 2000).

Tidal disruption can take place over hundreds of millions to billions of years, with the dwindling progenitor slowly losing stars that drift away along their slightly different orbits. The requirement that the stream stars originated in their progenitor provides a physical constraint that allows us to attempt to turn back the clock and find the smooth potential and any perturbations that can permit this inversion. This is perhaps the most useful property of stellar streams, giving us a means to peer into the dynamical past of a galaxy without making assumptions about its equilibrium, and so uncover properties of the global and local dark matter distribution that are otherwise invisible. Nevertheless, it should be borne in mind that tidal disruption typically produces ejecta with a spectrum of energies, meaning that stars move along the tidal stream pattern at different speeds, so that position along the pattern is not perfectly correlated with ejection time. Furthermore, epicyclic patterns are also present in streams, which can complicate their analysis (Küpper et al. 2008).

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
This promise of being able to map out the acceleration field of our Galaxy motivated us to develop a dedicated stream-detection algorithm, STREAMFINDER, with the intention to deploy it on the Gaia mission catalogs (Gaia Collaboration et al. 2016). STREAMFINDER is effectively a friend-finding algorithm, with a distance in the parameter space of observables defined so as to make objects on similar orbits and with similar stellar populations appear close together. The procedure is presented in detail in Malhan & Ibata (2018), Malhan et al. (2018, hereafter Paper I), and Ibata et al. (2021b, hereafter Paper II). Paper I applied the algorithm to the Gaia DR2 catalog, based on 22 months of astrometric observations, while Paper II extended the search to the Gaia EDR3 catalog, with 33 months of observations. In Martin et al. (2022), the STREAMFINDER sources were crossmatched with the Pristine survey catalog (Starkenburg et al. 2017; Martin et al. 2023), providing metallicity estimates for the stars and hence better discrimination against contamination, which allowed us to lower the detection threshold and so find further stream candidates. A compendium of Galactic stellar streams was compiled by Mateu (2023), and a bibliography of previous searches is given therein.

Our aim in this contribution is to present the STREAMFINDER results for the full Gaia DR3 data set (Gaia Collaboration et al. 2023), which supersedes EDR3. The new analysis makes use of the radial velocity measurements for 33 million stars provided by Gaia DR3, which we complement with measurements from other large spectroscopic surveys. We will also present measurements from our own dedicated spectroscopic follow-up with the Very Large Telescope (VLT)/UVES and Isaac Newton Telescope/Intermediate Dispersion Spectrograph (INT/IDS) as well as our earlier CFHT/ESPaDOnS measurements, providing the community the full Gaia DR3 stream catalog as uncovered by our algorithm. The main scientific goal of the present contribution is to conduct a detailed analysis of the constraints that these observations impose on the large-scale mass distribution in the Milky Way.

The outline of this paper is as follows. In Section 2, we begin by explaining the small alterations we made to the algorithm for this new Gaia DR3 analysis, while Section 3 briefly reviews our follow-up spectroscopic observations, and Section 4 presents the new sky maps and detections. In Section 5, we begin the modeling of the streams, assuming that the Milky Way is axisymmetric and that its center defines the origin of an inertial frame around which the streams orbit. Section 6 generalizes the model to allow the Local Group galaxies to move under the influence of each other, and the streams will be modeled as dissolving systems. In Section 7, we explain how we use test-particle simulations to correct simple orbit integrations, and thus greatly speed up the exploration of the parameters of the model. The results of the stream-fitting analysis are presented in Section 8. In Section 9, we briefly summarize some salient streams, and finally draw our conclusions in Section 10.

## 2. STREAMFINDER Alterations

Here, we use a search algorithm version very similar to that presented in Paper II. As before, we aim to find a low-contrast stream superimposed on a contaminating population of normal

Milky Way stars, and so we adopt a log-likelihood function of

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

where $\eta$ is the fraction of stars in the stream, $\mathcal{P}_{\text{stream}}(\theta)$ is the probability distribution function (PDF) of the stream model with parameters $\theta$, and $\mathcal{P}_{\text{cont}}$ is the PDF of the contaminants. The stream model is a model of a stellar population along a simple orbit that has been smeared out with a Gaussian

$$ \begin{aligned} \mathcal{P}_{\text{stream}}(\theta) &= \mathcal{P}_{\text{length}} \times \mathcal{P}_{\text{width}} \times \mathcal{P}_{\text{LF}} \times \mathcal{P}_{\text{color}} \times \mathcal{P}_{\mu} \times \mathcal{P}_{\varpi} \\ & \quad \times \mathcal{P}_{\nu} \\ &= \mathcal{P}_{\text{stream}}^{\text{Paper II}}(\theta) \times \mathcal{P}_{\nu}, \end{aligned} \quad (2) $$

where $\mathcal{P}_{\text{length}}$ is a uniform distribution along an $L = 20^{\circ}$ stream fragment, $\mathcal{P}_{\text{width}}$ is a Gaussian width of 50 pc dispersion, $\mathcal{P}_{\text{LF}}$ is the probability of being drawn from the luminosity function of the adopted simple stellar population (SSP) model (we adopt the PARSEC isochrone models, Bressan et al. 2012), $\mathcal{P}_{\text{color}}$ is a probability of the measured $G_{\text{BP}} - G_{\text{RP}}$ color of the star at the observed $G$ magnitude, and $\mathcal{P}_{\mu}$ and $\mathcal{P}_{\varpi}$ are the Gaussian probabilities of the observed astrometry given an assumed 5 km s$^{-1}$ intrinsic stream velocity dispersion, and a 50 pc line-of-sight dispersion. (The proper motion probability term $\mathcal{P}_{\mu}$ is a two-dimensional Gaussian PDF incorporating the proper motion uncertainties in R.A. and decl. and their correlation cross-terms.) Up to here, all is identical to the model described in Paper II. However, since the DR3 catalog provides many line-of-sight velocity measurements, we now also include a simple Gaussian probability $\mathcal{P}_{\nu}$ that the star has a line-of-sight velocity within 5 km s$^{-1}$ of the orbit. Most stars in the DR3 catalog do not have measured velocities, and for those, we set $\mathcal{P}_{\nu} = 1$.

In an earlier contribution we used a much wider (500 pc) spatial template in order to detect the stream of the Sagittarius dwarf galaxy and any similar structures (Ibata et al. 2020a). In contrast, our decision in the present work to adopt a characteristic 50 pc width spread and 5 km s$^{-1}$ velocity spread for the template was made to facilitate the detection of globular cluster streams. In the Harris (2010) catalog of Milky Way globular clusters, 79% have tidal radii <50 pc (their Gaussian width is, of course, much smaller), and 78% have a central velocity dispersion <5 km s$^{-1}$. Streams with much wider spatial or velocity dispersions will have their significance downweighted by our template choice; in future work, we intend to overcome this limitation by using a wider range of templates.

Including velocities in the likelihood function of the stream model means that we also have to include them in the contamination model $\mathcal{P}_{\text{cont}}$. We decided to retain the previously fitted, completely empirical, Gaussian mixture model (GMM) description $\mathcal{P}_{\text{cont}}^{\text{Paper II}}(\alpha, \delta, \mu_{\alpha}, \mu_{\delta}, G, G_{\text{BP}} - G_{\text{RP}}, \varpi)$ but include a new term $\mathcal{P}_{\text{GMM}}$ as follows:

$$ \begin{aligned} \mathcal{P}_{\text{cont}}(\alpha, \delta, \mu_{\alpha}, \mu_{\delta}, G, G_{\text{BP}} - G_{\text{RP}}, \varpi, \nu) \\ &= \mathcal{P}_{\text{cont}}^{\text{Paper II}}(\alpha, \delta, \mu_{\alpha}, \mu_{\delta}, G, G_{\text{BP}} - G_{\text{RP}}, \varpi) \\ & \quad \times \mathcal{P}_{\text{GMM}}(\nu|\alpha, \delta, G, G_{\text{BP}} - G_{\text{RP}}, \varpi). \end{aligned} \quad (3) $$

Note that the final term ignores the correlations of line-of-sight velocity with proper motion, which we felt was a necessary simplifying assumption given the relative scarcity of the velocity information compared to the astrometric information.
This promise of being able to map out the acceleration field of our Galaxy motivated us to develop a dedicated stream-detection algorithm, STREAMFINDER, with the intention to deploy it on the Gaia mission catalogs (Gaia Collaboration et al. 2016). STREAMFINDER is effectively a friend-finding algorithm, with a distance in the parameter space of observables defined so as to make objects on similar orbits and with similar stellar populations appear close together. The procedure is presented in detail in Malhan & Ibata (2018), Malhan et al. (2018, hereafter Paper I), and Ibata et al. (2021b, hereafter Paper II). Paper I applied the algorithm to the Gaia DR2 catalog, based on 22 months of astrometric observations, while Paper II extended the search to the Gaia EDR3 catalog, with 33 months of observations. In Martin et al. (2022), the STREAMFINDER sources were crossmatched with the Pristine survey catalog (Starkenburg et al. 2017; Martin et al. 2023), providing metallicity estimates for the stars and hence better discrimination against contamination, which allowed us to lower the detection threshold and so find further stream candidates. A compendium of Galactic stellar streams was compiled by Mateu (2023), and a bibliography of previous searches is given therein.

Our aim in this contribution is to present the STREAMFINDER results for the full Gaia DR3 data set (Gaia Collaboration et al. 2023), which supersedes EDR3. The new analysis makes use of the radial velocity measurements for 33 million stars provided by Gaia DR3, which we complement with measurements from other large spectroscopic surveys. We will also present measurements from our own dedicated spectroscopic follow-up with the Very Large Telescope (VLT)/UVES and Isaac Newton Telescope/Intermediate Dispersion Spectrograph (INT/IDS) as well as our earlier CFHT/ESPaDOnS measurements, providing the community the full Gaia DR3 stream catalog as uncovered by our algorithm. The main scientific goal of the present contribution is to conduct a detailed analysis of the constraints that these observations impose on the large-scale mass distribution in the Milky Way.

The outline of this paper is as follows. In Section 2, we begin by explaining the small alterations we made to the algorithm for this new Gaia DR3 analysis, while Section 3 briefly reviews our follow-up spectroscopic observations, and Section 4 presents the new sky maps and detections. In Section 5, we begin the modeling of the streams, assuming that the Milky Way is axisymmetric and that its center defines the origin of an inertial frame around which the streams orbit. Section 6 generalizes the model to allow the Local Group galaxies to move under the influence of each other, and the streams will be modeled as dissolving systems. In Section 7, we explain how we use test-particle simulations to correct simple orbit integrations, and thus greatly speed up the exploration of the parameters of the model. The results of the stream-fitting analysis are presented in Section 8. In Section 9, we briefly summarize some salient streams, and finally draw our conclusions in Section 10.

## 2. STREAMFINDER Alterations

Here, we use a search algorithm version very similar to that presented in Paper II. As before, we aim to find a low-contrast stream superimposed on a contaminating population of normal

Milky Way stars, and so we adopt a log-likelihood function of

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

where $\eta$ is the fraction of stars in the stream, $\mathcal{P}_{\text{stream}}(\theta)$ is the probability distribution function (PDF) of the stream model with parameters $\theta$, and $\mathcal{P}_{\text{cont}}$ is the PDF of the contaminants. The stream model is a model of a stellar population along a simple orbit that has been smeared out with a Gaussian

$$ \begin{aligned} \mathcal{P}_{\text{stream}}(\theta) &= \mathcal{P}_{\text{length}} \times \mathcal{P}_{\text{width}} \times \mathcal{P}_{\text{LF}} \times \mathcal{P}_{\text{color}} \times \mathcal{P}_{\mu} \times \mathcal{P}_{\varpi} \\ & \quad \times \mathcal{P}_{\nu} \\ &= \mathcal{P}_{\text{stream}}^{\text{Paper II}}(\theta) \times \mathcal{P}_{\nu}, \end{aligned} \quad (2) $$

where $\mathcal{P}_{\text{length}}$ is a uniform distribution along an $L = 20^{\circ}$ stream fragment, $\mathcal{P}_{\text{width}}$ is a Gaussian width of 50 pc dispersion, $\mathcal{P}_{\text{LF}}$ is the probability of being drawn from the luminosity function of the adopted simple stellar population (SSP) model (we adopt the PARSEC isochrone models, Bressan et al. 2012), $\mathcal{P}_{\text{color}}$ is a probability of the measured $G_{\text{BP}} - G_{\text{RP}}$ color of the star at the observed $G$ magnitude, and $\mathcal{P}_{\mu}$ and $\mathcal{P}_{\varpi}$ are the Gaussian probabilities of the observed astrometry given an assumed 5 km s$^{-1}$ intrinsic stream velocity dispersion, and a 50 pc line-of-sight dispersion. (The proper motion probability term $\mathcal{P}_{\mu}$ is a two-dimensional Gaussian PDF incorporating the proper motion uncertainties in R.A. and decl. and their correlation cross-terms.) Up to here, all is identical to the model described in Paper II. However, since the DR3 catalog provides many line-of-sight velocity measurements, we now also include a simple Gaussian probability $\mathcal{P}_{\nu}$ that the star has a line-of-sight velocity within 5 km s$^{-1}$ of the orbit. Most stars in the DR3 catalog do not have measured velocities, and for those, we set $\mathcal{P}_{\nu} = 1$.

In an earlier contribution we used a much wider (500 pc) spatial template in order to detect the stream of the Sagittarius dwarf galaxy and any similar structures (Ibata et al. 2020a). In contrast, our decision in the present work to adopt a characteristic 50 pc width spread and 5 km s$^{-1}$ velocity spread for the template was made to facilitate the detection of globular cluster streams. In the Harris (2010) catalog of Milky Way globular clusters, 79% have tidal radii <50 pc (their Gaussian width is, of course, much smaller), and 78% have a central velocity dispersion <5 km s$^{-1}$. Streams with much wider spatial or velocity dispersions will have their significance downweighted by our template choice; in future work, we intend to overcome this limitation by using a wider range of templates.

Including velocities in the likelihood function of the stream model means that we also have to include them in the contamination model $\mathcal{P}_{\text{cont}}$. We decided to retain the previously fitted, completely empirical, Gaussian mixture model (GMM) description $\mathcal{P}_{\text{cont}}^{\text{Paper II}}(\alpha, \delta, \mu_{\alpha}, \mu_{\delta}, G, G_{\text{BP}} - G_{\text{RP}}, \varpi)$ but include a new term $\mathcal{P}_{\text{GMM}}$ as follows:

$$ \begin{aligned} \mathcal{P}_{\text{cont}}(\alpha, \delta, \mu_{\alpha}, \mu_{\delta}, G, G_{\text{BP}} - G_{\text{RP}}, \varpi, \nu) \\ &= \mathcal{P}_{\text{cont}}^{\text{Paper II}}(\alpha, \delta, \mu_{\alpha}, \mu_{\delta}, G, G_{\text{BP}} - G_{\text{RP}}, \varpi) \\ & \quad \times \mathcal{P}_{\text{GMM}}(\nu|\alpha, \delta, G, G_{\text{BP}} - G_{\text{RP}}, \varpi). \end{aligned} \quad (3) $$

Note that the final term ignores the correlations of line-of-sight velocity with proper motion, which we felt was a necessary simplifying assumption given the relative scarcity of the velocity information compared to the astrometric information.
This $P_{GMM}$ term is calculated from the available velocities in large 5°.6 × 5°.6 tiles of a zenithal equal area projection of the sky.

The bulk of the radial velocity data that we provide to the algorithm are the 33 million stars with radial velocity measurements published as part of Gaia DR3. We complemented these measurements with other large public spectroscopic surveys, crossmatching the full Gaia DR3 survey with the APOGEE-2 survey (Majewski et al. 2017), the GALAH DR3 survey (Buder et al. 2018), the LAMOST DR7 survey (Cui et al. 2012), the Radial Velocity Experiment DR5 (Kunder et al. 2017), the Sloan Digital Sky Survey (SDSS)/Segue survey (Yanny et al. 2009), the Gaia-ESO survey (Randich et al. 2022), and the S5 survey (Li et al. 2019). We also included measurements of Palomar 5 stars by Odenkirchen et al. (2009) and Ibata et al. (2017a), as well as our previous velocity measurements of STREAMFINDER targets with CFHT/ESPaDOnS and the New Technology Telescope (NTT)/ EFOSC2 and VLT/UVES presented in Paper II. For those sources where more than one measurement was available, we adopted the value with the lowest uncertainty.

A further difference with respect to Paper II is that we now include stars from the Gaia DR3 RR Lyrae variable catalog (Clementini et al. 2022). Although we used the VizieR database to recover any measured line-of-sight velocity of these stars (where available) we decided to ignore this information in the STREAMFINDER analysis, as RR Lyrae are very strong radial velocity variables, and we could not reliably phase-correct the heterogeneous archival velocity measurements, so for RR Lyrae we set $P_v = 1$ in Equation (2). For these stars, we use the Muraveva et al. (2018) metallicity G-band absolute magnitude relation to predict a distance given the SSP trial metallicity, and set $P_{LF} = 1$ in Equation (2).

The algorithm proceeds star by star through the Gaia survey, using each star as a reference to launch trial orbits within the adopted fiducial potential (Dehnen & Binney 1998). The potential is actually going to be finely constrained with the detected streams in the following sections, but see Malhan & Ibata (2018) on the minor influence that a wrong potential in the first step has on the stream detection itself. The measured astrometry and photometry of the reference star are used to generate plausible trial orbits, and the assumed stellar population model provides the link between photometry and distance. When the line-of-sight velocity is available, it is used; otherwise, we scan through the missing information as explained in Papers I and II. For each trial orbit of the star, we find (and record) all the neighbors within the (10°) search radius that can be associated within a chosen threshold of the stream model. We define k to be the number of stars found in this way around a trial orbit centered on a particular star.

In many fields of astronomy, local measures of significance are used to define a detection, as the sources typically appear small (e.g., a star or distant galaxy) compared to the field of view of the instrument or survey, and there is a well-defined background. However, in the present situation, we are interested in detecting structures that subtend very long arcs on the sky with a wide range of (typically) foreground contamination. For streams we therefore seek a global statistic to quantify detection. The log-likelihood of Equation (1) can provide such a measure, but we would need to sum over the full Gaia survey (with more than $10^9$ stars) in order to compare models of different structures. Note also that we would need to

render the model more complex with each added stream; doing this in the Bayesian framework would require us to calculate the Bayesian evidence for each additional candidate structure to verify that it merits being added. This would be prohibitively computationally expensive. To simplify the calculation, we break up the log-likelihood into two terms $\ln \mathcal{L} = \ln \mathcal{L}_1 + \ln \mathcal{L}_2$, where $\ln \mathcal{L}_1$ corresponds to the parameter space region $R_1$ containing k stars (within, say $3\sigma$ of the model in the orbital and stellar populations parameters) and $\ln \mathcal{L}_2$ is for the complementary region $R_2$ with $n-k$ stars, where n is the number of stars in $R_1 + R_2$. Now for region $R_1$:

$$ \ln \mathcal{L}_1 = \sum_{j=1}^{k} \ln (\eta \mathcal{P}_{j,\text{model}} + (1 - \eta) \mathcal{P}_{j,\text{cont}}) - \ln \mathcal{P}_{j,\text{cont}}, \quad (4) $$

where the final term provides a convenient zero-point so that $\ln \mathcal{L}_1 = 0$ if no stream is present. In region $R_2$ we make the simplification that the stream model probability is zero, so the corresponding log-likelihood is

$$ \begin{aligned} \ln \mathcal{L}_2 &\approx \sum_{j=1}^{n-k} \ln((1 - \eta)\mathcal{P}_{j,\text{cont}}) - \ln \mathcal{P}_{j,\text{cont}} \\ &= (n - k)\ln(1 - \eta), \end{aligned} \quad (5) $$

where again we have included a zero-point choice so that $\ln \mathcal{L}_2 = 0$ when $\eta = 0$. Adding Equations (4) and (5) gives the total log-likelihood of

$$ \ln \mathcal{L} = (n - k)\ln(1 - \eta) + \sum_{j=1}^{k} \ln \left( 1 - \eta + \eta \frac{\mathcal{P}_{j,\text{model}}}{\mathcal{P}_{j,\text{cont}}} \right), \quad (6) $$

which allows a straightforward statistical comparison of the detected features over the full sky, and has the convenient property of being zero when $\eta = 0$.

Previously, in Papers I and II, we selected the stream solution $\eta$ with the highest likelihood according to Equation (6), which is a choice that works well at high Galactic latitude, but becomes problematic at low latitude where the number of possible associations becomes very large. We remind the reader that our algorithm is not finding the best single model fit to the global Gaia data, but rather our procedure checks whether a large-scale stream model centered at the phase-space position of each of the billions of stars in the survey is statistically significant or not. Hence, a given star may be considered a member of many different proposed groups. This is, of course, very different from the usual fitting problems in astronomy, where we fit the parameters of a model to a data set, or in detection problems, where we say fit galaxy models to a large number of spatially independent regions of a pixellated image. For the present work we therefore decided to use the logarithm of the likelihood per star of the group as our selection statistic

$$ \ln \hat{\mathcal{L}} = \ln \mathcal{L} - k, \quad (7) $$

where k is again the number of stars associated with the trial orbit.<fnref n="12" /> Optimizing this statistic forces the algorithm to find the

<footnote n="12">¹² In the present context, the statistic in Equation (7) is reminiscent of the Akaike information criterion, $AIC = 2m - 2 \ln \mathcal{L}$, which is a statistical estimator designed to optimize the balance between goodness of fit and model complexity (m is the number of parameters of the model). Because of the fact that stars can be reused in different stream solutions, the higher k is on average, the more potential overlapping models there will be, and so the higher the complexity of the global model will be.</footnote>
This $P_{GMM}$ term is calculated from the available velocities in large 5°6 × 5°6 tiles of a zenithal equal area projection of the sky.

The bulk of the radial velocity data that we provide to the algorithm are the 33 million stars with radial velocity measurements published as part of Gaia DR3. We complemented these measurements with other large public spectroscopic surveys, crossmatching the full Gaia DR3 survey with the APOGEE-2 survey (Majewski et al. 2017), the GALAH DR3 survey (Buder et al. 2018), the LAMOST DR7 survey (Cui et al. 2012), the Radial Velocity Experiment DR5 (Kunder et al. 2017), the Sloan Digital Sky Survey (SDSS)/Segue survey (Yanny et al. 2009), the Gaia-ESO survey (Randich et al. 2022), and the S5 survey (Li et al. 2019). We also included measurements of Palomar 5 stars by Odenkirchen et al. (2009) and Ibata et al. (2017a), as well as our previous velocity measurements of STREAMFINDER targets with CFHT/ESPaDOnS and the New Technology Telescope (NTT)/ EFOSC2 and VLT/UVES presented in Paper II. For those sources where more than one measurement was available, we adopted the value with the lowest uncertainty.

A further difference with respect to Paper II is that we now include stars from the Gaia DR3 RR Lyrae variable catalog (Clementini et al. 2022). Although we used the VizieR database to recover any measured line-of-sight velocity of these stars (where available) we decided to ignore this information in the STREAMFINDER analysis, as RR Lyrae are very strong radial velocity variables, and we could not reliably phase-correct the heterogeneous archival velocity measurements, so for RR Lyrae we set $P_v = 1$ in Equation (2). For these stars, we use the Muraveva et al. (2018) metallicity G-band absolute magnitude relation to predict a distance given the SSP trial metallicity, and set $P_{LF} = 1$ in Equation (2).

The algorithm proceeds star by star through the Gaia survey, using each star as a reference to launch trial orbits within the adopted fiducial potential (Dehnen & Binney 1998). The potential is actually going to be finely constrained with the detected streams in the following sections, but see Malhan & Ibata (2018) on the minor influence that a wrong potential in the first step has on the stream detection itself. The measured astrometry and photometry of the reference star are used to generate plausible trial orbits, and the assumed stellar population model provides the link between photometry and distance. When the line-of-sight velocity is available, it is used; otherwise, we scan through the missing information as explained in Papers I and II. For each trial orbit of the star, we find (and record) all the neighbors within the (10°) search radius that can be associated within a chosen threshold of the stream model. We define k to be the number of stars found in this way around a trial orbit centered on a particular star.

In many fields of astronomy, local measures of significance are used to define a detection, as the sources typically appear small (e.g., a star or distant galaxy) compared to the field of view of the instrument or survey, and there is a well-defined background. However, in the present situation, we are interested in detecting structures that subtend very long arcs on the sky with a wide range of (typically) foreground contamination. For streams we therefore seek a global statistic to quantify detection. The log-likelihood of Equation (1) can provide such a measure, but we would need to sum over the full Gaia survey (with more than $10^9$ stars) in order to compare models of different structures. Note also that we would need to

render the model more complex with each added stream; doing this in the Bayesian framework would require us to calculate the Bayesian evidence for each additional candidate structure to verify that it merits being added. This would be prohibitively computationally expensive. To simplify the calculation, we break up the log-likelihood into two terms $\ln \mathcal{L} = \ln \mathcal{L}_1 + \ln \mathcal{L}_2$, where $\ln \mathcal{L}_1$ corresponds to the parameter space region $R_1$ containing k stars (within, say $3\sigma$ of the model in the orbital and stellar populations parameters) and $\ln \mathcal{L}_2$ is for the complementary region $R_2$ with $n-k$ stars, where n is the number of stars in $R_1 + R_2$. Now for region $R_1$:

$$ \ln \mathcal{L}_1 = \sum_{j=1}^{k} \ln (\eta P_{j,\text{model}} + (1-\eta) P_{j,\text{cont}}) - \ln P_{j,\text{cont}}, \quad (4) $$

where the final term provides a convenient zero-point so that $\ln \mathcal{L}_1 = 0$ if no stream is present. In region $R_2$ we make the simplification that the stream model probability is zero, so the corresponding log-likelihood is

$$ \begin{aligned} \ln \mathcal{L}_2 &\approx \sum_{j=1}^{n-k} \ln((1-\eta)P_{j,\text{cont}}) - \ln P_{j,\text{cont}} \\ &= (n-k)\ln(1-\eta), \end{aligned} \quad (5) $$

where again we have included a zero-point choice so that $\ln \mathcal{L}_2 = 0$ when $\eta = 0$. Adding Equations (4) and (5) gives the total log-likelihood of

$$ \ln \mathcal{L} = (n-k)\ln(1-\eta) + \sum_{j=1}^{k} \ln \left( 1 - \eta + \eta \frac{P_{j,\text{model}}}{P_{j,\text{cont}}} \right), \quad (6) $$

which allows a straightforward statistical comparison of the detected features over the full sky, and has the convenient property of being zero when $\eta = 0$.

Previously, in Papers I and II, we selected the stream solution $\eta$ with the highest likelihood according to Equation (6), which is a choice that works well at high Galactic latitude, but becomes problematic at low latitude where the number of possible associations becomes very large. We remind the reader that our algorithm is not finding the best single model fit to the global Gaia data, but rather our procedure checks whether a large-scale stream model centered at the phase-space position of each of the billions of stars in the survey is statistically significant or not. Hence, a given star may be considered a member of many different proposed groups. This is, of course, very different from the usual fitting problems in astronomy, where we fit the parameters of a model to a data set, or in detection problems, where we say fit galaxy models to a large number of spatially independent regions of a pixellated image. For the present work we therefore decided to use the logarithm of the likelihood per star of the group as our selection statistic

$$ \ln \hat{\mathcal{L}} = \ln \mathcal{L} - k, \quad (7) $$

where k is again the number of stars associated with the trial orbit.<fnref n="12" /> Optimizing this statistic forces the algorithm to find the

<footnote n="12">¹² In the present context, the statistic in Equation (7) is reminiscent of the Akaike information criterion, $AIC = 2m - 2 \ln \mathcal{L}$, which is a statistical estimator designed to optimize the balance between goodness of fit and model complexity (m is the number of parameters of the model). Because of the fact that stars can be reused in different stream solutions, the higher k is on average, the more potential overlapping models there will be, and so the higher the complexity of the global model will be.</footnote>
orbit and associated stream fraction $\eta$ that maximizes the number of high-likelihood associations. By adjusting the acceptance threshold, stream maps built with this $\ln \hat{L}$ statistic can be selected to be very similar at high Galactic latitude to those built previously using the $\ln L$ statistic. However, close to the plane of the Galaxy, the new maps do not show spurious solutions that possess a large number of low probability members (i.e., that are composed mostly of stars at a high number of standard deviations from the orbit and stellar populations model).

### 3. Spectroscopic Observations

We used the VLT/UVES spectrograph (D'Oderico et al. 2000) to follow up on selected STREAMFINDER sources detected in the Gaia EDR3 and DR3 catalogs. These runs comprise runs 105.20AL.001 (2.5 nights in visitor mode), 110.246A.001 (40 hr service in service mode), and 111.2517.001 (3.6 nights in visitor mode). Our instrumental setup uses the DIC2 dichroic beamsplitter in the "437+760" setting, covering the wavelength ranges 3730–4990 and 5650–9460 Å. To reduce read noise, we binned the CCD in $2 \times 2$ pixel blocks, which in conjunction with a 1"/0 wide slit, yields a spectral resolution of approximately 40,000. Exposure times were selected on a star-by-star basis to reach signal-to-noise ratio ~3–5 for the fainter stars in the sample, so as to measure their radial velocities, but we set a minimum exposure time of 5 minutes. For the brighter stars, this minimum exposure time also allowed some elemental abundances to be measured. The spectra were reduced with the EsoReflex pipeline using daytime calibration arc lamps and flat-field images, resulting in extracted wavelength-calibrated one-dimensional spectra.

We also secured observations with the IDS long-slit spectrograph on the 2.5 m Isaac Newton Telescope in several runs over the course of 2022. Bright northern hemisphere stream stars were targeted with typically ~1 hr exposures at $G = 16$ mag. The instrument was configured with the RED+2 detector, the R1200R grating with a central wavelength set to 8500 Å, a 1" wide slit, and the GG495 order sorting filter.

The radial velocities of the target stars were measured with the IRAF `fxcor` algorithm, using the bright and relatively metal-poor star HD 182572 as a radial velocity standard. The UVES spectra were of sufficient quality and resolution to obtain excellent radial velocity measurements with $<1 \text{ km s}^{-1}$ uncertainty for stars to $G = 18$ mag, while the IDS spectra resulted in velocity uncertainties of ~10 km s$^{-1}$ at $G = 16$ mag.

### 4. Gaia DR3 STREAMFINDER Maps

We ran STREAMFINDER using seven different SSP templates with metallicity $[\text{Fe/H}] = -2.2, -1.9, -1.7, -1.5, -1.3, -1.1, -0.7$, and a common age 12.5 Gyr, as we are primarily aiming to discover halo structures. We allowed the solutions to explore the distance range of [1, 100] kpc. As in Paper II, the Gaia catalog was limited to sources with $G_0 < 20$ mag and V-band extinction $A_V < 2$ mag, so as to mitigate against variations in extinction and survey depth due to the satellite's scanning law.<fnref n="13" /> The full sky was processed

with the exception of circular regions around satellites; the radius of these masked regions was set to twice the tidal radius for globular clusters and seven times the half-light radius for satellite galaxies, as detailed in Paper II.

Figures 1 and 2 show the resulting sky maps of the detections of 24,540 (non-RR Lyrae) sources (listed in Table 1) with $\ln \hat{L} > 19.8$ (note that with the standard log-likelihood, the value $\ln L = 19.8$ corresponds to a $6\sigma$ threshold) colored according to proper motion along the Galactic longitude and latitude directions, respectively. The corresponding distance solutions are shown in Figure 3, colored according to the distance modulus. Figure 4 displays the 2178 (non-RR Lyrae) sources with measured line-of-sight velocities, while Figure 5 shows the mean metallicity of the streams (listed in Table 3) as calculated from the mean of the spectroscopic metallicity measurements (or photometric metallicities, for those structures without spectroscopic measurements). These maps display a rich tapestry of crisscrossing streams that possess coherent behavior in proper motion, line-of-sight velocity, and distance. In Figure 6, we have selected a contrasting color scheme to allow easier visual disambiguation of the overlapping structures that can be discerned in Figures 1–5.

To facilitate subsequent analyses, we create a local coordinate system on a per stream basis by fitting a great circle to each structure using a least squares criterion (we consider that the uncertainties in sky position are negligible). For the streams with known globular cluster progenitors, we require the fit to intersect the center of the globular cluster, and we also use the R.A. of the cluster to define the zero-point of the great circle coordinates $(\phi_1, \phi_2)$ along and perpendicular to the stream. The coordinate $\phi_1$ is oriented so that the stream's velocity is approximately parallel to it,<fnref n="14" /> while $\phi_2$ is in the direction perpendicular to the great circle. For stream #48 (Orphan-Chenab) we chose the zero-point proposed by Koposov et al. (2019), while for stream #53 (GD-1) we fully adopted the $(\phi_1, \phi_2)$ coordinate system of Koposov et al. (2010). The poles and zero-point of the $(\phi_1, \phi_2)$ coordinates are provided in Table 3. We also estimate the stellar mass of each stream, using the metallicity listed in Table 3 together with the fitted stream distance models (calculated as explained in Section 5), to correct for the missing stars in the stellar population, assuming that the survey is complete to $G_0 = 20$ mag. These stellar masses are also listed in Table 3, but we stress that they are rough approximations and involve many assumptions, including that we have identified the full extent of each stream and that only a single stellar population is present.

Figures 7–9 display the individual stream structures in an array of three columns. The first column shows the line-of-sight velocity as a function of the angle $\phi_1$ along each stream. The proper motion along the $\phi_1$ direction is displayed in the middle column, while the last column displays the color–magnitude diagram of the population. The same color scheme as in Figure 6 is employed.

The nature of the detections can be appreciated better by first considering their dynamical properties. Next, we use the streams to understand the global properties of the Milky Way, which will then give us access to their orbits, and then we return to a discussion of individual objects in Section 9.

<footnote n="13">¹³ To correct for interstellar extinction, we assume that all the dust is in the foreground and interpolate between the pixels of the Schlegel et al. (1998) maps, adopting the Schlafly & Finkbeiner (2011) recalibration with $R_V = 3.1$.</footnote>

<footnote n="14">¹⁴ The stream velocity vector need not be perfectly parallel to $\phi_1$ due to projection and reflex motion.</footnote>
orbit and associated stream fraction $\eta$ that maximizes the number of high-likelihood associations. By adjusting the acceptance threshold, stream maps built with this $\ln \hat{L}$ statistic can be selected to be very similar at high Galactic latitude to those built previously using the $\ln L$ statistic. However, close to the plane of the Galaxy, the new maps do not show spurious solutions that possess a large number of low probability members (i.e., that are composed mostly of stars at a high number of standard deviations from the orbit and stellar populations model).

### 3. Spectroscopic Observations

We used the VLT/UVES spectrograph (D'Oderico et al. 2000) to follow up on selected STREAMFINDER sources detected in the Gaia EDR3 and DR3 catalogs. These runs comprise runs 105.20AL.001 (2.5 nights in visitor mode), 110.246A.001 (40 hr service in service mode), and 111.2517.001 (3.6 nights in visitor mode). Our instrumental setup uses the DIC2 dichroic beamsplitter in the "437+760" setting, covering the wavelength ranges 3730–4990 and 5650–9460 Å. To reduce read noise, we binned the CCD in 2 × 2 pixel blocks, which in conjunction with a 1"/0 wide slit, yields a spectral resolution of approximately 40,000. Exposure times were selected on a star-by-star basis to reach signal-to-noise ratio ~3–5 for the fainter stars in the sample, so as to measure their radial velocities, but we set a minimum exposure time of 5 minutes. For the brighter stars, this minimum exposure time also allowed some elemental abundances to be measured. The spectra were reduced with the EsoReflex pipeline using daytime calibration arc lamps and flat-field images, resulting in extracted wavelength-calibrated one-dimensional spectra.

We also secured observations with the IDS long-slit spectrograph on the 2.5 m Isaac Newton Telescope in several runs over the course of 2022. Bright northern hemisphere stream stars were targeted with typically ~1 hr exposures at $G = 16$ mag. The instrument was configured with the RED+2 detector, the R1200R grating with a central wavelength set to 8500 Å, a 1" wide slit, and the GG495 order sorting filter.

The radial velocities of the target stars were measured with the IRAF `fxcor` algorithm, using the bright and relatively metal-poor star HD 182572 as a radial velocity standard. The UVES spectra were of sufficient quality and resolution to obtain excellent radial velocity measurements with <1 km s⁻¹ uncertainty for stars to $G = 18$ mag, while the IDS spectra resulted in velocity uncertainties of ~10 km s⁻¹ at $G = 16$ mag.

### 4. Gaia DR3 STREAMFINDER Maps

We ran STREAMFINDER using seven different SSP templates with metallicity [Fe/H] = −2.2, −1.9, −1.7, −1.5, −1.3, −1.1, −0.7, and a common age 12.5 Gyr, as we are primarily aiming to discover halo structures. We allowed the solutions to explore the distance range of [1, 100] kpc. As in Paper II, the Gaia catalog was limited to sources with $G_0 < 20$ mag and V-band extinction $A_V < 2$ mag, so as to mitigate against variations in extinction and survey depth due to the satellite's scanning law.<fnref n="13" /> The full sky was processed

with the exception of circular regions around satellites; the radius of these masked regions was set to twice the tidal radius for globular clusters and seven times the half-light radius for satellite galaxies, as detailed in Paper II.

Figures 1 and 2 show the resulting sky maps of the detections of 24,540 (non-RR Lyrae) sources (listed in Table 1) with $\ln \hat{L} > 19.8$ (note that with the standard log-likelihood, the value $\ln L = 19.8$ corresponds to a $6\sigma$ threshold) colored according to proper motion along the Galactic longitude and latitude directions, respectively. The corresponding distance solutions are shown in Figure 3, colored according to the distance modulus. Figure 4 displays the 2178 (non-RR Lyrae) sources with measured line-of-sight velocities, while Figure 5 shows the mean metallicity of the streams (listed in Table 3) as calculated from the mean of the spectroscopic metallicity measurements (or photometric metallicities, for those structures without spectroscopic measurements). These maps display a rich tapestry of crisscrossing streams that possess coherent behavior in proper motion, line-of-sight velocity, and distance. In Figure 6, we have selected a contrasting color scheme to allow easier visual disambiguation of the overlapping structures that can be discerned in Figures 1–5.

To facilitate subsequent analyses, we create a local coordinate system on a per stream basis by fitting a great circle to each structure using a least squares criterion (we consider that the uncertainties in sky position are negligible). For the streams with known globular cluster progenitors, we require the fit to intersect the center of the globular cluster, and we also use the R.A. of the cluster to define the zero-point of the great circle coordinates ($\phi_1, \phi_2$) along and perpendicular to the stream. The coordinate $\phi_1$ is oriented so that the stream's velocity is approximately parallel to it,<fnref n="14" /> while $\phi_2$ is in the direction perpendicular to the great circle. For stream #48 (Orphan-Chenab) we chose the zero-point proposed by Koposov et al. (2019), while for stream #53 (GD-1) we fully adopted the ($\phi_1, \phi_2$) coordinate system of Koposov et al. (2010). The poles and zero-point of the ($\phi_1, \phi_2$) coordinates are provided in Table 3. We also estimate the stellar mass of each stream, using the metallicity listed in Table 3 together with the fitted stream distance models (calculated as explained in Section 5), to correct for the missing stars in the stellar population, assuming that the survey is complete to $G_0 = 20$ mag. These stellar masses are also listed in Table 3, but we stress that they are rough approximations and involve many assumptions, including that we have identified the full extent of each stream and that only a single stellar population is present.

Figures 7–9 display the individual stream structures in an array of three columns. The first column shows the line-of-sight velocity as a function of the angle $\phi_1$ along each stream. The proper motion along the $\phi_1$ direction is displayed in the middle column, while the last column displays the color–magnitude diagram of the population. The same color scheme as in Figure 6 is employed.

The nature of the detections can be appreciated better by first considering their dynamical properties. Next, we use the streams to understand the global properties of the Milky Way, which will then give us access to their orbits, and then we return to a discussion of individual objects in Section 9.

<footnote n="13">¹³ To correct for interstellar extinction, we assume that all the dust is in the foreground and interpolate between the pixels of the Schlegel et al. (1998) maps, adopting the Schlafly & Finkbeiner (2011) recalibration with $R_V = 3.1$.</footnote>

<footnote n="14">¹⁴ The stream velocity vector need not be perfectly parallel to $\phi_1$ due to projection and reflex motion.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Four different projections in Galactic coordinates of the 24,540 STREAMFINDER sources listed in Table 1 with $\ln \hat{\mathcal{L}} > 19.8$ (more restrictive than a $6\sigma$ threshold). The color of the points encodes the proper motion $\mu_t$ along the Galactic longitude direction, as measured by Gaia in the DR3 catalog.</figcaption></figure>

**5. Stream Fitting with Corrected Orbits in an Isolated Milky Way**

We begin our modeling efforts by considering the Milky Way to be an isolated axisymmetric galaxy, whose center defines the origin of an inertial frame around which its satellites will orbit. In this context, we seek to find the most likely Galactic mass model that admits progenitor orbits that fit the stream data simultaneously. As we have mentioned before, stellar streams do not precisely delineate the orbital path of their progenitors. To overcome this complication, we will proceed in an iterative manner to find plausible functions $\Delta_{\Theta}(\phi_1)$ in the derived mass model that correct the offset between the stream and the progenitor orbit at $\phi_1$ for each observable phase-space coordinate $\Theta$.

While we could have tried to fit the whole data set of 87 streams and $>24,500$ stars, to find the best mass model, we restrict ourselves in the present contribution to a conservative sample containing velocity-confirmed members and RR Lyrae. We further impose that the individual streams should be long (>10° in length), containing at least five line-of-sight velocity measurements with a clear linear trend. Since we most probably only detect a small stream segment of the closest streams, we also select those structures at a heliocentric distance $>2.5$ kpc (this limit is somewhat arbitrary, and in future work, we will reexamine in detail the orbits of the closer streams). We set aside the Orphan stream for the initial analysis because although it has excellent data, it is probably affected by the Large Magellanic Cloud (LMC; Erkal et al. 2019). We have
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Four different projections in Galactic coordinates of the 24,540 STREAMFINDER sources listed in Table 1 with $\ln \hat{\mathcal{L}} > 19.8$ (more restrictive than a $6\sigma$ threshold). The color of the points encodes the proper motion $\mu_t$ along the Galactic longitude direction, as measured by Gaia in the DR3 catalog.</figcaption></figure>

**5. Stream Fitting with Corrected Orbits in an Isolated Milky Way**

We begin our modeling efforts by considering the Milky Way to be an isolated axisymmetric galaxy, whose center defines the origin of an inertial frame around which its satellites will orbit. In this context, we seek to find the most likely Galactic mass model that admits progenitor orbits that fit the stream data simultaneously. As we have mentioned before, stellar streams do not precisely delineate the orbital path of their progenitors. To overcome this complication, we will proceed in an iterative manner to find plausible functions $\Delta_{\Theta}(\phi_1)$ in the derived mass model that correct the offset between the stream and the progenitor orbit at $\phi_1$ for each observable phase-space coordinate $\Theta$.

While we could have tried to fit the whole data set of 87 streams and $>24,500$ stars, to find the best mass model, we restrict ourselves in the present contribution to a conservative sample containing velocity-confirmed members and RR Lyrae. We further impose that the individual streams should be long (>10° in length), containing at least five line-of-sight velocity measurements with a clear linear trend. Since we most probably only detect a small stream segment of the closest streams, we also select those structures at a heliocentric distance $>2.5$ kpc (this limit is somewhat arbitrary, and in future work, we will reexamine in detail the orbits of the closer streams). We set aside the Orphan stream for the initial analysis because although it has excellent data, it is probably affected by the Large Magellanic Cloud (LMC; Erkal et al. 2019). We have
<figure><img src="image_5.png" /><figcaption>**Figure 2.** As in Figure 1, but showing the measured proper motion $\mu_b$ of the sources along the Galactic longitude direction.</figcaption></figure>

previously shown (Ibata et al. 2019a) that the Fimbulthul stream possesses a complex structure due to the internal rotation of its progenitor, the massive cluster ω Cen, and so we decided not to model this structure either. The sample of 29 long and thin streams we fit here are flagged “1” in the column titled “sample” in Table 3. The data set comprises 1397 data points: 1098 stars with measured line-of-sight velocities, five globular clusters (with phase-space parameters derived from the compendia in Vasiliev 2019 and Baumgardt & Vasiliev 2021), and 294 Gaia RR Lyrae (for which we ignore any measured line-of-sight velocities due to concerns over their variability).

We constructed a similar Markov Chain Monte Carlo (MCMC) orbit-fitting algorithm to that presented in Ibata

et al. (2018), but altered to allow for simultaneous adjustment to many streams. The MCMC driver algorithm we use was originally developed by Ibata et al. (2011), based on the *walkers* scheme proposed by Goodman & Weare (2010). The *walkers* algorithm (that simulates a population of agents sampling the likelihood at different parameter space locations) is very efficient as it is affine-invariant and adapts automatically to the multidimensional covariance of the problem.

The Milky Way potential is generated via multipole expansion from analytical double-exponential density models (to represent the disk components) or double-power-law density models (to represent the spheroidal components), following Dehnen & Binney (1998). We allow each disk component to have up to three parameters: central surface
<figure><img src="image_5.png" /><figcaption>**Figure 2.** As in Figure 1, but showing the measured proper motion $\mu_b$ of the sources along the Galactic longitude direction.</figcaption></figure>

previously shown (Ibata et al. 2019a) that the Fimbulthul stream possesses a complex structure due to the internal rotation of its progenitor, the massive cluster ω Cen, and so we decided not to model this structure either. The sample of 29 long and thin streams we fit here are flagged “1” in the column titled “sample” in Table 3. The data set comprises 1397 data points: 1098 stars with measured line-of-sight velocities, five globular clusters (with phase-space parameters derived from the compendia in Vasiliev 2019 and Baumgardt & Vasiliev 2021), and 294 Gaia RR Lyrae (for which we ignore any measured line-of-sight velocities due to concerns over their variability).

We constructed a similar Markov Chain Monte Carlo (MCMC) orbit-fitting algorithm to that presented in Ibata

et al. (2018), but altered to allow for simultaneous adjustment to many streams. The MCMC driver algorithm we use was originally developed by Ibata et al. (2011), based on the *walkers* scheme proposed by Goodman & Weare (2010). The *walkers* algorithm (that simulates a population of agents sampling the likelihood at different parameter space locations) is very efficient as it is affine-invariant and adapts automatically to the multidimensional covariance of the problem.

The Milky Way potential is generated via multipole expansion from analytical double-exponential density models (to represent the disk components) or double-power-law density models (to represent the spheroidal components), following Dehnen & Binney (1998). We allow each disk component to have up to three parameters: central surface
<figure><img src="image_6.png" /><figcaption>Figure 3. As in Figure 1, but showing the distance solutions, displayed as distance modulus.</figcaption></figure>

density $\Sigma_d$, a scale length $h_R$ and scale height $h_z$

$$ \rho_d(R, z) = \frac{\Sigma_d}{2h_z} \exp\left(-\frac{R}{h_R} - \frac{|z|}{h_z}\right). \qquad (8) $$

In contrast we allow the spheroidal components to have six parameters: central density $\rho_0$, density flattening $q_m$, inner power-law slope $-\beta$, outer power-law slope $-\gamma$, scale radius $r_0$, and truncation radius $r_t$

$$ \rho_s(s) = \rho_0 \left(\frac{s}{r_0}\right)^{-\gamma} \left(1 + \frac{s}{r_0}\right)^{\gamma-\beta} e^{-s^2/r_t^2}, \qquad (9) $$

where $s = \sqrt{R^2 + z^2/q_m^2}$ is an ellipsoidal coordinate. Each orbit obviously requires six parameters; however, for those

streams with a known progenitor, we consider their center to be perfectly determined, leaving four free parameters, while for all others, we anchor the orbit at $\phi_1 = 0^\circ$ (except for GD-1, as discussed below), leaving five free parameters.

Hence, a full parameter exploration of the problem would have in excess of 150 parameters and would therefore be very challenging for an MCMC approach. However, the orbits can be fit independently of one another, which motivated our choice to implement a nested MCMC search. The outer MCMC loop simply proposes the Galaxy model parameters (and the components of the Solar peculiar motion vector, if so desired). The inner loop adopts those Galaxy model parameters and attempts to fit each stream independently. We implemented an option for two algorithm variants for the inner loop. The first option is another independent MCMC search with 500
<figure><img src="image_6.png" /><figcaption>Figure 3. As in Figure 1, but showing the distance solutions, displayed as distance modulus.</figcaption></figure>

density $\Sigma_d$, a scale length $h_R$ and scale height $h_z$

$$ \rho_d(R, z) = \frac{\Sigma_d}{2h_z} \exp\left(-\frac{R}{h_R} - \frac{|z|}{h_z}\right). \qquad (8) $$

In contrast we allow the spheroidal components to have six parameters: central density $\rho_0$, density flattening $q_m$, inner power-law slope $-\beta$, outer power-law slope $-\gamma$, scale radius $r_0$, and truncation radius $r_t$

$$ \rho_s(s) = \rho_0 \left(\frac{s}{r_0}\right)^{-\gamma} \left(1 + \frac{s}{r_0}\right)^{\gamma-\beta} e^{-s^2/r_t^2}, \qquad (9) $$

where $s = \sqrt{R^2 + z^2/q_m^2}$ is an ellipsoidal coordinate. Each orbit obviously requires six parameters; however, for those

streams with a known progenitor, we consider their center to be perfectly determined, leaving four free parameters, while for all others, we anchor the orbit at $\phi_1 = 0^\circ$ (except for GD-1, as discussed below), leaving five free parameters.

Hence, a full parameter exploration of the problem would have in excess of 150 parameters and would therefore be very challenging for an MCMC approach. However, the orbits can be fit independently of one another, which motivated our choice to implement a nested MCMC search. The outer MCMC loop simply proposes the Galaxy model parameters (and the components of the Solar peculiar motion vector, if so desired). The inner loop adopts those Galaxy model parameters and attempts to fit each stream independently. We implemented an option for two algorithm variants for the inner loop. The first option is another independent MCMC search with 500
<figure><img src="image_7.png" /><figcaption>**Figure 4.** As in Figure 1, but showing the line-of-sight velocity measurements of the stars in the sample. Stars without such measurements are marked with small gray points.</figcaption></figure>

iterations, where every 25 iterations of the inner loop, the step size is adjusted to try to maintain a 25% acceptance ratio (independently for each stream). The best solution and step size (for each stream) are recorded as if they were parameters of the outer loop and reused in the next run of the inner loop. The second option is to use a Nelder-Mead *downhill simplex* search (Press et al. 1992). We adopted a log-likelihood tolerance of $10^{-4}$, which gave robust results in our tests (each such search requires ≈1000 likelihood evaluations). The Nelder-Mead method can be sensitive to the chosen starting point of the search, so we reran it with different randomly drawn positions until the best likelihood solution is found again to within the chosen tolerance. We found that both the MCMC and downhill simplex options for the inner loop give similar results, but the

latter is computationally cheaper to reach a given accuracy, while the former allows one to estimate confidence intervals on the stream parameters. The results shown below were calculated using the downhill simplex as the optimization method of the inner loop.<fnref n="15" />

The outer loop uses 96 MCMC walkers, and we adjust the stretch factor of the proposed walker steps dynamically (see Goodman & Weare 2010) every 25 iterations, aiming to obtain a 25% acceptance ratio. We run the algorithm for at least 200,000 *burn-in* iterations, discard that data, and then continue

<footnote n="15">¹⁵ Although we first employ the MCMC method in the inner loop to obtain an estimate of the uncertainties of the stream model parameters, which are then reused in the generation of the (Gaussian) random starting points for the downhill simplex method.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** As in Figure 1, but showing the line-of-sight velocity measurements of the stars in the sample. Stars without such measurements are marked with small gray points.</figcaption></figure>

iterations, where every 25 iterations of the inner loop, the step size is adjusted to try to maintain a 25% acceptance ratio (independently for each stream). The best solution and step size (for each stream) are recorded as if they were parameters of the outer loop and reused in the next run of the inner loop. The second option is to use a Nelder-Mead *downhill simplex* search (Press et al. 1992). We adopted a log-likelihood tolerance of $10^{-4}$, which gave robust results in our tests (each such search requires ≈1000 likelihood evaluations). The Nelder-Mead method can be sensitive to the chosen starting point of the search, so we reran it with different randomly drawn positions until the best likelihood solution is found again to within the chosen tolerance. We found that both the MCMC and downhill simplex options for the inner loop give similar results, but the

latter is computationally cheaper to reach a given accuracy, while the former allows one to estimate confidence intervals on the stream parameters. The results shown below were calculated using the downhill simplex as the optimization method of the inner loop.<fnref n="15" />

The outer loop uses 96 MCMC walkers, and we adjust the stretch factor of the proposed walker steps dynamically (see Goodman & Weare 2010) every 25 iterations, aiming to obtain a 25% acceptance ratio. We run the algorithm for at least 200,000 *burn-in* iterations, discard that data, and then continue

<footnote n="15">¹⁵ Although we first employ the MCMC method in the inner loop to obtain an estimate of the uncertainties of the stream model parameters, which are then reused in the generation of the (Gaussian) random starting points for the downhill simplex method.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** As in Figure 1, but showing the mean spectroscopic metallicity $\langle[\text{Fe/H}]\rangle$ (or photometric metallicity when spectroscopic values are not available) of the streams, as listed in Table 3.</figcaption></figure>

for a further 500,000 iterations (i.e., for a total of 500,000 × 1000 evaluations of the combined likelihood of the 29 streams, with each orbit being sampled at a resolution of up to 10,000 points). For every step of the outer loop, the Galaxy parameters and the stream parameters of that step are recorded as the MCMC chain, which is the desired final output of the algorithm.<fnref n="16" /> We consider that the MCMC exploration

converged when the parameter statistics of the first and second half of the 500,000 iteration-long chains were statistically identical. One such optimization of the Galaxy and stream parameters typically takes ≈5000 CPU hr.

The log-likelihood objective function we choose to fit to the sample attempts to model the expected combination of thin streams, each enveloped by a wider cocoon population and a small amount of contamination. The adopted objective function is therefore

$$ \ln \mathcal{L} = \sum_{\text{data}} \ln[ (1 - \zeta - \chi)\mathcal{P}_{\text{stream}}^{\text{fit}} + \zeta\mathcal{P}_{\text{cocoon}} + \chi\mathcal{P}_{\text{contam}} ] + \ln \mathcal{L}_{\text{ancillary}}, \qquad (10) $$

<footnote n="16">¹⁶ Due to the low cost of the individual orbit integrations, which nevertheless require a non-negligible amount of memory to store their paths, this nested procedure proved to be very difficult to parallelize efficiently with standard parallelization libraries such as OPENMMPMPI. Consequently, we developed a custom parallelization scheme using persistent jobs (eliminating overheads due to memory allocation and data ingestion), one per machine core, that wait to be handed the different input parameters farmed out by the master process that handles the MCMC exploration.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** As in Figure 1, but showing the mean spectroscopic metallicity $\langle [Fe/H] \rangle$ (or photometric metallicity when spectroscopic values are not available) of the streams, as listed in Table 3.</figcaption></figure>

for a further 500,000 iterations (i.e., for a total of 500,000 × 1000 evaluations of the combined likelihood of the 29 streams, with each orbit being sampled at a resolution of up to 10,000 points). For every step of the outer loop, the Galaxy parameters and the stream parameters of that step are recorded as the MCMC chain, which is the desired final output of the algorithm.<fnref n="16" /> We consider that the MCMC exploration

converged when the parameter statistics of the first and second half of the 500,000 iteration-long chains were statistically identical. One such optimization of the Galaxy and stream parameters typically takes ≈5000 CPU hr.

The log-likelihood objective function we choose to fit to the sample attempts to model the expected combination of thin streams, each enveloped by a wider cocoon population and a small amount of contamination. The adopted objective function is therefore

$$ \ln \mathcal{L} = \sum_{\text{data}} \ln[(1 - \zeta - \chi)\mathcal{P}_{\text{stream}}^{\text{fit}} + \zeta\mathcal{P}_{\text{cocoon}} + \chi\mathcal{P}_{\text{contam}}] + \ln \mathcal{L}_{\text{ancillary}}, \quad (10) $$

<footnote n="16">¹⁶ Due to the low cost of the individual orbit integrations, which nevertheless require a non-negligible amount of memory to store their paths, this nested procedure proved to be very difficult to parallelize efficiently with standard parallelization libraries such as OPENMMPMPI. Consequently, we developed a custom parallelization scheme using persistent jobs (eliminating overheads due to memory allocation and data ingestion), one per machine core, that wait to be handed the different input parameters farmed out by the master process that handles the MCMC exploration.</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 6. Gallery of the stellar streams in the same projections as in Figures 1–5. The color of the points shows our chosen stream identification number in a modulo 8 representation.</figcaption></figure>

where $\zeta$ and $\chi$ are the fraction of stars in the cocoon and contamination, respectively. We conservatively set $\chi = 0.01$. In our examination of the GD-1 stream in Section 9.4 below, we deduce $\zeta = 0.27$ for that system; however, to be conservative, we set $\zeta = 0.1$ for all streams in our MCMC analysis. The

stream model

$$P_{\text{stream}}^{\text{fit}} = P_{\text{width}} \times P_{\text{color}} \times P_{\mu} \times P_{\infty} \times P_{\nu} \quad (11)$$

is exactly the same stream probability criterion as used in STREAMFINDER (Equation (2)) but ignores the luminosity
<figure><img src="image_9.png" /><figcaption>Figure 6. Gallery of the stellar streams in the same projections as in Figures 1–5. The color of the points shows our chosen stream identification number in a modulo 8 representation.</figcaption></figure>

where $\zeta$ and $\chi$ are the fraction of stars in the cocoon and contamination, respectively. We conservatively set $\chi = 0.01$. In our examination of the GD-1 stream in Section 9.4 below, we deduce $\zeta = 0.27$ for that system; however, to be conservative, we set $\zeta = 0.1$ for all streams in our MCMC analysis. The

stream model

$$P_{\text{stream}}^{\text{fit}} = P_{\text{width}} \times P_{\text{color}} \times P_{\mu} \times P_{\infty} \times P_{\nu} \quad (11)$$

is exactly the same stream probability criterion as used in STREAMFINDER (Equation (2)) but ignores the luminosity
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Kinematic profiles and color–magnitude relations of stellar streams. The left column of panels shows the heliocentric line-of-sight velocity as a function of angle $\phi_1$ along the streams, while the middle column displays proper motion $\mu_{\phi_1}$ along the $\phi_1$ direction. The color–magnitude properties of the stellar populations can be seen in the right column of panels.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Kinematic profiles and color–magnitude relations of stellar streams. The left column of panels shows the heliocentric line-of-sight velocity as a function of angle $\phi_1$ along the streams, while the middle column displays proper motion $\mu_{\phi_1}$ along the $\phi_1$ direction. The color–magnitude properties of the stellar populations can be seen in the right column of panels.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 8. Continuation of Figure 7.</figcaption></figure>

function information (as we are assuming at this point that the
stars have been appropriately selected). Also, in contrast to
STREAMFINDER, the orbital track is corrected by the function

$\Delta_{\Theta}(\phi_1)$ introduced above, so as to better approximate the
corresponding stream (these functions will be discussed in
more detail in Section 7 below). The cocoon model is assumed
<figure><img src="image_11.png" /><figcaption>Figure 8. Continuation of Figure 7.</figcaption></figure>

function information (as we are assuming at this point that the
stars have been appropriately selected). Also, in contrast to
STREAMFINDER, the orbital track is corrected by the function

$\Delta_{\Theta}(\phi_1)$ introduced above, so as to better approximate the
corresponding stream (these functions will be discussed in
more detail in Section 7 below). The cocoon model is assumed
<figure><img src="image_12.png" /><figcaption>Figure 9. Continuation of Figure 7.</figcaption></figure>

to be identical to the stream model, except that we widen the model dispersions by a factor of 5 in line-of-sight velocity (approximately the ratio inferred in Section 9.4 for GD-1), proper motion, distance, and width on the sky (i.e., equivalent to an intrinsic 25 km s⁻¹ velocity dispersion and an intrinsic 250 pc stream width). In the present context, the contamination is clearly a very complex multidimensional function and further compounding the complexity of our model does not seem like a good idea. So instead, we adopted a pragmatic contamination model, inspired by the *good-and-bad data model* of Sivia & Skilling (2006; see their Section 8.3.2). According to their model, the *good* fraction has Gaussian residuals $R = (F - D)/\sigma$ between the data *D* and the model prediction *F* given the expected dispersion *σ* (due to the model spread and the

measurement uncertainties), while for the *bad* fraction, these residuals are simply scaled down by some factor τ (Sivia & Skilling 2006 use the symbol γ). In this way, we choose the contamination model to be exactly the same as the stream model (not the cocoon), but with τ = 10, effectively inflating the σ of the stream model by that factor (for bright stars with small measurement uncertainties this corresponds to a velocity dispersion of 50 km s⁻¹ and a spatial dispersion of 0.5 kpc). This renders the objective function (Equation (10)) much less sensitive to the presence of outliers, even though the assumed contamination fraction (1%) is quite small.

The ln $L_{\text{ancillary}}$ term in Equation (10) is used to include constraints from ancillary data sets, all in the form of simple $\chi^2$-like terms. So, for instance, when we force the rotation
<figure><img src="image_12.png" /><figcaption>Figure 9. Continuation of Figure 7.</figcaption></figure>

to be identical to the stream model, except that we widen the model dispersions by a factor of 5 in line-of-sight velocity (approximately the ratio inferred in Section 9.4 for GD-1), proper motion, distance, and width on the sky (i.e., equivalent to an intrinsic 25 km s⁻¹ velocity dispersion and an intrinsic 250 pc stream width). In the present context, the contamination is clearly a very complex multidimensional function and further compounding the complexity of our model does not seem like a good idea. So instead, we adopted a pragmatic contamination model, inspired by the *good-and-bad data model* of Sivia & Skilling (2006; see their Section 8.3.2). According to their model, the *good* fraction has Gaussian residuals $R = (F - D)/\sigma$ between the data *D* and the model prediction *F* given the expected dispersion *σ* (due to the model spread and the

measurement uncertainties), while for the *bad* fraction, these residuals are simply scaled down by some factor τ (Sivia & Skilling 2006 use the symbol γ). In this way, we choose the contamination model to be exactly the same as the stream model (not the cocoon), but with τ = 10, effectively inflating the σ of the stream model by that factor (for bright stars with small measurement uncertainties this corresponds to a velocity dispersion of 50 km s⁻¹ and a spatial dispersion of 0.5 kpc). This renders the objective function (Equation (10)) much less sensitive to the presence of outliers, even though the assumed contamination fraction (1%) is quite small.

The ln $L_{\text{ancillary}}$ term in Equation (10) is used to include constraints from ancillary data sets, all in the form of simple $\chi^2$-like terms. So, for instance, when we force the rotation
curve of the proposed potential models $v_c^{\text{model}}$ to be consistent with Galactic rotation curve measurements $v_c^{\text{Eilers}}$ and associated uncertainties $\delta v_c^{\text{Eilers}}$ (Eilers et al. 2019), we include

$$ -\sum_{\text{data}} \frac{1}{2} \left( \frac{v_c^{\text{Eilers}} - v_c^{\text{model}}}{\delta v_c^{\text{Eilers}}} \right)^2 \quad (12) $$

in $\ln L_{\text{ancillary}}$. We limit the Eilers et al. (2019) rotation curve data to $R < 15$ kpc, as the outer stellar disk is both substantially warped and flared beyond that radius (see, e.g., Momany et al. 2006).

Further constraints on the radial mass distribution are available thanks to studies of the terminal velocity of H I in the inner Galaxy. The high-resolution profiles of McClure-Griffiths & Dickey (2007, 2016) are conveniently measured with respect to the local standard of rest, so we can directly compare them to the predictions of the mass models

$$ v_{\text{term}}^{\text{model}}(\ell) = \text{sgn}(\sin(\ell))v_c(R|\sin(\ell)|) - v_c(R_{\odot})\sin(\ell). \quad (13) $$

The H I terminal velocity curves as a function of angle show substantially localized bumps on scales of a few degrees (see, e.g., Figure 1 of McClure-Griffiths & Dickey 2016), implying correlated motion that our large-scale mass model cannot possibly reproduce. To mitigate against this, we decided to average the profiles over 2° intervals. Also, following McMillan (2017), we retain only those measurements with $|\sin(\ell)| > 0.5$ to avoid regions that are heavily affected by the presence of the Galactic bar. We assume an uncertainty of 6.2 km s⁻¹ of the terminal velocity measurements, as derived from the scatter between the measurements and the linear model fit by McClure-Griffiths & Dickey (2016).

We similarly add the datum of $71 \pm 6 M_{\odot} \text{ pc}^{-2}$ for the vertical force at $z = 1.1$ kpc above the Galactic midplane measured by Kuijken & Gilmore (1991). We also include the stellar density measurements toward the north Galactic cap of Ibata et al. (2017b) decomposed into thin and thick disk components according to photometric metallicity; the sample at Galactic latitude $b > 70°$ is used (shown in their Figure 12(f)) truncated to $z < 5$ kpc as their thin disk profile is noisy at larger $z$. As indicated in Ibata et al. (2017b), the stellar number density profiles are transformed into relative mass density profiles with simulations that use PARSEC isochrone models and take into account the selection function of their survey. They adopt an average age of 5 and 10 Gyr for the thin and thick disks, respectively, and note that the difference in the computed density corrections between these two ages is only 20%, implying that their measured relative vertical density profiles are not very sensitive to those reasonable assumptions on population age. The Ibata et al. (2017b) study chose not to constrain the normalization of the density profiles, so in the present work the thin and thick disk components are fit to those vertical profiles data while allowing for a free multiplicative constant.

Neither the Eilers et al. (2019) rotation curve measurements nor the long streams we have selected probe the inner regions of the Milky Way. We, therefore, chose for the present work to simply adopt the Galactic bulge component of the model previously fit by McMillan (2017). The truncation radius of the halo component is fixed at $r_t = 1000$ kpc, far beyond the expected virial radius of the Milky Way.

In our explorations we realized that the data set provides constraints on the peculiar velocity of the Sun that seem plausible given published studies, but that these constraints are somewhat degenerate with the freedom we allow to models of the potential. We defer an investigation into this issue to a future contribution, and therefore, fix the peculiar velocity of the Sun to the values found by Schönrich et al. (2010). Likewise, we also fix the distance to the Galactic Center to $R_{\odot} = 8.178$ kpc (GRAVITY Collaboration et al. 2019), and ignore the small uncertainties on this value.

We calculate the $M_{200}$ mass of the models, defined as the mass out to a radius of $r_{200}$ where the average interior density of the model is 200 times the critical density $\rho_c = 3H^2/8\pi G$ of the universe, and we assume a Hubble parameter of $H = 70$ km s⁻¹ Mpc⁻¹. Clearly, the calculated $M_{200}$ value involves a long extrapolation from the region constrained by the chosen sample of 29 streams. Nevertheless, we consider $M_{200}$ to be a useful constraint on the models, and require the solutions to be in the range of [0, 2] × 10¹² $M_{\odot}$. The upper bound on $M_{200}$ of $2 \times 10^{12} M_{\odot}$ corresponds to half of the Local Group mass of $4 \times 10^{12} M_{\odot}$ (95% upper limit) estimated by Peñarrubia et al. (2014). The other bounds are that the masses of all components, as well as the halo density flattening parameter $q_{m,h}$, should be positive.<fnref n="17" /> The thin and thick disk scale lengths are required to be in the range of [1, 10] kpc, so we impose that the thin disk scale height should be in the range of [0.05, 0.5] kpc and the thick disk scale height should be in the range of [0.05, 5] kpc. Finally, we require that the thick disk scale height should be larger than that of the thin disk.

In our first modeling efforts we considered what we thought was the simplest useful Galaxy model, containing only a fixed bulge, a halo, and a single disk component. However, we noticed that if we left out the Ibata et al. (2017b) vertical density constraint, the resulting best solutions required a disk component with an unrealistically large scale height (~1 kpc). This result is perhaps not surprising in retrospect since the selected streams mostly reside in halo regions at high extraplanar distance, where the effect of the thin disk is less pronounced. So, evidently, a minimal Galaxy model for our purposes also entails including a thick disk component. Because our data provide poor constraints on the potential at low Galactic latitude, we also adopt the two gas disk models proposed by McMillan (2017), which include an atomic gas disk of scale height of 85 pc and mass of $1.1 \times 10^{10} M_{\odot}$, as well as a molecular gas disk of scale height 45 pc and mass of $1.1 \times 10^{9} M_{\odot}$.

Before showing our results with this model, we next present a refined—but computationally much more costly—method that uses test particles to follow tidal dissolution in a moving Local Group potential. This model is also used to calculate the correction functions $\Delta_{\Theta}$ employed in the orbit fitting described above.

## 6. Stream Fitting with Test Particles in the Local Group

Contrary to the assumptions underlying the models presented in the previous section, the Milky Way is, of course, not an isolated galaxy, but a member of the Local Group, the other major components of this grouping being the Andromeda and

<footnote n="17">¹⁷ These bounds are implemented in the MCMC procedure simply by subtracting a very large number from $\ln L$ if the limits are crossed.</footnote>
curve of the proposed potential models $v_c^{\text{model}}$ to be consistent with Galactic rotation curve measurements $v_c^{\text{Eilers}}$ and associated uncertainties $\delta v_c^{\text{Eilers}}$ (Eilers et al. 2019), we include

$$-\sum_{\text{data}} \frac{1}{2} \left( \frac{v_c^{\text{Eilers}} - v_c^{\text{model}}}{\delta v_c^{\text{Eilers}}} \right)^2 \quad (12)$$

in $\ln L_{\text{ancillary}}$. We limit the Eilers et al. (2019) rotation curve data to $R < 15$ kpc, as the outer stellar disk is both substantially warped and flared beyond that radius (see, e.g., Momany et al. 2006).

Further constraints on the radial mass distribution are available thanks to studies of the terminal velocity of H I in the inner Galaxy. The high-resolution profiles of McClure-Griffiths & Dickey (2007, 2016) are conveniently measured with respect to the local standard of rest, so we can directly compare them to the predictions of the mass models

$$v_{\text{term}}^{\text{model}}(\ell) = \text{sgn}(\sin(\ell))v_c(R|\sin(\ell)|) - v_c(R_{\odot})\sin(\ell). \quad (13)$$

The H I terminal velocity curves as a function of angle show substantially localized bumps on scales of a few degrees (see, e.g., Figure 1 of McClure-Griffiths & Dickey 2016), implying correlated motion that our large-scale mass model cannot possibly reproduce. To mitigate against this, we decided to average the profiles over 2° intervals. Also, following McMillan (2017), we retain only those measurements with $|\sin(\ell)| > 0.5$ to avoid regions that are heavily affected by the presence of the Galactic bar. We assume an uncertainty of 6.2 km s⁻¹ of the terminal velocity measurements, as derived from the scatter between the measurements and the linear model fit by McClure-Griffiths & Dickey (2016).

We similarly add the datum of $71 \pm 6 M_{\odot} \text{ pc}^{-2}$ for the vertical force at $z = 1.1$ kpc above the Galactic midplane measured by Kuijken & Gilmore (1991). We also include the stellar density measurements toward the north Galactic cap of Ibata et al. (2017b) decomposed into thin and thick disk components according to photometric metallicity; the sample at Galactic latitude $b > 70°$ is used (shown in their Figure 12(f)) truncated to $z < 5$ kpc as their thin disk profile is noisy at larger $z$. As indicated in Ibata et al. (2017b), the stellar number density profiles are transformed into relative mass density profiles with simulations that use PARSEC isochrone models and take into account the selection function of their survey. They adopt an average age of 5 and 10 Gyr for the thin and thick disks, respectively, and note that the difference in the computed density corrections between these two ages is only 20%, implying that their measured relative vertical density profiles are not very sensitive to those reasonable assumptions on population age. The Ibata et al. (2017b) study chose not to constrain the normalization of the density profiles, so in the present work the thin and thick disk components are fit to those vertical profiles data while allowing for a free multiplicative constant.

Neither the Eilers et al. (2019) rotation curve measurements nor the long streams we have selected probe the inner regions of the Milky Way. We, therefore, chose for the present work to simply adopt the Galactic bulge component of the model previously fit by McMillan (2017). The truncation radius of the halo component is fixed at $r_t = 1000$ kpc, far beyond the expected virial radius of the Milky Way.

In our explorations we realized that the data set provides constraints on the peculiar velocity of the Sun that seem plausible given published studies, but that these constraints are somewhat degenerate with the freedom we allow to models of the potential. We defer an investigation into this issue to a future contribution, and therefore, fix the peculiar velocity of the Sun to the values found by Schönrich et al. (2010). Likewise, we also fix the distance to the Galactic Center to $R_{\odot} = 8.178$ kpc (GRAVITY Collaboration et al. 2019), and ignore the small uncertainties on this value.

We calculate the $M_{200}$ mass of the models, defined as the mass out to a radius of $r_{200}$ where the average interior density of the model is 200 times the critical density $\rho_c = 3H^2/8\pi G$ of the universe, and we assume a Hubble parameter of $H = 70$ km s⁻¹ Mpc⁻¹. Clearly, the calculated $M_{200}$ value involves a long extrapolation from the region constrained by the chosen sample of 29 streams. Nevertheless, we consider $M_{200}$ to be a useful constraint on the models, and require the solutions to be in the range of $[0, 2] \times 10^{12} M_{\odot}$. The upper bound on $M_{200}$ of $2 \times 10^{12} M_{\odot}$ corresponds to half of the Local Group mass of $4 \times 10^{12} M_{\odot}$ (95% upper limit) estimated by Peñarrubia et al. (2014). The other bounds are that the masses of all components, as well as the halo density flattening parameter $q_{m,h}$, should be positive.<fnref n="17" /> The thin and thick disk scale lengths are required to be in the range of [1, 10] kpc, so we impose that the thin disk scale height should be in the range of [0.05, 0.5] kpc and the thick disk scale height should be in the range of [0.05, 5] kpc. Finally, we require that the thick disk scale height should be larger than that of the thin disk.

In our first modeling efforts we considered what we thought was the simplest useful Galaxy model, containing only a fixed bulge, a halo, and a single disk component. However, we noticed that if we left out the Ibata et al. (2017b) vertical density constraint, the resulting best solutions required a disk component with an unrealistically large scale height (~1 kpc). This result is perhaps not surprising in retrospect since the selected streams mostly reside in halo regions at high extraplanar distance, where the effect of the thin disk is less pronounced. So, evidently, a minimal Galaxy model for our purposes also entails including a thick disk component. Because our data provide poor constraints on the potential at low Galactic latitude, we also adopt the two gas disk models proposed by McMillan (2017), which include an atomic gas disk of scale height of 85 pc and mass of $1.1 \times 10^{10} M_{\odot}$, as well as a molecular gas disk of scale height 45 pc and mass of $1.1 \times 10^{9} M_{\odot}$.

Before showing our results with this model, we next present a refined—but computationally much more costly—method that uses test particles to follow tidal dissolution in a moving Local Group potential. This model is also used to calculate the correction functions $\Delta_{\Theta}$ employed in the orbit fitting described above.

## 6. Stream Fitting with Test Particles in the Local Group

Contrary to the assumptions underlying the models presented in the previous section, the Milky Way is, of course, not an isolated galaxy, but a member of the Local Group, the other major components of this grouping being the Andromeda and

<footnote n="17">¹⁷ These bounds are implemented in the MCMC procedure simply by subtracting a very large number from $\ln L$ if the limits are crossed.</footnote>
Triangulum galaxies, approximately 800 kpc away, as well as the LMC, the Small Magellanic Cloud, and the Sagittarius dwarf galaxy that reside within the halo of our Galaxy. Predicting the paths of stellar streams in the Local Group requires some consideration of the combined gravitational effect of these six bodies, and perhaps additional mass concentrations depending on the particular path of the stream progenitor.

Yet, running full *N*-body simulations of the Local Group with a resolution sufficient to fit the observed stellar streams is far beyond current computational abilities, so we need to make some simplifying steps to tackle the problem. Here, we will assume that the major galaxies can each be modeled as moving potentials with a mass distribution that does not vary in time when viewed from their center of mass. Since this assumption becomes worse further back in time, we will only consider integrating backward for at most 5 Gyr.

We integrate the path of the Local Group members using a symplectic leapfrog scheme with direct summation of the forces from all other bodies. The major galaxies (the Milky Way, Andromeda, and the LMC) are assumed to be truncated at their (static) virial radius. This allows us to treat the bodies as point masses when considering their effect on the other Local Group galaxies that are beyond the virial radius. When including M31 or M33 (whose influence we neglect in the present work), we also take into account a dark energy term, using an additional acceleration term $\ddot{x} = H_0^2 \Omega_\Lambda x$, where $x$ in this case is the radial vector from the center of mass of the Local Group.

The dynamics of the LMC and Sagittarius galaxies are particularly difficult because we expect their orbits to have been substantially affected by dynamical friction. Here, we use the Chandrasekhar dynamical friction model to integrate the massive bodies backward in time, and adopt the formula of Correa Magnus & Vasiliev (2022) to approximate the Coulomb logarithm. We have confirmed with several spot-tests that the Correa Magnus & Vasiliev (2022) method predicts well the behavior in *N*-body simulations of the spiraling-in path of the LMC when the Milky Way is represented by a live Navarro–Frenk–White (NFW) halo.

The orbits of the stream progenitors are integrated backward in time for a time $T$ in this dynamic Local Group potential, starting from a trial present-day phase-space coordinate. At the initial time $-T$, we instantiate a set of $N = 10,000$ massless test particles that follow a spherical Plummer model (parameters: mass, scale radius).<fnref n="18" /> We then integrate this system forward in time until the present day. However, instead of following the full *N*-body evolution of the progenitor, we approximate tidal dissolution as follows. We use the analytic potential of the progenitor to approximate the evolving system, so that the forces on the test particles are the sum of the forces from the Local Group galaxies, plus the force due to the progenitor. If a particle moves beyond 10 scale radii from the progenitor, it is considered to be lost, and we decrease the mass of the progenitor by $1/N$. In contrast, we assume that the Plummer scale radius $a$ does not change during the tidal disruption. The advantage of this scheme is that it has saved us from having to follow the internal *N*-body dynamics of the progenitor, and it can also be parallelized trivially since the test

particles are independent of each other ($O(N)$ computational complexity). Our method is inspired by those of Mastrobuono-Battisti et al. (2012) and Ferrone et al. (2023), but we implement mass loss in the progenitor, as well as massive perturbers. Mastrobuono-Battisti et al. (2012) demonstrated that test-particle integration in nonvarying host globular clusters can reproduce very similar stream extent and morphology as full *N*-body simulations, but for the majority of the STREAMFINDER streams, we need to be able to account for the full dissolution of the progenitor.

The validity of our simplifying approximation of fixing the Plummer scale radius can be judged from Figure 10, where we show a separate *N*-body experiment of the evolution of a plausible GD-1 stream progenitor. For this test, we use the Galactic potential model #1 of Dehnen & Binney (1998) and integrate therein a Plummer model for 2 Gyr. Over this period of time, where the system loses almost all its mass (red profile), the half-mass-radius $r_h = (2^{2/3} - 1)^{-1/2}a$ (blue profile) remains close to constant. An equivalent result is discussed in Errani et al. (2024) for collision-less systems that follow an exponential density profile, suggesting that our simplifying assumption of $a = \text{const}$ is not unique to our choice of progenitor density profile. Note, however, that our idealized setup does not capture the effects of potential internal collisional dynamics on the progenitor structure.

We use exactly the same optimization machinery and the same objective function for the stream fits as for the orbit fits in Section 5 (except, of course, that the correction functions $\Delta\Theta$ are ignored). Unfortunately, this procedure of integrating test particles is substantially more computationally expensive to optimize: most importantly, we have to follow $N$ times more orbits, and there are two more parameters per stream (the progenitor mass and scale radius of the Plummer model), and the global potential requires seven additional parameters (mass and initial phase-space position) for each Local Group galaxy that we aim to include the potential of. While it would be possible to reduce $N$, there needs to be a sufficient number of tidally disrupted particles to model a stream, and our experiments showed that below $N \approx 10^4$, the stream model could become too noisy to allow a reliable comparison to data.

In the present contribution, we will only attempt a limited exploration of the parameter space of this complex model, and defer a full stream-fitting analysis with this machinery to a future contribution.

## 7. Orbit-Stream Correction Functions

The test-particle model introduced in the previous section is immediately useful to us as a means to calculate the $\Delta\Theta(\phi_1)$ correction functions needed for the corrected orbit fitting described in Section 5. We started by running the test-particle model optimization procedure with the Milky Way mass parameters fixed at the values found by McMillan (2017), and with all other Local Group galaxies set to zero mass. For the five systems with known globular cluster progenitors (Gjöll/NGC 3201, Fjörm/M68, NGC 5466, Palomar 5, and M5, i.e., streams #34, #47, #57, #61, and #62) in our sample of 29 streams, we fixed the mass and half-mass parameters at the values measured by Baumgardt & Hilker (2018), and also fixed the line-of-sight velocities at the values given by Vasiliev (2019); the distances and proper motions were allowed to vary within the uncertainties set by the measurements of Baumgardt & Vasiliev (2021) and Vasiliev (2019), respectively. For GD-1,

<footnote n="18">¹⁸ For streams that have dwarf galaxy progenitors an exponentially truncated NFW model would be more appropriate (Errani & Navarro 2021), and would be very simple to implement into this machinery.</footnote>
Triangulum galaxies, approximately 800 kpc away, as well as the LMC, the Small Magellanic Cloud, and the Sagittarius dwarf galaxy that reside within the halo of our Galaxy. Predicting the paths of stellar streams in the Local Group requires some consideration of the combined gravitational effect of these six bodies, and perhaps additional mass concentrations depending on the particular path of the stream progenitor.

Yet, running full *N*-body simulations of the Local Group with a resolution sufficient to fit the observed stellar streams is far beyond current computational abilities, so we need to make some simplifying steps to tackle the problem. Here, we will assume that the major galaxies can each be modeled as moving potentials with a mass distribution that does not vary in time when viewed from their center of mass. Since this assumption becomes worse further back in time, we will only consider integrating backward for at most 5 Gyr.

We integrate the path of the Local Group members using a symplectic leapfrog scheme with direct summation of the forces from all other bodies. The major galaxies (the Milky Way, Andromeda, and the LMC) are assumed to be truncated at their (static) virial radius. This allows us to treat the bodies as point masses when considering their effect on the other Local Group galaxies that are beyond the virial radius. When including M31 or M33 (whose influence we neglect in the present work), we also take into account a dark energy term, using an additional acceleration term $\ddot{x} = H_0^2 \Omega_\Lambda x$, where $x$ in this case is the radial vector from the center of mass of the Local Group.

The dynamics of the LMC and Sagittarius galaxies are particularly difficult because we expect their orbits to have been substantially affected by dynamical friction. Here, we use the Chandrasekhar dynamical friction model to integrate the massive bodies backward in time, and adopt the formula of Correa Magnus & Vasiliev (2022) to approximate the Coulomb logarithm. We have confirmed with several spot-tests that the Correa Magnus & Vasiliev (2022) method predicts well the behavior in *N*-body simulations of the spiraling-in path of the LMC when the Milky Way is represented by a live Navarro–Frenk–White (NFW) halo.

The orbits of the stream progenitors are integrated backward in time for a time $T$ in this dynamic Local Group potential, starting from a trial present-day phase-space coordinate. At the initial time $-T$, we instantiate a set of $N = 10,000$ massless test particles that follow a spherical Plummer model (parameters: mass, scale radius).<fnref n="18" /> We then integrate this system forward in time until the present day. However, instead of following the full *N*-body evolution of the progenitor, we approximate tidal dissolution as follows. We use the analytic potential of the progenitor to approximate the evolving system, so that the forces on the test particles are the sum of the forces from the Local Group galaxies, plus the force due to the progenitor. If a particle moves beyond 10 scale radii from the progenitor, it is considered to be lost, and we decrease the mass of the progenitor by $1/N$. In contrast, we assume that the Plummer scale radius $a$ does not change during the tidal disruption. The advantage of this scheme is that it has saved us from having to follow the internal *N*-body dynamics of the progenitor, and it can also be parallelized trivially since the test

particles are independent of each other ($O(N)$ computational complexity). Our method is inspired by those of Mastrobuono-Battisti et al. (2012) and Ferrone et al. (2023), but we implement mass loss in the progenitor, as well as massive perturbers. Mastrobuono-Battisti et al. (2012) demonstrated that test-particle integration in nonvarying host globular clusters can reproduce very similar stream extent and morphology as full *N*-body simulations, but for the majority of the STREAMFINDER streams, we need to be able to account for the full dissolution of the progenitor.

The validity of our simplifying approximation of fixing the Plummer scale radius can be judged from Figure 10, where we show a separate *N*-body experiment of the evolution of a plausible GD-1 stream progenitor. For this test, we use the Galactic potential model #1 of Dehnen & Binney (1998) and integrate therein a Plummer model for 2 Gyr. Over this period of time, where the system loses almost all its mass (red profile), the half-mass-radius $r_h = (2^{2/3} - 1)^{-1/2}a$ (blue profile) remains close to constant. An equivalent result is discussed in Errani et al. (2024) for collision-less systems that follow an exponential density profile, suggesting that our simplifying assumption of $a = \text{const}$ is not unique to our choice of progenitor density profile. Note, however, that our idealized setup does not capture the effects of potential internal collisional dynamics on the progenitor structure.

We use exactly the same optimization machinery and the same objective function for the stream fits as for the orbit fits in Section 5 (except, of course, that the correction functions $\Delta\Theta$ are ignored). Unfortunately, this procedure of integrating test particles is substantially more computationally expensive to optimize: most importantly, we have to follow $N$ times more orbits, and there are two more parameters per stream (the progenitor mass and scale radius of the Plummer model), and the global potential requires seven additional parameters (mass and initial phase-space position) for each Local Group galaxy that we aim to include the potential of. While it would be possible to reduce $N$, there needs to be a sufficient number of tidally disrupted particles to model a stream, and our experiments showed that below $N \approx 10^4$, the stream model could become too noisy to allow a reliable comparison to data.

In the present contribution, we will only attempt a limited exploration of the parameter space of this complex model, and defer a full stream-fitting analysis with this machinery to a future contribution.

## 7. Orbit-Stream Correction Functions

The test-particle model introduced in the previous section is immediately useful to us as a means to calculate the $\Delta\Theta(\phi_1)$ correction functions needed for the corrected orbit fitting described in Section 5. We started by running the test-particle model optimization procedure with the Milky Way mass parameters fixed at the values found by McMillan (2017), and with all other Local Group galaxies set to zero mass. For the five systems with known globular cluster progenitors (Gjöll/NGC 3201, Fjörm/M68, NGC 5466, Palomar 5, and M5, i.e., streams #34, #47, #57, #61, and #62) in our sample of 29 streams, we fixed the mass and half-mass parameters at the values measured by Baumgardt & Hilker (2018), and also fixed the line-of-sight velocities at the values given by Vasiliev (2019); the distances and proper motions were allowed to vary within the uncertainties set by the measurements of Baumgardt & Vasiliev (2021) and Vasiliev (2019), respectively. For GD-1,

<footnote n="18">¹⁸ For streams that have dwarf galaxy progenitors an exponentially truncated NFW model would be more appropriate (Errani & Navarro 2021), and would be very simple to implement into this machinery.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Tidal dissolution of an *N*-body GD-1 progenitor model, showing the approximate constancy of the half-mass–radius as the system loses almost all its mass. For this experiment, a Plummer sphere of $10^3$ particles with a scale radius *a* = 50 pc and a mass of $3 \times 10^4 M_\odot$ was evolved using the falcon (Dehnen 2002) *N*-body integrator. A softening length of 1 pc was adopted.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Difference between the line-of-sight velocity correction functions $\Delta_v(\phi_1) = v_{\text{stream}} - v_{\text{orbit}}$ corresponding to the most likely Milky Way model and that resulting from imposing $\beta_h = 3$. It can be seen that the correction functions differ by only ≈1 km s⁻¹. The points are plotted at the $\phi_1$ locations of the stars in the 29 selected streams, and we again use the color map of Figure 6.</figcaption></figure>

we set the progenitor position to $\alpha = 157^\circ 74$, as will be justified below in Section 8. For the remaining 23 streams, we unfortunately have no constraint on the position of the progenitor, so it is set at the estimated midpoint of the stream ($\phi_1 = 0^\circ$). The masses of the structures are set to the stellar mass estimates listed in Table 3, but we conservatively set a floor of $2 \times 10^4 M_\odot$ to these values (approximately the mass of our estimate of the Palomar 5 progenitor). For these streams, we also assume a progenitor Plummer scale length of 30 pc. The correction function $\Delta_\Theta(\phi_1)$ is then calculated independently for each observable $\Theta$ (i.e., sky position, distance

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Comparison of the GD-1 data (stream #53, red points) with a best-fit model (black points) in a potential model that is constrained to have fixed halo power-law parameters of $\gamma_h = 1$ and $\beta_h = 3$, fixed scale heights of 0.3 and 0.9 kpc, respectively, for the thin and thick disk components, together with a fixed thick disk scale length of 3.6 kpc (all these parameters are also fixed in the McMillan 2017 model). The model fits the halo central density, halo scale parameter, the disk scale length, and the masses and relative fractions of the thin and thick disks. Although all GD-1 data are displayed, only those stars with measured radial velocity (blue points) were used in the fit. The GD-1 model has an initial mass of $3 \times 10^4 M_\odot$, a Plummer scale length of 50 pc, and is integrated over 2.5 Gyr. To aid visual comparison, the sky positions (e) include a version of the model shifted by −5° (gray) as well as the data shifted by +5° (magenta) in $\phi_2$. The gap in the data visible in all panels at $\phi_1 \approx 68^\circ$ is a consequence of our choice in running STREAMFINDER to excise Gaia data around dense sources (in this case, M67).</figcaption></figure>

modulus, line-of-sight velocity, and proper motions) as a fourth-order polynomial fit to the best-fit stream minus the corresponding progenitor orbit. This $\Delta_\Theta(\phi_1)$ is subsequently used in the corrected orbit-fitting algorithm to find an improved Milky Way mass model. We repeated this procedure a further two times, rederiving $\Delta_\Theta(\phi_1)$ for the best mass model, and rerunning the MCMC search of Section 5. We stopped at the third iteration as the derived Milky Way parameters were statistically identical to those in the second iteration.

The procedure we have devised makes several simplifications to speed up the calculations, including the use of a single correction function $\Delta_\Theta(\phi_1)$ for each stream during an MCMC run. Ideally, a new correction function should be calculated for each stream for every set of Milky Way and stream parameters
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Tidal dissolution of an *N*-body GD-1 progenitor model, showing the approximate constancy of the half-mass–radius as the system loses almost all its mass. For this experiment, a Plummer sphere of $10^3$ particles with a scale radius *a* = 50 pc and a mass of $3 \times 10^4 M_\odot$ was evolved using the falcon (Dehnen 2002) *N*-body integrator. A softening length of 1 pc was adopted.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Difference between the line-of-sight velocity correction functions $\Delta_v(\phi_1) = v_{\text{stream}} - v_{\text{orbit}}$ corresponding to the most likely Milky Way model and that resulting from imposing $\beta_h = 3$. It can be seen that the correction functions differ by only ≈1 km s⁻¹. The points are plotted at the $\phi_1$ locations of the stars in the 29 selected streams, and we again use the color map of Figure 6.</figcaption></figure>

we set the progenitor position to $\alpha = 157^\circ 74$, as will be justified below in Section 8. For the remaining 23 streams, we unfortunately have no constraint on the position of the progenitor, so it is set at the estimated midpoint of the stream ($\phi_1 = 0^\circ$). The masses of the structures are set to the stellar mass estimates listed in Table 3, but we conservatively set a floor of $2 \times 10^4 M_\odot$ to these values (approximately the mass of our estimate of the Palomar 5 progenitor). For these streams, we also assume a progenitor Plummer scale length of 30 pc. The correction function $\Delta_\Theta(\phi_1)$ is then calculated independently for each observable $\Theta$ (i.e., sky position, distance

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Comparison of the GD-1 data (stream #53, red points) with a best-fit model (black points) in a potential model that is constrained to have fixed halo power-law parameters of $\gamma_h = 1$ and $\beta_h = 3$, fixed scale heights of 0.3 and 0.9 kpc, respectively, for the thin and thick disk components, together with a fixed thick disk scale length of 3.6 kpc (all these parameters are also fixed in the McMillan 2017 model). The model fits the halo central density, halo scale parameter, the disk scale length, and the masses and relative fractions of the thin and thick disks. Although all GD-1 data are displayed, only those stars with measured radial velocity (blue points) were used in the fit. The GD-1 model has an initial mass of $3 \times 10^4 M_\odot$, a Plummer scale length of 50 pc, and is integrated over 2.5 Gyr. To aid visual comparison, the sky positions (e) include a version of the model shifted by −5° (gray) as well as the data shifted by +5° (magenta) in $\phi_2$. The gap in the data visible in all panels at $\phi_1 \approx 68^\circ$ is a consequence of our choice in running STREAMFINDER to excise Gaia data around dense sources (in this case, M67).</figcaption></figure>

modulus, line-of-sight velocity, and proper motions) as a fourth-order polynomial fit to the best-fit stream minus the corresponding progenitor orbit. This $\Delta_\Theta(\phi_1)$ is subsequently used in the corrected orbit-fitting algorithm to find an improved Milky Way mass model. We repeated this procedure a further two times, rederiving $\Delta_\Theta(\phi_1)$ for the best mass model, and rerunning the MCMC search of Section 5. We stopped at the third iteration as the derived Milky Way parameters were statistically identical to those in the second iteration.

The procedure we have devised makes several simplifications to speed up the calculations, including the use of a single correction function $\Delta_\Theta(\phi_1)$ for each stream during an MCMC run. Ideally, a new correction function should be calculated for each stream for every set of Milky Way and stream parameters
<table><caption>**Table 1**
The First 10 Rows of the STREAMFINDER Catalog of 24,540 Stars Detected in Gaia DR3</caption><thead><tr><th>DR3 ID<br>(1)</th><th>α<br>(deg)<br>(2)</th><th>δ<br>(deg)<br>(3)</th><th>ω<br>(mas)<br>(4)</th><th>μ<sub>α</sub><br>(mas yr<sup>-1</sup>)<br>(5)</th><th>μ<sub>δ</sub><br>(mas yr<sup>-1</sup>)<br>(6)</th><th>G<sub>0</sub><br>(mag)<br>(7)</th><th>(G<sub>BP</sub> - G<sub>RP</sub>)<sub>0</sub><br>(mag)<br>(8)</th><th>ν<sub>h</sub><br>(km s<sup>-1</sup>)<br>(9)</th><th>δν<sub>h</sub><br>(km s<sup>-1</sup>)<br>(10)</th><th>s<br>(11)</th><th>s<sub>ID</sub><br>(12)</th></tr></thead><tbody><tr><td>2746748683780680704</td><td>0.131771</td><td>8.022269</td><td>-0.523</td><td>-0.338</td><td>1.390</td><td>19.35</td><td>0.74</td><td>...</td><td>...</td><td>...</td><td>1</td></tr><tr><td>2746696319539384320</td><td>0.328548</td><td>7.886247</td><td>0.045</td><td>-0.683</td><td>1.425</td><td>19.48</td><td>0.59</td><td>...</td><td>...</td><td>...</td><td>1</td></tr><tr><td>2745893989584959488</td><td>0.714213</td><td>7.207192</td><td>0.090</td><td>-0.553</td><td>1.405</td><td>15.83</td><td>1.04</td><td>-115.03</td><td>0.35</td><td>13</td><td>1</td></tr><tr><td>2745895810651091456</td><td>0.714858</td><td>7.301568</td><td>-0.001</td><td>-0.439</td><td>1.322</td><td>18.75</td><td>0.83</td><td>...</td><td>...</td><td>...</td><td>1</td></tr><tr><td>2745605054248759296</td><td>0.992495</td><td>6.646490</td><td>-0.413</td><td>-0.760</td><td>1.801</td><td>19.04</td><td>0.68</td><td>...</td><td>...</td><td>...</td><td>1</td></tr><tr><td>2745591005412035584</td><td>1.135276</td><td>6.569216</td><td>-0.055</td><td>-0.251</td><td>1.656</td><td>18.00</td><td>0.89</td><td>...</td><td>...</td><td>...</td><td>1</td></tr><tr><td>2742045419713596416</td><td>1.655054</td><td>5.419648</td><td>0.186</td><td>-0.531</td><td>1.219</td><td>19.22</td><td>0.75</td><td>...</td><td>...</td><td>...</td><td>1</td></tr><tr><td>2742028759537330176</td><td>1.740425</td><td>5.219596</td><td>0.047</td><td>-0.572</td><td>1.522</td><td>19.59</td><td>0.71</td><td>...</td><td>...</td><td>...</td><td>1</td></tr><tr><td>2740941544398892032</td><td>2.436206</td><td>4.142149</td><td>0.545</td><td>-0.651</td><td>1.321</td><td>18.60</td><td>0.79</td><td>...</td><td>...</td><td>...</td><td>1</td></tr><tr><td>2740973155358226432</td><td>2.751825</td><td>4.308954</td><td>-0.239</td><td>-0.391</td><td>1.247</td><td>19.20</td><td>0.68</td><td>...</td><td>...</td><td>...</td><td>1</td></tr></tbody></table>

Note. Column (1) provides the Gaia DR3 identification of the star, columns (2)–(6) list the EDR3 equatorial coordinates α and δ, parallax ω and proper motions μ<sub>α</sub>(*cos(δ)), μ<sub>δ</sub>. The extinction-corrected magnitude G<sub>0</sub>, and color (G<sub>BP</sub> - G<sub>RP</sub>)<sub>0</sub> used in STREAMFINDER are listed in columns (7) and (8). Columns (9) and (10) list the best-measured heliocentric line-of-sight velocity, as derived from the corresponding source “s” in column (11). In column (11), the source identifications “s” are: 1 = APOGEE, 2 = GALAH, 3 = Gaia RVS, 4 = LAMOST, 5 = S5, 6 = SDSS, 7 = BOSS, 8 = ESPaDOnS (this work), 9 = AAOmega (from Ibata et al. 2017a), 10 = FLAMES (from Ibata et al. 2017a), 11 = UVES (from Odenkirchen et al. 2009), 12 = EFOSC (this work), 13 = UVES (this work), 14 = INT (this work), 15 = Yuan et al. (2022b), 16 = Li et al. (2021), 17 = Caldwell et al. (2020), 18 = T. S. Li et al. (2018), 19 = Koposov et al. (2018), 20 = G.-W. Li et al. (2018), 21 = Simon et al. (2020), 22 = Walker et al. (2015), 23 = VizieR, 24 = GES, 25 = DESI. Finally, column (12) provides a unique stream identification label. (This table is available in its entirety in machine-readable form.)

examined in the MCMC search.<fnref n="19" /> In Figure 11, we show the difference between the correction functions in the best-fit model, and in the best model when imposing a halo outer power-law parameter of β<sub>h</sub> = 3 (so as to be consistent with an NFW profile). In these two quite different Milky Way models, the correction functions are very similar, within ≈1 km s<sup>-1</sup> in line-of-sight velocity (which is the easiest observable to interpret when comparing streams, as there is no distance dependence).

## 8. Stream-fitting Results

We begin by attempting to fit a static axisymmetric model to the Galaxy, ignoring any stellar stream data; to this end, we use the objective function in Equation (10), with the first term on the right-hand side set to zero (i.e., we only fit to the ancillary data). For this first experiment, we fix the halo power-law parameters to γ<sub>h</sub> = 1 and β<sub>h</sub> = 3 so as to reproduce the expected behavior of dark matter halos (Navarro et al. 1997), set the halo matter flattening to q<sub>m,h</sub> = 1, set the thin and thick disk scale heights to 0.3 and 0.9 kpc respectively, and set the thick disk scale length to 3.6 kpc; these are all choices that were adopted by McMillan (2017). The main difference with respect to that study is that we are using the Eilers et al. (2019) rotation curve measurements. We find a thin disk scale length of h<sub>R,d</sub> = 2.22 ± 0.07(2.53 ± 0.14) kpc, a halo scale radius of r<sub>0,h</sub> = 12.4 ± 2.0(19.0 ± 4.9) kpc, a total stellar mass of M<sub>*</sub> = 4.63 ± 0.33(5.43 ± 0.57) × 10<sup>10</sup> M<sub>☉</sub>, M<sub>200</sub> = 0.86 ± 0.09(1.3 ± 0.3) × 10<sup>12</sup> M<sub>☉</sub>, and a mass interior to 50 kpc of M<sub>R&lt;50</sub> = 0.41 ± 0.02 × 10<sup>12</sup> M<sub>☉</sub>, where the values in brackets are those reported by McMillan (2017), when available. Figure 6 of McMillan (2017) shows their posterior PDF of M<sub>R&lt;50</sub> peaking at a considerably higher mass of M<sub>R&lt;50</sub> ≈ 0.51 × 10<sup>12</sup> M<sub>☉</sub>. Our aim with this test was to obtain a baseline for the

following stream analysis, and to explore how the improved data changes these inferences.

We now repeat this experiment, but add in the GD-1 stream, which is the structure for which we have the most numerous data (323 stars with line-of-sight velocity measurements and three RR Lyrae). This time we allow the halo mass flattening parameter q<sub>m,h</sub> to also be free. With this data set, we find h<sub>R,d</sub> = 2.16 ± 0.08 kpc, q<sub>m,h</sub> = 0.86 ± 0.05, r<sub>0,h</sub> = 10.7 ± 1.6 kpc, M<sub>*</sub> = 4.24 ± 0.35 × 10<sup>10</sup> M<sub>☉</sub>, M<sub>200</sub> = 0.77 ± 0.07 × 10<sup>12</sup> M<sub>☉</sub> and M<sub>R&lt;50</sub> = 0.39 ± 0.01 × 10<sup>12</sup> M<sub>☉</sub>, and the corresponding stream fit is shown in Figure 12. While the stream model captures the large-scale behavior of the kinematics of GD-1, it possesses some small, yet systematic, deviations from the proper motion trend in the data, especially in the range of 60° < φ<sub>1</sub> < 90°. The density distribution of the model particles (panel (e)) is also not a good match to the data, given that the model produces two high-density peaks at φ<sub>1</sub> ≈ 0° and φ<sub>1</sub> ≈ 70°, and a small remnant (at φ<sub>1</sub> = 29.7) still remains, but is not present in reality. However, the large peaks in the data at φ<sub>1</sub> ≈ 10° and ≈ 50°, which we ascribed in an earlier contribution (Ibata et al. 2020b) to epicyclic motion, can be seen to have counterparts in the model. Given this match, and also given that φ<sub>1</sub> = 32.4 corresponds to the stellar mass-weighted midpoint of the stream,<fnref n="20" /> we consider α = 157.74 to be currently the best estimate for the position of the center of the GD-1 progenitor. Henceforth, we fix the present-day position of GD-1's progenitor to that location.

We continue the procedure of rendering the model more realistic, adding in the five streams in our sample that have known globular cluster progenitors (Gjöll/NGC 3201, Fjörm/ M68, NGC 5466, Palomar 5, and M5). The phase-space parameters of these objects, and in particular their distances, are very well constrained, and in particular their distances, are very well constrained, which makes them especially powerful probes of the Galactic potential compared to a similarly sized sample of streams. In addition to the previously

<footnote n="19">¹⁹ Of course, this would be pointless in practice, as it would be more accurate and more computationally efficient to simply use the test-particle machinery developed in Section 6 and forgo the corrected orbit-fitting approach altogether.</footnote>

<footnote n="20">²⁰ We calculate this midpoint by assuming that the survey is complete to G<sub>0</sub> = 20 mag, and correct for missing stars given the distance model and the metallicity listed in Table 3. However, we did not correct for the masked region around M67.</footnote>
<table><caption>**Table 1**
The First 10 Rows of the STREAMFINDER Catalog of 24,540 Stars Detected in Gaia DR3</caption><thead><tr><th>DR3 ID<br>(1)</th><th>α<br>(deg)<br>(2)</th><th>δ<br>(deg)<br>(3)</th><th>ω<br>(mas)<br>(4)</th><th>μ<sub>α</sub><br>(mas yr<sup>-1</sup>)<br>(5)</th><th>μ<sub>δ</sub><br>(mas yr<sup>-1</sup>)<br>(6)</th><th>G<sub>0</sub><br>(mag)<br>(7)</th><th>(G<sub>BP</sub> - G<sub>RP</sub>)<sub>0</sub><br>(mag)<br>(8)</th><th>ν<sub>h</sub><br>(km s<sup>-1</sup>)<br>(9)</th><th>δν<sub>h</sub><br>(km s<sup>-1</sup>)<br>(10)</th><th>s<br>(11)</th><th>s<sub>ID</sub><br>(12)</th></tr></thead><tbody><tr><td>2746748683780680704</td><td>0.131771</td><td>8.022269</td><td>-0.523</td><td>-0.338</td><td>1.390</td><td>19.35</td><td>0.74</td><td>...</td><td>...</td><td>...</td><td>1</td></tr><tr><td>2746696319539384320</td><td>0.328548</td><td>7.886247</td><td>0.045</td><td>-0.683</td><td>1.425</td><td>19.48</td><td>0.59</td><td>...</td><td>...</td><td>...</td><td>1</td></tr><tr><td>2745893989584959488</td><td>0.714213</td><td>7.207192</td><td>0.090</td><td>-0.553</td><td>1.405</td><td>15.83</td><td>1.04</td><td>-115.03</td><td>0.35</td><td>13</td><td>1</td></tr><tr><td>2745895810651091456</td><td>0.714858</td><td>7.301568</td><td>-0.001</td><td>-0.439</td><td>1.322</td><td>18.75</td><td>0.83</td><td>...</td><td>...</td><td>...</td><td>1</td></tr><tr><td>2745605054248759296</td><td>0.992495</td><td>6.646490</td><td>-0.413</td><td>-0.760</td><td>1.801</td><td>19.04</td><td>0.68</td><td>...</td><td>...</td><td>...</td><td>1</td></tr><tr><td>2745591005412035584</td><td>1.135276</td><td>6.569216</td><td>-0.055</td><td>-0.251</td><td>1.656</td><td>18.00</td><td>0.89</td><td>...</td><td>...</td><td>...</td><td>1</td></tr><tr><td>2742045419713596416</td><td>1.655054</td><td>5.419648</td><td>0.186</td><td>-0.531</td><td>1.219</td><td>19.22</td><td>0.75</td><td>...</td><td>...</td><td>...</td><td>1</td></tr><tr><td>2742028759537330176</td><td>1.740425</td><td>5.219596</td><td>0.047</td><td>-0.572</td><td>1.522</td><td>19.59</td><td>0.71</td><td>...</td><td>...</td><td>...</td><td>1</td></tr><tr><td>2740941544398892032</td><td>2.436206</td><td>4.142149</td><td>0.545</td><td>-0.651</td><td>1.321</td><td>18.60</td><td>0.79</td><td>...</td><td>...</td><td>...</td><td>1</td></tr><tr><td>2740973155358226432</td><td>2.751825</td><td>4.308954</td><td>-0.239</td><td>-0.391</td><td>1.247</td><td>19.20</td><td>0.68</td><td>...</td><td>...</td><td>...</td><td>1</td></tr></tbody></table>

Note. Column (1) provides the Gaia DR3 identification of the star, columns (2)–(6) list the EDR3 equatorial coordinates α and δ, parallax ω and proper motions μ<sub>α</sub>(*cos(δ)), μ<sub>δ</sub>. The extinction-corrected magnitude G<sub>0</sub>, and color (G<sub>BP</sub> - G<sub>RP</sub>)<sub>0</sub> used in STREAMFINDER are listed in columns (7) and (8). Columns (9) and (10) list the best-measured heliocentric line-of-sight velocity, as derived from the corresponding source “s” in column (11). In column (11), the source identifications “s” are: 1 = APOGEE, 2 = GALAH, 3 = Gaia RVS, 4 = LAMOST, 5 = S5, 6 = SDSS, 7 = BOSS, 8 = ESPaDOnS (this work), 9 = AAOmega (from Ibata et al. 2017a), 10 = FLAMES (from Ibata et al. 2017a), 11 = UVES (from Odenkirchen et al. 2009), 12 = EFOSC (this work), 13 = UVES (this work), 14 = INT (this work), 15 = Yuan et al. (2022b), 16 = Li et al. (2021), 17 = Caldwell et al. (2020), 18 = T. S. Li et al. (2018), 19 = Koposov et al. (2018), 20 = G.-W. Li et al. (2018), 21 = Simon et al. (2020), 22 = Walker et al. (2015), 23 = VizieR, 24 = GES, 25 = DESI. Finally, column (12) provides a unique stream identification label. (This table is available in its entirety in machine-readable form.)

examined in the MCMC search.<fnref n="19" /> In Figure 11, we show the difference between the correction functions in the best-fit model, and in the best model when imposing a halo outer power-law parameter of β<sub>h</sub> = 3 (so as to be consistent with an NFW profile). In these two quite different Milky Way models, the correction functions are very similar, within ≈1 km s<sup>-1</sup> in line-of-sight velocity (which is the easiest observable to interpret when comparing streams, as there is no distance dependence).

## 8. Stream-fitting Results

We begin by attempting to fit a static axisymmetric model to the Galaxy, ignoring any stellar stream data; to this end, we use the objective function in Equation (10), with the first term on the right-hand side set to zero (i.e., we only fit to the ancillary data). For this first experiment, we fix the halo power-law parameters to γ<sub>h</sub> = 1 and β<sub>h</sub> = 3 so as to reproduce the expected behavior of dark matter halos (Navarro et al. 1997), set the halo matter flattening to q<sub>m,h</sub> = 1, set the thin and thick disk scale heights to 0.3 and 0.9 kpc respectively, and set the thick disk scale length to 3.6 kpc; these are all choices that were adopted by McMillan (2017). The main difference with respect to that study is that we are using the Eilers et al. (2019) rotation curve measurements. We find a thin disk scale length of h<sub>R,d</sub> = 2.22 ± 0.07(2.53 ± 0.14) kpc, a halo scale radius of r<sub>0,h</sub> = 12.4 ± 2.0(19.0 ± 4.9) kpc, a total stellar mass of M<sub>*</sub> = 4.63 ± 0.33(5.43 ± 0.57) × 10<sup>10</sup> M<sub>☉</sub>, M<sub>200</sub> = 0.86 ± 0.09(1.3 ± 0.3) × 10<sup>12</sup> M<sub>☉</sub>, and a mass interior to 50 kpc of M<sub>R&lt;50</sub> = 0.41 ± 0.02 × 10<sup>12</sup> M<sub>☉</sub>, where the values in brackets are those reported by McMillan (2017), when available. Figure 6 of McMillan (2017) shows their posterior PDF of M<sub>R&lt;50</sub> peaking at a considerably higher mass of M<sub>R&lt;50</sub> ≈ 0.51 × 10<sup>12</sup> M<sub>☉</sub>. Our aim with this test was to obtain a baseline for the

following stream analysis, and to explore how the improved data changes these inferences.

We now repeat this experiment, but add in the GD-1 stream, which is the structure for which we have the most numerous data (323 stars with line-of-sight velocity measurements and three RR Lyrae). This time we allow the halo mass flattening parameter q<sub>m,h</sub> to also be free. With this data set, we find h<sub>R,d</sub> = 2.16 ± 0.08 kpc, q<sub>m,h</sub> = 0.86 ± 0.05, r<sub>0,h</sub> = 10.7 ± 1.6 kpc, M<sub>*</sub> = 4.24 ± 0.35 × 10<sup>10</sup> M<sub>☉</sub>, M<sub>200</sub> = 0.77 ± 0.07 × 10<sup>12</sup> M<sub>☉</sub> and M<sub>R&lt;50</sub> = 0.39 ± 0.01 × 10<sup>12</sup> M<sub>☉</sub>, and the corresponding stream fit is shown in Figure 12. While the stream model captures the large-scale behavior of the kinematics of GD-1, it possesses some small, yet systematic, deviations from the proper motion trend in the data, especially in the range of 60° < φ<sub>1</sub> < 90°. The density distribution of the model particles (panel (e)) is also not a good match to the data, given that the model produces two high-density peaks at φ<sub>1</sub> ≈ 0° and φ<sub>1</sub> ≈ 70°, and a small remnant (at φ<sub>1</sub> = 29.7) still remains, but is not present in reality. However, the large peaks in the data at φ<sub>1</sub> ≈ 10° and ≈ 50°, which we ascribed in an earlier contribution (Ibata et al. 2020b) to epicyclic motion, can be seen to have counterparts in the model. Given this match, and also given that φ<sub>1</sub> = 32.4 corresponds to the stellar mass-weighted midpoint of the stream,<fnref n="20" /> we consider α = 157.74 to be currently the best estimate for the position of the center of the GD-1 progenitor. Henceforth, we fix the present-day position of GD-1's progenitor to that location.

We continue the procedure of rendering the model more realistic, adding in the five streams in our sample that have known globular cluster progenitors (Gjöll/NGC 3201, Fjörm/ M68, NGC 5466, Palomar 5, and M5). The phase-space parameters of these objects, and in particular their distances, are very well constrained, and in particular their distances, are very well constrained, which makes them especially powerful probes of the Galactic potential compared to a similarly sized sample of streams. In addition to the previously

<footnote n="19">¹⁹ Of course, this would be pointless in practice, as it would be more accurate and more computationally efficient to simply use the test-particle machinery developed in Section 6 and forgo the corrected orbit-fitting approach altogether.</footnote>

<footnote n="20">²⁰ We calculate this midpoint by assuming that the survey is complete to G<sub>0</sub> = 20 mag, and correct for missing stars given the distance model and the metallicity listed in Table 3. However, we did not correct for the masked region around M67.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Correlations between Milky Way parameters in a six-component Galaxy model with fixed bulge and gaseous disks, but free halo (four parameters, with fixed $\beta_h = 3$), thin disk (three parameters), and thick disk (three parameters). The useful quantity $M_{R<50}$ is not a fitting parameter, but is derived from the mass model determined by the other parameters. This mass model determines the potential in which 29 stream structures are fit. The red vertical line in each histogram marks the mode of the distribution. Using this mode, each distribution is cut into two subsamples, and the green lines show the 68th percentile lowest (left of the peak) and 68th percentile highest point (right of the peak) of the two subsamples. The contour lines in the correlation panels encompass $0.5\sigma$, $1\sigma$, $1.5\sigma$ and $2\sigma$ confidence regions.</figcaption></figure>

fitted Galactic parameters, the scale heights of the thin and thick disk components, the scale length of the thick disk, and the halo inner power exponent are also allowed to vary. Compared to the least constrained model we present below, only the halo outer power-law parameter remains fixed at

$\beta_h = 3$. The resulting fit has $h_{R,d} = 2.31 \pm 0.14$ kpc, $h_{z,d} = 345 \pm 9$ pc, $h_{R,t} = 2.02 \pm 0.42$ kpc, $h_{z,t} = 859 \pm 25$ pc, $\gamma_h = 0.89 \pm 0.19$, $q_{m,h} = 0.75 \pm 0.03$, $r_{0,h} = 22.1 \pm 2.8$ kpc, $M_* = 5.60 \pm 0.51 \times 10^{10} M_\odot$, $M_{200} = 1.42 \pm 0.25 \times 10^{12} M_\odot$, and $M_{R<50} = 0.50 \pm 0.03 \times 10^{12} M_\odot$. The most striking difference
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Correlations between Milky Way parameters in a six-component Galaxy model with fixed bulge and gaseous disks, but free halo (four parameters, with fixed $\beta_h = 3$), thin disk (three parameters), and thick disk (three parameters). The useful quantity $M_{R<50}$ is not a fitting parameter, but is derived from the mass model determined by the other parameters. This mass model determines the potential in which 29 stream structures are fit. The red vertical line in each histogram marks the mode of the distribution. Using this mode, each distribution is cut into two subsamples, and the green lines show the 68th percentile lowest (left of the peak) and 68th percentile highest point (right of the peak) of the two subsamples. The contour lines in the correlation panels encompass $0.5\sigma$, $1\sigma$, $1.5\sigma$ and $2\sigma$ confidence regions.</figcaption></figure>

fitted Galactic parameters, the scale heights of the thin and thick disk components, the scale length of the thick disk, and the halo inner power exponent are also allowed to vary. Compared to the least constrained model we present below, only the halo outer power-law parameter remains fixed at

$\beta_h = 3$. The resulting fit has $h_{R,d} = 2.31 \pm 0.14$ kpc, $h_{z,d} = 345 \pm 9$ pc, $h_{R,t} = 2.02 \pm 0.42$ kpc, $h_{z,t} = 859 \pm 25$ pc, $\gamma_h = 0.89 \pm 0.19$, $q_{m,h} = 0.75 \pm 0.03$, $r_{0,h} = 22.1 \pm 2.8$ kpc, $M_* = 5.60 \pm 0.51 \times 10^{10} M_\odot$, $M_{200} = 1.42 \pm 0.25 \times 10^{12} M_\odot$, and $M_{R<50} = 0.50 \pm 0.03 \times 10^{12} M_\odot$. The most striking difference
<figure><img src="image_17.png" /><figcaption>Figure 14. As in Figure 13, but here the halo outer power-law parameter $\beta_h$ is also allowed to vary, for a total of 11 free parameters of the Galactic mass model.</figcaption></figure>

to the previous fit is that the stream kinematics require the thick disk scale length to be much shorter than the 3.6 kpc measured by Jurić et al. (2008) and assumed by McMillan (2017); with a shorter thick disk scale length the halo scale radius $r_{0,h}$ becomes longer.

Exactly the same procedure is repeated, but this time using the full sample of 29 streams. As before, the Galactic potential model has a fixed bulge and gaseous disks, but a free halo, thin disk, and thick disk. The results of the MCMC exploration with the halo outer power-law parameter set to $\beta_h = 3$ are shown in Figure 13, while in Figure 14, $\beta_h$ is also allowed to vary.

Table 2 provides the summary statistics of the two experiments. We note that some of the marginalized distributions are clearly highly skewed, so the mean and rms values can be misleading. So instead, we report the distribution mode, which we find as follows: using the binning scheme shown in Figures 13 and 14 (30 bins between the minimum and maximum values), we identify the peak bin, and refine the position of the mode assuming that the distribution is locally quadratic over the region of the peak and the two adjacent bins. The marginalized distributions are then divided into two subsets, containing the data with lower and higher values than the mode, and in each
<figure><img src="image_17.png" /><figcaption>**Figure 14.** As in Figure 13, but here the halo outer power-law parameter $\beta_h$ is also allowed to vary, for a total of 11 free parameters of the Galactic mass model.</figcaption></figure>

to the previous fit is that the stream kinematics require the thick disk scale length to be much shorter than the 3.6 kpc measured by Jurić et al. (2008) and assumed by McMillan (2017); with a shorter thick disk scale length the halo scale radius $r_{0,h}$ becomes longer.

Exactly the same procedure is repeated, but this time using the full sample of 29 streams. As before, the Galactic potential model has a fixed bulge and gaseous disks, but a free halo, thin disk, and thick disk. The results of the MCMC exploration with the halo outer power-law parameter set to $\beta_h = 3$ are shown in Figure 13, while in Figure 14, $\beta_h$ is also allowed to vary.

Table 2 provides the summary statistics of the two experiments. We note that some of the marginalized distributions are clearly highly skewed, so the mean and rms values can be misleading. So instead, we report the distribution mode, which we find as follows: using the binning scheme shown in Figures 13 and 14 (30 bins between the minimum and maximum values), we identify the peak bin, and refine the position of the mode assuming that the distribution is locally quadratic over the region of the peak and the two adjacent bins. The marginalized distributions are then divided into two subsets, containing the data with lower and higher values than the mode, and in each
<figure><img src="image_18.png" /><figcaption>**Figure 15.** The total $M_{200}$ masses of the Milky Way solutions shown in Figures 13 (green) and 14 (blue).</figcaption></figure>

subset, we respectively find the position of the datum that is the 68.27% lowest and the 68.27% highest, and consider these values as the $1\sigma$ confidence intervals.

The halo mass flattening turns out to be very well constrained with $q_{m,h} = 0.749^{+0.026}_{-0.030}(0.709^{+0.033}_{-0.026})$, where the value in brackets corresponds to the fit in which $\beta_h$ is allowed to vary. Interestingly the (the negative of) the inner power-law slope of the halo is close to the universal NFW profile found in cosmological dark matter simulations. Nominally, we find a shallower outer power-law profile $\beta_h = 2.53^{+0.42}_{-0.16}$, but one that is consistent with the NFW profile at $1\sigma$. The values of the halo scale radius $r_{0,h} = 14.7^{+4.7}_{-1.0}(14.7^{+4.4}_{-3.3})$ kpc are closely consistent, although we note that the maximum-likelihood solutions are larger: $r_{0,h} = 21.4(19.5)$ kpc. We find that $r_{0,h}$ is strongly correlated with $\beta_h$, such that larger-scale radii require steeper outer power-law slopes.

The MCMC chain revealed strong correlations between the disk scale lengths and the central disk surface densities, and we also found that the combined mass of the disks $M_{d+t} = 4.20^{+0.44}_{-0.53}(4.34^{+0.50}_{-0.65}) \times 10^{10} M_{\odot}$ is much better constrained than that of the individual components. Likewise, we found that the local halo density $\rho_{h,\odot}$ is much better constrained than the central halo density $\rho_{0,h}$. These realizations motivated our use of the combination of parameters shown in Figures 13 and 14. (The central surface density of the disk can be recovered simply using the fact that the total mass of an exponential disk is $M = 2\pi\Sigma h_R^2$, while the central density of the halo is related to the local density via Equation (9).) For the thin and thick disk components, respectively, the scale heights are found to be $h_{z,d} = 0.347^{+0.007}_{-0.010}(0.345^{+0.009}_{-0.010})$ kpc and $h_{z,t} = 0.856^{+0.030}_{-0.020}(0.861^{+0.022}_{-0.025})$ kpc with scale lengths of $h_{R,d} = 2.172^{+0.184}_{-0.079}(2.198^{+0.153}_{-0.062})$ kpc and $h_{R,t} = 1.62^{+0.72}_{-0.13}(1.77^{+0.86}_{-0.17})$ kpc. The mass within $R=50$ kpc of these two models is $M_{R<50} = 0.458^{+0.030}_{-0.026}(0.490^{+0.024}_{-0.033}) \times 10^{12} M_{\odot}$. Determining $M_{200}$, however, requires a substantial extrapolation out of the region where we have data. The distribution of $M_{200}$ masses of these Milky Way models as derived from the MCMC chains is

<caption>**Table 2**</caption>

Statistics of the Parameters of the MCMC Chains Shown in Figures 13 and 14 Derived from Fitting a Flexible Galactic Model with the Sample of 29 Streams

<table><thead><tr><th>Quantity</th><th>β<sub>h</sub> = 3</th><th>L<sub>max</sub><sup>βh=3</sup></th><th>β<sub>h</sub> = free</th><th>L<sub>max</sub><sup>βh=free</sup></th></tr></thead><tbody><tr><td>ρ<sub>h,⊙</sub> [M<sub>☉</sub> pc<sup>-3</sup>/1000]</td><td>11.4 ± 0.7</td><td>11.4</td><td>11.6<sup>+0.7</sup><sub>-0.8</sub></td><td>11.2</td></tr><tr><td>q<sub>m,h</sub></td><td>0.749<sup>+0.026</sup><sub>-0.030</sub></td><td>0.735</td><td>0.709<sup>+0.033</sup><sub>-0.026</sub></td><td>0.727</td></tr><tr><td>γ<sub>h</sub></td><td>0.97<sup>+0.17</sup><sub>-0.21</sub></td><td>1.03</td><td>0.99<sup>+0.26</sup><sub>-0.30</sub></td><td>1.14</td></tr><tr><td>β<sub>h</sub></td><td>3</td><td>3</td><td>2.53<sup>+0.42</sup><sub>-0.16</sub></td><td>2.54</td></tr><tr><td>r<sub>0,h</sub> [kpc]</td><td>14.7<sup>+4.7</sup><sub>-1.0</sub></td><td>21.4</td><td>14.7<sup>+4.4</sup><sub>-3.3</sub></td><td>19.5</td></tr><tr><td>M<sub>d+t</sub>[10<sup>10</sup> M<sub>☉</sub>]</td><td>4.20<sup>+0.44</sup><sub>-0.53</sub></td><td>4.17</td><td>4.34<sup>+0.50</sup><sub>-0.65</sub></td><td>4.27</td></tr><tr><td>M<sub>t</sub>/M<sub>d+t</sub></td><td>0.095<sup>+0.143</sup><sub>-0.014</sub></td><td>0.201</td><td>0.09<sup>+0.10</sup><sub>-0.01</sub></td><td>0.14</td></tr><tr><td>h<sub>R,d</sub>[kpc]</td><td>2.172<sup>+0.184</sup><sub>-0.079</sub></td><td>2.241</td><td>2.198<sup>+0.153</sup><sub>-0.062</sub></td><td>2.251</td></tr><tr><td>h<sub>z,d</sub>[kpc]</td><td>0.347<sup>+0.007</sup><sub>-0.010</sub></td><td>0.348</td><td>0.345<sup>+0.009</sup><sub>-0.010</sub></td><td>0.349</td></tr><tr><td>h<sub>R,t</sub>[kpc]</td><td>1.62<sup>+0.72</sup><sub>-0.13</sub></td><td>1.74</td><td>1.77<sup>+0.86</sup><sub>-0.17</sub></td><td>2.07</td></tr><tr><td>h<sub>z,t</sub>[kpc]</td><td>0.856<sup>+0.030</sup><sub>-0.020</sub></td><td>0.858</td><td>0.861<sup>+0.022</sup><sub>-0.025</sub></td><td>0.860</td></tr><tr><td>M<sub>R<50</sub></td><td>0.458<sup>+0.030</sup><sub>-0.026</sub></td><td>0.477</td><td>0.490<sup>+0.024</sup><sub>-0.033</sub></td><td>0.518</td></tr><tr><td>M<sub>200</sub></td><td>1.09<sup>+0.19</sup><sub>-0.14</sub></td><td>1.24</td><td>1.70<sup>+0.20</sup><sub>-0.49</sub></td><td>2.00</td></tr><tr><td>M<sub>d</sub>[10<sup>10</sup> M<sub>☉</sub>]</td><td>3.44<sup>+0.44</sup><sub>-0.54</sub></td><td>3.34</td><td>3.52<sup>+0.54</sup><sub>-0.51</sub></td><td>3.65</td></tr><tr><td>M<sub>t</sub></td><td>0.40<sup>+0.61</sup><sub>-0.09</sub></td><td>0.84</td><td>0.39<sup>+0.43</sup><sub>-0.07</sub></td><td>0.62</td></tr></tbody></table>

**Note.** The first two sets of values use a fixed value $\beta_h = 3$ for the halo outer power-law parameter, while the final two sets of columns are for the experiment where $\beta_h$ is allowed to be free. The derived model masses are listed on the final four rows of the table. The $L_{\max}$ columns show the parameter values of the highest likelihood datum in the MCMC chains.

shown in Figure 15. The model with unconstrained $\beta_h$ is clearly affected by the imposed limit of $M_{200} < 2 \times 10^{12} M_{\odot}$. For this reason, we prefer the $\beta_h = 3$ model, although most of the other properties of the two models are consistent.

The fitted orbits of the progenitors of the 29 streams integrated in the best-fit potential from Figure 14 are shown in Figure 16. These orbits can be seen to approximate quite closely the observational data, demonstrating that the $\Delta\Theta$ correction functions are not required to be large for these systems. Similarly, in Figure 17, we compare these orbits to the simulated streams resulting from the test-particle integrations in the same potential, which again shows a close correspondence. Our fitting procedure has clearly managed to encapsulate the large-scale properties of the kinematic, and hence, of the dynamical behavior of this sample of streams.

## 9. Gaia DR3 Detections

Many of the structures displayed in Figures 1–9 were presented previously in Papers I and II, but 28 are new detections (the objects labeled “New” in Table 3 and stream #11 emanating from the Hydrus dwarf galaxy, Koposov et al. 2018). We have verified that they do not correspond to previously reported detections in the compendium of Mateu (2023) 21 or in the phase-space analysis of Dodd et al. (2023).

As may be expected, the stream-like structures discovered in each new Gaia data release generally have a smaller number of member stars than those found in earlier releases, simply because more populated streams were easier to find first. And despite our efforts in trying to obtain follow-up spectroscopy, in 16 of the new streams we have secured fewer than five radial velocity measurements. However, some of the new features

<footnote n="21">21 We note that stream #2 corresponds to the “SGP-S” structure discovered by Yang et al. (2022) that is not present in the Mateu (2023) catalog.</footnote>
<figure><img src="image_18.png" /><figcaption>Figure 15. The total $M_{200}$ masses of the Milky Way solutions shown in Figures 13 (green) and 14 (blue).</figcaption></figure>

subset, we respectively find the position of the datum that is the 68.27% lowest and the 68.27% highest, and consider these values as the $1\sigma$ confidence intervals.

The halo mass flattening turns out to be very well constrained with $q_{m,h} = 0.749^{+0.026}_{-0.030}(0.709^{+0.033}_{-0.026})$, where the value in brackets corresponds to the fit in which $\beta_h$ is allowed to vary. Interestingly the (the negative of) the inner power-law slope of the halo is close to the universal NFW profile found in cosmological dark matter simulations. Nominally, we find a shallower outer power-law profile $\beta_h = 2.53^{+0.42}_{-0.16}$, but one that is consistent with the NFW profile at $1\sigma$. The values of the halo scale radius $r_{0,h} = 14.7^{+4.7}_{-1.0}(14.7^{+4.4}_{-3.3})$ kpc are closely consistent, although we note that the maximum-likelihood solutions are larger: $r_{0,h} = 21.4(19.5)$ kpc. We find that $r_{0,h}$ is strongly correlated with $\beta_h$, such that larger-scale radii require steeper outer power-law slopes.

The MCMC chain revealed strong correlations between the disk scale lengths and the central disk surface densities, and we also found that the combined mass of the disks $M_{d+t} = 4.20^{+0.44}_{-0.53}(4.34^{+0.50}_{-0.65}) \times 10^{10} M_{\odot}$ is much better constrained than that of the individual components. Likewise, we found that the local halo density $\rho_{h,\odot}$ is much better constrained than the central halo density $\rho_{0,h}$. These realizations motivated our use of the combination of parameters shown in Figures 13 and 14. (The central surface density of the disk can be recovered simply using the fact that the total mass of an exponential disk is $M = 2\pi\Sigma h_R^2$, while the central density of the halo is related to the local density via Equation (9).) For the thin and thick disk components, respectively, the scale heights are found to be $h_{z,d} = 0.347^{+0.007}_{-0.010}(0.345^{+0.009}_{-0.010})$ kpc and $h_{z,t} = 0.856^{+0.030}_{-0.020}(0.861^{+0.022}_{-0.025})$ kpc with scale lengths of $h_{R,d} = 2.172^{+0.184}_{-0.079}(2.198^{+0.153}_{-0.062})$ kpc and $h_{R,t} = 1.62^{+0.72}_{-0.13}(1.77^{+0.86}_{-0.17})$ kpc. The mass within $R=50$ kpc of these two models is $M_{R<50} = 0.458^{+0.030}_{-0.026}(0.490^{+0.024}_{-0.033}) \times 10^{12} M_{\odot}$. Determining $M_{200}$, however, requires a substantial extrapolation out of the region where we have data. The distribution of $M_{200}$ masses of these Milky Way models as derived from the MCMC chains is

<caption>Table 2</caption>

Statistics of the Parameters of the MCMC Chains Shown in Figures 13 and 14 Derived from Fitting a Flexible Galactic Model with the Sample of 29 Streams

<table><thead><tr><th>Quantity</th><th>β<sub>h</sub> = 3</th><th>L<sub>max</sub><sup>βh=3</sup></th><th>β<sub>h</sub> = free</th><th>L<sub>max</sub><sup>βh=free</sup></th></tr></thead><tbody><tr><td>ρ<sub>h,⊙</sub> [M<sub>☉</sub> pc<sup>-3</sup>/1000]</td><td>11.4 ± 0.7</td><td>11.4</td><td>11.6<sup>+0.7</sup><sub>-0.8</sub></td><td>11.2</td></tr><tr><td>q<sub>m,h</sub></td><td>0.749<sup>+0.026</sup><sub>-0.030</sub></td><td>0.735</td><td>0.709<sup>+0.033</sup><sub>-0.026</sub></td><td>0.727</td></tr><tr><td>γ<sub>h</sub></td><td>0.97<sup>+0.17</sup><sub>-0.21</sub></td><td>1.03</td><td>0.99<sup>+0.26</sup><sub>-0.30</sub></td><td>1.14</td></tr><tr><td>β<sub>h</sub></td><td>3</td><td>3</td><td>2.53<sup>+0.42</sup><sub>-0.16</sub></td><td>2.54</td></tr><tr><td>r<sub>0,h</sub> [kpc]</td><td>14.7<sup>+4.7</sup><sub>-1.0</sub></td><td>21.4</td><td>14.7<sup>+4.4</sup><sub>-3.3</sub></td><td>19.5</td></tr><tr><td>M<sub>d+t</sub>[10<sup>10</sup> M<sub>☉</sub>]</td><td>4.20<sup>+0.44</sup><sub>-0.53</sub></td><td>4.17</td><td>4.34<sup>+0.50</sup><sub>-0.65</sub></td><td>4.27</td></tr><tr><td>M<sub>t</sub>/M<sub>d+t</sub></td><td>0.095<sup>+0.143</sup><sub>-0.014</sub></td><td>0.201</td><td>0.09<sup>+0.10</sup><sub>-0.01</sub></td><td>0.14</td></tr><tr><td>h<sub>R,d</sub>[kpc]</td><td>2.172<sup>+0.184</sup><sub>-0.079</sub></td><td>2.241</td><td>2.198<sup>+0.153</sup><sub>-0.062</sub></td><td>2.251</td></tr><tr><td>h<sub>z,d</sub>[kpc]</td><td>0.347<sup>+0.007</sup><sub>-0.010</sub></td><td>0.348</td><td>0.345<sup>+0.009</sup><sub>-0.010</sub></td><td>0.349</td></tr><tr><td>h<sub>R,t</sub>[kpc]</td><td>1.62<sup>+0.72</sup><sub>-0.13</sub></td><td>1.74</td><td>1.77<sup>+0.86</sup><sub>-0.17</sub></td><td>2.07</td></tr><tr><td>h<sub>z,t</sub>[kpc]</td><td>0.856<sup>+0.030</sup><sub>-0.020</sub></td><td>0.858</td><td>0.861<sup>+0.022</sup><sub>-0.025</sub></td><td>0.860</td></tr><tr><td>M<sub>R<50</sub></td><td>0.458<sup>+0.030</sup><sub>-0.026</sub></td><td>0.477</td><td>0.490<sup>+0.024</sup><sub>-0.033</sub></td><td>0.518</td></tr><tr><td>M<sub>200</sub></td><td>1.09<sup>+0.19</sup><sub>-0.14</sub></td><td>1.24</td><td>1.70<sup>+0.20</sup><sub>-0.49</sub></td><td>2.00</td></tr><tr><td>M<sub>d</sub></td><td>3.44<sup>+0.44</sup><sub>-0.54</sub></td><td>3.34</td><td>3.52<sup>+0.54</sup><sub>-0.51</sub></td><td>3.65</td></tr><tr><td>M<sub>t</sub></td><td>0.40<sup>+0.61</sup><sub>-0.09</sub></td><td>0.84</td><td>0.39<sup>+0.43</sup><sub>-0.07</sub></td><td>0.62</td></tr></tbody></table>

Note. The first two sets of values use a fixed value $\beta_h = 3$ for the halo outer power-law parameter, while the final two sets of columns are for the experiment where $\beta_h$ is allowed to be free. The derived model masses are listed on the final four rows of the table. The $L_{\max}$ columns show the parameter values of the highest likelihood datum in the MCMC chains.

shown in Figure 15. The model with unconstrained $\beta_h$ is clearly affected by the imposed limit of $M_{200} < 2 \times 10^{12} M_{\odot}$. For this reason, we prefer the $\beta_h = 3$ model, although most of the other properties of the two models are consistent.

The fitted orbits of the progenitors of the 29 streams integrated in the best-fit potential from Figure 14 are shown in Figure 16. These orbits can be seen to approximate quite closely the observational data, demonstrating that the $\Delta\Theta$ correction functions are not required to be large for these systems. Similarly, in Figure 17, we compare these orbits to the simulated streams resulting from the test-particle integrations in the same potential, which again shows a close correspondence. Our fitting procedure has clearly managed to encapsulate the large-scale properties of the kinematic, and hence, of the dynamical behavior of this sample of streams.

## 9. Gaia DR3 Detections

Many of the structures displayed in Figures 1–9 were presented previously in Papers I and II, but 28 are new detections (the objects labeled “New” in Table 3 and stream #11 emanating from the Hydrus dwarf galaxy, Koposov et al. 2018). We have verified that they do not correspond to previously reported detections in the compendium of Mateu (2023) 21 or in the phase-space analysis of Dodd et al. (2023).

As may be expected, the stream-like structures discovered in each new Gaia data release generally have a smaller number of member stars than those found in earlier releases, simply because more populated streams were easier to find first. And despite our efforts in trying to obtain follow-up spectroscopy, in 16 of the new streams we have secured fewer than five radial velocity measurements. However, some of the new features

<footnote n="21">21 We note that stream #2 corresponds to the “SGP-S” structure discovered by Yang et al. (2022) that is not present in the Mateu (2023) catalog.</footnote>
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Orbital fits to the sample of 29 long streams integrated in the best potential model with $\beta_h$ allowed to vary. The large and small four-cornered markers show, respectively, the positions of the center of globular clusters and RR Lyrae. The filled circles are stars in the streams with measured line-of-sight velocities. For each parameter on the ordinate axes, we show the corresponding behavior of the best-fit orbits. The color scheme of Figure 6 is used here.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Orbital fits to the sample of 29 long streams integrated in the best potential model with $\beta_h$ allowed to vary. The large and small four-cornered markers show, respectively, the positions of the center of globular clusters and RR Lyrae. The filled circles are stars in the streams with measured line-of-sight velocities. For each parameter on the ordinate axes, we show the corresponding behavior of the best-fit orbits. The color scheme of Figure 6 is used here.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 17.** As in Figure 16, but comparing the progenitor orbits to the test-particle stream simulations.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 17. As in Figure 16, but comparing the progenitor orbits to the test-particle stream simulations.</figcaption></figure>
<table><caption>Table 3
Streams Detected in Gaia DR3</caption><thead><tr><th>s<sub>ID</sub></th><th>Name</th><th>Sample</th><th>α<sub>0</sub></th><th>α<sub>pole</sub></th><th>δ<sub>pole</sub></th><th>n</th><th>n<sub>v</sub></th><th>M<sub>*</sub><br>(M<sub>☉</sub>)</th><th>⟨[Fe/H]⟩</th><th>References</th><th>Remarks</th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th><th>(4)</th><th>(5)</th><th>(6)</th><th>(7)</th><th>(8)</th><th>(9)</th><th>(10)</th><th>(11)</th><th>(12)</th></tr></thead><tbody><tr><td>1</td><td>C-20</td><td>1</td><td>0.561</td><td>275.038</td><td>-30.9</td><td>29</td><td>8</td><td>1200</td><td>-2.93</td><td>(1)</td><td>...</td></tr><tr><td>2</td><td>SGP-S</td><td>0</td><td>6.705</td><td>273.23</td><td>-5.703</td><td>36</td><td>1</td><td>530</td><td>-2.22</td><td>(2)</td><td>...</td></tr><tr><td>3</td><td>NGC 288</td><td>0</td><td>13.188</td><td>101.35</td><td>1.686</td><td>173</td><td>7</td><td>2200</td><td>-1.32</td><td>(3)</td><td>...</td></tr><tr><td>4</td><td>ATLAS</td><td>1</td><td>20.185</td><td>75.623</td><td>47.653</td><td>208</td><td>80</td><td>19,000</td><td>-2.2</td><td>(4)</td><td>...</td></tr><tr><td>5</td><td>Phoenix</td><td>1</td><td>24.341</td><td>309.372</td><td>12.542</td><td>58</td><td>26</td><td>3000</td><td>-2.7</td><td>(5)</td><td>...</td></tr><tr><td>6</td><td>Kwando</td><td>0</td><td>24.930</td><td>112.058</td><td>9.835</td><td>138</td><td>62</td><td>40,000</td><td>-2.29</td><td>(6)</td><td>Very wide, rough [Fe/H] estimate</td></tr><tr><td>7</td><td>New-1</td><td>0</td><td>27.124</td><td>296.02</td><td>-1.042</td><td>11</td><td>3</td><td>420</td><td>-1.85</td><td>(2)</td><td>...</td></tr><tr><td>8</td><td>New-2</td><td>0</td><td>34.127</td><td>53.727</td><td>-85.092</td><td>90</td><td>18</td><td>73</td><td>-2.07</td><td>(7)</td><td>...</td></tr><tr><td>9</td><td>C-13</td><td>0</td><td>35.678</td><td>311.207</td><td>-6.984</td><td>100</td><td>5</td><td>750</td><td>-0.25</td><td>(8)</td><td>Velocities unclear, strange CMD</td></tr><tr><td>10</td><td>Gaia-12</td><td>1</td><td>37.357</td><td>147.363</td><td>43.851</td><td>46</td><td>5</td><td>1200</td><td>-3.28</td><td>(2)</td><td>...</td></tr><tr><td>11</td><td>Hydrus</td><td>0</td><td>37.389</td><td>63.965</td><td>9.619</td><td>39</td><td>14</td><td>7700</td><td>-3.48</td><td>(2)</td><td>New stream of known structure</td></tr><tr><td>12</td><td>NGC 1261</td><td>0</td><td>48.068</td><td>18.701</td><td>31.13</td><td>173</td><td>23</td><td>8100</td><td>-1.27</td><td>(3)</td><td>Velocities messy, orbit apocenter</td></tr><tr><td>13</td><td>NGC 1261a</td><td>0</td><td>48.068</td><td>153.876</td><td>-10.473</td><td>136</td><td>4</td><td>9700</td><td>-1.27</td><td>(3)</td><td>Related to N1261?</td></tr><tr><td>14</td><td>NGC 1261b</td><td>0</td><td>48.068</td><td>84.128</td><td>32.219</td><td>70</td><td>12</td><td>5100</td><td>-1.27</td><td>(3)</td><td>Related to N1261?</td></tr><tr><td>15</td><td>Indus</td><td>0</td><td>56.206</td><td>38.115</td><td>19.901</td><td>1021</td><td>96</td><td>24,000</td><td>-2.09</td><td>(9)</td><td>...</td></tr><tr><td>16</td><td>NGC 1851</td><td>0</td><td>78.528</td><td>151.161</td><td>19.287</td><td>105</td><td>11</td><td>3400</td><td>-1.18</td><td>(3)</td><td>...</td></tr><tr><td>17</td><td>New-3</td><td>0</td><td>82.663</td><td>254.437</td><td>-72.963</td><td>120</td><td>10</td><td>120</td><td>-2.71</td><td>(10)</td><td>Two clumps?</td></tr><tr><td>18</td><td>New-4</td><td>1</td><td>90.988</td><td>51.229</td><td>45.719</td><td>29</td><td>8</td><td>730</td><td>-1.76</td><td>(2)</td><td>...</td></tr><tr><td>19</td><td>Leiptr</td><td>1</td><td>91.505</td><td>152.738</td><td>39.307</td><td>412</td><td>37</td><td>3000</td><td>-2.17</td><td>(2)</td><td>...</td></tr><tr><td>20</td><td>C-12</td><td>0</td><td>94.401</td><td>289.858</td><td>-59.422</td><td>929</td><td>13</td><td>14,000</td><td>-1.34</td><td>(10)</td><td>Two overlapping structures?</td></tr><tr><td>21</td><td>New-5</td><td>0</td><td>98.455</td><td>249.927</td><td>-87.804</td><td>13</td><td>9</td><td>16</td><td>-1.02</td><td>(2)</td><td>Two structures?</td></tr><tr><td>22</td><td>New-6</td><td>1</td><td>99.776</td><td>106.87</td><td>-14.425</td><td>146</td><td>12</td><td>340</td><td>-2.18</td><td>(7)</td><td>...</td></tr><tr><td>23</td><td>NGC 2298</td><td>0</td><td>102.248</td><td>137.971</td><td>48.3</td><td>125</td><td>1</td><td>2000</td><td>-1.92</td><td>(3)</td><td>...</td></tr><tr><td>24</td><td>C-25</td><td>0</td><td>112.659</td><td>324.848</td><td>44.684</td><td>121</td><td>25</td><td>1300</td><td>-2.3</td><td>(1)</td><td>Related to C-11?</td></tr><tr><td>25</td><td>C-11</td><td>0</td><td>117.062</td><td>349.03</td><td>24.975</td><td>112</td><td>9</td><td>1100</td><td>-2.91</td><td>(1)</td><td>Velocities unclear</td></tr><tr><td>26</td><td>NGC 2808</td><td>0</td><td>138.013</td><td>171.607</td><td>21.534</td><td>292</td><td>10</td><td>7300</td><td>-1.14</td><td>(3)</td><td>...</td></tr><tr><td>27</td><td>New-7</td><td>0</td><td>139.625</td><td>318.111</td><td>66.353</td><td>37</td><td>11</td><td>110</td><td>-2.23</td><td>(7)</td><td>Velocities unclear</td></tr><tr><td>28</td><td>New-8</td><td>0</td><td>142.991</td><td>237.5</td><td>6.097</td><td>7</td><td>1</td><td>3</td><td>-1.99</td><td>(7)</td><td>...</td></tr><tr><td>29</td><td>New-9</td><td>0</td><td>144.479</td><td>318.071</td><td>58.997</td><td>10</td><td>7</td><td>38</td><td>-1.95</td><td>(7)</td><td>...</td></tr><tr><td>30</td><td>C-10</td><td>0</td><td>145.970</td><td>55.813</td><td>5.518</td><td>158</td><td>5</td><td>1800</td><td>-0.91</td><td>(8)</td><td>Thick RGB</td></tr><tr><td>31</td><td>New-10</td><td>0</td><td>147.682</td><td>61.538</td><td>17.134</td><td>47</td><td>7</td><td>100</td><td>-0.99</td><td>(10)</td><td>...</td></tr><tr><td>32</td><td>New-11</td><td>0</td><td>153.299</td><td>234.748</td><td>-11.648</td><td>18</td><td>6</td><td>110</td><td>-2.03</td><td>(7)</td><td>Velocities unclear</td></tr><tr><td>33</td><td>Gaia-10</td><td>1</td><td>153.312</td><td>156.208</td><td>-74.386</td><td>141</td><td>26</td><td>7600</td><td>-1.4</td><td>(1)</td><td>...</td></tr><tr><td>34</td><td>Gjoll</td><td>1</td><td>154.403</td><td>158.539</td><td>43.545</td><td>607</td><td>40</td><td>2100</td><td>-1.63</td><td>(3)</td><td>...</td></tr><tr><td>35</td><td>C-9</td><td>0</td><td>155.889</td><td>115.351</td><td>-66.055</td><td>183</td><td>35</td><td>630</td><td>-0.72</td><td>(1)</td><td>Several clumps?</td></tr><tr><td>36</td><td>C-24</td><td>0</td><td>158.186</td><td>228.19</td><td>-12.131</td><td>244</td><td>27</td><td>2700</td><td>-0.93</td><td>(1)</td><td>Velocities unclear</td></tr><tr><td>37</td><td>New-12</td><td>0</td><td>161.861</td><td>84.236</td><td>14.306</td><td>89</td><td>3</td><td>580</td><td>-1.85</td><td>(2)</td><td>...</td></tr><tr><td>38</td><td>New-13</td><td>0</td><td>162.209</td><td>264.091</td><td>-13.901</td><td>47</td><td>1</td><td>250</td><td>-1.19</td><td>(2)</td><td>High velocity</td></tr><tr><td>39</td><td>Slidr</td><td>1</td><td>166.076</td><td>92.169</td><td>-64.021</td><td>330</td><td>34</td><td>840</td><td>-1.7</td><td>(1)</td><td>...</td></tr><tr><td>40</td><td>New-14</td><td>0</td><td>170.248</td><td>69.207</td><td>-26.38</td><td>18</td><td>4</td><td>28</td><td>-0.95</td><td>(2)</td><td>Compact</td></tr><tr><td>41</td><td>Ylgr</td><td>1</td><td>174.004</td><td>256.584</td><td>17.379</td><td>919</td><td>20</td><td>11,000</td><td>-2.09</td><td>(7)</td><td>...</td></tr><tr><td>42</td><td>Sylgr</td><td>0</td><td>176.261</td><td>273.255</td><td>-60.454</td><td>256</td><td>33</td><td>700</td><td>-2.92</td><td>(11)</td><td>High-velocity scatter</td></tr><tr><td>43</td><td>New-15</td><td>0</td><td>178.134</td><td>66.201</td><td>58.999</td><td>152</td><td>23</td><td>960</td><td>-2.01</td><td>(7)</td><td>...</td></tr><tr><td>44</td><td>Gaia-7</td><td>1</td><td>181.295</td><td>78.993</td><td>-35.727</td><td>282</td><td>15</td><td>1300</td><td>-0.88</td><td>(10)</td><td>...</td></tr><tr><td>45</td><td>Gaia-8</td><td>1</td><td>186.342</td><td>305.604</td><td>-48.129</td><td>453</td><td>18</td><td>3600</td><td>-1.77</td><td>(10)</td><td>...</td></tr><tr><td>46</td><td>New-16</td><td>0</td><td>188.285</td><td>102.418</td><td>5.746</td><td>44</td><td>3</td><td>170</td><td>-1.62</td><td>(2)</td><td>...</td></tr><tr><td>47</td><td>Fjorm</td><td>1</td><td>189.867</td><td>103.086</td><td>19.455</td><td>297</td><td>29</td><td>1700</td><td>-2.23</td><td>(3)</td><td>...</td></tr><tr><td>48</td><td>Orphan</td><td>0</td><td>191.105</td><td>68.052</td><td>-15.708</td><td>867</td><td>247</td><td>130,000</td><td>-1.9</td><td>(12)</td><td>Massive system</td></tr><tr><td>49</td><td>Gaia-1</td><td>1</td><td>192.56</td><td>287.833</td><td>-35.545</td><td>200</td><td>28</td><td>1100</td><td>-1.8</td><td>(7)</td><td>...</td></tr><tr><td>50</td><td>C-23</td><td>1</td><td>195.564</td><td>122.347</td><td>-16.78</td><td>29</td><td>6</td><td>310</td><td>-2.36</td><td>(1)</td><td>...</td></tr><tr><td>51</td><td>LMS-1</td><td>0</td><td>198.23</td><td>43.633</td><td>62.255</td><td>358</td><td>29</td><td>15,000</td><td>-2.09</td><td>(13)</td><td>Massive and wide system</td></tr><tr><td>52</td><td>C-22</td><td>1</td><td>199.904</td><td>56.016</td><td>68.974</td><td>39</td><td>11</td><td>470</td><td>-2.05</td><td>(8)</td><td>...</td></tr><tr><td>53</td><td>GD-1</td><td>1</td><td>200</td><td>214.928</td><td>-29.881</td><td>1468</td><td>323</td><td>14,000</td><td>-2.49</td><td>(1)</td><td>...</td></tr><tr><td>54</td><td>Fimbulthul</td><td>0</td><td>201.697</td><td>304.173</td><td>-22.859</td><td>3724</td><td>29</td><td>16,000</td><td>-1.53</td><td>(3)</td><td>Massive system</td></tr><tr><td>55</td><td>Fimbulthul-S</td><td>0</td><td>201.697</td><td>16.148</td><td>1.199</td><td>1734</td><td>25</td><td>4400</td><td>-1.53</td><td>(10)</td><td>Massive system</td></tr><tr><td>56</td><td>New-17</td><td>0</td><td>210.104</td><td>147.375</td><td>-50.685</td><td>13</td><td>0</td><td>29</td><td>-2.35</td><td>(2)</td><td>...</td></tr><tr><td>57</td><td>NGC5466</td><td>1</td><td>211.364</td><td>161.594</td><td>-50.069</td><td>43</td><td>6</td><td>1900</td><td>-1.98</td><td>(3)</td><td>...</td></tr><tr><td>58</td><td>New-18</td><td>0</td><td>214.041</td><td>302.311</td><td>-1.888</td><td>107</td><td>1</td><td>280</td><td>-0.94</td><td>(2)</td><td>...</td></tr><tr><td>59</td><td>Gaia-6</td><td>1</td><td>214.600</td><td>268.951</td><td>-41.694</td><td>145</td><td>14</td><td>1800</td><td>-1.53</td><td>(1)</td><td>...</td></tr><tr><td>60</td><td>New-19</td><td>0</td><td>221.839</td><td>294.781</td><td>-22.535</td><td>40</td><td>4</td><td>100</td><td>-0.79</td><td>(2)</td><td>Compact, velocities unclear</td></tr><tr><td>61</td><td>Pal-5</td><td>1</td><td>229.022</td><td>320.298</td><td>-54.042</td><td>129</td><td>69</td><td>17,000</td><td>-1.41</td><td>(3)</td><td>...</td></tr><tr><td>62</td><td>M5</td><td>1</td><td>229.638</td><td>325.223</td><td>63.075</td><td>91</td><td>11</td><td>710</td><td>-1.29</td><td>(3)</td><td>...</td></tr></tbody></table>
<table><caption>Table 3
Streams Detected in Gaia DR3</caption><thead><tr><th>s<sub>ID</sub></th><th>Name</th><th>Sample</th><th>α<sub>0</sub></th><th>α<sub>pole</sub></th><th>δ<sub>pole</sub></th><th>n</th><th>n<sub>v</sub></th><th>M<sub>*</sub><br>(M<sub>c</sub>)</th><th>⟨[Fe/H]⟩</th><th>References</th><th>Remarks</th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th><th>(4)</th><th>(5)</th><th>(6)</th><th>(7)</th><th>(8)</th><th>(9)</th><th>(10)</th><th>(11)</th><th>(12)</th></tr></thead><tbody><tr><td>1</td><td>C-20</td><td>1</td><td>0.561</td><td>275.038</td><td>-30.9</td><td>29</td><td>8</td><td>1200</td><td>-2.93</td><td>(1)</td><td>...</td></tr><tr><td>2</td><td>SGP-S</td><td>0</td><td>6.705</td><td>273.23</td><td>-5.703</td><td>36</td><td>1</td><td>530</td><td>-2.22</td><td>(2)</td><td>...</td></tr><tr><td>3</td><td>NGC 288</td><td>0</td><td>13.188</td><td>101.35</td><td>1.686</td><td>173</td><td>7</td><td>2200</td><td>-1.32</td><td>(3)</td><td>...</td></tr><tr><td>4</td><td>ATLAS</td><td>1</td><td>20.185</td><td>75.623</td><td>47.653</td><td>208</td><td>80</td><td>19,000</td><td>-2.2</td><td>(4)</td><td>...</td></tr><tr><td>5</td><td>Phoenix</td><td>1</td><td>24.341</td><td>309.372</td><td>12.542</td><td>58</td><td>26</td><td>3000</td><td>-2.7</td><td>(5)</td><td>...</td></tr><tr><td>6</td><td>Kwando</td><td>0</td><td>24.930</td><td>112.058</td><td>9.835</td><td>138</td><td>62</td><td>40,000</td><td>-2.29</td><td>(6)</td><td>Very wide, rough [Fe/H] estimate</td></tr><tr><td>7</td><td>New-1</td><td>0</td><td>27.124</td><td>296.02</td><td>-1.042</td><td>11</td><td>3</td><td>420</td><td>-1.85</td><td>(2)</td><td>...</td></tr><tr><td>8</td><td>New-2</td><td>0</td><td>34.127</td><td>53.727</td><td>-85.092</td><td>90</td><td>18</td><td>73</td><td>-2.07</td><td>(7)</td><td>...</td></tr><tr><td>9</td><td>C-13</td><td>0</td><td>35.678</td><td>311.207</td><td>-6.984</td><td>100</td><td>5</td><td>750</td><td>-0.25</td><td>(8)</td><td>Velocities unclear, strange CMD</td></tr><tr><td>10</td><td>Gaia-12</td><td>1</td><td>37.357</td><td>147.363</td><td>43.851</td><td>46</td><td>5</td><td>1200</td><td>-3.28</td><td>(2)</td><td>...</td></tr><tr><td>11</td><td>Hydrus</td><td>0</td><td>37.389</td><td>63.965</td><td>9.619</td><td>39</td><td>14</td><td>7700</td><td>-3.48</td><td>(2)</td><td>New stream of known structure</td></tr><tr><td>12</td><td>NGC 1261</td><td>0</td><td>48.068</td><td>18.701</td><td>31.13</td><td>173</td><td>23</td><td>8100</td><td>-1.27</td><td>(3)</td><td>Velocities messy, orbit apocenter</td></tr><tr><td>13</td><td>NGC 1261a</td><td>0</td><td>48.068</td><td>153.876</td><td>-10.473</td><td>136</td><td>4</td><td>9700</td><td>-1.27</td><td>(3)</td><td>Related to N1261?</td></tr><tr><td>14</td><td>NGC 1261b</td><td>0</td><td>48.068</td><td>84.128</td><td>32.219</td><td>70</td><td>12</td><td>5100</td><td>-1.27</td><td>(3)</td><td>Related to N1261?</td></tr><tr><td>15</td><td>Indus</td><td>0</td><td>56.206</td><td>38.115</td><td>19.901</td><td>1021</td><td>96</td><td>24,000</td><td>-2.09</td><td>(9)</td><td>...</td></tr><tr><td>16</td><td>NGC 1851</td><td>0</td><td>78.528</td><td>151.161</td><td>19.287</td><td>105</td><td>11</td><td>3400</td><td>-1.18</td><td>(3)</td><td>...</td></tr><tr><td>17</td><td>New-3</td><td>0</td><td>82.663</td><td>254.437</td><td>-72.963</td><td>120</td><td>10</td><td>120</td><td>-2.71</td><td>(10)</td><td>Two clumps?</td></tr><tr><td>18</td><td>New-4</td><td>1</td><td>90.988</td><td>51.229</td><td>45.719</td><td>29</td><td>8</td><td>730</td><td>-1.76</td><td>(2)</td><td>...</td></tr><tr><td>19</td><td>Leipt</td><td>1</td><td>91.505</td><td>152.738</td><td>39.307</td><td>412</td><td>37</td><td>3000</td><td>-2.17</td><td>(2)</td><td>...</td></tr><tr><td>20</td><td>C-12</td><td>0</td><td>94.401</td><td>289.858</td><td>-59.422</td><td>929</td><td>13</td><td>14,000</td><td>-1.34</td><td>(10)</td><td>Two overlapping structures?</td></tr><tr><td>21</td><td>New-5</td><td>0</td><td>98.455</td><td>249.927</td><td>-87.804</td><td>13</td><td>9</td><td>16</td><td>-1.02</td><td>(2)</td><td>Two structures?</td></tr><tr><td>22</td><td>New-6</td><td>1</td><td>99.776</td><td>106.87</td><td>-14.425</td><td>146</td><td>12</td><td>340</td><td>-2.18</td><td>(7)</td><td>...</td></tr><tr><td>23</td><td>NGC 2298</td><td>0</td><td>102.248</td><td>137.971</td><td>48.3</td><td>125</td><td>1</td><td>2000</td><td>-1.92</td><td>(3)</td><td>...</td></tr><tr><td>24</td><td>C-25</td><td>0</td><td>112.659</td><td>324.848</td><td>44.684</td><td>121</td><td>25</td><td>1300</td><td>-2.3</td><td>(1)</td><td>Related to C-11?</td></tr><tr><td>25</td><td>C-11</td><td>0</td><td>117.062</td><td>349.03</td><td>24.975</td><td>112</td><td>9</td><td>1100</td><td>-2.91</td><td>(1)</td><td>Velocities unclear</td></tr><tr><td>26</td><td>NGC 2808</td><td>0</td><td>138.013</td><td>171.607</td><td>21.534</td><td>292</td><td>10</td><td>7300</td><td>-1.14</td><td>(3)</td><td>...</td></tr><tr><td>27</td><td>New-7</td><td>0</td><td>139.625</td><td>318.111</td><td>66.353</td><td>37</td><td>11</td><td>110</td><td>-2.23</td><td>(7)</td><td>Velocities unclear</td></tr><tr><td>28</td><td>New-8</td><td>0</td><td>142.991</td><td>237.5</td><td>6.097</td><td>7</td><td>1</td><td>3</td><td>-1.99</td><td>(7)</td><td>...</td></tr><tr><td>29</td><td>New-9</td><td>0</td><td>144.479</td><td>318.071</td><td>58.997</td><td>10</td><td>7</td><td>38</td><td>-1.95</td><td>(7)</td><td>...</td></tr><tr><td>30</td><td>C-10</td><td>0</td><td>145.970</td><td>55.813</td><td>5.518</td><td>158</td><td>5</td><td>1800</td><td>-0.91</td><td>(8)</td><td>Thick RGB</td></tr><tr><td>31</td><td>New-10</td><td>0</td><td>147.682</td><td>61.538</td><td>17.134</td><td>47</td><td>7</td><td>100</td><td>-0.99</td><td>(10)</td><td>...</td></tr><tr><td>32</td><td>New-11</td><td>0</td><td>153.299</td><td>234.748</td><td>-11.648</td><td>18</td><td>6</td><td>110</td><td>-2.03</td><td>(7)</td><td>Velocities unclear</td></tr><tr><td>33</td><td>Gaia-10</td><td>1</td><td>153.312</td><td>156.208</td><td>-74.386</td><td>141</td><td>26</td><td>7600</td><td>-1.4</td><td>(1)</td><td>...</td></tr><tr><td>34</td><td>Gjoll</td><td>1</td><td>154.403</td><td>158.539</td><td>43.545</td><td>607</td><td>40</td><td>2100</td><td>-1.63</td><td>(3)</td><td>...</td></tr><tr><td>35</td><td>C-9</td><td>0</td><td>155.889</td><td>115.351</td><td>-66.055</td><td>183</td><td>35</td><td>630</td><td>-0.72</td><td>(1)</td><td>Several clumps?</td></tr><tr><td>36</td><td>C-24</td><td>0</td><td>158.186</td><td>228.19</td><td>-12.131</td><td>244</td><td>27</td><td>2700</td><td>-0.93</td><td>(1)</td><td>Velocities unclear</td></tr><tr><td>37</td><td>New-12</td><td>0</td><td>161.861</td><td>84.236</td><td>14.306</td><td>89</td><td>3</td><td>580</td><td>-1.85</td><td>(2)</td><td>...</td></tr><tr><td>38</td><td>New-13</td><td>0</td><td>162.209</td><td>264.091</td><td>-13.901</td><td>47</td><td>1</td><td>250</td><td>-1.19</td><td>(2)</td><td>High velocity</td></tr><tr><td>39</td><td>Slidr</td><td>1</td><td>166.076</td><td>92.169</td><td>-64.021</td><td>330</td><td>34</td><td>840</td><td>-1.7</td><td>(1)</td><td>...</td></tr><tr><td>40</td><td>New-14</td><td>0</td><td>170.248</td><td>69.207</td><td>-26.38</td><td>18</td><td>4</td><td>28</td><td>-0.95</td><td>(2)</td><td>Compact</td></tr><tr><td>41</td><td>Ylgr</td><td>1</td><td>174.004</td><td>256.584</td><td>17.379</td><td>919</td><td>20</td><td>11,000</td><td>-2.09</td><td>(7)</td><td>...</td></tr><tr><td>42</td><td>Sylgr</td><td>0</td><td>176.261</td><td>273.255</td><td>-60.454</td><td>256</td><td>33</td><td>700</td><td>-2.92</td><td>(11)</td><td>High-velocity scatter</td></tr><tr><td>43</td><td>New-15</td><td>0</td><td>178.134</td><td>66.201</td><td>58.999</td><td>152</td><td>23</td><td>960</td><td>-2.01</td><td>(7)</td><td>...</td></tr><tr><td>44</td><td>Gaia-7</td><td>1</td><td>181.295</td><td>78.993</td><td>-35.727</td><td>282</td><td>15</td><td>1300</td><td>-0.88</td><td>(10)</td><td>...</td></tr><tr><td>45</td><td>Gaia-8</td><td>1</td><td>186.342</td><td>305.604</td><td>-48.129</td><td>453</td><td>18</td><td>3600</td><td>-1.77</td><td>(10)</td><td>...</td></tr><tr><td>46</td><td>New-16</td><td>0</td><td>188.285</td><td>102.418</td><td>5.746</td><td>44</td><td>3</td><td>170</td><td>-1.62</td><td>(2)</td><td>...</td></tr><tr><td>47</td><td>Fjorm</td><td>1</td><td>189.867</td><td>103.086</td><td>19.455</td><td>297</td><td>29</td><td>1700</td><td>-2.23</td><td>(3)</td><td>...</td></tr><tr><td>48</td><td>Orphan</td><td>0</td><td>191.105</td><td>68.052</td><td>-15.708</td><td>867</td><td>247</td><td>130,000</td><td>-1.9</td><td>(12)</td><td>Massive system</td></tr><tr><td>49</td><td>Gaia-1</td><td>1</td><td>192.56</td><td>287.833</td><td>-35.545</td><td>200</td><td>28</td><td>1100</td><td>-1.8</td><td>(7)</td><td>...</td></tr><tr><td>50</td><td>C-23</td><td>1</td><td>195.564</td><td>122.347</td><td>-16.78</td><td>29</td><td>6</td><td>310</td><td>-2.36</td><td>(1)</td><td>...</td></tr><tr><td>51</td><td>LMS-1</td><td>0</td><td>198.23</td><td>43.633</td><td>62.255</td><td>358</td><td>29</td><td>15,000</td><td>-2.09</td><td>(13)</td><td>Massive and wide system</td></tr><tr><td>52</td><td>C-22</td><td>1</td><td>199.904</td><td>56.016</td><td>68.974</td><td>39</td><td>11</td><td>470</td><td>-2.05</td><td>(8)</td><td>...</td></tr><tr><td>53</td><td>GD-1</td><td>1</td><td>200</td><td>214.928</td><td>-29.881</td><td>1468</td><td>323</td><td>14,000</td><td>-2.49</td><td>(1)</td><td>...</td></tr><tr><td>54</td><td>Fimbulthul</td><td>0</td><td>201.697</td><td>304.173</td><td>-22.859</td><td>3724</td><td>29</td><td>16,000</td><td>-1.53</td><td>(3)</td><td>Massive system</td></tr><tr><td>55</td><td>Fimbulthul-S</td><td>0</td><td>201.697</td><td>16.148</td><td>1.199</td><td>1734</td><td>25</td><td>4400</td><td>-1.53</td><td>(10)</td><td>Massive system</td></tr><tr><td>56</td><td>New-17</td><td>0</td><td>210.104</td><td>147.375</td><td>-50.685</td><td>13</td><td>0</td><td>29</td><td>-2.35</td><td>(2)</td><td>...</td></tr><tr><td>57</td><td>NGC5466</td><td>1</td><td>211.364</td><td>161.594</td><td>-50.069</td><td>43</td><td>6</td><td>1900</td><td>-1.98</td><td>(3)</td><td>...</td></tr><tr><td>58</td><td>New-18</td><td>0</td><td>214.041</td><td>302.311</td><td>-1.888</td><td>107</td><td>1</td><td>280</td><td>-0.94</td><td>(2)</td><td>...</td></tr><tr><td>59</td><td>Gaia-6</td><td>1</td><td>214.600</td><td>268.951</td><td>-41.694</td><td>145</td><td>14</td><td>1800</td><td>-1.53</td><td>(1)</td><td>...</td></tr><tr><td>60</td><td>New-19</td><td>0</td><td>221.839</td><td>294.781</td><td>-22.535</td><td>40</td><td>4</td><td>100</td><td>-0.79</td><td>(2)</td><td>Compact, velocities unclear</td></tr><tr><td>61</td><td>Pal-5</td><td>1</td><td>229.022</td><td>320.298</td><td>-54.042</td><td>129</td><td>69</td><td>17,000</td><td>-1.41</td><td>(3)</td><td>...</td></tr><tr><td>62</td><td>M5</td><td>1</td><td>229.638</td><td>325.223</td><td>63.075</td><td>91</td><td>11</td><td>710</td><td>-1.29</td><td>(3)</td><td>...</td></tr></tbody></table>
<table><caption>**Table 3**
(Continued)</caption><thead><tr><th>s<sub>ID</sub><br>(1)</th><th>Name<br>(2)</th><th>Sample<br>(3)</th><th>α<sub>0</sub><br>(4)</th><th>α<sub>pole</sub><br>(5)</th><th>δ<sub>pole</sub><br>(6)</th><th>n<br>(7)</th><th>n<sub>v</sub><br>(8)</th><th>M<sub>*</sub><br>(M<sub>c</sub>)<br>(9)</th><th>⟨[Fe/H]⟩<br>(10)</th><th>References<br>(11)</th><th>Remarks<br>(12)</th></tr></thead><tbody><tr><td>63</td><td>Kshir</td><td>1</td><td>230.196</td><td>229.22</td><td>-21.937</td><td>141</td><td>7</td><td>2200</td><td>-1.83</td><td>(7)</td><td>...</td></tr><tr><td>64</td><td>Svol</td><td>1</td><td>231.711</td><td>28.918</td><td>60.957</td><td>234</td><td>16</td><td>2000</td><td>-1.98</td><td>(1)</td><td>...</td></tr><tr><td>65</td><td>Gaia-9</td><td>1</td><td>236.238</td><td>169.091</td><td>-13.272</td><td>233</td><td>25</td><td>950</td><td>-2.21</td><td>(1)</td><td>...</td></tr><tr><td>66</td><td>Ophiuchus</td><td>0</td><td>242.576</td><td>18.564</td><td>-79.685</td><td>391</td><td>72</td><td>4800</td><td>-1.98</td><td>(3, 14)</td><td>...</td></tr><tr><td>67</td><td>NGC6101</td><td>0</td><td>246.451</td><td>227.979</td><td>16.954</td><td>94</td><td>3</td><td>4800</td><td>-1.98</td><td>(3)</td><td>...</td></tr><tr><td>68</td><td>M92</td><td>0</td><td>259.281</td><td>232.737</td><td>-43.4</td><td>202</td><td>23</td><td>2100</td><td>-2.31</td><td>(3)</td><td>...</td></tr><tr><td>69</td><td>NGC6397</td><td>0</td><td>265.175</td><td>268.821</td><td>35.811</td><td>1207</td><td>17</td><td>2500</td><td>-2.02</td><td>(3)</td><td>...</td></tr><tr><td>70</td><td>Gaia-11</td><td>0</td><td>267.612</td><td>195.778</td><td>-15.134</td><td>84</td><td>3</td><td>1100</td><td>-1.19</td><td>(2)</td><td>Strange CMD, velocities unclear</td></tr><tr><td>71</td><td>Hrid</td><td>0</td><td>280.228</td><td>187.604</td><td>2.034</td><td>666</td><td>29</td><td>2000</td><td>-1.13</td><td>(7)</td><td>...</td></tr><tr><td>72</td><td>C-7</td><td>1</td><td>283.419</td><td>76.064</td><td>-35.517</td><td>239</td><td>16</td><td>1500</td><td>-1.52</td><td>(10)</td><td>Two superposed streams</td></tr><tr><td>73</td><td>New-20</td><td>0</td><td>314.963</td><td>50.439</td><td>18.37</td><td>12</td><td>4</td><td>8</td><td>-1.76</td><td>(15)</td><td>...</td></tr><tr><td>74</td><td>New-21</td><td>0</td><td>316.569</td><td>243.966</td><td>-56.299</td><td>762</td><td>3</td><td>7000</td><td>-0.87</td><td>(2)</td><td>High velocity</td></tr><tr><td>75</td><td>Phlegethon</td><td>1</td><td>322.720</td><td>55.023</td><td>-8.319</td><td>632</td><td>41</td><td>1700</td><td>-2.19</td><td>(8)</td><td>...</td></tr><tr><td>76</td><td>NGC7089</td><td>0</td><td>323.363</td><td>235.973</td><td>72.865</td><td>15</td><td>2</td><td>300</td><td>-1.65</td><td>(3)</td><td>...</td></tr><tr><td>77</td><td>NGC7099</td><td>0</td><td>325.092</td><td>79.501</td><td>-44.028</td><td>54</td><td>3</td><td>680</td><td>-2.27</td><td>(3)</td><td>...</td></tr><tr><td>78</td><td>New-22</td><td>0</td><td>325.847</td><td>211.35</td><td>-34.781</td><td>55</td><td>10</td><td>580</td><td>-0.98</td><td>(2)</td><td>Velocities unclear</td></tr><tr><td>79</td><td>New-23</td><td>0</td><td>329.825</td><td>188.254</td><td>66.646</td><td>44</td><td>6</td><td>120</td><td>-2.25</td><td>(8)</td><td>...</td></tr><tr><td>80</td><td>New-24</td><td>0</td><td>334.183</td><td>59.385</td><td>20.639</td><td>56</td><td>7</td><td>170</td><td>-0.64</td><td>(16)</td><td>Compact, velocities unclear</td></tr><tr><td>81</td><td>New-25</td><td>0</td><td>335.959</td><td>50.194</td><td>23.447</td><td>121</td><td>12</td><td>18,000</td><td>-1.68</td><td>(2)</td><td>Velocities unclear, superposed streams?</td></tr><tr><td>82</td><td>New-26</td><td>0</td><td>341.875</td><td>23.599</td><td>-51.695</td><td>43</td><td>4</td><td>300</td><td>-0.17</td><td>(2)</td><td>Velocities unclear</td></tr><tr><td>83</td><td>New-27</td><td>0</td><td>342.224</td><td>80.278</td><td>-5.92</td><td>37</td><td>8</td><td>410</td><td>-2.08</td><td>(2)</td><td>Velocities unclear</td></tr><tr><td>84</td><td>NGC7492</td><td>0</td><td>347.111</td><td>89.005</td><td>-34.892</td><td>26</td><td>1</td><td>3300</td><td>-1.78</td><td>(3)</td><td>...</td></tr><tr><td>85</td><td>C-19</td><td>1</td><td>354.356</td><td>81.45</td><td>-6.346</td><td>46</td><td>12</td><td>3200</td><td>-3.58</td><td>(1)</td><td>...</td></tr><tr><td>86</td><td>Jhelum</td><td>1</td><td>355.231</td><td>5.299</td><td>37.216</td><td>986</td><td>160</td><td>17,000</td><td>-2.12</td><td>(9)</td><td>...</td></tr><tr><td>87</td><td>Tuc-3</td><td>0</td><td>359.150</td><td>176.704</td><td>-30.396</td><td>102</td><td>34</td><td>12,000</td><td>-2.49</td><td>(17)</td><td>...</td></tr></tbody></table>

**Note.** Columns (1) and (2) list the identification number *s*<sub>ID</sub> and name of the stream; a value of “1” in column (3) marks whether the stream is part of the chosen sample of 29 streams that were used to constrain the Galactic potential. Columns (4)–(6) list, respectively, the zero-point in R.A. and position of the pole of the coordinate system used to derive the *φ*<sub>1</sub> and *φ*<sub>2</sub> stream coordinates. Columns (7) and (8) report the number of stars found in the structure, and the number with measured radial velocity, respectively. Column (9) reports the stellar mass estimate, computed assuming that each stream sample is complete to *G*<sub>0</sub> = 20 mag, and that its stellar population follows a PARSEC model isochrone of age 12.5 Gyr with the metallicity value given in column (10). Column (11) lists the source of the metallicity information. Some brief comments on the nature of the structure are given in column (12). We note that stream #6 (“Kwando,” Grillmair 2017) appears to be the same structure as the “Cetus/Palca” stream (Yuan et al. 2022a). Note also that one star (Gaia ID 3267948604442696448) in stream #8 (“New-2”) is a member of group #4 of Dodd et al. (2023), but closer inspection shows that these two structures are different, with stream #8 showing narrow spatial and kinematic coherence.

**References.** (1) Martin et al. (2022); (2) Martin et al. (2023); (3) Harris (2010); (4) Li et al. (2021); (5) Wan et al. (2020); (6) Grillmair (2017); (7) SDSS; (8) LAMOST; (9) Ji et al. (2020); (10) This work; (11) Roederer & Gnedin (2019); (12) Koposov et al. (2023); (13) Yuan et al. (2020); (14) Caldwell et al. (2020); (15) Gaia BP/RP; (16) Andrae et al. (2023); (17) Li et al. (2019).

were harder to discover simply because they are more diffuse than those found earlier. For instance, streams 8, 17, and 22 are examples of very metal-poor structures ([Fe/H] = -2.07, -2.71, and -2.18, respectively) that subtend a long angle over the sky (62°, 69°, and 80°, respectively) are nearby (≈1.2, ≈1.5, and ≈1.25 kpc distant, respectively) and display a well-defined velocity gradient.

We will discuss the new findings in detail in a future contribution, but here, we nevertheless briefly mention some interesting highlights, either because the new stream is particularly interesting, or because the new data or analysis conveys some new insight into the previously known structure. We will discuss these highlighted streams in order of their *s*<sub>ID</sub> index.

### 9.1. Stream #17 (New-3)

This newly discovered structure has 120 candidate members, including 10 radial velocity measurements and 22 with excellent parallaxes (ω/δω > 10), which are all

within 2 kpc of the Sun. A chemical abundance analysis with the MyGIsFOS code (Sbordone et al. 2014) of the stars observed at high resolution with the UVES instrument will be presented elsewhere. However, we note here that the two well-measured stars (with Gaia IDs: 3089847099636770560, 3074553030332697344) have [Fe/H] metallicities of -2.75 and -2.67, respectively, which places this system close to the most metal-poor known.

Another surprising property of this stream is its orbit. We use the machinery described in Section 5 to fit the phase-space data, using again only the velocity-confirmed members. The resulting progenitor path is displayed in Figure 18 in relation to the Sun (yellow circle), the Galactic Center, and the positions of the RR Lyrae members of the Sagittarius dwarf galaxy in Gaia DR3, which we have selected using the parameter filter described in Ibata et al. (2020a). It therefore seems highly probable that stream #17 formed in the early universe as a globular cluster satellite of the Sagittarius dwarf, and was brought into the Galaxy with it.