Shiva and Shakti: Presumed Proto-Galactic Fragments in the Inner Milky Way

Using Gaia Data Release 3 astrometry and spectroscopy, we study two new substructures in the orbit–metallicity space of the inner Milky Way: Shakti and Shiva. They were identified as two confined, high-contrast overdensities in the (L z , E) distribution of bright (G < 16) and metal-poor (−2.5 < [M/H] < − 1.0) stars. Both have stellar masses of M ⋆ ≳ 107 M ⊙, and are distributed on prograde orbits inside the solar circle in the Galaxy. Both structures have an orbit-space distribution that points toward an accreted origin; however, their abundance patterns—from APOGEE—are such that are conventionally attributed to an in situ population. These seemingly contradictory diagnostics could be reconciled if we interpret the abundances [Mg/Fe], [Al/Fe], [Mg/Mn] versus [Fe/H] distribution of their member stars merely as a sign of rapid enrichment. This would then suggest one of two scenarios. Either these prograde substructures were created by some form of resonant orbit trapping of the field stars by the rotating bar; a plausible scenario proposed by Dillamore et al. Or, Shakti and Shiva were protogalactic fragments that formed stars rapidly and coalesced early, akin to the constituents of the poor old heart of the Milky Way, just less deep in the Galactic potential and still discernible in orbit space.


INTRODUCTION
One focus of modern astrophysics is to understand how massive galaxies, such as our Milky Way, formed and evolved.To address this broad theme with observations of the present epoch Universe, it has proven powerful to take the Milky Way itself as a galaxy model organism for which to map -and then eventually understand -the detailed age-orbit-abundance distribution of its stars.This constrains when and on which orbits its stars formed, how the stellar birth material became gradually enriched, and which dynamical processes -e.g.mergers or secular evolution -shaped the Galactic structure.The importance of various processes must depend both on the epoch and on the Galactocentric radius (or potential well depth) at which they occur.
Parsing the overall stellar orbit-abundance distribution into an extensive set of individual components has been a long tradition, as each of them may represent a discrete event or regime in our Galaxy's formation history.In particular, it has become established to categorise each component or substructure as either born from gas already within the Milky Way's dominant potential well (dubbed in-situ) or born within a (at the time) distinct subhalo (or satellite galaxy) that was later accreted into the Milky Way.
Among old (≥ 8 Gyrs) stellar populations in the Milky Way, the prototypical in-situ population is the α-enhanced or thick disk.This inference is based on the disk-like kinematics that reflects the settling of the stars' birth gas in the Milky Way's main potential well (e.g.Birnboim & Dekel 2003;Stern et al. 2021), starting t age ∼ 12.5 Gyr ago (Xiang & Rix 2022).When looking at these stars in the (L z , E) plane of orbital parameters, their distribution follows the line of circular, in-plane orbits (see Fig. 1), though not exactly owing to their substantive velocity dispersion.And this assessment of an in-situ origin is based on the abundance patterns: these stars are enhanced in [Mg/Fe] and [Mg/Mn] and show a rapid rise of [Al/Fe] with [Fe/H].All of these are taken as signatures of rapid enrichment that presum-ably requires high densities and a deep potential well (Hawkins et al. 2016;Belokurov & Kravtsov 2022).
Among the stellar halo components, the Gaia-Sausage-Enceladus (henceforth GSE) population (Helmi et al. 2018;Belokurov et al. 2018) exemplifies an accreted population (see Fig. 1; right panel).The very low angular momentum and large apocenters of their orbits cannot be readily explained by an in-situ origin, and the low-α abundance patterns resemble those of the 'surviving' dwarfs satellites of the Milky Way and point toward slow enrichment in a more diffuse potential well (Wyse 2016;Hasselquist et al. 2021).
Yet, when looking in detail at the very early phases of our Galaxy's evolution (say ≳ 11 − 12 Gyr ago), this in-situ vs. accreted distinction may become practically ambiguous and conceptually subtle.For instance, there is the "Splash" population (see Fig. 1), identified by Bonaca et al. (2017);Di Matteo et al. (2019); Belokurov et al. (2020) that has been termed the "metal-rich insitu" halo: its orbits are radial (hence halo-like, favouring and accreted origin), but the population is metal-rich with [M/H] > ∼ − 1.0 and α-enhanced, pointing towards an in-situ scenario.But both properties can be understood if its population arose from stars that were kicked out of the early, nascent in-situ disk into radial orbits during an early collision with some massive merger.
Recent studies of ancient stars in the inner Galaxy (Kruijssen et al. 2019;Belokurov & Kravtsov 2022;Rix Two new substructures in the Milky Way et al. 2022;Horta et al. 2022;Conroy et al. 2022) have identified a dynamically hot population of tightly-bound stars, namely the Poor Old Heart of the Milky Way (or Aurora), henceforth POH (see Fig. 1).Most of these stars show the chemical signatures of rapid enrichment (i.e., these stars are enhanced in [Mg/Fe] and [Mg/Mn] and show a rapid rise of [Al/Fe] with [Fe/H]) that are conventionally attributed to an in-situ origin.Simulations of galaxy formation (e.g., Renaud et al. 2021) suggest that this component originated presumably from the coalescence of several proto-Galactic fragments that were of comparable, though not equal, mass.Whether at such early phases, it is useful to attribute the term insitu for only one of these presumed protogalactic fragments (Belokurov & Kravtsov 2022) is not clear (cf.Rix et al. 2022;Conroy et al. 2022).Nonetheless, the abundance pattern of these stars may serve as a reminder that they only imply rapid enrichment (which may require high star formation rates, deep potential wells and high densities), rather than as a direct indicator of the dynamical pre-history.

