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

HAL
open science

# The Global Dynamical Atlas of the Milky Way Mergers: Constraints from Gaia EDR3-based Orbits of Globular Clusters, Stellar Streams, and Satellite Galaxies

Khyati Malhan, Rodrigo A Ibata, Sanjib Sharma, Benoit Famaey, Michele Bellazzini, Raymond Carlberg, Richard D'souza, Zhen Yuan, Nicolas Martin, Guillaume Thomas

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

► To cite this version:

Khyati Malhan, Rodrigo A Ibata, Sanjib Sharma, Benoit Famaey, Michele Bellazzini, et al.. The Global Dynamical Atlas of the Milky Way Mergers: Constraints from Gaia EDR3-based Orbits of Globular Clusters, Stellar Streams, and Satellite Galaxies. The Astrophysical Journal, 2022, 926 (2), pp.107. 10.3847/1538-4357/ac4d2a. hal-03585659

HAL Id: hal-03585659

https://hal.science/hal-03585659

Submitted on 24 Feb 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

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

CrossMark

# The Global Dynamical Atlas of the Milky Way Mergers: Constraints from Gaia EDR3-based Orbits of Globular Clusters, Stellar Streams, and Satellite Galaxies

Khyati Malhan¹,<fnref n="9" /> Ⓓ, Rodrigo A. Ibata² Ⓓ, Sanjib Sharma³,⁴ Ⓓ, Benoit Famaey² Ⓓ, Michele Bellazzini⁵ Ⓓ, Raymond G. Carlberg⁶ Ⓓ, Richard D’Souza⁷ Ⓓ, Zhen Yuan² Ⓓ, Nicolas F. Martin¹,² Ⓓ, and Guillaume F. Thomas⁸ Ⓓ

¹Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117, Heidelberg, Germany; kmalhan07@gmail.com

²Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France

³Sydney Institute for Astronomy, School of Physics, The University of Sydney, NSW 2006, Australia

⁴ARC Centre of Excellence for All Sky Astrophysics in Three Dimensions (ASTRO-3D), Australia

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

⁶Department of Astronomy & Astrophysics, University of Toronto, Toronto, ON M5S 3H4, Canada

⁷Vatican Observatory, Specola Vaticana, V-00120, Vatican City State, Italy

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

Received 2021 December 11; revised 2022 January 14; accepted 2022 January 18; published 2022 February 17

## Abstract

The Milky Way halo was predominantly formed by the merging of numerous progenitor galaxies. However, our knowledge of this process is still incomplete, especially in regard to the total number of mergers, their global dynamical properties and their contribution to the stellar population of the Galactic halo. Here, we uncover the Milky Way mergers by detecting groupings of globular clusters, stellar streams, and satellite galaxies in action (J) space. While actions fully characterize the orbits, we additionally use the redundant information on their energy (E) to enhance the contrast between the groupings. For this endeavor, we use Gaia EDR3-based measurements of 170 globular clusters, 41 streams, and 46 satellites to derive their J and E. To detect groups, we use the ENLINK software, coupled with a statistical procedure that accounts for the observed phase-space uncertainties of these objects. We detect a total of N = 6 groups, including the previously known mergers Sagittarius, Cetus, Gaia-Sausage/Enceladus, LMS-1/Wukong, Arjuna/Sequoia/I’itoi, and one new merger that we call Pontus. All of these mergers, together, comprise 62 objects (≈25% of our sample). We discuss their members, orbital properties, and metallicity distributions. We find that the three most-metal-poor streams of our galaxy—“C-19” ([Fe/H] = -3.4 dex), “Sylgr” ([Fe/H] = -2.9 dex), and “Phoenix” ([Fe/H] = -2.7 dex)—are associated with LMS-1/Wukong, showing it to be the most-metal-poor merger. The global dynamical atlas of Milky Way mergers that we present here provides a present-day reference for galaxy formation models.

*Unified Astronomy Thesaurus concepts: Globular star clusters (656); Milky Way formation (1053); Milky Way stellar halo (1060); Dwarf galaxies (416); Stellar streams (2166); Galaxy formation (595); Galaxy structure (622)*

## 1. Introduction

The stellar halo of the Milky Way was predominantly formed by the merging of numerous progenitor galaxies (Ibata et al. 1994; Helmi et al. 1999; Chiba & Beers 2000; Majewski et al. 2003; Bell et al. 2008; Newberg et al. 2009; Nissen & Schuster 2010; Belokurov et al. 2018; Helmi et al. 2018; Koppelman et al. 2019a; Matsuno et al. 2019; Myeong et al. 2019; Naidu et al. 2020; Yuan et al. 2020b), and this observation appears consistent with the ΛCDM-based models of galaxy formation (e.g., Bullock & Johnston 2005; Pillepich et al. 2018). However, challenging questions remain, for instance: How many progenitor galaxies actually merged with our galaxy? What were the initial physical properties of these merging galaxies, including their stellar and dark matter masses, their stellar population, and their chemical composition (e.g., their [Fe/H] distribution function)? Which objects among the observed population of globular clusters, stellar

streams, and satellite galaxies in the Galactic halo were accreted inside these mergers? Answering these questions is important to understand the hierarchical buildup of our galaxy and thereby to inform galaxy formation models.

It was recently proposed that a significant fraction of the Milky Way’s stellar halo (~95% of the stellar population) resulted from the merging of ≈9–10 progenitor galaxies. This scenario is suggested by Naidu et al. (2020), who identified these mergers by selecting “overdensities” in the chemodynamical space of ~5700 giant stars (these giants lie within 50 kpc from the Galactic center). Many of their selections were based on the knowledge of the previously known mergers. Here, our motivation is also to find the mergers of our galaxy but using a different approach from theirs. First, our objective is to be able to detect these mergers using the data (and not select them) while being agnostic about the previously claimed mergers of our galaxy. Second, we aim for a procedure that is possibly reproducible in cosmological simulations. Finally, we use a very different sample of halo objects, comprising only of globular clusters, stellar streams, and satellite galaxies.

The Milky Way halo harbors a large population of globular clusters (Harris 2010; Vasiliev & Baumgardt 2021), stellar streams (Ibata et al. 2021; Li et al. 2021a), and satellite galaxies (McConnachie & Venn 2020; Battaglia et al. 2022), and these

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

<footnote n="9">⁹ Humboldt Fellow and IAU Gruber Fellow.</footnote>

<footnote>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.</footnote>
objects represent the most ancient and metal-poor structures of our galaxy (e.g., Harris 2010; Kirby et al. 2013; Helmi 2020).<fnref n="10" /> A majority of halo streams are the tidal remnants of either globular clusters or very low-mass satellites (Ibata et al. 2021; Li et al. 2021a, see Section 2.1). For all of these halo objects, a significant fraction of their population is expected to have been brought into the Galactic halo inside massive progenitor galaxies (e.g., Deason et al. 2015; Kruijssen et al. 2019; Carlberg 2020). This implies that these objects, which today form part of our Galactic halo, can be used to trace their progenitor galaxies (e.g., Malhan et al. 2019b; Massari et al. 2019; Bonaca et al. 2021). Consequently, this knowledge also has direct implications on the long-standing question—which of the globular clusters (or streams) were initially formed within the stellar halo (and represent an in situ population) and which were initially formed in different progenitors that only later merged into the Milky Way (and represent an ex situ population). Therefore, using these halo objects as tracers of mergers is also important to understand their own origin and birth sites.

In this regard, using these halo objects also provides a powerful means to detect even the most-metal-poor mergers of the Milky Way. This is important, for instance, to understand the origin of the metal-poor population of the stellar halo (e.g., Komiya et al. 2010; Sestito et al. 2021) and also to constrain the formation scenarios of the metal-deficient globular clusters inside high-redshift galaxies (e.g., Forbes et al. 2018b). In the context of stellar halos, we currently lack knowledge of the origin of the “metallicity floor” for globular clusters; this has recently been pushed down from $[Fe/H] = -2.5$ dex (Harris 2010) to $[Fe/H] = -3.4$ dex (Martin et al. 2022a). In fact, the stellar halo harbors several very metal-poor globular clusters (e.g., Simpson 2018) and also streams (e.g., Roederer & Gnedin 2019; Wan et al. 2020; Martin et al. 2022a). These observations raise the question: Did these metal-poor objects originally form in the Milky Way itself or were they accreted inside the merging galaxies? Moreover, such globular-cluster-merger associations also allow us to understand the cluster formation processes inside those protogalaxies that had formed in the early universe (e.g., Freeman & Bland-Hawthorn 2002; Frebel & Norris 2015).

Our underlying strategy to identify the Milky Way’s mergers is as follows. For each halo object we compute their three actions *J* and then detect those “groups” that clump together tightly in *J* space.<fnref n="11" /> However, to enhance the contrast between groups, we additionally use the redundant information on their energy *E* as this allows us to separate the groups even more confidently.

The motivation behind this strategy can be explained as follows. First, imagine a progenitor galaxy (that is yet to be merged with the Milky Way) containing its own population of globular clusters, satellite galaxies, and streams,<fnref n="12" /> along with its population of stars. Upon merging with the Milky Way, the

progenitor galaxy will get tidally disrupted and deposit its contents into the Galactic halo. If the tidal disruption occurs slowly, the stars of the merging galaxy will themselves form a vast stellar stream in the Galactic halo (e.g., this is the case for the Sagittarius merger; Ibata et al. 2020; Vasiliev & Belokurov 2020). However, if the disruption occurs rapidly, then the stars will quickly get phase-mixed and no clear signature of the stream will be visible (e.g., this is expected for the Gaia-Sausage/Enceladus merger; Belokurov et al. 2018; Helmi et al. 2018). In either case, the member objects of the progenitor galaxy (i.e., its member globular clusters, satellites, and streams), which are now inside the Galactic halo, will possess very similar values of actions *J*. This is because the dynamical quantities *J* are conserved for a very long time, if the potential of the primary galaxy changes adiabatically. The Milky Way’s potential likely evolved adiabatically (e.g., Cardone & Sereno 2005) and, therefore, those objects that merge inside the same progenitor galaxy are expected to remain tightly clumped in the *J* space of the Milky Way, even long after they have been tidally removed from their progenitor. While *E* is not by itself an adiabatic invariant, objects that merge together are expected to occupy a small subset of the energy space. Hence, even though actions fully characterize the orbits, *E* is useful as a redundant “weight” to enhance the contrast between different groups. Moreover, because the mass of the merging galaxies ($M_{\text{halo}} \lesssim 10^{9-11} M_{\odot}$; e.g., Robertson et al. 2005; D’Souza & Bell 2018) are typically much smaller than that of our galaxy ($M_{\text{halo}} \sim 10^{12} M_{\odot}$; e.g., Karukes et al. 2020), the merged objects are expected to occupy only a small volume of the (*J*, *E*) space. Therefore, detecting tightly clumped groups of halo objects in *J* and *E* space potentially provides a powerful means to detect the past mergers that contributed to the Milky Way’s halo.

This strategy for detecting mergers has now become feasible in the era of the ESA/Gaia mission (Gaia Collaboration et al. 2016) because the precision of this astrometric data set allows one to compute reasonably accurate ($J$, $E$) values for a very large population of halo objects. In particular, the excellent Gaia EDR3 data set (Gaia Collaboration et al. 2021; Lindegren et al. 2021) has provided the means to obtain very precise phase-space measurements for an enormously large number of globular clusters (e.g., Vasiliev & Baumgardt 2021), stellar streams (e.g., Ibata et al. 2021; Li et al. 2021a), and satellite galaxies (e.g., McConnachie & Venn 2020), and we use these measurements in the present study.

Before proceeding further, we note that some recent studies have also analyzed energies and angular momenta of globular clusters (e.g., Massari et al. 2019) and streams (e.g., Bonaca et al. 2021). However, the objective of those studies was largely to associate these objects with the previously known mergers of the Milky Way.<fnref n="13" /> Here, our objective is fundamentally different, namely—to detect the Milky Way’s mergers by being completely agnostic about the previously hypothesized groupings of mergers and accretions.

This paper is arranged as follows. In Section 2, we describe the data used for the halo objects and explain our method to compute their actions and energy values. In Section 3, we present our procedure for detecting the mergers by finding “groups of objects” in (*J*, *E*) space. In Section 4, we analyze the detected mergers for their member objects, their dynamical

<footnote n="10">¹⁰ While globular clusters and satellite galaxies represent two very different categories of stellar systems, streams do not represent a third category as they are produced from either globular clusters or satellite galaxies. Streams differ from the other two objects only in terms of their dynamical evolution, in the sense that streams are much more dynamically evolved.</footnote>

<footnote n="11">¹¹ Conceptually, *J* represents the amplitude of an object’s orbit along different directions. For instance, in cylindrical coordinates, $\mathbf{J} = (J_R, J_\phi, J_z)$, where $J_\phi$ represents the z component of angular momentum ($\equiv L_z$), and $J_R$ and $J_z$ describe the extent of oscillations in cylindrical radius and z directions, respectively.</footnote>

<footnote n="12">¹² The population of streams can originate from the tidal stripping of the member globular clusters and/or satellites inside the progenitor galaxy (e.g., Carlberg 2018).</footnote>

<footnote n="13">¹³ An exception is the study by Myeong et al. (2019), who used the Gaia DR2 measurements of globular clusters and hypothesized the “Sequoia” merger.</footnote>
properties, and their [Fe/H] distribution function. In Section 5, we discuss the properties of a specific candidate merger. Additionally, in Section 6, we find several physical connections between streams and other objects (based on the similarity of their orbits and [Fe/H]). Finally, we discuss our findings and conclude in Section 7.

## 2. Computing Actions and Energy of Globular Clusters, Stellar Streams, and Satellite Galaxies

To compute the (J, E) of an object, we require (1) data of the complete 6D phase-space measurements of that object, i.e., its 2D sky position (α, δ), heliocentric distance (D⊙) or parallax (π), 2D proper motion (μ*α* ≡ μα cos δ, μδ), and line-of-sight velocity (vlos), and (2) a Galactic potential model that suitably represents the Milky Way. Below, Section 2.1 describes the phase-space measurements of *n* = 257 objects and Section 2.2 details the adopted Galactic potential model and our procedure for computing the (J, E) quantities.

### 2.1. Data

For globular clusters, we obtain their phase-space measurements from the Vasiliev & Baumgardt (2021) catalog. This catalog provides, for 170 globular clusters, their phase-space measurements, and we use the observed heliocentric coordinates (i.e., α, δ, D⊙, μ*α*, μδ, vlos) along with the associated uncertainties. Vasiliev & Baumgardt (2021) derive the 4D astrometric measurement (α, δ, μ*α*, μδ) of globular clusters using the Gaia EDR3 data set, while the parameters D⊙ and vlos are based on a combination of Gaia EDR3 and other surveys.

For satellite galaxies, we obtain their phase-space measurements from the McConnachie & Venn (2020) catalog. This catalog provides data in heliocentric coordinates format similar to that described above but for the satellite galaxies. From this catalog, we use only those objects that lie within a distance of D⊙ < 250 kpc (equivalent to the virial radius of the Milky Way; Correa Magnus & Vasiliev 2021), yielding a sample of 44 objects. In McConnachie & Venn (2020), the uncertainties on each component of the proper motion are only the observational uncertainties, and therefore, we add in quadrature a systematic uncertainty of 0.033 mas yr⁻¹ to each component of proper motion (A. McConnachie, private communication). While inspecting this catalog, we found that it lacks the proper motion measurements of two other satellites of the Milky Way, namely Bootes III (Grillmair 2009) and the Sagittarius dSph (Ibata et al. 1994). For Bootes III, we obtain its Gaia DR2–based proper motion from Carlin & Sand (2018). For Sagittarius, we use the Vasiliev & Belokurov (2020) catalog that provides Gaia DR2–based proper motions for this dwarf. From this, we compute the median and uncertainty for the Sagittarius dSph as (μ*α*, μδ) = (-2.67 ± 0.45, -1.40 ± 0.40) mas yr⁻¹. Our final sample comprises 46 satellite galaxies.

For stellar streams, we acquire their phase-space measurements primarily from the Ibata et al. (2021) catalog, but we also use some other public stream catalogs (as described below). We first use the Ibata et al. (2021) catalog that contains those streams detected in the Gaia DR2 and EDR3 data sets using the STREAMFINDER algorithm (Malhan & Ibata 2018; Malhan et al. 2018; Ibata et al. 2019b). In this catalog, all the stream stars possess the 5D astrometric measurements on (α, δ, π, μ*α*, μδ), along with their observational uncertainties, as listed in the EDR3 catalog. However, most of these stream stars lack spectroscopic

vlos measurements; this is because (to date) Gaia has provided vlos for only very bright stars with $G \lesssim 12$ mag. Therefore, to obtain the missing vlos measurements, we use various available spectroscopic surveys and also rely on the data from our own follow-up spectroscopic campaigns. These spectroscopic measurements are already presented in Ibata et al. (2021) for the streams “Pal 5” (originally discovered by Odenkirchen et al. 2001), “GD-1” (Grillmair & Dionatos 2006), “Orphan” (Grillmair 2006; Belokurov et al. 2007), “Atlas” (Shipp et al. 2018), “Gaia-1” (Malhan et al. 2018), “Phlegethon” (Ibata et al. 2018), “Slidr” (Ibata et al. 2019b), “Ylgr” (Ibata et al. 2019b), “Leiptr” (Ibata et al. 2019b), “Svöl” (Ibata et al. 2019b), “Gjöll” (Ibata et al. 2019b, the stream of NGC 3201; Hansen et al. 2020; Palau & Miralda-Escude 2021), “Fjörm” (Ibata et al. 2019b, the stream of NGC 4590/M68; Palau & Miralda-Escude 2019), “Sylgr” (Ibata et al. 2019b, the low-metallicity stream with [Fe/H] = -2.92 dex; Roederer & Gnedin 2019), “Fimbulthul” (stream of the ω Centauri cluster, Ibata et al. 2019a), “Kshir” (Malhan et al. 2019a), “M92” (Thomas et al. 2020; Sollima 2020), “Hríd,” “C-7,” “C-3,” “Gunnthrà,” and “NGC 6397.” This spectroscopic campaign suggests that $\ge$85% of the Ibata et al. (2021) sample stars are bona fide stream members.

Ibata et al. (2021) also detected other streams, namely “Indus” (Shipp et al. 2018), “Jhelum” (Shipp et al. 2018), “NGC 5466” (Belokurov et al. 2006; Grillmair & Johnson 2006), “M 5” (Grillmair 2019), “Phoenix” (Balbinot et al. 2016, the low-metallicity globular cluster stream with [Fe/H] = -2.7 dex; Wan et al. 2020), “Gaia-6,” “Gaia-9,” “Gaia-10,” “Gaia-12,” “NGC 7089.” For these streams, we obtain their vlos measurements in this study by cross-matching their stars with various public spectroscopic catalogs, namely SDSS/Segue (Yanny et al. 2009), LAMOST DR7 (Zhao et al. 2012), APOGEE DR16 (Majewski et al. 2017), S5 DR1 (Li et al. 2019), and our own spectroscopic data (that we have collected from our follow-up campaigns; Ibata et al. 2021).

Finally, we include additional streams into our analysis from some of the public stream catalogs. From Malhan et al. (2021) we take the data for the “LMS-1” stream (a recently discovered dwarf galaxy stream, Yuan et al. 2020a). We use the Yuan et al. (2021) catalog for the streams “Palca” (Shipp et al. 2018), “C-20” (Ibata et al. 2021), and “Cetus” (Newberg et al. 2009). For “Ophiuchus” (Bernard et al. 2014), we use those stars provided in the Caldwell et al. (2020) catalog that possess membership probabilities >0.8. We also include those streams provided in the S5 DR1 survey (Li et al. 2019) but were not detected by Ibata et al. (2021). These include “Elqui,” “AliqaUma,” and “Chenab.” Lastly, we also include into our analysis the “C-19” stream (the most-metal-poor globular cluster stream known to date with [Fe/H] ≈ -3.4 dex; Martin et al. 2022a). For all these streams, we use the Gaia EDR3 astrometry.

In summary, we use a total of 41 stellar streams for this study. This stream sample comprises *n* = 9192 Gaia EDR3 stars, of which 1485 possess spectroscopic vlos measurements. The parallaxes of all the stream stars are corrected for the global parallax zero point in Gaia EDR3 using the Lindegren et al. (2021) value.

Figure 1 shows phase-space measurements of all the *n* = 257 objects considered in our study. In this plot, the distance of a given stream corresponds to the inverse of the uncertainty-weighted average mean parallax value of its member stars.
<figure><img src="image_6.png" /><figcaption>**Figure 1.** The Galactic maps showing phase-space measurements of *n* = 257 halo objects used in our study, namely 170 globular clusters (denoted by “star” markers), 41 stellar streams (denoted by “dot” markers), and 46 satellite galaxies (denoted by “square” markers). From panel (a) to (d), these objects are colored by their heliocentric distances (*D*<sub>☉</sub>), line-of-sight velocities (v<sub>los</sub>), proper motion in the Galactic longitude direction (μ<sub>ℓ</sub><sup>*</sup>), and proper motion in the Galactic latitude direction (μ<sub>b</sub>), respectively. In panel (d), we only show streams with their names and do not plot other objects to avoid crowding.</figcaption></figure>

## 2.2. Computing Actions and Energy of the Halo Objects

To compute the orbits of the objects in our sample, we adopt the Galactic potential model of McMillan (2017). This is a static and axisymmetric model comprising a bulge, disk components, and an NFW (Navarro–Frenk–White) halo. For this potential model, the total galactic mass within the galactocentric distance $r_{\text{gal}} < 20$ kpc is $2.2 \times 10^{11} M_{\odot}$, $r_{\text{gal}} < 50$ kpc is $4.9 \times 10^{11} M_{\odot}$, and $r_{\text{gal}} < 100$ kpc is $8.1 \times 10^{11} M_{\odot}$. Another model that is often used to represent the Galactic potential is MWPotential2014 of Bovy (2015), and this model (on average) is ~1.5 times lighter than the McMillan (2017) model. For our study, we prefer the McMillan (2017) model because (1) the predicted velocity curve of this model is more consistent with the measurements of the Milky Way (e.g., Bovy 2020; Nitschai et al. 2021), and (2) we find that all halo objects in this mass model possess $E < 0$ (i.e., their orbits are bound); however, in the case of MWPotential2014 we infer that 34 clusters and all the satellite galaxies possess $E \ge 0$ (i.e., their orbits are unbound). To set the McMillan (2017) potential model and to compute (J, E) and other orbital parameters, we make use of the galpy module (Bovy 2015). Moreover, to transform the heliocentric phase-space measurements of the objects into the Galactocentric frame (which is required for computing orbits), we adopt the Sun’s Galactocentric distance

from Gravity Collaboration et al. (2018) and the Sun’s galactic velocity from Reid et al. (2014) and Schonrich et al. (2010).

