The Poor Old Heart of the Milky Way

The initial inner Galaxy sample was devised via the following ADQL query:

• select ∗ from gaiadr3.gaia_source
• where
• parallax < 100.∗power(10.,0.2∗(0.9 - (phot_g_mean_mag---1.5∗(bp_rp-1.))))
• and parallax < 1.
• and abs(b)<30 and (l<30 or l > 330.)
• and bp_rp between 1.0 and 3.5
• and phot_bp_mean_mag < 15.5
• and has_xp_continuous='true'

The first condition on the parallax selects stars whose absolute magnitude is more luminous than M_G = 0.9, aimed at selecting giant stars at luminosities that just include the RC; the -1.5*(bp_rp-1.) term reflects an approximately dereddened magnitude. The query is articulated in such a way that zero or negative parallaxes are accommodated. The second parallax condition only excludes the immediate "foreground" at D < 1 kpc. The conditions on l and b select a rectangular region of 30° around the GC. The lower limit on bp_rp eliminates most luminous hot stars for which [M/H] estimates are not feasible (but not all in the presence of reddening). The condition on phot_bp_mean_mag is designed to select objects bright enough in the blue part (BP) of the XP spectra for robust [M/H] estimates; it still allows one to select RC stars at the distance of the GC in the presence of moderate reddening. The final line in the selection ensures that continuous XP spectra are available. This query yields 2.1 million objects.

While about 75% of this sample has RVS velocities in Gaia DR3, only a small fraction of them have metallicities or abundances derived from the RVS spectra (Recio-Blanco et al. 2022). This leaves the route of estimating [M/H] from the low-resolution XP spectra of these objects. Consequently, our objective is to train a machine-learning algorithm that can precisely, accurately, and robustly predict [M/H] from Gaia DR3 data.