3
But also the inference path from a particular orbit distribution to, say, an accreted origin may not be straightforward or unique.For instance, as mentioned above, the "Splash" population has halo-like orbits and appear to favour an accreted origin; but, it actually comprises metal-rich stars with [M/H] > ∼ − 1.0 that suggest it was the earliest part of the α-enhanced disk whose stars were dramatically stirred during an early, quite massive merger.Likewise, the "R" feature indicated in Fig. 1, which represents ridge-like overdensity of stars, may lend itself to the interpretation as an accreted population.But Dillamore et al. (2023) showed that this feature comprises both metal-rich and metal-poor stars and their orbits lie somewhere between disk-and halo-like.Given that these stars' azimuthal orbit frequency coincides with that of the Galactic bar, this orbit-space substructure can be very plausibly attributed to the trapping of the (inner) halo stars in bar resonances.
In this study, our primary goal is to examine the (inner) Milky Way's (L z , E)-[M/H] space, using the recent Gaia DR3-based XGBoost dataset, (Rix et al. 2022;Andrae et al. 2023, A23), to look for previously unrecognised substructures.In doing so, we identified two orbitabundance space overdensities that we call Shakti and Shiva, and show these possess very intriguing properties that send "mixed messages" about their in-situ or accreted origin.To this end, Section 2 describes the data used.Section 2.2 details our method to compute the orbital parameters of the stars.In section 3, we analyze the (L z , E) distribution of stars in different slices of [M/H], thereby identifying Shakti and Shiva.Section 4 and Section 5 present the characterization of these stellar populations.In Section 6, we lay out possible origin scenarios for these two populations.We finally conclude in Section 7.

XGBoost dataset
For this, one of the most extensive and well-vetted data sets is the recent Gaia DR3-based XGBoost catalogue (Rix et al. 2022;Andrae et al. 2023, A23).A23 provides data-driven XGboost estimates of stellar metallicity [M/H] XP , surface gravity log(g) XP and effective temperature T eff,XP for 175 million sources.These are based on their Gaia DR3 low-resolution BP/RP (or, XP) spectra, after being trained on stellar parameters from APOGEE that have been augmented by a set of very metal-poor stars.Although all of these stars have Gaia DR3's 5D astrometry, only brighter sources (with G RVS < ∼ 14 mag) also have Gaia RVS velocities.A23 combines one of the most extensive data set with 6D phase-space with [M/H] that are robust, precise, and vetted against a number of external catalogs, such as APOGEE DR17 (Abdurro'uf et al. 2021), GALAH DR3 (Buder et al. 2021), LAMOST DR7 (Zhao et al. 2012) and SEGUE (Yanny et al. 2009).In our analysis, we restrict ourselves to only those stars with phot g mean mag < 16, 3800 K ≤ T eff,XP < 4800 K and log(g) XP < 3.5, because in this range XGBoost's [M/H] estimates are most robust and precise (with δ[M/H] < 0.1).Even with these quality cuts, this A23 subset contains ∼ 30 times more stars compared to APOGEE DR17 and LAMOST DR7, ∼ 70 times more stars compared to GALAH DR3 and ∼ 5000 times more stars than in SEGUE.
We apply additional quality to ensure good Gaia astrometry (similar to Malhan 2022): (1) parallax over error≥ 5; this quality cut ensures that our resulting sample contains well-measured parallaxes with relative parallax error of ≤ 20%.Effectively, it also ensures that ϖ > 0 and that the estimated heliocentric distances D ⊙ (= 1/ϖ) are physical.We correct for the parallax zero-point of each star as ϖ corrected = ϖ observed −(−0.017)(Lindegren et al. 2021).
(2) phot bp rp excess factor< 1.7; this moderate cut removes stars with strong G BP − G RP color excess that occurs in the very crowded regions of the Galaxy (for normal stars, phot bp rp excess factor≈ 1, Arenou et al. 2018).( 3) ruwe< 1.4; this value has been prescribed as a quality cut for "good" astrometric solutions by Lindegren et al. (2021).This results in a sample of 5, 915, 741 stars.Finally, we excise all the stars that lie within five tidal radii of the globular clusters (listed in Harris 2010) in the 2D projection space.This removes 116, 017 stars.
Our final 'base' dataset comprises n = 5, 799, 724 stars, and its overview is provided in Fig. 2. We correct for dust extinction using Schlegel et al. (1998) maps and assuming the extinction ratios A G /A V = 0.86117, A BP /A V = 1.06126,A RP /A V = 0.64753 (as listed on the web interface of the PARSEC isochrones, Bressan et al. 2012).Henceforth, all magnitudes and colors refer to these extinction-corrected values.