To compute the (J, E) values of globular clusters, we do the following. For a given globular cluster, we sample 1000 orbits using the mean and the uncertainty on its phase-space measurement. For that particular cluster, this provides an (J, E) distribution of 1000 data points, and this distribution represents the uncertainty in the derived (J, E) value for that cluster. Note that this (J, E) uncertainty, for a given cluster, reflects its uncertainty on the phase-space measurement. This orbit-sampling procedure is repeated for all globular clusters, and for each cluster we retain their respective (J, E) distribution. This (J, E) distribution is a vital information, and we subsequently use this while detecting the mergers (as shown in Section 3). The resulting (J, E) distribution of all the globular clusters is shown in Figure 2, where each object is effectively represented by a distribution of 1000 points.

We analyze actions in cylindrical coordinates, i.e., in the $J \equiv (J_R, J_\phi, J_z)$ system, where $J_\phi$ corresponds to the z component of angular momentum (i.e., $J_\phi \equiv L_z$) and negative $J_\phi$ represents prograde motion (i.e., rotational motion in the direction of the Galactic disk). Similarly, components $J_R$ and $J_z$ describe the extent of oscillations in cylindrical radius and z directions, respectively. Figure 2 shows these globular clusters in (1) the “projected action space,” represented by a diagram of $J_\phi/J_{\text{total}}$
<figure><img src="image_7.png" /></figure>

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

<figure><img src="image_9.png" /><figcaption>**Figure 2.** Action-energy (J, E) space of the Milky Way showing the globular clusters (top panels), stellar streams (middle panels), and satellite galaxies (bottom panels). Each object can be seen as a “cloud” of 1000 Monte Carlo representations of its orbit (see Section 2.2). In each row, the left panel corresponds to the projected action-space map, where the horizontal axis is $J_\phi/J_{\text{tot}}$ and the vertical axis is $(J_z - J_R)/J_{\text{tot}}$ with $J_{\text{tot}} = J_R + J_z + |J_\phi|$. In these panels, the points are colored by the total energy of their orbits (E). The right panels show the z component of the angular momentum ($J_\phi \equiv L_z$) vs. E, and the points are colored by the orthogonal component of their angular momenta ($L_\perp = \sqrt{L_x^2 + L_y^2}$).</figcaption></figure>
versus $(J_z - J_R)/J_{total}$, where $J_{total} = J_R + J_z + |J_\phi|$, and (2) the $J_\phi$ versus E space. The reason for using the projected action space is that this plot is effective in separating objects that lie along circular, radial, and in-plane orbits, and it is considered to be superior to other commonly used kinematic spaces (e.g., Lane et al. 2022). We also use the orthogonal component of the angular momentum $L_\perp = \sqrt{L_x^2 + L_y^2}$ for representation. Note that even though $L_\perp$ is not fully conserved in an axisymmetric potential, it still serves as a useful quantity for orbital characterization (e.g., Bonaca et al. 2021). Along with retrieving the (J, E) values, we also retrieve other orbital parameters (e.g., $r_{apo}$, $r_{peri}$, eccentricity—these values are used at a later stage for the analysis of the detected mergers).

To compute the (J, E) values of satellite galaxies, we use exactly the same orbit-sampling procedure as described above for globular clusters. The corresponding (J, E) distribution is shown in Figure 2.

To compute the (J, E) values of stellar streams, we follow the orbit-fitting procedure; this approach is more sophisticated than the above-described orbit-sampling procedure and more suitable for stellar streams. That is, we obtain (J, E) solutions of a given stream by fitting orbits to the phase-space measurements of all its member stars (e.g., Koposov et al. 2010). This procedure ensures that the resulting orbit solution provides a reasonable representation of the entire stream structure and also that the resulting (J, E) values are precise.<fnref n="14" /> We use this method only for narrow and dynamically cold streams (that make up most of our stream sample), but for the other broad and dynamically hot streams we rely on the orbit-sampling procedure (see further below). To carry out the orbit fitting of streams, we follow the same procedure as described in Malhan et al. (2021). Briefly, we survey the parameter space using our own Metropolis–Hastings-based Markov Chain Monte Carlo (MCMC) algorithm, where the log-likelihood of each member star *i* is defined as

$$ \ln \mathcal{L}_i = -\ln((2\pi)^{5/2} \sigma_{\text{sky}} \sigma_{\varpi} \sigma_{\mu_\alpha} \sigma_{\mu_\delta} \sigma_{v_{\text{los}}}) + \ln N - \ln D, \quad (1) $$

where

$$ \begin{aligned} N &= \prod_{j=1}^{5} (1 - e^{-R_j^2/2}), \qquad D = \prod_{j=1}^{5} R_j^2, \\ R_1^2 &= \frac{\theta_{\text{sky}}^2}{\sigma_{\text{sky}}^2}, \qquad R_2^2 = \frac{(\varpi^d - \varpi^o)^2}{\sigma_{\varpi}^2}, \\ R_3^2 &= \frac{(\mu_\alpha^d - \mu_\alpha^o)^2}{\sigma_{\mu_\alpha}^2}, \qquad R_4^2 = \frac{(\mu_\delta^d - \mu_\delta^o)^2}{\sigma_{\mu_\delta}^2}, \\ R_5^2 &= \frac{(v_{\text{los}}^d - v_{\text{los}}^o)^2}{\sigma_{v_{\text{los}}}^2}. \end{aligned} \quad (2) $$

Here, $\theta_{\text{sky}}$ is the on-sky angular difference between the orbit and the data point, $\varpi^d$, $\mu_\alpha^d$, $\mu_\delta^d$ and $v_{\text{los}}^d$ are the measured data parallax, proper motion, and line-of-sight velocity, with the corresponding orbital model values marked with “o.” The Gaussian dispersions $\sigma_{\text{sky}}$, $\sigma_{\varpi}$, $\sigma_{\mu_\alpha}$, $\sigma_{\mu_\delta}$, $\sigma_{v_{\text{los}}}$ are the sum in quadrature of the intrinsic dispersion of the model and the observational uncertainty of each data point. The particular

reason for adopting this “conservative formulation” of the log-likelihood function (Sivia 1996) is to lower the contribution from outliers that could be contaminating the stream data. Furthermore, in a given stream, we set all those stars that lack spectroscopic measurements to $v_{\text{los}} = 0 \text{ km s}^{-1}$ with a $10^4 \text{ km s}^{-1}$ Gaussian uncertainty. While undertaking this orbit-fitting procedure for a given stream, we chose to anchor the orbit solutions at a fixed R.A. value (which was approximately halfway along the stream), while leaving all the other parameters to be varied. We do this because without setting an anchor, the solution would have wandered over the full length of the stream. The success of such a procedure in fitting streams has been demonstrated before (Malhan & Ibata 2019; Malhan et al. 2021). This procedure works well for most of the streams, as the final MCMC chains are converged and the resulting best-fit orbits provide good representations to the phase-space structures of all these streams.

The above orbit-fitting procedure was carried out for the majority of streams; however, for a subset of them we considered it better to instead adopt the orbit-sampling procedure. This subset includes LMS-1, Orphan, Fimbulthul, Cetus, Svol, NGC 6397, Ophiuchus, C-3, Gaia-6, and Chenab. The orbit-sampling procedure means that we no longer use Equation (1) (this ensures that the resulting orbit provides a reasonable fit to the entire stream structure), but instead, we simply sample orbits using directly the phase-space measurements of the individual member stars (this does not guarantee an orbit-fit to the entire stream structure). The reason for adopting this scheme for LMS-1, Cetus (which are dwarf galaxy streams), and Fimbulthul (which is the stream of the massive $\omega$ Cen cluster) is that these are dynamically hot and physically broad streams, and the aforementioned orbit-fitting procedure would have underestimated their dispersions in the derived (J, E) quantities. Similarly, Ophiuchus also appears to possess a broader dispersion in $v_{\text{los}}$ space (~10–15 km s⁻¹; see Figure 10 of Caldwell et al. 2020). For Orphan, which is a stream with a “twisted” shape (due to perturbation by the LMC; Erkal et al. 2019), we deemed it better to sample its orbits (Li et al. 2021a also adopt a similar procedure to compute the orbit of Orphan). For the remaining streams, although they did appear narrow and linear in the $(\alpha, \delta)$ and $(\mu_\alpha^*, \mu_\delta)$ spaces, it was difficult to visualize this linearity in $v_{\text{los}}$ space. This was primarily because these streams lack enough spectroscopic measurements so that a clear stream signal can be visible in $v_{\text{los}}$ space. Therefore, it was difficult to apply the orbit-fitting procedure for them, and we resort to the orbit-sampling procedure. For all of these streams, the sampling in $\alpha, \delta, \mu_\alpha^*, \mu_\delta$ and $v_{\text{los}}$ was performed directly using the measurements and the associated uncertainties. However, to sample over the distance parameter in a given stream, we computed the average distance (and the uncertainty) using the uncertainty-weighted average mean parallax of the member stars.

The above orbit-fitting and orbit-sampling schemes generate the MCMC chains for the orbital parameters of all 41 streams, and for each stream, we randomly sample 1000 steps (this we do after rejecting the burn-in phase). These sampled values are shown in Figure 2. Note that for most of the streams, their (J, E) dispersions are much smaller than those of globular clusters and satellite galaxies. This is because the orbits of streams are much more precisely constrained (because we employ the above orbit-fitting procedure). The derived orbital properties of our streams are provided in Tables 1 and 2.

<footnote n="14">¹⁴ By employing the orbit-fitting procedure, we are assuming that the entire phase-space structure of a stream can be well represented by an orbit. Although streams do not strictly delineate orbits (Sanders & Binney 2013), our assumption is still reasonable as far as the scope of this study is concerned.</footnote>
<table><caption>**Table 1**

Constrained Heliocentric Parameters of Stellar Streams. For Each Stream, the Following Values Represent the Posterior Distribution at the Stream's "Anchor" Point (i.e., at a Fixed R.A. Value)</caption><thead><tr><th>Stream</th><th>No. of Gaia Sources</th><th>No. of v<sub>los</sub> sources</th><th>R.A. (deg)</th><th>Decl. (deg)</th><th>D<sub>☉</sub> (kpc)</th><th>μ<sup>*</sup><sub>☉</sub> (mas yr<sup>-1</sup>)</th><th>μ<sub>δ</sub> (mas yr<sup>-1</sup>)</th><th>v<sub>los</sub> (km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Gjoll</td><td>102</td><td>35</td><td>82.1</td><td>-13.95<sup>+0.35</sup><sub>-0.36</sub></td><td>3.26<sup>+0.03</sup><sub>-0.03</sub></td><td>23.58<sup>+0.08</sup><sub>-0.09</sub></td><td>-23.7<sup>+0.06</sup><sub>-0.05</sub></td><td>78.73<sup>+2.36</sup><sub>-1.84</sub></td></tr><tr><td>Leiptr</td><td>237</td><td>67</td><td>89.11</td><td>-28.37<sup>+0.2</sup><sub>-0.27</sub></td><td>7.39<sup>+0.07</sup><sub>-0.07</sub></td><td>10.59<sup>+0.03</sup><sub>-0.04</sub></td><td>-9.9<sup>+0.03</sup><sub>-0.04</sub></td><td>194.22<sup>+2.23</sup><sub>-1.86</sub></td></tr><tr><td>Hrid</td><td>233</td><td>24</td><td>280.51</td><td>33.3<sup>+0.75</sup><sub>-0.6</sub></td><td>2.75<sup>+0.1</sup><sub>-0.07</sub></td><td>-5.88<sup>+0.11</sup><sub>-0.08</sub></td><td>20.08<sup>+0.21</sup><sub>-0.19</sub></td><td>-238.77<sup>+3.3</sup><sub>-5.52</sub></td></tr><tr><td>Pal5</td><td>48</td><td>29</td><td>229.65</td><td>0.26<sup>+0.1</sup><sub>-0.13</sub></td><td>20.16<sup>+0.24</sup><sub>-0.33</sub></td><td>-2.75<sup>+0.03</sup><sub>-0.02</sub></td><td>-2.68<sup>+0.02</sup><sub>-0.02</sub></td><td>-57.03<sup>+1.08</sup><sub>-1.04</sub></td></tr><tr><td>Gaia-1</td><td>106</td><td>8</td><td>190.96</td><td>-9.16<sup>+0.15</sup><sub>-0.1</sub></td><td>5.57<sup>+0.16</sup><sub>-0.1</sub></td><td>-14.39<sup>+0.04</sup><sub>-0.04</sub></td><td>-19.72<sup>+0.03</sup><sub>-0.04</sub></td><td>214.91<sup>+3.5</sup><sub>-2.16</sub></td></tr><tr><td>Ylgr</td><td>699</td><td>32</td><td>173.82</td><td>-22.31<sup>+0.22</sup><sub>-0.3</sub></td><td>9.72<sup>+0.16</sup><sub>-0.14</sub></td><td>-0.44<sup>+0.04</sup><sub>-0.03</sub></td><td>-7.65<sup>+0.05</sup><sub>-0.04</sub></td><td>317.86<sup>+2.83</sup><sub>-3.05</sub></td></tr><tr><td>Fjorm</td><td>182</td><td>28</td><td>251.89</td><td>65.38<sup>+0.24</sup><sub>-0.22</sub></td><td>6.42<sup>+0.16</sup><sub>-0.14</sub></td><td>3.92<sup>+0.07</sup><sub>-0.08</sub></td><td>3.1<sup>+0.06</sup><sub>-0.06</sub></td><td>-25.37<sup>+1.89</sup><sub>-2.19</sub></td></tr><tr><td>Kshir</td><td>55</td><td>16</td><td>205.88</td><td>67.25<sup>+0.13</sup><sub>-0.17</sub></td><td>9.57<sup>+0.08</sup><sub>-0.08</sub></td><td>-7.67<sup>+0.04</sup><sub>-0.04</sub></td><td>-3.92<sup>+0.04</sup><sub>-0.05</sub></td><td>-249.88<sup>+2.62</sup><sub>-2.92</sub></td></tr><tr><td>Gunnthra</td><td>61</td><td>8</td><td>284.22</td><td>-73.49<sup>+0.23</sup><sub>-0.14</sub></td><td>2.83<sup>+0.12</sup><sub>-0.13</sub></td><td>-15.83<sup>+0.11</sup><sub>-0.13</sub></td><td>-24.04<sup>+0.15</sup><sub>-0.17</sub></td><td>132.26<sup>+6.23</sup><sub>-4.97</sub></td></tr><tr><td>Slidr</td><td>181</td><td>29</td><td>160.05</td><td>10.22<sup>+0.43</sup><sub>-0.41</sub></td><td>2.99<sup>+0.11</sup><sub>-0.09</sub></td><td>-24.6<sup>+0.08</sup><sub>-0.08</sub></td><td>-6.65<sup>+0.06</sup><sub>-0.06</sub></td><td>-87.98<sup>+3.44</sup><sub>-3.17</sub></td></tr><tr><td>M92</td><td>84</td><td>9</td><td>259.89</td><td>43.08<sup>+0.2</sup><sub>-0.2</sub></td><td>8.94<sup>+0.2</sup><sub>-0.18</sub></td><td>-5.15<sup>+0.05</sup><sub>-0.05</sub></td><td>-0.63<sup>+0.06</sup><sub>-0.04</sub></td><td>-140.66<sup>+6.28</sup><sub>-7.53</sub></td></tr><tr><td>NGC 3201</td><td>388</td><td>4</td><td>152.46</td><td>-46.32<sup>+0.11</sup><sub>-0.08</sub></td><td>4.99<sup>+0.01</sup><sub>-0.02</sub></td><td>8.87<sup>+0.02</sup><sub>-0.02</sub></td><td>-2.22<sup>+0.02</sup><sub>-0.02</sub></td><td>489.63<sup>+3.36</sup><sub>-3.82</sub></td></tr><tr><td>Atlas</td><td>46</td><td>10</td><td>25.04</td><td>-29.81<sup>+0.1</sup><sub>-0.1</sub></td><td>19.93<sup>+0.76</sup><sub>-0.75</sub></td><td>0.04<sup>+0.02</sup><sub>-0.02</sub></td><td>-0.89<sup>+0.02</sup><sub>-0.02</sub></td><td>-85.65<sup>+1.48</sup><sub>-1.58</sub></td></tr><tr><td>C-7</td><td>120</td><td>10</td><td>287.15</td><td>-50.17<sup>+0.16</sup><sub>-0.14</sub></td><td>6.77<sup>+0.28</sup><sub>-0.21</sub></td><td>-13.79<sup>+0.07</sup><sub>-0.06</sub></td><td>-12.38<sup>+0.06</sup><sub>-0.07</sub></td><td>55.05<sup>+1.5</sup><sub>-2.51</sub></td></tr><tr><td>Palca</td><td>24</td><td>24</td><td>36.57</td><td>-36.15<sup>+0.33</sup><sub>-0.31</sub></td><td>12.31<sup>+1.68</sup><sub>-1.44</sub></td><td>0.9<sup>+0.02</sup><sub>-0.02</sub></td><td>-0.23<sup>+0.04</sup><sub>-0.04</sub></td><td>106.32<sup>+2.62</sup><sub>-2.54</sub></td></tr><tr><td>Sylgr</td><td>165</td><td>19</td><td>179.68</td><td>-2.44<sup>+0.27</sup><sub>-0.4</sub></td><td>3.77<sup>+0.07</sup><sub>-0.11</sub></td><td>-13.98<sup>+0.12</sup><sub>-0.14</sub></td><td>-12.9<sup>+0.09</sup><sub>-0.1</sub></td><td>-184.8<sup>+1.548</sup><sub>-8.15</sub></td></tr><tr><td>Gaia-9</td><td>286</td><td>15</td><td>233.27</td><td>60.42<sup>+0.04</sup><sub>-0.11</sub></td><td>4.68<sup>+0.08</sup><sub>-0.09</sub></td><td>-12.49<sup>+0.11</sup><sub>-0.12</sub></td><td>6.37<sup>+0.14</sup><sub>-0.08</sub></td><td>-359.86<sup>+4.65</sup><sub>-4.11</sub></td></tr><tr><td>Gaia-10</td><td>90</td><td>9</td><td>161.47</td><td>15.17<sup>+0.14</sup><sub>-0.14</sub></td><td>13.32<sup>+0.34</sup><sub>-0.28</sub></td><td>-4.14<sup>+0.05</sup><sub>-0.05</sub></td><td>-3.15<sup>+0.04</sup><sub>-0.04</sub></td><td>289.64<sup>+2.75</sup><sub>-3.32</sub></td></tr><tr><td>Gaia-12</td><td>38</td><td>1</td><td>41.05</td><td>16.45<sup>+0.13</sup><sub>-0.13</sub></td><td>15.71<sup>+1.29</sup><sub>-1.03</sub></td><td>5.84<sup>+0.05</sup><sub>-0.05</sub></td><td>-5.66<sup>+0.07</sup><sub>-0.06</sub></td><td>-303.83<sup>+22.55</sup><sub>-15.07</sub></td></tr><tr><td>Indus</td><td>454</td><td>45</td><td>340.12</td><td>-60.58<sup>+0.1</sup><sub>-0.1</sub></td><td>14.96<sup>+0.19</sup><sub>-0.16</sub></td><td>3.59<sup>+0.03</sup><sub>-0.03</sub></td><td>-4.89<sup>+0.02</sup><sub>-0.03</sub></td><td>-49.15<sup>+2.45</sup><sub>-3.68</sub></td></tr><tr><td>Jhelum</td><td>972</td><td>246</td><td>351.95</td><td>-51.74<sup>+0.08</sup><sub>-0.08</sub></td><td>11.39<sup>+0.13</sup><sub>-0.15</sub></td><td>7.23<sup>+0.04</sup><sub>-0.04</sub></td><td>-4.37<sup>+0.04</sup><sub>-0.03</sub></td><td>-1.29<sup>+2.6</sup><sub>-3.11</sub></td></tr><tr><td>Phoenix</td><td>35</td><td>19</td><td>23.96</td><td>-50.01<sup>+0.24</sup><sub>-0.24</sub></td><td>16.8<sup>+0.33</sup><sub>-0.36</sub></td><td>2.72<sup>+0.03</sup><sub>-0.03</sub></td><td>-0.07<sup>+0.03</sup><sub>-0.03</sub></td><td>45.92<sup>+1.63</sup><sub>-1.58</sub></td></tr><tr><td>NGC5466</td><td>62</td><td>4</td><td>214.41</td><td>26.84<sup>+0.12</sup><sub>-0.11</sub></td><td>14.09<sup>+0.27</sup><sub>-0.25</sub></td><td>-5.64<sup>+0.03</sup><sub>-0.03</sub></td><td>-0.72<sup>+0.03</sup><sub>-0.02</sub></td><td>95.04<sup>+7.4</sup><sub>-5.91</sub></td></tr><tr><td>M5</td><td>139</td><td>5</td><td>206.96</td><td>13.5<sup>+0.15</sup><sub>-0.14</sub></td><td>7.44<sup>+0.12</sup><sub>-0.11</sub></td><td>3.5<sup>+0.03</sup><sub>-0.04</sub></td><td>-8.76<sup>+0.04</sup><sub>-0.04</sub></td><td>-42.97<sup>+3.33</sup><sub>-3.83</sub></td></tr><tr><td>C-20</td><td>34</td><td>9</td><td>359.81</td><td>8.63<sup>+0.16</sup><sub>-0.16</sub></td><td>18.11<sup>+1.45</sup><sub>-1.39</sub></td><td>-0.58<sup>+0.03</sup><sub>-0.03</sub></td><td>1.44<sup>+0.02</sup><sub>-0.02</sub></td><td>-116.87<sup>+1.46</sup><sub>-1.44</sub></td></tr><tr><td>C-19</td><td>34</td><td>8</td><td>355.28</td><td>28.82<sup>+0.63</sup><sub>-1.17</sub></td><td>18.04<sup>+0.55</sup><sub>-0.53</sub></td><td>1.25<sup>+0.03</sup><sub>-0.03</sub></td><td>-2.74<sup>+0.03</sup><sub>-0.05</sub></td><td>-193.48<sup>+2.61</sup><sub>-2.52</sub></td></tr><tr><td>Elqui</td><td>4</td><td>4</td><td>19.77</td><td>-42.36<sup>+0.3</sup><sub>-0.29</sub></td><td>51.41<sup>+4.64</sup><sub>-7.04</sub></td><td>0.33<sup>+0.02</sup><sub>-0.03</sub></td><td>-0.49<sup>+0.02</sup><sub>-0.02</sub></td><td>15.86<sup>+8.82</sup><sub>-20.58</sub></td></tr><tr><td>AliqaUma</td><td>5</td><td>5</td><td>34.08</td><td>-33.97<sup>+0.31</sup><sub>-0.34</sub></td><td>21.48<sup>+2.32</sup><sub>-1.2</sub></td><td>0.24<sup>+0.02</sup><sub>-0.03</sub></td><td>-0.79<sup>+0.03</sup><sub>-0.03</sub></td><td>-42.33<sup>+2.29</sup><sub>-2.23</sub></td></tr><tr><td>Phlegethon</td><td>365</td><td>41</td><td>319.89</td><td>-32.07<sup>+0.43</sup><sub>-0.37</sub></td><td>3.29<sup>+0.05</sup><sub>-0.05</sub></td><td>-3.97<sup>+0.09</sup><sub>-0.09</sub></td><td>-37.66<sup>+0.08</sup><sub>-0.09</sub></td><td>15.9<sup>+4.97</sup><sub>-6.12</sub></td></tr><tr><td>GD-1</td><td>811</td><td>216</td><td>160.02</td><td>45.9<sup>+0.25</sup><sub>-0.19</sub></td><td>8.06<sup>+0.07</sup><sub>-0.07</sub></td><td>-6.75<sup>+0.04</sup><sub>-0.03</sub></td><td>-10.88<sup>+0.04</sup><sub>-0.05</sub></td><td>-101.83<sup>+2.05</sup><sub>-2.47</sub></td></tr></tbody></table>