Estimating [M/H] based on machine learning involves four aspects: (1) it involves the choice of the training sample, for which we adopt SDSS DR17 (APOGEE; Abdurro'uf et al. 2022) since its abundances are well validated, it covers the inner Galaxy, and contains mainly giants, many with high extinctions; (2) it involves the choice of the data features on which to train the prediction; these can be derived from the XP spectra but may also entail, as we discuss below, external data available across the sample, such as near-infrared photometry from AllWISE (Cutri et al. 2021); (3) it involves the choice of the machine-learning algorithm, for which we adopt an extreme gradient boosting algorithm (Chen & Guestrin 2016, hereafter XGBoost); XGBoost is one of the best algorithms currently available, being straightforward and computationally inexpensive to train, and it can outperform other methods such as deep learning (e.g., Grinsztajn et al. 2022); and (4) finally we need to validate our [M/H] prediction, which then speaks to the quality of the data, the training sample, and the algorithm.

The choice of the correct data features to train and predict [M/H] from the XP spectra is not trivial. One obvious option is to directly use the 2 × 55 modified Gauss–Hermite coefficients that describe XP spectra. But our numerical experimentation implies that other or additional data features may yield more precise and robust results with a machine-learning approach, ¹² for several pragmatic reasons. Much of the information about narrow metal line features is contained in high-order coefficients, which are often noisy. We have prior information on which parts of the spectrum are highly diagnostic in [M/H] estimates and which ones are not; we need not ask an algorithm to learn this. Any estimate of [M/H] requires an implicit estimate of the effective temperature, T_eff, which is covariant with [M/H]; and T_eff is also highly covariant with reddening. Therefore, providing additional information that helps break that degeneracy is precious. The near-infrared photometry from AllWISE (Cutri et al. 2021) exists across the sky and proves powerful in this respect.

In practice, we follow the approach of Montegriffo et al. (2022) to calculate synthetic photometry from the XP spectra in a wide range of mostly narrowband filters, in particular those filters with established utility in identifying metal-poor stars: Strömgren filters (Strömgren 1966), JPAS+ filters (Marín-Franch et al. 2012), and Pristine H&K filters (Starkenburg et al.2017). On this basis, we use XGBoost to train the prediction of [M/H] using the entire SDSS DR17 APOGEE of giants with G_BP < 15.5 (∼230,000 stars). Details are given in the Appendix.

Training on a particular data set and using a model that is discriminative, rather than generative, raises several issues. First, we (obviously) tie our results to the metallicity scale of SDSS DR17. Second, estimates of [M/H] may be drawn into the support of the training set. In particular, XGBoost will not extrapolate [M/H] significantly beyond the APOGEE DR17 range, as it is a tree-based method that segments the feature space and assigns a mean label to each such segment. Finally, generative data-driven spectral models, such as the Cannon (Ness et al. 2015), will presumably degrade more gracefully toward low signal-to-noise (S/N). Here we have addressed the last point implicitly by restricting the sample to G_BP < 15.5.

This results in [M/H] estimates for about 2 million giants toward the inner Galaxy, 1.58 million of which have RVS velocities, and 1.25 million have both RVS velocities and ϖ/σ_ϖ ≥ 5, which is our minimal condition for an orbit estimate. For subsequent analysis, we have eliminated 2% of the sample as reddened hot stars, which can be readily recognized by their position in the m₁–bp_rp color plane, where m₁ ≡ v − 2b + y is the metallicity-sensitive Strömgren filter combination (see Strömgren 1966, for details). Specifically, we eliminate sources with m₁ + 0.16 (bp_rp − 1) < 0 and bp_rp < 2.7.

We can explore and validate the precision and accuracy of these [M/H] estimates through various comparisons with analogous estimates from (higher-resolution) spectroscopic surveys. For this validation, we consider SDSS APOGEE (Abdurro'uf et al. 2022), the Large Sky Area Multi-Object Fibre Spectroscopic Telescope (LAMOST; here Xiang et al. 2019), Data Release 3 (DR3) of GALactic Archaeology with HERMES (GALAH; Buder et al. 2021), and Gaia's GSP-Spec (Recio-Blanco et al. 2022).

We start by considering the ∼17,000 stars in common between APOGEE and our sample, shown in Figure 1 (top panel). This panel of Figure 1 illustrates that, for this sample, XGBoost provides a remarkably precise (≲ 0.1 in the median), accurate, and robust [M/H] estimate across −2.5 < [M/H] < 0.5. It should be noted that XGBoost was trained on the full all-sky APOGEE sample, so this is a validation on a small subset of the full training set, not a pure cross-validation. The accuracy of the [M/H] prediction (i.e., lack of systematic offset) is therefore by construction. It is particularly remarkable that basically all stars predicted to be metal-poor from [M/H]_XP are indeed metal-poor according to [M/H]_APOGEE: a selection of, for example, metal-poor objects by [M/H]_XP appears to be nearly pure. The bottom panel of Figure 1 shows that these [M/H] estimates remain unbiased, at a level of a few percent, even in the presence of substantial dust extinction, A_K = 0.3, corresponding to about three magnitudes of A_V .

Zoom In Zoom Out Reset image size

Appendix A.2 details a true cross-validation of our [M/H]_XP estimates against three external data sets: LAMOST, GALAH, and GSP-Spec. Broadly speaking, the cross-validation affirms that the [M/H]_XP estimates are precise and robust (enabling pure [M/H]-based sample selection). Our estimates have a rms difference of ≈0.1 dex compared to these data sets, with minimal bias. Most importantly, the performance remains unbiased toward the metal-poor end, affirming our confidence in selecting metal-poor stars based on XP metallicities. We refer to Appendix A.2 for more detailed validation of the [M/H] estimates.

For sample members with suitable 6D phase-space information, we calculate their orbits. The sky positions and proper motions are available for all of them, radial velocities and good parallax measurements only for a good fraction of the sample. In practice, we take all stars that have RVS velocities with δ v_RVS < 5 km s⁻¹ and ϖ/σ_ϖ > 5, about 1.25 million objects across all [M/H].

To compute orbits and compute orbital properties, we use a four-component Milky Way mass model consisting of spherical Hernquist nucleus and bulge components, an (approximate) exponential disk component (Smith et al. 2015), and a spherical Navarro–Frenk–White halo component. We adopt a radial scale length h_R = 2.6 kpc and scale height h_z = 300 pc (Bland-Hawthorn & Gerhard 2016), and fit for the masses and scale radii of the nucleus, bulge, and halo components using the same compilation of Milky Way enclosed mass measurements as used to define the MilkyWayPotential in gala (Price-Whelan 2017; Price-Whelan et al. 2020), with the additional constraint of having a circular velocity at the solar position v_c (R₀) = 229 km s⁻¹ (Eilers et al. 2019). We compute actions using the "Stäckel Fudge" (Binney 2012; Sanders 2012) as implemented in galpy (Bovy 2015). To compute the orbital eccentricity, pericenter, and apocenter values, we numerically integrate the orbits with a time step of 0.5 Myr across four times the radial orbit period (estimated by the computed orbital frequencies from the action solver).

The [M/H] estimates for 1.5 million stars toward the GC, along with orbits for the RVS subsample of 1.25 million stars, are available as supplementary material to this article and are hosted online ¹³ (Andrae et al. 2022a).

Figure 2 compares the on-sky density distribution of metal-rich stars with [M/H] > −0.4 (left panel) to metal-poor ones with [M/H] < −1.5 (middle panel). The metal-poor stars exhibit a striking central concentration relative to the metal-rich sample, suggesting a spheroidal population in the inner Galaxy. The projected sky density of metal-poor sample members in Figure 2 (middle panel) is dramatically altered by the high dust extinction toward the GC (right panel). The fine-scale structure in the projected [M/H] < −1.5 stellar distribution can mostly be explained by dust extinction, as the right-hand panel of Figure 2 shows. This is expected as we require G_BP < 15.5 for sample selection.

Zoom In Zoom Out Reset image size

We do not have precise parallaxes for the entire sample, which complicates the analysis of the spatial structure. For the overall sample, the median parallax precision is ϖ/σ_ϖ ∼ 12, with 10% of the sample having ϖ/σ_ϖ ≲ 4, and 1% of the sample with negative parallaxes. For the metal-poor sample members with [M/H] < −1.5, which turn out to be more distant than the sample mean; the median precision is only ϖ/σ_ϖ ∼ 5. Therefore, the most robust way to present the distance distribution of objects in a directional cone may be to consider their parallaxes. Given that the objects of particular interest are at ∼8 kpc distance, it matters that we apply a Gaia parallax zero-point correction, for which we adopt 0.02 mas (Lindegren et al. 2021).

Figure 3 shows ${n}_{* }\left(\varpi \,| \,[{\rm{M}}/{\rm{H}}]\right)$ in five [M/H] bins, with the minimal distance (maximal parallax) of 1 kpc at the right edge and the parallax expected for the GC indicated by the dotted line (1/8.2 kpc, or 0.12 mas; Gravity Collaboration et al. 2021). This figure shows that the distance distribution of giants toward the GC depends dramatically on [M/H]. Metal-rich stars ([M/H] > −0.4) are distributed throughout the disk with heliocentric distances D_⊙ = 1−5 kpc (or 0.2–1.0 mas). In contrast, metal-poor populations become increasingly more concentrated toward the GC. Indeed, the most metal-poor stars are predominantly clustered around the GC's parallax, revealing that there is a centrally concentrated very metal-poor population in the Milky Way.

Zoom In Zoom Out Reset image size

Figures 2 and 3 qualitatively confirm the picture of a metal-poor population that is very concentrated toward the GC. The fact that there are almost no [M/H] < −1.5 stars at very low latitudes shows that almost all lie behind the strong dust extinction 2–3 kpc inward of the Sun. Rigorous modeling of these stars' spatial distribution of these stars density distribution (e.g., Rix et al. 2021) is beyond the scope of this paper, with the severe dust extinction and crowding effects presenting formidable challenges. However, we can model the parallax distribution (Figure 3) in the central regions less affected by dust to quantify the radial extent of this metal-poor population.

We assume that the 3D distribution is a spherical Gaussian as a function of Galactocentric radius R ${\rho }_{\mathrm{Gauss}}(R,{\sigma }_{{R}_{\mathrm{GC}}})$, with a width of ${\sigma }_{{R}_{\mathrm{GC}}}$. As long as we can see to the far side of the GC (see Figure 3), we should expect the parallax (ϖ) distribution in a cone of dΩ in the direction (l,b) to be

$$\begin{eqnarray}&&n(\varpi )=d{\rm{\Omega }}\,d\varpi \,{\varpi }^{-4}\,{\rho }_{\mathrm{Gauss}}\left(R(\varpi ,l,b),{\sigma }_{{{\rm{R}}}_{\mathrm{GC}}}\right).\end{eqnarray} \tag{ 1 }$$

Phrasing this model in terms of the ϖ allows the model to deal with vanishing (or even negative) observed parallaxes. We find that the maximum likelihood of the observed parallaxes (and their uncertainties) for the [M/H] < −1.5 stars implies an extent of ${\sigma }_{{R}_{\mathrm{GC}}}\sim 2.7\,\mathrm{kpc}$. Interestingly, this radius corresponds to an angle of 18° at the distance of the GC, which seems very plausible in light of the projected density distribution (Figure 2; middle panel).

We now turn our attention to the [M/H] distribution, n_*([M/H]), in the inner Galaxy, which is illustrated in Figure 4. The first remarkable feature of the distribution is the sheer number of objects, most apparent from the cumulative distribution in Figure 4: there are >4000 stars with an estimated [M/H] < −2, ∼18,000 with [M/H] < −1.5, and almost 100,000 stars with [M/H] < −1, i.e., below the abundance floor of the old, α-enhanced disk. This is an order of magnitude more objects than previously published metal-poor samples ([M/H] < −1) of the inner Galaxy (e.g., Arentsen et al. 2020a, 2020b).

Zoom In Zoom Out Reset image size

We have already presented two arguments that the low-[M/H] portion of n_*([M/H]) is not severely contaminated by spurious [M/H] estimates from some of the far more numerous stars of higher (true) metallicity: the validation with APOGEE and the parallax, or distance, distribution, which is most peaked for the most metal-poor objects. The [M/H] distribution itself provides additional evidence. It rises very steeply from [M/H] = −2.5 to −2.2. The distribution then follows a power law of slope $d(\mathrm{log}{n}_{* })/d[{\rm{M}}/{\rm{H}}]\approx 1$ over a wide range of metallicities to [M/H] ≈ −0.9, where it steepens toward a peak at [M/H] = −0.6, beyond which n_*([M/H]) starts dropping toward the highest metallicities present, [M/H] ≈ 0.5.

The bottom panel of Figure 4 shows the metallicity distribution with stars of an estimated Galactocentric distance ${\hat{R}}_{\mathrm{GC}}\lt 4\,\mathrm{kpc}$ in green, and the stars of the intervening disk in gray. As giants, in particular red clump stars, are of comparable luminosity and color over a wide range of metallicities, the differential selection effects across different [M/H] should be modest. Hence, the [M/H] distribution of stars R_GC < 4 kpc in Figure 4 should be relatively unbiased. The [M/H] distribution of course varies with position in the Galaxy (both ∣z∣ and R_GC), as Figure 4 illustrates.

The decrease in n_*([M/H]) for stars at the metal-rich end ([M/H] > 0) with R_GC < 4 kpc (green) may seem surprising, because we know from APOGEE (e.g., Eilers et al. 2022) that the inner Galaxy is teeming with high-metallicity stars. This decrease in n_*([M/H]) is straightforwardly explained by Figure 2, which shows that in the inner Galaxy we are almost completely lacking sample members with ∣Z∣ < 700 pc; any thinner population will inevitably be absent from our sample due to dust extinction. The APOGEE sample has a different MDF in part because its Two Micron All Sky Survey–selected stars better sample the low-∣Z∣, high-[M/H] population.

3.2.1. The Mass of the Central Proto-galactic Stellar Population

We now proceed to see what the observed metallicity distribution might tell us about the early central enrichment history.

The bulk of the metal-poor population, below the minimal [M/H] of the old disk at [M/H] ∼ −1, follows $d(\mathrm{log}{n}_{* })/d[{\rm{M}}/{\rm{H}}]\approx 1$: the cumulative number of stars N_*(<Z) below a (linear) metallicity Z grows linearly with Z. This seems to be a natural slope for n_*([M/H]) in the early phases of a self-enriching, in situ system, as shown by the following argument, which builds on the models of Weinberg et al. (2017). Under fairly general conditions, the evolution of the interstellar medium's (ISM's) metallicity in a simple, fully mixed system can be described by

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{Z}}{{dt}}={y}_{Z}\,{\dot{M}}_{* }-(1+\eta -r){\dot{M}}_{* }Z,\end{eqnarray} \tag{ 2 }$$