Orbits
To compute the orbital parameters of stars (such as L z , E, etc), we need to transform their heliocentric phase-space measurements into Galactocentric Cartesian coordinates.This is achieved using the Sun's Galactocentric location as R ⊙ = 8.112 kpc (Gravity Collaboration et al. 2018) and the Sun's Galactocentric velocity as (v x,⊙ , v y,⊙ , v z,⊙ ) = (12.9,245.6, 7.78) km s −1 (Drimmel & Poggio 2018).Next, the Galactic potential is set up using the McMillan (2017) model, which is a static and axisymmetric model comprising a bulge, disk components and an NFW halo.With this potential, we compute the orbital parameters, using the galpy module (Bovy 2015).and c).It also shows clearly some of the well-established low-metallicity 'sub-structures', such as the metal-rich in-situ halo/Splash (Bonaca et al. 2017;Di Matteo et al. 2019;Belokurov et al. 2020, visible in the panel c), GSE merger debris (Belokurov et al. 2018;Helmi et al. 2018, visible in the panels d and e) and POH (Belokurov & Kravtsov 2022;Rix et al. 2022, visible in the panels d and e).
Panels d and e of Fig. 3 (which together correspond to the metal-poor regime −2.0 ≤ [M/H] XP < −1.0) reveal the presence of three prominent substructures that are well confined in (L z , E) space.One of them with net L z ≈ 0 corresponds to the known GSE structure.The other two prograde substructures at (L z , E) ∼ (−1000 kpc km s −1 , −1.7×10 5 km 2 s −2 ) and at (L z , E) ∼ (−300 kpc km s −1 , −1.9 × 10 5 km 2 s −2 ) have previously not been observed as such prominent overdensities or appreciated as possible substructures.They are the focus of the present work, and for brevity, we denote them as Shakti and Shiva respectively 1 .The fact that these substructures are observed here for the first time with such high contrast, owes to the novel richness of the XGBoost dataset.
In the remaining study, we characterize the stellar populations of Shakti and Shiva and aim to understand their nature and origin (i.e., in-situ or accreted).At this juncture, it is important to note that the (L z , E) location of Shakti is similar to two previously identified substructures: the ridge-like component pointed out in Dillamore et al. (2023), who attribute them to halo stars that got spun up as they got trapped in the resonances of the Galactic bar (see our Fig.1); and also a sparse population of metal-poor disk stars (e.g., Carollo et al. 2019;Naidu et al. 2020).To what extent these substructures are identical, overlapping or connected with Shakti is a question we address in Section 6.

CHARACTERIZATION OF SHAKTI
Below, Section 4.1 describes our procedure to select the Shakti population, Section 4.2 details its chemodynamical properties and in Section 4.3 we compute its stellar mass.

Selection of Shakti stars
We select Shakti stars drawing on simple visual inspection of the n L z , E | [M/H] XP distribution, inspired by the high-contrast appearance among low-metallicity stars (see panel d of Fig. 3).Initially, we had tried for automated and systematic selection by employing some detection algorithms (namely, ENLINK, Sharma & Johnston 2009, and Gaussian mixture model, Pedregosa et al. 2011) and also experimented with different data sub-spaces (e.g., the 2D (L z , E) space, the 3D action (J R , J ϕ , J z ) space, the 4D (J, E) space, the 3D (L z , E, [M/H] XP ) space).We did not find that any of these served our specific purposes better than selection "by hand".

