Evidence for Multiple Accretion Events in the Gaia-Sausage/Enceladus Structures

The studies of the stellar halo of the Milky Way (MW), which is thought to be assembled via multiple hierarchical mergers (White & Frenk 1991), provide valuable clues to its formation and evolutionary history, as its long dynamical timescale preserves a fossil record of past accretion events (Bland-Hawthorn & Gerhard 2016). The advent of the Gaia Data Releases (Gaia DRs; Gaia Collaboration et al. 2016, 2018), which provide precise astrometric information for many millions of stars, has dramatically expanded the detail of our view of the MW's accretion history. For instance, the combination of Gaia and large spectroscopic survey data has enabled the detection of distinctive accretion signatures from the Gaia-Sausage (GS; Belokurov et al. 2018) and Gaia-Enceladus (GE; Helmi et al. 2018). Other small-scale accretion events have also been discovered (Myeong et al. 2018, 2019; Koppelman et al. 2019; Naidu et al. 2020; Necib et al. 2020; Yuan et al. 2020; Horta et al. 2021; Re Fiorentin et al. 2021).

According to the study of Belokurov et al. (2018), who used a sample of main-sequence (MS) stars from the Sloan Digital Sky Survey (SDSS; York et al. 2000) along with Gaia DR1 astrometry, the GS structure exhibits strong radial anisotropy and a mildly prograde motion of 20 ∼ 30 km s⁻¹ for stars with −1.7 < [Fe/H] < −1.0. Based on cosmological zoom-in simulations of the formation of the stellar halo, they argued that the GS is the result of an accretion event of a massive dwarf galaxy with orbital eccentricity e > 0.9. From a sample of nearby MS and blue horizontal-branch stars with SDSS DR9 (Ahn et al. 2012) spectroscopy and Gaia DR2 proper motions, Deason et al. (2018) demonstrated that the GS stars with −1.5 < [Fe/H] < −1.0 have eccentricities e > 0.9 in the range of $10\lt {r}_{\max }$ (apogalactic distance) <30 kpc, and an average ${r}_{\max }$ agrees with the break radius of the MW stellar halo (Deason et al. 2013).

From an analysis of the kinematics, chemical abundances, ages, and spatial distributions of disk and halo stars in Gaia DR2 with available spectroscopy from the Apache Point Observatory Galactic Evolution Experiment (APOGEE) DR14 (Majewski et al. 2017; Abolfathi et al. 2018), Helmi et al. (2018) identified that GE stars primarily occupy the inner-halo region (IHR), and exhibit a small retrograde net motion. They also suggested that the GE stars are debris from a massive dwarf galaxy, inferred from their relatively low [α/Fe] and a large spread in metallicity. Mackereth et al. (2019) additionally found a robust correlation between the chemical abundances and the orbital eccentricities of local halo stars using the APOGEE DR14 data. Their results indicated that the majority of local halo stars have highly radial orbits with e > 0.7 (see also Mackereth & Bovy 2020), and relatively low abundances of [Mg/Fe], [Al/Fe], and [Ni/Fe], which are comparable to those of stars observed in massive surviving satellite galaxies of the MW.

The general consensus of the above studies is that both the GS and GE stars exhibit eccentric and radial orbits in the inner halo, and they are possibly the remnants of a single massive (∼10⁹ M_⊙) disrupted galaxy. Nonetheless, closer inspection has suggested that they have different rotational motions, and that the GS stars appear to have e > 0.9, higher than the GE stars (Belokurov et al. 2018; Deason et al. 2018; Helmi et al. 2018; Mackereth et al. 2019). Thus, it is still an open question whether or not the GS and GE stars have originated from a single accretion event.

Recent studies have demonstrated that the dynamical properties of likely progenitors for different substructures can be a powerful tool to distinguish one accretion episode from another. For example, studies of minor-merging simulations (Read et al. 2008; Villalobos & Helmi 2008; Jean-Baptiste et al. 2017; Karademir et al. 2019) showed that the key ingredients of merged galaxies to understand the orbits of stripped stars are their orbital eccentricity, inclination, and inferred mass. In particular, the orbital eccentricities of disrupted stars largely remain unchanged after accretion (Mackereth et al. 2019), and the majority of these stars orbit with the same orbital inclinations as of their parent galaxies (Re Fiorentin et al. 2015).