where M_Z is the mass in metals, y_Z is the initial-mass-function-averaged yield, ${\dot{M}}_{* }$ is the star formation rate, Z ≡ M_Z /M_g is the metal mass fraction of the star-forming gas, $\eta \equiv {\dot{M}}_{\mathrm{out}}/{\dot{M}}_{* }$ is the mass loading factor of a gaseous outflow, and r ≈ 0.4 represents the recycling of metals that are formed into stars but quickly returned when the stars die. Equation (2) treats enrichment as instantaneous, which should be a good approximation at low metallicities where core-collapse supernovae and massive star winds are the dominant sources. It also assumes that outflows have the same metallicity as the ambient ISM; if outflows instead consist of a fraction, f, of the supernova ejecta plus entrained ISM, then the effect is to reduce y_Z by a factor (1 − f) so that it reflects only metals retained by the ISM. Importantly for the case at hand, if Z is well below the equilibrium metallicity Z_eq ≡ y_Z /(1 + η − r) then the entire second term of Equation (2) can be neglected.

In this low-metallicity limit, Equation (2) is simply ${\dot{M}}_{Z}={y}_{Z}{\dot{M}}_{* }$, and for a metallicity-independent yield, its time integral implies M_Z = y_Z M_* and thus Z = y_Z M_*/M_g . The log-slope of the MDF is

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{log}{n}_{* }}{d[{\rm{M}}/{\rm{H}}]}\,=\,\displaystyle \frac{d\mathrm{log}{M}_{* }}{d\mathrm{log}Z}\,=\,\displaystyle \frac{Z}{{M}_{* }}\cdot \displaystyle \frac{{\dot{M}}_{* }}{\dot{Z}}\,=\,\displaystyle \frac{{y}_{Z}}{\dot{Z}{\tau }_{* }},\end{eqnarray} \tag{ 3 }$$