Note. This anchor is defined during the orbit-fitting procedure. From left to right, the columns provide the stream's name, number of Gaia EDR3 sources in the stream, number of sources with spectroscopic line-of-sight velocities ($v_{los}$), R.A. (which acts as the anchor point in our orbit-fitting procedure), decl., heliocentric distance ($D_{☉}$), proper motions ($\mu_{\phi}^*$, $\mu_{\delta}^*$), and $v_{los}$. The quoted values are medians of the sampled posterior distributions and the corresponding uncertainties represent their 16th and 84th percentiles. Only those streams for which the orbit-fitting procedure was employed are listed (see Section 2.2).

## 2.3. A Qualitative Analysis of the Orbits

As a passing analysis, we qualitatively examine some basic orbital properties of globular clusters, satellite galaxies, and stellar streams. The knowledge gained from this analysis allows us to put our final results in some context.

For globular clusters, we find that ~70% of them move along prograde orbits (i.e., their $J_\phi < 0$), 18% move along polar orbits (i.e., their orbital planes are inclined almost perpendicularly to the Galactic disk plane, with $\phi \ge 75^\circ$), ~12% have their orbits nearly confined to the Galactic plane (i.e., their $\phi \le 20^\circ$), 11% have disk-like orbits (i.e., both prograde and in plane), 1% have in-plane and retrograde orbits, and 10% have highly eccentric orbits (with ecc > 0.8). This excess of prograde globular clusters could be indicating that the Galactic halo itself initially had an excess of prograde clusters or it may owe to the possible spinning of the dark matter halo (e.g., Obreja et al. 2022).

For satellites, we find that ~60% of them move along prograde orbits, ~10% have highly eccentric orbits, and ~45% move along polar orbits (most of these “polar” satellites belong

to the “Vast Plane of Satellites” structure; see Pawlowski et al. 2021). None of the satellites move in the disk plane; this could be because satellites on coplanar orbits are expected to be destroyed quickly compared to those on polar orbits (e.g., Penarrubia et al. 2002). The satellites possess quite high energies and angular momenta compared to the globular clusters (and also stellar streams, as we note below). The high $E$ values of satellites suggest that many of them are not ancient inhabitants of the Milky Way but have only recently arrived into our galaxy (perhaps ~4 Gyr ago; e.g., Hammer et al. 2021).

For stellar streams, we find that 55% of them move along prograde orbits, 22% move along polar orbits, and 5% possess highly eccentric orbits. Some of these polar streams are LMS-1, C-19, Sylgr, Jhelum, Elqui, Gaia-10, Ophiuchus, and Hrid. None of the streams orbit in the disk plane. Our inference on the prograde distribution of streams is somewhat consistent with the study of Panithanpaisal et al. (2021), who analyzed FIRE 2 cosmological simulations and found that Milky Way–mass galaxies should have an even distribution of streams on prograde and retrograde orbits.
<table><caption>**Table 2**
Actions, Energies, Orbital Parameters, and Metallicities of the Stellar Streams</caption><thead><tr><th>Stream</th><th>(J<sub>R</sub>, J<sub>0</sub>, J<sub>2</sub>)<br>(kpc km s<sup>-1</sup>)</th><th>Energy<br>(×10<sup>2</sup> km s<sup>-2</sup>)</th><th>r<sub>peri</sub><br>(kpc)</th><th>r<sub>apo</sub><br>(kpc)</th><th>z<sub>max</sub><br>(kpc)</th><th>ecc.</th><th>[Fe/H]<br>(dex)</th></tr></thead><tbody><tr><td>LMS-1</td><td>(255<sup>+239</sup><sub>-149</sub>, -627<sup>+183</sup><sub>-232</sub>, 2514<sup>+383</sup><sub>-263</sub>)</td><td>-1227<sup>+65</sup><sub>-39</sub></td><td>10.8<sup>+2.5</sup><sub>-1.8</sub></td><td>20.6<sup>+3.7</sup><sub>-1.9</sub></td><td>20.2<sup>+3.6</sup><sub>-2.0</sub></td><td>0.32<sup>+0.11</sup><sub>-0.12</sub></td><td>-2.1 ± 0.4</td></tr><tr><td>Gjoll</td><td>(783<sup>+73</sup><sub>-65</sub>, 2782<sup>+60</sup><sub>-60</sub>, 274<sup>+7</sup><sub>-6</sub>)</td><td>-1152<sup>+19</sup><sub>-18</sub></td><td>8.5<sup>+0.1</sup><sub>-0.1</sub></td><td>27.4<sup>+1.4</sup><sub>-1.2</sub></td><td>10.8<sup>+0.5</sup><sub>-0.5</sub></td><td>0.52<sup>+0.01</sup><sub>-0.01</sub></td><td>-1.78</td></tr><tr><td>Leiptr</td><td>(1455<sup>+133</sup><sub>-119</sub>, 4128<sup>+77</sup><sub>-73</sub>, 378<sup>+11</sup><sub>-10</sub>)</td><td>-933<sup>+19</sup><sub>-18</sub></td><td>12.3<sup>+0.1</sup><sub>-0.1</sub></td><td>45.1<sup>+2.3</sup><sub>-2.1</sub></td><td>17.6<sup>+0.9</sup><sub>-0.9</sub></td><td>0.57<sup>+0.01</sup><sub>-0.01</sub></td><td>...</td></tr><tr><td>Hrid</td><td>(1642<sup>+110</sup><sub>-79</sub>, 78<sup>+54</sup><sub>-47</sub>, 83<sup>+6</sup><sub>-5</sub>)</td><td>-1319<sup>+30</sup><sub>-22</sub></td><td>1.1<sup>+0.0</sup><sub>-0.0</sub></td><td>22.0<sup>+1.4</sup><sub>-1.0</sub></td><td>7.0<sup>+0.9</sup><sub>-0.5</sub></td><td>0.9<sup>+0.0</sup><sub>-0.0</sub></td><td>-1.1</td></tr><tr><td>Pal5</td><td>(282<sup>+19</sup><sub>-17</sub>, -744<sup>+61</sup><sub>-44</sub>, 1357<sup>+42</sup><sub>-31</sub>)</td><td>-1385<sup>+11</sup><sub>-15</sub></td><td>6.9<sup>+0.3</sup><sub>-0.4</sub></td><td>15.8<sup>+0.2</sup><sub>-0.3</sub></td><td>14.7<sup>+0.2</sup><sub>-0.3</sub></td><td>0.39<sup>+0.02</sup><sub>-0.02</sub></td><td>-1.35 ± 0.06</td></tr><tr><td>Orphan</td><td>(959<sup>+978</sup><sub>-271</sub>, -3885<sup>+405</sup><sub>-1017</sub>, 1199<sup>+484</sup><sub>-213</sub>)</td><td>-949<sup>+175</sup><sub>-64</sub></td><td>15.6<sup>+3.8</sup><sub>-2.1</sub></td><td>41.2<sup>+23.6</sup><sub>-6.2</sub></td><td>26.4<sup>+16.9</sup><sub>-4.8</sub></td><td>0.48<sup>+0.09</sup><sub>-0.06</sub></td><td>-1.85 ± 0.53</td></tr><tr><td>Gaia-1</td><td>(3638<sup>+1307</sup><sub>-632</sub>, 2678<sup>+89</sup><sub>-57</sub>, 997<sup>+93</sup><sub>-58</sub>)</td><td>-794<sup>+90</sup><sub>-55</sub></td><td>8.2<sup>+0.1</sup><sub>-0.0</sub></td><td>67.6<sup>+18.3</sup><sub>-8.9</sub></td><td>45.7<sup>+13.6</sup><sub>-6.5</sub></td><td>0.78<sup>+0.04</sup><sub>-0.03</sub></td><td>-1.36</td></tr><tr><td>Fimbulthul</td><td>(202<sup>+109</sup><sub>-78</sub>, 427<sup>+244</sup><sub>-588</sub>, 95<sup>+197</sup><sub>-44</sub>)</td><td>-1847<sup>+73</sup><sub>-16</sub></td><td>2.4<sup>+0.8</sup><sub>-0.7</sub></td><td>7.2<sup>+0.4</sup><sub>-0.3</sub></td><td>2.4<sup>+3.5</sup><sub>-0.7</sub></td><td>0.51<sup>+0.12</sup><sub>-0.13</sub></td><td>-1.36 to -1.8</td></tr><tr><td>Ylgr</td><td>(205<sup>+40</sup><sub>-35</sub>, 2766<sup>+68</sup><sub>-556</sub>, 556<sup>+30</sup><sub>-25</sub>)</td><td>-1219<sup>+20</sup><sub>-19</sub></td><td>11.5<sup>+0.1</sup><sub>-0.1</sub></td><td>20.7<sup>+1.2</sup><sub>-1.1</sub></td><td>11.2<sup>+0.8</sup><sub>-0.7</sub></td><td>0.29<sup>+0.02</sup><sub>-0.02</sub></td><td>-1.87</td></tr><tr><td>Fjorm</td><td>(831<sup>+60</sup><sub>-62</sub>, -2332<sup>+24</sup><sub>-24</sub>, 877<sup>+43</sup><sub>-40</sub>)</td><td>-1123<sup>+15</sup><sub>-15</sub></td><td>9.1<sup>+0.1</sup><sub>-0.1</sub></td><td>29.1<sup>+1.1</sup><sub>-1.1</sub></td><td>19.7<sup>+1.0</sup><sub>-1.0</sub></td><td>0.52<sup>+0.01</sup><sub>-0.01</sub></td><td>-2.2</td></tr><tr><td>Kshir</td><td>(18<sup>+6</sup><sub>-4</sub>, 2755<sup>+60</sup><sub>-53</sub>, 491<sup>+14</sup><sub>-13</sub>)</td><td>-1268<sup>+12</sup><sub>-11</sub></td><td>13.4<sup>+0.2</sup><sub>-0.2</sub></td><td>16.0<sup>+0.6</sup><sub>-0.5</sub></td><td>8.2<sup>+0.3</sup><sub>-0.3</sub></td><td>0.09<sup>+0.01</sup><sub>-0.01</sub></td><td>-1.78</td></tr><tr><td>Cetus</td><td>(815<sup>+513</sup><sub>-317</sub>, -2416<sup>+84</sup><sub>-1064</sub>, 2287<sup>+1282</sup><sub>-954</sub>)</td><td>-1000<sup>+124</sup><sub>-64</sub></td><td>14.7<sup>+7.2</sup><sub>-4.5</sub></td><td>35.9<sup>+9.9</sup><sub>-3.7</sub></td><td>30.2<sup>+10.9</sup><sub>-4.9</sub></td><td>0.45<sup>+0.14</sup><sub>-0.1</sub></td><td>-2.0</td></tr><tr><td>Svol</td><td>(97<sup>+94</sup><sub>-32</sub>, -1501<sup>+384</sup><sub>-248</sub>, 224<sup>+107</sup><sub>-54</sub>)</td><td>-1566<sup>+89</sup><sub>-61</sub></td><td>5.9<sup>+0.6</sup><sub>-0.6</sub></td><td>10.0<sup>+2.8</sup><sub>-1.0</sub></td><td>5.0<sup>+0.9</sup><sub>-0.8</sub></td><td>0.28<sup>+0.07</sup><sub>-0.05</sub></td><td>-1.98 ± 0.10</td></tr><tr><td>Gunntha</td><td>(69<sup>+14</sup><sub>-7</sub>, 852<sup>+67</sup><sub>-77</sub>, 218<sup>+30</sup><sub>-34</sub>)</td><td>-1765<sup>+31</sup><sub>-28</sub></td><td>4.2<sup>+0.3</sup><sub>-0.4</sub></td><td>7.2<sup>+0.3</sup><sub>-0.2</sub></td><td>3.8<sup>+0.4</sup><sub>-0.4</sub></td><td>0.27<sup>+0.03</sup><sub>-0.02</sub></td><td>...</td></tr><tr><td>Slidr</td><td>(1076<sup>+217</sup><sub>-149</sub>, -1358<sup>+23</sup><sub>-23</sub>, 1831<sup>+126</sup><sub>-102</sub>)</td><td>-1086<sup>+41</sup><sub>-32</sub></td><td>8.7<sup>+0.1</sup><sub>-0.1</sub></td><td>32.3<sup>+3.5</sup><sub>-2.5</sub></td><td>29.1<sup>+3.4</sup><sub>-2.4</sub></td><td>0.58<sup>+0.03</sup><sub>-0.02</sub></td><td>-1.8</td></tr><tr><td>M92</td><td>(361<sup>+9</sup><sub>-9</sub>, 181<sup>+39</sup><sub>-41</sub>, 544<sup>+83</sup><sub>-70</sub>)</td><td>-1639<sup>+12</sup><sub>-10</sub></td><td>3.0<sup>+0.2</sup><sub>-0.1</sub></td><td>10.7<sup>+0.2</sup><sub>-0.2</sub></td><td>9.9<sup>+0.5</sup><sub>-0.6</sub></td><td>0.56<sup>+0.01</sup><sub>-0.01</sub></td><td>-2.16 ± 0.05</td></tr><tr><td>NGC 6397</td><td>(75<sup>+5</sup><sub>-5</sub>, -586<sup>+15</sup><sub>-29</sub>, 222<sup>+23</sup><sub>-14</sub>)</td><td>-1851<sup>+11</sup><sub>-6</sub></td><td>3.4<sup>+0.1</sup><sub>-0.1</sub></td><td>6.4<sup>+0.1</sup><sub>-0.1</sub></td><td>3.7<sup>+0.2</sup><sub>-0.1</sub></td><td>0.3<sup>+0.01</sup><sub>-0.01</sub></td><td>...</td></tr><tr><td>NGC 3201</td><td>(975<sup>+48</sup><sub>-48</sub>, 2860<sup>+32</sup><sub>-32</sub>, 296<sup>+5</sup><sub>-5</sub>)</td><td>-1110<sup>+10</sup><sub>-10</sub></td><td>8.5<sup>+0.0</sup><sub>-0.0</sub></td><td>30.5<sup>+0.8</sup><sub>-0.8</sub></td><td>12.3<sup>+0.3</sup><sub>-0.3</sub></td><td>0.56<sup>+0.01</sup><sub>-0.01</sub></td><td>...</td></tr><tr><td>Ophiuchus</td><td>(507<sup>+387</sup><sub>-202</sub>, -160<sup>+34</sup><sub>-41</sub>, 1192<sup>+91</sup><sub>-98</sub>)</td><td>-1490<sup>+130</sup><sub>-84</sub></td><td>3.9<sup>+0.3</sup><sub>-0.4</sub></td><td>14.2<sup>+4.9</sup><sub>-2.7</sub></td><td>14.1<sup>+4.9</sup><sub>-2.6</sub></td><td>0.58<sup>+0.11</sup><sub>-0.1</sub></td><td>-1.80 ± 0.09</td></tr><tr><td>Atlas</td><td>(757<sup>+39</sup><sub>-35</sub>, -1817<sup>+34</sup><sub>-33</sub>, 2093<sup>+100</sup><sub>-86</sub>)</td><td>-1061<sup>+12</sup><sub>-12</sub></td><td>11.7<sup>+0.3</sup><sub>-0.3</sub></td><td>32.4<sup>+1.0</sup><sub>-0.9</sub></td><td>28.6<sup>+1.1</sup><sub>-1.0</sub></td><td>0.47<sup>+0.01</sup><sub>-0.01</sub></td><td>-2.22 ± 0.03</td></tr><tr><td>C-7</td><td>(1059<sup>+397</sup><sub>-215</sub>, 706<sup>+17</sup><sub>-25</sub>, 728<sup>+107</sup><sub>-69</sub>)</td><td>-1319<sup>+100</sup><sub>-65</sub></td><td>3.5<sup>+0.0</sup><sub>-0.0</sub></td><td>21.0<sup>+2.2</sup><sub>-2.8</sub></td><td>18.1<sup>+5.1</sup><sub>-2.8</sub></td><td>0.72<sup>+0.05</sup><sub>-0.04</sub></td><td>...</td></tr><tr><td>C-3</td><td>(142<sup>+538</sup><sub>-64</sub>, 468<sup>+1110</sup><sub>-1185</sub>, 872<sup>+910</sup><sub>-510</sub>)</td><td>-1571<sup>+338</sup><sub>-111</sub></td><td>5.7<sup>+2.0</sup><sub>-1.2</sub></td><td>10.0<sup>+13.2</sup><sub>-1.4</sub></td><td>8.7<sup>+12.5</sup><sub>-1.3</sub></td><td>0.35<sup>+0.18</sup><sub>-0.1</sub></td><td>...</td></tr><tr><td>Palca</td><td>(91<sup>+37</sup><sub>-24</sub>, -1830<sup>+28</sup><sub>-29</sub>, 1076<sup>+138</sup><sub>-128</sub>)</td><td>-1300<sup>+30</sup><sub>-28</sub></td><td>10.8<sup>+0.4</sup><sub>-0.3</sub></td><td>16.5<sup>+1.5</sup><sub>-1.3</sub></td><td>12.7<sup>+1.5</sup><sub>-1.4</sub></td><td>0.21<sup>+0.03</sup><sub>-0.03</sub></td><td>-2.02 ± 0.23</td></tr><tr><td>Sylgr</td><td>(602<sup>+141</sup><sub>-202</sub>, -702<sup>+28</sup><sub>-24</sub>, 2220<sup>+94</sup><sub>-153</sub>)</td><td>-1192<sup>+36</sup><sub>-61</sub></td><td>8.7<sup>+0.0</sup><sub>-0.0</sub></td><td>24.6<sup>+2.4</sup><sub>-3.7</sub></td><td>23.8<sup>+2.4</sup><sub>-3.7</sub></td><td>0.48<sup>+0.03</sup><sub>-0.06</sub></td><td>-2.92 ± 0.06</td></tr><tr><td>Gaia-6</td><td>(125<sup>+51</sup><sub>-71</sub>, 907<sup>+342</sup><sub>-229</sub>, 557<sup>+288</sup><sub>-204</sub>)</td><td>-1593<sup>+129</sup><sub>-70</sub></td><td>6.0<sup>+1.4</sup><sub>-1.8</sub></td><td>9.5<sup>+3.1</sup><sub>-0.3</sub></td><td>6.9<sup>+3.6</sup><sub>-0.9</sub></td><td>0.3<sup>+0.08</sup><sub>-0.1</sub></td><td>-1.16</td></tr><tr><td>Gaia-9</td><td>(393<sup>+37</sup><sub>-38</sub>, 1928<sup>+47</sup><sub>-61</sub>, 852<sup>+22</sup><sub>-20</sub>)</td><td>-1255<sup>+15</sup><sub>-18</sub></td><td>8.7<sup>+0.1</sup><sub>-0.1</sub></td><td>20.8<sup>+0.8</sup><sub>-0.9</sub></td><td>14.7<sup>+0.6</sup><sub>-0.6</sub></td><td>0.41<sup>+0.01</sup><sub>-0.01</sub></td><td>-1.94</td></tr><tr><td>Gaia-10</td><td>(2189<sup>+66</sup><sub>-64</sub>, 287<sup>+43</sup><sub>-43</sub>, 1542<sup>+155</sup><sub>-140</sub>)</td><td>-1051<sup>+20</sup><sub>-18</sub></td><td>4.3<sup>+0.4</sup><sub>-0.3</sub></td><td>37.7<sup>+1.6</sup><sub>-1.4</sub></td><td>37.2<sup>+1.6</sup><sub>-1.5</sub></td><td>0.8<sup>+0.01</sup><sub>-0.01</sub></td><td>-1.4</td></tr><tr><td>Gaia-12</td><td>(9834<sup>+5378</sup><sub>-4378</sub>, 7340<sup>+608</sup><sub>-801</sub>, 794<sup>+80</sup><sub>-107</sub>)</td><td>-433<sup>+98</sup><sub>-142</sub></td><td>18.5<sup>+0.9</sup><sub>-1.2</sub></td><td>194.3<sup>+96.8</sup><sub>-75.0</sub></td><td>83.0<sup>+42.9</sup><sub>-32.5</sub></td><td>0.83<sup>+0.05</sup><sub>-0.08</sub></td><td>-2.6</td></tr><tr><td>Indus</td><td>(99<sup>+75</sup><sub>-79</sub>, -1121<sup>+35</sup><sub>-36</sub>, 2211<sup>+61</sup><sub>-49</sub>)</td><td>-1232<sup>+17</sup><sub>-15</sub></td><td>12.6<sup>+0.2</sup><sub>-0.1</sub></td><td>18.9<sup>+1.0</sup><sub>-0.9</sub></td><td>17.8<sup>+0.9</sup><sub>-0.8</sub></td><td>0.2<sup>+0.02</sup><sub>-0.02</sub></td><td>-1.96 ± 0.41</td></tr><tr><td>Jhelum</td><td>(594<sup>+49</sup><sub>-56</sub>, -356<sup>+19</sup><sub>-17</sub>, 2557<sup>+62</sup><sub>-72</sub>)</td><td>-1193<sup>+17</sup><sub>-22</sub></td><td>8.7<sup>+0.2</sup><sub>-0.2</sub></td><td>24.5<sup>+1.1</sup><sub>-1.3</sub></td><td>24.3<sup>+1.1</sup><sub>-1.3</sub></td><td>0.48<sup>+0.01</sup><sub>-0.01</sub></td><td>-1.83 ± 0.34</td></tr><tr><td>Phoenix</td><td>(107<sup>+10</sup><sub>-10</sub>, -1563<sup>+32</sup><sub>-32</sub>, 1578<sup>+56</sup><sub>-62</sub>)</td><td>-1259<sup>+10</sup><sub>-12</sub></td><td>11.7<sup>+0.4</sup><sub>-0.5</sub></td><td>18.1<sup>+0.3</sup><sub>-0.3</sub></td><td>15.6<sup>+0.3</sup><sub>-0.3</sub></td><td>0.22<sup>+0.01</sup><sub>-0.01</sub></td><td>-2.70 ± 0.06</td></tr><tr><td>NGC5466</td><td>(1769<sup>+144</sup><sub>-114</sub>, 619<sup>+38</sup><sub>-35</sub>, 1373<sup>+93</sup><sub>-80</sub>)</td><td>-1098<sup>+29</sup><sub>-24</sub></td><td>4.8<sup>+0.3</sup><sub>-0.2</sub></td><td>33.7<sup>+2.3</sup><sub>-1.8</sub></td><td>31.8<sup>+2.7</sup><sub>-1.7</sub></td><td>0.75<sup>+0.01</sup><sub>-0.01</sub></td><td>...</td></tr><tr><td>M5</td><td>(1366<sup>+69</sup><sub>-57</sub>, -353<sup>+20</sup><sub>-18</sub>, 931<sup>+46</sup><sub>-40</sub>)</td><td>-1246<sup>+17</sup><sub>-14</sub></td><td>3.4<sup>+0.1</sup><sub>-0.1</sub></td><td>24.8<sup>+1.0</sup><sub>-0.8</sub></td><td>23.6<sup>+1.1</sup><sub>-0.9</sub></td><td>0.76<sup>+0.01</sup><sub>-0.01</sub></td><td>-1.34 ± 0.05</td></tr><tr><td>C-20</td><td>(1329<sup>+526</sup><sub>-350</sub>, -3042<sup>+131</sup><sub>-167</sub>, 3823<sup>+503</sup><sub>-484</sub>)</td><td>-800<sup>+65</sup><sub>-59</sub></td><td>20.8<sup>+1.3</sup><sub>-1.3</sub></td><td>58.5<sup>+12.0</sup><sub>-8.8</sub></td><td>52.4<sup>+11.3</sup><sub>-8.5</sub></td><td>0.47<sup>+0.05</sup><sub>-0.04</sub></td><td>-2.44</td></tr><tr><td>NGC7089</td><td>(800<sup>+713</sup><sub>-270</sub>, -638<sup>+414</sup><sub>-455</sub>, 359<sup>+99</sup><sub>-123</sub>)</td><td>-1504<sup>+307</sup><sub>-104</sub></td><td>2.9<sup>+1.7</sup><sub>-0.7</sub></td><td>14.7<sup>+12.7</sup><sub>-3.0</sub></td><td>10.9<sup>+7.5</sup><sub>-3.8</sub></td><td>0.71<sup>+0.06</sup><sub>-0.06</sub></td><td>...</td></tr><tr><td>C-19</td><td>(383<sup>+53</sup><sub>-47</sub>, -210<sup>+48</sup><sub>-46</sub>, 2712<sup>+258</sup><sub>-253</sub>)</td><td>-1232<sup>+21</sup><sub>-21</sub></td><td>9.3<sup>+1.0</sup><sub>-1.0</sub></td><td>21.6<sup>+0.5</sup><sub>-0.5</sub></td><td>21.6<sup>+0.6</sup><sub>-0.6</sub></td><td>0.4<sup>+0.04</sup><sub>-0.03</sub></td><td>-3.38 ± 0.06</td></tr><tr><td>Elqui</td><td>(2072<sup>+543</sup><sub>-602</sub>, -273<sup>+166</sup><sub>-191</sub>, 4324<sup>+321</sup><sub>-352</sub>)</td><td>-868<sup>+27</sup><sub>-35</sub></td><td>12.1<sup>+1.8</sup><sub>-2.0</sub></td><td>54.0<sup>+4.6</sup><sub>-6.2</sub></td><td>53.9<sup>+4.6</sup><sub>-6.3</sub></td><td>0.64<sup>+0.07</sup><sub>-0.08</sub></td><td>-2.22 ± 0.37</td></tr><tr><td>Chenab</td><td>(2469<sup>+463</sup><sub>-286</sub>, -4062<sup>+601</sup><sub>-509</sub>, 3735<sup>+338</sup><sub>-338</sub>)</td><td>-690<sup>+52</sup><sub>-53</sub></td><td>22.0<sup>+2.7</sup><sub>-2.7</sub></td><td>81.0<sup>+12.8</sup><sub>-10.4</sub></td><td>69.1<sup>+10.6</sup><sub>-8.2</sub></td><td>0.58<sup>+0.01</sup><sub>-0.01</sub></td><td>-1.78 ± 0.34</td></tr><tr><td>AliqaUma</td><td>(738<sup>+138</sup><sub>-71</sub>, -1838<sup>+96</sup><sub>-64</sub>, 2025<sup>+223</sup><sub>-126</sub>)</td><td>-1067<sup>+30</sup><sub>-18</sub></td><td>11.6<sup>+0.4</sup><sub>-0.3</sub></td><td>31.9<sup>+2.6</sup><sub>-1.5</sub></td><td>27.9<sup>+3.0</sup><sub>-1.6</sub></td><td>0.47<sup>+0.03</sup><sub>-0.02</sub></td><td>-2.30 ± 0.06</td></tr><tr><td>Phlegethon</td><td>(815<sup>+120</sup><sub>-93</sub>, 1882<sup>+39</sup><sub>-37</sub>, 231<sup>+11</sup><sub>-10</sub>)</td><td>-1272<sup>+32</sup><sub>-28</sub></td><td>5.5<sup>+0.0</sup><sub>-0.0</sub></td><td>22.1<sup>+1.8</sup><sub>-1.4</sub></td><td>9.4<sup>+0.9</sup><sub>-0.7</sub></td><td>0.6<sup>+0.02</sup><sub>-0.02</sub></td><td>-1.96 ± 0.05</td></tr><tr><td>GD-1</td><td>(164<sup>+35</sup><sub>-28</sub>, 2952<sup>+66</sup><sub>-61</sub>, 938<sup>+23</sup><sub>-22</sub>)</td><td>-1153<sup>+15</sup><sub>-14</sub></td><td>14.1<sup>+0.1</sup><sub>-0.1</sub></td><td>23.0<sup>+1.1</sup><sub>-1.0</sub></td><td>14.8<sup>+0.7</sup><sub>-0.7</sub></td><td>0.24<sup>+0.02</sup><sub>-0.02</sub></td><td>-2.24 ± 0.21</td></tr></tbody></table>