Retrograde motions of stars in the outer halo of a galaxy have been shown by merging simulations to arise from dwarf parent galaxies on low-inclination retrograde orbits (Bignone et al. 2019). Two massive dwarf galaxies with different orbital directions on high-eccentricity orbits (Murante et al. 2010) can also produce stars with retrograde motions. Observations also confirm these predictions by various simulations. Helmi et al. (2018) found, from isolated simulations of minor mergers (e.g., Villalobos & Helmi 2009), that the retrograde motion of the GE stars is similar to that of a retrograde encounter with a low-inclination dwarf. In addition, Simion et al. (2019) verified that the Hercules-Aquila cloud and Virgo overdensity are dominated by stars on highly eccentric orbits, which are commensurate with the kinematic and orbital properties of the GS stars. They also showed that both of the diffuse debris clouds associated with this structure have high orbital-inclination trajectories. These studies imply that dwarf galaxies with different orbital inclination leave distinct dynamical signatures in their disrupted stars; thus the orbital inclination of their disrupted stars can be used to trace the orbital properties of their progenitors.

Following the above reasoning, in this Letter we first identify two groups of stars—low orbital-inclination (LOI) and high orbital-inclination (HOI) stars—among stars with high eccentricity (e > 0.7) and [Fe/H] < −1.0, the typical properties of the GS and GE (GSE, hereafter) stars, and report on the distinct chemical and dynamical signatures of LOI and HOI groups, providing strong evidence for multiple accretion events involved in the formation of the GSE structures.

This Letter is organized as follows. We explain our sample selection in Section 2, and calculations of velocity components and orbital parameters of our sample stars in Section 3. In Section 4, we examine kinematic and orbital properties of LOI and HOI stars in the IHR, and present evidence for accretion events that are distinct from the GSE structures, as well as the identification of retrograde motions associated with stars in the outer-halo region (OHR). Section 5 discusses the implications of our findings; a summary follows in Section 6.

We have collected a sample of stars with available medium-resolution (R ∼ 2000) spectra from SDSS DR12 (Alam et al. 2015), which includes objects from the legacy SDSS program, the Sloan Extension for Galactic Understanding and Exploration (SEGUE; Yanny et al. 2009), and the Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al. 2013), covering the extinction-corrected magnitude and dereddened color ranges 14.0 < g₀ < 20.0 and 0.0 < (g − r)₀ < 1.2, respectively.

Using the SEGUE Stellar Parameter Pipeline (Allende Prieto et al. 2008; Lee et al. 2008a, 2008b), we estimated stellar atmospheric parameters (T_eff, $\mathrm{log}\,g$, and [Fe/H]) for each star. In order to obtain reliable stellar-parameter estimates for MS and MS turnoff (MSTO) stars, we restrict our analysis to stars with average spectral signal-to-noise ratios (S/Ns) greater than 10.0, 4400 ≤ T_eff ≤ 7000 K, and $\mathrm{log}\,g\geqslant 3.5$. For stars with multiple observations, we chose the star with the highest S/N, and removed stars with apparently defective spectra.

Proper motions for stars with errors less than 1.0 mas yr⁻¹ were obtained through cross-matching with Gaia DR2 (Gaia Collaboration et al. 2018). Radial velocities were adopted from the SDSS pipeline; these have a precision of 5–20 km s⁻¹, depending on the S/N of the spectrum, and negligible zero-point errors (Yanny et al. 2009). For stellar-distance estimates, we employed the methods of Beers et al. (2000, 2012), as our program stars are mostly too faint to have reliable parallaxes available from Gaia DR2. Their reported uncertainty is on the order of 15%–20%, as verified by comparing our derived distances with Gaia DR2 distances based on parallaxes with relative errors less than 10% (Kim et al. 2019; Lee et al. 2019).

