A Swing of the Pendulum: The Chemodynamics of the Local Stellar Halo Indicate Contributions from Several Radial Merger Events

It is currently our understanding that the Milky Way (MW) halo was primarily built up by accretion events (Helmi 2020). Streams of stars stripped from satellites that are currently in the process of falling into the MW (Newberg & Carlin 2016) provide evidence of this accretion history. Within the last several years, it has become apparent that the MW halo contains merger debris from at least one major merger event that is no longer identifiable in stellar density alone, but also requires velocity information of stars and/or compact objects in the inner and outer MW halo (Belokurov et al. 2018; Helmi et al. 2018; Donlon et al. 2019, 2020; Kruijssen et al. 2020; Naidu et al. 2020; Horta et al. 2021; Naidu et al. 2021; Malhan et al. 2022).

The primary motivation for study of MW halo substructure is cosmological: the present-day configuration of the MW allows us to test our theories for galaxy formation and the early Universe. One of the first models that attempted to explain how the MW formed was the monolithic collapse model (Eggen et al. 1962), which proposed that large galaxies formed out of freefall collapse of massive, uniform gas clouds.

The figurative "pendulum" of public thought swung the other way when Searle & Zinn (1978) proposed a competing hierarchical model of halo formation, as the generally accepted standard model of Λ-CDM cosmology (Dekel & Silk 1986; Frenk et al. 1988; Navarro et al. 1996) favors a stratified collapse model (White & Rees 1978; White & Frenk 1991). In this model, many small pockets of star formation form out of the collapse of gas clouds. These clusters are gravitationally attracted to one another, and the accretion of lower-mass objects onto larger halos causes the build-up of progressively more massive objects over time. The popular opinion of how the MW's halo was formed transitioned away from a single collapse and toward a picture of continuous accretion.

In order to provide constraints for cosmological models, it is useful to identify these accretion events within our Galaxy. Numerical models of halo formation were examined by Helmi & White (1999), who concluded that a disrupting satellite would leave "few obvious asymmetries [...] in the distribution of particles in configuration space," and that identification of halo substructure would require three-dimensional velocity measurements. The same year, Helmi et al. (1999) identified one such accretion event in halo stars, based on their kinematics. They estimated that the local Solar region is populated by 300–500 low-mass stellar streams, which would cause the stellar halo to appear phase-mixed even if each individual accretion event was not actually phase-mixed.

Again, the pendulum swung back with the discovery that the Sagittarius spheroidal dwarf galaxy (Sgr dSph; Ibata et al. 1994) had a tidal stream that was apparent in the density of stars (Yanny et al. 2000; Ibata et al. 2001; Newberg et al. 2002; Majewski et al. 2003). With the realization that it was possible to detect stellar halo substructure using only the positions of stars, tidal streams became a popular way to study the halo. Presently, we know of 70+ stellar streams populating the Galactic halo (a compilation of known streams is available in the galstreams Python package; Mateu 2017), although new streams are constantly being identified and characterized (e.g., Ibata et al. 2021; Martin et al. 2022), and some streams that were previously thought to be independent may be combined into the same structure after further analysis (e.g., Li et al. 2021).

At the same time, it also became apparent that there were MW halo structures that did not appear to be stellar streams, but were nonetheless apparent in stellar density. These include the Virgo Overdensity (VOD; Vivas et al. 2001; Newberg et al. 2002; Jurić et al. 2008; Donlon et al. 2019) and the Hercules Aquila Cloud (Belokurov et al. 2007; Gryncewicz et al. 2021), which were both identified as overdensities in the locations of halo stars. These structures, along with the existence of the MW halo's stellar break at ∼25 kpc (Deason et al. 2011), were interpreted as either the apocenter pileup of a few accreted satellites, or a single major merger.

This "ancient last major merger" hypothesis became especially popular. Decades ago, it was proposed that the MW thick disk was the result of a major merger heating the MW's proto-disk (Quinn & Goodman 1986; Velazquez & White 1999). This merger event would have to be sizeable, compared to the proto-MW, in order to induce a sufficiently large change in the velocity dispersion of the proto-disk stars (e.g., Grand et al. 2020).

The large age of thick disk stars meant that the merger event had to have been accreted early on in the MW's history (at least 8 Gyr ago), as the end of thick disk star formation would mark the time of the merger event. It was also argued by Deason et al. (2013) that the stellar break in the MW halo could be explained by a major accretion event 6–9 Gyr ago, and Deason et al. (2015) found that the ratio of blue straggler stars to blue horizontal branch stars in the MW halo favored an accretion history consisting of a few more-massive dwarfs rather than many smaller objects. This argument was followed up by Deason et al. (2018), who concluded that the metallicity and location of the stellar break could only be explained by a scenario involving the accretion of a single massive dwarf progenitor. The stage had been set; all that remained was to locate this "last major merger."

With the advent of the Gaia survey (Gaia Collaboration et al. 2016), parallaxes and proper motions of billions of stars within several kpc of the Sun became available, which made it possible to perform large-scale kinematic and dynamical analysis on the local Solar region. Multiple groups rapidly identified the ancient last major merger of the MW. Belokurov et al. (2018) found that the local stellar halo stars predominantly had orbits with high radial anisotropy, and claimed that these very radial stars composed the debris of the ancient last major merger; they named this merger the "Gaia Sausage." Simultaneously, Helmi et al. (2018) identified a population of local halo stars based on their kinematics and chemical abundances, and claimed that they composed the debris of the "ancient last major merger," which they named the "Gaia-Enceladus" structure. This structure is typically referred to as the Gaia Sausage/Enceladus (GSE) structure.

It was natural at this time to associate these halo stars with the ancient last major merger concept that had occupied the literature for several decades. Because the ancient last major merger was also thought to dominate the stellar halo, it was also reasonable to associate the VOD and HAC with this merger event (Simion et al. 2019). Notably, Naidu et al. (2021) were able to recover the general shape of the stellar halo with a simulation of a very massive (2.5:1) disky progenitor that fell into the MW on a retrograde orbit around 8 Gyr ago.

In addition to the dynamical arguments for the existence of the GSE, numerous independent studies of spectroscopic data (Helmi et al. 2018; Vincenzo et al. 2019; Feuillet et al. 2020; Naidu et al. 2020; Hasselquist et al. 2021; Buder et al. 2022) found that the metallicity distribution function (MDF) for GSE stars is consistent with the expected MDF of a massive dwarf galaxy (Kirby et al. 2013). Namely, the MDF of GSE stars has a single peak near [Fe/H] ∼ −1.2, with an extended metal-poor tail.

However, some evidence argues against the ancient major merger scenario. It was found that the formation of a thick disk does not require a major merger event; if the MW's proto-disk were lumpy rather than uniform, it would be possible to generate a disk system that is bimodal in both chemical abundances and velocity dispersion (Clarke et al. 2019; Amarante et al. 2020; Beraldo e Silva et al. 2020). Further, Beraldo e Silva et al. (2021) showed that the thin and thick disks likely formed at the same time early on in the MW's history. These studies call into question the requirement that the radial stars in the local Solar region must have been accreted at early times in order to puff up the thick disk.

The inner halo structure known as the "Splash" (Belokurov et al. 2020) initially supported an ancient last major merger scenario, because it was thought to originate from heating of the inner MW proto-disk stars to orbits with large vertical heights, probably with the same merger that heated the thick disk. However, this structure also can be explained if the early MW contained a lumpy proto-disk (Amarante et al. 2020), or instead could have formed within gas outflows that were thrown out from the disk and then fell back into the inner MW (Yu et al. 2020).

The morphology of these radial stars in phase space also called the large age of the structure into question. Donlon et al. (2019) found that the VOD could be recreated using an N-body simulation of a dwarf galaxy on a radial orbit that began disrupting 2 Gyr ago, but they were not able to create overdensities in the halo with long integration times. They named the merger event that was responsible for creating the VOD the "Virgo Radial Merger" (VRM). Following those results, Donlon et al. (2020) found that the VOD and HAC regions contained shell structure, which is an indicator of a radial merger event. They found that, if the progenitor of the VOD and HAC had been disrupted more than 5 Gyr ago, then we would not expect to observe any shell structure in the MW. Based on this phase-mixing argument, as well as reverse orbit integration of the shell structures, Donlon et al. (2020) argued that the progenitor of the VOD and HAC collided with the Galactic center within the last 3 Gyr (although Donlon et al. 2020 used a spherical halo potential for their orbit integrations, which might artificially decrease the amount of time for which large-scale asymmetries remain visible, based on the findings of Han et al. 2022b). Additionally, Grunblatt et al. (2021) computed asteroseismic ages for several GSE stars, and found that some stars had been formed as recently as 6 Gyr ago; this called the proposed infall time of the GSE progenitor into question, but allowed for a more recent merger time.

The recent time of the VRM appeared to be in tension with other published results. First, the thick disk could not have been heated 3 Gyr ago. That is okay if the thick disk was heated by another mechanism; the more recent merger would interact with a more massive MW, and should not cause a recent disk heating event (see the mass/velocity scales of Grand et al. 2020). The VRM could instead have produced the observed increases in star formation rate in the MW disk around 2.5–2.7 Gyr ago (Snaith et al. 2014; Mor et al. 2019), although these starbursts could instead be caused by the passage of the Sagittarius Dwarf Galaxy (Sgr dSph; Ruiz-Lara et al. 2020). One might worry about a recent merger time, given that the GSE/VRM stars are predominantly old; however, the peak ages of these stars are not necessarily tied to the time that their progenitor was accreted. The best example of this is that the Sgr dSph stars appear to have been formed before 8 Gyr ago (de Boer et al. 2015; Hasselquist et al. 2021), even though the dwarf is still bound today. For a detailed discussion of these points, we direct the reader to Section 7 of Donlon et al. (2020).

Since the initial discoveries of the ancient last major merger, there have been numerous works that propose the existence of different phase-mixed merger events within the MW halo. These include the strongly retrograde "Thamnos" structure (Koppelman et al. 2019), the "Kraken" structure that was predicted via populations of globular clusters (Kruijssen et al. 2020), the "Inner Galaxy Structure" (also called "Heracles") that was identified in halo stars (Horta et al. 2021), multiple retrograde structures identified in halo stars, (Naidu et al. 2020), and the "Pontus" structure identified in integral clustering of MW globular clusters and dwarf galaxies (Malhan et al. 2022). These different halo structures all make up the phase-mixed stellar halo, and it is not likely that they all come from a single massive merger event. However, none of these structures are thought to dominate the stellar halo near the location of the Sun.

Of particular importance to this work are the halo structures identified in the Solar neighborhood by Donlon et al. (2022). While previous work called into question the arguments for an ancient major merger or its association with halo structures like the Virgo Overdensity, Donlon et al. (2022) found multiple radial structures in a volume of halo stars that should have been dominated by only the GSE.