Note. From left to right, the columns provide the stream's name, action components (J, energy (E), pericentric distance ($r_{\text{peri}}$), apocentric distance ($r_{\text{apo}}$), maximum height of the orbit from the Galactic plane ($z_{\text{max}}$), eccentricity (ecc), and [Fe/H] measurements (most of these are spectroscopic and a few are photometric). The (J, E) and other orbital parameters are derived in this study. The quoted orbital parameter values are the medians of the sampled posterior distributions and the corresponding uncertainties reflect their 16th and 84th percentiles. The [Fe/H] values of Gaia-6, Gaia-9, Gaia-10, and Gaia-12 correspond to the median of the spectroscopic sample that we obtained in this study. The other streams with spectroscopic [Fe/H] include LMS-1 (its value is taken from Malhan et al. 2021), Gjöll, Ylgr, Slidr, Fjörm (Ibata et al. 2019b), Jhelum, Chenab, Elqui, Ophiuchus, Orphan, Palca, Indus (Li et al. 2021a), Fimbulthul (Ibata et al. 2019a), Gaia-1, C-2 and Hrid (Malhan et al. 2020), Cetus (Yam et al. 2013), Sylgr (Roederer & Gnedin 2019), GD-1 (Malhan & Ibata 2019), Kshir (Malhan et al. 2019a), C-20 (Yuan et al. 2021), Pal 5 (Ishigaki et al. 2016), Atlas, and AliqaUma (Li et al. 2021b). Streams with photometric [Fe/H] include Phlegethon, Svöl, M92, and M5 (their values are taken from Martin et al. 2022b).

As a final passing analysis, and not to deviate too much from the prime objective of the paper, we quickly compare the distribution of the orbital phase and eccentricity of all the halo objects (as shown in Figure 13 of Appendix A). First, we observe a pileup of objects at the pericenter and at the apocenter, and this is more prevalent for globular clusters and

stellar streams and not so much for satellite galaxies. Particularly, in the case of streams, we note that more objects are piled-up at the pericenter than at the apocenter. This effect points toward our inefficiency in detecting those streams that, at the present day, could be close to their apocenters (at distances $D \gtrsim 30$ kpc). This inefficiency, in part, is also
because of Gaia's limiting magnitude at $G \sim 21$. Our result is different from that of Li et al. (2021a), who find that more streams (in their sample of 12 streams) are piled-up at the apocenter. Second, we find that most of the objects (be it clusters, streams, or satellites) have eccentricities $e \approx 0.5$, and it is rare for the objects to possess very radial orbits ($e \approx 1$) or very circular orbits ($e \approx 0$). This last inference, with regard to streams, is consistent with that of Li et al. (2021a).

In summary, we now possess (J, E) information for a total of $n = 257$ halo objects of the Milky Way (as shown in Figure 2). In the next section, we process the entire (J, E) data to detect groups of objects (i.e., mergers). Therefore, at this stage, it is important to clarify that some of the objects are being counted twice in our data set. These objects include those systems that have counterparts both in the globular cluster catalog and the stream catalog. For instance, a subset of these objects include Pal 5, NGC 3201, ω Centauri and M5. One possible way to proceed would be to remove their counterparts from either of the catalogs. However, there could be many other streams in our catalog that could be physically associated to other globular clusters (e.g., see Section 6) or even to other streams (e.g., Orphan-Chenab, Koposov et al. 2019; Palca-Cetus, Chang et al. 2020; AliqaUma-Atlas, Li et al. 2021b), and it is a difficult task to separate these plausible associations. We therefore consider it to be less biased to proceed with all of the detected structures. Prior associations will be discussed in our final grouping analysis.

### 3. Detecting Groups of Objects in (J, E) Space

To search for the Milky Way mergers, we essentially process the data shown in Figure 2 and detect groups of objects that tightly clump together in the (J, E) space. To detect these groups, we employ the ENLINK software (Sharma & Johnston 2009) and couple it with a statistical procedure that accounts for the uncertainties in the (J, E) values of every object. Below, we first briefly describe the working of ENLINK and then our procedure to detect groups.

#### 3.1. Description of ENLINK

ENLINK is a density-based hierarchical group-finding algorithm that detects groups of any shape and density in a multidimensional data set. This software employs nonparametric methods to find groups, i.e., it makes no assumptions about the number of groups being identified or their form. These functionalities of ENLINK are particularly useful for our study because a priori we neither know the number of groups (i.e., number of mergers) that are present in the (J, E) data set, nor the shapes of these groups (because objects that accrete inside the same merging galaxy can realize extended/irregular ellipsoidal shapes in (J, E) space; e.g., Wu et al. 2022).

To detect groups in the data set, ENLINK does not use the typical Euclidean metric, but builds a locally adaptive Mahalanobis (LAM) metric. The importance of this metric can be explained as follows. Generally speaking, the task of finding groups in a given data set ultimately boils down to computing “distances” between different data points. Then, those data points that lie at smaller distances from each other form part of the same group. In a scenario where the correlation between different dimensions of the data set are zero or negligible, one can simply adopt the Euclidean metric to compute these distances. In this case, the distance between two data points $x_i$ and $x_j$ is given by

$s^2(x_i, x_j) = (x_i - x_j)^T \cdot (x_i - x_j)$, where $x$ is a 1D matrix whose length equals the dimension of the data set. However, in real data sets, correlations between different dimensions are nonzero. Particularly in our case, one expects significant correlation in the space constructed with J and E dimensions. Therefore, to find groups in such a correlated data set, one effectively requires a multivariate equivalent of the Euclidean distance. This is the importance of LAM that ENLINK employs, because the Mahalanobis distance is the distance between a point and a distribution (and not between two data points). At its heart, ENLINK uses the LAM metric, where the distance between the two data points (under discrete approximation) is defined as

$$s^2(x_i, x_j) = |\Sigma(x_i, x_j)|^{1/d} \cdot (x_i - x_j)^T \cdot \Sigma^{-1}(x_i, x_j) \cdot (x_i - x_j), \quad (3)$$

where $d$ is the dimension of the data, $\Sigma$ is the covariance matrix, $\Sigma(x_i, x_j) = 0.5[\Sigma(x_i) + \Sigma(x_j)]$ and $\Sigma^{-1}(x_i, x_j) = 0.5[\Sigma^{-1}(x_i) + \Sigma^{-1}(x_j)]$.

The above formula can be intuitively understood as follows. Consider the term $(x_i - x_j)^T \cdot \Sigma^{-1}$. Here, $(x_i - x_j)$ is the distance between two data points. This is then multiplied by the inverse of the covariance matrix $\Sigma$ (or divided by the covariance matrix). So, this is essentially a multivariate equivalent of the regular standardization $y = (x - \mu)/\sigma$. The effect of dividing by covariance is that if the $x$ values in the data set are strongly correlated, then the covariance will be high and dividing by this large covariance will reduce the distance. On the other hand, if the $x$ are not correlated, then the covariance is small and the distance is not reduced by much. The overall workings and implementation of ENLINK are detailed in Sharma & Johnston (2009), and this software has also been previously applied to various data sets (e.g., Sharma et al. 2010; Wu et al. 2022).

#### 3.2. Applying ENLINK

To detect groups, we work in the four-dimensional space of $x_i \equiv (J_{R,i}, J_{\phi,i}, J_{z,i}, E_i)$, where $i$ represents a given halo object and the units of J and E are kpc km s⁻¹ and km² s⁻², respectively. The reason for working with both J and E quantities is that their combined information allowed us to detect several groups (as we show below). Initially, we operated ENLINK only in the three-dimensional space of J. However, this resulted in the detection of the Sagittarius group (Ibata et al. 1994; Bellazzini et al. 2020) (although with unusual membership of objects), the Arjuna/Sequoia/I'toi group (Naidu et al. 2020; Bonaca et al. 2021), and one to two other very low-significance groups. At first, this may seem odd that ENLINK requires the additional (redundant) E information to find high-significance groups because J fully characterized the orbits and the parameter E brings no additional dynamical information. However, this oddity relates to the uncertainties on J and E. For instance, the relative uncertainties on $(J_{R,i}, J_{\phi,i}, J_{z,i})$ for all the objects in our sample (on average) are (12%, 17%, 9%), while the relative uncertainty on E is only 2%. Therefore, ENLINK prefers these precise values of E in addition to J as this helps it to easily distinguish between different groups.

The ENLINK parameters that we use are neighbors, min_cluster_size, min_peak_height, cluster_separation_method, density_method, and gme-tric. neighbors is the “smoothing” that is used to compute a local density for each data point because ENLINK
first estimates the density and then finds groups in the density field. To search for groups in a $d$-dimensional data set, ENLINK requires $\text{neighbors} \ge (d+1)$. In our case, $d=4$ (3 components of $J$ and $E$) and therefore we set $\text{neighbors}=5$. Second, we set $\text{min_cluster_size}=5$. This is because it is difficult to find groups smaller than the smoothing length (i.e., we satisfy the $\text{min_cluster_size} \ge \text{neighbors}$ condition of ENLINK). $\text{min_peak_height}$ can be thought of as the signal-to-noise ratio of the detected groups, and we set $\text{min_peak_height}=3.0$. For the parameters $\text{cluster\_separation_method}$ and $\text{gmetric}$, we adopt the default values (i.e., 0). Further, we set $\text{density_method}=\text{sbr}$ as this uses an adaptive metric to detect groups. We also tried different metric definitions, but these gave very similar results to those we obtained from the above parameter setting.<fnref n="15" /> Our experimentation with various parameter settings makes us confident that we are detecting robust groups.

Before unleashing ENLINK onto the $(J, E)$ data set, we couple it with a statistical procedure that accounts for the dispersion in the $(J, E)$ values of the objects (these dispersions are visible in Figure 2). This is important because ENLINK itself does not account for the dispersion associated with each data point. This statistical procedure can be explained as follows. Fundamentally, we want to compute a “group probability” ($P_{\text{Group}}$) for each halo object, such that this probability is higher for those objects that belong to the groups detected by ENLINK. To compute this $P_{\text{Group}}$ value, we undertake an iterative procedure.

In the first iteration, each halo object is represented by a single $(J, E)$ value that is sampled from its MCMC chain (we obtained these MCMC chains in Section 2.2). At this stage, the total number of $(J, E)$ data points equals the total number of objects (i.e., 257). After this, we process this $(J, E)$ data using ENLINK. An attribute that ENLINK returns is a 1D array labels. labels has the same length as the number of input data points, and it stores the grouping information. That is, all the elements in labels possess integer values in the range 1 to $n$, where $n$ is the total number of groups detected by ENLINK, and elements that form part of the same group receive the same values. Furthermore, elements for which $\text{labels}=1$ correspond to those objects that form part of the largest group. For all the objects with $\text{labels} \ge 2$, we explicitly set their probability of group membership at iteration $i$ to be $P_{\text{Group},i}=1$. Among objects with $\text{labels}=1$, we accept only those objects that possess $\text{density} > 99\%$ percentile and set their $P_{\text{Group},i}=1$, while the remaining low $\text{density}$ objects are set as $P_{\text{Group},i}=0$.<fnref n="16" /> In the next iteration, a new set of $(J, E)$ values is sampled and the above procedure is repeated. Note that in this new iteration, the input $(J, E)$ data has changed, and therefore, the same object can now belong to a different group and thus receive a different labels value

and a different $P_{\text{Group},i}$ value. We iterate this procedure 1000 times. This produces, for each halo object, a one-dimensional array (of length 1000) that contains a combination of 0 s or 1 s. For each halo object, we take the average of this array and this we interpret as the group probability $P_{\text{Group}}$ of that object. The $P_{\text{Group}}$ parameter can be defined as the probability of an object belonging to a group in $(J, E)$ space. Indeed, those halo objects that lie in denser regions of $(J, E)$ space—i.e., objects whose $(J, E)$ distributions overlap significantly—will possess higher $P_{\text{Group}}$ values.

Figure 3 shows the $(J, E)$ distribution of the halo objects as a function of the computed $P_{\text{Group}}$ values. In this figure, each object is represented by the median of its $(J, E)$ distribution. It can be seen that different objects possess different $P_{\text{Group}}$ values. We also note that objects with high $P_{\text{Group}}$ values lie in denser regions of $(J, E)$ space, suggesting that our procedure of detecting groups has worked as desired. In Figure 3, one can already visually identify many possible groupings—comprising those objects that possess high $P_{\text{Group}}$ values and that appear well separated from other groups.

### 3.3. Detecting High-significance Groups

Due to the relatively large $(J, E)$ uncertainties, the ENLINK algorithm’s output of proposed groupings varies considerably over the 1000 random iterations described above. This means that the proposed groups cannot be immediately used to identify the Milky Way’s mergers.

Therefore, we proceed by first defining a threshold value $P_{\text{Threshold}}$, such that objects with $P_{\text{Group}} \ge P_{\text{Threshold}}$ belong to high-significance groups, and this corresponds to a likely detection. To find a suitable $P_{\text{Threshold}}$ value, we follow a pragmatic approach. We repeat the above analysis of computing the $P_{\text{Group}}$ values of all the halo objects, except this time we use a “randomized” version of our real $(J, E)$ data. This randomized data is artificially created, where each object is first assigned a random orbital pole and then its new $(J, E)$ values are computed. These randomized $(J, E)$ data are shown in Figure 14 in Appendix B. Such a randomization procedure erases any plausible correlations between the objects in $(J, E)$ space. For the resulting PDF of the new $P_{\text{Group}}$ values (that is shown in Figure 16), its 90 percentile limit motivates setting a threshold at $P_{\text{Threshold}}=0.3$ for a $2\sigma$ detection. This procedure may seem convoluted, but it is required by the astrometric uncertainties which project in a complicated, nonlinear way into $(J, E)$ space (hence the usual techniques of error propagation would not have been appropriate). This method of finding the $P_{\text{Threshold}}$ value is detailed in Appendix B. Consequently, for the real $P_{\text{Group}}$ values (shown in Figure 4), all those objects that possess $P_{\text{Group}} \ge P_{\text{Threshold}}$ are considered as high-significance groups.

The selection $P_{\text{Group}} \ge P_{\text{Threshold}}$ yields 108 objects (42% of the total $n=257$ objects), and these are shown in Figure 5. These objects include 81 globular clusters, 25 stellar streams, and 2 satellite galaxies. This figure also shows different objects being linked by straight lines. This “link,” between two given objects, represents the frequency with which these objects were classified as members of the same group (as per the procedure described in Section 3.2). Thicker links imply higher frequency. In Figure 5, these links are pruned by removing those cases where two objects resulted in the same group in less than (approximately) one third (300/1000) of the realizations. Due to this pruning, a couple of objects can be seen without

<footnote n="15">¹⁵ For example, instead of using the adaptive metric, we defined a constant metric using the uncertainties on $(J, E)$ by setting $\text{gmetric}=2$ and using the custom\_metric parameter. We made this test because as we are dealing with a very low number of data points (only 257 points), we wanted to ensure that the detected groups are robust and are not noise driven. However, in this case we found similar results to those with the original ENLINK setting.</footnote>

<footnote n="16">¹⁶ The reason that we make such a distinction for the objects in the labels = 1 group is that a majority of objects in this largest group are those that could not be associated with any “well-defined group” by ENLINK (these represent the background objects). However, even in this group, some of the high-density objects may still represent a real merger. Therefore, in order to consider these potential objects of interest, we accept only those objects that satisfy the threshold density criteria.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 3.** (J, E) distribution of the halo objects as a function of their group probability $P_{\text{Group}}$ (see Section 3). In this plot, each object is represented by the median of its (J, E) distribution (which is shown in Figure 2). The globular clusters are denoted by “stars,” streams by “circles,” and satellite galaxies by “squares.” Note that objects with higher $P_{\text{Group}}$ values lie in the denser regions of this (J, E) space. Such objects with high $P_{\text{Group}}$ values, which also clump together in (J, E) space, form part of the same group. In panel (b), we label all the high-significance groups.</figcaption></figure>

any links, even though they satisfy the condition $P_{\text{Group}} \ge P_{\text{Threshold}}$, and it is therefore difficult to associate them with one unique group. The power of Figure 5 is that such a representation automatically reveals the detection of several independent groups.

Figure 5 shows that we have detected nine high-significance groups and the properties of these groups are discussed below.

## 4. Analyzing the Detected Groups

We detect a total of nine distinct groups at $\ge 2\sigma$ significance. Among these, we interpret $N=6$ groups as the mergers of the Milky Way because the remaining three actually contain the

in situ population of the Milky Way (see below). The merger groups comprise 62 halo objects (25% of the total 257 objects considered in our study), including 35 globular clusters, 25 streams, and 2 satellite galaxies. For each of the merger groups, we analyze the objects' memberships (which are also summarized in Table 3), their (J, E) properties (which result from Figure 5), orbital parameters as a function of [Fe/H] (see Figure 6), [Fe/H] distribution function (MDF; see Figure 7), other orbital parameters (see Figure 8), and also estimate the masses of the corresponding progenitor galaxies.