where the last equality introduces the star formation efficiency (SFE) timescale ${\tau }_{* }\equiv {M}_{g}/{\dot{M}}_{* }$. Using

$$\begin{eqnarray}&&\dot{Z}=\displaystyle \frac{d}{{dt}}({M}_{Z}/{M}_{g})={y}_{Z}\displaystyle \frac{{\dot{M}}_{* }}{{M}_{g}}-{y}_{Z}\displaystyle \frac{{M}_{* }{\dot{M}}_{g}}{{M}_{g}^{2}}\end{eqnarray} \tag{ 4 }$$

gives the end result:

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{log}{n}_{* }}{d[{\rm{M}}/{\rm{H}}]}\,=\,{\left(\displaystyle \frac{\dot{Z}{\tau }_{* }}{{y}_{Z}}\right)}^{-1}\,=\,{\left[1-\left(\displaystyle \frac{{M}_{* }}{{M}_{g}}\right)\left(\displaystyle \frac{{\tau }_{* }{\dot{M}}_{g}}{{M}_{g}}\right)\right]}^{-1}.\end{eqnarray} \tag{ 5 }$$

The combination $({\tau }_{* }{\dot{M}}_{g}/{M}_{g})$ corresponds to the fractional change of the gas reservoir mass over one SFE timescale.