Given the distances, radial velocities, and proper motions adopted for our sample of stars, we derived their space velocity components in a spherical coordinate system. For these calculations, we adopted V_LSR = 236 km s⁻¹ (Kawata et al. 2019) for the rotation of the local standard of rest (LSR), a solar position of R_⊙ = 8.2 kpc (Bland-Hawthorn & Gerhard 2016) in the disk plane from the Galactic center, and a vertical distance of Z_⊙ = 20.8 pc (Bennet & Bovy 2019) from the midplane. The solar peculiar motion with respect to the LSR was assumed to be (U, V, W)_⊙ = (−11.10, 12.24, 7.25) km s⁻¹ (Schönrich et al. 2010), where the velocity components U, V, and W are positive in the direction toward the Galactic anticenter, Galactic rotation, and north Galactic pole, respectively. In our adopted system, a star with V_ϕ > 0 km s⁻¹ has a prograde motion; retrograde rotation is indicated by V_ϕ < 0 km s⁻¹. Stars with V_r > 0 km s⁻¹ move away from the Galactic center, and stars with V_θ > 0 km s⁻¹ move toward the south Galactic pole.

We also made use of a Galactocentric Cartesian reference frame, denoted by (X, Y, Z), where the axes are positive in orientation toward the Sun, Galactic rotation, and north Galactic pole, respectively. In addition, we introduced an angle (α) between the orientation of the total angular momentum vector and the L_Z axis, and a position angle (Θ) of L_⊥ measured from the negative L_X axis to the positive L_Y -axis direction of the angular momentum component in the X–Y plane. These are defined by:

$$\begin{eqnarray*}&&\alpha ={\cos }^{-1}\left(\displaystyle \frac{{L}_{Z}}{L}\right)\,\mathrm{and}\,{\rm{\Theta }}={\tan }^{-1}\left(-\displaystyle \frac{{L}_{Y}}{{L}_{X}}\right),\end{eqnarray*}$$

where $L=\sqrt{{L}_{\perp }^{2}+{L}_{Z}^{2}}$ is the total angular momentum, ${L}_{\perp }=\sqrt{{L}_{X}^{2}+{L}_{Y}^{2}};$ L_X , L_Y , and L_Z are the X, Y, and Z components of the angular momentum, and are positive along the positive X-axis, the positive Y-axis, and the negative Z-axis directions, respectively. Schematic representations for α and Θ are shown in Figure 1. In this notation, stars with L_Z > 0 have prograde orbits and α < 90°. Retrograde orbits have α > 90°, and the inclination angle (i) of their orbital plane increases as α approaches 90°. For prograde motions, the inclination angle i = α, whereas i = 180° − α for retrograde orbits.

Zoom In Zoom Out Reset image size

We employed an analytic Stäckel-type potential (see Chiba & Beers 2000; Kim et al. 2019 for details) in order to calculate the orbital parameters of our sample stars, including the perigalactic distance (${r}_{\min }$, the minimum distance of an orbit from the Galactic center), apogalactic distance (${r}_{\max }$, the maximum distance of an orbit from the Galactic center), and stellar orbital eccentricity ($e=({r}_{\max }-{r}_{\min })/({r}_{\max }+{r}_{\min })$), as well as ${Z}_{\max }$ (the maximum distance of a stellar orbit above or below the Galactic plane).

Uncertainties on the derived kinematic and orbital values are obtained from 100 Monte Carlo simulations with adopted uncertainties of 20% in distance, and quoted uncertainties in the radial velocity and proper motions, assuming Gaussian error distributions. Prior to examining the metallicity distribution function (MDF) of our program stars, we derived simple selection functions to correct for the target-selection bias. The selection function is defined by the fraction of the spectroscopically targeted stars among the photometrically available targets in bins of 0.2 and 0.05 mag for a color–magnitude diagram of r₀ and g₀ − r₀, separately for each SDSS plug-plate, as described in Lee et al. (2019).

Even though we obtained a total number of N = 328, 102 stars with valid orbital parameters, in this study we focus on the lower metallicity stars with high-eccentricity orbits, in the range of [Fe/H] < −1.0 and e > 0.7, in order to reduce contamination from disk stars heated by the GSE (Belokurov et al. 2020), and include the GSE stars according to Mackereth et al. (2019), resulting in a total of 68,776 stars to be explored. Note that all of our selected stars are also included in the GSE even using the selection criteria by Naidu et al. (2020).

In this section, we analyze the kinematic and orbital properties of stars with highly eccentric and radial orbits in two regions: the IHR ($15\leqslant {r}_{\max }\lt 30\,\mathrm{kpc}$) and the OHR (${r}_{\max }$ ≥ 30 kpc). In the IHR, we search for any differences in the dynamical properties between the LOI and HOI stars, while in the OHR, we explore the mean rotational velocity as a function of mean ${r}_{\max }$, and the stellar MDFs for the LOI and HOI populations.