For each group discussed below, we also make a comparison between our object membership and those proposed in previous studies. Therefore, to also facilitate this comparison visually,
<figure><img src="image_11.png" /><figcaption>**Figure 4.** PDF of the group probability $P_{\text{Group}}$ of all the halo objects (see Section 3.3). The vertical line represents the $P_{\text{Threshold}}$ value and all the objects with $P_{\text{Group}} \ge P_{\text{Threshold}}$ belong to high-significance groups. The value quoted in the top-right corner is the median of the distribution and the corresponding uncertainties reflect the 16th and 84th percentiles.</figcaption></figure>

we provide Figure 17 in Appendix C. This figure is constructed by adopting the object-merger associations from other studies (specifically from Massari et al. 2019; Myeong et al. 2019; Kruijssen et al. 2020; Bonaca et al. 2021). We explicitly note that our results are not based on Figure 17, and we use it purely for comparison with our Figure 5.

## 4.1. Group 1: The Sagittarius Merger

The first group we detect is a high-energy and prograde overdensity in (J, E) space. Its member objects possess dynamical properties in the range $E \sim [-0.91, -0.79] \times 10^5 \text{ km}^2 \text{s}^{-2}$, $J_R \sim [1205, 2090] \text{ km s}^{-1} \text{kpc}$, $J_\phi \sim [-2115, -265] \text{ km s}^{-1} \text{kpc}$, $J_z \sim [3285, 5350] \text{ km s}^{-1} \text{kpc}$, $L_\perp \sim [4565, 6835] \text{ km s}^{-1} \text{kpc}$, eccentricity $\sim [0.5, 0.6]$, $r_{\text{peri}} \sim [11, 22] \text{ kpc}$, $r_{\text{apo}} \sim [45, 60] \text{ kpc}$, and $\phi \sim [66^\circ, 86^\circ]$, where $\phi$ defines how “polar” the merger group is. This highly polar group represents the previously known Sagittarius merger (Ibata et al. 1994; Majewski et al. 2003; Bellazzini et al. 2020).

We find that eight objects belong to this group: six globular clusters (namely, Pal 12, Whiting 1, Terzan 7, Terzan 8, Arp 2, NGC 6715/M54), one stream (namely, Elqui), and one satellite (namely, the Sagittarius dSph itself). Our globular cluster member list is similar to those previously reported by other studies (e.g., Massari et al. 2019; Bellazzini et al. 2020; Forbes 2020). We note that our Sagittarius group lacks NGC 2419 as its member, but previous studies have advocated for this association based on the fact that this cluster also lies within the phase-space distribution of the Sagittarius stream (e.g., Sohn et al. 2018; Bellazzini et al. 2020). A possible reason why our analysis does not identify a strong association between NGC 2419 and Sagittarius could be due to this cluster’s large (J, E) uncertainties that arise due to its large observational uncertainties (because it is a very distant cluster, $D_\odot \approx 83$ kpc). On the other hand, our stream-Sagittarius association is completely different from that of Bonaca et al. (2021). Bonaca et al. (2021) found five stream-Sagittarius associations by comparing the $(J_\phi, E)$ values of their stream

sample to the $(J_\phi, E)$ distribution of the mergers previously found by Naidu et al. (2021), but their stream member list does not include Elqui.<fnref n="17" /> In fact, we find that most of their Sagittarius stream members actually belong to the Cetus group (see below). Moreover, given that the Elqui stream is produced from a low-mass dwarf galaxy (Li et al. 2021a), this further suggests that Elqui was likely the satellite dwarf galaxy of the progenitor Sagittarius galaxy (i.e., of the Sagittarius dSph galaxy itself).

We use the above-listed member objects of Sagittarius and analyze their [Fe/H]. The [Fe/H] measurements of streams are taken from Table 2 and for globular clusters we rely on the Harris (2010) catalog. Figure 6 shows the orbital properties of the objects as a function of their [Fe/H]. One can notice that the member objects of Sagittarius possess varied metallicities, and this is consistent with previous studies (e.g., Massari et al. 2019; Bellazzini et al. 2020). To quantify this [Fe/H] distribution, we also construct the MDF shown in Figure 7. This MDF has a median [Fe/H] = -0.85 dex and spans a wide range from -2.22 dex to -0.32 dex.

For the progenitor Sagittarius galaxy, we determine its halo mass ($M_{\text{halo}}$) and stellar mass ($M_*$) as follows. We first determine $M_{\text{halo}}$ using the globular-cluster-to-halo-mass relation (Hudson et al. 2014) and then convert this $M_{\text{halo}}$ to $M_*$ using the stellar-to-halo-mass relation (Read et al. 2017). To this end, we use the masses of the individual globular clusters from Baumgardt et al. (2019). The combined masses of the clusters provide $M_{\text{halo}} \sim 4.5 \times 10^{10} M_\odot$, and this further implies $M_* \sim 13 \times 10^7 M_\odot$. These mass values are similar to those found by previous studies (e.g., Gibbons et al. 2017; Niederste-Ostholt et al. 2012).

Note that such a method provides a very rough estimate of the mass values and is not very accurate. This is because (1) both the Hudson et al. (2014) and Read et al. (2017) relations have some scatter that we do not account for. (2) Such a method makes a strong assumption that the present-day observations of the globular-cluster-to-halo-mass relation and stellar-to-halo-mass relation do not evolve with redshift. Because our estimates are not corrected for redshift, they provide an overestimate of the actual progenitor mass (at merging time). (3) On the other hand, such a mass estimation technique uses knowledge of only member globular clusters and not member stellar streams (some of which could be produced from globular cluster themselves). Therefore, this may underestimate the actual progenitor mass. (4) Such a method does not account for other objects that in principle could belong to the merger groups but were actually not identified by our study. For instance, the globular cluster AM 4 has also been previously linked to the Sagittarius group (Forbes 2020), but we do not identify it here. In view of these limitations, we note that this method provides an approximate value on the mass of the progenitor galaxy (at the time of merging).

## 4.2. Group 2: The Cetus Merger

This group is the most prograde among all the detected groups and possesses dynamical properties in the range $E \sim [-1.09, -0.93] \times 10^5 \text{ km}^2 \text{s}^{-2}$, $J_R \sim [605, 1075] \text{ km s}^{-1} \text{kpc}$, $J_\phi \sim [-2700, -1360] \text{ km s}^{-1} \text{kpc}$, $J_z \sim [1835, 2820] \text{ km s}^{-1} \text{kpc}$,

<footnote n="17">¹⁷ Our stream sample contains all the streams that Bonaca et al. (2021) associated with Sagittarius, except for “Turranburra.”</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 5.** (*J*, *E*) distribution of the groups detected in our study. The plot shows several independent groups that comprise those objects with high probabilities (i.e., *P*<sub>Group</sub> ≥ *P*<sub>Threshold</sub>; see Section 3.3). The left panel shows *J*<sub>φ</sub> vs. *J*<sub>R</sub>, and the objects are colored by their *J*<sub>z</sub> values. The right panel shows *J*<sub>φ</sub> vs. *E*, colored by *L*<sub>⊥</sub>. The gray points are all the remaining objects with *P*<sub>Group</sub> < *P*<sub>Threshold</sub>. The straight lines between any two objects indicate the frequency of these objects being members of the same group—the thicker the line, the higher is this frequency. These lines are colored using the same scheme described above. Such a representation automatically reveals several independent groups.</figcaption></figure>

*L*<sub>⊥</sub> ~ [2905, 4635] km s<sup>-1</sup> kpc, eccentricity ~ [0.4, 0.6], *r*<sub>peri</sub> ~ [8, 16] kpc, *r*<sub>apo</sub> ~ [31, 43] kpc, and *φ* ~ [50°, 65°]. It corresponds to the previously known Cetus merger (Newberg et al. 2009; Yuan et al. 2019). Inspecting Figure 5, it can be seen that Cetus is situated in the vicinity of the Sagittarius group. However, these two groups overall possess quite different *J*<sub>R</sub> and *J*<sub>φ</sub> components and different orbital properties, and they can also be distinguished on the basis of their [Fe/H] properties (the Cetus members are overall more metal poor than the Sagittarius members).

We find that six objects belong to this group: four streams (namely, Cetus itself, Slidr, Atlas, AliqaUma), one cluster (namely, NGC 5824), and one satellite galaxy (Willman 1). Among the stream member list, AliqaUma and Atlas were recently associated with the Cetus stream by Li et al. (2021a). On the other hand, Bonaca et al. (2021) associated most of these streams with the Sagittarius group. Bonaca et al. (2021) found three other streams to be associated with Cetus, but these streams are not present in our data sample.<fnref n="18" /> Furthermore, we could not find the streams C-20 and Palca as members of this group, but their associations have been suggested by previous studies (e.g., Chang et al. 2020; Li et al. 2021a; Yuan et al. 2021). As for the globular-cluster–Cetus association, NGC 5824 has been previously linked with the Cetus stream by various studies on the basis that this cluster lies within the phase-space distribution of the Cetus stream (e.g., Newberg et al. 2009; Yuan et al. 2019; Chang et al. 2020). However, other studies indicate that NGC 5824 is associated with the Sagittarius group (Massari et al. 2019; Forbes 2020).

Surprisingly, some of the previous studies do not mention the Cetus group in their analysis. For instance, Massari et al. (2019) made a selection in *L*<sub>z</sub>–*L*<sub>⊥</sub> space to identify the Sagittarius globular clusters and found that this integral-of-motion space also

contains NGC 5824, so they assigned it to the Sagittarius group. On the other hand, Forbes (2020) identify their merger groups by combining the orbit information of globular clusters from Massari et al. (2019) and the ages and [Fe/H] from Kruijssen et al. (2019). Also, they guide their analysis with the previously known cluster–merger memberships from Massari et al. (2019). A possible reason that these studies could not identify Cetus is because they were analyzing only globular clusters, and the Cetus group (likely) contains only one such object—NGC 5824 itself. However, we are able to detect Cetus because we have combined the globular cluster information with that of streams and satellites, and the Cetus group clearly contains many streams. As for the satellite–Cetus association, this is the first time that Willman 1 has been associated with this group (to the best of our knowledge). It could be that Willman 1, which is an ultrafaint dwarf galaxy (Willman et al. 2005), is actually the remnant of the progenitor Cetus galaxy (in other words, the remnant of the Cetus stream). This scenario is also supported by the fact that the [Fe/H] of Willman 1 (≈−2.1 dex; McConnachie & Venn 2020) is very similar to that of the Cetus stream (see Table 2).

In Figure 5, one can see two additional objects that lie close to the Cetus group, namely the globular cluster NGC 4590/M68 and the stream Fjörm (which is the stream produced from NGC 4590/M68). These two objects have very similar (*J*<sub>R</sub>, *J*<sub>φ</sub>, *E*) values to those of the Cetus group but possess lower *J*<sub>z</sub> values, rendering this association rather tentative. We note that NGC 4590 was previously associated with the Helmi substructure by Massari et al. (2019), Forbes (2020), and Kruijssen et al. (2020) and with the Canis Major progenitor galaxy (Martin et al. 2004) by Kruijssen et al. (2019). On the other hand, Fjörm was previously linked with Sagittarius by Bonaca et al. (2021).

<footnote n="18">¹⁸ These streams are Willka Yaku, Triangulum, and Turbio.</footnote>
<table><caption>**Table 3**
Various Groups Detected in Our Study along with Their Member Globular Clusters, Stellar Streams, and Satellite Galaxies</caption><thead><tr><th>Merger/ In Situ Group</th><th>No. of Members</th><th>Member Globular Clusters</th><th>Member Stellar Streams</th><th>Member Satellite Galaxies</th></tr></thead><tbody><tr><td>Sagittarius (Section 4.1)</td><td>8</td><td>Pal 12, Whiting 1, Terzan 7, Terzan 8, NGC 6715/M54, Arp 2</td><td>Elqui</td><td>Sagittarius dSph</td></tr><tr><td>Cetus (Section 4.2)</td><td>6–8</td><td>NGC 5824<br>tentative:<br>NGC 4590/M68 [stream:Fjörm]</td><td>Cetus [stream of Cetus],<br>Slidr, Atlas, AliqaUma<br>Fjörm [stream of NGC 4590/M68]</td><td>Willman 1</td></tr><tr><td>Gaia-Sausage/Enceladus (Section 4.3)</td><td>16–18</td><td>NGC 7492, NGC 6229, NGC 6584, NGC 5634, IC 1257, NGC 1851, NGC 2298, NGC 4147, NGC 1261, NGC 6981/M72, NGC 1904/M79, NGC 7089/M2 [stream:NGC 7089], NGC 5904/M5 [stream: M5]<br>tentative:<br>NGC 6864/M75</td><td>C-7, Hrid,<br>M 5 [stream of NGC 5904/M5]<br>NGC 7089 [stream of NGC 7089/M2]</td><td></td></tr><tr><td>Arjuna/Sequoia/I’itoi (Section 4.4)</td><td>9</td><td>NGC 6101,<br>NGC 3201 [streams: NGC 3201, Gjöll]</td><td>GD-1, Phlegethon, Gaia-9, Kshir,<br>NGC 3201 [stream of NGC 3201],<br>Gjöll [stream of NGC 3201], Ylgr</td><td></td></tr><tr><td>LMS-1/Wukong (Section 4.5)</td><td>11</td><td>NGC 5272/M3, NGC 5053, NGC 5024/M53, Pal 5 [stream: Pal 5]</td><td>LMS-1 [stream of LMS-1/Wukong],<br>C-19, Sylgr, Phoenix, Indus,<br>Jhelum, Pal 5 [stream of Pal 5]</td><td></td></tr><tr><td>Pontus (Section 4.6)</td><td>8–10</td><td>NGC 288, NGC 5286, NGC 7099/M30<br>NGC 6205/M13, NGC 6779/M56, NGC 6341/M92 [stream: M92], NGC 362<br>tentative:<br>NGC 6864/M75</td><td>M92 [stream of NGC 6341/M92]<br>NGC 7089 [stream of NGC 7089/M2]</td><td></td></tr><tr><td>Candidate merger (not detected, but selected, Section 5)</td><td>5</td><td>NGC 5466 [stream: NGC 5466], NGC 7492</td><td>Gaia-10,<br>NGC 5466 [stream of NGC 5466]</td><td>Tucana III</td></tr><tr><td>Galactic disk (Section 4.7)</td><td>6–7</td><td>Pal 10, NGC 6838/M71, NGC 6356, IC 1276/Pal 7, Pal 11, NGC 104/47 Tuc tentative:<br>NGC 7078/M15</td><td></td><td></td></tr><tr><td>Galactic Bulge (Section 4.7)</td><td>28</td><td>Terzan 2/HP 3, 1636-283/ESO452, Gran 1, Djorg 2/ESO 456, NGC 6453, NGC 6401, NGC 6304, NGC 6256, NGC 6325, Pal 6, Terzan 6/HP 5, Terzan 1/HP 2, NGC 6528, NGC 6522, NGC 6626/M28, Terzan 9, Terzan 5 11, NGC 6355, NGC 6638, NGC 6624, NGC 6266/M62, NGC 6642, NGC 6380/Ton1, NGC 6717/Pal9, NGC 6558, NGC 6342, HP 1/BH 229, NGC 6637/M69</td><td></td><td></td></tr><tr><td>Galactic Bulge/disk/low energy (Section 4.7)</td><td>11</td><td>Terzan 3, NGC 6569, NGC 6366, NGC 6139, BH 261/AL 3, NGC 6171/M107, Pal 8, Lynga 7/BH 184, NGC 6316, FSR 1716, NGC 6441</td><td></td><td></td></tr></tbody></table>

**Note.** This table is based on the associations that are visible in Figure 5. The detailed properties of these groups are described in Section 4. From left to right the columns provide the name of the group, total number of halo objects that are members of this group, names of the member globular clusters, names of the member streams, and names of the member satellite galaxies. In case of globular clusters, we provide in brackets the names of their streams (if any present in our sample). Likewise, in the case of stellar streams, we provide in brackets the names of their parent globular clusters (in case the parent globular cluster of the stream is known).

From Figure 6, it can be seen that all the Cetus member objects possess similar [Fe/H] values. We find that the MDF of this group has a median at [Fe/H] = −2.05 dex and spans a

very narrow range from −2.3 dex to −1.8 dex. Using the mass of NGC 5824, we estimate the mass of the progenitor Cetus galaxy as $M_{\text{halo}} \sim 2 \times 10^{10} M_{\odot}$, and this in turn provides
<figure><img src="image_13.png" /><figcaption>**Figure 6.** The orbital parameters of the halo objects as a function of their [Fe/H] values. We consider the [Fe/H] of only those objects that possess $P_{\text{Group}} \ge P_{\text{Threshold}}$, and the remaining objects are shown as gray points. The left plot shows $J_{\phi}$ vs. E, and the right plot shows $r_{\text{peri}}$ vs. $r_{\text{apo}}$. The LMS-1/Wukong group has a minimum of [Fe/H] = -3.4 dex.</figcaption></figure>

$M_* \sim 3 \times 10^7 M_{\odot}$. Note that these mass values likely represent a severe underestimation of the true (infall) mass of the progenitor Cetus galaxy, because this group contains many streams whose masses we have not accounted for (because we do not possess that information).

### 4.3. Group 3: The Gaia–Sausage/Enceladus Merger

This group represents the largest of all the mergers that we detect here. The member objects of this group possess relatively low values in $|J_{\phi}|$ and $L_{\perp}$, implying that they lie along radial orbits. This group possesses dynamical properties in the range $E \sim [-1.44, -1.16] \times 10^5 \text{ km}^2 \text{s}^{-2}$, $J_R \sim [935, 2075] \text{ km s}^{-1} \text{ kpc}$, $J_{\phi} \sim [-715, 705] \text{ km s}^{-1} \text{ kpc}$, $J_z \sim [85, 1505] \text{ km s}^{-1} \text{ kpc}$, $L_{\perp} \sim [85, 1520] \text{ km s}^{-1} \text{ kpc}$, eccentricity $\sim [0.7, 0.9]$, $r_{\text{peri}} \sim [1, 4] \text{ kpc}$, $r_{\text{apo}} \sim [16, 30] \text{ kpc}$, and $\phi \sim [27^{\circ}, 85^{\circ}]$. This group represents the Gaia–Sausage/Enceladus merger (Belokurov et al. 2018; Helmi et al. 2018; Myeong et al. 2018; Massari et al. 2019).

We find that 16 objects belong to this group: 3 streams (namely, C-7, M5, Hrid) and 13 globular clusters (namely, NGC 7492, NGC 6229, NGC 6584, NGC 5634, NGC 5904/M5, NGC 2298, NGC 4147, NGC 1261, NGC 6981/M72, NGC 7089/M2, IC 1257, NGC 1904/M79, NGC 1851). There exist two additional objects close to this group, namely the globular cluster NGC 6864/M75 and the stream NGC 7089, but their association is not very strong (because of their slightly lower $J_R$ values).

These streams–Gaia–Sausage/Enceladus associations are reported here for the first time. Unlike Bonaca et al. (2021), we do not find the streams Ophiuchus and Fimbulthul to be members of this group. As for the globular-cluster–Gaia–Sausage/Enceladus associations, our list contains half of those 10 clusters that were previously associated with this group by Myeong et al. (2018). However, more recent studies have attributed a large number of globular clusters to the Gaia–Sausage/Enceladus merger. For instance, Massari et al. (2019) associated $\approx 32$ globular clusters to this merger, although some

of their associations were tentative. They found these association by making hard cuts in the $(J_{\phi}, E, L_{\perp})$ space that were previously used by Helmi et al. (2018) to select the Gaia-Sausage/Enceladus stellar debris. Massari et al. (2019) further supported their associations by arguing that the resulting globular clusters show a tight age–metallicity relation (AMR). We show the AMR of our Gaia–Sausage/Enceladus globular clusters in Figure 9 that (visually) appears to be tighter than Figure 4 of Massari et al. (2019). The study of Forbes (2020), which is based on the analysis of Massari et al. (2019), found 28 globular-cluster–Gaia–Sausage/Enceladus associations. We find that some of these additional globular clusters, which have recently been linked with Gaia–Sausage/Enceladus by other studies, likely belong to a different merger group (see Section 4.6).

The halo objects associated with the Gaia–Sausage/Enceladus merger span a very wide range in [Fe/H] from $-2.4$ dex to $-1.1$ dex, with the median of the MDF located at [Fe/H] = $-1.6$. This large spread in MDF supports the scenario that Gaia–Sausage/Enceladus was a massive galaxy. We estimate the mass of the progenitor galaxy to be $M_{\text{halo}} \sim 10 \times 10^{10} M_{\odot}$ and $M_* \sim 50 \times 10^7 M_{\odot}$. This mass estimate is consistent with those found by previous studies from chemical evolution models (e.g., Fernández-Alvar et al. 2018; Helmi et al. 2018), counts of metal-poor and highly eccentric stars (e.g., Mackereth & Bovy 2020), and the mass–metallicity relation (e.g., Naidu et al. 2020).

### 4.4. Group 4: The Arjuna/Sequoia/I’itoi Merger

This group is highly retrograde and its member objects possess dynamical properties in the range $E \sim [-1.27, -1.02] \times 10^5 \text{ km}^2 \text{s}^{-2}$, $J_R \sim [20, 1190] \text{ km s}^{-1} \text{ kpc}$, $J_{\phi} \sim [1880, 2955] \text{ km s}^{-1} \text{ kpc}$, $J_z \sim [230, 940] \text{ km s}^{-1} \text{ kpc}$, $L_{\perp} \sim [980, 2530] \text{ km s}^{-1} \text{ kpc}$, eccentricity $\sim [0.1, 0.6]$, $r_{\text{peri}} \sim [5, 14] \text{ kpc}$, $r_{\text{apo}} \sim [15, 37] \text{ kpc}$, and $\phi \sim [25^{\circ}, 46^{\circ}]$. We refer to this group as the Arjuna/Sequoia/I’itoi group, because it likely comprises objects that actually resulted from three independent mergers: Sequoia (Myeong et al. 2019), Arjuna,
<figure><img src="image_14.png" /><figcaption>**Figure 7.** The metallicity distribution function (MDF) of different groups detected in our study. This MDF is constructed using the [Fe/H] measurements of globular clusters and streams that belong to different groups. The LMS-1/Wukong group has a minimum of [Fe/H] = −3.4 dex.</figcaption></figure>

and I'toi (Naidu et al. 2020). This understanding comes from Naidu et al. (2020), who performed a chemodynamical analysis of stars and proposed that at this (E, J<sub>φ</sub>) location, there exist three different (but somewhat overlapping) stellar populations: a metal-rich population whose MDF peaks at [Fe/H] ≈ −1.2 dex (namely Arjuna), another one whose MDF peaks at [Fe/H] ≈ −1.6 dex (namely Sequoia), and the most metal-poor among these whose MDF peaks at [Fe/H] ≈ −2 dex (namely I'toi). Here, we analyze this detected group as a single merger because our detected grouping contains only a handful of objects and therefore it is difficult to detect

any plausible subgroups within this Arjuna/Sequoia/I'toi group.