Equation (5) shows that an MDF log-slope ≈ 1 is a generic result in the low-metallicity regime (Z ≪ Z_eq), arising whenever the gas mass is constant (${\dot{M}}_{g}=0$) or more generally when the star-to-gas ratio, M_*/M_g , is small enough that the second term in Equation (4) can be neglected. Figure 4 exhibits this generic slope over the range −2 < [M/H] < −1. The steeper slope at [M/H] < −2 could plausibly arise because the gas reservoir is small but rapidly growing, with ${\tau }_{* }{\dot{M}}_{g}/{M}_{g}\gt 1$. The steeper slope at −1 < [M/H] < −0.6 may reflect rapid gas accretion (high ${\dot{M}}_{g}/{M}_{g}$) coinciding with the onset of the old α-enhanced disk (e.g., Belokurov & Kravtsov 2022; Conroy et al.2022; Xiang & Rix 2022). At still higher metallicities our neglect of sink terms in Equation (2) becomes a poor approximation.

We now turn to the orbit distribution of the sample at hand (see Section 2.4), which can tell us which of the metal-poor stars selected in the inner Galaxy remain confined to the central regions of the Milky Way. The orbit distribution can tell us which stars are just passers-through near their pericenter but on orbits that take them into the outer halo. In particular, we might expect to see the pericenter members of accreted halo components like GSE, with stars on highly eccentric orbits that take them to R_apo > 10 kpc. We can also delineate how early—or how metal-poor, if we use abundances to estimate relative stellar ages— the population develops net rotation and how rapidly the kinematics change with [M/H] toward near-circular motion in the plane of the Galaxy (see Arentsen et al. 2020a; Belokurov & Kravtsov 2022; Conroy et al. 2022).

We start by considering the orbit distribution, n_*(e, R_apo), in terms of the orbital eccentricity, e, and the apocenter, R_apo. We initially focus on the [M/H] range −1.9 < [M/H] < −1.2 since this encompasses the bulk of known GSE stars. This distribution of 28,000 stars is shown in Figure 5. The density in the (e, R_apo) plane is affected and limited by two experimental aspects. First, our initial Gaia query selects almost exclusively stars with R_GC < 8 kpc, and hence all stars whose pericenter distance exceeds 8 kpc are excluded; this boundary is indicated by the thick gray line. Except at e < 0.2 the n_*(e, R_apo) distribution falls off toward large R_apo well before this limit. There may be two reasons for the seeming dearth of stars with R_apo < 2 kpc near the center. First, severe dust extinction simply obscures much of the central kiloparsec (Figure 2). Second, any uncertainties in the 6D phase-space coordinates, especially uncertainties in the distance, are far more likely to increase the inferred R_apo than decrease it, as it is a positive definite quantity. These aspects deserve careful modeling, which is beyond the scope of this paper.

Zoom In Zoom Out Reset image size

The distribution in Figure 5 can be characterized by two regimes: most stars form a nearly flat distribution in eccentricity from e = 0.1 to e = 0.8, but with a narrow range in R_apo, with 〈R_apo〉 ≈ 4 kpc. This implies that most stars form an approximately isotropic distribution that stays confined to the inner Galaxy and hence will not appear in surveys of halo stars that focused on larger R_GC and off the Galactic plane. The second regime is that of e ≥ 0.75 and R_apo ≳ 10 kpc. The orbits of these stars fully match the expectations of the pericenter members of the GSE stars, the accreted component that dominates the halo population between 10 and 30 kpc (Belokurov et al. 2018a; Helmi et al. 2018; Naidu et al. 2021). Indeed, recent work by Belokurov et al. (2022) has shown quite clearly that there are likely GSE members with pericenters within 3 kpc. There is tantalizing but inconclusive evidence from Figure 5 that the distribution of GSE stars might extend to small R_apo and lower e. Figure 5 indicates an excess of stars near (e, R_apo) ≈ (0.7, 4 kpc) that appears to continue down from the accreted GSE sequence. Here it would become important to disentangle GSE from other proposed accreted structures in the inner Galaxy (e.g., Kruijssen et al. 2019; Horta et al. 2021; see Section 3.5).