4.1. Definition of LOI and HOI Populations

The top panel of Figure 2 presents the logarithmic number density of the selected program stars in the Θ versus α plane. In this map, we clearly see distinct kinematic features. Stars with relatively HOI of 55° < α < 125° are mostly concentrated around Θ = 90° and 270°, and have the highest density at α = 90°, whereas stars with relatively LOI of α ≤ 55° or α ≥ 125° are distributed over all ranges of Θ. Moreover, the α histogram (black color) for our program stars, shown in the bottom panel of Figure 2, is well-fit with the sum (cyan line) of three different Gaussian distributions, which represent one (blue) for the HOI and two (red) for the LOI population. By assembling these characteristics together, we define LOI stars as those in the range of α ≤ 55° or α ≥ 125°; HOI stars are those in the range of 55° < α < 125°. The condition of 55° < α < 125° is equivalent to ${L}_{\perp }^{2}/{L}^{2}\gt {\sin }^{2}55^\circ =0.67$, or ${L}_{\perp }\gt \pm \tan 55^\circ {L}_{Z}=\pm 1.43{L}_{Z}$ in the L_Z versus L_⊥ space.

Zoom In Zoom Out Reset image size

4.2. Separation of GSE Stars into Low and High Orbital-inclination Groups

First, in order to select GSE stars having highly eccentric and radial orbits in the IHR, we plot stars in the ranges of $15\leqslant {r}_{\max }\lt 30\,\mathrm{kpc}$ (Deason et al. 2018) and [Fe/H] < −1.0 in the plane of ${V}_{{\rm{r}}}^{2}/{V}^{2}$ versus orbital eccentricity, as shown in the top panel of Figure 3. In this figure, cyan, magenta, and orange dots represent stars with e > 0.7 on radial and tangential orbits, and stars with e ≤ 0.7, respectively. Then, we choose GSE members (cyan dots in the panel), defined by e > 0.7 (Mackereth et al. 2019) and ${V}_{{\rm{r}}}^{2}/{V}^{2}$ > 0.6 (Belokurov et al. 2018), where ${V}_{{\rm{r}}}^{2}/{V}^{2}$ > 0.6 is equal to β > 0.67. The velocity anisotropy parameter (β) of each star is defined as β = 1–${V}_{{\rm{t}}}^{2}/2{V}_{{\rm{r}}}^{2}$ (see Binney & Tremaine 2008; Elias et al. 2020).

Zoom In Zoom Out Reset image size

Next, we divide the GSE members into two groups of stars: the LOI and HOI subsamples, by following the definitions described in Section 4.1. The bottom panel of Figure 3 shows the ${r}_{\max }$ distribution of the LOI (red) and HOI (blue) populations, respectively. As can be seen, there are nearly equal numbers of stars for the LOI and HOI populations as a function of ${r}_{\max }$. This indicates that their progenitor galaxies are likely different, as it is unlikely to have equal numbers of LOI and HOI stars if they are accreted from a single dwarf galaxy.

4.3. Inner-halo Region

The orange solid box in the top-right panel of Figure 4 indicates the GE stars identified in Helmi et al. (2018). We note in the panel that the GE stars in the box are mostly made up of LOI and HOI stars with retrograde motions. Furthermore, we can observe that the HOI stars (blue histogram) have mostly e > 0.9 in the bottom right panel of Figure 4, which is characteristic of the GS stars reported by Deason et al. (2018). Thus, we can confirm that most of our selected LOI (red dots) and HOI (blue dots) stars follow typical properties of the GSE stars in the E_tot versus L_Z plane, although we applied slightly different selection criteria for the GSE member stars from the original works (e.g., Helmi et al. 2018). As a result, our selection criteria for the GSE members do not alter the interpretation of the dynamical properties of the genuine GSE stars.

Zoom In Zoom Out Reset image size