We find that nine objects belong to this group: seven streams (namely, Phlegethon, Gaia-9, NGC 3201, Gjöll, GD-1, Kshir, and Ylgr) and two globular clusters (namely, NGC 6101 and NGC 3201). These two globular clusters were previously associated with Sequoia by Myeong et al. (2019), although they associated five additional clusters with this group that we do not identify here. To discover the Sequoia group, Myeong et al. (2019) applied a “Friends-of-Friends” grouping algorithm to the projected action space containing only globular clusters (essentially, they applied their algorithm to the Gaia DR2 version of the top-left panel shown in Figure 2). With this, they found a group of globular clusters (which they named Sequoia) whose combined dynamical properties ranged from E = [−2.2, −0.97] × 10<sup>5</sup> km<sup>2</sup> s<sup>−2</sup>, J<sub>R</sub> ~ [54, 1400] km s<sup>−1</sup> kpc, J<sub>φ</sub> = [250, 3210] km s<sup>−1</sup> kpc, and J<sub>z</sub> = [66, 800] km s<sup>−1</sup> kpc (they used the same Galactic potential model as ours). This dynamical range is larger than the range we infer for the Arjuna/Sequoia/I'toi group, especially in J<sub>φ</sub> and E. Moreover, it could be due to their wide E selection that even the low-energy cluster NGC 6401 ends up in their Sequoia group; we note that NGC 6401 possesses such a low energy (Massari et al. 2019) that it likely belongs to the Galactic bulge (see below). On the other hand, our two Arjuna/Sequoia/I'toi globular clusters were previously associated with both Sequoia and Gaia–Sausage/Enceladus by Massari et al. (2019); we note that their (J<sub>φ</sub>, E) selection is motivated by the results of Myeong et al. (2019). But Massari et al. (2019) also found five additional member clusters for Sequoia, and many of these are not present in the Myeong et al. (2019) selection. The Sequoia group found by Forbes (2020) is very similar to that of Massari et al. (2019), likely because the selection of the former study is based on the latter. As for the stream–Arjuna/Sequoia/I’toi associations, Myeong et al. (2019) analyzed the (J, E) of only GD-1 and argued against its association. However, Bonaca et al. (2021) favored an association of GD-1 with Arjuna/Sequoia/I’toi, along with those of Phlegethon, Gjöll, and Ylgr.

The MDF of this group spans a wide range from −2.24 dex to −1.56 dex, with the median of [Fe/H] ~ −1.78 dex. Interestingly, this [Fe/H] median is similar to the [Fe/H] of the Kshir stream (Malhan et al. 2019a). Kshir is a broad stream that moves in the Milky Way along a very similar orbit to that of GD-1, and this observation encouraged Malhan et al. (2019a) to propose that Kshir is likely the stellar stream produced from the tidal stripping of the merging galaxy that brought in GD-1. If true, this indicates that Kshir is likely the stream of the progenitor Arjuna/Sequoia/I’toi galaxy, perhaps that of the Sequoia galaxy (given the similarity in their [Fe/H]).

Using only the member globular clusters of Arjuna/Sequoia/I'toi and not the streams, we estimate the mass of the progenitor galaxy to be M<sub>halo</sub> ~ 1.2 × 10<sup>10</sup> M<sub>☉</sub> and M<sub>*</sub> ~ 1.5 × 10<sup>7</sup> M<sub>☉</sub>. Note that the actual masses should be higher than these computed masses because this group contains several streams (as compared to globular clusters) whose masses we could not account for. Interestingly, these mass values are similar to those derived by Myeong et al. (2019) for the Sequoia merger, using similar techniques, even though we could not identify many of their globular clusters as members of our Arjuna/Sequoia/I'toi group.
<figure><img src="image_15.png" /><figcaption>**Figure 8.** Comparing the orbital properties of those objects that belong to different groups. The left panel shows $r_{\text{peri}}$ vs. eccentricity, and the right panel shows $r_{\text{peri}}$ vs. $r_{\text{apo}}$. We use the same color for all those objects that belong to the same group. The gray points represent all the objects in our sample. The objects that form part of the same group are clumped in this space (in addition to being clumped in (J, E) space; see Figure 5).</figcaption></figure>

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

### 4.5. Group 5: The LMS-1/Wukong Merger

This group has a slight prograde motion and its member objects are very tightly clumped in (J, E) space. It possesses dynamical properties in the range $E \sim [-1.41, -1.19] \times 10^5 \text{ km}^2 \text{s}^{-2}$, $J_R \sim [100, 605] \text{ km s}^{-1} \text{kpc}$, $J_\phi \sim [-1560, -210] \text{ km s}^{-1} \text{kpc}$, $J_z \sim [875, 2710] \text{ km s}^{-1} \text{kpc}$, $L_\perp \sim [1400, 3085] \text{ km s}^{-1} \text{kpc}$, eccentricity $\sim [0.2, 0.5]$, $r_{\text{peri}} \sim [5, 13] \text{ kpc}$, $r_{\text{apo}} \sim [15, 25] \text{ kpc}$, and $\phi \sim [58^\circ, 85^\circ]$. This polar group corresponds to the low-mass-stream-1 (LMS-1)/Wukong merger (Yuan et al. 2020a; Naidu et al. 2020; Malhan et al. 2021).

We find that 11 objects belong to this group: 7 streams (LMS-1 itself, Phoenix, Pal 5, C-19, Indus, Sylgr, and Jhelum) and 4 globular clusters (namely NGC 5272/M3, NGC 5053, Pal 5, and NGC 5024/M53). With regard to the stream–LMS-1/Wukong associations, Phoenix, Indus, Jhelum, and Sylgr were tentatively associated with this merger by Bonaca et al. (2021). The association of Indus was also favored by Malhan et al. (2021); however, this study had argued against Jhelum’s association. Our result here could be different from Malhan et al. (2021) because here we are using different data for streams; that consequentially results in different (J, E) solutions. Furthermore, because Indus and Jhelum are tidal debris of dwarf galaxies (Li et al. 2021a), this indicates that they were likely the satellite dwarf galaxies of the progenitor LMS-1/Wukong galaxy (also see Malhan et al. 2021). As for the globular cluster–LMS-1/Wukong associations, Koppelman et al. (2019b) and Massari et al. (2019) associated NGC 5024, NGC 5053, and NGC 5272 with the Helmi substructure (Helmi et al. 1999). Koppelman et al. (2019b), in particular, supported the association of four additional clusters with the Helmi substructure, but we find many of these clusters as part of the Gaia–Sausage/Enceladus group. As for Massari et al. (2019), they could only associate the clusters with those merger groups that were known at that time, but LMS-1/Wukong was detected after their study by Yuan et al. (2020a) and Naidu et al. (2020). On the other hand, recent studies have shown that there indeed exists a strong association of NGC 5024 and NGC 5053 with the LMS-1/Wukong group (Yuan et al. 2020a; Naidu et al. 2020; Malhan et al. 2021), based on the fact that these clusters lie within the phase-space distribution of the LMS-1 stream

(e.g., Yuan et al. 2020a; Malhan et al. 2021). Another recent study by Wan et al. (2020) advocates for a dynamical connection between Phoenix, Pal 5, and NGC 5053. In summary, our analysis supports these recent studies and makes a stronger case that all of these objects are associated with the LMS-1/Wukong merger.

We find that LMS-1/Wukong is the most metal-poor merger of the Milky Way because this group contains the three most metal-poor streams of our galaxy, namely C-19, Sylgr, and Phoenix (see their [Fe/H] values in Table 2). Overall, this group has a wide MDF ranging from −3.38 dex to −1.41 dex with the median of [Fe/H] $\sim -2$ dex. We note that this median is similar to the metallicity of the LMS-1 stream (Malhan et al. 2021). Using the masses of the globular clusters, we estimate the progenitor galaxy’s mass to be $M_{\text{halo}} \sim 2.7 \times 10^{10} \text{M}_\odot$ and $M_* \sim 5.5 \times 10^{10} \text{M}_\odot$. These mass values are higher than those reported in Malhan et al. (2021) because here we find a higher number of globular-cluster–LMS-1/Wukong associations.

### 4.6. Group 6: Discovery of the Pontus Merger

We detect a new group that possesses low energy and is slightly retrograde. Its dynamical properties are in the range $E \sim [-1.72, -1.56] \times 10^5 \text{ km}^2 \text{s}^{-2}$, $J_R \sim [245, 725] \text{ km s}^{-1} \text{kpc}$, $J_\phi \sim [-5, 470] \text{ km s}^{-1} \text{kpc}$, $J_z \sim [115, 545] \text{ km s}^{-1} \text{kpc}$, $L_\perp \sim [390, 865] \text{ km s}^{-1} \text{kpc}$, eccentricity $\sim [0.5, 0.8]$, $r_{\text{peri}} \sim [1, 3] \text{ kpc}$, $r_{\text{apo}} \sim [8, 13] \text{ kpc}$, and $\phi \sim [54^\circ, 89^\circ]$. We refer to this group as Pontus.<fnref n="19" /> We find that eight objects belong in this group: one stream (namely, M92) and seven clusters (namely, NGC 288, NGC 5286, NGC 7099/M30, NGC 6205/M13, NGC 6341/M92, NGC 6779/M56, and NGC 362). There exist two additional objects close to this group, namely the globular cluster NGC 6864/M75 and the stream NGC 7089, but their association was not very strong (because of their slightly higher $J_R$ and slightly lower $J_\phi$ values).

Pontus lies close to Gaia–Sausage/Enceladus in (J, E) space (although the two groups possess very different $J_R$ values) and essentially all of Pontus’s globular clusters (which we mentioned above) have been previously associated with

<footnote n="19">¹⁹ In Greek mythology, “Pontus” (meaning “the Sea”) is the name of one of the first children of the deity Gaia.</footnote>
<figure><img src="image_17.png" /><figcaption>**Figure 9.** Comparing properties of those objects that belong to the groups Gaia-Sausage/Enceladus and Pontus. The left panel shows the age-metallicity relationship of the globular clusters belonging to these groups; the [Fe/H] and age values are taken from Kruijssen et al. (2019). The right panel shows the velocity behavior of these objects in spherical polar coordinates, namely radial $V_r$ and azimuthal $V_\phi$. The filled ellipses represent the 1.5$\sigma$ confidence contour and the Gaussians represent the mean and the standard deviation in the $V_\phi$ components of the member objects of these groups.</figcaption></figure>

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

Gaia–Sausage/Enceladus(Massari et al. 2019). Given this potential overlapping between the two groups, it is natural to ask: Do these groups in fact represent different merging events, or is it that our procedure has fragmented the large Gaia–Sausage/Enceladus group into two pieces? We argue that fragmentation cannot be the reason, otherwise the neighboring Sagittarius and Cetus groups should also be regarded as a single group, as these latter groups are much closer to each other in (J, E) space compared to the former groups. Similarly, even the neighboring LMS-1/Wukong and Gaia–Sausage/Enceladus groups could be distinguished by our detection procedure. To understand the nature of Pontus and Gaia–Sausage/Enceladus groups, we use their member objects and analyze their dynamical properties and the AMR.

First, we find that the objects belonging to these two groups possess different dynamical properties. The average eccentricity of Pontus objects is smaller than that of the Gaia-Sausage/Enceladus objects (see Figure 8), while we note that the eccentricity range of our Gaia-Sausage/Enceladus group is similar to that of Myeong et al. (2018). This implies that the orbits of Gaia-Sausage/Enceladus objects are more radial than those of the Pontus objects (this can also be discerned by comparing their $J_R$ values). Also, the average $r_{apo}$ of Pontus objects is smaller than that of the Gaia-Sausage/Enceladus objects. Furthermore, we also compare the velocity behavior of their member objects in spherical polar coordinates, namely radial $V_r$ and azimuthal $V_\phi$ (see Figure 9). The motivation for adopting this particular coordinate system comes from Belokurov et al. (2018), who used a similar system to originally identify the “sausage” structure. From Figure 9, we note that both the Pontus and Gaia-Sausage/Enceladus distributions are stretched along the $V_r$ direction (implying radial orbits), although their $V_\phi$ components (on average) differ by ≈60 km s⁻¹. This implies, as noted above, that Pontus objects are more retrograde than the Gaia-Sausage/Enceladus objects. Also, Gaia-Sausage/Enceladus objects possess a larger dispersion in $V_\phi$ compared to Pontus objects.

Moreover, in Figure 9, the AMR for the globular clusters belonging to these two groups also appear quite different, especially the age difference of their metal-poor clusters (⟨−1.5 dex) is ≳1 Gyr. In view of this investigation, we conclude that Pontus and Gaia-Sausage/Enceladus represent two distinct and independent merging events: Gaia-Sausage/Enceladus comprising slightly younger globular clusters than those present in Pontus and Gaia-Sausage/Enceladus’s objects possessing overall different dynamical properties compared to Pontus’s objects.

Similarly, we argue that the Pontus group is also different than the Thamnos substructures identified by Koppelman et al. (2019a). Koppelman et al. (2019a) suggested that at the location (E, $J_\phi$) ~ (−1.65 × 10⁵ km² s⁻², 1500 km s⁻¹ kpc) and ~ (−1.75 × 10⁵ km² s⁻², 900 km s⁻¹ kpc), there lie two substructures, namely Thamnos 1 and Thamnos 2 (see their Figure 2). Motivated by their selection, Naidu et al. (2020) selected Thamnos stars around a small region of (E, $J_\phi$) ~ (−1.75 × 10⁵ km² s⁻², 500 km s⁻¹ kpc) (see their Figure 23). Given that these (E, $J_\phi$) locations for Thamnos are different from those of Pontus, we argue that Pontus is independent of Thamnos. Moreover, the metallicity of Thamnos 2 members is different from that of Pontus members. This we argue by inspecting Figure 2 of Koppelman et al. (2019a), which shows that the metallicity of Thamnos 2 stars range from [Fe/H] ~ −1.4 dex to −1.1 dex, and this is different from the metallicity range of Pontus (see below). We also note that a few of the Pontus member clusters were previously tentatively associated with the Canis Major progenitor galaxy (Kruijssen et al. 2019).

The MDF of Pontus spans a range from −2.3 dex to −1.3 dex with a median of [Fe/H] = −1.7 dex. We estimate the mass of the progenitor Pontus galaxy to be $M_{halo}$ ~ 5 × 10¹⁰ M☉ and $M_*$ ~ 15 × 10⁷ M☉.

### 4.7. The in Situ Groups 7, 8, and 9

We detect three additional groups; however, their locations in (E, $J_\phi$) space indicates that they do not represent any merger
but actually belong to the in situ population of the Milky Way—the population of the Galactic disk and the Galactic bulge. This we infer based on the fact that the member objects of these groups possess low E, low $r_{apo}$, and high [Fe/H]—as expected from the in situ globular cluster population (see Figure 6).

The first of these groups possess dynamical properties in the range $E \sim [-1.77, -1.66] \times 10^5 \text{ km}^2 \text{s}^{-2}$, $J_R \sim [15, 135] \text{ km s}^{-1} \text{kpc}$, $J_\phi \sim [-1315, -960] \text{ km s}^{-1} \text{kpc}$, $J_z \sim [5, 235] \text{ km s}^{-1} \text{kpc}$, $L_\perp \sim [115, 730] \text{ km s}^{-1} \text{kpc}$, eccentricity $\sim [0.1, 0.4]$, $r_{\text{peri}} \sim [3, 6] \text{ kpc}$, $r_{\text{apo}} \sim [7, 9] \text{ kpc}$, and $\phi \sim [5^\circ, 36^\circ]$. Given these dynamical properties, especially low eccentricity (implying circular orbits), low $\phi$ value (implying that the objects orbit close to the Galactic plane), and the values of $r_{\text{peri}}$ and $r_{\text{apo}}$ being similar to stars in the Galactic disk, we interpret this as the disk group. This group contains six globular clusters (their names are provided in Table 3). There exists one additional cluster, NGC 7078/M15, that lies close to this group in (J, E) space but we do not identify it as a strong associate. The member globular clusters are metal rich and the corresponding MDF ranges from −0.8 dex to −0.1 dex with a median of [Fe/H] $\sim -0.65$ dex (see Figure 7). It is interesting to note that this MDF minimum is consistent with the results of Zinn (1985), who found [Fe/H] $\ge -0.8$ as the threshold between the disk and halo clusters (they inferred this simply on the basis of the bimodality of the [Fe/H] distribution of the globular clusters). All of our globular clusters, including the very metal-poor NGC 7078/M15, were previously associated with the disk by Massari et al. (2019); although they associated a total of 26 clusters to the disk, we could not identify all of these objects.

The second group possesses the lowest energy among all the detected groups, implying that its member objects orbit deep in the potential of the Milky Way—close to the Galactic center. The members of this group possess the dynamical properties in the range $E \sim [-2.55, -2.22] \times 10^5 \text{ km}^2 \text{s}^{-2}$, $J_R \sim [5, 140] \text{ km s}^{-1} \text{kpc}$, $J_\phi \sim [-380, 150] \text{ km s}^{-1} \text{kpc}$, $J_z \sim [0, 275] \text{ km s}^{-1} \text{kpc}$, $L_\perp \sim [10, 405] \text{ km s}^{-1} \text{kpc}$, eccentricity $\sim [0.1, 0.8]$, $r_{\text{peri}} \sim [0, 2] \text{ kpc}$, $r_{\text{apo}} \sim [1, 4] \text{ kpc}$, and $\phi \sim [1^\circ, 89^\circ]$. This group comprises 28 globular clusters (their names are provided in Table 3). Given these dynamical properties, especially very low values of E, $r_{\text{peri}}$, and $r_{\text{apo}}$, and that the objects are spherically distributed (as we note from the range of the $\phi$ parameter), we interpret this as the Galactic bulge group. The member objects span a wide range in [Fe/H], ranging from −1.5 dex to −0.1 dex with a median at [Fe/H] $\sim -1.0$ dex. We confirm that several of these objects have been associated with the Galactic bulge by Massari et al. (2019), although they associated a total of 36 globular clusters to the bulge. A few of our member objects were interpreted by Massari et al. (2019) as “unassociated objects with low energy,” but Forbes (2020) interpreted these objects as those accreted inside the Koala progenitor galaxy. We further note that some of our clusters have also been interpreted as the bulge objects by Horta et al. (2020) on the basis of their high alpha-element abundances and high [Fe/H] values.

We detect a third group that is slightly prograde and possesses E values between that of the bulge and disk groups (see Figure 5). Its dynamical properties lie in the range $E \sim [-2.17, -1.92] \times 10^5 \text{ km}^2 \text{s}^{-2}$, $J_R \sim [10, 110] \text{ km s}^{-1} \text{kpc}$, $J_\phi \sim [-630, -250] \text{ km s}^{-1} \text{kpc}$, $J_z \sim [55, 190] \text{ km s}^{-1} \text{kpc}$, $L_\perp \sim [215, 485] \text{ km s}^{-1} \text{kpc}$, eccentricity $\sim [0.2, 0.4]$, $r_{\text{peri}} \sim [1, 3] \text{ kpc}$, $r_{\text{apo}} \sim [3, 6] \text{ kpc}$, and $\phi \sim [18^\circ, 56^\circ]$. This group contains 11 globular clusters (these are listed in Table 3). The metallicity of these objects range from [Fe/H] $\sim -1.65$ dex to −0.4 dex with a

median of [Fe/H] = −0.7 dex. For this group, while its dynamical properties appear consistent with those of the disk (e.g., low eccentricity, the $r_{apo}$ range, and low $\phi$ values), its relatively lower [Fe/H] value appears more consistent with that of the bulge. This makes it challenging to associate these objects with either disk or bulge. Perhaps Massari et al. (2019) were also in a similar conundrum that they interpreted some of these objects as the bulge clusters, some as disk clusters, and others simply as “low-energy objects.” Among our member objects, NGC 6441 was tentatively associated with the Kraken merger by Kruijssen et al. (2020). We argue that our objects likely do not belong to Kraken because our objects possess slightly negative $J_\phi$ values (on average), while Kraken objects have an average $J_\phi \sim 0$ (that we observed from Figure 17 in Appendix C). Also, a few of our objects have been interpreted as either bulge or simply low-energy clusters by Horta et al. (2020) on the basis of their chemical compositions.

## 5. A Candidate Merger

All the above-mentioned groups were detected at $\ge 2\sigma$ significance following the ENLINK procedure described in Section 3. However, during our multiple ENLINK runs (while we were initially experimenting with different parameters), we noticed a particular group that comprised five objects whose $P_{\text{Group}}$ values were fluctuating close to the detection threshold. This is likely due to the ENLINK parameter min_cluster_size that we set to 5 (see Section 3.2), thus making it difficult for ENLINK to detect groups containing $\le 5$ objects. Motivated by the possibility that this group may represent an actual merger, we discuss its properties below.

This group possesses a slight retrograde motion and relatively high energy (as shown in Figure 10). The dynamical properties of this group lie in the range $E \sim [-1.2, -0.94] \times 10^5 \text{ km}^2 \text{s}^{-2}$, $J_R \sim [1385, 2525] \text{ km s}^{-1} \text{kpc}$, $J_\phi \sim [130, 915] \text{ km s}^{-1} \text{kpc}$, $J_z \sim [1125, 2095] \text{ km s}^{-1} \text{kpc}$, $L_\perp \sim [1435, 2795] \text{ km s}^{-1} \text{kpc}$, eccentricity $\sim [0.7, 0.8]$, $r_{\text{peri}} \sim [3, 7] \text{ kpc}$, $r_{\text{apo}} \sim [27, 47] \text{ kpc}$, and $\phi \sim [69^\circ, 85^\circ]$. The group comprises two globular clusters (namely, NGC 5466 and NGC 7492), two stellar streams (NGC 5466 and Gaia-10), and one dwarf galaxy (Tucana III). The $\phi$ parameter indicates that the member objects have very “polar” orbits. Particularly for Tucana III, Gaia-10, and NGC 5466, we note that their orbital planes are very similar; this further lends credence to their possible association. The MDF of this group ranges from [Fe/H] = −2.4 dex to −1.4 dex. These minima and maxima are set by Tucana III (Simon et al. 2017) and Gaia-10, respectively. NGC 5466 has [Fe/H] $\sim -1.98$ dex (Lamb et al. 2015) and NGC 7492 has [Fe/H] $\sim -1.8$ dex (Cohen & Melendez 2005). We further compare the stellar population of these objects in terms of their color-magnitude distributions (CMDs), and this is shown in Figure 11. These objects possess strikingly similar CMDs, despite their differences in [Fe/H].<fnref n="20" /> In summary, the similarities in the stellar population of these objects, together with their coincidence in (J, E) space, indicate that these objects were perhaps born at the same time inside the same progenitor galaxy.

Previous studies have associated NGC 5466 and NGC 7492 with the Sequoia and Gaia–Sausage/Enceladus groups, respectively (Massari et al. 2019; Forbes 2020). On the other