We now ask to what extent different metal-poor stars in the inner Galaxy have net angular momentum. In general, more metal-rich stars are generally expected and observed to have a more coherent sense of angular momentum. We explore this in the two panels of Figure 6, which show how close, or far, stars at a given [M/H] and R_apo are from being an ensemble of stars on circular, in-plane orbits. This is quantified via the ratio of the angular momentum, which is the azimuthal action J_ϕ , to the total action ${J}_{\mathrm{tot}}=\sqrt{{J}_{\phi }^{2}+{J}_{R}^{2}+{J}_{z}^{2}}$. The top panel shows that the increase of $\overline{{J}_{\phi }/{J}_{\mathrm{tot}}}$ with [M/H] depends somewhat on the orbit size: At a given [M/H] orbits with smaller R_apo are farther away from the coplanar quasi-circular orbits than at R_apo ∼ 5 kpc.

Zoom In Zoom Out Reset image size

The bottom panel of Figure 6 shows $\overline{{J}_{\phi }/{J}_{\mathrm{tot}}}([{\rm{M}}/{\rm{H}}])$, marginalized over 3 < R_apo/kpc < 7: starting already at [M/H] < −2 there is some, albeit small, net $\overline{{J}_{\phi }/{J}_{\mathrm{tot}}}([{\rm{M}}/{\rm{H}}])$, which increases gradually to ∼0.3 at [M/H] = −1.3. Among more metal-rich populations, $\overline{{J}_{\phi }/{J}_{\mathrm{tot}}}([{\rm{M}}/{\rm{H}}])$ steeply rise with rising [M/H] to [M/H] ∼ −0.7, after which it levels off in the regime dominated by cold rotation. This analysis implies that some modest net rotation is present in the population well below [M/H] < −1.5, affirming and extending the findings of Arentsen et al. (2020a), Belokurov & Kravtsov (2022), and Conroy et al.(2022).

We now enlist element abundances for our sample stars that overlap with APOGEE to explore which parts of the sample appear accreted from a chemical enrichment perspective (as opposed to proto-Galactic, or in situ). Such abundance-based arguments for the origin of stars are based on the idea that enrichment works differently in (sub)halos of different mass (e.g., Hawkins et al. 2015; Das et al. 2020; Horta et al. 2021; Belokurov & Kravtsov 2022): the production of α elements and aluminum aided by high star formation intensities in the deeper potential wells of the proto-Galaxy, compared to the shallower potential of lower-mass satellite galaxies. However, these abundance diagnostics of the stars' formation environment may be less discriminating at lower metallicities (e.g., Conroy et al. 2022).

We remove stars belonging to globular clusters and apply a few basic quality cuts to the APOGEE abundance data based on their uncertainties and flags. We plot the resulting chemical distributions, [X/Fe] versus [Fe/H], with ${\rm{X}}=\left\{\alpha ,\mathrm{Al},\mathrm{Mn}\right\}$, as a function of eccentricity in Figure 7, by median apocenter radius. The presumed boundaries between in situ and accreted are indicated by the yellow dashed lines.

Zoom In Zoom Out Reset image size

Figure 7 shows clearly that at modest eccentricities (left and middle columns), most stars belong to the α-enhanced sequence, implying a proto-Galactic origin. At high eccentricities (right columns), a prominent low-α sequence with large apocenters appears. These stars on eccentric orbits with large apocenters appear to be remnants of the accreted GSE merger based on their chemistry and kinematics (e.g., Mackereth et al. 2019; Bonaca et al. 2020; Naidu et al. 2020; Hasselquist et al. 2021; Horta et al. 2022). Note that among the stars of highly eccentric orbits there is a portion whose abundances point toward proto-Galactic origin (right column). Remarkably, they almost exclusively have apocenters ≲5 kpc, very much like their cousins on less eccentric orbits, but very unlike the chemically identified GSE debris. The population of eccentric stars attributable to the proto-Galaxy has been seen before. It has been dubbed "Aurora" by Belokurov & Kravtsov (2022) and has been explored at even lower metallicities by Conroy et al. (2022). Our sample probes to considerably lower R_GC than either of these works, and reveals the bulk and full extent of the metal-poor proto-Galaxy at the heart of the Milky Way.

It is apparent from Figures 5 and 7 that both orbits and abundances can help to differentiate between presumed proto-Galactic stars and stars accreted later from a lower-mass satellite. A differentiation by chemistry alone appears increasingly difficult at low [M/H] (Conroy et al. 2022). But the combination of both aspects, summarized in Figure 7, can make such a differentiation more convincing.