The top-left panel of Figure 4 shows that the LOI and HOI subsamples have essentially identical distributions of V_r. However, the V_ϕ distribution (bottom left panel) and the e-distribution (bottom right panel) tell an entirely different story—the eccentricity distributions exhibit completely different behaviors between the LOI (red) and HOI (blue) populations. The LOI stars (red) are well-separated into two subgroups of prograde and retrograde motion in the V_ϕ distribution, and the LOI group mainly occupies the range of e < 0.9, while the majority of stars in the HOI population have e > 0.9. These properties provide clear evidence that these stars have experienced different accretion episodes. The HOI stars (blue) have almost zero net rotational velocity, suggesting that they share a common progenitor.

We can use other orbital parameters, e.g., ${Z}_{\max }$ and ${r}_{\max }$, to explore more the discrete dynamical signatures of the LOI and HOI populations. To accomplish this, we first introduce an angle, $\phi ={\tan }^{-1}({Z}_{\max }/\sqrt{{r}_{\max }^{2}-{Z}_{\max }^{2}})$. As illustrated in the top panel of Figure 5, which shows the distribution in ${Z}_{\max }$ versus $\sqrt{{r}_{\max }^{2}-{Z}_{\max }^{2}}$ for LOI (red dots) and HOI (blue dots) stars on highly eccentric and radial orbits, this angle is measured from the axis of $\sqrt{{r}_{\max }^{2}-{Z}_{\max }^{2}}$ to a line connecting the origin to a star in the coordinate plane of $\sqrt{{r}_{\max }^{2}-{Z}_{\max }^{2}}$ and ${Z}_{\max }$. Generally, this angle (ϕ) increases with increasing ${Z}_{\max }$ at a given ${r}_{\max }$ (or E_tot).

Zoom In Zoom Out Reset image size

The middle and bottom panels of Figure 5 show the distributions of the ϕ angles for LOI (red) and HOI (blue) stars, in the two regions of ${r}_{\max }$: all ${r}_{\max }$ (middle panel) and $15\leqslant {r}_{\max }\lt 30\,\mathrm{kpc}$ (bottom panel). These panels immediately indicate that the HOI stars have mostly higher ϕ values than those of the LOI stars, both over all ${r}_{\max }$ and $15\leqslant {r}_{\max }\lt 30\,\mathrm{kpc}$ ranges. Once again, this discrepancy suggests a dynamical distinction between the two populations. One can assume that if both HOI and LOI stars were accreted from a relatively massive dwarf galaxy, as claimed in the literature (e.g., Belokurov et al. 2018; Helmi et al. 2018), no matter what the orbital inclination they have, both groups of stars should exhibit similar ϕ distributions, regardless of their ${r}_{\max }$.

4.4. Outer-halo Region

We now extend our search for diverse dynamical signatures between LOI and HOI stars to the OHR, the region with ${r}_{\max }$ ≥ 30 kpc. Once again, we only consider stars with [Fe/H] < −1.0 on highly eccentric (e > 0.7) orbits, and apply the selection criteria for LOI and HOI stars described in Section 4.1 to the OHR. First, we search for evidence of retrograde motions and dissimilar MDFs for stars in the OHR. The left plot of top panels in Figure 6 shows the profiles of mean rotation velocity, as a function of ${r}_{\max }$, for the full sample of stars (brown), and for stars with ${V}_{{\rm{r}}}^{2}/{V}^{2}\gt 0.6$ (black). The middle plot shows the V_ϕ distributions for stars with ${r}_{\max }$ ≥ 30 kpc. The right plot represents the MDFs for the stars with ${r}_{\max }$ ≥ 30 kpc and V_ϕ < 0 km s⁻¹. The black and magenta histogram in the middle and right plots in each panel denote stars with ${V}_{{\rm{r}}}^{2}/{V}^{2}\gt 0.6$ (radial-motion dominated) and ${V}_{{\rm{r}}}^{2}/{V}^{2}\leqslant 0.6$ (tangential-motion dominated), respectively. Plots in the top panels are for the LOI population, whereas bottom panels are the same as in the top panels, but for the HOI population.

Zoom In Zoom Out Reset image size