Using local dwarf stars identified by Kim & Lépine (2022) combined with proper motions from Gaia EDR3, Donlon et al. (2022) identified VRM stars using the criteria that they are relatively metal-rich, are slightly counter-rotating from the disk, and have large radial velocities as viewed from the Galactic center (v_R ). In addition, Donlon et al. (2022) identified two new structures in these dwarf stars. The first, named "Cronus," was identified as a relatively metal-rich structure with low V_R , and is slightly corotating with the MW disk. The second, named "Nereus," was identified as a metal-poor structure with no net rotational velocity, and has v_R somewhere between those of Cronus and the VRM. Notably, all of these structures are made up of stars that would typically be classified as belonging to the GSE structure. Donlon et al. (2022) conclude that the combination of the VRM, Nereus, and Cronus structures cannot be explained by a single ancient, massive merger event, and that the structure commonly referred to as the GSE is really a combination of all of these components.

In this work, we classify local Solar neighborhood halo stars using their chemical abundances as well as their dynamical properties. We find that these local halo stars cluster into four components of the stellar halo, which correspond to the VRM, Nereus, and Cronus components of the halo that were identified in Donlon et al. (2022), as well as Thamnos. We again conclude that the inner halo stars that are commonly attributed to the GSE merger event actually belong to a combination of these four halo components. Our results suggest the VRM was only recently accreted by the MW, Cronus was accreted early on in the MW's history, and Nereus and Thamnos are made up of multiple smaller radial merger events that may have a range of accretion times.

We believe that we are experiencing yet another shift in our understanding of the MW's inner stellar halo. At first, it was thought that galaxies formed under monolithic collapse, until hierarchical collapse was shown to describe the Universe. In the past, the literature was certain that stellar streams would only be found kinematically—until a stream was discovered using the positions of stars. The popular understanding of the MW's inner stellar halo has transitioned from a belief that the local stellar halo was built up from 300 to 500 structures, to a belief that the majority of the MW's inner stellar halo was deposited by a single merger event. Now, with substantial evidence that the local stellar halo was actually built up from several merger events, we find that the pendulum is once again swinging, and that our understanding of the stellar halo must now allow for several mergers contributing to the accreted halo stars in the local Solar neighborhood.

In this section, we describe the process of obtaining and cleaning the data that we analyze in this work.