<footnote n="20">²⁰ The reason that Tucana III’s CMD appears scattered is because we construct the CMD using the photometry from Gaia EDR3, and Gaia has a limiting magnitude at $G \sim 21$.</footnote>
<table><caption>**Table 4**
Actions, Energies, Orbital Parameters, and [Fe/H] of Those Globular Clusters/Satellite Galaxies that are Physically Connected to Stellar Streams (See Section 6)</caption><thead><tr><th>Globular Cluster/<br>Satellite Galaxy</th><th>(J<sub>R</sub>, J<sub>φ</sub>, J<sub>z</sub>)<br>(kpc km s<sup>-1</sup>)</th><th>Energy<br>(km<sup>2</sup> s<sup>-2</sup>)</th><th>r<sub>peri</sub><br>(kpc)</th><th>r<sub>apo</sub><br>(kpc)</th><th>z<sub>max</sub><br>(kpc)</th><th>Eccentricity</th><th>[Fe/H]<br>(dex)</th></tr></thead><tbody><tr><td>NGC 3201</td><td>(787<sup>+29</sup><sub>-27</sub>, 2724<sup>+21</sup><sub>-18</sub>, 276<sup>+5</sup><sub>-4</sub>)</td><td>-115824<sup>+752</sup><sub>-724</sub></td><td>8.3<sup>+0.0</sup><sub>-0.0</sub></td><td>27.1<sup>+0.5</sup><sub>-0.5</sub></td><td>10.8<sup>+0.2</sup><sub>-0.2</sub></td><td>0.53<sup>+0.01</sup><sub>-0.01</sub></td><td>-1.59</td></tr><tr><td>NGC 4590</td><td>(765<sup>+23</sup><sub>-25</sub>, -2359<sup>+14</sup><sub>-15</sub>, 824<sup>+12</sup><sub>-12</sub>)</td><td>-113703<sup>+611</sup><sub>-668</sub></td><td>9.1<sup>+0.1</sup><sub>-0.0</sub></td><td>28.1<sup>+0.4</sup><sub>-0.5</sub></td><td>18.5<sup>+0.3</sup><sub>-0.3</sub></td><td>0.51<sup>+0.0</sup><sub>-0.01</sub></td><td>-2.23</td></tr><tr><td>NGC 5139</td><td>(157<sup>+5</sup><sub>-5</sub>, 540<sup>+10</sup><sub>-10</sub>, 136<sup>+6</sup><sub>-6</sub>)</td><td>-185068<sup>+253</sup><sub>-258</sub></td><td>2.7<sup>+0.1</sup><sub>-0.1</sub></td><td>6.9<sup>+0.0</sup><sub>-0.0</sub></td><td>3.0<sup>+0.1</sup><sub>-0.1</sub></td><td>0.44<sup>+0.01</sup><sub>-0.01</sub></td><td>-1.53</td></tr><tr><td>Kim 3</td><td>(623<sup>+4925</sup><sub>-526</sub>, -569<sup>+800</sup><sub>-1117</sub>, 2095<sup>+3421</sup><sub>-1872</sub>)</td><td>-117699<sup>+58531</sup><sub>-40035</sub></td><td>9.4<sup>+5.7</sup><sub>-3.7</sub></td><td>24.5<sup>+90.9</sup><sub>-14.0</sub></td><td>21.8<sup>+91.7</sup><sub>-16.6</sub></td><td>0.5<sup>+0.29</sup><sub>-0.24</sub></td><td>-1.6<sup>+0.45</sup><sub>-0.30</sub></td></tr><tr><td>Willman 1</td><td>(837<sup>+306</sup><sub>-262</sub>, -2432<sup>+758</sup><sub>-490</sub>, 2726<sup>+1775</sup><sub>-892</sub>)</td><td>-93587<sup>+8055</sup><sub>-7180</sub></td><td>16.0<sup>+4.4</sup><sub>-2.6</sub></td><td>42.4<sup>+7.1</sup><sub>-6.9</sub></td><td>37.3<sup>+10.1</sup><sub>-8.6</sub></td><td>0.44<sup>+0.07</sup><sub>-0.08</sub></td><td>-2.1</td></tr></tbody></table>

Note. The orbital parameters are derived in this study and the values represent the medians of the sampled posterior distributions and the corresponding uncertainties reflect their 16th and 84th percentiles. The [Fe/H] measurements are taken from the literature. We take the [Fe/H] of all the globular clusters from the Harris (2010) catalog, except for Kim 3 (we use Kim et al. 2016). The [Fe/H] value of the satellite galaxies are taken from the McConnachie & Venn (2020) catalog.

hand, Gaia-10 (which is being analyzed here for the first time) and Tucana III have not been previously associated with any merger group.

From Figure 10, one notices that the globular cluster NGC 5466 and its stream have different values of (J, E), although they should ideally have very similar orbits. Upon closer inspection we find that the computed orbit of the NGC 5466 stream fits the phase-space data nicely; however, its orbital distance solution has an offset of ~1.5 kpc from that of the globular cluster. This offset in distance changes the resulting (J, E) solutions. This offset arises because we are constraining the distance of the streams using only the parallaxes (see Section 2.2), but the uncertainties on the parallaxes are too large to properly constrain the strong distance gradient that is present particularly in this stream (Jensen et al. 2021).

## 6. Physical Connections between Streams and Other Objects

Until this point in the analysis, the question that we have tried to answer is: Which set of globular clusters, streams, and satellites were accreted inside which merger galaxies? In this section, our objective is slightly different from the main objective of this paper. Here, we want to investigate: Which stream was produced by the tidal stripping of which parent system (i.e., which globular cluster or satellite galaxy)? At this juncture of our analysis, exploring this question is relatively straightforward because we possess the orbital information of all halo objects. Therefore, all we require to do is to compare the orbits of streams with those of globular clusters and satellites to find plausible physical connections among them.<fnref n="21" /> Below, we first briefly describe the motivation behind this investigation and then explain our method of identifying these physical connections.

A stellar stream is generally pictured as two tidal tails emerging from its parent system (be it a globular cluster or a satellite galaxy, e.g., Dehnen et al. 2004). While this is true for a few streams in our sample (e.g., Pal 5, Fjörm, NGC 3201, NGC 5466), a majority of them are observed without their parent systems (≈75% of the streams considered in this study). The two main reasons for this could be (1) if the parent system

of a stream has been completely dissolved due to tidal stripping, then it cannot be observed at the present day, and (2) if the parent system is spatially separated from the stream, then it becomes difficult to recognize any physical connection between the two. The latter scenario is especially possible for those streams that merged inside their progenitor galaxies because mergers can deposit their stellar content in disparate regions of phase space (e.g., Jean-Baptiste et al. 2017). It is this second scenario that we want to explore here. Note that this knowledge of a stream's parent object is also crucial for the stream's N-body dynamical modeling (e.g., Thomas et al. 2016; Bonaca et al. 2019).

Our strategy to find parent objects of streams is straightforward. We compare the (J, E) and [Fe/H] values of streams with those of globular clusters and satellites. Thus, a stream is considered to be physically associated with an object if (1) their (J, E) values differ by <2σ, and (2) their [Fe/H] values are similar. Through this (J, E)–[Fe/H] comparison, we identify four stream–globular-cluster pairs and one stream–satellite pair. The former pairs include Gjöll–NGC 3201 (originally noted in Hansen et al. 2020; Palau & Miralda-Escude 2021), Fimbulthul–NGC 5139/ω Centauri (Ibata et al. 2019a), Fjörm–NGC 4590/M 68 (Palau & Miralda-Escude 2019), and Ophiuchus–Kim 3 (we note that the large dispersion in the (J, E) values of Kim 3 renders this association tentative). The latter pair includes Cetus–Willman 1. The details of these parent objects are provided in Table 4, and these can be compared with the parameters of the streams in Table 2.

We find additional associations based only on the (J, E) values because a majority of streams lack [Fe/H] measurements. In this case, we consider only those associations for which the (J, E) values differ by <1σ. These pairs include C3–NGC 288.

We also examine if there are any streams that are possibly connected to other streams, implying that these physically separated structures actually represent different segments of the same stream. To this end, we consider only those associations for which the (J, E) values differ by <1σ and find the following associations: Atlas-AliqaUma (originally noted in Li et al. 2021b), C3–Gaia-6, and Fimbulthul–Gunnthrà. This last association suggests that Gunnthrà is the leading tidal arm of the ω Centauri cluster and Fimbulthul is the trailing arm; in a way, we complete the picture proposed by Ibata et al. (2019a).

## 7. Discussion and Conclusion

We implemented an objective search strategy and detected N = 6 apparently independent mergers that our galaxy has

<footnote n="21">²¹ This is a very simple method of finding parent systems of streams, and it does not guarantee that the resulting connections represent the reality (at least not in all the cases). For confirmation, such methods should be complemented with the [Fe/H] comparison of the stream and also proper N-body modeling; however, these are beyond the scope of this paper. It is for these reasons that we proceeded in Section 3 with both the globular cluster and stellar stream counterparts, as we deemed it less biased.</footnote>
<figure><img src="image_19.png" /><figcaption>**Figure 10.** (*J*, *E*) distribution of the objects belonging to the candidate merger (see Section 5). These objects are shown along with all the other objects in our sample. Similar to Figure 2, each object is shown as a “cloud” of 1000 Monte Carlo representations of its orbit.</figcaption></figure>

suffered during its long history. We achieve this by searching for statistically significant groups of halo objects in (*J*, *E*) space. Our data set comprises *n* = 257 objects (namely, 170 globular clusters, 41 stellar streams, and 46 satellite galaxies). The (*J*, *E*) values of these objects are derived using Gaia EDR3-based catalogs. To detect the mergers in (*J*, *E*) space, we used the ENLINK group-finding software, coupled with our statistical procedure that accounts for the uncertainties on the derived (*J*, *E*) values of independent objects (that arises due to the measurement uncertainties on the phase-space properties of the objects). All of the mergers are detected at ≥2σ confidence, and together, they comprise 62 objects (25% of the total sample size), including 35 globular clusters, 25 streams, and 2 satellite galaxies.

We successfully recovered many of the previously known mergers (namely, Sagittarius, Cetus, Gaia–Sausage/Enceladus, Arjuna/Sequoia/Γ’itoi, and LMS-1/Wukong).<fnref n="22" /> Moreover, we discovered a new merger (which we call Pontus) and a possible merger candidate. Section 4 details the overall properties of these mergers—their member objects (Table 3), their (*J*, *E*) distribution (Figure 5), their orbital properties as a function of the [Fe/H] of member objects (Figure 6), their [Fe/H] distribution function (Figure 7), their other orbital properties (Figure 8), and the masses of their progenitor galaxies.

Below, we first discuss some of the key results of our study, then consider how our analysis advances our knowledge about the formation of the Milky Way halo, and then finally explore what information this study provides to galaxy formation models.

## 7.1. Key Results

First, we report the detection of a new merger, Pontus (see Section 4.6). This is a retrograde merger (i.e., *J*<sub>φ</sub> > 0) that

possesses relatively low energy (see Figure 5). Its low energy could indicate an early accretion, but this is hard to confirm at this stage. We note that some of the globular clusters in Pontus are quite metal poor and old (Figure 9). While Pontus is detected by our analysis, we find another group through a visual inspection of the (*J*, *E*) data, and this represents a new candidate merger (Section 5).

We find that the LMS-1/Wukong merger (Yuan et al. 2020a; Naidu et al. 2020) represents the most metal-poor merger of the Milky Way, with its MDF having a minimum of [Fe/H] ≈ −3.4 dex (Figure 7, Section 4.5). The three most metal-poor streams—C-19 ([Fe/H] = −3.38 ± 0.06 dex, Martin et al. 2022a), Sylgr ([Fe/H] = −2.92 ± 0.06, Roederer & Gnedin 2019), and Phoenix ([Fe/H] = −2.70 ± 0.06, Wan et al. 2020)—belong to this merger. Thus, the progenitor of LMS-1/Wukong was probably an ancient and extremely metal-poor protogalaxy that formed in the early universe. Furthermore, most of the member objects of LMS-1/Wukong possess low mass. The exceptions are the globular clusters NGC 5024 (M ~ 4.5 × 10<sup>5</sup>M<sub>☉</sub>) and NGC 5272/M3 (M ~ 4 × 10<sup>5</sup>M<sub>☉</sub>, Baumgardt et al. 2019), but we note that both of them are Type I clusters,<fnref n="23" /> and they are also metal poor (Harris 2010; Boberg et al. 2016). Moreover, NGC 5024 has also been confirmed to possess a stellar population similar to that of the LMS-1 stream (Malhan et al. 2021). These results, at some level, lend credence to the possibility that massive globular clusters can form inside metal-poor protogalaxies. In the future, it will be interesting to combine the chemical properties of the member clusters, and also streams, of the LMS-1/Wukong group (e.g., their [Fe/H], alpha abundances, etc.) to explore their formation inside the progenitor LMS-1/Wukong galaxy. In a broader context, this will also be useful to understand the formation of the most-metal-deficient

<footnote n="22">²² The mergers that we could not detect are Helmi (Helmi et al. 1999), Thamnos 1 and 2 (Koppelman et al. 2019a), and Kraken/Koala (Kruijssen et al. 2020; Forbes 2020).</footnote>

<footnote n="23">²³ Type I clusters correspond to those that exhibit a simpler stellar population and where first and second generations of stars are not easily distinguishable (likely because they host only one generation of stars, Milone et al. 2017).</footnote>
<figure><img src="image_20.png" /><figcaption>**Figure 11.** Comparing the color–magnitude distribution of the objects belonging to the candidate merger (see Section 5). We use the Gaia EDR3 photometry that is corrected for extinction using the Schlafly et al. (2014) dust-extinction map. The absolute magnitude ($M_{G0}$) is obtained by correcting the dereddened magnitudes for the distance of the objects. The distances of globular clusters are taken from Vasiliev & Baumgardt (2021), those of Tucana III are taken from McConnachie & Venn (2020) and for the stream from Table 1. The quoted [Fe/H] values correspond to the spectroscopic measurements of these objects (see text).</figcaption></figure>

globular clusters inside those ancient protogalaxies that formed
in the early universe.

Furthermore, in the context of the hierarchical buildup of the Milky Way halo, the distribution function of these accreted mergers can be qualitatively constrained by examining Figure 5 (and also Figure 10). For instance, there is an even distribution of mergers along the prograde and retrograde directions. A second interesting feature that we note in Figure 5(a) is that all the prograde mergers possess higher $J_z$ values and all the retrograde mergers possess lower $J_z$ values. Moreover, three of these mergers have their orbits inclined at $\phi \sim 90^\circ$ to the Galactic plane, i.e., they possess polar orbits structures in the Milky Way halo. These mergers include Sagittarius (Section 4.1), Cetus (Section 4.2), and LMS-1/Wukong (see Section 4.5). The presence of these “polar” mergers could be indicating that the polar orbits in the Milky Way are indeed stable (see Penarrubia et al. 2002), and therefore, those mergers that would have accreted along polar orbits would remain polar (although their orbits may still precess). We also note a lack of mergers that orbit close to the Galactic plane. A possible reason could be that objects accreted along such orbits are quickly phase-mixed, and this renders their detection rather difficult with present methods.

Our analysis also allows us to (roughly) constrain the accreted stellar mass of the Milky Way halo. This parameter depends critically on the most massive mergers, such as those we detect here. To compute this parameter we add up the stellar mass values of the independent mergers, which yields a Galactic stellar halo mass of $M_{*,MW \text{ halo}} \sim 10^9 M_{\odot}$. This value is

smaller than that obtained by some recent studies (e.g.,
~1.5 × 10⁹ M☉ from Deason et al. 2019), but we caution that
our value is based only on the member globular clusters in each
merger and not the member stellar streams, and the inclusion of
streams can change this result.

We also find various physical associations between stellar streams and other objects (i.e., globular clusters and satellite galaxies; see Section 6). Most of these associations are based on the similarities in their (*J*, *E*) and [Fe/H] values; however, a few of them are proposed based only on (*J*, *E*) comparisons. To confirm these latter connections, future [Fe/H] and stellar population comparisons are necessary.

## 7.2. What about the Remaining Halo Objects That are Not Identified as Members of Any Detected Group?

A large fraction of objects in our sample are not identified as members of any detected groups—be it the merger groups or the in situ groups (i.e., Galactic disk and bulge). Quantitatively, this is 58% of the sample of 257 objects at a $2\sigma$ threshold. We consider below some possible interpretations of the low-significance groupings.

With regard to those objects that possess a low $P_{\text{Group}}$ value and high energy (with, say, $E \gtrsim -1.6 \times 10^5 \text{ km}^2 \text{ s}^{-2}$), there are three possible scenarios: (1) Some of these objects were accreted inside the same mergers that we detected here, but our method could not confidently associate them with these mergers. This could owe either to the inefficiency of our method or large phase-space uncertainties of these objects, or both. (2) Some of these objects were accreted inside other mergers that we could not detect here (e.g., the Helmi structure; Helmi et al. 1999). A possible reason why we did not detect these mergers could be that these are highly phase-mixed substructures and their member objects are largely dispersed in (J, E) space. (3) A majority of these objects may have accreted inside a multitude of low-mass mergers, which perhaps had masses in the range $M_{\text{halo}} \lesssim 10^{7-9} M_{\odot}$. Because such low-mass galaxies typically harbor only ~1–2 globular clusters (e.g., Forbes et al. 2018a), their detection is not possible using our method.

To explore these possible scenarios, we consider Figure 12 that shows all those objects with $P_{\text{Group}} \ge 0.20$. This time we find 78 objects (30% of our total sample) representing potential merger objects. Lowering this threshold also reveals the candidate merger that we described in Section 5. But it is clearly difficult to interpret the reality of these low-significance objects.

With regard to those objects that possess a low $P_{\text{Group}}$ value and low energy (say, with $E \lesssim -1.6 \times 10^5 \text{ km}^2 \text{ s}^{-2}$), there are two possible scenarios: (1) All of these objects indeed represent the in situ population (i.e., disk objects or bulge objects); however, we could only identify those objects that are very tightly clumped in (J, E) space. (2) Some of these objects were probably accreted into the Milky Way but we could not identify them as members of high-significance groups. This may include the Kraken/Koala group (Kruijssen et al. 2020; Forbes 2020), as shown in Figure 17.

In particular, most of the satellite galaxies did not form part of any group according to our analysis. This could be because these satellites have only recently arrived in the Milky Way halo (as indicated by their high *E* values in Figure 2) and are themselves are undergoing merging (see Hammer et al. 2021).
<figure><img src="image_21.png" /><figcaption>**Figure 12.** Same as Figure 5, but lowering the probability threshold value $P_{Threshold}$ to 0.2 (see Section 7.2).</figcaption></figure>

Naturally, each satellite independently represents a (minor) merger event.

### 7.3. Caveats and Limitations

Our method cannot detect those mergers that initially only possessed a population of stars and no globular clusters or streams (e.g., low-mass galaxies with $M_{halo} \lesssim 10^{7-9} M_{\odot}$). Perhaps this is the reason that we could not detect the Thamnos substructure (Koppelman et al. 2019a) and/or the Helmi substructure (Helmi et al. 1999)<fnref n="24" />. Such mergers can only be identified using the phase-space information of the stars (e.g., Naidu et al. 2020). Hence, our study is capable of only putting a lower limit on the number of Milky Way mergers (in this case, $N=6$).

To detect mergers using ENLINK, we had to supplement the *J* information with the (redundant) *E* information. The additional information of *E* acts as a “weight,” thus enhancing the contrast between different detected mergers. As we noted above, this is because the relative uncertainties on *E* are smaller than those on *J*, and thus ENLINK uses this additional (precise) *E* information to detect groups. Future Gaia data sets will deliver much more precise phase-space measurements of the halo objects, and this in turn will reduce the uncertainties on the derived *J* quantities and it may eventually become possible to detect mergers by using only the *J* information.

Our strategy of detecting mergers relies on the derived (*J*, *E*) values of the halo objects, and these parameters depend on the choice of the Galactic potential. Admittedly, we tried only one potential model (McMillan 2017), but this is a realistic, static and axisymmetric model. On the other hand, there is now mounting evidence that the potential of the Milky Way halo is time dependent and asymmetric. This is due to the ongoing accretion of the massive LMC system that has perturbed the Galactic halo (especially in the outer halo regions $r_{gal} \gtrsim 30$ kpc;

e.g., Garavito-Camargo et al. 2021; Petersen & Penarrubia 2021). Studies indicate that this perturbation is significant enough to alter the orbits of halo objects. In future, it will be important to also try time-dependent potential models to integrate the orbits of the halo objects (e.g., Correa Magnus & Vasiliev 2021), compute *J* values along these orbital paths, and then detect the mergers.

Moreover, we also ignored the effect of the bar together with all the other (almost) unknowable time-dependent asymmetries that the galaxy has possessed over its long history, as well as the effect of dynamical friction. The effect of the bar should not have a significant effect on the detected mergers, first of all because its quadrupole potential drops off quickly with distance (as $r^{-3}$), and it will therefore solely affect (weakly) groups that have small pericenters, like Gaia–Sausage/Enceladus and Pontus. Resonances with the bar should not be an issue for the halo orbits under consideration, and because the streams have a relatively narrow range of energies, while the stream path itself will not be the same in the barred potential as in an axisymmetric one, the stream will not be disrupted. In any case, objects that remained close to each other in (*J*, *E*) space, as computed within the axisymmetric background potential, can certainly still be associated with each other. A likely much more important issue that we also ignored is dynamical friction. For instance, Villalobos & Helmi (2008) simulated the accretion of a merger with a mass of $10^{11} M_{\odot}$ and showed that the mass of the merger, which remains bound, decreases mostly during pericentric passages of the merger. For this reason, we cannot in principle formally exclude that Gaia–Sausage/Enceladus and Pontus represent different stages of disruption of one single massive accretion event. However, Pontus has a lower orbital energy than Gaia–Sausage/Enceladus, meaning it should then have been stripped later from the common progenitor galaxy, and while low-mass mergers are known to circularize under dynamical friction, high-mass ones tend to radialize (Vasiliev et al. 2021). Therefore, with a progenitor mass of $\sim 5 \times 10^{10} M_{\odot}$, Pontus is safely in the radialization regime, which would therefore make

<footnote n="24">²⁴ Although the Helmi substructure was previously (unambiguously) associated with five globular clusters by Massari et al. (2019), we find all of these clusters to be associated with one of the detected mergers.</footnote>
its lower eccentricity than Gaia–Sausage/Enceladus an argument against the scenario of these two groups being different stages of disruption of one single massive accretion event. Moreover, we recall that the two groups have a different AMR, which reinforces the hypothesis that these are in fact the remnants of distinct accretion events.