Inspection of the left plot of the top panels shows that the LOI stars with strong radial motions (black) exhibit stronger retrograde motions in the OHR than in the IHR, and drops down to V_ϕ ∼ −70 km s⁻¹ at ${r}_{\max }$ ∼ 75 kpc. This behavior can also be inferred from the V_ϕ distribution (black histogram) in the middle plot, which presents more stars with retrograde motion than in prograde motion. By way of contrast, the HOI group of stars with strong radial motion (black) in the left plot of bottom panels exhibits a small retrograde motion, V_ϕ ∼ −15 km s⁻¹, in the OHR, and does not change with increasing ${r}_{\max }$ up to ∼ 85 kpc. The middle plot indicates that their V_ϕ distribution (black histogram) appears to be symmetric around V_ϕ ∼ 0 km s⁻¹. As a result, we realize that the behaviors in the V_ϕ distributions for both LOI and HOI stars with ${V}_{{\rm{r}}}^{2}/{V}^{2}\gt 0.6$ in the OHR are not much different in the IHR. One more interesting aspect is that we generally see a steeper gradient (brown) of V_ϕ over ${r}_{\max }$ in the OHR in the left plot of top and bottom panels. This is driven by the stars on the tangential orbits (${V}_{{\rm{r}}}^{2}/{V}^{2}\leqslant 0.6$) with retrograde motion.

We now examine the MDF of stars with retrograde motions in the OHR. From the right plot of top and bottom panels of Figure 6, one can clearly see different MDFs between the radial- and tangential-dominated samples for the LOI and HOI stars. The peaks of the MDFs for LOI stars on highly radial and tangential orbits are at [Fe/H] = −1.5 and −1.7, respectively, while the peaks of the HOI stars are at [Fe/H] = −1.3 and −1.9, respectively. Converting to the stellar masses of their progenitors by the mass–metallicity relationship for dwarf galaxies (Kirby et al. 2013), we obtained M_* ∼ 2 × 10⁶ M_⊙ and M_* ∼ 4 × 10⁵ M_⊙ for LOI stars, and M_* ∼ 3 × 10⁶ M_⊙ and M_* ∼ 3 × 10⁵ M_⊙ for HOI stars, respectively. Similarly, the stellar mass for the progenitor of the LOI stars with prograde motion on highly radial orbits is M_* ∼ 3 × 10⁶ M_⊙. Consequently, the discrepancies in the MDFs and progenitor masses suggest that the progenitor of the stars on highly radial orbits experienced different star formation histories from that of those on tangential orbits, for both the LOI and HOI subsamples.

In the previous section, we showed that the LOI and HOI stars in the IHR exhibit different distributions in V_ϕ , e, and the ϕ angle. The LOI stars include more objects with e < 0.9 and small ϕ angles, and they are well-separated into stars with prograde and retrograde motions. On the other hand, the HOI stars have more objects with e > 0.9 and high ϕ angles, and they exhibit a single V_ϕ -distribution with a peak at V_ϕ ∼ 0 km s⁻¹.

For stars in the OHR, the objects with tangential orbits among the LOI population exhibit stronger retrograde motions than the HOI population. Considering the stars with retrograde motions in the OHR, the tangential-dominated sample exhibits a more metal-poor distribution than the radial-dominated one for both the LOI and HOI populations.

How do we interpret our findings in terms of the accretion history of the Galactic halo? Let us first consider the dynamical signatures of merging galaxies predicted from various numerical simulations. It is well-known that tidally stripped stars follow the orbital trajectories of their parent galaxies when they are fully accreted into the MW (e.g., Quinn & Goodman 1986; Van den Bosch et al. 1999; Amorisco 2017)—the accreted stars preserve the orbital eccentricity (Mackereth et al. 2019) and inclination (Re Fiorentin et al. 2015) of their parent galaxies. Numerical simulations also predict stars with retrograde motions in the OHR of a galaxy, accreted not only from the merging of two massive dwarf galaxies with different orbital directions on highly eccentric and low-inclination orbits (Murante et al. 2010), owing to the dissimilar efficacy of dynamical friction (Quinn & Goodman 1986), but also from one merged dwarf galaxy on a retrograde orbit of low inclination (Bignone et al. 2019). In terms of the spatial distribution of accreted stars, stars of a low-mass dwarf galaxy under the influence of very weak dynamical friction are tidally stripped off in the outer region of the MW, but its high eccentricity causes its stars to be deposited in the inner region of the MW (Karademir et al. 2019), whereas, due to stronger dynamical friction, a more-massive dwarf galaxy loses its stars in the inner region of the MW (Amorisco 2017).