We constructed a data set consisting of stars from SDSS Data Release 17 (Blanton et al. 2017; Abdurro'uf et al. 2022) that were in the APOGEE survey (Majewski et al. 2017), and stars from GALAH DR3 (Kos et al. 2017; Buder et al. 2021). The data from these surveys were supplemented by Gaia EDR3 data (Gaia Collaboration et al. 2016, 2021).

For both of the APOGEE and GALAH samples, the six-dimensional phase-space information for each star was derived using astrometric positions, distances obtained from Gaia EDR3 parallaxes, Gaia EDR3 proper motions, and line-of-sight velocities from the corresponding spectroscopic survey.

Stars from both catalogs were required to satisfy the Gaia EDR3 astrometric cuts astrometric_excess_noise <0.25, which removes stars that contain disagreements between observations of a source and the Gaia astrometric model (Lindegren et al. 2012), and phot_bp_rp_excess_factor <1.5, which removes stars that have contamination issues in the Gaia BP and RP photometry (Riello et al. 2018).

The Gaia EDR3 parallaxes were adjusted using the gaiadr3_zero-point Python package (Lindegren et al. 2021). Additionally, it was required that the parallax for each star be positive before zero-point offset correction, and that the relative error in the parallax was less than 15%. We restricted our data to be within 5 kpc of the Sun, as distances derived from Gaia parallaxes beyond this range can become unreliable.

This work uses a right-handed Galactocentric Cartesian coordinate system. The Sun is located at (−8.21, 0, 0.25) kpc (consistent with Jurić et al. 2008; Gravity Collaboration et al. 2019), where positive X points from the Sun toward the Galactic center, positive Y points in the direction of the Sun's motion, and $\hat{Z}=\hat{X}\times \hat{Y}.$ The circular velocity at the Sun's position is 233.1 km s⁻¹, and the Sun's motion with respect to the LSR is (11.1, 15.17, 7.25) km s⁻¹ (Schönrich et al. 2010, except v_y is slightly larger than their reported value; see Buder et al. 2022 for details). This kinematic frame was used for the GALAH DR3 value-added dynamical quantity catalog GALAH_DR3_VAC_dynamics_v2, and we choose to use it for the APOGEE data as well for consistency. However, it should be noted that this choice of frame can have a substantial effect on the measured dynamical quantities of stars, particularly L_z (see Appendix A). We take the sign of the azimuthal velocity v_ϕ to be positive in the direction of disk motion. Note, however, that the sign of L_z is negative for prograde motion.

The GALAH_DR3_VAC_dynamics_v2 catalog includes several orbital parameters, including the maximum orbital height above the plane ${z}_{\max }$, maximum (r_apo) and minimum (r_peri ) distances from the Galactic center, and orbital eccentricity e. It also includes energies and actions that are computed in the axisymmetric McMillan2017 potential (McMillan 2017) using the galpy package (Bovy 2015). These actions were computed using a Staekel fudge algorithm with a focus of 0.45. We took these dynamical values directly from the GALAH_DR3_VAC_dynamics_v2 catalog for the GALAH stars, except we negated the GALAH values of L_z (=J_ϕ ) to match the definition of our coordinate system. Dynamical quantities were calculated for the APOGEE stars using the actionAngleStaeckel class in galpy, with the same settings as the GALAH DR3 dynamical quantity calculations (Buder et al. 2022).

We queried all available chemical abundance data for each star, which included 23 abundances for APOGEE stars and 35 abundances for GALAH stars. However, for some abundances there were only data for a few or no stars in our final data set, so they could not be included in our analysis. For APOGEE stars, the [P/Fe] and [Cu/Fe] abundances were not included (for a total of 21 APOGEE usable abundances); for GALAH stars, the [C/Fe], [ScII/Fe], [CrII/Fe], [Sr/Fe], [Rb/Fe], [Mo/Fe], and [Ru/Fe] abundances were not included (for a total of 28 usable GALAH abundances).

Of particular interest are stellar [Mg/Mn] and [Na/Fe] abundances. [Mg/Mn] is often preferred to [Mg/Fe] in other works utilizing GALAH abundances (e.g., Hawkins et al. 2015; Das et al. 2020; Buder et al. 2022), due to the fact that Mg and Mn probe orthogonal nucleosynthesis processes (Ting et al. 2012). Buder et al. (2022) provide detailed reasoning for using [Mg/Mn] and [Na/Fe] in their analysis, and they found success with using those abundances to identify accreted halo stars.

In order to ensure that our data set did not contain stars that are members of compact structures, we removed all stars that were located within 03 of any cluster that is listed in the Catalog of Parameters for Milky Way Globular Clusters (Harris 1996) and is located within 10 kpc of the Sun. In the case of the Omega Centauri globular cluster, all stars within 05 of the center of the cluster were removed instead, due to the large apparent size of the object.

Finally, in order to remove the majority of disk stars from our data set, we required that all stars had velocities that differed by at least 180 km s⁻¹ from the LSR. After this cut, there were still many thick disk stars present in the data, which were removed via chemical abundance cuts (Section 3).

Additional quality cuts were made for both APOGEE stars and GALAH stars, which are outlined in the following subsections. Figure 1 shows the spatial distribution of the cleaned data set. Stars from the two surveys are located in different volumes, due to the on-sky and distance distributions of stars in each survey.

Zoom In Zoom Out Reset image size

2.1. APOGEE Data Selection and Cleaning

2.2. GALAH Data Selection and Cleaning

This section contains our procedure for selecting thick disk and Splash stars in our data set based on their kinematics and chemical abundances, as well as a brief discussion about the properties of these stars in our data compared to other literature. This selection of thick disk & Splash stars will be removed from our sample in order to analyze the properties of the accreted stellar halo.

The kinematic cuts made in the previous section removed the majority of thin disk stars, which have velocities very similar to that of the LSR. However, thick disk stars can have velocities that are different enough from the LSR that they are not removed by our kinematic cut. Additionally, Belokurov et al. (2020) showed that there was a population of metal-rich inner halo stars that appears to smoothly extend stars with thick disk [Fe/H] values down to somewhat retrograde v_ϕ values, which they called the "Splash"; these stars will similarly not be removed from our velocity cut.

3.1. Selecting Thick Disk/Splash Stars

Das et al. (2020) showed that, in APOGEE data, disk stars appeared to separate from halo stars in plots of [Al/Fe] versus [Fe/H]. We plot our APOGEE data in [Al/Fe] versus [Fe/H] in the top row of Figure 2, except we color each star by its velocity with respect to the LSR:

$$\begin{eqnarray}&&{V}_{\mathrm{LSR}}=\sqrt{{v}_{R}^{2}+{({v}_{\phi }-233.1)}^{2}+{v}_{z}^{2}}.\end{eqnarray} \tag{ 1 }$$

Thick disk/Splash stars should have smaller values of V_LSR than halo stars that are not associated with the disk. There is a clear clump of stars with small V_LSR at ([Fe/H], [Al/Fe]) = (−0.6, 0.3) in the APOGEE data. Additionally, there are two trails of stars leading from the low-V_LSR clump down and to the left; these are high-α disk stars (Belokurov & Kravtsov 2022), which we also want to remove from our data set. These in situ stars were selected using the criteria

$$\begin{eqnarray}&&\begin{array}{l}[\mathrm{Al}/\mathrm{Fe}]\gt -0.5\ [\mathrm{Fe}/{\rm{H}}]-0.6,\mathrm{and}\\ [\mathrm{Al}/\mathrm{Fe}]\gt \ 0.\end{array}\end{eqnarray} \tag{ 2 }$$

This cut is shown with dashed black lines in the top left panel of Figure 2. This cut classified 1686 stars (73.3%) of the cleaned APOGEE sample as "in situ" stars, and the remaining 615 (26.7%) APOGEE stars were classified as accreted halo stars.

Zoom In Zoom Out Reset image size

Das et al. (2020) also used [Mg/Mn], [C/Fe], and [N/Fe] to identify stars that belonged to the thick and thin disk. We chose to avoid requiring our APOGEE stars to have low-uncertainty [Mn/Fe] and [N/Fe] information, because using those abundances greatly restricts the number of viable halo stars. While most of our APOGEE halo stars have high-quality [C/Fe] abundances, we were not able to separate disk and halo stars using these [C/Fe] abundances.

Chemical cuts were also used in the GALAH stars to isolate thick disk/Splash stars. However, we chose to use [Na/Fe] rather than [Al/Fe] for this task because there are 5828 [Na/Fe] abundances in our data set, but only 3937 GALAH [Al/Fe] abundances. Buder et al. (2022) also chose to use [Na/Fe] to determine candidate GSE stars in GALAH data instead of [Al/Fe], for the same reason. There is a clump of stars with small V_LSR at ([Fe/H], [Na/Fe]) = (−0.6, 0.15) in the GALAH data. Our selection of thick disk/Splash stars in the GALAH data is

$$\begin{eqnarray}&&\begin{array}{l}[\mathrm{Na}/\mathrm{Fe}]\gt -0.6\ [\mathrm{Fe}/{\rm{H}}]-0.6,\mathrm{and}\\ [\mathrm{Fe}/{\rm{H}}]\gt -1.\end{array}\end{eqnarray} \tag{ 3 }$$

This cut is shown as dashed black lines in the top right panel of Figure 2. This cut classified 3,413 stars (58.6%) of the cleaned GALAH sample as "in situ" stars, and the remaining 2,415 (41.4%) cleaned GALAH stars were classified as accreted halo stars.

The ratio of "in situ" stars to accreted stars is substantially different for the APOGEE and GALAH samples. This may be due in part to our explicit removal of the α-rich disk stars in APOGEE data, which we were not able to do in the GALAH data. However, the main difference is likely the selection functions of the surveys: APOGEE is designed to obtain spectra of any red (primarily giant) stars, while GALAH contains a larger fraction of accreted stars because its targets were chosen with Galactic archeology in mind.

3.2. Properties of Thick Disk/Splash Stars

The middle row of Figure 2 shows v_⊥ versus v_ϕ for halo stars (green) and thick disk/Splash stars (purple), where

$$\begin{eqnarray}&&{v}_{\perp }=\sqrt{{v}_{R}^{2}+{v}_{z}^{2}}.\end{eqnarray} \tag{ 4 }$$

The thick disk/Splash stars preferentially hug the boundary of the velocity cut in the bottom right of each panel, where V_LSR is the smallest. However, there are some thick disk/Splash stars that appear kinematically similar to the halo, and have large v_⊥ and small or negative v_ϕ . Note that there are halo stars located at the same v_⊥ and v_ϕ as the thick disk/Splash stars. If one uses a more restrictive velocity cut instead of a chemical abundance cut to isolate halo stars, there will still be thick disk/Splash star contamination at low and negative v_ϕ , and the data will lose a number of halo stars with positive v_ϕ .

The bottom row of Figure 2 shows v_ϕ and [Fe/H] for halo stars and thick disk/Splash stars, as well as the original boundary drawn over Splash stars by Belokurov et al. (2020). Surprisingly, we find that the thick disk/Splash stars that we select with our chemical abundance cuts are shifted to lower [Fe/H] than was suggested by the Belokurov et al. (2020) selection box. Horta et al. (2021) also found Splash stars with [Fe/H] < −0.8, and Donlon et al. (2022) claimed that there may be thick disk/Splash dwarf stars with [Fe/H] values as low as −1.

Amarante et al. (2020) claimed that a population of Splash stars can be generated by disk star+gas dynamics alone, without the need for a major merger event. However, that work was not able to exactly match the observed Splash star data; Belokurov et al. (2020) placed the lower v_ϕ boundary for Splash stars at −150 km s⁻¹, but the Amarante et al. (2020) model struggled to create Splash stars with v_ϕ < −50 km s⁻¹. In our data, the observed thick disk/Splash star distribution appears to terminate near v_ϕ = −50 km s⁻¹, which matches the cutoff of the Amarante et al. (2020) model and is therefore consistent with a scenario in which the Splash star population was not created by a merger event.

Additionally, we do not see a bend in the v_ϕ -[Fe/H] distribution of thick disk/Splash stars above v_ϕ > −100 km s⁻¹; this is because our kinematic cuts removed the majority of disk stars, which dramatically outnumber Splash stars. Because Belokurov et al. (2020) drew their Splash star region over a row-normalized histogram, they were not able to identify Splash stars in the rows that contained disk stars. In our data, after the majority of disk stars have been removed, the thick disk/Splash star population actually extends up to v_ϕ > 200 km s⁻¹.

In this section, we explore the expected chemodynamic distributions that would be generated from a single massive, ancient merger event. We then compare this model with observed data and show that the energies and abundances of local halo stars are not consistent with this scenario. These observed chemodynamic trends in local halo stars can instead be explained by the local halo containing debris from multiple merger events.

4.1. Predictions for Chemodynamic Trends Due to a Single Massive, Ancient Merger Event

Based on the current literature surrounding the GSE, we can surmise two generally accepted ideas about its progenitor dwarf galaxy: (i) it was massive, and (ii) it was accreted early on in the MW's history (Belokurov et al. 2018; Helmi et al. 2018; Gallart et al. 2019; Bonaca et al. 2020; Kruijssen et al. 2020).

A massive progenitor informs us that the progenitor probably sustained enough star formation to generate a negative radial abundance gradient in its member stars. For example, Mercado et al. (2021) showed that, in FIRE simulations, dwarf galaxies that assembled their stars early enough to have created the GSE have strong negative radial metallicity gradients. While Maiolino & Mannucci (2019) noted that there have been observations of positive metallicity gradients in galaxies, these are only weakly positive and are probably transient phenomena due to recent star formation, and therefore they are less likely to be good representatives of the ancient GSE merger. A massive progenitor (compared to the Milky Way at the time of infall) will also lose orbital energy due to dynamical friction.

The second point tells us that, regardless of the progenitor's initial state, we should expect its debris to be well-mixed at the present day. With these assumptions in mind, we can build a qualitative model for the chemodynamic distributions of a single massive, ancient merger's member stars at the present day.

Figure 3 shows two different scenarios for the accretion of a massive, ancient merger onto an MW-like galaxy. The two scenarios differ based on whether the outer regions of the progenitor become unbound before or after the progenitor dwarf galaxy becomes sufficiently radialized and collides with the host galaxy.

Zoom In Zoom Out Reset image size

In the first scenario, the outer portion of the progenitor with small [Fe/H] and large [α/Fe] is stripped before the bound core of the progenitor collides with the host galaxy. Once the outer portions are stripped from the progenitor, they can no longer lose energy and angular momentum via dynamical friction, so these quantities are locked in for the low-[Fe/H] stars. The bound progenitor core will continue falling in until it collides with the host galaxy, and will progressively lose more energy and angular momentum. At the present day, we expect to see a radial abundance gradient in the dwarf member stars, which are now located in the host galaxy. The low-[Fe/H] stars from the progenitor will have nonzero angular momentum and higher energy than the stars from the progenitor with large [Fe/H], which will have low energy and little-to-no angular momentum.

In the second scenario, the outer portion of the progenitor remains bound until the progenitor collides with the host galaxy. In this case, both the low- and high-[Fe/H] stars in the progenitor have the same angular momentum and energy on average. However, the low-[Fe/H] portion of the progenitor has a larger internal velocity dispersion and larger range of positions than the high-[Fe/H] portion of the progenitor. As a result, after the collision, the low-[Fe/H] stars will have a larger range of energy and angular momentum values with respect to the host galaxy. Additionally, the larger velocity dispersion of the low-[Fe/H] stars gives them somewhat higher energies after accretion.

In reality, the infall of a massive dwarf galaxy will be described by some mixture of these two scenarios. Some material may be stripped from the progenitor early on, and the progenitor may still have an abundance gradient at the time of impact.

If the majority of the stars in the local stellar halo were deposited by a single massive, ancient merger event, those stars should follow the chemodynamic trends shown in Figure 3. However, this is not the case; in the following subsections, we show that the observed chemodynamic trends in the local stellar halo are not consistent with either of these scenarios, or a mixture of scenarios.

4.2. Energy and Abundance Trends

Figure 4 shows the distributions of [Fe/H] and [α/Fe] abundances as functions of each star's orbital energy for the accreted APOGEE and GALAH stars that have orbital eccentricity > 0.75. This is a reasonable criterion for selecting GSE member stars (e.g., Naidu et al. 2020); thus, this figure only contains stars which are likely to belong to the GSE. As one looks at stars with progressively higher energies, the average [Fe/H] of those stars increases. Similarly, the stars with higher energy have lower average [α/Fe] values.

Zoom In Zoom Out Reset image size

In the rightmost column of Figure 4, we show the predicted energy-[Fe/H] and energy-[α/Fe] trends from Figure 3 in gray. On top of this expected trend, we show the sigma-clipped mean of the APOGEE and GALAH data as red and blue dashed lines, respectively. The observed data do not follow the predicted trend for [Fe/H] or [α/Fe]. While our model predicts a trend with a negative slope in [Fe/H] versus energy, both the APOGEE and GALAH data show trends with positive slopes. Similarly, our model predicts a positively sloped trend in [α/Fe] versus energy, but the APOGEE and GALAH data show negatively sloping trends.

Similar energy-abundance trends are also present in other abundances than [Fe/H] and [α/Fe]. Figure 5 shows heatmaps of orbital energy as a function of various chemical abundances in the APOGEE and GALAH data. Overall, stars with larger values of [X/Fe] abundances have lower energy than stars with low [X/Fe] abundances, which is the opposite of the expected trend based on our single merger model. Other works have previously found a connection between the chemical abundances of halo stars and their dynamics; our results are similar to those of Das et al. (2020) and Ness et al. (2022), who also used [Mg/Fe] and [Al/Fe] in order to isolate distinct halo structures in APOGEE data, as well as those of Buder et al. (2022), who used [Na/Fe] and [Mg/Mn] in their analysis of GALAH stars.

Zoom In Zoom Out Reset image size

We conclude that the energy/abundance trends in the observed data are inconsistent with the predicted chemodynamic trends of the GSE; the observed data can only be explained if the progenitor galaxy had a strong positive metallicity gradient, which is unlikely for an early, massive merger. However, it is possible to explain the observed chemodynamic trends if the local solar region contains debris from several distinct radial merger events. Dwarf galaxies have a range of chemical abundances, and when they are accreted, the energy of their constituent material will depend on many factors such as the mass, orbital parameters, and infall time of the progenitor dwarf galaxy. If, for example, a low-[Fe/H], high-[α/Fe] dwarf galaxy fell in at early times (so that its material would have low-energy at the present day), and a high-[Fe/H], low-[α/Fe] dwarf galaxy fell in recently (so that its material would have high energy), then that could explain the observed energy-abundance trends. Notably, this scenario is consistent with the MW accretion history proposed by Donlon et al. (2022).

4.3.  L_z -Abundance Trends

The chemodynamic trends observed in Section 4 suggest that it should be possible to separate different halo structures in the solar neighborhood using multidimensional combinations of dynamical quantities and chemical abundances. In this section, we investigate this idea further, using the Gaussian Mixture Model (GMM) implementation from the scikit-learn python package (Pedregosa et al. 2011) to split the data into multidimensional Gaussian components. We then analyze the identified GMM components, and match them to known halo substructure based on their chemodynamic properties.

GMMs can struggle to correctly identify the correct number and shapes of components in a sample if they are fitting quantities that are correlated with each other. This is because, if two quantities are individually normally distributed but share some correlation, then their joint distribution is no longer expected to be a multidimensional normal distribution. Notably, this is true for orbital energy E and L_z , which are both functions of v_ϕ . In order to remove the correlation between the quantities in our fitting method, we chose to replace orbital energy with a "pseudo-energy" $\widetilde{E}$, defined as

$$\begin{eqnarray}&&\widetilde{E}=\displaystyle \frac{{{\boldsymbol{v}}}^{2}-{v}_{\phi }^{2}}{2}+{\rm{\Phi }}({\boldsymbol{x}}).\end{eqnarray} \tag{ 5 }$$

Provided that v_ϕ is small compared to the total velocity of each star, this quantity will be roughly an integral of the motion.

We took the quantities $\widetilde{E}$, L_z , [Fe/H], [α/Fe], [Na/Fe], and [Al/Fe] for each star in the APOGEE sample, and then optimized a single-component GMM over this data. We then calculated the Bayesian Information Criteria (BIC; Schwarz 1978) for this fit using the built-in scikit-learn BIC function for GMMs.

Once we had the BIC for the single-component GMM, we fit a two-component GMM on the same data and recalculated the BIC for the two-component fit. The BIC is a heuristic that is designed to evaluate the trade-off between the improvement in the likelihood score of a model fit with more parameters and the number of parameters that were added to that new model. A decrease in the BIC indicates a better model. We repeated this process for up to seven GMM components, and then compared the BIC of each model. The BIC was minimized by the fit with two components, so we adopt the two-component model as our best fit. The results of this fit are shown in Figure 6.

Zoom In Zoom Out Reset image size

We then repeated this procedure for the GALAH sample, except that [Al/Fe] was replaced by [Mg/Mn], because fewer GALAH stars have high-quality [α/Fe] measurements. The GMM model with four components had the lowest BIC for the GALAH data. The results of the four-component GMM fit for the GALAH data are shown in Figure 7.

Zoom In Zoom Out Reset image size

Each star in the APOGEE and GALAH samples was then assigned to one of the fit components in the corresponding GMM fit. The probability p_k that a star with parameters s belongs to the kth Gaussian component is defined as

$$\begin{eqnarray}&&{p}_{k}=\displaystyle \frac{{f}_{k}({\boldsymbol{s}})}{\sum _{i}^{4}{f}_{i}({\boldsymbol{s}})},\end{eqnarray} \tag{ 6 }$$

where f_k ( s ) is the evaluation of the multivariate normal distribution of the kth component at s .

This will separate the stars into the components from which they are actually associated only if the components are not overlapping in the six dimensions that we use to separate them. If the components are overlapping, then we only recover the underlying distributions, not the actual associations for each star. Here, we only need the distribution of each component in our six chemodynamical parameters to identify the properties of each component and the structure(s) to which each component belongs.

To explore this idea further, we assigned a "confidence" value to each star, defined as the maximum p_k that a given star had for all components. This confidence value can be interpreted as the likelihood that each given star was assigned to the correct GMM component, assuming that the GMM components are known exactly. For the APOGEE data, the mean confidence value is 90.3% and the median confidence value is 96.9%. For the GALAH data, the mean confidence value is 79.3% and the median confidence value is 83.3%. We tried removing all stars with confidence values below an arbitrary threshold value (75% and 90%) from the data when performing our analysis, but the overall high confidence values for the stars in the model meant that removing these stars did not substantially change any results.

Appendix B provides an assessment of the algorithm's ability to identify individual dwarf galaxies from a mixed distribution, based on a fit of known MW dwarf galaxy data. The algorithm is able to accurately separate the mixed data into the constituent dwarf galaxies, so we are confident that the algorithm is able to identify distinct dwarf galaxy components in the nearby stellar halo data as well.

5.1. Identifying the Separated Components

5.2. Differences between the APOGEE and GALAH Fits

Our GMM fits only utilize two kinematic and four chemical quantities for each star. However, our data contain many more dynamical quantities and chemical abundances than those used in the fitting procedure. In this section, we look at additional properties of the four components of the stellar halo.

The means and standard deviations of each dynamical quantity and chemical abundance for the identified halo structures are given in Table 1.

Table 1. Means and Standard Deviations of Kinematic, Dynamical, and Chemical Properties of Stars Assigned to Each Merger Event, Separated by Spectroscopic Survey

	Virgo Radial Merger		Cronus	Nereus		Thamnos
Quantity	APOGEE	GALAH	GALAH	APOGEE	GALAH	GALAH
∣vR ∣ (km/s)	156 ± 89	199 ± 94	91 ± 61	136 ± 95	147 ± 97	126 ± 80
vϕ (km/s)	7 ± 43	4 ± 44	54 ± 63	9 ± 91	24 ± 90	−55 ± 80
vz (km/s)	3 ± 92	−8 ± 85	7 ± 86	1 ± 93	−2 ± 84	5 ± 76
Eccentricity	0.73 ± 0.20	0.91 ± 0.09	0.74 ± 0.18	0.69 ± 0.23	0.79 ± 0.18	0.74 ± 0.20
${z}_{\max }$ (kpc)	7.0 ± 4.8	5.3 ± 5.7	3.1 ± 1.7	7.5 ± 13.4	4.4 ± 4.6	3.6 ± 2.5
rm (kpc)	7.6 ± 2.8	7.4 ± 3.7	4.5 ± 1.5	8.3 ± 7.6	6.8 ± 4.4	5.9 ± 2.1
rapo (kpc)	13.1 ± 4.7	14.3 ± 7.2	7.7 ± 2.2	14.1 ± 14.4	12.2 ± 8.2	10.3 ± 3.7
rperi (kpc)	2.2 ± 1.9	0.6 ± 0.7	1.3 ± 1.2	2.6 ± 2.3	1.4 ± 1.4	1.5 ± 1.3
Energy (105 km2 s−2)	−1.56 ± 0.16	−1.57 ± 0.20	−1.82 ± 0.14	−1.57 ± 0.23	−1.66 ± 0.24	−1.69 ± 0.14
Lz (kpc km s−1)	−64 ± 368	−24 ± 298	−347 ± 383	−64 ± 744	−146 ± 610	371 ± 528
JR (kpc km s−1)	631 ± 368	1023 ± 583	358 ± 157	677 ± 970	715 ± 630	547 ± 319
Jz (kpc km s−1)	391 ± 440	174 ± 219	174 ± 152	350 ± 416	174 ± 209	153 ± 170
[Fe/H]	−1.14 ± 0.19	−1.01 ± 0.21	−1.24 ± 0.25	−1.54 ± 0.25	−1.43 ± 0.36	−1.34 ± 0.23
[α/Fe]	0.18 ± 0.06	0.14 ± 0.08	0.24 ± 0.06	0.24 ± 0.07	0.27 ± 0.15	0.24 ± 0.10
[Li/Fe]	⋯	1.40 ± 0.71	1.37 ± 0.48	⋯	1.99 ± 0.67	1.61 ± 0.62
[C/Fe]	−0.34 ± 0.16	⋯	⋯	−0.28 ± 0.27	⋯	⋯
[CI/Fe]	−0.20 ± 0.20	⋯	⋯	−0.10 ± 0.26	⋯	⋯
[N/Fe]	0.14 ± 0.13	⋯	⋯	0.24 ± 0.26	⋯	⋯
[O/Fe]	0.31 ± 0.12	0.50 ± 0.22	0.62 ± 0.22	0.37 ± 0.12	0.67 ± 0.32	0.64 ± 0.20
[Na/Fe]	0.05 ± 0.35	−0.21 ± 0.12	−0.12 ± 0.09	0.05 ± 0.59	−0.01 ± 0.12	−0.15 ± 0.14
[Mg/Fe]	0.17 ± 0.08	0.10 ± 0.10	0.23 ± 0.09	0.25 ± 0.08	0.15 ± 0.16	0.20 ± 0.11
[Al/Fe]	−0.21 ± 0.06	−0.06 ± 0.21	0.13 ± 0.16	−0.24 ± 0.12	0.32 ± 0.38	0.20 ± 0.24
[Si/Fe]	0.19 ± 0.06	0.12 ± 0.11	0.24 ± 0.11	0.23 ± 0.07	0.28 ± 0.19	0.25 ± 0.14
[S/Fe]	0.37 ± 0.14	⋯	⋯	0.38 ± 0.24	⋯	⋯
[K/Fe]	0.18 ± 0.16	0.09 ± 0.15	0.17 ± 0.16	0.21 ± 0.32	0.14 ± 0.20	0.14 ± 0.14
[Ca/Fe]	0.17 ± 0.08	0.20 ± 0.11	0.28 ± 0.09	0.19 ± 0.18	0.29 ± 0.15	0.27 ± 0.11
[Sc/Fe]	⋯	0.07 ± 0.10	0.10 ± 0.09	⋯	0.10 ± 0.13	0.10 ± 0.10
[Ti/Fe]	−0.06 ± 0.11	0.20 ± 0.16	0.26 ± 0.14	−0.05 ± 0.21	0.37 ± 0.29	0.28 ± 0.21
[TiII/Fe]	0.19 ± 0.23	0.30 ± 0.13	0.37 ± 0.12	0.12 ± 0.29	0.39 ± 0.17	0.36 ± 0.13
[V/Fe]	−0.03 ± 0.37	0.05 ± 0.33	−0.01 ± 0.26	0.10 ± 0.42	0.17 ± 0.56	0.00 ± 0.36
[Cr/Fe]	−0.20 ± 0.30	−0.14 ± 0.13	−0.15 ± 0.13	−0.10 ± 0.38	−0.06 ± 0.22	−0.13 ± 0.15
[Mn/Fe]	−0.34 ± 0.12	−0.33 ± 0.10	−0.37 ± 0.09	−0.34 ± 0.20	−0.30 ± 0.19	−0.37 ± 0.11
[Co/Fe]	−0.21 ± 0.34	0.07 ± 0.45	0.20 ± 0.45	−0.10 ± 0.48	0.50 ± 0.75	0.24 ± 0.54
[Ni/Fe]	−0.05 ± 0.12	−0.16 ± 0.12	−0.12 ± 0.13	−0.06 ± 0.11	−0.09 ± 0.25	−0.13 ± 0.15
[Cu/Fe]	⋯	−0.49 ± 0.14	−0.40 ± 0.16	⋯	−0.41 ± 0.29	−0.48 ± 0.18
[Zn/Fe]	⋯	0.16 ± 0.18	0.23 ± 0.18	⋯	0.23 ± 0.19	0.22 ± 0.18
[Y/Fe]	⋯	0.12 ± 0.27	0.25 ± 0.30	⋯	0.12 ± 0.33	0.18 ± 0.30
[Zr/Fe]	⋯	0.44 ± 0.47	0.55 ± 0.38	⋯	0.67 ± 0.69	0.64 ± 0.47
[Ba/Fe]	⋯	0.31 ± 0.34	0.36 ± 0.37	⋯	0.25 ± 0.42	0.31 ± 0.36
[La/Fe]	⋯	0.40 ± 0.33	0.39 ± 0.31	⋯	0.59 ± 0.50	0.42 ± 0.32
[Ce/Fe]	−0.19 ± 0.28	−0.05 ± 0.32	−0.10 ± 0.32	−0.23 ± 0.39	0.13 ± 0.46	−0.04 ± 0.35
[Nd/Fe]	⋯	0.56 ± 0.22	0.48 ± 0.23	⋯	0.64 ± 0.35	0.53 ± 0.24
[Sm/Fe]	⋯	0.26 ± 0.31	0.19 ± 0.33	⋯	0.43 ± 0.49	0.22 ± 0.32
[Eu/Fe]	⋯	0.49 ± 0.15	0.40 ± 0.20	⋯	0.52 ± 0.27	0.47 ± 0.28

Note. Bold font indicates that the quantity was used in the GMM fit to identify the components.

Download table as:  ASCIITypeset image

6.1. Dynamics

In addition to the dynamical parameters that are available in the GALAH value-added catalogs, we introduce an additional dynamical parameter,

$$\begin{eqnarray}&&{r}_{m}=\displaystyle \frac{{r}_{\mathrm{apo}}+{r}_{\mathrm{peri}}}{2},\end{eqnarray} \tag{ 7 }$$

which has previously been used to characterize possible VRM stars by Zhao & Chen (2021). For structures that predominantly contain stars on radial orbits, this quantity can be more useful than either the apogalacticon or perigalacticon distances on their own; this is because apogalacticon distance is functionally identical to total orbital energy for radial orbits, and perigalacticon distance is usually set by the radial orbit's eccentricity. Additionally, r_m and ${z}_{\max }$ appear to contain information about orbital resonances with the Galactic bar (Moreno et al. 2015; Zhao & Chen 2021).

Figure 8 contains the distributions of v_z , orbital eccentricity, ${z}_{\max }$, and r_m for all of the GMM components. In the GALAH data, the VRM has larger orbital eccentricity, ${z}_{\max }$, and r_m values than the other components. Cronus has the lowest r_m values, which indicates that Cronus stars are typically located closer to the center of the Galaxy than the other three components. In the APOGEE data, the VRM has a slightly larger eccentricity, r_m , and ${z}_{\max }$ than Nereus, although the values are also consistent with being the same for both structures. In the GALAH data, Cronus also has a net positive v_z , while the VRM stars have a slightly negative v_z , and Nereus and Thamnos stars appear to have no net v_z overall. In the APOGEE data, the VRM and Nereus stars both have near zero overall v_z . In dynamical equilibrium, the v_z distribution of a population is expected to be symmetric around zero, so the offsets in v_z for the VRM and Cronus stars either indicate selection effects or that their stars are not in equilibrium.

Zoom In Zoom Out Reset image size

Figure 9 contains the energy and L_z angular momentum data for each sample, split up into the different GMM components. In the APOGEE sample, the VRM has no clear net rotation, but Nereus has a prograde L_z . The energies of the two components are roughly the same, although the VRM appears to have slightly higher energy than Nereus. In the GALAH sample, Cronus has a net prograde rotation, Nereus has a slight prograde rotation, the VRM has a slight retrograde rotation, and Thamnos has a strong retrograde rotation. Cronus has the lowest energy, followed by Thamnos, then Nereus, and the VRM has the highest energy. In both samples, the L_z and E of the local stellar halo components agree with the analysis of Donlon et al. (2022), provided that the overall L_z shift is taken into account (Appendix A).

Zoom In Zoom Out Reset image size

The Sequoia structure (Myeong et al. 2019) can be seen in Figure 9 at very retrograde L_z and high $\widetilde{E}$. It is not identified as a distinct halo component in either sample; in the APOGEE data, it is assigned to the Nereus component, and in the GALAH data it is assigned to a mixture of the Nereus and Thamnos components. As the Sequoia is a known halo substructure that our algorithm is not able to identify, it may be the case that there are other distinct structures that our procedure is not capable of separating from other components of the local stellar halo.

6.2. Chemical Abundances

The large number of chemical abundances available to us in the APOGEE and GALAH survey data allow us to obtain detailed chemical abundance distributions for each GMM component. Figure 10 contains histograms of each chemical abundance that was available in the APOGEE data, and Figure 11 contains histograms of all the available GALAH chemical abundances. These histograms are split up by component, which shows the differences in chemical composition of each local halo structure.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In both the APOGEE and GALAH data, the VRM is characterized by relatively low abundances of most elements, except for [Mn/Fe] (particularly in the APOGEE data) and possibly [Nd/Fe] and [Eu/Fe]. This could be due in part to the high [Fe/H] abundance of the VRM, which would appear to correspondingly lower the apparent amount of other elements with respect to Fe.

In general, Cronus stars have larger [X/Fe] chemical abundances than VRM stars. This is especially true for [O/Fe], [Na/Fe], [Mg/Fe], [Al/Fe], [Si/Fe], [Ca/Fe], [Ti/Fe], [TiII/Fe], [Co/Fe], [Cu/Fe], [Y/Fe], and [Zr/Fe] compared to VRM stars. As O, Mg, Si, and Ca are α elements, the high abundances of these elements in Cronus stars is reflected in Cronus stars' high [α/Fe] abundances.

Cronus' [Fe/H] content is similar to that of the larger classical dwarf galaxies such as the Fornax dSph or Leo I (Kirby et al. 2013). The high α-element abundances of Cronus stars indicate that the progenitor of Cronus terminated star formation before Type-Ia supernovae elements could enrich the gas within the progenitor. In contrast, stars in the Fornax dSph have roughly solar levels of [Mg/Fe] and [Si/Fe], which suggests that it sustained multiple epochs of bursty star formation (Hendricks et al. 2014). One way to explain Cronus stars having both high [Fe/H] and high [α/Fe] is if the Cronus progenitor was artificially quenched early on in its star formation history. It is plausible that this quenching could have been due to the Cronus progenitor being accreted by the MW.

In order to generate an [Fe/H] of ∼−1.5 while still maintaining a high [α/Fe] abundance, the progenitor of Cronus was likely more massive than the classical MW dSphs. However, there is a somewhat smaller number of Cronus stars in the local Solar region than VRM stars. This could be due to Cronus' location deeper in the MW's potential well; Cronus stars are preferentially distributed closer to the Galactic center (e.g., the bottom panel of Figure 8), whereas VRM stars are distributed throughout the inner halo. It could also be because the VRM debris not phased-mixed (Donlon et al. 2020), so the density of VRM stars in the local stellar halo could be higher than the average density of VRM stars in the MW.

As Mn is primarily produced in Type-Ia supernovae (Kobayashi et al. 2020), the large values of [Mn/Fe] in VRM stars indicate that the VRM population continued forming stars long enough to experience chemical enrichment via SNe Ia. When a dwarf galaxy is accreted, it is expected that its star formation is quenched; therefore, the VRM was probably only accreted recently, because we expect that it had a lengthy star formation history that was not quenched long ago. In contrast, the large [α/Fe] abundance of Cronus compared to its [Fe/H] abundance indicates that it may have been accreted earlier in the MW's history than the VRM, as it is plausible that its star formation was quenched when it was accreted. The [Fe/H] abundances of VRM stars indicate that the progenitor of the VRM was probably similar in size to the larger classical dSphs, such as Fornax and Leo I.

Cronus stars have higher abundances of [Y/Fe] and [Zr/Fe] and lower abundances of [Nd/Fe], [Sm/Fe], and [Eu/Fe] compared to VRM stars. [Y/Fe], and [Zr/Fe] are s-process elements, which are thought to be formed primarily in asymptotic giant branch stars (Busso et al. 1999). [Nd/Fe], [Sm/Fe], and [Eu/Fe] are r-process elements, which are primarily generated in neutron star mergers (Kobayashi et al. 2020). Interestingly, the VRM and Cronus stars appear to have similar abundance distributions of [Ba/Fe], [La/Fe], and [Ce/Fe], which are also neutron-capture elements. The differences in r-process and s-process element abundances in the VRM and Cronus stars likely contain information about the star formation histories of the two structures.

In this section, we use the chemical abundances and actions of stars in order to explore whether (i) radial action cuts produce an unbiased sample of the local stellar halo (specifically the stars that are typically associated with the GSE), and (ii) whether the components we identified in this work are chemically consistent with a single massive dwarf galaxy. We find that the typical radial action cuts used in the literature predominantly select VRM stars, and that the chemical evolution paths of our halo components indicate that their progenitors had different origins, although some of the components may be related to each other.

7.1. Radial Action Cuts and Local Stellar Halo Populations

In order to further explore the dynamics of our halo star data, we made use of the radial orbital action J_R . The radial action can be thought of as the magnitude of a star's orbital oscillation toward or away from the Galactic center; near the Sun, the radial action is essentially equal to the part of a star's energy that comes from its radial motion. As the potential energy of all stars in the local solar region is roughly similar, it follows that J_R is roughly proportional to ${v}_{R}^{2},$ or equivalently $\sqrt{{J}_{R}}\propto | {v}_{R}| .$

Within the local stellar halo, GSE debris is characterized as having a large spread of radial velocity and low spread in rotational velocity (Belokurov et al. 2018). As a result, it will have a large range of ∣v_R ∣ values, and therefore a large range of $\sqrt{{J}_{R}}$ values. Despite this, it is common to isolate GSE debris by only using stars with low L_z and $\sqrt{{J}_{R}}$ greater than some value, usually 30 (kpc km s⁻¹)^1/2 (Feuillet et al. 2020; Hasselquist et al. 2021; Buder et al. 2022).

The top panel of Figure 12 shows histograms of $\sqrt{{J}_{R}}$ for each of the GALAH GMM components. Cronus stars dominate the data with $\sqrt{{J}_{R}}\lt 25$ (kpc km s⁻¹)^1/2, and VRM stars dominate the data with larger $\sqrt{{J}_{R}}$ values. The "GSE" stars that are isolated by a cut of $\sqrt{{J}_{R}}\gt 30$ (kpc km s⁻¹)^1/2 are really VRM stars, and therefore they are not representative of the chemical and kinematic properties of the local stellar halo as a whole. The often-used cut at $\sqrt{{J}_{R}}$ = 30 (kpc km s⁻¹)^1/2 is higher than the observed transition from Cronus to VRM stars at $\sqrt{{J}_{R}}$ = 25 (kpc km s⁻¹)^1/2. We find it unlikely that using a larger $\sqrt{{J}_{R}}$ cut provides any meaningful advantage; the high cut unnecessarily removes VRM stars from the sample, and contamination from Nereus and Thamnos stars will occur even in data with $\sqrt{{J}_{R}}\gt 40$ (kpc km s⁻¹)^1/2.

Zoom In Zoom Out Reset image size

It may not be so surprising that a high-$\sqrt{{J}_{R}}$ selection produces a sample that is not representative of the entire local stellar halo. A high-$\sqrt{{J}_{R}}$ selection will only select stars that are in the double-lobed velocity structure in the local stellar halo (Donlon et al. 2022), rather than the entire "sausage" velocity structure. This selection ignores stars that are located near (v_r , v_ϕ ) = (0, 0), which make up a substantial fraction of halo stars. If only the stars in the double-lobed velocity structure are considered to be GSE stars, then there must be some other origin for the stars with low ∣v_r ∣, and the GSE stars cannot make up the majority of the local stellar halo.

7.2. Using Radial Action to Isolate Chemical Evolution Paths

While populations of stars are often expected to clump together in chemical abundances, in reality stars from a forming (dwarf) galaxy will populate a CEP (Andrews et al. 2017). The CEP encodes information about the total mass, star formation history, and chemical abundances of its host galaxy. This is because a given (dwarf) galaxy contains stars from multiple epochs of its star formation, which will each have different chemical abundances due to the overall chemical evolution of the (dwarf) galaxy's gas reserves over time. Typically, CEPs are drawn on plots of an α-element, such as Mg or Si, against [Fe/H]; this allows one to roughly track the ratio of SNe Ia and SNe II elements in a population over time.

The bottom panel of Figure 12 shows plots of [Mg/Fe] versus [Fe/H] for the GALAH sample, split up into stars with $\sqrt{{J}_{R}}$ above and below 25 (kpc km s⁻¹)^1/2. The CEP for the high-$\sqrt{{J}_{R}}$ stars sits at a lower [Mg/Fe] than the low-$\sqrt{{J}_{R}}$ stars, which confirms that the $\sqrt{{J}_{R}}$ cut is isolating two different chemical populations in the local stellar halo.

Does the existence of two CEPs in the local stellar halo require multiple progenitors? If the progenitor of the GSE was sufficiently massive, then it may have multiple populations of stars (e.g., disk & halo) that formed in environments with different chemical abundances. This could hypothetically give the populations different CEPs. Then, the halo stars from the GSE progenitor may have been stripped earlier on than the GSE progenitor disk stars, which could potentially give them different kinematics (i.e., Figure 3).

There are multiple problems with this hypothesis. Figure 17 of Hasselquist et al. (2021) shows the stars from the Sgr dSph, split up into stars that are still bound to the dwarf, and stars that reside in the Sagittarius Stream tidal tails. It is expected that the present-day stream stars were primarily located in the outer regions of the Sgr dSph, and would have formed in a different chemical environment than the stars within the presently bound region of the Sgr dSph. However, the stream stars and bound stars are located on the same CEP: the stream stars are located on the left half of the CEP, and the bound stars are located on the right. Based on this, there is no reason to expect that any dwarf galaxy progenitor would have formed multiple different CEPs like we see in the data. One could imagine that the Sgr dSph was just not massive enough to create multiple CEPs, but this is also unlikely, because the LMC data in Hasselquist et al. (2021) do not contain multiple CEPs. However, it is once again easy to explain the existence of multiple CEPs if the local stellar halo contains debris from multiple merger events.

Figure 13 shows CEPs of the GALAH GMM components for various elements versus [Fe/H]. The CEPs of each structure are distinct from one another, except for Cronus and Thamnos, which have extremely similar CEPs for each element. This could suggest that Cronus and Thamnos are related in some way, despite their different dynamics. It is plausible that Cronus and Thamnos were accreted at roughly the same time, given their similar energies, and they could have been components of a group infall event like the LMC and its satellites today. This idea is supported by the fact that the LMC and SMC, which are known to be accreting via group infall, have extremely similar CEPs (Hasselquist et al. 2021). Additionally, the VRM and Nereus CEPs are not identical, but they are similar to each other. It is also possible that the progenitors of the VRM and Nereus were also associated with a group infall event, given their similar CEPs and energies.

Zoom In Zoom Out Reset image size

In this section, we explore the hypothesis that the identified Nereus and Thamnos components are each built up from more than one progenitor. We find that their chemical abundance distributions are more clustered than the chemical abundance distributions of the VRM and Cronus, and that there are visible clumps of stars in the Nereus and Thamnos chemical abundance distributions. These findings are consistent with Nereus and Thamnos being built up from multiple progenitor dwarf galaxies.

The large [Fe/H] abundances of the VRM and Cronus indicate that their progenitors were probably fairly massive, as it is not possible for low-mass dwarf galaxies to maintain star formation for long enough to enrich their gas with SNe Ia elements. In contrast, the lower [Fe/H] and larger [α/Fe] abundances of Nereus and Thamnos suggest that their progenitors had lower masses than the progenitors of the VRM and Cronus. Because the VRM is not phase-mixed (Donlon et al. 2020), and it is possible that Cronus is also in disequilibrium (see discussion of Figure 8), we cannot simply use the star counts of the local stellar halo to determine the ratios of the components' stellar masses.

One peculiar trend in the data is that Nereus and Thamnos stars have larger dispersion for many chemical abundances than VRM and Cronus stars. This large dispersion could be explained if Nereus and Thamnos are actually composed of the debris from many smaller merger events, rather than a single large progenitor such as the VRM and Cronus.

If Nereus or Thamnos were built up of smaller merger events, then each individual progenitor would have formed through its own star formation history in its own environment, and therefore each would have its own distinct abundance distribution compared to each of the other small progenitors. Each progenitor's abundance distribution would have a small intrinsic scatter compared to the overall population because of the progenitor's small mass. This would cause the abundance distributions of Nereus/Thamnos to contain many small clumps, each one corresponding to an individual progenitor. Conversely, if Nereus/Thamnos were formed from a single massive progenitor, then one would expect a single large distribution without any smaller abundance clumps.

One way of evaluating whether any of our halo components contain deeper substructure is by quantitatively measuring the amount of clustering, or "clumpiness," of the stellar abundances within each component. To this end, we utilize a form of Kullback–Leibler Divergence (KLD; Kullback & Leibler 1951; Fabricius et al. 2021), which is defined as

$$\begin{eqnarray}&&{\rm{KLD}}={\int }_{A}P({\boldsymbol{x}})\mathrm{log}\left(\displaystyle \frac{P({\boldsymbol{x}})}{Q({\boldsymbol{x}})}\right)\ {d}^{2}x,\end{eqnarray} \tag{ 8 }$$

where P( x ) is the two-dimensional joint probability density function, and Q( x ) = Π_i P_i (x_i ), which is the product of the marginalized one-dimensional probability density functions. The KLD is larger for more clustered distributions and smaller for less clustered distributions, which allows us to measure how clustered a given two-dimensional distribution is.

For each halo component, we numerically computed the KLD for each combination of abundances within the APOGEE and GALAH data sets. The results of this are shown in Figure 14; in the APOGEE data, the KLD of Nereus abundance distributions are much larger than the KLD for VRM stars, and in the GALAH data the mean KLD of Nereus is larger than the KLD of Thamnos, which is larger than the KLDs for Cronus and the VRM. This indicates that the abundance distributions of Nereus and Thamnos stars are more clustered than the abundance distributions of VRM and Cronus stars. This supports our hypothesis that Nereus and Thamnos are actually built up of many small accretion events, rather than a merger event with a single progenitor. Because Nereus has an even higher KLD than Thamnos, it is reasonable to believe that it is built up from more progenitors than Thamnos.

Zoom In Zoom Out Reset image size

This clumpiness is apparent in abundance diagrams; one good example is [Al/Fe] versus [Fe/H], which is shown in Figure 15. While the majority of VRM and Cronus stars cluster into a single group, the Nereus and Thamnos stars clump into several apparent substructures. It is easiest to explain the multiple peaks in the chemical abundance distributions of Nereus and Thamnos stars if they are built up from the material of multiple progenitors that formed independently.

Zoom In Zoom Out Reset image size

These results are in tentative agreement with the results of Koppelman et al. (2019), who found that Thamnos has two distinct populations. While they did not find any evidence that the different populations in Thamnos were required to be from different progenitors, we believe that this new evidence warrants further consideration of the idea that Thamnos is not built up from a single merger event.

Determining the statistically significant number of abundance clumps within Nereus and Thamnos, the relative size of each clump, and if/how these clumps correspond to different progenitors will require detailed analysis. These questions have important implications for the accretion history of the MW, including potentially constraining the number of small merger events with radial orbits that happen at various times, the average (stellar) mass of these minor merger events, and at what times these minor mergers happen.

In this section, we compare our conclusions with those of other literature. Importantly, we show that we are able to explain the results of Nissen & Schuster (2010) as a combination of VRM and Cronus stars, rather than separate populations of accreted and in situ stars.

9.1. Nissen & Schuster (2010)

One major work that analyzed potential populations of halo stars was that of Nissen & Schuster (2010) (hereafter "NS10"), who used the kinematics and [Mg/Fe], [Ni/Fe], and [Na/Fe] abundances of halo stars in order to show that there were two distinct populations of stars in the local stellar halo. They found that there was a population of stars with low [α/Fe] abundances and net retrograde motion, as well as a chemically distinct population of stars with high [α/Fe] abundances and net prograde motion.

These two halo populations have since been interpreted as accreted (low-α) and "in situ" (high-α) populations within the halo (Belokurov & Kravtsov 2022; Buder et al. 2022), and the NS10 low-α stars have even been used as a guide for selecting stars that belong to the GSE (Buder et al. 2022).

However, we offer an alternative interpretation of the NS10 data. First, we remove all NS10 stars that were classified as belonging to the thick disk, as well as any NS10 stars that are above the Splash/thick disk chemical cut given in Equation (3). Then, instead of splitting the NS10 stars into a high- and a low-[α/Fe] population, we split the stars based on whether they are located above or below the line

$$\begin{eqnarray}&&[\mathrm{Ni}{\rm{/}}\mathrm{Fe}]=1.05\ [\mathrm{Na}{\rm{/}}\mathrm{Fe}]+0.08.\end{eqnarray} \tag{ 9 }$$

This separation is shown in the top left panel of Figure 16. It is apparent in the top middle panel of Figure 16 that these two separate populations appear to trace out two distinct chemical evolution paths. These two CEPs are stacked one on top of the other, which is reminiscent of the bottom panel of Figure 12. This connection is confirmed in the top right panel of Figure 16, where the two CEPs have different $| {v}_{R}| \propto \sqrt{{J}_{R}}$ distributions.

Zoom In Zoom Out Reset image size

These two selections are consistent with the VRM and Cronus stars, as we showed earlier that the VRM and Cronus stars are split in $\sqrt{{J}_{R}}.$ The similarity between the two CEPs in the NS10 data and the GALAH VRM and Cronus components is shown in the bottom row of Figure 16. The [Ni/Fe] abundances of NS10 stars in the bottom row of Figure 16 are adjusted by −0.04 dex, to account for the mean offset between NS10 abundances and GALAH data (Buder et al. 2022). There does not appear to be a corresponding group of Nereus or Thamnos stars in NS10, because NS10 only collected data for stars with [Fe/H] > −1.6, and Nereus and Thamnos stars have low [Fe/H].

9.2. Hierarchical Clustering

9.3. Other Literature

Helmi et al. (2018) (hereafter H+18) was one of the two papers that are credited with the discovery of the GSE structure. H+18 provided the first [α/Fe]-[Fe/H] diagram of APOGEE GSE stars in the local solar neighborhood. Figure 17 shows a comparison of the H+18 data with our APOGEE GMM halo components. Overall, the H+18 data appear very similar to our VRM component, except that their stars extend up to higher [α/Fe] and lower [Fe/H] than the VRM. The population of Nereus stars with −2.0 < [Fe/H] < −1.5 and 0.15 < [α/Fe] < 0.3 are mostly missing from the H+18 data, although the very [α/Fe]-rich Nereus stars are present in the data. These differences could be due to more APOGEE data having been released since H+18, or it could be due to the more restrictive kinematic cuts made by H+18.

Zoom In Zoom Out Reset image size

Das et al. (2020) found a "blob" of stars in APOGEE data, which they claimed was composed of accreted material based on their chemical abundances. This blob has been associated with the GSE in later works (for example, Buder et al. 2022), and the chemical abundances of this blob are consistent with those of the stellar halo components identified in this work. Das et al. (2020) found that stars that are located toward the bottom right of a [Mg/Fe]-[Fe/H] plot (where we expect most of our VRM stars to be located) have smaller ages on average than stars up and to the left on a [Mg/Fe]-[Fe/H] plot (where Cronus, Nereus, and Thamnos stars are found). This supports our claim that the progenitor of the VRM was accreted recently and was able to support an extensive period of star formation.

Buder et al. (2022) also use a GMM procedure to identify populations of stars in GALAH data; however, their GMM models only fit using chemical abundances, without any kinematic information. They find that only two of these components separate meaningfully from the others: an accreted component they identify as the GSE, and an "in situ" component with high [Na/Fe] abundances. This accreted component has chemical abundances similar to those of the VRM, Cronus, and Nereus. Interestingly, in one of their GMM fits (panel "e" in their Figure 6), their GMM appears to fit two distinct low-[Fe/H] and high-[Fe/H] components to the accreted stars. The low-[Fe/H] component probably corresponds to Nereus/Thamnos, while the high-[Fe/H] component probably contains the VRM and Cronus. In general, the Buder et al. (2022) GMM procedure is not able to separate the VRM, Cronus, Nereus, and Thamnos components from one another (though the stated goal of Buder et al. (2022) was to identify accreted stars within the halo, and not to identify substructures within those accreted stars). This indicates that the orbital dynamics of the local halo stars provide essential information for disentangling the multiple components of the local stellar halo.

Buder et al. (2022) discusses the differences between the typical dynamical cut that is used to isolate GSE stars versus a chemical cut that was used to isolate likely GSE stars. They find that there is some overlap between the two selections, but that the dynamical cut and the chemical cut select different samples of stars. Buder et al. (2022) interpret that this is because the overlap between the dynamical cut and the chemical cut provides a sample of GSE stars with low contamination and relatively low completeness, whereas the other two cuts have high contamination but higher completeness. We argue that these cuts are in reality selecting different structures within halo stars: the chemical cuts select the majority of the accreted material in the halo, but the kinematic cuts are specifically isolating high-$\sqrt{{J}_{R}}$ stars (which primarily belong to the VRM debris), so it makes sense that the two selections do not produce identical samples.

Han et al. (2022a) find that the MW's stellar halo is doubly broken, with one stellar break located at r = 12 kpc, and another located at r = 28 kpc. They interpret these distances as pericenters in the orbit of the GSE progenitor, where large amounts of material would have been stripped from the dwarf galaxy. We offer an alternative interpretation of this result; the inner stellar break could be due to apocenter pileup of Cronus material, while the outer stellar break could be due to apocenter pileup of VRM stars. This is consistent with our results that Cronus stars have lower energy than VRM stars, and that the mean orbital distance of Cronus stars is smaller than that of VRM stars. Thamnos and Nereus material likely came from multiple lower-mass progenitors, so it is unlikely that they would produce strong stellar breaks in the halo.

In the last few years, there has been an explosion of literature surrounding the identification of the "last major merger," which has since been named the Gaia-Enceladus or Gaia Sausage (GSE) merger event. This literature is based on the central idea that a galaxy with a sizeable mass compared to the proto-MW collided with our Galaxy between 6 and 10 Gyr ago, which had many profound effects on the MW that can still be seen today.

We argue two points in this work, which offer an alternative viewpoint to the GSE literature: (i) A single massive merger event at early times is not able to explain the observed chemodynamics of the local stellar halo stars, and (ii) The local stellar halo stars can be separated into (at least) four components, which are the Virgo Radial Merger (VRM), Nereus, Cronus, and Thamnos. We claim that all of the phenomena laid out in this paper can be explained by the scenario where many dwarf galaxies were accreted onto the Milky Way on more-or-less radial orbits over the Galaxy's lifetime, and are now mostly phased-mixed within the inner halo (except for the VRM).

In order to analyze local stellar halo stars, we removed the thick disk/Splash stars from our sample via a chemical abundance cut. These thick disk/Splash stars appear shifted to lower [Fe/H] and higher v_ϕ compared to the selection box in Belokurov et al. (2020). This group of high-v_ϕ Splash stars does not appear to bend into the typical disk distribution because we made kinematic cuts that removed stars with disk kinematics from the data.

Our claim that a single, massive, early merger is inconsistent with observations is based on a number of chemodynamic trends within APOGEE and GALAH data. We create qualitative models for the infall of a single massive, ancient merger event with an internal metallicity gradient. In this scenario, the accreted stars with lower-[Fe/H] debris would have larger energy and L_z with respect to the host galaxy than the accreted high-[Fe/H] stars. In both surveys, we show that halo stars with lower (higher) energy have chemical abundances that trend the opposite way compared to our single merger model. Additionally, the L_z -abundance trends of the observed data would require a prograde orbit of the GSE progenitor, which is not consistent with the results of other studies. These trends cannot be explained by our model of the accretion of a single massive, ancient dwarf galaxy such as the GSE. However, it is simple to explain these chemodynamic trends if the local stellar halo actually contains debris from multiple merger events.

We show that stars with high $\sqrt{{J}_{R}}$ and stars with low $\sqrt{{J}_{R}}$ lie on different chemical evolution paths. Additionally, the CEPs of the four identified halo components are distinct, except for Cronus and Thamnos, which may be related in some way. The existence of multiple chemical evolution paths, each with different kinematics, is not consistent with a single, massive, early merger event. Even if one assumes that the GSE contributes the majority of stars in the local stellar halo, this result indicates that the kinematic selection of $\sqrt{{J}_{R}}\gt 30$ (kpc km s⁻¹)^1/2, which is often used to select GSE stars, does not produce a representative sample of the accreted stars near the Sun. We are able to explain the existence of several chemical evolution paths by the existence of debris from multiple merger events in the local stellar halo.

Our claim that the stars in the local stellar halo can be separated into (at least) four components is based on a Bayesian Gaussian mixture model regression algorithm, which finds that the best model fit for the GALAH local stellar halo stars contains four components. The GMM algorithm preferred two components for the APOGEE data due to the different footprints of the two surveys, which limits the energy range of APOGEE stars. We interpret these components as inner halo substructure, and each component corresponds to a particular merger event (or set of merger events). The identified substructures are summarized here:

• 1.  
  VRM: A high-[Fe/H], low-[α/Fe] substructure with large energy and a small retrograde rotation. The high [Fe/H], low [α/Fe], and large energy of this structure suggest that it was likely accreted recently. The VRM is closest to the typical characterization of the GSE, as it has high orbital eccentricity, creates a double-lobed v_r -v_ϕ velocity structure in the solar neighborhood, and has large $\sqrt{{J}_{R}}$.
• 2.  
  Cronus: A high-[Fe/H], high-[α/Fe] substructure with low energy and prograde rotation. The large [α/Fe] values of the Cronus stars indicate an early time of accretion, which is also consistent with Cronus stars' low energy.
• 3.  
  Nereus: A low-[Fe/H], high-[α/Fe] substructure with intermediate energy, and little-to-no net rotation. Nereus stars form several clumps in abundance space, which indicates that Nereus is probably made up of material from several smaller progenitors that formed independently and then fell in at relatively early times.
• 4.  
  Thamnos: This component has a relatively low energy, relatively low [Fe/H], and contains stars on strongly retrograde orbits. Thamnos stars' chemical abundances are also clumpy and have a wide dispersion, which is in agreement with the idea that Thamnos came from more than one progenitor dwarf galaxy.

We compare our results to previous analyses of halo stars. In general, these studies assume that all accreted stars arise from a single structure. However, the data from Nissen & Schuster (2010) is remarkably consistent with our multiple-progenitor picture when split into chemical evolution paths rather than high and low [α/Fe].

Some of the typical cuts made on the halo to select GSE stars, such as selecting all stars with orbital eccentricity >0.8, are really selecting a combination of VRM, Cronus, and Nereus stars. Conversely, some selections, such as $\sqrt{{J}_{R}}\gt 30$ (kpc km s⁻¹)^−1/2, are preferentially selecting for only one of these components. Treating these selections as identical makes sense when one assumes that the entire local stellar halo is composed of debris from a single merger event, but in reality these various selections correspond to different mixtures of substructures within the local stellar halo, which can result in different observed properties depending on how the data was selected.

The Milky Way's stellar halo was originally considered to be smooth in density; now, we know that the stellar halo contains a multitude of tidal streams, which provide evidence of many accretion events over a large span of time. It makes sense that the radial portion of the stellar halo would be formed in a similar way, with many accretion events over the course of the Milky Way's history, rather than a single merger event depositing the majority of the material in the inner halo. The picture of a single ancient, massive merger event dominating the MW's inner halo is no longer sufficient to explain the wealth of chemodynamic substructure near the Sun, and the figurative pendulum must swing back to a view where the local stellar halo was built up over time from multiple merger events.

This work was supported by NSF grant AST 19-08653; contributions made by Manit Limlamai; and the NASA/NY Space Grant.

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

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics ∣Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

This work made use of the Third Data Release of the GALAH Survey (Buder et al. 2022). The GALAH Survey is based on data acquired through the Australian Astronomical Observatory, under programs: A/2013B/13 (The GALAH pilot survey); A/2014A/25, A/2015A/19, A2017A/18 (The GALAH survey phase 1); A2018A/18 (Open clusters with HERMES); A2019A/1 (Hierarchical star formation in Ori OB1); A2019A/15 (The GALAH survey phase 2); A/2015B/19, A/2016A/22, A/2016B/10, A/2017B/16, A/2018B/15 (The HERMES-TESS program); and A/2015A/3, A/2015B/1, A/2015B/19, A/2016A/22, A/2016B/12, A/2017A/14 (The HERMES K2-follow-up program). We acknowledge the traditional owners of the land on which the AAT stands, the Gamilaraay people, and pay our respects to elders past and present. This paper includes data that has been provided by AAO Data Central (datacentral.org.au).

Software: gaiadr3_zeropoint (Lindegren et al. 2021), McMillan2017 (McMillan 2017), galpy (Bovy 2015), scikit-learn (Pedregosa et al. 2011), astropy (Astropy Collaboration et al. 2013, 2018)

The z-component of angular momentum can be computed as

$$\begin{eqnarray}&&{L}_{z}={{xv}}_{y}-{{yv}}_{x}.\end{eqnarray} \tag{ A1 }$$

The calculated Cartesian Galactocentric positions (x, y) and velocities (v_x , v_y ) of each star are functions of distance from the Sun, and are also dependent on the choice of the Sun's phase-space coordinates. Therefore, it follows that L_z will have a dependence on these values.

Donlon et al. (2022) (hereafter D+22) used a uniform parallax zero-point offset of −0.017 mas to correct the Gaia EDR3 parallaxes; in this work, we use the parallax zero-point offset calibration provided by Lindegren et al. (2021). These different choices of parallax zero-point offsets will result in slightly different calculated distances to each star, which will change the computed L_z value to that star (see Helmi et al. 2018 for an example of this).

Similarly, D+22 use a v_LSR of 220 km s⁻¹, a Solar position of 8 kpc from the Galactic center, and a Solar reflex motion of (10.1, 4.0, 6.7) km s⁻¹. This work uses a v_LSR of 233.1 km s⁻¹, a Solar position of 8.21 kpc from the Galactic center, and a Solar reflex motion of (1.1, 15.17, 7.25) km s⁻¹. These choices of kinematic frame will also generate differences in the measured L_z values of each star.

It is helpful to determine the impact of this choice of frame on measured values of L_z in order to see how the results of this work compare to the results of D+22. To this end, we define the change in angular momentum,

$$\begin{eqnarray}&&{\rm{\Delta }}{L}_{z}={L}_{z,{\rm{D+22}}}-{L}_{z,{\rm{This}}\,{\rm{Work}}},\end{eqnarray} \tag{ A2 }$$

as well as the change in measured distance,

$$\begin{eqnarray}&&{\rm{\Delta }}d={d}_{{\rm{D+22}}}-{d}_{{\rm{This}}\,{\rm{Work}}}.\end{eqnarray} \tag{ A3 }$$

Figure 18 shows the relationship between these values. In general, the larger Δd, the larger the corresponding ΔL_z . Note the overall offset of ΔL_z ∼ 200 kpc km s⁻¹ at Δd = 0: this is because the change in the chosen Solar phase-space coordinates produces an offset in the measured L_z .

Zoom In Zoom Out Reset image size

The variation in L_z is also dependent on other parameters, such as position on the sky. This is shown in Figure 19, where it is clear that the measured change in L_z is smallest near the Galactic center, and is largest in the regions near the Galactic poles and the anticenter. The value of ΔL_z also varies based on the apparent magnitude of each star, because the parallax zero-point calibration has a complex dependence on the apparent magnitude of each object (Figure 3 of Lindegren et al. 2021).

Zoom In Zoom Out Reset image size

D+22 found that the VRM had a net L_z of ∼400 kpc km s⁻¹, whereas Nereus had a net L_z of ∼0 kpc km s⁻¹. It is important to make sure that this difference in L_z is physical and not simply due to the choice of kinematic frame and parallax zero-point offset. Figure 20 shows the cumulative distribution functions of ΔL_z split up for VRM and Nereus stars. There is a slightly smaller average ΔL_z value for Nereus stars, i.e., the L_z values of Nereus stars in D+22 were too small compared to the L_z values of VRM stars, when comparing to L_z values in the frame used in this work. This is likely due to the fact that Nereus stars are on average dimmer than VRM stars (shown in the bottom panel of Figure 20), so Nereus stars will have smaller distances (and also smaller Δd) on average compared to VRM stars, and therefore a smaller ΔL_Z as per Figure 18. However, this effect is small (only on the order of 10 kpc km s⁻¹) and cannot explain the ∼400 kpc km s⁻¹ difference between the two populations measured in D+22.

Zoom In Zoom Out Reset image size

We conclude that the choice of frame does have a sizeable effect on measured L_Z values, but it is not able to explain the differences in L_z values of VRM and Nereus stars observed in D+22. The variations between the measured L_z of the VRM and Nereus in this work and D+22 could be due to multiple factors: there is a large dispersion in the measured L_z of these structures in this work, which could "wash out" differences between the two substructures; this work primarily studies giant stars, whereas D+22 studies dwarf stars, and the measured L_z of the two populations could be different in some way; or the actual L_z distribution of the merger debris for Nereus and the VRM varies sizeably on the sky, and the regions of the sky that are probed by this work are different enough from those analyzed in D+22 that there is a measured difference in the L_z of each merger event.

In order to evaluate whether our GMM algorithm is able to recover individual dwarf galaxies as components, we fit our model on a data set in which it is known to which dwarf galaxy each star belongs.

For this exercise, we used the data from Hasselquist et al. (2021) (hereafter H+21), by matching their Table 2 to the APOGEE DR17 allStar catalog. Kinematic parameters were calculated using the APOGEE line-of-sight velocities and the Gaia EDR3 proper motions and parallax values provided in the catalog. This yielded a data set with 6D kinematics and APOGEE chemical abundances for 8005 stars, where H+21 had assigned each star to one of the following structures: the LMC, the SMC, the Sgr dSph plus its tidal tails, Fornax (Fnx), or the GSE.

Figure 21 shows the pseudo-energy $\widetilde{E}$ versus L_z for stars in the data set with positive parallaxes and relative parallax uncertainties smaller than 100%. We only show GSE stars with relative parallax uncertainties below 15%, because GSE stars can be found closer to the Sun than the other structures in the data. The distributions of kinematic quantities for the LMC, SMC, and Sgr stars are strongly non-Gaussian. This is probably because there are substantial distance uncertainties for most stars; less than 10% of stars in the data have relative parallax uncertainties below 15%, and large uncertainties in distance will cause corresponding uncertainties in the measured values of $\widetilde{E}$ and L_z . It is difficult to determine more accurate distances to these stars than those calculated from their Gaia parallaxes, because they do not share a common spectral type and they do not belong to a single population.

Zoom In Zoom Out Reset image size

Our GMM algorithm only works well for quantities that are not correlated with one another, which does not appear to be the case for $\widetilde{E}$ and L_z of the H+21 data due to their large distance uncertainties. However, the measured chemical abundances of the APOGEE data are of high quality. In order to utilize the chemical abundances of these stars, we generated a mock set of kinematic data that are guaranteed to have no correlation between $\widetilde{E}$ and L_z . We assume that if we had accurate kinematic information for all of the stars in the data set, the $\widetilde{E}$–L_z distributions of a given dwarf galaxy would be roughly Gaussian. Each structure in the H+21 data was assigned a two-dimensional normal distribution of $\widetilde{E}$ and L_z , where the mean values and variances were chosen to be roughly analogous to those shown in Figure 21. Then, each star in the data set was given random $\widetilde{E}$ and L_z values; these values are pulled from the distribution of the structure to which that star is assigned.

The GMM algorithm was then run on the mock kinematic values for the data, along with their measured chemical abundances. We followed the procedure of Section 5, i.e., we used the [Fe/H], [α/Fe], [Na/Fe], and [Al/Fe] abundances of the data to fit the model. The best-fit model contained six components: one component corresponded to each of the H+21 structures, plus an additional "background" component with a large dispersion in kinematic quantities and chemical abundances but a small number of stars (less than 5% of the stars in the data set). Then, following Section 5, each star was assigned to a GMM component.

Figure 22 shows the mock kinematic distributions of the structures in H+21, as well as the mock kinematic distributions of the fit GMM components. The GMM algorithm does a good job of recovering the underlying components of the distribution; however, this may not be particularly surprising or interesting considering a GMM was used to generate the mock kinematic data. The main difference between the true distributions and the fit distributions is that the algorithm incorrectly assigned several of the LMC stars to the SMC distribution.

Zoom In Zoom Out Reset image size

Figure 23 shows the chemical distributions of the H+21 structures compared to the fit distributions. Overall, the fit model does a good job of recovering the true distributions of the data. As seen in the mock kinematic data, a number of stars that actually belong to the LMC are assigned to the SMC by the model. The model also struggled to associate a small number of [Fe/H]-poor, [α/Fe]-rich stars in Sgr with the Sgr component (green); this may be because the Sgr component appears to have multiple peaks in [Fe/H] and [α/Fe].

Zoom In Zoom Out Reset image size

We also ran the GMM algorithm on the H+21 data using the measured kinematic quantities, in order to explore how robust the model is for non-Gaussian distributions. The best-fit model had the correct number of components, and it was able to assign 68% of the stars to the correct components, compared to the mock kinematics, where the algorithm was able to correctly assign 78% of stars to the correct component. It is possible that some stars in the H+21 data do not actually belong to the dwarf galaxy to which they were assigned, which would reduce our algorithm's ability to correctly identify to which component a given star belongs.

We are confident that the GMM algorithm that we present in this work is able to recover distributions of stars that belong to one or more dwarf galaxies, with fairly high accuracy, even when those distributions are complicated and overlapping.