Horta et al. (2021) selected "accreted" stars in these APOGEE chemical spaces as evidence of a past "Heracles" merger buried deep within the Milky Way, distinguished from GSE on the basis of a bimodal distribution in total orbital energy. This claimed galaxy bears resemblance to the previously proposed Kraken/Koala mergers (Kruijssen et al. 2019; Forbes 2020; Kruijssen et al. 2020), and these names may well refer to the same event. However, Lane et al. (2022) pointed out that the apparent bimodality in orbital energy used by Horta et al. (2021) to select Heracles/Kraken/Koala might be an artifact of the spatial selection function of APOGEE. Furthermore, Heracles/Kraken/Koala and the proto-Galactic Aurora are virtually indistinguishable in all chemical spaces (Horta et al. 2022; Myeong et al. 2022; Naidu et al. 2022).

Taken together, these lines of evidence suggest that the stars attributed to Heracles, Kraken, Koala, and Aurora, along with the numerous metal-poor stars identified here, are different orbit and metallicity regimes of one proto-Galactic stellar population, now forming the metal-poor, ancient heart of the Milky Way. Whereas past studies have mainly identified the eccentric tail of this population that traverses the solar neighborhood or off the Galactic plane (Belokurov & Kravtsov 2022; Conroy et al. 2022; Myeong et al. 2022), we now show more directly that the bulk of this population resides in the bulge itself, mostly within a few kiloparsecs of the GC. Indeed, this matches theoretical predictions from cosmological zoom-in simulations, which predict a spheroidal distribution of proto-Galaxy stars predominantly concentrated around the GC (e.g., El-Badry et al. 2018; Wetzel et al. 2022; Figure 13 in Belokurov & Kravtsov 2022).

The picture presented here appears to obviate the need for an ex situ merger like Heracles/Kraken/Koala, at least based on the orbits and chemistry of field stars. This does not necessarily imply that these mergers did not occur; rather, it emphasizes the—conceptual and practical—ambiguity between in situ and accreted stars during the earliest phases of the Milky Way, which led us to use the collective term proto-Galaxy. This view is supported by cosmological zoom-in simulations, which affirm the difficulty in distinguishing accreted galaxies from the Milky Way proto-Galaxy at early times (e.g., Renaud et al. 2021a, 2021b; Orkney et al. 2022).