The dynamical atlas of the Milky Way mergers that we present here provides a global view of the galaxy formation in action. Thus, our study contributes to the initial steps of unraveling the full hierarchical buildup of our galaxy and also to understanding the origin of the globular clusters and stellar streams of the Milky Way halo. This endeavor of detecting mergers using the (J, E) quantities is only possible due to the amazingly rich phase-space information that the ESA/Gaia mission has provided, and this places us in a very exciting position to disentangle the merging events in the Milky Way halo. The community now looks forward to future Gaia data releases and upcoming spectroscopic surveys (e.g., WEAVE, Dalton et al. 2014, 4MOST, de Jong et al. 2019, SDSS-V), as the combination of these data sets will provide a gold mine of information in terms of phase space, metallicity, elemental abundances, and ages of stars. With such a wealth of information, we will be able to explore the “temporal” aspect of Galactic archeology (e.g., Feuillet et al. 2021) by building an understanding of the “chronological” merging history of the Milky Way.

We thank the referee for helpful comments and suggestions. K.M. acknowledges support from the Alexander von Humboldt Foundation at Max-Planck-Institut für Astronomie, Heidelberg. K.M. is also grateful to the IAU's Gruber Foundation Fellowship Programme for their financial support. R.I., B.F., Z.Y., and N.F.M. acknowledge funding from the Agence Nationale de la Recherche (ANR project ANR-18-CE31-0006, ANR-18-CE31-0017 and ANR-19-CE31-0017), from CNRS/INSU through the Programme National Galaxies et Cosmologie, and from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement No. 834148). M.B. acknowledges the financial support to this research by INAF, through the Mainstream grant 1.05.01.86.22 assigned to the project “Chemo-dynamics of globular clusters: the Gaia revolution” (P.I. E. Pancino). G.T. acknowledges support from the Agencia Estatal de Investigación (AEI) of the Ministerio de Ciencia e Innovación (MCINN) under grant FJC2018-037323-I and from the financial support through the grant (AEI/FEDER, UE) AYA2017-89076-P, as well as by the Ministerio de Ciencia,

Innovación y Universidades (MCIU), through the State Budget and by the Consejería de Economía, Industria, Comercio y Conocimiento of the Canary Islands Autonomous Community, through the Regional Budget.

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

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

SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofísica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.

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

## Appendix A

### Comparing the Orbital Phase and Eccentricity of the Halo Objects

As a passing analysis in Section 2.3, we compared the distribution of the orbital phase (f) and eccentricity (e) of globular clusters, stellar streams, and satellite galaxies. This comparison is shown in Figure 13.

<figure><img src="image_22.png" /><figcaption>**Figure 13.** Analyzing the orbital phase (f) and eccentricity (e) of globular clusters, stellar streams, and satellite galaxies. Panel (a) compares the f distribution of various objects, panel (b) compares the e distribution (where medians of different distributions are also quoted), and panel (c) shows f vs. e.</figcaption></figure>
Figure 13(a) compares the *f* distribution of different types of objects. We define the orbital phase of an object (in the radial direction) as $f = \frac{r_{gal} - r_{peri}}{r_{apo} - r_{peri}}$, where $r_{gal}$ is the present-day galactocentric distance of the object, and $r_{peri}$ and $r_{apo}$ represent its pericenter and apocenter distances, respectively. This definition of *f* comes from Fritz et al. (2018) and Li et al. (2021a). Per this definition, *f* = 0 implies that the object is at its pericenter and *f* = 1 implies that the object is at its apocenter. The *f* distribution in Figure 13(a) shows the piling-up of objects at both the pericenter and the apocenter. This effect is prevalent for globular clusters and stellar streams, but not so much for satellite galaxies. In particular, in the case of streams, we note that more objects are piled-up at the pericenter than at the apocenter. This effect points toward our inefficiency in detecting those streams that may currently be located at larger distances ($D_{\odot} \gtrsim 30$ kpc) close to their apocenters. This inefficiency, in part, could be due to Gaia's limiting magnitude at $G \sim 21$. We note that this inference on the *f* distribution of streams is different than that of Li et al. (2021a), who found that more streams are piled-up at the apocenter (likely due to a different selection strategy).

Figure 13(b) compares the *e* distribution of different types of objects. The medians of the *e* distribution for globular clusters, streams, and satellites are provided in panel (b). The median of the eccentricity distribution lies at *e* ≈ 0.5, which is true for all three types of objects. This implies that, in general, very radial orbits (with *e* ≈ 1) or very circular orbits (with *e* ≈ 0) are rare occurrences. This inference, in particular for streams, is consistent with that of Li et al. (2021a).

Figure 13(c) shows *e* as a function of *f*, and we infer no obvious trend from this plot.

## Appendix B
Determining the $P_{Threshold}$ Value

In Section 3, we required a suitable probability threshold value $P_{Threshold}$, such that all the halo objects with $P_{Group} \ge P_{Threshold}$ can be considered to belong to high-significance groups. To set this $P_{Threshold}$ value, we artificially construct a “randomized” version of the real (J, E) data as follows.

For a given halo object, we first randomly set its orbital pole and then recompute its (J, E) corresponding to this new orbit (taking the real (J, E) uncertainties into account). This is repeated for all the halo objects. This randomization erases any real correlations between the objects in (J, E) space. For a given halo object, this transformation from the real coordinate system to this new coordinate system is done as follows. From the 3D position vector and 3D velocity vector (x, v), we draw (isotropically) a random value (x', v'). After performing this transformation on all the objects, we compute the new (J, E) values. This new randomized (J, E) distribution for all the objects is shown in Figure 14. We note that, as a consequence of this transformation, the *E* parameter of a given object only slightly changes; however, the *J* parameters are significantly altered.

Next, we find groups in this randomized (J, E) distribution using exactly the same procedure as described in Section 3. Figure 15 shows the probabilities $P_{Group}$ for all these randomized objects in the (J, E) space. Initially, we were expecting that all the objects would receive very low $P_{Group}$ values (close to zero), because all the correlations must be erased and hence no groups would be found. While this is true for most of the objects in Figure 15, however, a few objects located at ($E, J_{\phi}$) ~ (-1.25 × 10⁵ km² s⁻², 0 kpc km s⁻¹) possess strikingly high $P_{Group}$ values. In other words, ENLINK has detected one significant group at this ($E, J_{\phi}$) location. Intriguingly, this ($E, J_{\phi}$) location coincides with the dynamical location of the Gaia–Sausage/Enceladus merger. We repeated the above (J, E) randomization procedure multiple times; however, a group was always detected at this ($E, L_z$) location. We thus tentatively suggest that this compromises the reality of those objects that have been linked with the Gaia–Sausage/Enceladus merger (both in our study and in previous studies). We defer the exploration of this problem to future work.

Figure 16 shows the PDF of the computed probabilities $P_{Group}$ of all randomized objects, and this serves as the background model for our main analysis. From this PDF, we use its 90th percentile to set the threshold probability at $P_{Threshold} = 0.3$ (rounded off from 0.295) for a 2σ detection.
Randomised population of globular clusters in (E, J) space

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

(a)

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

(b)

Randomised population of streams in (E, J) space

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

(a)

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

(b)

Randomised population of satellite galaxies in (E, J) space

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

(a)

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

(b)

<caption>**Figure 14.** Same as Figure 2, but for the artificially randomized (J, E) distribution drawn from the sample.</caption>
<figure><img src="image_29.png" /></figure>

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

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

<figure><img src="image_32.png" /><figcaption>**Figure 15.** Same as Figure 3, but for the artificially randomized (*J*, *E*) distribution.</figcaption></figure>
<figure><img src="image_33.png" /><figcaption>**Figure 16.** Same as Figure 4, but for the artificially randomized (*J*, *E*) distribution of objects. The vertical line shows the 90th percentile of this PDF, and we use the corresponding value to set the desired *P*<sub>Threshold</sub> value. The other quoted value is the median of the distribution, and the corresponding uncertainties reflect its 16th and 84th percentiles.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>**Figure 17.** (*J*, *E*) space showing groupings of objects as found by previous studies. Different groups are shown with different colors, the globular clusters are shown with the “star” markers, and stellar streams are shown with the “circle” markers.</figcaption></figure>

## Appendix C
### Groupings of Previous Studies

We construct Figure 17, which shows groupings of objects in the (*J*, *E*) space from previous studies. We do this so that the previous results can be compared with our result shown in Figure 5.

To construct Figure 17, we use our derived (*J*, *E*) values of globular clusters and streams, but use the associations that are described by previous studies. To this end, the globular cluster groupings are adopted from Table 1 of Massari et al. (2019) for the groups Sagittarius, Helmi, Gaia–Sausage/Enceladus, and Bulge; Table 1 of Myeong et al. (2019) for Sequoia; and Table 1 of Kruijssen et al. (2020) for Kraken. The stream groupings are

adopted from Table 1 of Bonaca et al. (2021). From these previous studies, we specifically highlight only those objects that have been unambiguously linked to just one group and avoid highlighting those objects that have been linked to multiple groups; the only exception is the Arjuna/Sequoia/I'toi group.

## ORCID iDs

* Khyati Malhan https://orcid.org/0000-0002-8318-433X
* Rodrigo A. Ibata https://orcid.org/0000-0002-3292-9709
* Sanjib Sharma https://orcid.org/0000-0002-0920-809X
* Benoit Famaey https://orcid.org/0000-0003-3180-9825
* Michele Bellazzini https://orcid.org/0000-0001-8200-810X
* Raymond G. Carlberg https://orcid.org/0000-0002-7667-0081
Richard D'Souza https://orcid.org/0000-0001-9269-8167

Zhen Yuan https://orcid.org/0000-0002-8129-5415

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

Guillaume F. Thomas https://orcid.org/0000-0002-2468-5521

## References

Balbinot, E., Yanny, B., Li, T. S., et al. 2016, ApJ, 820, 58

Battaglia, G., Taibi, S., Thomas, G. F., & Fritz, T. K. 2022, A&A, 657, A54

Baumgardt, H., Hilker, M., Sollima, A., & Bellini, A. 2019, MNRAS, 482, 5138

Bell, E. F., Zucker, D. B., Belokurov, V., et al. 2008, ApJ, 680, 295

Bellazzini, M., Ibata, R., Malhan, K., et al. 2020, A&A, 636, A107

Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, MNRAS, 478, 611

Belokurov, V., Evans, N. W., Irwin, M. J., Hewett, P. C., & Wilkinson, M. I. 2006, ApJL, 637, L29

Belokurov, V., Evans, N. W., Irwin, M. J., et al. 2007, ApJ, 658, 337

Bernard, E. J., Ferguson, A. M. N., Schlafly, E. F., et al. 2014, MNRAS, 443, L84

Boberg, O. M., Friel, E. D., & Vesperini, E. 2016, ApJ, 824, 5

Bonaca, A., Hogg, D. W., Price-Welan, A. M., & Conroy, C. 2019, ApJ, 880, 38

Bonaca, A., Naidu, R. P., Conroy, C., et al. 2021, ApJL, 909, L26

Bovy, J. 2015, ApJS, 216, 29

Bovy, J. 2020, arXiv:2012.02169

Bullock, J. S., & Johnston, K. V. 2005, ApJ, 635, 931

Caldwell, N., Bonaca, A., Price-Welan, A. M., Sesar, B., & Walker, M. G. 2020, AJ, 159, 287

Cardone, V. F., & Sereno, M. 2005, A&A, 438, 545

Carlberg, R. G. 2018, ApJ, 861, 69

Carlberg, R. G. 2020, ApJ, 889, 107

Carlin, J. L., & Sand, D. J. 2018, ApJ, 865, 7

Chang, J., Yuan, Z., Xue, X.-X., et al. 2020, ApJ, 905, 100

Chiba, M., & Beers, T. C. 2000, AJ, 119, 2843

Cohen, J. G., & Melendez, J. 2005, AJ, 129, 1607

Correa Magnus, L., & Vasiliev, E. 2021, arXiv:2110.00018

Dalton, G., Trager, S., Abrams, D. C., et al. 2014, Proc. SPIE, 9147, 91470L

de Jong, R. S., Agertz, O., Berbel, A. A., et al. 2019, Msngr, 175, 3

Deason, A. J., Belokurov, V., & Sanders, J. L. 2019, MNRAS, 490, 3426

Deason, A. J., Wetzel, A. R., Garrison-Kimmel, S., & Belokurov, V. 2015, MNRAS, 453, 3568

Dehnen, W., Odenkirchen, M., Grebel, E. K., & Rix, H.-W. 2004, AJ, 127, 2753

D'Souza, R., & Bell, E. F. 2018, MNRAS, 474, 5300

Erkal, D., Belokurov, V., Laporte, C. F. P., et al. 2019, MNRAS, 487, 2685

Fernández-Alvar, E., Carigi, L., Schuster, W. J., et al. 2018, ApJ, 852, 50

Feuillet, D. K., Sahlholdt, C. L., Feltzing, S., & Casagrande, L. 2021, MNRAS, 508, 1489

Forbes, D. A. 2020, MNRAS, 493, 847

Forbes, D. A., Read, J. I., Gieles, M., & Collins, M. L. M. 2018a, MNRAS, 481, 5592

Forbes, D. A., Bastian, N., Gieles, M., et al. 2018b, RSPSA, 474, 20170616

Frebel, A., & Norris, J. E. 2015, ARA&A, 53, 631

Freeman, K., & Bland-Hawthorn, J. 2002, ARA&A, 40, 487

Fritz, T. K., Battaglia, G., Pawlowski, M. S., et al. 2018, A&A, 619, A103

Gaia Collaboration, Brown, A. G., Vallenari, A., et al. 2021, A&A, 649, A1

Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1

Garavito-Camargo, N., Besla, G., Laporte, C. F. P., et al. 2021, ApJ, 919, 109

Gibbons, S. L. J., Belokurov, V., & Evans, N. W. 2017, MNRAS, 464, 794

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

Grillmair, C. J. 2006, ApJL, 645, L37

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

Grillmair, C. J. 2019, ApJ, 884, 174

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

Hammer, F., Wang, J., Pawloski, M., et al. 2021, ApJ, 922, 93

Hansen, T. T., Riley, A. H., Strigari, L. E., et al. 2020, ApJ, 901, 23

Harris, W. E. 2010, arXiv:1012.3224

Helmi, A. 2020, ARA&A, 58, 205

Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Natur, 563, 85

Helmi, A., White, S. D. M., de Zeeuw, P. T., & Zhao, H. 1999, Natur, 402, 53

Horta, D., Schiavon, R. P., Mackereth, J. T., et al. 2020, MNRAS, 493, 3363

Hudson, M. J., Harris, G. L., & Harris, W. E. 2014, ApJL, 787, L5

Ibata, R., Bellazzini, M., Thomas, G., et al. 2020, ApJL, 891, L19

Ibata, R., Malhan, K., Martin, N., et al. 2021, ApJ, 914, 123

Ibata, R. A., Bellazzini, M., Malhan, K., Martin, N., & Bianchini, P. 2019a, NatAs, 3, 667

Ibata, R. A., Gilmore, G., & Irwin, M. J. 1994, Natur, 370, 194

Ibata, R. A., Malhan, K., & Martin, N. F. 2019b, ApJ, 872, 152

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

Ishigaki, M. N., Hwang, N., Chiba, M., & Aoki, W. 2016, ApJ, 823, 157

Jean-Baptiste, I., Di Matteo, P., Haywood, M., et al. 2017, A&A, 604, A106

Jensen, J., Thomas, G., McConnachie, A. W., et al. 2021, MNRAS, 507, 1923

Karukes, E. V., Benito, M., Iocco, F., Trotta, R., & Geringer-Sameth, A. 2020, JCAP, 05, 033

Kim, D., Jerjen, H., Mackey, D., Da Costa, G. S., & Milone, A. P. 2016, ApJ, 820, 119

Kirby, E. N., Cohen, J. G., Guhathakurta, P., et al. 2013, ApJ, 779, 102

Komiya, Y., Habe, A., Suda, T., & Fujimoto, M. Y. 2010, ApJ, 717, 542

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

Koposov, S. E., Belokurov, V., Li, T. S., et al. 2019, MNRAS, 485, 4726

Koppelman, H. H., Helmi, A., Massari, D., Price-Welan, A. M., & Starkenburg, T. K. 2019a, A&A, 631, L9

Koppelman, H. H., Helmi, A., Massari, D., Roelenga, S., & Bastian, U. 2019b, A&A, 625, A5

Kruijssen, J. M. D., Pfeffer, J. L., Reina-Campos, M., Crain, R. A., & Bastian, N. 2019, MNRAS, 486, 3180

Kruijssen, J. M. D., Pfeffer, J. L., Chevance, M., et al. 2020, MNRAS, 498, 2472

Lamb, M. P., Venn, K. A., Shetrone, M. D., Sakari, C. M., & Pritzl, B. J. 2015, MNRAS, 448, 42

Lane, J. M. D., Bovy, J., & Mackereth, J. T. 2022, MNRAS, 510, 5119

Li, T. S., Koposov, S. E., Zucker, D. B., et al. 2019, MNRAS, 490, 3508

Li, T. S., Ji, A. P., Pace, A. B., et al. 2021a, arXiv:2110.06950

Li, T. S., Koposov, S. E., Erkal, D., et al. 2021b, ApJ, 911, 149

Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2

Mackereth, J. T., & Bovy, J. 2020, MNRAS, 492, 3631

Majewski, S. R., Skrutskie, M. F., Weinberg, M. D., & Ostheimer, J. C. 2003, ApJ, 599, 1082

Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, AJ, 154, 94

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

Malhan, K., & Ibata, R. A. 2019, MNRAS, 486, 2995

Malhan, K., Ibata, R. A., Carlberg, R. G., et al. 2019a, ApJL, 886, L7

Malhan, K., Ibata, R. A., Carlberg, R. G., Valluri, M., & Freese, K. 2019b, ApJ, 881, 106

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

Malhan, K., Ibata, R. A., & Martin, N. F. 2020, arXiv:2012.05271

Malhan, K., Yuan, Z., Ibata, R. A., et al. 2021, ApJ, 920, 51

Martin, N. F., Ibata, R. A., Bellazzini, M., et al. 2004, MNRAS, 348, 12

Martin, N. F., Venn, K. A., Aguado, D. S., et al. 2022a, Natur, 601, 45

Martin, N. F., Ibata, R. A., Starkenburg, E., et al. 2022b, arXiv:2201.01310

Massari, D., Koppelman, H. H., & Helmi, A. 2019, A&A, 630, L4

Matsuno, T., Aoki, W., & Suda, T. 2019, ApJL, 874, L35

McConnachie, A. W., & Venn, K. A. 2020, RNAAS, 4, 229

McMillan, P. J. 2017, MNRAS, 465, 76

Milone, A. P., Piotto, G., Renzini, A., et al. 2017, MNRAS, 464, 3636

Myeong, G. C., Evans, N. W., Belokurov, V., Sanders, J. L., & Koposov, S. E. 2018, ApJL, 863, L28

Myeong, G. C., Vasiliev, E., Iorio, G., Evans, N. W., & Belokurov, V. 2019, MNRAS, 488, 1235

Naidu, R. P., Conroy, C., Bonaca, A., et al. 2020, ApJ, 901, 48

Naidu, R. P., Conroy, C., Bonaca, A., et al. 2021, ApJ, 923, 92

Newberg, H. J., Yanny, B., & Willett, B. A. 2009, ApJL, 700, L61

Niederer-Ostholt, M., Belokurov, V., & Evans, N. W. 2012, MNRAS, 422, 207

Nissen, P. E., & Schuster, W. J. 2010, A&A, 511, L10

Nitschai, M. S., Eilers, A.-C., Neumayer, N., Cappellari, M., & Rix, H.-W. 2021, ApJ, 916, L12

Obreja, A., Buck, T., & Macciò, A. V. 2022, A&A, 657, A15

Odenkirchen, M., Grebel, E. K., Rockosi, C. M., et al. 2001, ApJL, 548, L165

Palau, C. G., & Miralda-Escude, J. 2019, MNRAS, 488, 1535

Palau, C. G., & Miralda-Escude, J. 2021, MNRAS, 504, 2727

Panithanpaisal, N., Sanderson, R. E., Wetzel, A., et al. 2021, ApJ, 920, 10

Pawlowski, M. S., Oria, P.-A., Taibi, S., Famaey, B., & Ibata, R. 2021, arXiv:2111.05358

Penarrubia, J., Kroupa, P., & Boily, C. M. 2002, MNRAS, 333, 779

Petersen, M. S., & Penarrubia, J. 2021, NatAs, 5, 251

Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077
Read, J. I., Iorio, G., Agertz, O., & Fraternali, F. 2017, MNRAS, 467, 2019

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

Robertson, B., Bullock, J. S., Font, A. S., Johnston, K. V., & Hernquist, L. 2005, ApJ, 632, 872

Roederer, I. U., & Gnedin, O. Y. 2019, ApJ, 883, 84

Sanders, J. L., & Binney, J. 2013, MNRAS, 433, 1813

Schlafly, E. F., Green, G., Finkbeiner, D. P., et al. 2014, ApJ, 789, 15

Schonrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829

Sestito, F., Buck, T., Starkenburg, E., et al. 2021, MNRAS, 500, 3750

Sharma, S., & Johnston, K. V. 2009, ApJ, 703, 1061

Sharma, S., Johnston, K. V., Majewski, S. R., et al. 2010, ApJ, 722, 750

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

Simon, J. D., Li, T. S., Drlica-Wagner, A., et al. 2017, ApJ, 838, 11

Simpson, J. D. 2018, MNRAS, 477, 4565

Sivia, D. 1996, Data Analysis: A Bayesian Tutorial (Clarendon Press: Oxford Science Publications)

Sohn, S. T., Watkins, L. L., Fardal, M. A., et al. 2018, ApJ, 862, 52

Sollima, A. 2020, MNRAS, 495, 2222

Thomas, G. F., Ibata, R., Famaey, B., Martin, N. F., & Lewis, G. F. 2016, MNRAS, 460, 2711

Thomas, G. F., Jensen, J., McConnachie, A., et al. 2020, ApJ, 902, 89

Vasiliev, E., & Baumgardt, H. 2021, MNRAS, 505, 5978

Vasiliev, E., & Belokurov, V. 2020, MNRAS, 497, 4162

Vasiliev, E., Belokurov, V., & Evans, W. 2021, arXiv:2108.00010

Villalobos, Á., & Helmi, A. 2008, MNRAS, 391, 1806

Wan, Z., Lewis, G. F., Li, T. S., et al. 2020, Natur, 583, 768

Willman, B., Blanton, M. R., West, A. A., et al. 2005, AJ, 129, 2692

Wu, Y., Valluri, M., Panithanpaisal, N., et al. 2022, MNRAS, 509, 5882

Yam, W., Carlin, J. L., Newberg, H. J., et al. 2013, ApJ, 776, 133

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

Yuan, Z., Chang, J., Beers, T. C., & Huang, Y. 2020a, ApJL, 898, L37

Yuan, Z., Smith, M. C., Xue, X.-X., et al. 2019, ApJ, 881, 164

Yuan, Z., Myeong, G. C., Beers, T. C., et al. 2020b, ApJ, 891, 39

Yuan, Z., Malhan, K., Sestito, F., et al. 2021, arXiv:2112.05775

Zhao, G., Zhao, Y., Chu, Y., Jing, Y., & Deng, L. 2012, RAA, 12, 723

Zinn, R. 1985, ApJ, 293, 424