Orbital and chemical properties of Shakti
Fig. 4 shows the distribution of the Shakti stars in the orbital space, spatial coordinates and chemical space.These distributions are quantified in Table 1, with uncertainties derived from bootstrap re-sampling.The basic properties can be summarized as follows.
Shakti stars orbit the Milky Way within the Galactocentric distance range of D gal ∼ 2−10 kpc (Fig. 4, panel b); they reach z max ∼ 1 − 6 kpc away from the Galactic plane, and have moderate eccentricities of e ∼ 0.1 − 0.6.Furthermore, these stars are well phase-mixed in the Galaxy (as shown in the R − z and ℓ − b distributions).
We now turn to understanding Shakti's detailed abundance patterns.We initially had selected Shakti members via n L z , E | [M/H] XP .Now we aim to construct the detailed element ratios p [X 1 /X 2 ]|Shakti member that can serve as diagnostics to differentiate Galactic populations (e.g., Hawkins et al. 2015;Belokurov & Kravtsov 2022).Since Gaia DR3 cannot deliver these2 , we cross-match the Gaia-identified Shakti sample members with SDSS/APOGEE DR17 (Abdurro'uf et al.The reason behind this is unclear, however, but it may be pointing towards a single, chemically simple origin for the Shakti stars, rather than a mix of different populations.While CMD-based analyses are generally quite useful to study stellar populations, here this is not very helpful for two reaons: all Shakti stars were chosen to be metal-poor, and hence by implication old, population of giant stars (see Section 2 and panels e and f of Figure 2).Second, for most members the  Note-The "median" represents the median of the distribution and the corresponding uncertainties represent their 16 and 84 percentiles.The "spread" is the inter-quartile difference between the 16 and 84 percentiles and the corresponding uncertainties represent their 16 and 84 percentiles.
parallaxes are not precise enough to enable meaninful age estimates.
How the above properties of Shakti relate to its nature and origin is discussed below in Section 6 (together with that of Shiva).

Shakti's Stellar Mass
Obviously, the sample members that we attribute to Shakti are a highly biased subset of that structure's stars, biased in luminosity and orbital phase.
If Shakti reflected an accreted satellite, we could employ the much-used mass-metallicity relation for dwarf galaxies (Kirby et al. 2013).Adopting a mean [Fe/H] APOGEE for Shakti would imply M ⋆ > ∼ 10 7 M ⊙ , as this mean metallicity resembles that of the known WLM dwarf galaxy, which has [Fe/H]≈ −1.28 and M ⋆ ∼ 1.7 × 10 7 M ⊙ (Zhang et al. 2012;McConnachie 2012).It is important to note that this [Fe/H]-M ⋆ relation possesses significant scatter and also a redshift dependence (i.e., higher masses at higher redshifts at fixed [Fe/H], e.g., Zahid et al. 2014;Ma et al. 2016).
Alternatively, we can simply count Shakti's stars, accounting for the fact that we only consider giants; and that we will not see its members at all orbital phases.Indeed, Frankel et al. (2023) suggested to exploit the fact that for phase-mixed structures the distribution in all three orbital angles Θ is expected to be uniform.If we can estimate the "angle incompleteness" λ −1 sample , then we can estimate Shakti's total stellar mass, M total ⋆ via where N sample is the number of (giant) star sample members that we attribute to Shakti (Fig. 4a), λ sample is the factor accounting for the sample incompleteness3 and "M ⋆ /giant in sample" is the stellar population mass per giant star (at a given absolute magnitude M G0 < M G0,max ).These parameters are explained and computed in Appendix A. In this way, we arrive at a stellar mass estimate for Shakti of M total ⋆ ∼ 0.7 × 10 7 M ⊙ .

CHARACTERIZATION OF SHIVA
For characterizing the Shiva population, we follow the same approach as that described above for Shakti (in Section 4), except with slight modifications that are most suitable to study this substructure.
The selection of Shiva stars is again guided by the appearance of a prominent overdensity in panel d of Fig. 5 shows the distribution of Shiva stars in the orbital space, spatial coordinates and the chemical space, and these distributions are quantified in Table 1.These stars orbit the Milky Way within the Galactocentric distance of D gal ∼ 0 − 8 kpc (see r peri − r apo distribution), with the amplitudes perpendicular to the Galactic plane ranging between z max ∼ 1 − 6 kpc and possess a moderately high range of eccentricities with e ∼ 0.2 − 1.0.These stars are phase-mixed in the Galaxy (as shown in the R − z and ℓ − b distributions).To compute Shiva's stellar mass, if we base an estimate on the mass-metallicity relation of dwarf galaxies (Kirby et al. 2013), we arrive at M ⋆ > ∼ 10 7 M ⊙ : e.g., similar to the WLM known WLM dwarf galaxy, which has [Fe/H]≈ −1.28 and M ⋆ ∼ 1.7 × 10 7 M ⊙ (Zhang et al. 2012;McConnachie 2012).On the other hand, if we count the Shiva member stars and correct them for (orbital) angle incompleteness and the sample's luminosity cut, we find M ⋆ ∼ 2.5 × 10 7 M ⊙ , somewhat comparable to that of the POH substructure (Rix et al. 2022)