We want to derive [M/H] from XP information for objects that still have significant flux in the blue (G_BP < 15.5) yet may be reddened by several magnitudes. To help break any degeneracies between extinction and T_eff, we include near-infrared photometry from AllWISE (Cutri et al. 2021). We adopt SDSS DR17 APOGEE abundances (Abdurro'uf et al. 2022) as a training sample, since it covers the inner disk and contains many giants stars with high extinctions.

For the machine-learning algorithm, we choose the extreme gradient boosting algorithm (Chen & Guestrin 2016, hereafter XGBoost), as it is computationally inexpensive to train and can outperform other algorithms, e.g., deep learning (e.g., Grinsztajn et al. 2022). To improve the [M/H] predictions for our sample (consisting of red giant stars by construction), we also restrict the training sample from SDSS DR17 sample to only giants, via $\mathrm{log}g\lt 3.5$.

This leaves us with a suitable choice of data features, for which we explored a number of options. Our starting point was the set of XP coefficients (Carrasco et al. 2021; De Angeli et al.2022), which we normalized by the Gaia apparent G flux. In Table 1, we quote the test errors for using the different numbers of XP coefficients. The performance is overall very good but depends only weakly on how many coefficients we use. This suggests that also the high-order coefficients contain useful information, and, at least for this application and the APOGEE DR17 sample, a truncation of XP coefficients is neither required nor helpful in our case.

Given the results in Montegriffo et al. (2022), we then attempted to only use photometry synthesized from XP spectra using GaiaXPy ¹⁴ (Ruz-Mieres 2022) and combine this with AllWISE near-infrared photometry. Specifically, we synthesize photometry from the systems of Pristine, StromgrenStd, Jpas, Jplus, SkyMapper, and HstAcswfc. As we can see from Table 1, using only such colors has a similar performance as only using the XP coefficients, despite including AllWISE photometry. However, a substantial improvement is obtained when combining all XP coefficients with these colors.

In the end, we used GaiaXPy to compute the following synthetic photometry and combined it with all 110 XP coefficients in the XGBoost training on SDSS DR17 (APOGEE):

$$\begin{eqnarray*}&&\begin{array}{r}G-{W}_{1}\\ {G}_{\mathrm{BP}}-{G}_{\mathrm{RP}}\\ {G}_{\mathrm{BP}}-{W}_{2}\\ {W}_{1}-{W}_{2}\\ {\mathrm{CaHK}}_{\mathrm{Pristine}}-{W}_{1}\\ {\mathrm{CaHK}}_{\mathrm{Pristine}}-{v}_{\mathrm{StromgrenStd}}\\ {\mathrm{CaHK}}_{\mathrm{Pristine}}-{b}_{\mathrm{StromgrenStd}}\\ {\mathrm{CaHK}}_{\mathrm{Pristine}}-{y}_{\mathrm{StromgrenStd}}\\ {\mathrm{Jplus}}_{g}-{\mathrm{Jplus}}_{i}\\ {\mathrm{Jplus}}_{0395}-{\mathrm{CaHK}}_{\mathrm{Pristine}}\\ {\mathrm{Jplus}}_{0515}-{\mathrm{CaHK}}_{\mathrm{Pristine}}\\ {\mathrm{Jplus}}_{0861}-{\mathrm{CaHK}}_{\mathrm{Pristine}}\\ {G}_{\mathrm{RP}}-{\mathrm{Jplus}}_{i}\\ {v}_{\mathrm{StromgrenStd}}-2\cdot {b}_{\mathrm{StromgrenStd}}+{y}_{\mathrm{StromgrenStd}}\\ {\mathrm{CaHK}}_{\mathrm{Pristine}}-2\cdot {\mathrm{Jplus}}_{0410}+{\mathrm{Jplus}}_{0430}\\ \end{array}\end{eqnarray*}$$

We emphasize that the G_BP − W₂ color has the longest wavelength leverage and therefore is highly affected by extinction, whereas the W₁ − W₂ color is comparatively insensitive to extinction.

As is evident from Table 1, this final configuration has slightly worse performance than the best configuration. However, the median absolute errors are similar, that is, most sources get similar results. The only noteworthy difference is in the rms error, suggesting that the final configuration may have a few more outliers. Still, the results are better than what we can achieve from coefficients alone.

Figure 8 illustrates a cross-validation comparison of our results to other spectroscopic surveys. Such cross-validation is easy to interpret if one can assume that the external data set represents "ground truth." We applied a SNR_G > 25 to the LAMOST sample and required that there were no quality (concern) flags in any of the GALAH spectra. The left panel shows that the [M/H]_XP estimates are accurate for LAMOST, which itself is tied to the APOGEE [M/H] scale, with a slight increase in scatter toward low [M/H]. The comparison with GALAH (middle panel) affirms the precision and purity of the low-[M/H]_XP estimates. Note that GALAH has a different scaling between [M/H] and [Fe/H] than APOGEE and LAMOST, which may explain the offset of the estimates for low-[M/H] stars from the 1-to-1 line. In the parent sample without stringent GALAH quality cuts, there are many stars for which GALAH implies that they are [M/H] < −1, but the [M/H]_XP estimate disagrees. The position of those stars in the dereddened color–magnitude plane suggests that they indeed have [M/H] > −1. In other words, for these cases [M/H]_XP appears to be the more robust estimate. The comparison with GSP-Spec in the right panel of Figure 8 confirms the quality of GSP-Spec's [M/H] estimates. But it also shows that the GSP-Spec values do not cover low [M/H] (by design; see Recio-Blanco et al. 2022). It is possible to get GSP-Spec values for a wider range of metallicities by relaxing quality cuts, but this results in a more contaminated sample with systematic grid effects.

Zoom In Zoom Out Reset image size

Because a star's effective temperature (T_eff) can be very helpful in sample selection and characterization, we also estimate T_eff using the same XP features and XGBoost trained on SDSS DR17 APOGEE T_eff, as illustrated in Figure 9. The XGBoost predictions based on XP spectra and WISE photometry match the APOGEE values within a median ΔT_eff = 32 K, and with a mean difference of only 3 K. Figure 9 compares our [M/H] estimates to APOGEE as a function of Gaia G_BP magnitude, showing the slight degradation toward the faintest sources in the sample.

Zoom In Zoom Out Reset image size

We also compared our [M/H] estimates to the extensive set of photometric estimates for metal-poor stars from Chiti et al. (2021) based on SkyMapper Data Release 2 (DR2/SM2) photometry. We found a mean [M/H] offset of 0.28 dex, and for stars with [M/H]_SM2 ∼ −1.5 a central 68% interval of $p\left(\,{[{\rm{M}}/{\rm{H}}]}_{\mathrm{SM}2}-{[{\rm{M}}/{\rm{H}}]}_{\mathrm{XP}}\,\right)$ of 1 dex. Clearly, the cross-validation of our XP-based [M/H] estimates and SkyMapper DR2 shows far more scatter than our cross-validation with external spectroscopic data sets. This illustrates that for this sample XP-based [M/H] estimates outperform those of SM2.