The Global Dynamical Atlas of the Milky Way Mergers: Constraints from Gaia EDR3– based Orbits of Globular Clusters, Stellar Streams, and Satellite Galaxies Khyati Malhan1,9 , Rodrigo A. Ibata2 , Sanjib Sharma3,4 , Benoit Famaey2 , Michele Bellazzini5 , Raymond G. Carlberg6 , Richard D’Souza7 , Zhen Yuan2 , Nicolas F. Martin1,2 , and Guillaume F. Thomas8 1 Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117, Heidelberg, Germany; kmalhan07@gmail.com 2 Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France 3 Sydney Institute for Astronomy, School of Physics, The University of Sydney, NSW 2006, Australia 4 ARC Centre of Excellence for All Sky Astrophysics in Three Dimensions (ASTRO-3D), Australia 5 INAF—Osservatorio di Astrofisica e Scienza dello Spazio, via Gobetti 93/3, I-40129 Bologna, Italy 6 Department of Astronomy & Astrophysics, University of Toronto, Toronto, ON M5S 3H4, Canada 7 Vatican Observatory, Specola Vaticana, V-00120, Vatican City State, Italy 8 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 composi- tion (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 chemody- namical 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 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 https://doi.org/10.3847/1538-4357/ac4d2a © 2022. The Author(s). Published by the American Astronomical Society. 9 Humboldt Fellow and IAU Gruber Fellow. 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. 1 objects represent the most ancient and metal-poor structures of our galaxy (e.g., Harris 2010; Kirby et al. 2013; Helmi 2020).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.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,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 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 (Mhalo  109−11Me; e.g., Robertson et al. 2005; D’Souza & Bell 2018) are typically much smaller than that of our galaxy (Mhalo ∼ 10 12Me; 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.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 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. 11 Conceptually, J represents the amplitude of an object’s orbit along different directions. For instance, in cylindrical coordinates, J ≡ (JR, Jf, Jz), where Jf represents the z component of angular momentum (≡Lz), and JR and Jz describe the extent of oscillations in cylindrical radius and z directions, respectively. 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). 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. 2 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 Malhan et al. 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 (De) or parallax (ϖ), 2D proper motion ( cos , m m d m º a a d * ), 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 measure- ments from the Vasiliev & Baumgardt (2021) catalog. This catalog provides, for 170 globular clusters, their phase-space measurements, and we use the observed heliocentric coordi- nates (i.e., α, δ, De, ma*, μδ, vlos) along with the associated uncertainties. Vasiliev & Baumgardt (2021) derive the 4D astrometric measurement (α, δ, ma* , μδ) of globular clusters using the Gaia EDR3 data set, while the parameters De and vlos are based on a combination of Gaia EDR3 and other surveys. For satellite galaxies, we obtain their phase-space measure- ments 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 De < 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−1 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 (ma*, μδ)= (−2.67± 0.45,−1.40± 0.40)mas yr −1. 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 (α, δ, ϖ, ma*, μδ), 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 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 mea- surements 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” (Grill- mair 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 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), “M5” (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 measure- ments 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. 3 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 Malhan et al. 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 rgal < 20 kpc is 2.2 × 10 11Me, rgal < 50 kpc is 4.9× 1011Me, and rgal < 100 kpc is 8.1 × 10 11Me. 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… 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≡ (JR, Jf, Jz) system, where Jf corresponds to the z component of angular momentum (i.e., Jf ≡ Lz) and negative Jf represents prograde motion (i.e., rotational motion in the direction of the Galactic disk). Similarly, components JR and Jz 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 Jf/Jtotal 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 (De), line-of-sight velocities (vlos), proper motion in the Galactic longitude direction ( ℓm*), and proper motion in the Galactic latitude direction direction (μb), respectively. In panel (d), we only show streams with their names and do not plot other objects to avoid crowding. 4 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 Malhan et al. 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 Jf/Jtot and the vertical axis is (Jz – JR)/Jtot with Jtot = JR + Jz + |Jf|. 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 (Jf ≡ Lz) vs. E, and the points are colored by the orthogonal component of their angular momenta (L L L x y 2 2 = + ^ ). 5 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 Malhan et al. versus (Jz − JR)/Jtotal, where Jtotal = JR + Jz + |Jf|, and (2) the Jf 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 L L x y 2 2 = + ^ for representation. Note that even though L⊥ 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., rapo, rperi, 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.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 (( ) ) ( )  N D ln ln 2 ln ln , 1 i v 5 2 sky los p s s s s s = - + - v m m a d where ( ) ( ) ( ) ( ) ( ) ( ) N e D R R R R R R v v 1 , , , , , , . 2 j R j j v 1 5 2 1 5 2 1 2 sky 2 sky 2 2 2 d o 2 2 3 2 d o 2 2 4 2 d o 2 2 5 2 los d los o 2 2 j 2 los   q s v v s m m s m m s s = - = = = - = - = - = - v a a m d d m = - = a d Here, θsky is the on-sky angular difference between the orbit and the data point, , , d d d v m m a d and vlos 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 , , , , v sky los s s s s s v m m a d 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 vlos = 0 km s −1 with a 104 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 measure- ments 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 ω 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 vlos space (∼10–15 km s −1; 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 (α, δ) and (ma*, μδ) spaces, it was difficult to visualize this linearity in vlos space. This was primarily because these streams lack enough spectroscopic measurements so that a clear stream signal can be visible in vlos 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 α, δ, ma*, μδ, and vlos 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. 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. 6 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 Malhan et al. 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 Jf < 0), 18% move along polar orbits (i.e., their orbital planes are inclined almost perpendi- cularly to the Galactic disk plane, with f … 75°), ∼12% have their orbits nearly confined to the Galactic plane (i.e., their f „ 20°), 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 Hrìd. 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 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) Stream No. of Gaia No. of R.A. Decl. De ma* μδ vlos Sources vlos sources (deg) (deg) (kpc) (mas yr−1) (mas yr−1) (km s−1) Gjoll 102 35 82.1 13.95 0.36 0.35 - - + 3.26 0.03 0.03 - + 23.58 0.09 0.08 - + 23.7 0.05 0.06 - - + 78.73 1.84 2.36 - + Leiptr 237 67 89.11 28.37 0.27 0.2 - - + 7.39 0.07 0.07 - + 10.59 0.04 0.03 - + 9.9 0.04 0.03 - - + 194.22 1.86 2.23 - + Hrid 233 24 280.51 33.3 0.6 0.75 - + 2.75 0.07 0.1 - + 5.88 0.08 0.11 - - + 20.08 0.19 0.21 - + 238.77 5.52 3.3 - - + Pal5 48 29 229.65 0.26 0.13 0.1 - + 20.16 0.33 0.24 - + 2.75 0.02 0.03 - - + 2.68 0.02 0.02 - - + 57.03 1.04 1.08 - - + Gaia-1 106 8 190.96 9.16 0.1 0.15 - - + 5.57 0.1 0.16 - + 14.39 0.04 0.04 - - + 19.72 0.04 0.03 - - + 214.91 2.16 3.5 - + Ylgr 699 32 173.82 22.31 0.3 0.22 - - + 9.72 0.14 0.16 - + 0.44 0.03 0.04 - - + 7.65 0.04 0.05 - - + 317.86 3.05 2.83 - + Fjorm 182 28 251.89 65.38 0.22 0.24 - + 6.42 0.14 0.16 - + 3.92 0.08 0.07 - + 3.1 0.06 0.06 - + 25.37 2.19 1.89 - - + Kshir 55 16 205.88 67.25 0.17 0.13 - + 9.57 0.08 0.08 - + 7.67 0.04 0.04 - - + 3.92 0.05 0.04 - - + 249.88 2.92 2.62 - - + Gunnthra 61 8 284.22 73.49 0.14 0.23 - - + 2.83 0.13 0.12 - + 15.83 0.13 0.11 - - + 24.04 0.17 0.15 - - + 132.26 4.97 6.23 - + Slidr 181 29 160.05 10.22 0.41 0.43 - + 2.99 0.09 0.11 - + 24.6 0.08 0.08 - - + 6.65 0.06 0.06 - - + 87.98 3.17 3.44 - - + M92 84 9 259.89 43.08 0.2 0.2 - + 8.94 0.18 0.2 - + 5.15 0.05 0.05 - - + 0.63 0.04 0.06 - - + 140.66 7.53 6.28 - - + NGC 3201 388 4 152.46 46.32 0.08 0.11 - - + 4.99 0.02 0.01 - + 8.87 0.02 0.02 - + 2.22 0.02 0.02 - - + 489.63 3.82 3.36 - + Atlas 46 10 25.04 29.81 0.1 0.1 - - + 19.93 0.75 0.76 - + 0.04 0.02 0.02 - + 0.89 0.02 0.02 - - + 85.65 1.58 1.48 - - + C-7 120 10 287.15 50.17 0.14 0.16 - - + 6.77 0.21 0.28 - + 13.79 0.06 0.07 - - + 12.38 0.07 0.06 - - + 55.05 2.51 1.5 - + Palca 24 24 36.57 36.15 0.31 0.33 - - + 12.31 1.44 1.68 - + 0.9 0.02 0.02 - + 0.23 0.04 0.04 - - + 106.32 2.54 2.62 - + Sylgr 165 19 179.68 2.44 0.4 0.27 - - + 3.77 0.11 0.07 - + 13.98 0.14 0.12 - - + 12.9 0.1 0.09 - - + 184.8 8.15 15.48 - - + Gaia-9 286 15 233.27 60.42 0.11 0.04 - + 4.68 0.09 0.08 - + 12.49 0.12 0.11 - - + 6.37 0.08 0.14 - + 359.86 4.11 4.65 - - + Gaia-10 90 9 161.47 15.17 0.14 0.14 - + 13.32 0.28 0.34 - + 4.14 0.05 0.05 - - + 3.15 0.04 0.04 - - + 289.64 3.32 2.75 - + Gaia-12 38 1 41.05 16.45 0.13 0.13 - + 15.71 1.03 1.29 - + 5.84 0.05 0.05 - + 5.66 0.06 0.07 - - + 303.83 15.07 22.55 - - + Indus 454 45 340.12 60.58 0.1 0.1 - - + 14.96 0.16 0.19 - + 3.59 0.03 0.03 - + 4.89 0.03 0.02 - - + 49.15 3.68 2.45 - - + Jhelum 972 246 351.95 51.74 0.08 0.08 - - + 11.39 0.15 0.13 - + 7.23 0.04 0.04 - + 4.37 0.03 0.04 - - + 1.29 3.11 2.6 - - + Phoenix 35 19 23.96 50.01 0.24 0.24 - - + 16.8 0.36 0.33 - + 2.72 0.03 0.03 - + 0.07 0.03 0.03 - - + 45.92 1.58 1.63 - + NGC5466 62 4 214.41 26.84 0.11 0.12 - + 14.09 0.25 0.27 - + 5.64 0.03 0.03 - - + 0.72 0.02 0.03 - - + 95.04 5.91 7.4 - + M5 139 5 206.96 13.5 0.14 0.15 - + 7.44 0.11 0.12 - + 3.5 0.04 0.03 - + 8.76 0.04 0.04 - - + 42.97 3.83 3.33 - - + C-20 34 9 359.81 8.63 0.16 0.16 - + 18.11 1.39 1.45 - + 0.58 0.03 0.03 - - + 1.44 0.02 0.02 - + 116.87 1.44 1.46 - - + C-19 34 8 355.28 28.82 1.17 0.63 - + 18.04 0.53 0.55 - + 1.25 0.03 0.03 - + 2.74 0.05 0.03 - - + 193.48 2.52 2.61 - - + Elqui 4 4 19.77 42.36 0.29 0.3 - - + 51.41 7.04 4.64 - + 0.33 0.03 0.02 - + 0.49 0.02 0.02 - - + 15.86 20.38 8.82 - + AliqaUma 5 5 34.08 33.97 0.34 0.31 - - + 21.48 1.2 2.32 - + 0.24 0.03 0.02 - + 0.79 0.03 0.03 - - + 42.33 2.23 2.29 - - + Phlegethon 365 41 319.89 32.07 0.37 0.43 - - + 3.29 0.05 0.05 - + 3.97 0.09 0.09 - - + 37.66 0.09 0.08 - - + 15.9 6.12 4.97 - + GD-1 811 216 160.02 45.9 0.19 0.25 - + 8.06 0.07 0.07 - + 6.75 0.03 0.04 - - + 10.88 0.05 0.04 - - + 101.83 2.47 2.05 - - + 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 (vlos), R.A. (which acts as the anchor point in our orbit-fitting procedure), decl., heliocentric distance (De), proper motions ( , m m a d * ), and vlos. 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). 7 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 Malhan et al. 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 De  30 kpc). This inefficiency, in part, is also Table 2 Actions, Energies, Orbital Parameters, and Metallicities of the Stellar Streams Stream (JR, Jf, Jz) Energy rperi rapo zmax ecc. [Fe/H] (kpc km s−1) (×102 km2 s−2) (kpc) (kpc) (kpc) (dex) LMS-1 ( ) 255 , 627 , 2514 149 239 232 183 263 383 - - + - + - + 1227 39 65 - - + 10.8 1.8 2.5 - + 20.6 1.9 3.7 - + 20.2 2.0 3.6 - + 0.32 0.12 0.11 - + −2.1 ± 0.4 Gjoll ( ) 783 , 2782 , 274 65 73 60 60 6 7 - + - + - + 1152 18 19 - - + 8.5 0.1 0.1 - + 27.4 1.2 1.4 - + 10.8 0.5 0.5 - + 0.52 0.01 0.01 - + −1.78 Leiptr ( ) 1455 , 4128 , 378 119 133 73 77 10 11 - + - + - + 933 18 19 - - + 12.3 0.1 0.1 - + 45.1 2.1 2.3 - + 17.6 0.9 0.9 - + 0.57 0.01 0.01 - + L Hrid ( ) 1642 , 78 , 83 79 110 47 54 5 6 - + - + - + 1319 22 30 - - + 1.1 0.0 0.0 - + 22.0 1.0 1.4 - + 7.0 0.5 0.9 - + 0.9 0.0 0.0 - + −1.1 Pal5 ( ) 282 , 744 , 1357 17 19 44 61 51 42 - - + - + - + 1385 15 11 - - + 6.9 0.4 0.3 - + 15.8 0.3 0.2 - + 14.7 0.3 0.2 - + 0.39 0.02 0.02 - + −1.35 ± 0.06 Orphan ( ) 959 , 3885 , 1199 271 978 1017 405 213 484 - - + - + - + 949 64 175 - - + 15.6 2.1 3.8 - + 41.2 6.2 23.6 - + 26.4 4.8 16.9 - + 0.48 0.06 0.09 - + −1.85 ± 0.53 Gaia-1 ( ) 3638 , 2678 , 997 632 1307 57 89 58 93 - + - + - + 794 55 90 - - + 8.2 0.0 0.1 - + 67.6 8.9 18.3 - + 45.7 6.5 13.6 - + 0.78 0.03 0.04 - + −1.36 Fimbulthul ( ) 202 , 427 , 95 78 109 588 244 44 197 - + - + - + 1847 16 73 - - + 2.4 0.7 0.8 - + 7.2 0.3 0.4 - + 2.4 0.7 3.5 - + 0.51 0.13 0.12 - + −1.36 to −1.8 Ylgr ( ) 205 , 2766 , 556 35 40 66 68 25 30 - + - + - + 1219 19 20 - - + 11.5 0.1 0.1 - + 20.7 1.1 1.2 - + 11.2 0.7 0.8 - + 0.29 0.02 0.02 - + −1.87 Fjorm ( ) 831 , 2332 , 877 62 60 24 24 40 43 - - + - + - + 1123 15 15 - - + 9.1 0.1 0.1 - + 29.1 1.1 1.1 - + 19.7 1.0 1.0 - + 0.52 0.01 0.01 - + −2.2 Kshir ( ) 18 , 2755 , 491 4 6 53 60 13 14 - + - + - + 1268 11 12 - - + 13.4 0.2 0.2 - + 16.0 0.5 0.6 - + 8.2 0.3 0.3 - + 0.09 0.01 0.01 - + −1.78 Cetus ( ) 815 , 2416 , 2287 317 513 1064 841 954 1282 - - + - + - + 1000 64 124 - - + 14.7 4.5 7.2 - + 35.9 3.7 9.9 - + 30.2 4.9 10.9 - + 0.45 0.1 0.14 - + −2.0 Svol ( ) 97 , 1501 , 224 32 94 248 384 54 107 - - + - + - + 1566 61 89 - - + 5.9 0.6 0.6 - + 10.0 1.0 2.8 - + 5.0 0.8 0.9 - + 0.28 0.05 0.07 - + −1.98 ± 0.10 Gunnthra ( ) 69 , 852 , 218 7 14 77 67 34 30 - + - + - + 1765 28 31 - - + 4.2 0.4 0.3 - + 7.2 0.2 0.3 - + 3.8 0.4 0.4 - + 0.27 0.02 0.03 - + L Slidr ( ) 1076 , 1358 , 1831 149 217 23 23 102 126 - - + - + - + 1086 32 41 - - + 8.7 0.1 0.1 - + 32.3 2.5 3.5 - + 29.1 2.4 3.4 - + 0.58 0.02 0.03 - + −1.8 M92 ( ) 361 , 181 , 544 9 9 41 39 70 83 - + - + - + 1639 10 12 - - + 3.0 0.1 0.2 - + 10.7 0.2 0.2 - + 9.9 0.6 0.5 - + 0.56 0.01 0.01 - + −2.16 ± 0.05 NGC 6397 ( ) 75 , 586 , 222 5 5 29 15 14 23 - - + - + - + 1851 6 11 - - + 3.4 0.1 0.1 - + 6.4 0.1 0.1 - + 3.7 0.1 0.2 - + 0.3 0.01 0.01 - + L NGC 3201 ( ) 975 , 2860 , 296 48 48 33 32 5 5 - + - + - + 1110 10 10 - - + 8.5 0.0 0.0 - + 30.5 0.8 0.8 - + 12.3 0.3 0.3 - + 0.56 0.01 0.01 - + L Ophiuchus ( ) 507 , 160 , 1192 202 387 41 34 98 91 - - + - + - + 1490 84 130 - - + 3.9 0.4 0.3 - + 14.2 2.7 4.9 - + 14.1 2.6 4.9 - + 0.58 0.1 0.11 - + −1.80 ± 0.09 Atlas ( ) 757 , 1817 , 2093 35 39 33 34 86 100 - - + - + - + 1061 11 12 - - + 11.7 0.3 0.3 - + 32.4 0.9 1.0 - + 28.6 1.0 1.1 - + 0.47 0.01 0.01 - + −2.22 ± 0.03 C-7 ( ) 1059 , 706 , 728 215 397 25 17 69 107 - + - + - + 1319 65 100 - - + 3.5 0.0 0.0 - + 21.0 2.8 5.2 - + 18.1 2.8 5.1 - + 0.72 0.04 0.05 - + L C-3 ( ) 142 , 468 , 872 64 538 1185 1110 510 910 - + - + - + 1571 111 338 - - + 5.7 1.2 2.0 - + 10.0 1.4 13.2 - + 8.7 1.3 12.5 - + 0.35 0.1 0.18 - + L Palca ( ) 91 , 1830 , 1076 24 37 29 28 128 138 - - + - + - + 1300 28 30 - - + 10.8 0.3 0.4 - + 16.5 1.3 1.5 - + 12.7 1.4 1.5 - + 0.21 0.03 0.03 - + −2.02 ± 0.23 Sylgr ( ) 602 , 702 , 2220 202 141 24 28 153 94 - - + - + - + 1192 61 36 - - + 8.7 0.0 0.0 - + 24.6 3.7 2.4 - + 23.8 3.7 2.4 - + 0.48 0.06 0.03 - + −2.92 ± 0.06 Gaia-6 ( ) 125 , 907 , 557 71 51 229 342 204 288 - + - + - + 1593 70 129 - - + 6.0 1.8 1.4 - + 9.5 0.3 3.1 - + 6.9 0.9 3.6 - + 0.3 0.1 0.08 - + −1.16 Gaia-9 ( ) 393 , 1928 , 852 38 37 61 47 20 22 - + - + - + 1255 18 15 - - + 8.7 0.1 0.1 - + 20.8 0.9 0.8 - + 14.7 0.6 0.6 - + 0.41 0.01 0.01 - + −1.94 Gaia-10 ( ) 2189 , 287 , 1542 64 66 43 43 140 155 - + - + - + 1051 18 20 - - + 4.3 0.3 0.4 - + 37.7 1.4 1.6 - + 37.2 1.5 1.6 - + 0.8 0.01 0.01 - + −1.4 Gaia-12 ( ) 9834 , 7340 , 794 4378 5378 801 608 107 80 - + - + - + 433 142 98 - - + 18.5 1.2 0.9 - + 194.3 75.0 96.8 - + 83.0 32.5 42.9 - + 0.83 0.08 0.05 - + −2.6 Indus ( ) 99 , 1121 , 2211 19 25 36 35 49 61 - - + - + - + 1232 15 17 - - + 12.6 0.1 0.2 - + 18.9 0.9 1.0 - + 17.8 0.8 0.9 - + 0.2 0.02 0.02 - + −1.96 ± 0.41 Jhelum ( ) 594 , 356 , 2557 56 49 17 19 72 62 - - + - + - + 1193 22 17 - - + 8.7 0.2 0.2 - + 24.5 1.3 1.1 - + 24.3 1.3 1.1 - + 0.48 0.01 0.01 - + −1.83 ± 0.34 Phoenix ( ) 107 , 1563 , 1578 10 11 32 35 62 56 - - + - + - + 1259 12 10 - - + 11.7 0.5 0.4 - + 18.1 0.3 0.3 - + 15.6 0.3 0.3 - + 0.22 0.01 0.01 - + −2.70 ± 0.06 NGC5466 ( ) 1769 , 619 , 1373 114 144 35 38 80 93 - + - + - + 1098 24 29 - - + 4.8 0.2 0.3 - + 33.7 1.8 2.3 - + 31.8 1.7 2.2 - + 0.75 0.01 0.01 - + L M5 ( ) 1366 , 353 , 931 57 69 18 20 40 46 - - + - + - + 1246 14 17 - - + 3.4 0.1 0.1 - + 24.8 0.8 1.0 - + 23.6 0.9 1.1 - + 0.76 0.01 0.01 - + −1.34 ± 0.05 C-20 ( ) 1329 , 3042 , 3823 350 526 167 131 484 503 - - + - + - + 800 59 65 - - + 20.8 1.3 1.3 - + 58.5 8.8 12.0 - + 52.4 8.5 11.3 - + 0.47 0.04 0.05 - + −2.44 NGC7089 ( ) 800 , 638 , 359 270 713 555 414 123 99 - - + - + - + 1504 104 307 - - + 2.9 0.7 1.2 - + 14.7 3.0 12.7 - + 10.9 3.8 7.5 - + 0.71 0.06 0.06 - + L C-19 ( ) 383 , 210 , 2712 47 53 46 48 253 258 - - + - + - + 1232 21 21 - - + 9.3 1.0 1.0 - + 21.6 0.5 0.5 - + 21.6 0.6 0.5 - + 0.4 0.03 0.04 - + −3.38 ± 0.06 Elqui ( ) 2072 , 273 , 4324 602 543 191 166 352 321 - - + - + - + 868 35 27 - - + 12.1 2.0 1.8 - + 54.0 6.2 4.6 - + 53.9 6.3 4.6 - + 0.64 0.08 0.07 - + −2.22 ± 0.37 Chenab ( ) 2469 , 4062 , 3735 286 463 509 601 338 341 - - + - + - + 690 53 52 - - + 22.0 2.7 2.2 - + 81.0 10.4 12.8 - + 69.1 8.2 10.6 - + 0.58 0.01 0.01 - + −1.78 ± 0.34 AliqaUma ( ) 738 , 1838 , 2025 71 138 64 96 126 223 - - + - + - + 1067 18 30 - - + 11.6 0.3 0.4 - + 31.9 1.5 2.6 - + 27.9 1.6 3.0 - + 0.47 0.02 0.03 - + −2.30 ± 0.06 Phlegethon ( ) 815 , 1882 , 231 93 120 37 39 10 11 - + - + - + 1272 28 32 - - + 5.5 0.0 0.0 - + 22.1 1.4 1.8 - + 9.4 0.7 0.9 - + 0.6 0.02 0.02 - + −1.96 ± 0.05 GD-1 ( ) 164 , 2952 , 938 28 35 61 66 22 23 - + - + - + 1153 14 15 - - + 14.1 0.1 0.1 - + 23.0 1.0 1.1 - + 14.8 0.7 0.7 - + 0.24 0.02 0.02 - + −2.24 ± 0.21 Note. From left to right, the columns provide the stream’s name, action components (J), energy (E), pericentric distance (rperi), apocentric distance (rapo), maximum height of the orbit from the Galactic plane (zmax), 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 Hríd (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). 8 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 Malhan et al. because of Gaia’s limiting magnitude at G ∼ 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 ≈ 0.5, and it is rare for the objects to possess very radial orbits (e ≈ 1) or very circular orbits (e ≈ 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 nonpara- metric 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 Mahala- nobis (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 xi and xj is given by ( ) ( ) ( ) x x x x x x s , . , i j i j T i j 2 = - - 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 Mahalonobis 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 ( ) ∣ ( )∣ ( ) ( ) ( ) ( ) x x x x x x x x x x s , , . . , . , 3 i j i j i j i j i j d T 2 1 1 = S - S - - where d is the dimension of the data, Σ is the covariance matrix, Σ(xi, xj) = 0.5[Σ(xi) + Σ(xj)] and Σ −1(xi, xj) = 0.5[Σ−1(xi) + Σ −1(xj)]. The above formula can be intuitively understood as follows. Consider the term ( ) x x . i j T 1 - S- . Here, (xi − xj) is the distance between two data points. This is then multiplied by the inverse of the covariance matrix Σ (or divided by the covariance matrix). So, this is essentially a multivariate equivalent of the regular standardization y = (x − μ)/σ. 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 xi ≡ (JR,i, Jf,i, Jz,i, Ei), where i represents a given halo object and the units of J and E are kpc km s−1 and km2 s−2, 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’itoi 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 (JR,i, Jf,i, Jz,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 9 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 Malhan et al. first estimates the density and then finds groups in the density field. To search for groups in a d-dimensional data set, ENLINK requires neighbors … (d+ 1). In our case, d= 4 (3 components of J and E) and therefore we set neighbors = 5. Second, we set min_cluster_size = 5. This is because it is difficult to find groups smaller than the smoothing length (i.e., we satisfy the min_cluster_size … neighbors condition of ENLINK). min_peak_height can be thought of as the signal-to-noise ratio of the detected groups, and we set min_peak_height = 3.0. For the parameters cluster_ separation_method and gmetric, we adopt the default values (i.e., 0). Further, we set density_method = 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.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” (PGroup) for each halo object, such that this probability is higher for those objects that belong to the groups detected by ENLINK. To compute this PGroup 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 labels = 1 correspond to those objects that form part of the largest group. For all the objects with labels … 2, we explicitly set their probability of group membership at iteration i to be PGroup,i = 1. Among objects with labels = 1, we accept only those objects that possess density > 99th percentile and set their PGroup,i = 1, while the remaining low density objects are set as PGroup,i = 0. 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 PGroup,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 PGroup of that object. The PGroup 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 PGroup values. Figure 3 shows the (J, E) distribution of the halo objects as a function of the computed PGroup 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 PGroup values. We also note that objects with high PGroup 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 PGroup 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 PThreshold, such that objects with PGroup … PThreshold belong to high-significance groups, and this corresponds to a likely detection. To find a suitable PThreshold value, we follow a pragmatic approach. We repeat the above analysis of comput- ing the PGroup 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 PGroup values (that is shown in Figure 16), its 90 percentile limit motivates setting a threshold at PThreshold = 0.3 for a 2σ 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 PThreshold value is detailed in Appendix B. Consequently, for the real PGroup values (shown in Figure 4), all those objects that possess PGroup … PThreshold are considered as high-significance groups. The selection PGroup … PThreshold 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 15 For example, instead of using the adaptive metric, we defined a constant metric using the uncertainties on (J, E) by setting 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. 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 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. 10 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 Malhan et al. any links, even though they satisfy the condition PGroup … PThreshold, 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 …2σ 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 3. (J, E) distribution of the halo objects as a function of their group probability PGroup (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 PGroup values lie in the denser regions of this (J, E) space. Such objects with high PGroup 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. 11 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 Malhan et al. 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 ∼ [−0.91, −0.79] × 105 km2 s−2, JR ∼ [1205, 2090] km s −1 kpc, Jf ∼ [−2115, −265] km s−1 kpc, Jz ∼ [3285, 5350] km s −1 kpc, L⊥ ∼ [4565, 6835] km s−1 kpc, eccentricity ∼[0.5, 0.6], rperi ∼ [11, 22] kpc, rapo ∼ [45, 60] kpc, and f∼ [66°, 86°], where f 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, De ≈ 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 (Jf, E) values of their stream sample to the (Jf, E) distribution of the mergers previously found by Naidu et al. (2021), but their stream member list does not include Elqui.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 (Mhalo) and stellar mass (M*) as follows. We first determine Mhalo using the globular-cluster-to-halo-mass rela- tion (Hudson et al. 2014) and then convert this Mhalo 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 Mhalo ∼ 4.5 × 10 10Me, and this further implies M*∼ 13 × 10 7Me. 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∼ [−1.09, −0.93]× 105 km2 s−2, JR ∼ [605, 1075] km s−1 kpc, Jf ∼ [−2700, −1360] km s−1 kpc, Jz ∼ [1835, 2820] km s −1 kpc, Figure 4. PDF of the group probability PGroup of all the halo objects (see Section 3.3). The vertical line represents the PThreshold value and all the objects with PGroup … PThreshold 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. 17 Our stream sample contains all the streams that Bonaca et al. (2021) associated with Sagittarius, except for “Turranburra.” 12 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 Malhan et al. L⊥ ∼ [2905, 4635] km s −1 kpc, eccentricity∼ [0.4, 0.6], rperi ∼ [8, 16] kpc, rapo ∼ [31, 43] kpc, and f∼ [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 JR and Jf 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.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 Lz−L⊥ 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 (JR, Jf, E) values to those of the Cetus group but possess lower Jz values, rendering this association rather tentative. We note that NGC 4590 was previously associated with the Helmi sub- structure 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). 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., PGroup … PThreshold; see Section 3.3). The left panel shows Jf vs. JR, and the objects are colored by their Jz values. The right panel shows Jf vs. E, colored by L⊥. The gray points are all the remaining objects with PGroup < PThreshold. 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. 18 These streams are Willka Yaku, Triangulum, and Turbio. 13 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 Malhan et al. 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 Mhalo ∼ 2× 10 10Me, and this in turn provides Table 3 Various Groups Detected in Our Study along with Their Member Globular Clusters, Stellar Streams, and Satellite Galaxies Merger/ No. of Member Member Member In Situ Group Members Globular Clusters Stellar Streams Satellite Galaxies Sagittarius 8 Pal 12, Whiting 1, Terzan 7, Terzan 8, Elqui Sagittarius dSph (Section 4.1) NGC 6715/M54, Arp 2 Cetus 6–8 NGC 5824 Cetus [stream of Cetus], Willman 1 (Section 4.2) Slidr, Atlas, AliqaUma tentative: NGC 4590/M68 [stream:Fjörm] Fjörm [stream of NGC 4590/M68] Gaia-Sausage/Enceladus 16–18 NGC 7492, NGC 6229, NGC 6584, C-7, Hrìd, (Section 4.3) NGC 5634, IC 1257, NGC 1851, M 5 [stream of NGC 5904/M5] NGC 2298, NGC 4147, NGC 1261, NGC 6981/M72, NGC 1904/M79, NGC 7089/M2 [stream:NGC 7089], NGC 5904/M5 [stream: M5] tentative: NGC 6864/M75 NGC 7089 [stream of NGC 7089/M2] Arjuna/Sequoia/I’itoi 9 NGC 6101, GD-1, Phlegethon, Gaia-9, Kshir, (Section 4.4) NGC 3201 [streams: NGC 3201, Gjöll] NGC 3201 [stream of NGC 3201], Gjöll [stream of NGC 3201], Ylgr LMS-1/Wukong 11 NGC 5272/M3, NGC 5053, LMS-1 [stream of LMS-1/Wukong], (Section 4.5) NGC 5024/M53, C-19, Sylgr, Phoenix, Indus, Pal 5 [stream: Pal 5] Jhelum, Pal 5 [stream of Pal 5] Pontus 8–10 NGC 288, NGC 5286, NGC 7099/M30 M92 [stream of NGC 6341/M92] (Section 4.6) NGC 6205/M13, NGC 6779/M56, NGC 6341/M92 [stream: M92], NGC 362 tentative: NGC 6864/M75 NGC 7089 [stream of NGC 7089/M2] Candidate merger 5 NGC 5466 [stream: NGC 5466], Gaia-10, Tucana III (not detected, but NGC 7492 NGC 5466 [stream of NGC 5466] selected, Section 5) Galactic disk 6-7 Pal 10, NGC 6838/M71, NGC 6356, (Section 4.7) IC 1276/Pal 7, Pal 11, NGC 104/47 Tuc tentative: NGC 7078/M15 Galactic Bulge 28 Terzan 2/HP 3, 1636-283/ESO452, Gran 1, (Section 4.7) 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 Galactic Bulge/ 11 Terzan 3, NGC 6569, NGC 6366, NGC 6139, disk/low energy BH 261/AL 3, NGC 6171/M107, Pal 8, (Section 4.7) Lynga 7/BH 184, NGC 6316, FSR 1716, NGC 6441 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). 14 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 Malhan et al. M*∼ 3× 10 7Me. 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 |Jf| and L⊥, implying that they lie along radial orbits. This group possesses dynamical properties in the range E∼ [−1.44, −1.16]× 105 km2 s−2, JR ∼ [935, 2075] km s −1 kpc, Jf ∼ [−715, 705] km s −1 kpc, Jz ∼ [85, 1505] km s −1 kpc, L⊥ ∼ [85, 1520] km s−1 kpc, eccentricity∼ [0.7, 0.9], rperi ∼ [1, 4] kpc, rapo ∼ [16, 30] kpc, and f∼ [27°, 85°]. 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, Hrìd) 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 JR 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≈32 globular clusters to this merger, although some of their associations were tentative. They found these association by making hard cuts in the (Jf, E, L⊥) 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 Mhalo ∼ 10× 10 10Me and M*∼ 50× 10 7Me. 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 ∼ [−1.27, −1.02] × 105 km2 s−2, JR ∼ [20, 1190] km s−1 kpc, Jf ∼ [1880, 2955] km s−1 kpc, Jz ∼ [230, 940] km s −1 kpc, L⊥ ∼ [980, 2530] km s−1 kpc, eccentricity ∼ [0.1, 0.6], rperi ∼ [5, 14] kpc, rapo ∼ [15, 37] kpc, and f ∼ [25°, 46°]. 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 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 PGroup … PThreshold, and the remaining objects are shown as gray points. The left plot shows Jf vs. E, and the right plot shows rperi vs. rapo. The LMS-1/ Wukong group has a minimum of [Fe/H] = −3.4 dex. 15 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 Malhan et al. and I’itoi (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, Jf) 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’itoi). 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’itoi 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]× 105 km2 s−2, JR ∼ [54, 1400] km s −1 kpc, Jf = [250, 3210] km s−1 kpc, and Jz = [66, 800] km s −1 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’itoi group, especially in Jf 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’itoi globular clusters were pre- viously associated with both Sequoia and Gaia–Sausage/ Enceladus by Massari et al. (2019); we note that their (Jf, 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’itoi 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’itoi, 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’itoi galaxy, perhaps that of the Sequoia galaxy (given the similarity in their [Fe/H]). Using only the member globular clusters of Arjuna/ Sequoia/I’itoi and not the streams, we estimate the mass of the progenitor galaxy to be Mhalo ∼ 1.2 × 10 10Me and M*∼ 1.5 × 10 7Me. 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’itoi group. 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. 16 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 Malhan et al. 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∼ [−1.41, −1.19] × 105 km2 s−2, JR ∼ [100, 605] km s −1 kpc, Jf ∼ [−1560, − 210] km s−1 kpc, Jz ∼ [875, 2710] km s−1 kpc, L⊥ ∼ [1400, 3085] km s−1 kpc, eccentricity ∼ [0.2, 0.5], rperi ∼ [5, 13] kpc, rapo ∼ [15, 25] kpc, and f ∼ [58°, 85°]. This polar group corr- esponds 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/Wukongwas 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/ Wukongmerger. 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] ∼−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 Mhalo ∼ 2.7 × 10 10Me and M*∼ 5.5 × 10 7Me. 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∼ [−1.72, −1.56]× 105 km2 s−2, JR ∼ [245, 725] km s −1 kpc, Jf ∼ [−5, 470] km s −1 kpc, Jz ∼ [115, 545] km s −1 kpc, L⊥ ∼ [390, 865] km s−1 kpc, eccentricity∼ [0.5, 0.8], rperi ∼ [1, 3] kpc, rapo ∼ [8, 13] kpc, and f∼ [54°, 89°]. We refer to this group as Pontus.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 add- itional objects close to this group, namely the globular cluster NGC 6864/M75 and the stream NGC 7089, but their associa- tion was not very strong (because of their slightly higher JR and slightly lower Jf values). Pontus lies close to Gaia–Sausage/Enceladus in (Jf, E) space (although the two groups possess very different JR values) and essentially all of Pontus’s globular clusters (which we mentioned above) have been previously associated with Figure 8. Comparing the orbital properties of those objects that belong to different groups. The left panel shows rperi vs. eccentricity, and the right panel shows rperi vs. rapo. 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). 19 In Greek mythology, “Pontus” (meaning “the Sea”) is the name of one of the first children of the deity Gaia . 17 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 Malhan et al. 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 eccentri- city 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 JR values). Also, the average rapo 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 Vr and azimuthal Vf (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 Vr direction (implying radial orbits), although their Vf components (on average) differ by ≈60 km s−1. 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 Vf 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, Jf)∼ (–1.65× 10 5 km2 s−2, 1500 km s−1 kpc) and∼(–1.75× 105 km2 s−2, 900 km s−1 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, Jf)∼ (–1.75× 10 5 km2 s−2, 500 km s−1 kpc) (see their Figure 23). Given that these (E, Jf) 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 Mhalo ∼ 5× 1010Me and M*∼ 15 × 10 7Me. 4.7. The in Situ Groups 7, 8, and 9 We detect three additional groups; however, their locations in (E, Jf) space indicates that they do not represent any merger 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 Vr and azimuthal Vθ. The filled ellipses represent the 1.5σ confidence contour and the Gaussians represent the mean and the standard deviation in the Vf components of the member objects of these groups. 18 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 Malhan et al. 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 rapo, 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 ∼ [−1.77, −1.66] × 105 km2 s−2, JR ∼ [15, 135] km s−1 kpc, Jf ∼ [−1315, −960] km s −1 kpc, Jz ∼ [5, 235] km s−1 kpc, L⊥ ∼ [115, 730] km s −1 kpc, eccentricity ∼[0.1, 0.4], rperi ∼ [3, 6] kpc, rapo ∼ [7, 9] kpc, and f ∼ [5°, 36°]. Given these dynamical properties, especially low eccentricity (implying circular orbits), low f value (implying that the objects orbit close to the Galactic plane), and the values of rperi and rapo 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]∼−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]… −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∼ [−2.55, −2.22]× 105 km2 s−2, JR ∼ [5, 140] km s−1 kpc, Jf ∼ [−380, 150] km s −1 kpc, Jz ∼ [0, 275] km s−1 kpc, L⊥ ∼ [10, 405] km s −1 kpc, eccentricity∼ [0.1, 0.8], rperi ∼ [0, 2] kpc, rapo ∼ [1, 4] kpc, and f∼ [1°, 89°]. This group comprises 28 globular clusters (their names are provided in Table 3). Given these dynamical properties, especially very low values of E, rperi, and rapo, and that the objects are spherically distributed (as we note from the range of the f 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]∼−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∼ [−2.17, −1.92]× 105 km2 s−2, JR ∼ [10, 110] km s −1 kpc, Jf ∼ [−630, −250] km s−1 kpc, Jz ∼ [55, 190] km s −1 kpc, L⊥ ∼ [215, 485] km s−1 kpc, eccentricity∼ [0.2, 0.4], rperi ∼ [1, 3] kpc, rapo ∼ [3, 6] kpc, and f∼ [18°, 56°]. This group contains 11 globular clusters (these are listed in Table 3). The metallicity of these objects range from [Fe/H]∼−1.65 dex to −0.4 dex with a median of [Fe/H]=−0.7 dex. For this group, while its dyna- mical properties appear consistent with those of the disk (e.g., low eccentricity, the rapo range, and low f 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 Jf values (on average), while Kraken objects have an average Jf ∼ 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 …2σ 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 PGroup 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 „ 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∼ [−1.2, −0.94]× 105 km2 s−2, JR ∼ [1385, 2525] km s−1 kpc, Jf ∼ [130, 915] km s−1 kpc, Jz ∼ [1125, 2095] km s −1 kpc, L⊥ ∼ [1435, 2795] km s −1 kpc, eccentricity∼ [0.7, 0.8], rperi ∼ [3, 7] kpc, rapo ∼ [27, 47] kpc, and f∼ [69°, 85°]. 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 f 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]∼−1.98 dex (Lamb et al. 2015) and NGC 7492 has [Fe/H]∼−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].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 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 ∼ 21. 19 The Astrophysical Journal, 926:107 (30pp), 2022 February 20 Malhan et al.