Following the above reasoning, we can infer that our LOI stars on highly eccentric and radial orbits may be accreted from two massive dwarf galaxies with highly eccentric and low-inclination orbits in prograde and retrograde motions, respectively. Stars stripped from them might drive the prograde motion in the IHR and retrograde motion in the OHR, due to different dynamical friction efficacy arising from the different orbital directions. On the other hand, the HOI stars on highly eccentric and radial orbits are regarded as tidally stripped ones from one massive dwarf galaxy merged on an extremely eccentric and high-inclination orbit under the influence of dynamical friction.

Meanwhile, the LOI and HOI stars with retrograde motions on highly eccentric and tangential orbits have the peaks of their MDFs located in the more metal-poor regime in the OHR than the radial-dominated stars. In addition, many of those stars have higher E_tot than stars on radial orbits. These aspects lead us to conclude that they may be stripped in the OHR of the MW, due to weak dynamical friction and self-gravity from low-mass dwarf galaxies that are on low- and high-inclination orbits with high eccentricity and retrograde motion at high E_tot, resulting in stripped stars that have large ${r}_{\max }$. Taken as a whole, accounting for the dynamical characteristics of our LOI and HOI stars require at least five different accretion episodes.

To demonstrate that the LOI and HOI stars we analyze are different entities from other known substructures, we compared our two groups (LOI and HOI) of stars to those on retrograde orbits identified in the literature (Koppelman et al. 2019; Myeong et al. 2019; Naidu et al. 2020) in the L_Z –E_tot diagram shown in Figure 7. The figure is drawn with the stars with ${r}_{\max }$ < 200 kpc, calculated using the McMillan (2017) potential for easy comparison with those in the literature. Substructures of stars in the LOI and HOI subsamples we have identified are represented in red and blue colors, respectively. The Sequoia event, with 0.5 < e ≤ 0.7, ${V}_{{\rm{r}}}^{2}/{V}^{2}\leqslant 0.45$, and ${r}_{\max }$ ≥ 20 kpc, and Thamnos 1, with e ≤ 0.3 and ${r}_{\max }$ < 12 kpc, are displayed with orange and magenta colors, respectively. The cyan color marks subsamples with e ≤ 0.5 and ${r}_{\max }$ ≥ 20 kpc, similar to those reported by Naidu et al. (2020). From inspection, our LOI and HOI stars on retrograde motions are clearly well-separated from the other known substructures.

Zoom In Zoom Out Reset image size

Using metallicities, radial velocities, and distances from SDSS DR12 and proper motions from Gaia DR2, we have presented an analysis of the kinematic and orbital properties for MS and MSTO stars with eccentricities greater than 0.7 and [Fe/H] < −1.0, which are the typical properties of the GSE stars, after separating them into LOI and HOI subsamples.

LOI and HOI stars on highly eccentric and radial orbits exhibit different dynamical characteristics. LOI stars mostly have e < 0.9 and lower ${Z}_{\max }$, whereas HOI stars have e > 0.9 and higher ${Z}_{\max }$. Moreover, LOI stars are separated into two groups with prograde and retrograde motions, and exhibit prograde motions in the inner halo and retrograde motions in the outer halo. Accordingly, they are regarded as stars accreted from two massive dwarf galaxies with prograde and retrograde orbits of low inclination under the influence of different dynamical friction due to different orbital directions. On the other hand, HOI stars globally have a symmetric distribution of rotational velocity near zero, although they exhibit a small retrograde motion of V_ϕ ∼ −15 km s⁻¹ in the outer halo. These stars are considered to be stripped from a massive dwarf galaxy on a orbit of high inclination, based on their MDF with a peak of [Fe/H] = −1.3.

In addition, our analysis indicates that at least two low-mass progenitors are required to explain the distinct MDFs and dynamical properties between the LOI and HOI stars that are on retrograde motion with highly eccentric and tangential orbits in the OHR.

We thank an anonymous referee for a careful review of this paper, which has improved the clarity of its presentation. Y.S.L. acknowledges support from the National Research Foundation (NRF) of Korea grant funded by the Ministry of Science and ICT (NRF-2018R1A2B6003961 and NRF-2021R1A2C1008679). T.C.B. acknowledges partial support for this work from grant PHY 14-30152; Physics Frontier Center/JINA Center for the Evolution of the Elements (JINA-CEE), awarded by the U.S. National Science Foundation.

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

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