ON THE "NATURE" AND ORIGIN OF SHAKTI AND SHIVA
As we show below, Shakti and Shiva populations possess quite unconventional, intriguing properties that have not been observed previously in any substructure.This renders it quite challenging to understand the true nature and origin of these substructures.Below, we list down their properties and simultaneously interpret them to understand which of the following origin scenarios they possibly favour -(a) in-situ/disk origin (i.e., corresponding to the metal-rich disk population), (b) "proto-Galactic" origin (i.e., representing Aurora/POHlike metal-poor population, Belokurov & Kravtsov 2022;Rix et al. 2022), (c) dynamically "heated" disk origin (i.e., corresponding to the stars that were kicked out of the disk during an early collision with some massive merger, Bonaca et al. 2017;Di Matteo et al. 2019), (d) "spun up" halo population (i.e., representing those halo stars that got trapped in a ridge-like component in resonances with the Galactic bar, Dillamore et al. 2023) or (e) accreted origin (i.e., representing a population that was accreted inside a foreign galaxy, much like the GSE).
The following discussion is summarised in Table 2.
In doing so, we also require to compare the chemical abundances of Shakti and Shiva samples with those of GSE (which represents a prototypical massive merger) and POH.For this, we construct the GSE sample using our base data and follow Limberg et al. ( 2022 6 and are interesting for the following reasons. First, accreted metal-poor stars (with [Fe/H]< −1.0) are presumed to stem from mergers that are less massive and less dense than the Milky Way (e.g., Helmi 2020).Consequently, such accreted stars should show chemical signatures of slower enrichment inside the low-mass progenitors (Wyse 2016;Hawkins et al. 2015), which is reflected in their low [Al/Fe] ( < ∼ −0.2) and with practically all the stars possessing [Al/Fe]< 0. This latter point is based on Hasselquist et al. (2021) who analysed the five most massive dwarf galaxies accreted by the Milky Way -LMC, SMC, Sagittarius, Fornax and GSE; also see Horta et al. (2022).
Second, metal-poor stars whose [α/Fe], [Al/Fe] and [Mg/Mn] abundances point towards rapid enrichment, most likely would have happened within a massive and dense progenitor.Recently, Belokurov & Kravtsov (2022) analysed the APOGEE DR17 dataset, and dubbed the entire population of metalpoor ([Fe/H] < ∼ − 1.0) and Al-rich ([Al/Fe]> −0.1) stars as Aurora, attributing them to correspond to the Milky Way pre-disk in-situ population that initially surrounded the Milky Way in a kinematically hot spheroidal configuration (and some of which later settled into the disk configuration, Sestito et al. 2020), but at the present-day is sparsely distributed within the solar radius.This population may overlap with the POH identified by Rix et al. (2022) as a dense central concentration of metal-poor, α-enhanced stars.Rix et al. (2022) suggest that these stars need not come from a single in-situ progenitor, but may come from a few distinct dense and massive progenitors that coalesced at high-redshift (T > ∼ 12.5 Gyr ago) to form the proto-Galaxy.
To examine whether Shakti's and Shiva's abundance patterns look similar to the accreted-GSE or to Aurora/POH we perform Kolmogorov-Smirnov (KS) test for the null hypothesis that the given two samples are drawn from the same distribution.Shakti and Shiva samples are statistically very similar to that of POH, but very different than that of GSE (the null hypothesis can be rejected at > 99% level).Furthermore, these two substructures are even more enhanced in higher-order elements than LMC and SMC (for instance, compare our Fig.6 with Figure 5 of Hasselquist et al. 2021).Moreover, the abundances of Shakti and Shiva are also very similar to each other, suggesting a similar progenitor or formation history.
In sum, the fact that Shakti and Shiva are [Fe/H]poor but show abundance patterns indicative of rapid enrichment, similar to Aurora/POH, favours a "proto-Galactic" origin.Perhaps Shakti and Shiva were very early distinct entities that were simply more rapidly enriched than GSE (and also than LMC and SMC), perhaps they were more dense; at least at the time of accretion.
• Both Shakti and Shiva populations are observed as high-contrast, tightly-bound overdensities in the (L z , E) space.Furthermore, both possess prograde motion, with Shakti possessing a much higher |L z | value (similar to that of the disk at the given E) and Shiva possessing a lower |L z | value.Moreover, they are more tightly bound in the Milky Way potential (E) than GSE but less tighly bound that the POH.These properties can be seen in Fig. 3.
Observations of such prominent overdensities in the (L z , E) space have conventionally been considered unambiguous telltale signs of accreted populations (e.g., Helmi et al. 2018;Naidu et al. 2020;Yuan et al. 2020), that were once born in a potential well that was distinct from the (current) Milky Way.
While Shakti's high angular momentum makes it tempting to connect it with the in-situ/disk origin, this scenario seems implausible given its [Fe/H]< −1.0.Secondly, while its (L z , E) location is similar to that of Dillamore et al. (2023)'s ridge-like component, its compact "blob"-like appearance in (L z , E) space argues against them being a ridge of former halo stars that became trapped in resonances with the Galactic bar.The feature that Dillamore et al. ( 2023) identified had a lower contrast than Shakti.If Shakti were part of a resonance-driven ridge in orbit space, we would expect the chemical properties to be similar to that of the halo stars, in particular to GSE (because the majority of halo stars presumably originated from GSE alone, Helmi 2020).But Shakti is chemically different from GSE, as shown by the abundance data.Furthermore, we note that Shakti's location in orbit space is broadly similar to a postulated "metal-poor disk" (e.g., Carollo et al. 2019;Naidu et al. 2020).However, this feature is much more extended and of low contrast than Shakti.
If Shakti and Shiva were once distinct protogalactic fragments, they must have been incorporated ≳ 12 Gyr ago; at least their star formation must have preceded the oldest and most metal-poor part of the disk (Xiang & Rix 2022).This is qualitatively supported by the idea of taking the binding energy (E) as a proxy for the accretion time, which places them between the POH ( > ∼ 12.5 Gyr, Rix et al. 2022) and GSE ( > ∼ 11 Gyr, Belokurov et al. 2018;Helmi et al. 2018).
• Both Shakti and Shiva populations are distributed within D gal < ∼ 8 − 10 kpc of the Milky Way, appear phase-mixed, and possess moderately large z max and eccentricity values.These properties can be seen in Figures 4 and 5.This indicates that the stars lie on dynamically hot orbits, hotter than the oldest part of the (α-enhanced) disk.This may be consistent with both the accreted and the "heated" disk scenarios.However, the latter scenario seems unlikely, given that Shakti and Shiva are observed as high-contrast orbit-substructures (and not sparse distribution of stars, such as that observed in the Splash substructure).Even their [M/H]-deficiency opposes any plausible connection with the disk.
In view of the above discussion, the following scenario appears most consistent with the totality of the observational evidence: Shakti and Shiva represent two of those early, massive progenitors that coalseced at high redshift (perhaps T > ∼ 12 Gyr ago), perhaps the last event form the proto-Galaxy, before disk formation commenced.These progenitors enriched more rapidly, and were likely even more dense and massive than GSE (and also LMC and SMC), at least back then at the time of accretion.This conjecture is supported by the combination of their unconventional properties.
To conclude the discussion, we ask whether Shakti and Shiva have any population of globular clusters dynamically "associated" with them?This is an interesting question because recent studies have shown that it is common for massive mergers to bring in populations of such objects (e.g., Massari et al. 2019;Myeong et al. 2019;Malhan et al. 2022).To this end, it seems reasonable to link those Milky Way globular clusters that possess orbits and metallicities similar to these substructures.For the globular clusters, we obtain their orbits and metallicities from Malhan et al. (2022) and Harris (2010), respectively.Next, we consider only those globular clusters that lie within the same (L z , E) polygon and [M/H] range as that described to select the substructure stars.Shakti is found to be linked with only n = 1 cluster (namely, NGC 6235) and Shiva with n = 7 clusters (namely, NGC 5986, NGC 6121/M 4, NGC 6218/M 12, NGC 6287, NGC 6333/M 9, NGC 6397,NGC 6809/M 55).

CONCLUSION AND DISCUSSION
We have identified two new stellar structures in the Milky Way, seen as prominent, confined overdensities in the (L z , E) distribution of bright (G < 16), metal-poor (−2.5 < [M/H] < −1.0) stars.This was enabled by using the novel XGBoost dataset (Rix et al. 2022;Andrae et al. 2023) that provides robust and precise metallicity estimates (based on Gaia DR3's BP/RP spectra) that can be combined with 6D phase-space measurements (from Gaia DR3 itself) for millions of stars; practically for the highest number of stars ever.More detailed chemical abundance patterns were enabled by SDSS DR17.
Shakti and Shiva populations possess an unconventional combination of orbital and abundance properties that have not been observed previously.This is because it has become possible only recently to combine Gaia and SDSS DR17 datasets to disentangle different Galactic populations (e.g., based on the stellar orbits, [M/H] and [Al/Fe]), even in the inner Galaxy and at old ages, plausibly unearthing ancient (T > ∼ 12.5 Gyr old) populations that formed the proto-Galaxy (Belokurov & Kravtsov 2022;Rix et al. 2022).There may be a preponderance of observational clues that argue for Shakti and Shiva representing two of early, massive (M ⋆ ∼ 10 7 M ⊙ ) -or at least dense -fragments that coalesced at high redshifts (perhaps z > ∼ 5 or T > ∼ 12 Gyr ago) to form the proto-Galaxy.They are less tightly bound than Aurora/POH and therefore have a better chance of appearing as distinct overdensities in orbit-space.
However, there are important aspects in which Shakti and Shiva match expectations for ancient stars that gradually picked up angular momentum by resonant orbit trapping caused by the Galactic bar; a scenario proposed recently by (Dillamore et al. 2023).This scenario predicts ridges in the E − L z plane towards the prograde direction.And for the observationally inferred bar speed, these ridges should lie at orbital frequencies (or binding energies E) broadly similar to those of Shakti and Shiva Ṡuch orbital trapping should affect stars of all metallicities, and we see a fairly broad metallicity distribution in both substructures.However, quantitatively, this scenario does not match our observations: Shakti's compact distribution in L z does not at all look like an extended ridge as observed in the Dillamore et al. (2023) study.Nor does it explain the abundance pattern mismatch between Shakti and GSE, which is the presumed low-L z precursor in this scenario.Whether these inconsistencies, which argue against a spun-up-by-bar scenario, can be reconciled remains to be seen.
In closing, we put Shakti and Shiva in the dynamical context of the various other sub structures detected in the past.Fig. 7 shows different substructures of the Milky Way4 , most of which have been discovered using the Gaia dataset, and their detailed chemical abundances are generally studied using the SDSS, LAM-OST or GALAH surveys.Excitingly, Gaia's excellent all-sky 5D astrometry and metallicities are now being complemented by SDSS-V (Kollmeier et al. 2017), in particular also for Shakti and Shiva.In the near future, these spectroscopic efforts will be complemented by measurements from the upcoming WEAVE (Dalton et al. 2014) and4MOST (de Jong et al. 2019) surveys that will map the northern and southern hemispheres, respectively.Furthermore, observations by the Rubin Observatory (fka LSST, Ivezić et al. 2019) will eventually provide a means to determine precise distances of the halo stars.These surveys together will provide 6D phase-space measurements and abundances of stars out to > ∼ 50 kpc, thereby allowing us to create a highresolution chemo-dynamical map of different Galactic populations, giving us once in a lifetime opportunity to unravel the complete formation history of our Galaxy.
Malhan & Rix ulation mass per giant star (at a given absolute magnitude M G0 < M G0,max ).It is this parameter whose estimate relies on the above assumption."M ⋆ /giant in sample" is the stellar population mass per giant star (at the absolute magnitude M G0 < M G0,max ).Below, we detail our procedure to compute the stellar mass.
A.1.Shakti First, we compute the parameter λ sample .The fact that the sample is incomplete can be realised from Fig. 8 which shows that the distribution of stars in the Θ space is not uniform.This is because we can only observe giant stars in our quadrant of the Galaxy.However, based on our above assumption for a phase-mixed population, the expectation is that this Θ distribution would have been uniform, had we been able to observe all the Shakti stars in the entire Galaxy.Therefore, we require to estimate the mean occupation number N , which is the average number of stars in a given Θ voxel.For this, we manually select the Θ sub-space where the distribution appears to be most complete (as shown in Fig. 8): θ R ∈ [0.8, 3.1], θ ϕ ∈ [3.0, 3.8], θ z ∈ [3.8, 5.4] and set the voxel bin size as θ bin = 0.5 radians.The resulting cumulative distribution of Shakti stars is shown in Fig. 8. Upon comparing this curve with the Poisson curves, we find the average stellar occupancy in a given voxel as N = 2.6.This value needs to be multiplied with the total number of voxels, i.e., 2π/θ bin = 13 3 = 2197, which gives the total number of stars after correcting for sample incompleteness as N rmcorrected = 5712.This implies λ sample = N corrected /N observed = 5712/1719 = 3.3.
Second, we compute the parameter "M ⋆ /giant in sample" following the same approach described in Rix et al. (2022).For this, we use the globular cluster NGC 362, which serves as a template for an old, [M/H]= −1.25 population.Beyond NGC 362's half-mass radius of 1.2 ′ (Baumgardt & Hilker 2018), we find about 95 giants with M G0 < 0 (this magnitude limit is set by the faint limit of the Shakti stars after the above Θ sub-space selection).Given NGC 362's total mass of 252, 000 M ⊙ (Baumgardt & Hilker 2018), this implies one giant (at M G0 < 0) per 1300 M ⊙ of the total stellar population mass.
Finally, the above λ sample (= 3.3) and M ⋆ /giant in sample(= 1300 M ⊙ ) values are substituted in equation A1 to obtain the stellar mass as M total ⋆ ≈ 0.7 × 10 7 M ⊙ .

Figure 2 .
Figure 2. Overview of the overall dataset used in this analysis, without any kinematic or [M/H]XP cuts.Panel (a): Distribution of stars in the Galactic coordinates.The legend on the top-left mentions the total number of stars in our dataset.Panel (b): Distribution of stars in the Galactocentric x − y plane.In this Cartesian system, the Galactic center (denoted by "+") lies at the origin and the Sun (denoted by a yellow dot) is at (x, y, z) = (−8.122,0, 0) kpc.Panel (c): Same as panel b, but in the Galactic R − z plane.Panel (d): Velocity behaviour of the sources in spherical polar coordinates, namely radial vr and azimuthal v ϕ .Panel (e): Color-magnitude distribution (CMD), where each star is coloured by its [M/H]XP value.Note that the y-axis is the absolute magnitude.Panel (f ): Effective temperature (T eff,XP ) vs surface gravity (log(g)XP).Panel (g): Metallicity distribution function (MDF).
Fig. 3 shows the (L z , E) orbit distribution of these stars in different [M/H] XP slices, i.e. n L z , E | [M/H] XP .Foremost, this Figure illustrates the known fact that the stellar orbit distribution depends dramatically on their metallicity, with the Galactic disk(s) dominating at [M/H] XP > −1 (see panels band c).It also shows clearly some of the well-established low-metallicity 'sub-structures', such as the metal-rich in-situ halo/Splash(Bonaca et al. 2017;Di Matteo et al. 2019;Belokurov et al. 2020, visible in the panel c), GSE merger debris(Belokurov et al. 2018;Helmi et al. 2018, visible in the panels d and e) and POH(Belokurov & Kravtsov 2022; Rix et al. 2022, visible  in the panels d and e).Panels d and e of Fig. 3 (which together correspond to the metal-poor regime −2.0 ≤ [M/H] XP < −1.0) Figure 3. (Lz, E) distribution of stars in different slices of metallicity [M/H]XP.Panel (a): (Lz, E) distribution of the entire dataset is shown as a scatter plot.Panel (b) -(f ): (Lz, E) distribution in different [M/H]XP slices is shown as 2D histogram plots.[M/H]XP decreases from panel (b) to (f).In every panel, the top text mentions the [M/H]XP slice range and the top-left legend mentions the number of stars within this [M/H]XP slice.In a given panel, the 2D histogram is of (n bin × n bin ) bins, where n bin = min(100, 0.5 √ n data ) and n data is the number of stars within that [M/H]XP slice.The labelled substructures correspond to those that we visually identify as overdensities.To select the Shakti population, we use the base data and restrict ourselves in the metallicity range −2.5 ≤ [M/H] XP < −1.0 and then retain only those stars that lie inside a 2D (L z , E) polygon described by four vertices as [−1400, −1.66], [−800, −1.66], [−600, −1.72], [−1200, −1.72]; the units are in [ kpc km s −1 , ×10 5 km 2 s −2 ].This polygon best describes this population.The [M/H] range is in concordance with the [M/H] slices in which we identify Shakti in Fig. 3, and we especially probe the low-[M/H] range to gather even the very metal-poor stars.This selection results in N sample = 1719 stars.

Figure 5 .
Figure 5. Chemo-dynamical properties of the Shiva sub-structure.Same as Fig. 4, but for the Shiva stars.Next, we find Shiva's detailed chemical abundance properties in an analogous fashion to that of Shakti.In this case, with cross-matching of the Gaia-identified Shiva sample members with APOGEE, we find n = 77 stars possessing [Fe/H], [Mg/Fe] and [Al/Fe] measurements, and of these, only 68 stars possess [Mg/Mn] measurements.We find it is relatively metal-poor, with a median of [Fe/H]≈ −1.25, it is quite rich in Mg, Mn and Al, as listed in Table 1.Similarly to Shakti, Shiva's population also has 50% of stars with [Al/Fe]> 0 and shows a remarkably tight relation in [Fe/H] -[Al/Fe].
whose stars are distributed in the inner D gal ∼ 3.5 kpc of the Galaxy and they possess [M/H]< −1.0.

Figure 7 .
Figure 7. Milky Way's action-energy (J, E) space showing stars from different substructures.Each data point corresponds to a Gaia star.The left panel shows the (Lz, E) distribution and the right panel shows the 'projected action' space, where the horizontal axis is J ϕ /Jtot and the vertical axis is (Jz-JR)/Jtot with Jtot = JR + Jz + |J ϕ |.

Table 1 .
Orbital and chemical properties of Shakti and Shiva stellar populations.

Table 2 .
Arguments on the nature and origin of Shakti and Shiva.