© The Authors 2025

1 Introduction

A key objective in Galactic astronomy is to construct a cohesive formation history of the Milky Way down to its earliest times. More broadly, we aim to understand the extent to which a galaxy retains its formation history over cosmic timescales and to investigate the physical processes shaping disc galaxies using the Milky Way as an example (i.e. Galactic archaeology). The Milky Way is the perfect cosmic backyard for this because it offers detailed unparalleled, star-by-star 6D kinematics and chemical information.

The rotationally supported disc contains most of the stars in the Galaxy. It is well established that the disc consists of multiple components or populations, distinguished by their chemistry, kinematics, spatial extent, and age (Norris et al. 1985; Gilmore & Reid 1983; Chiba & Beers 2000; Nissen & Schuster 2010; Bovy et al. 2012; Haywood et al. 2013; Hayden et al. 2015). Structurally, there is evidence of thin and thick discs with different scale heights. Chemically, there is a distinct dichotomy and bimodality in the abundance of α-elements relative to iron ([α/Fe]); the Milky Way disc manifests a high-α and low-α sequence at fixed [Fe/H] (Hayden et al. 2015; Imig et al. 2023). Generally, the high-α population is older and has a larger scale height than the low-α population. While the connection between high-α and thick disc and low-α and thin disc is fairly strong in the solar vicinity, this connection weakens significantly at larger galactocentric distances (e.g. Hayden et al. 2015).

On the other hand, the Galactic halo is crucial for understanding the formation and evolution of our Galaxy, as it holds the remnants of ancient cosmic events and provides a window into the processes that shaped the Milky Way. Our current understanding of Galactic stellar halo formation suggests a dual process: (i) gas accretion from cosmic filaments that drives secular evolution and in situ star formation, and (ii) the accretion of various mass building blocks, which contribute their baryonic material and dark matter to the larger host galaxy, adding to its mass as these building blocks are consumed (see recent reviews by Helmi 2020; Deason & Belokurov 2024). Thanks to the combination of Gaia astrometry and extensive spectroscopic surveys, our understanding of stars on halo-like orbits has advanced significantly in recent years. A notable finding was the discovery of stars with chemistry identical to the high-α disc but with halo-like orbits (e.g. Bonaca et al. 2017; Koppelman et al. 2018; Helmi et al. 2018; Belokurov et al. 2020). One interpretation of these stars is that they were born in situ and later dynamically heated to halo-like orbits, as some simulations predicted (e.g. Grand et al. 2020). The age distribution of these in situ halo or hot thick disc stars is cut off at lookback times of 8-11 Gyr (Gallart et al. 2019; Xiang & Rix 2022), suggesting that the heating event, possibly a merger, occurred at z ∼ 1−2 (see also Montalbán et al. 2021). Numerous other accreted systems have also been identified as part of the accreted stellar halo (Ibata et al. 1994; Helmi & White 1999; Belokurov et al. 2018; Helmi et al. 2018; Koppelman et al. 2019; Myeong et al. 2019; Yuan et al. 2020; Horta et al. 2021; Viswanathan et al. 2024a). The most significant of these is the Gaia-Enceladus-Sausage (GES) system, which contributes most of the accreted halo within ∼ 6−30 kpc of the Galactic centre (Naidu et al. 2020). There is also evidence for another major building block in the inner Galaxy (Heracles: Horta et al. 2021), whose stellar mass could possibly have been as massive as the GES (Horta & Schiavon 2025). It has been speculated that there may be a link between Heracles and a grouping of globular cluster populations classified in terms of their orbital and/or age-metallicity properties (Massari et al. 2019; Kruijssen et al. 2019; Forbes 2020). Based on the chemical-dynamical properties of this population and a comparison of these with cosmological simulations such as FIRE (Horta et al. 2024), it is likely that Heracles coalesced with the primordial Galaxy before the GES (z ≳ 2). However, deciphering if such populations arise from one single building block or multiple blocks is still an open question.

While the characterisation and origin of the main disc and halo populations are becoming clearer, the very earliest epochs of the Galaxy and the emergence of the disc remain poorly understood. The physical mechanism behind disc formation remains a subject of active debate. Theoretically, disc galaxies are thought to form through the dissipative collapse of gas within dark matter haloes, coupled with the conservation of angular momentum imparted during the early stages of large-scale structure formation (Eggen et al. 1962; Peebles 1969; White & Rees 1978; Fall & Efstathiou 1980; Ryden & Gunn 1987). Observationally, deep imaging and spectroscopy at high redshift (e.g. z>1−2) have revealed a population of rotationally supported, gas-rich discs already in place within the first few billion years of cosmic history (e.g. Förster Schreiber et al. 2009; Wisnioski et al. 2015; Übler et al. 2019), supporting the view that disc formation is an early and widespread outcome of galaxy evolution. In the context of the Milky Way, recent studies discover very and extremely metal-poor stars (VMP, [Fe/H]<−2.0, and EMP, [Fe/H]<−3.0) on disc-like orbits (Sestito et al. 2019, 2020; Di Matteo et al. 2020; Viswanathan et al. 2024b, 2025), driving the important question of when and how the Milky Way’s old disc formed. While there is a complex relation between the metallicity of a star and its birth time, typically metal-poor stars in the Milky Way are old, and therefore studying them can provide insights into the early epochs. The orbits and abundance distributions of old metal-poor stars today reflect the earliest phases of the Milky Way’s star formation and enrichment history (Beers & Christlieb 2005; Frebel & Norris 2015; Starkenburg et al. 2018; Lucey et al. 2019; Horta et al. 2021; Ardern-Arentsen et al. 2024; Horta et al. 2024; McCluskey et al. 2024). Belokurov & Kravtsov (2022) used [Al/Fe] from APOGEE (Majewski et al. 2017) to distinguish between in situ and accreted stars, identifying in situ stars down to [Fe/H] ∼−1.5. They found that the average tangential velocity of the in situ stars increased rapidly at [Fe/H] ∼−1. They interpreted this transition as the epoch when the Galaxy transitioned from a relatively disordered state to well-ordered rotation. Conroy et al. (2022) used abundances and ages from the H3 Survey, finding that this transition coincides with a non-monotonic rise in [α/Fe] abundances. This suggests a near-instantaneous change in star formation efficiency or gas inflow (see also Chen et al. 2024). At the metal-poor end, Rix et al. (2022) revealed a significant concentration of metal-poor stars ([M/H]<−1.5) near the Galactic centre, corroborating earlier studies of the very metal-poor component in the Galactic bulge region (e.g. Lucey et al. 2019; Arentsen et al. 2020). This also corroborates with other works such as Tumlinson (2010); Starkenburg et al. (2017); El-Badry et al. (2018); Horta et al. (2021); Ardern-Arentsen et al. (2024). The exact relationship between the proto-Galaxy and the metal-poor stars within the solar radius is still unclear, although Belokurov & Kravtsov (2023) provide evidence that the proto-Galaxy’s density follows a steep power law with galactocentric radius. Horta et al. (2024) used Milky Way analogues from high-resolution FIRE simulations, showing that the proto-Galactic populations in all their Milky Way analogues exhibit weak net rotation aligned with the present-day disc. Semenov et al. (2024) used the Illustris TNG50 simulations to interpret evidence of the Milky Way’s early disc formation. They suggest that rapid early mass growth was key to early disc formation. However, pin-pointing the exact driver is difficult, as both hot halo formation and gravitational potential steepening occur alongside disc emergence, indicating that a concentrated potential might be a result, not a cause, of disc formation (Dillamore et al. 2024).

Recently, Zhang et al. (2024) analyzed the Gaia DR3 RVS kinematic sample of the metal-poor population from Andrae et al. (2023a) using XGBoost algorithm on Gaia XP spectra in the Milky Way and identified two kinematic groups in velocity–[M/H] space: one a stationary accreted halo and the other with a net prograde rotation of 80 km s⁻¹ associating it with the proto-Galaxy. Chandra et al. (2024) used a similar Gaia RVS+XP sample with an additional α-separation to study the Milky Way evolution in a similar parameter space (orbital circularity–[Fe/H]) instead of velocities ([V_r, V_ϕ, V_z]−[Fe/H]). They found that the old high-α disc begins at [Fe/H] ∼−1 dex, with a dramatic transition from a proto-Galaxy to old disc.

For this work, we built on these studies in more detail and set out to study the transition of a proto-Galactic population to a high-α old disc using kinematics from Gaia DR3 RVS sample cross-matched with Gaia XP spectra [α/M] and [M/H] abundances from Andrae et al. (2023a) and Li et al. (2024). Section 2 presents the underlying data used in this work, and introduces the α-separation. In Section 3, we introduce the azimuthal velocity versus metallicity space where we study the evolution of the high-α disc. We also compare the results from this space with APOGEE chemical abundances to see if our inferences hold well with high-resolution spectroscopy. Section 4 presents a simple two-component Gaussian mixture model (GMM) corresponding to a disc-like and halo-like population using knots across the metallicity range to carry information between different metallicities, to interpret the transition between the proto-Galactic component to a rotation-dominated high-α disc. In Section 5, we repeat our analysis using orbital circularity instead of azimuthal velocity, as circularity is position independent. In Section 6, we discuss the conclusions derived in this work in the context of what we already know about the proto-Milky Way and the old high-α disc. We also compare our results with the recent literature and discuss some limitations and future prospects in this regard. Section 7 presents a summary and our conclusions.

2 Gaia DR3 XP+RVS Data

The Gaia DR3 catalogue provides low-resolution spectroscopy (XP) for approximately 200 million stars brighter than G=17.65, with radial velocities (RVS) available for a subset of around 30 million (Gaia Collaboration 2023; De Angeli et al. 2023; Katz et al. 2023). Subsequent studies have derived precise metallicity and α-abundance estimates from these spectra (Andrae et al. 2023a; Zhang et al. 2024; Martin, Starkenburg et al. 2024; Li et al. 2024). This dataset is exceptional due to its immense size, homogeneity, all-sky coverage, and relatively straightforward selection function. It is therefore an ideal resource to comprehensively investigate our Galaxy’s formation and evolution.

Andrae et al. (2023a) developed a powerful tool using a machine learning model called XGBoost – to estimate a star’s metallicity ([M/H]) from its low-resolution XP spectrum. They trained this model using a robust dataset from the Apache Point Observatory Galactic Evolution Experiment (APOGEE: Majewski et al. 2017), further enhanced with metal-poor stars from Li et al. (2022). The metallicities used for training are highly accurate, having been validated against established spectroscopic surveys. For more details on how they chose specific features and validated their catalogue, we refer directly to Andrae et al. (2023a). We leverage their metallicity data to create a refined sample of red giant branch (RGB) stars with line-of-sight velocities from Gaia DR3 radial velocity spectrometer (RVS). The following selection is made to the Andrae et al. (2023a) catalogue similar to their vetted RGB sample released as Table 2 in their work:

• We focused on bright stars (phot_g_mean_mag < 16) with highly reliable parallaxes (ϕ/σ_ϕ > 5). This ensures sufficient data quality for robust metallicity estimates.

• We excluded hot stars (teff_xgboost > 5200 K) as their measured metallicities can be misleadingly low.

• We applied a series of cuts based on colour and broadband magnitudes (logg_xgboost, M_{W 1}, G−W2, G_BP−W1) to select genuine red giants on the Hertzsprung-Russell (HR) diagram. These are listed here:

  • logg_xgboost < 3.5

  • M_{W 1} > −0.3−0.006 × (5500-teff_xgboost)

  • M_{W 1} > −0.01 × (5300-teff_xgboost), where M_{W 1} =W1+5 × log₁₀(π/100)

  • (G-W2) < 0.2+0.77 ×(G_BP-W1)

• We removed stars with a high probability (>90%) of belonging to globular clusters (GCs). This probability comes from a catalogue by Vasiliev & Baumgardt (2021).

This resulted in a sample of approximately 11 million RGB stars with Gaia XP metallicities and RVS radial velocities.

Li et al. (2024) created a catalogue of alpha-over-metallicity abundance ([α/M]) values derived from Gaia XP spectra using a neural network model. This model learns to use XP spectra to predict stellar labels but also uses and predicts the high-resolution APOGEE spectra that lead to the same stellar labels. This catalogue has also been cross-checked against existing surveys for accuracy, demonstrating a median error of only ∼ 0.05 dex in [α/M] for stars within our sample. We cross-matched our clean RGB sample with this catalogue, yielding a final sample of 9645972 RGB stars with metallicities and α-abundances from Gaia XP spectra, and full 6D phase space coordinates from Gaia astrometry and RVS line-of-sight velocities.

Figure 1 shows the logarithmic density distribution of [α/M] versus [M/H] of our final sample. While Figure 1 showcases a clear separation between high-α and low-α stars in terms of their chemical abundance, this distinction does not necessarily translate directly to their height above the Milky Way’s midplane. The chemical difference between these populations maps more strongly onto their radial distribution (distance from the Galactic Centre) than their vertical thickness (Hayden et al. 2015). Therefore, we leverage this chemical separation (high-α versus low-α) as a tool to differentiate the two chemically distinct disc populations and then analyze their spatial and kinematic properties (positions and velocities) independently.

Fig. 1 Logarithmic density of [α/M] vs [M/H]. The purple band represents the high- and low-α sequence separation defined in this work (see text for details). Stars in the purple band are excluded. The bulk of accreted last major merger (GES) is primarily restricted to the low-α population with our selection.

2.1 High- and low-α sequences

The purple band in Figure 1 shows the selection cut used to separate high-/low-α populations. Stars in the purple band are cut out to ensure that we have purer samples of high- and low-α stars. This is especially important for high-α stars since we use the high-α sample as the sample of stars tracing the evolution of the proto-Galaxy to high-α disc over metallicity and any contamination from the denser low-α disc can affect the purity of the high-α sample. Typically, the high-alpha sequence is attributed to stars formed in situ. However, in practise, the in situ versus accreted separation is not as simple as using an α-separation, the consequences of which are discussed in Section 6.4.

The selection is defined using the following equations: $$High-\alpha =\{\begin{array}{cc}[\mathrm{M}/\mathrm{H}]<-0.6\&[\alpha /\mathrm{M}]>0.28\hfill & \\ [\mathrm{M}/\mathrm{H}]\in [-0.6,0.125]\&[\alpha /\mathrm{M}]>-0.25\times [\mathrm{M}/\mathrm{H}]\hfill & +0.13\hfill \\ & \\ [\mathrm{M}/\mathrm{H}]>0.125\&[\alpha /\mathrm{M}]>0.1\hfill \end{array},$$(1) $$\text{Low}-\alpha =\{\begin{array}{cc}[\mathrm{M}/\mathrm{H}]<-0.8\&[\alpha /\mathrm{M}]<0.21\hfill & \\ [\mathrm{M}/\mathrm{H}]\in [-0.8,0.07]\&[\alpha /\mathrm{M}]<-0.21\times [\mathrm{M}/\mathrm{H}]\hfill & +0.045\hfill \\ & \\ [\mathrm{M}/\mathrm{H}]>0.07\&[\alpha /\mathrm{M}]<0.03\hfill \end{array}.$$(2)

This selection is somewhat different from the one implemented by Chandra et al. (2024) who use the same sample and abundances from Gaia XP spectra. The main difference between Chandra et al. (2024) and our selection are: (i) we use a more stringent selection for high-α to make sure that the bulk of accreted GES merger (Belokurov et al. 2018; Helmi et al. 2018) does not contaminate the high-α population (see Figure 1), and (ii) the purple band is also quantitatively larger to ensure a sample with the highest purity. We tested the validity of this separation by using the α-selection on APOGEE DR 17 data in the [M/H]−[α/M] space and using [Al/Fe]−[Mg/Mn] space to see if they occupy accreted, low-α disc or high-α disc stars regions as described by the tracks used in Horta et al. (2021). The contamination rate from this validation is as low as 5% and 8% for high- and low-α stars respectively. The high-α selection used here is stricter than the selection used in the literature, in order to remove as much accreted low-α stars as possible (mainly, GES, the last major merger). Even though we show the full sample in Figure 1, we restrict the analysis in the rest of this work to metallicities above −2.5.

2.2 Positions and kinematics

To estimate the distance to each RGB star, we use photo-geometric distance provided by Bailer-Jones et al. (2021), which is a Bayesian approach that incorporates both Gaia parallax measurements and a prior model of the Milky Way’s stellar density distribution. While using a simpler method based solely on zero-point-corrected parallax (Lindegren et al. 2021) yields similar results, the approach by Bailer-Jones et al. (2021) offers a more robust framework for estimating distances across stars with varying signal-to-noise ratios. To ensure that our astrometry was reliable, we focused on stars with a low renormalised unit weight error (ruwe < 1.4). This metric indicates good astrometric quality, potentially filtering out binary star systems. We adopted several standard assumptions: a Local Standard of Rest velocity (V_LSR) of 232 km s⁻¹, a distance of 8.2 kpc between the Sun and the Milky Way’s centre (GRAVITY Collaboration 2018), and the Sun’s peculiar motion of (U_⊙, V_⊙, W_⊙) =(11.1, 12.24, 7.25) km s⁻¹ (Schönrich et al. 2010). This allow us to calculate the positions and velocities for all the stars in our final sample.

Figure 2 shows the distribution of our sample in cylindrical galactocentric radius, R, versus height above the midplane, z, colour-coded by the mean metallicity value per pixel. The top panel of Figure 2 shows this distribution for all the stars in our sample, while bottom left and right panels show this distribution for high- and low-α samples, respectively. It is worth noting that the high-α stars probe a smaller range in R-z compared to low-α stars; most of the low-α stars at larger scale heights are expected to come from accretion events, and this is likely why the low-α sample covers a wider scale height. We see a strong and steep negative metallicity gradient with increasing z for the low-α sample; we also see this in the full sample, likely because it is dominated by stars in our low-α selection, while the high-α sample has a shallower negative metallicity gradient with increasing z. The high-α disc has a larger scale height compared to the low-α disc and we also see a clear disc flaring for the low-α disc (as seen for the Galactic thin disc) with a metallicity gradient towards larger R, reminiscent of radial migration (see e.g. Haywood et al. 2013; Hayden et al. 2015; Ratcliffe et al. 2024).

3 Milky Way populations in the [M/H]−v_ϕ plane

This section summarises the evolution of azimuthal velocity¹ versus metallicity for our high-/low-α samples. In examining this plane, we aim to gain an intuition of when and how quickly metal-poor populations (that for high-α likely contain components of the proto-Galaxy) go from a low net spin and kinematically hotter configuration to a more rotationally supported high-α disc population.

Fig. 2 Distribution of stars in cylindrical Galactic coordinates (R–z plane) colour-coded by their mean metallicities for all the stars in our sample (top), high-α selection (bottom left), and low-α selection (bottom right). We can see that dust near the Milky Way’s midplane significantly impacts our survey selection. We see a sharp negative metallicity gradient with respect to height above the disc plane in all stars and low-α stars. High- α stars have a shallower negative metallicity gradient.

3.1 Azimuthal velocity versus metallicity trends

Figure 3 illustrates the formation and evolution of the Milky Way using azimuthal velocity, v_ϕ, as a function of metallicity, [M/H], for our full sample (top), and our high-/low-α samples (bottom). For all panels, each [M/H] bin has a size 0.01 dex, and we show the sum normalised histogram of azimuthal velocity, plotted as a 2D column-normalised histogram² for all, high-, and low-α stars.

Although it is intuitive to read metallicity as an age sequence, we know that the age-metallicity relation (AMR) of the Milky Way is not monotonic (Bensby et al. 2014; Xiang & Rix 2022; Gallart et al. 2024; Xiang et al. 2024). This is true for all stars (top panel, which has a mix of in situ and accreted stars), and low-α stars (bottom right panel), as they are made of two different stellar populations (thin disc and accreted mergers). However, Xiang & Rix (2022) have shown that the stars in high-α sequence have a relatively consistent and monotonic AMR. They also show that the high-α disc reached solar metallicities at around ∼ 8 Gyr in lookback time. Therefore, the high-α panel should trace the first ∼5−6 Gyr of the Milky Way’s evolution, which is also proposed by Chandra et al. (2024) in their orbital circularity versus metallicity plane. All the 2D histogram in Figure 3 have 16th, 50th, and 84th percentile tracks of v_ϕ versus [M/H] overlaid. The median and percentile tracks look very similar if we restrict our analysis to solar neighbourhood stars (d<3 kpc), as the velocities are less position dependent in a small volume around the Sun. However, in order to preserve the number statistics, especially in the low-metallicity end, we use the entire sample in the rest of this paper instead of restricting only to solar neighbourhood stars. As an additional check, we use orbital circularity versus metallicity instead of azimuthal velocity, which is expected to be position independent in Section 5. We discuss the different panels separately in the following subsections.

Fig. 3 Column-normalised (by sum) 2D histogram of stars in the [M/H]–vϕ plane (azimuthal velocity vs metallicity) for all the stars (top), high-α selection (bottom left), and low-α selection (bottom right). The running median track is shown as dashed black line and the 16th and 84th percentile tracks are shown as black lines in all panels. The running median tracks for the low-α panel look more like a step-function, while high-α tracks is shallower, supporting the gradual monotonic increase, that can be interpreted as a gradual spin-up from old proto-Galactic populations to the present day high-α disc with respect to metallicities.

3.1.1 All stars

In Figure 3, top panel, we see strongly rotational (v_ϕ ∼ 200 km s⁻¹), kinematically cold thin (young, low-α) disc dominating higher metallicities down to [M/H] ∼−0.7. Below these metallicities, we see the kinematically hotter thick (pld, high-α) disc population (v_ϕ ∼ 180 km s⁻¹) which then has a disconnection at [M/H] ∼−1.0. Furthermore, it is possible to see the GES merger (and other accretion events) dominating radial orbits (v_ϕ ∼ 0 km s⁻¹) at metallicities below [M/H]<−1.0. This creates a step function behavior (also seen in the running median tracks) around the metallicities between −1.3 and −0.9 which has been reported in the literature as the spin-up phase, where we have the dramatic transition from what has been dubbed a “proto-Galactic” population with low net rotation to a rotation-supported disc component (Belokurov & Kravtsov 2022; Chandra et al. 2024; Zhang et al. 2024; Kurbatov et al. 2024).

3.1.2 Low-α stars

In Figure 3, bottom right panel, we see strongly rotation-supported low-α disc (v_ϕ ∼ 200 km s⁻¹) dominating higher metallicities ([M/H] > −1) and is almost fully disjoint to the accreted population at lower metallicities ([M/H]<−1) that are radial and present (almost) isotropic rotational velocities (v_ϕ ∼ 0 km s⁻¹). Here, we can clearly see a step function behaviour at [M/H] ∼−1.0, which demarks the transition from the accreted populations to the low-α disc. We see little-to-no evolution in the azimuthal velocity and velocity dispersions for the low-α disc stars with increasing metallicities. This is also because metallicity clearly does not trace age for low-α disc stars, but instead the birth radius of the star (e.g. Ratcliffe et al. 2024), i.e. metallicities for low-α stars cannot be read as a temporal axis.

Fig. 4 Column-normalised (by sum) 2D histogram of stars in the [M/H]−vϕ plane for APOGEE DR 17 giants (with the same high- and low-α separation as shown in Figure 1) for all the stars (right), high-α selection (left), and low-α selection (middle). The running median track is shown as dashed black line and the 16th and 84th percentile tracks are shown as black lines in all panels. The median, 16th and 84th percentile tracks for our Gaia XP sample are shown in grey for comparison. The tracks between Gaia XP and APOGEE samples show remarkable resemblance in high-α and small differences for low-α and all stars due to lower contamination in low-α selection in APOGEE (see Appendix A).

3.1.3 High- α stars

In Figure 3, bottom left panel, we see high-α stars with their azimuthal velocity evolving across the metallicity range monotonically. Both the 2D column-normalised histogram and the median tracks show the evolution of a less rotationally supported (kinematically hotter) metal-poor component (likely the proto-Galaxy), with a small but non-negligible mean azimuthal velocity value (v_ϕ ∼ 50 km s⁻¹, McCluskey et al. 2024; Semenov et al. 2024; Horta et al. 2024). As metallicity increases, we see the high-α sample increasing in v_ϕ towards more prograde orbits, gradually transitioning from a hotter less rotating component to a rotation-dominated high-α disc (v_ϕ ∼ 180 km s⁻¹). The velocity dispersion for high-α disc orbits becomes smaller with increasing metallicities. Therefore, Figure 3 likely illustrates the emergence of the high-α disc, revealing the evolution of proto-Galactic populations gradually spinning-up with increasing metallicities.

It is important to note that our [α/M] selection would not fully remove accreted stars in the high-α sample, especially in the metal-poor end of accretion events, [M/H]<−1.3, where GES and other debris like Heracles (Horta et al. 2021) appear. Therefore, accreted halo stars would still be present (although probably contributing a small fraction) in the metal-poor end. Additionally, the mean velocities and velocity dispersion that we show are present day velocity distributions. Therefore, Figure 3 shows the instantaneous velocity information for these populations, and not the velocity profile at formation; for these metal-poor (old) populations, these two could be drastically different due to dynamical heating over time, for example. It is also important to note that the hotter high-α disc (in situ halo or hot thick disc) population, is still present in the high-α sample, but smaller in number than the rotation-supported high-α disc stars. These “hot” high-α disc stars can be easily seen in row-normalised [M/H]−v_ϕ plane as shown in Appendix C.

From the bottom left panel of Figure 3, we find evidence in support of a metal-poor high-α population (likely part of the proto-Galaxy) gaining rotation at a slower pace across wider range of metallicities (between −2 and −0.7) that eventually settles into a high-α disc population.

3.2 Azimuthal velocity versus metallicity tracks using high-resolution APOGEE abundances

Figure 4 shows 2D column-normalised histograms of [M/H]−v_ϕ plane for high-α (left), low-α (middle), and all stars (right) for APOGEE DR17 stars with the same α-selection described by Equations (1) and (2). In APOGEE DR17, we select stars with log g>3.5, excluding potentially problematic data based on quality flags such as ASPCAPFLAG, and WARNING, or BAD flags defined in ASPCAP for T_eff, log g,[M/H], and [α/M]. The running 16th, 50th, and 84th percentile tracks for the APOGEE data are shown as black lines and the tracks from our final sample in Figure 3 are shown in grey in all three panels. We find very good resemblance between the track in the APOGEE and Gaia XP data, especially for high-α stars stars. It showcases that the scenario of a gradual spin-up of the proto-Galaxy to the high-α disc population (across metallicities) is not due to large errors in Gaia XP abundances, as it is also supported by a much more reliable high-resolution chemical abundance data from APOGEE. The low-α density distribution (and median tracks) are the ones that vary the most between the two samples. This would be expected if contamination from in situ centrally concentrated stars from the proto-Galaxy were present in our low-α sample, while APOGEE has a much cleaner and purer low-α selection, due to higher quality of the chemical abundances. This can also be seen by the low-α disc and halo populations being completely disjoint in the low-α 2D histograms. This is discussed more in detail in Appendix A using [M/H]−v_ϕ plane towards and away from the inner Galaxy and noting the differences for low-α 2D histograms in Figure 4. Lastly, the all stars panel in APOGEE also resembles what we see in our Gaia XP sample, with stars at lower metallicities dominated by accreted stars and proto-Galaxy populations with low net spin, and higher metallicities dominated by the kinematically hotter high-α disc and kinematically colder low-α disc at even higher metallicities.

Given these results, we are confident that our selection presented in Figure 1 is efficient to separate different stellar populations in the Milky Way, where we can analyse the [M/H]−v_ϕ plane. In the following section, we set out to model these data using a new mixture model, with the aim of quantifying the point at which metal-poor populations transition into the more rotationally supported disc populations; we also aim to assess the fraction of halo/disc populations at different metallicity bins.

4 Modelling the high-/low-α stars in the [M/H]−v_ϕ plane

To decipher the underlying structure of the data in [M/H]−v_ϕ space for both the high-/low-α samples, we create a simple model described by two Gaussian components in v_ϕ (one that is more disc-like and other that is more halo-like) the probabilities of which is conditioned on the value of the metallicity, thereby carrying over the two-component GMM information across metallicities. This model is implemented using numpyro (Phan et al. 2019; Bingham et al. 2019), a lightweight probabilistic programming library that provides a numpy backend for pyro.

We argue that this model is more well-suited for this problem than a standard GMM decomposition (even if the number of components is set free and determined using the Bayesian Information Criteria) in bins of [M/H], because the underlying Gaussians are independent of each other in between different [M/H] bins. Therefore, instead of manually linking different Gaussians across metallicity, our model P(v_ϕ |[M/H]) is able to connect the two Gaussians across the metallicities as a continuous conditional distribution. This is also advantageous as it allows us to follow and interpret the different Gaussian components and thus examine how they vary across metallicity space. An underlying assumption of our method is that we can model the low-α and high-α discs as a Gaussian distribution in v_ϕ space. To first order, this is a good approximation (see middle panels of Figure 9); it is also a good approximation for the stellar halo, as our sample is dominated by the debris of one major accreted satellite (i.e. GES).

Moreover, our model allows us to interpret how different stellar populations evolve in v_ϕ across [M/H]. It also enables us to measure the transition from a halo-like non-rotating population to a disc-like rotating one for both the high-/low-α samples. We reason that the two-component fit is also physically motivated as we know that the Galaxy has a halo and a disc in the high- and low-α regime. While simple in model form, we show here how this model is able to accurately capture the data.

In practice, the model takes in [M/H] and v_ϕ as input parameters, and uses priors on two (separate) Gaussian components to split the data into a halo-like and disc-like population across [M/H]. In detail, we choose 16 evenly spaced knots on a log scale³ across the metallicity space from −2.5<[M/H]<0.01. For each batch of data in these bins, our model computes the fraction of the data that can be modelled by one Gaussian compared to the other (i.e. the fractional contribution of each Gaussian); this allows us to quantify how much of the data is better fitted by a halo-like populations versus a disc-like one in each [M/H] bin. We set the priors of each Gaussian component as a normal distribution with 𝒩(μ, σ), where μ=0 km s⁻¹ and σ=150 km s⁻¹. For the priors on the weights (or fractional contribution) of each Gaussian at the locations of each [M/H]-knot in the spline, we use a Dirichlet prior with a constant concentration value of 0.5 for both the components⁴. The means, standard deviations, and relative weights for the two Gaussian components are then sampled using a Hamiltonian Monte Carlo (HMC) inference, using the No U-Turn Sampler (NUTS). We run the sampler for 200 warm-up steps, and 100 sampling steps. This results in the Markov chain Monte Carlo fits for high-α and low-α stars in our sample. We show the graphical model representation of this model in the left panel of Figure 5 and a list of model parameters and their functional forms are provided in the second column of Table 1.

When running this method on our high-/low-α samples, we find that the r_hat split Gelman Rubin diagnostic parameter is less than 1.1 for the majority of the [M/H] bins, indicating that the chains have converged. Moreover, the chains look stable, and n_eff, the number of effective samples, is at least 30 or more for the means, standard deviations, and relative weights. We do note however that for the lowest [M/H] bins, the r_hat values are higher for high-α stars (up to 1.72 for the knots in the metal-poor regime, even though the chains are converged) than for the low-α stars (oscillate below 1.1 regardless of the warm up steps and sample size for the MCMC chains). We tested varying the number of warm-up and sampling steps, but found that this did not affect our results. This result indicates that the samples we generate are a fair representation of the posterior probability distribution over the model parameters for the low-α stars. This is also the case for the metal-rich bins in the high-α sample. However, at lower [M/H] for high-α stars, our model requires more flexibility to generate a fair representation of posterior distributions. The low-α sample is comprised of two clearly distinct populations, the low-α disc and GES debris at different metallicities. Conversely, the high-α sample is comprised of stars that appear to follow a single trajectory across the entire metallicity range. Moreover, metal-poor high-α populations are likely an amalgamation of stellar populations (i.e. Heracles, old in situ material). Thus, our results hint that our fixed means two-Gaussian components are too simple to be able to capture the distribution of this data perfectly. We can see the results from this model in Figure 3. Here, the low-α stars cleanly separate into two distinct stellar populations, each of which dominates a different metallicity range; we note that the median tracks display almost a step function. At high (low) [M/H], the low-α disc (GES debris) dominates. This is not the case for the high-α star sample, that displays a monotonically increasing behaviour in v_ϕ for increasing metallicity. Therefore, a model with the means and standard deviations of the azimuthal velocity varying across the metallicities may be more suited for the high-α stars.

However, despite this limitation, it is useful to compare the results from our model for the high- and low-α population. Thus, we proceed with this simple two-component model for both high-/low-α populations to assess the fractional contribution of disc-like and halo-like components in each sample.

4.1 Interpretation of the frozen means model in the [M/H]–v on high/low-α stars

Figure 6 shows the fractional contribution of a halo-like component in orange and disc-like component in purple with the 1 σ uncertainty in the fractions as orange and purple shaded bands, respectively. The knots in metallicity chosen to run the model are shown as scatter points. The converged velocity means for high-/low-α stars are 188 km s⁻¹ and 223 km s⁻¹, respectively. These values match well the mean rotational velocity of the high-α (thick) and low-α (thin) discs. The corresponding velocity dispersions are 39 km s⁻¹ and 25 km s⁻¹ for high-/low-α discs, showing that the high-α disc is kinematically hotter than the low-α disc. Furthermore, the halo-like components have a more radial mean velocity and hotter velocity dispersion when compared to their respective high-/low-α disc samples. Here, the high-α halo-like component has a smaller velocity dispersion (74 km s⁻¹) than low-α stars (107 km s⁻¹). This could simply be due to the fact that the two-component fit with fixed mean velocities and standard deviations might not be the best representation of the underlying data for high-α stars. Conversely, it could also imply that the high-α halo component is kinematically colder than the low-α (GES) one.

The main difference we can see between the fractional contribution of disc-like and halo-like stars in high-/low-α sub-samples is that the low-α stars completely switch between halo and disc around a very narrow −1.1 and −0.8. This steprange in metallicities between function behaviour is clearly seen in the right panel of Figure 6. The small fall in halo-like fractions at lower metallicities could simply be due to noisy data in the metal-poor end. For the low-α stars, within the uncertainties in the weights, we can clearly see that the metal-poor end is fully composed of halo-like stars ([M/H ≤−1.1]) and the metal-rich end is fully composed of disc-like stars ([M/H] ≥−0.6). Conversely for high-α stars, we do not see such a steep turn-over between the halo-like and disc-like samples. In contrast, the disc-like component gradually increases with increasing metallicity in its relative fraction, while the halo-like component gradually decreases. The disc-like component is present with 18% relative contribution in the very metal-poor end ([M/H]<−2) in this simple model. This result favours a more gradual spin-up (with respect to metallicities) scenario. However, as discussed above, our model is not able to fully capture the distribution of the high-α stars in the metal-poor regime. To improve the underlying model to better capture the data in order to quantify the spin up of the high-α disc, we run a separate model in the following section that is able to let the mean velocity and velocity dispersion evolve over increasing metallicities.

Fig. 5 Graphical model representation of the evolving means model, modelling the conditional distribution P(vϕ |[M/H]) with a two-component GMM, where the means and standard deviations are fixed (left) and varying (right) with respect to the underlying metallicities. The evolving or varying means model is also implemented for circularity as the conditional distribution or P(η |[M/H]). The subscript 1 and 2 stands for the disc-like and halo-like components, respectively, with μ and σ as the Gaussian mean and standard deviations in vϕ or η, conditioned on the value of [M/H]. w stands for the relative contribution of the disc-like component with the relative contribution of the halo-like component defined as (1−w). A list of the model parameters and their priors and functional forms are given in Table 1.

4.2 Quantifying the evolution of the high-α disc with metallicity

In order to quantify the spin-up of the high-α disc population, we modify the simple two-component model from the previous section to allow the mean velocity and velocity dispersion to vary between the spline knots across the metallicity range. We fix the mean velocity of the halo component to be 0 km s⁻¹ across the metallicity knots (the value converged from the HMC run using fixed means and standard deviations). This is fixed in order to separate the metal-poor high-α in situ population (likely components of the proto-Galaxy, with a small but non-negligible prograde rotation) from the more radial halo population, as they heavily overlap. However, we let the standard deviation of the halo-like component vary across the metallicity knots with a truncated normal prior with a mean of 100 km s⁻¹ and standard deviation of 50 km s⁻¹, with lower and upper limits restricted between 50 and 200 km s⁻¹.

It is important to note that the high-α selection is not perfect, and still captures the metal-poor end of accreted mergers (like GES). Therefore, fitting the halo-like component is still important in our high-α sample. The mean and standard deviation of the velocity of the disc-like component are set to vary with increasing metallicity. This will allow us to measure the spin-up of a proto-Galactic population to a rotation-dominated high-α disc. The metallicity range used for high-α stars with this model is between −2.0 and +0.1. For the more metal-rich bins, we undersample the data to have an equal number of stars in each bin. The number of stars is set to the number in the lowest metallicity bin. This downsampling reduces the sample size for the most metal-rich bins, that in turn reduces the computation time. This also allows us to have equidistant knots, allowing the sampling of the metal-poor end as well as the metal-rich end (as opposed to the log-space knots used in the previous subsection, which sample the metal-poor end less than the metal-rich end). We choose 8 metallicity knots (equidistant in linear space) for the relative fraction and 9 metallicity knots (equidistant in linear space) for the velocity means and standard deviations. The mean velocity, velocity dispersion, and relative weights are computed for each knot in metallicity and spline interpolated for the data in-between the knots.

In order to measure the spin-up of the disc-like component, and to enable the information of the disc-like component to be carried across metallicity bins, we implement a Gaussian process such that every finite collection of the azimuthal velocities (measuring the mean azimuthal velocity of the disc-like component) indexed by its metallicity has a multivariate normal distribution described by a rational quadratic kernel function⁵ $${cov}_{\text{spin-up}}({[\mathrm{M}/\mathrm{H}]}_{\mathrm{i}},{[\mathrm{M}/\mathrm{H}]}_{\mathrm{j}})={\sigma }^2{(1+\frac{{({[\mathrm{M}/\mathrm{H}]}_{\mathrm{i}}-{[\mathrm{M}/\mathrm{H}]}_{\mathrm{j}})}^2}{2\alpha {\ell }^2})}^{-\alpha },$$(3) where σ² is the overall variance, ℓ is the characteristic length scale of the covariance function, that describes the range in [M/H] to which the covariance function typically holds, and α is the positive scale-mixture parameter of the covariance function (α > 0), that in simple terms, describes the curvature of the covariance function. This function models the covariance between each pair of ([M/H]_i,[M/H]_j). The standard deviation on the mean of the azimuthal velocity describing the disc-like component has a uniform prior between 50 and 200, the characteristic length scale has a uniform prior between 0 and 1, and the scale-mixture parameter has a uniform prior between 0 and 4. All these parameters are sampled with the MCMC along with the means, standard deviations, and weights of the disc-like and halo-like components across the metallicity bins. The mean azimuthal velocity of the disc-like component has a multivariate normal distribution prior centred at 120 km s⁻¹ described by the rational quadratic kernel function covariance matrix shown in equation (3). The mean azimuthal velocity is therefore defined as the mean function together with the kernel function that define the Gaussian process distribution of the azimuthal velocity of the disc-like component with varying metallicity. The standard deviation is described by a half normal prior centred at 150 km s⁻¹. The relative weight of the disc-like component is forced to be monotonically increasing using the positive_ordered_vector constraint on the argument, which automatically forces the halo-like component to be monotonically decreasing, removing noisy fits in the low-metallicity end.

The NUTS sampler with a HMC inference is run on the sampled data with the model described above for 100 warm-up steps, and 50 samples to generate from the Markov chain for the high-α stars. We show the graphical model representation of this model in the right panel of Figure 5 and a list of model parameters and their functional forms are provided in the third column of Table 1. The high-α sample chains look much more stable and converged with this evolving velocity means and standard deviations model, they have r_hat below 1.1 (between 0.98 and 1.03), and n_eff, the number of effective sample is at least 30 or more for the means, standard deviations and relative weights, and about 20 for the covariance function parameters. The converged characteristic length is 0.75 (showing larger correlation between the velocities at different metallicities), and the scale-mixture parameter is 1.95. In summary, these results suggest that this model provided posterior distributions that better describe the underlying high-α sample.

Figure 7 shows the converged parameters as a function of metallicity. As in Figure 6, the metallicity knots are shown as scatter points and 1 σ uncertainties on the converged parameters are shown as coloured bands. The disc-like component is shown in purple and the halo-like component is shown in orange. The left panel of Figure 7 shows the evolution of the mean azimuthal velocity and azimuthal velocity dispersion for the high-α disc-like component. Here, the trend gradually increases from v_ϕ ∼ 50 km s⁻¹ at [M/H] ∼−2 to v_ϕ ∼ 200 km s⁻¹ at [M/H] ∼ 0. This result can be interpreted as the high-α disc spinning up from a proto-Galactic population to a rotation dominated disc. While previous work have alluded to this in the literature (Chandra et al. 2024; Zhang et al. 2024), this is the first time that the evolution of the high-α proto-Galaxy-to-disc population has been quantitatively measured. In more detail, our results suggest that the proto-Galactic phase (pre-disc) lasts for approximately 0.3 dex (between −2 and −1.7 dex). After this, the proto-Galaxy populations gain azimuthal velocities as metallicity increases – in an almost linearly fashion – until they settle into a disc at around −0.5 dex. In summary, our results show that the so-called “spin-up” phase of the Galaxy happens gradually across a large range of [M/H], starting from metallicities as low as [M/H] ∼−1.7.

Furthermore, throughout this phase, the velocity dispersion of the disc-like component decreases with metallicity, which also highlights how as the population gains v_ϕ, it becomes less dispersion dominated. The velocity dispersion of the halo component (orange, that has fixed mean of v_ϕ=0 km s⁻¹) is also decreasing with increasing metallicities. This could be due to different substructures dominating different metallicities. For example, we know that the GES merger has a lower velocity dispersion (of about 50 km s⁻¹) compared to the rest of the halo resting at about 100 km s⁻¹ velocity dispersion (see Figure 3). The rise in standard deviation at higher metallicities (when the high-α disc is in place) is not physical, but is rather caused by the model trying to make a broad halo component to fit the asymmetric tail of the high-α disc. This is one of the disadvantages of using a Gaussian distribution.

The middle panel of Figure 7 shows the ratio of v_ϕ to σ_ϕ, measuring how rotationally supported or “discy” the stars are. This ratio is basically zero (extremely pressure-supported) for the halo-like component, as we set the mean of the halo-like component to be close to 0 km s⁻¹. However, the disc-like component has a clear rise in v_ϕ/σ_ϕ, reaching up to a ratio of 6 at solar metallicities. Moreover, the right panel of Figure 7 shows the relative contribution of halo-like (orange) and disclike (purple) components at different metallicities. The accreted halo contribution decrease quickly with increasing metallicities. It dominates the distribution only at the lowest metallicity bins, below [M/H] ≤−1.7. In the very metal-poor end ([M/H]<−2) of our high-α selection, the halo-like to disc-like component (i.e. halo-like to proto-Galaxy-like) ratio is 40%: 60%. It is worth noting that we do not call this component “disc”, but simply refer to this component as “disc-like” for the modelling purpose. This component captures the proto-Galaxy to spin-up phase to high-α disc. The disc-like component’s fractional contribution increases slowly from [M/H] ∼−2, with an approximately constant gradient, up to [M/H] ∼−0.5. Upon reaching this point, the disc-like component (i.e. high-α disc) dominates the sample. All these results are highly in favour of a gradual spin-up for the disc-like component with increasing metallicities. In terms of the evolution of the Galaxy, our results highly favour the scenario where a proto-Galaxy population with low (but non-negligible) v_ϕ profile spins up into a fully rotation-dominated high-α disc.

In Figure 8, we show the [M/H]−v_ϕ plane for high-α sample, with the sampled data in 2D column-normalised histogram on the left, the fitted model normalised with the integral under the curve equal to 1 (similar to column normalisation by sum) in the middle, and normalised residual (i.e. data – model) on the right. All the panels have the running 16th, 50th, and 84th percentile tracks from Figure 3 shown as black lines. We can clearly see that the sampled data follows the median tracks very well. We construct the model using the splines on the mean velocities, velocity dispersions, and fractional contributions for the two-component fits. The normalised residual is constructed by subtracting the probability density function of the sampled data with the model’s probability density function (PDF) in each cell (with bin sizes of 6.25 km s⁻¹ in v_ϕ and 0.04 in [M/H]). The residuals are scaled by the sum of sampled data and the model’s PDF. This way, the residual only goes from −1 to 1 in value. If the residual is less than 0 the model predicts more stars than the data shows, and if the residual is greater than 0 the data has more stars in that region than the model predicts.

When inspecting Figure 8, the first thing one catches by eye in the residuals (right panel) is the horizontal patch of red stars at higher metallicities (between −0.5<[M/H]<0). Here, the model predicts fewer stars than what is present in reality. We conjecture this is because the v_ϕ distribution of high-α disc stars is asymmetric probably caused by the mechanism of asymmetric drift and is strongly non-Gaussian with a long tail towards lower v_ϕ for this metallicity bin. Thus, the model is not able to capture well these stars. This asymmetric drift is stronger in high-α disc than the low-α disc. However, it is present in both populations (Anguiano et al. 2020). The same effect can also be seen in the model with frozen means and standard deviation for azimuthal velocity in the low-α sample. Due to this asymmetric tail towards lower velocities, we find that the model is underfitting the data by 22−26% at higher metallicities ([M/H] > −0.5), while the under/overfitting is as low as 2−5% at lower metallicities.

The evolving mean high-α model presented here models well the bulk of the stars (also represented by the percentile tracks). However, we note that it struggles to fit the stars in the periphery of the v_ϕ distribution (edge of the grid values in [M/H]−v_ϕ plane), mostly due to Poisson error. When compared to the residuals for the model with frozen means and standard deviation for azimuthal velocity in the high-α sample, the evolving means model has much smaller systematic effects. This is due to the frozen means model’s underlying assumption of frozen mean velocity and velocity dispersion. Therefore, this improved model with evolving means is a much better representation of the high-α sample. Furthermore, it is important to note that in the residuals, the model is fully represented by a gradual (almost linear) rise in azimuthal velocity over the entire range of metallicities. If the data/model were better represented by a rapid and more exponential growth of v_ϕ with respect to a narrow range of metallicity – as reported in the literature –, we would find large systematic effects in our residuals (that are not seen). The absence of such systematic effects is more supporting evidence that the metal-poor high-α (or proto-Galaxy population) gradually spins-up to a rotation-supported high-α disc over a wide range of metallicities.

To our knowledge, this model is a first attempt at a simplified, physically motivated, and easily interpretable representation of the azimuthal velocity versus metallicity plane in the high-α regime, capturing the evolution of the first 5-6 Gyr of the proto-Milky Way populations to the high-α disc. However, given the consequences of the simplified Gaussian distribution assumption, this model is a more qualititative representation in the metal-rich end ([M/H] > −0.5).

Fig. 6 Fractional contribution of disc-like (purple) and halo-like (orange) GMM components as a function of metallicity for high- and low-α stars. The scatter points show the location of the knots chosen to run our model. The bands are 1 σ uncertainties on the converged weights. These fractions are computed for a fixed Gaussian for halo and disc-like stars when the MCMC is converged. The converged azimuthal velocity and velocity dispersion is shown in the legend, shown in units of km s−1. The low-α model clearly shows two separate populations (thin disc and accreted) shown by the step function that could be misclassified as a rapid spin-up, while the high-α model shows a shallower monotonically increasing(decreasing) profile for the disc(halo) with [M/H].

Fig. 7 High- α model with evolving means and standard deviation for azimuthal velocity of the disc-like component. Left: Evolution of azimuthal velocity and azimuthal velocity dispersion as a function of metallicity, capturing a monotonic transition of a chaotic dispersion-dominated state (predominantly made of proto-Galactic populations) to a rotationally supported state with disc-like motion in purple; the velocity dispersion evolution of halo is shown in orange. Centre: Evolution of azimuthal velocity over azimuthal velocity dispersion, Vϕ/σϕ as a function of metallicity for the disc phase and halo component is shown in purple and orange, respectively. Right: Fractional contribution of disc phase (proto-Galaxy to high-α disc) (purple) and halo-like (orange) GMM components as a function of metallicity. In all the panels, the scatter points show the location of the knots chosen to run our model. The bands are 1 σ errors on the converged weights, velocities, and velocity dispersions. The fractions are computed for an evolving Gaussian for the disc phase and fixed Gaussian for the halo phase when the MCMC is converged.

Fig. 8 Distribution of sampled data column normalised by sum (left), model colour-coded by the probability density function (centre), and the normalised residual (right) in the [M/H]−vϕ plane (azimuthal velocity vs metallicity) for high-α evolving means model. The flat distribution of red stars at higher metallicities in the residuals show that the model is not able to fit the non-Gaussianity (long tail towards lower velocities) of the high-α disc. The systematic scatter in the periphery is due to random error in the data. However, for [M/H]<−0.7, the lack of strong systematic patterns, the low-amplitude, and the small scatter of the residuals in the central regions (with the major bulk of the data) verifies the validity of the model. The median, 16th and 84th percentile tracks for our data are shown in black lines to guide the eye towards the bulk of the data.

5 Evolution of orbital circularity with [M/H]

The orbital circularity metric, η, defined with respect to the present-day disc alignment, ranging from −1 (perfectly retrograde, in-plane) to 1 (perfectly prograde, in-plane), interprets a star’s orbital properties within the Galaxy. Intermediate values represent elliptical orbits (closer to 0, isotropic, radial). This quantitative measure of orbital shape is very useful in understanding the dynamical processes that shaped the formation of the high-α disc. We compute orbital circularities for our sample of stars using the gala software package (Price-Whelan 2017). We integrate the orbits of stars within a realistic Milky Way model for the Galactic potential, MilkyWayPotential2022 (Price-Whelan et al. 2022), which incorporates four crucial components: a dense central core, a surrounding bulge of stars, a flattened disc of stars and gas, and a vast, spherical dark matter halo. This model is also calibrated to match observations of the Milky Way’s rotation curve by Eilers et al. (2019) and a compilation of Milky Way’s total mass enclosed (see Hunt et al. 2022). Using the adopted gravitational potential, we calculate the vertical angular momentum (L_{z, circ}) and total energy (E) for a perfectly circular orbit. By interpolating the resulting L_{z, circ}(E) curve for each observed star’s total energy, we determine the orbital circularity as η=L_z/L_{z, circ}.

The evolution of orbital circularity with respect to increasing metallicities is shown as 2D column normalised (by sum) histograms in the top panels of Figure 9 for high-α, low-α, and all stars. All the 2D histograms have 16th, 50th, and 84th percentile tracks of η versus [M/H] overlaid. To compare the orbital circularity with the azimuthal velocities used in this work, the 1D histograms of azimuthal velocities in bins of metallicities are shown in the bottom panels of Figure 9 for the high-α, low-α, and all samples. We also show the 1D histograms of orbital circularity, η, in bins of metallicities in the second row of Figure 9. The use of orbital circularity is very important because, unlike v_ϕ, orbital circularity is position independent (especially at larger distances away from the solar neighbourhood). From the high-α panel of Figure 9, we can confirm the gradual evolution of circularity from a non-zero median circularity (η ∼ 0.1) metal-poor (halo-like) population to a rotation-dominated (η ∼ 0.9) disc-like component over a broad range of metallicities, ranging between −2.5<[M/H]<−0.7 (see the median tracks overlaid), if they are composed of a single stellar population. This is in favour of the gradual spin-up phase of the Milky Way’s high-α disc over increasing metallicities. This can also be seen in the 1D histograms of v_ϕ (fourth row) and η (second row), with the lighter colours (lower metallicities) having a bimodal distribution; this bimodal feature is likely attributed to the superposition of a halo-like population with no rotation and a proto-Galaxy-like population with small rotation (Horta et al. 2024). The bimodal distribution at lower metallicities is strikingly clear for orbital circularities as annotated in the second row of Figure 9. Furthermore, the low-α panel of Figure 9 reveals that the low-α stars are made of two (almost) disjoint populations: an accreted halo and low-α (thin) disc. The halo dominates the lower metallicity bins, whilst the low-α disc dominates the higher metallicities, as expected. This can be seen also in the median tracks (top panel of Figure 9), which delineate a trajectory similar to a step function. The low-α 1D histograms of v_ϕ also show that the accreted halo (lighter colours, lower metallicities) dominates at v_ϕ ∼ 0 km s⁻¹ from [M/H] ∼−2 to [M/H] ∼−1, without a decrease in the relative weight (i.e. the peak in the distribution stays relatively the same. After [M/H] ∼−1, the v_ϕ distribution rapidly changes into a highly rotating, disc dominated, population, that spans over 0.7 dex in metallicity (darker colours, higher metallicities).

In the all stars bottom panel of Figure 9, we see a combination of the high-α and low-α components: accreted halo+proto-Galactic populations dominating the lower metallicities, high-α disc populations emerging from −1<[M/H]<−0.7, and low-α disc populations taking over the distribution at metallicities higher than [M/H] ∼−0.6. The same picture can be deduced from the 1D histograms of v_ϕ, that trace the same distribution as the circularity. Moreover, the median tracks on the 2D histograms show a rapid rise in circularity at around [M/H] ∼−1.0, which is driven by a transition from the accreted (mostly GES) debris to the high-/low-α discs, similar to what is seen in the low-α sample. However, in the all stars panel, the jump is not as sharp due to the presence of high-α disc populations as well. The conclusions are therefore consistent with the use of orbital circularity or azimuthal velocity. Given the position dependence on the azimuthal velocity, the use of orbital circularity brings more confidence that the gradual rise in rotational support over increasing metallicities for the high-α population is bonafide.

This result is important. If one does not account for α-separation like done in this work, a conclusion of the disc spin up in the [M/H]−v_ϕ diagram could be interpreted as rapid, which we find is not the case (see Figure 9). The rapid transition from radial (v_ϕ ∼ 0 km s⁻¹) orbits to circular ones (v_ϕ ∼ 200 km s⁻¹) is caused by the presence of accreted populations and the low-α disc, and is only seen when examining low-α populations. On the contrary, when looking solely at high-α populations – which should trace directly the spin up of the old proto-Galaxy to the high-α disc –, it is immediately clear that the relation in this diagram is much more gradual with respect to increasing metallicities.

Previous work has attempted to look at the transition between hot/radial orbits and cold/circular ones using this α-selection (Chandra et al. 2024). Thus, it is important that we compare our results to theirs. We argue that the reason these [M/H]−η 2D column-normalised histograms look strikingly different (especially for the high-α samples) from the study by Chandra et al. (2024) is because of the way the 2D histograms are plotted (both studies use the same data and similar α-separation curves). Chandra et al. (2024) column-normalises their histograms by amplitude (tracing the mode of the distribution), while we column-normalise by their sum (tracing the underlying PDF). Column-normalisation by amplitude traces the mode of the distribution, which makes the whole 2D histograms more noisy (as the distribution is scaled by a factor of the standard deviation of the curve). Tracing the mode also means that in a bimodal distribution across a large range of metallicities (like for the high-α sample), the mode switches between one and the other rapidly within a small bin size in metallicity. This could lead to the data manifesting a sharp increase from η ∼ 0 to η ∼ 1, as seen in Chandra et al. (2024), when in reality the data shows a more gradual increase in circularity with respect to metallicities (this work). Therefore, it is important to know how different normalisation methods can give rise to different interpretations of the underlying data and choose the normalisation method that best represents the science question that needs solving. In our case, as we aim to understand how the low-[M/H] regime transitions into the high-[M/H] regime for both high-/low-α populations, we choose to plot the sum column-normalised distribution. Furthermore, the differences with Chandra et al. (2024)’s results also arise from a small difference in the α-separation. These normalisation and α-selection differences and their implications are discussed in detail in Appendix B.

Even though orbital circularity depends on the adopted Galactic potential, it is position independent as opposed to azimuthal velocity. This makes it valuable to model the circularity evolution across metallicities. Due to the non-Gaussian nature of orbital circularity (as it abruptly cuts at −1 and 1 for a perfectly circular retrograde and prograde orbits respectively), the models described earlier in this work using azimuthal velocities cannot be directly applied to orbital circularities. Therefore, the underlying Gaussians are modelled as folded normal distributions (folded at −1 and 1). The modelling is only performed on high-α stars because they trace the transition of the protoGalaxy to rotation-supported high-α disc more clearly as seen in Figure 9, while the low-α stars can be simply described by two separate components similar to the [M/H]−v_ϕ plane.

The halo-like component has a fixed mean of 0 while varying standard deviation with a truncated normal prior with a mean of 0.4 and standard deviation of 0.2, with lower/upper limits restricted between 0.2 and 0.7 across the metallicity knots. We use the same metallicity range as in Section 4.2, with 12 and 6 metallicity knots (equidistant in linear space and downsampled between each bins to have the same number of stars) for relative fractions and orbital circularity means and standard deviation. The definition of a covariance matrix to enable the information between each components across the metallicity is unchanged in this circularity model, with the standard deviation on the mean of the circularity having a uniform prior between 0.1 and 0.9. The mean orbital circularity of the disc-like component is described by a multivariate normal distribution centred at 0.3 and the standard deviation is described by a half normal prior centred at 0.3. The final distribution of both the components are converted to a folded normal distribution to account for the folding at −1 and 1.

The model is run with 250 warm-up steps and 50 sample chains with the chains converged and r_hat less than 1.1 and effective samples above 30. We show the graphical model representation of this model in the right panel of Figure 5 and a list of model parameters and their functional forms are provided in the fourth column of Table 1. Figure 10 shows the converged parameters as a function of metallicity after the spline interpolation between the metallicity knots. We see that the orbital circularity of the disc-like component is steadily increasing with increasing metallicities over a large range of metallicities, [M/H] ∼−1.7 to −1 and is still steadily increasing up to [M/H] ∼−0.5. However, the interpretation of metal-rich stars ([M/H] > −1.0) is trickier given that the long asymmetric tail of the disc-like population is poorly fit due to the assumption of the underlying distribution as a simple (folded) Gaussian distribution. For the metal-poor stars, we can more confidently say that the spin up phase is more extended over a large range of metallicities from [M/H] ∼ 1.7 to −1. The relative fractional contribution from the spin up (disc-like) and halo-like population also supports the slower evolution of circularity with respect to metallicities, as opposed to the rapid spin up shown in the literature (see e.g. Chandra et al. 2024; Kurbatov et al. 2024). The major difference between the azimuthal velocity and orbital circularity evolution from the modelling perspective is (i) the metallicity at which the halo and spin up component crossover in the relative contribution, which is much more metal-rich in circularity ([M/H] ∼−1.2) compared to velocity ([M/H] ∼−1.7), and (ii) the mean orbital circularity is more circular (η ∼ 0.57) at lower metallicities ([M/H]<−1.5, also seen in the 1D histogram in Figure 9) than previously reported and more circular than what mean azimuthal velocities show at the same metallicities (v_ϕ ∼ 50 km s⁻¹). Both of these differences can be explained due to the fact that circularity has less position dependence than velocities (given that our sample has stars outside the solar neighbourhood; 49% of our stars are at heliocentric distances larger than 3 kpc). This is because the mean azimuthal velocity reduces close to the inner Galaxy when compared to solar neighbourhood, making the mean azimuthal velocity smaller in value. This affects the metal-poor stars more, as metal-poor stars are centrally concentrated (Rix et al. 2022), especially in the high-α sample. The latter difference could also arise from the fact that in the literature, the metal-poor high-α stars are not modelled with both an accreted and in situ population, whereas we can clearly see a bimodal distribution in the 1D histograms in Figure 9, justifying our choice of modelling the leftover accreted halo stars in the high-α regime along with the proto-Milky Way-like population. Therefore, the evolution of orbital circularity over increasing metallicities shows that the transition from a slowly rotating population to high-α disc is more gradual and not as rapid at [M/H] ∼−1. However, an important caveat to mention is that the orbital circularity represents how circular the orbit of a star is compared to the present-day disc orientation, which does not necessarily trace the circularity at formation.

Fig. 9 Top: column-normalised (by sum) 2D histogram of stars in the [M/H]−η plane (orbital circularity vs metallicity, top) and [M/H]−vϕ plane (azimuthal velocity vs metallicity, middle bottom) for all the stars (right), high-α selection (left), and low-α selection (midde). The running median track is shown as dashed black line and the 16th and 84th percentile tracks are shown as black lines in all panels. The running median tracks for the all stars and low-α panels look more like a step-function, while the high-α tracks are shallower, supporting the more gradual spin-up phase with respect to metallicities. (middle top) Probability density functions (PDF) of η in bins of [M/H] for all the stars (right), high-α selection (left), and low-α selection (middle). (bottom) Probability density functions (PDF) of vϕ in bins of [M/H] for all the stars (right), high-α selection (left), and low-α selection (middle). We note the second peak with low net spin in high-α panel slowly gains rotation with increasing [M/H] and spins-up into the old high-α disc in both middle top and bottom panels.

Fig. 10 High- α model with evolving means and standard deviation for orbital circularity of the disc-like component. Left: Evolution of orbital circularity, and orbital circularity dispersion as a function of metallicity, capturing a gradual transition of a chaotic state (predominantly made of proto-Galactic populations) to a rotationally supported disc-like motion in purple and the orbital circularity dispersion evolution of halo is shown in orange. Right: Fractional contribution of disc phase (proto-Galaxy to high-α disc) (purple) and halo-like (orange) GMM components as a function of metallicity. In all the panels, the scatter points show the location of the knots chosen to run our model. The bands are 1 σ uncertainties on the converged weights, circularity, and circularity dispersions. The fractions are computed for an evolving Gaussian for the disc phase and fixed Gaussian for the halo phase when the MCMC is converged.

6 Discussion

In this section, we discuss the summary of our results and the implications of the gradual spin-up phase across a large range in metallicities. We furthermore perform a simple GMM decomposition in bins of metallicities, to support the gradual spin-up phase scenario. Finally, we present the limitations and future scope of this work.

6.1 Summary of results

In this work, we have set out to model the azimuthal velocity and orbital circularity evolution of high-/low-α stars across metallicity space using Gaia XP element abundances, DR3 astrometry, and RVS radial velocites, to understand the transition phase of the proto-Galactic population to the high-α disc. At first, we model the conditional distribution P(v_ϕ |[M/H]) using a two-component (disc-like and halo-like) GMM with frozen means and standard deviations for the azimuthal velocities for the high- and low-α stars. We find the inferred posterior matches better with the data for low-α stars than high-α stars. The fractional contribution from disc-like and halo-like components look like a step function at [M/H] ∼−1 for low-α stars, while the high-α stars have a relatively gradual rise over increasing metallicities in the fractional contribution (Figure 6). Secondly, given that the high-α stars are not modelled well by the frozen means model, we model the conditional distribution P(v_ϕ |[M/H]) using a two-component (disc-like and halo-like) GMM with evolving means and standard deviations for the azimuthal velocities for the high-α stars. From this exercise, we see that both the mean azimuthal velocity and fractional contribution of the disc-like component is gradually increasing starting from a v_ϕ ∼ 50 km s⁻¹ over increasing metallicities (−1.7<[M/H]<−1.0). Thirdly, we perform the same exercise for high-α stars in orbital circularity versus metallicity plane, [M/H]–η, given that the orbital circularity is position independent, as opposed to azimuthal velocity. From this exercise, we also see a gradual rise in orbital circularity starting from an η ∼ 0.57 over increasing metallicities for the disc-like component.

Using different flavours of these mixture models, we have provided several lines of evidence that the metal-poor high-α disc increases its average azimuthal velocity, orbital circularity and rotational support gradually and monotonically across a wide range of [M/H], spanning approximately −1.7<[M/H]< −1. These data favour the scenario of a gradual spin-up of the metal-poor high-α disc (likely the proto-Galaxy) to a rotationally supported high-α disc at higher [M/H]. On the contrary, due to the superposition of the GES debris and the low-α disc in the low-α sample, the transition from metal-poor (halo) populations to metal-rich (disc) ones is much sharper, appearing almost like a step-function at [M/H] ∼−1. Due to the GES debris dominating the metal-poor sample for all stars in our dataset, this yields a similar profile when inspecting the Gaia XP sample without any α selection. Thus, our results highlight the importance of the [α/M] selection for studying the azimuthal velocity evolution of the old Milky Way disc and to avoid the GES debris. This also suggests that the proto-Galactic debris gained rotation gradually over increasing metallicities, eventually settling into a high-α disc and the disc formation was not as rapid or dramatic across a narrow range of metallicities as previously reported.

6.2 On the spin up of the Milky Way disc

A comprehensive understanding of disc evolution necessitates accounting for the complex interplay between stellar metallicity, α-element enhancement, age, and spatial distribution. Foundational models such as the monolithic collapse scenario proposed by Eggen et al. (1962) described the Milky Way forming rapidly from a single collapsing gas cloud. While this framework was influential, it cannot fully account for the detailed chemical and structural features observed today-particularly the clear dichotomy between the thin and thick disc components and the chemical bimodality in the [α/Fe]−[Fe/H] plane. To address these limitations, models like the two-infall scenario (Chiappini et al. 1997) introduced sequential gas accretion episodes, offering a more dynamic and layered view of disc formation. As highlighted in recent studies (e.g. Spitoni et al. 2019), combining chemical abundances with kinematic, spatial, and age information has been essential for revealing the Milky Way’s more intricate evolutionary pathways. This broader context informs our interpretation of disc spin-up mechanisms and their connection to the Galaxy’s formation history.

At the earliest cosmic times, the systems that seeded the formation of the Milky Way are conjectured to be manifested as a chaotic, unstructured, “proto-Galactic” population (e.g. Mowla et al. 2024), shaped by continual merging of building blocks that gave birth to the Galaxy’s ancient stars (Tumlinson 2010; El-Badry et al. 2018; Horta et al. 2024; Semenov et al. 2024; Xiang et al. 2024). Evidence of stellar populations associated with the proto-Milky Way have been proposed in the literature (Lucey et al. 2019; Arentsen et al. 2020; Reggiani et al. 2020; Horta et al. 2021; Belokurov & Kravtsov 2022; Ardern-Arentsen et al. 2024), with a clear picture now emerging thanks to the vast Gaia data (Rix et al. 2022). Expectations from cosmological simulations corroborate these observational findings (e.g. Horta et al. 2024; McCluskey et al. 2024; Semenov et al. 2024); these suggest that at present day, the remnants of proto-Milky Way fragments should present weak but systematic net rotation with respect to. the Galactic disc (v_ϕ ∼ 50 km s⁻¹), should primarily be concentrated towards the innermost Galactic regions (r ≲ 10 kpc), and should host stars with high [α/Fe] abundance patterns, similar to the old Milky Way disc.

Ascertaining the transition from a dispersion dominated population to a rotationally supported one is key to unravel the genesis of the Galactic disc. α-enhanced stars trace an early phase of the Galaxy’s chemical evolution, having formed before the onset of significant iron enrichment from Type Ia supernovae-typically within the first ∼ 1 Gyr (Tinsley 1979; Matteucci & Greggio 1986; Raiteri et al. 1996). Their elevated [α/Fe] ratios reflect nucleosynthetic contributions from short-lived massive stars via core-collapse supernovae, prior to the delayed iron input from Type Ia events. This short enrichment timescale implies that the observed changes in angular momentum or spin among these populations reflect relatively rapid dynamical processes in the high-α disc, rather than a gradual evolution in Galactic timescales (Haywood et al. 2013; Minchev et al. 2013; Mackereth et al. 2019). However, such short dynamical timescales can be manifested across a wide range of metallicities. As such, tracking the kinematic evolution of these in situ populations provides a direct view of the internal processes that shaped the early high-α disc component. Published age determinations for stars in the high-[α/Fe] disc population consistently indicate that these stars are predominantly old, typically older than ∼ 10 Gyr (Martig et al. 2016; Silva Aguirre et al. 2018; Imig et al. 2023). This narrow and ancient age distribution places strong constraints on the timescales for both chemical enrichment and dynamical evolution within this population. In particular, it implies that the increase in spin observed among the high-α stars must have occurred rapidly – within the first 1−2 Gyr of Galactic evolution – consistent with enrichment dominated by core-collapse supernovae and early dynamical settling of the high-α disc. Due to the inability to currently measure stellar ages precisely for metal-poor stars, metallicity ([M/H]), a much more reliable stellar parameter estimate to determine, is often used instead. While [M/H] is not a direct tracer of age, the reason why [M/H] is used is because, under the assumption that a stellar population chemically evolves in a progressive manner, stars that form later are more enriched in metals than their previous generations; in other words, old populations tend to be metal-poor while younger populations tend to be metal-rich. Thus, [M/H] has been shown to be a useful metric to study the evolution of kinematic samples (e.g. Belokurov & Kravtsov 2022; Xiang & Rix 2022; Xiang et al. 2024). This is especially the case in the current era of large-scale stellar surveys. For example, the vast Gaia dataset of [M/H] (Andrae et al. 2023a) and now also [α/M] (Li et al. 2024) for over two million stars has given an unprecedented (qualitative) view on the transition from the proto-Milky Way to the high-α disc (e.g. Chandra et al. 2024; Zhang et al. 2024). However, most of these studies either: (i) examine the profile of azimuthal velocity, v_ϕ, and/or orbital circularity, η, with [M/H] agglomerating stars with different [α/M]. This is known to be an issue due to the inclusion of stellar halo populations in the sample, that have different [M/H] and kinematic distributions; (ii) model the transition of dispersion dominated to rotationally dominated populations in the [M/H]–velocities plane using disconnected GMMs (see Section 6.3), making it difficult to link components across [M/H]; (iii) do not provide a quantitative measure of how the Milky Way transitioned from being dispersion dominated to rotationally dominated.

The correlation between azimuthal velocity and metallicity in the thick disc has been established in several earlier studies. Using kinematically selected samples at high Galactic latitudes, Spagna et al. (2010) first identified a steady increase in rotation velocity with [Fe/H] among thick disc stars. This trend was later confirmed for chemically defined high-[α/Fe] populations by Kordopatis et al. (2017), who found a comparable v_ϕ−[Fe/H] gradient. Although these earlier analyses were constrained by smaller sample sizes and narrower metallicity ranges, they provided foundational evidence of a kinematic-chemical link within the thick disc. In addition, the presence of a metal-poor thick disc ([Fe/H] ≤−1.0 dex) is now well established (e.g. Ruchti et al. 2011; Kordopatis et al. 2011; Bonaca et al. 2017; Di Matteo et al. 2019), extending the high-α population to lower metallicities than those accessible in the last decades. Belokurov & Kravtsov (2022) used high-resolution spectroscopic data from the APOGEE survey to study the evolution of stars around the solar neighborhood they deem in situ in the metallicity-v_tan plane; these authors found that the Milky Way “spun up” rapidly between −1.3<[Fe/H]<−0.9 (see also Zhang et al. (2024) for a similar result using Gaia XP metallicities). Along those lines, Chandra et al. (2024) used the available [α/Fe] for the Gaia XP sample to study the profile of [M/H]−η, finding that the high-α disc emerges rapidly at [M/H] ∼−1.

Our analysis builds on this previous works by tracing these trends across a broader metallicity baseline and with more robust statistics, enabled by the larger, well-characterised stellar samples now available with Gaia DR3. In this work, we have found that when performing a similar [α/M] cut to Chandra et al. (2024), our transition from proto-Milky Way populations to the high-α disc (Figure 9) shows a much more gradual increase in circularity with respect to metallicities when compared to these previous works. In more detail, both when examining and modelling the [M/H]−v_ϕ and [M/H]−η planes, we have found that the spin up of the old Milky Way disc occurs across a much wider range in [M/H], namely between −1.7<[M/H]<−1. The reason for the discrepancy between our results and previous work is because: (i) Belokurov & Kravtsov (2022) do not go down to lower metallicities ([M/H]<−1.5) and have lower number statistics compared to our work. The 1D histograms of azimuthal velocities in bins of [M/H] in Belokurov & Kravtsov (2022) already shows preliminary hints of a gradual spin up (see their Figure 5); (ii) Zhang et al. (2024) perform a traditional GMM without any α-selection, resulting in the GES debris driving the rapid spin up inference at a narrow range in metallicities; (iii) Chandra et al. (2024) shows a qualitative view of [M/H]−η plane, while their high-α sample is contaminated by GES debris more than ours, and their column normalisation is performed by amplitude instead of sum. Therefore, the results presented in this work, in comparison with the recent literature, shows that the proto-Galaxy to high-α disc transition has likely been gradual over increasing metallicities and not as dramatic over a narrow range of metallicities as previously reported. While direct comparisons with high-redshift measurements of V/σ (e.g. Wisnioski et al. 2015) are limited by differences in methodology and redshift coverage, our results (V/σ ≍ 6 and 10 for high- and low-α disc at [M/H]=0 respectively) offer a z=0 anchor point for the rotational support of chemically distinct disc populations in the Milky Way.

We finish this section by raising a speculative, yet interesting point. In Figure 7 and Figure 10, we find that the profile corresponding to the “spin up” displays a bump in the profile at [M/H] ≍−0.8. We have checked that this bump is not artificially caused by our fitting methods. Interestingly, the location of this bump in [M/H] coincides with the end of the low-α accreted sample (primarily the GES merger debris). It is postulated that this population is the debris from a major merger in the Milky Way that dominates the local stellar halo (6 ≲ r ≲ 30 kpc: e.g. Belokurov et al. 2018; Helmi et al. 2018; Koppelman et al. 2018). This merger possibly catalyzed the formation of the in situ halo, or “hot high-α disc” comprising stars on halo-like orbits with chemistry similar to the high-α disc (Di Matteo et al. 2019; Bonaca et al. 2020; Belokurov et al. 2020). Thus, we speculate that this bump could be a sign of the impact of the GES with the old high-α Milky Way thick disc.

6.3 Gaussian mixture model on 3D velocities in bins of metallicities

In this subsection, we perform the traditional N-component GMM decomposition of high-α, low-α, and all stars subsamples in bins of [M/H] on the 3D velocity space in cyclindrical coordinates −v_r, v_ϕ, and v_z. We set out to perform this exercise to understand how this method compares to our method presented in Section 4.2.

Typically, within any GMM framework, each sub-population is characterised by three key parameters. The first is a weighting factor, which reflects the relative proportion of stars belonging to that particular group. The second is the mean velocity, representing the first moment of the velocity distribution, and the average velocity along each of the three dimensions for stars within that group. Finally, the second moment of the velocity, the velocity dispersion, captures the variation in velocities around the average for stars within that group. We use the same method outlined in Zhang et al. (2024) for the GMM decomposition, to be able to compare the effects of α-separation in this type of GMM decomposition. To perform the GMM analysis, we utilised a software package named pyGMMis (Melchior & Goulding 2016). This software incorporates the sophisticated “Extreme Deconvolution” technique developed by Bovy et al. (2011) to account for the inherent uncertainties in stellar velocity measurements and the covariances between the input parameters.

Our analysis focuses specifically on metal-poor stars. We restricted the GMM decomposition to stars with a metallicity range of −2.5 to −0.7, expressed as [M/H]. Below this we are hampered by low number statistics and above these metallicities we are affected by the non-Gaussian nature of the low-α and high-α discs distribution in velocity space that could compromise the accuracy of the GMM results. The [M/H] bins are spaced equidistantly on a log scale with seven(six) bins for the high-(low-) α samples. The all stars sample is fit using the same six bins as in the low-α sample⁶. Zhang et al. (2024) excluded stars located further than |z|>2.5 kpc from the Galactic plane, which is a cut that we do not use to be able to study the entire underlying data. We argue that as we are modelling the halo and disc populations separately, this cut is unnecessary. For each metallicity bin, we applied the GMM to the three-dimensional velocity space defined by the cylindrical galactocentric coordinates: radial velocity (v_r), azimuthal velocity (v_ϕ), and vertical velocity (v_z). Moreover, each star’s measurement uncertainty was incorporated as a covariance matrix with the diagonals representing the variances for each velocity component and the off-diagonals representing the covariance between a pair of velocities. To determine the optimal number of Gaussian components at each metallicity bin, we employed a standard approach based on the Bayesian Information Criterion (BIC), defined in the same way as Zhang et al. (2024). Here, lower BIC values indicate a better fit, balancing model complexity with data fidelity. To mitigate the risk of getting stuck in a local minima during the optimisation (sub-optimal solutions), we repeated the GMM fitting process 100 times with different starting conditions and recorded the BIC value for each trial. The fit with the lowest BIC value (in the median BIC curve) was considered the optimal solution, suggesting a high likelihood of finding the globally best fit. Increasing the complexity of the model by adding an extra component, while keeping the log-likelihood unchanged, raises the BIC value by ∼ 100. Given this order of magnitude, some N-component GMMs have very similar BIC values, especially in relatively metal-rich regimes, where the rotating component is fit with multiple Gaussians instead of 1. When BIC values are similar, we prefer models with fewer components, as they offer a more straightforward physical interpretation. This occurs for at least 1 bin in high-α and 2 bins in low-α and all stars sub-samples. We do not compute the uncertainty in the GMM parameters as Zhang et al. (2024) found them to be generally around ∼ 0.1 km s⁻¹. Thus, uncertainties in the fits can be treated as negligible.

In Figure 11, we show the ratio of azimuthal velocity to vertical velocity dispersion (V_ϕ/σ_z), which is commonly used as a measure of how rotation-supported or “disc-like” a sample of stars is. We restrict our analysis to the most prograde GMM component (largest mean v_ϕ) per [M/H] bin to trace the evolution of the Milky Way’s disc; these are also the component with an almost zero mean v_r and v_z. In addition, to ensure that our model captures the halo and disc components well, we verified that the other GMM components the model finds fit for a halo-like sub-structure (V_ϕ/σ_z ⋘ 1). However, as we are only interested in the evolution of the disc-like GMM, we restrict our analysis to the most prograde GMM component in this work in Figure 11.

From Figure 11, we can see that the high-α stars gradually increase in their V_ϕ/σ_z (in purple) towards higher metallicity bins while low-α (in orange) and all stars (in black) have a more rapid increase in V_ϕ/σ_z between two [M/H] bins, ∼−1.3 and ∼−1. The size of the scatter points are directly proportional to their relative weights in each bin. Therefore, our results suggest that not accounting for the α-separation could make the spin-up seem more rapid and exponential over a small range of metallicities. In our case, as we have been able to distinguish high-/low-α populations, we find that the high-α prograde GMMs favour a gradual spin-up of a metal-poor high-α population (likely the proto-Galaxy) to a metal-rich high-α disc over increasing metallicities. The GMMs from Zhang et al. (2024) show a two-component fit in the most metal-poor bin while our model fits only one component. This is due to the fact that our all stars sample is composed of all stars with reliable metallicities from Andrae et al. (2023b), while having an α estimate from Li et al. (2024), whereas Zhang et al. (2024) only used the metallicity estimates and therefore, we end up having far fewer VMP stars than theirs. Zhang et al. (2024)’s results also suggest a more rapid growth in V_ϕ/σ_ϕ over a shorter range in metallicity that is slightly more metal-poor than ours (−1.7<[M/H] <−1.3), which is explained by the difference in bin sizes and the scale height |z|<2.5 kpc cut. GMM decomposition can be highly dependent on the choice of the [M/H] bins. However, this GMM decomposition is performed in order to make a qualitative comparison between high- and low-α samples, and is not a one-to-one comparison to the literature results. The evolution of the velocity components and its dispersion focused on the high-α stars tracing the evolution of the proto-Galaxy is shown in Appendix C.

In summary, our results reveal that the high-α subsample shows a gradual rise in V_ϕ/σ_z, favouring a gradual spin-up phase for the Milky Way’s high-α disc over increasing metallicities.

Fig. 11 Vϕ/σz ratio (azimuthal velocity to vertical velocity dispersion) for the most prograde Gaussian component in each metallicity bin (spaced evenly on a log scale) from our GMM runs in Section 6.3. The square size indicates the relative weight of the most prograde Gaussian in each metallicity bin. This ratio is a unitless quantity indicating the level of rotational support within a disc, essentially quantifying how dynamical cold it is at z=0. The dashed line at Vϕ/σz=1 marks the boundary between pressure-supported (halo-like) vs rotation-dominated kinematics (disc-like).

6.4 Limitations and future scope

In this work, we use an α-separation on Gaia XP based [α/M] and [M/H] abundances to separate high-/low-α stars and use the high-α stars to trace the azimuthal velocity (v_ϕ) evolution of the old Milky Way disc over a range of metallicities (−2.5< [M/H]<0.1). The high-α selection implemented in this work is a simple piece-wise function (as described by equations (1) and (2)), and is supposed to select only in situ stars. However, at lower metallicities [M/H]<−1.5, accreted merger populations (e.g. Heracles, the high-alpha tail of the GES and possibly other merger remnants) overlap, so a simple α-cut can no longer establish a clean separation between in situ versus accreted population purely. This is modelled with our mixture model with evolving mean velocities and velocity dispersion for the spin-up phase and a background halo model that is isotropic as described in Sections 4.2 and 5. However, this model assumes a Gaussian distribution of azimuthal velocity at any given metallicity, which is not strictly true for the high-α disc, which has an asymmetric drift (long tail towards lower v_ϕ, see Figure 9). Because of this, the model tries to inflate the halo velocity dispersion at higher metallicities to fit this tail of the high-α disc which is nonphysical. This also makes the model yield strong systematic residuals when compared to the data in the metal-rich end, where the high-α disc displays a strong prograde profile. This is a consequence of a simple Gaussian assumption that the model is based on. In the future, it would be possible to extend this model by using a more sophisticated dynamically motivated distribution function that represents more accurately the Galaxy’s disc/halo in order to fully measure the evolution of velocities from the chaotic proto-Galactic state to an ordered high-α disc state with respect to metallicities. In doing so, we could also use the total velocity distributions to constrain the enclosed mass and in turn, measure the mass growth of the Milky Way from a proto-Galactic state at high-redshift to the old high-α disc state as a function of increasing metallicity. We reserve this exploration for a future work.

It is also important to note that the measured mean velocities and velocity dispersions with respect to decreasing metallicities are all present-day velocities and not the velocity they had at formation. From cosmological simulations, there have been clues that the Galaxy must have undergone post-formation dynamical heating which could change the net rotation that we see at present-day (McCluskey et al. 2024; Horta et al. 2024).

One other factor to keep in mind is the ability to precisely capture the subtle trends in [α/Fe] with our high-α sample, given our α separation. Recent studies such as Belokurov & Kravtsov (2022) and Conroy et al. (2022) used APOGEE and H 3 surveys, respectively, to show that the [α/Fe] ratio gradually declines between −3<[Fe/H]<−1.3, and then rises to a higher value to meet the high-α disc population. It is very likely that we are not catching this dip in [α/Fe] with our simple α-separation. However, this should not affect the main conclusions of this work.

It is also important to recognise that the high-[α/Fe] disc population is not a single-age cohort, but rather spans a range of formation times and chemical enrichment levels. In particular, stars with the highest [α/Fe] values (e.g. ≥ 0.3 dex) likely formed very early – within the first few hundred million years of Galactic evolution – while others in the high-α sequence formed later, albeit still early on Galactic timescales (Matteucci & Brocato 1990; Chiappini et al. 1997; Hasselquist et al. 2019). While our analysis focuses on the median azimuthal velocity and orbital circularity evolution of the high-α disc as a whole, disentangling how internal structure in [α/Fe] within this population contributes to the observed trends is an avenue for future exploration.

7 Conclusions and outlook

In this work, we set out to model the azimuthal velocity evolution of high- and low-α stars across metallicity space using Gaia XP element abundances, DR3 astrometry, and RVS radial velocities. By employing various mixture models, we have uncovered several lines of evidence that provide new insights into the kinematics of the metal-poor high-α disc and its evolution within the Milky Way.

Our analysis reveals that the metal-poor high-α disc shows a gradual and monotonic increase in its average azimuthal velocity distribution over a broad range of [M/H], spanning approximately −1.7<[M/H]<−1. This finding supports the scenario of a gradual spin-up of the metal-poor high-α disc, likely representing the proto-Galaxy, evolving into a rotationally supported high-α disc as [M/H] increases. This gradual transition underscores the dynamic evolution of the early Milky Way over a large range in metallicities, reflecting a continuous and progressive change in the rotational characteristics of the metal-poor high-α stellar population.

In contrast, the low-α sample presents a different picture due to the superposition of the GES and other accreted debris along with the low-α disc. The transition from metal-poor (halo) populations to metal-rich (disc) populations in this sample is much sharper, resembling a step-function at [M/H] ∼−1. The dominance of the GES debris in the metal-poor sample for all stars in our dataset results in a similar profile when inspecting the Gaia XP sample without any α-selection. This sharp transition highlights the distinct formation and evolutionary histories of these stellar populations compared to the more gradual evolution in rotation across metallicities observed in the high-α disc.

Our results emphasise the critical importance of [α/M] selection for studying the azimuthal velocity (or the orbital circularity) evolution of the old Milky Way disc. By distinguishing between high- and low-α populations, we can better understand the complex interplay between different stellar populations and their contributions to the overall dynamics of the Galaxy. This distinction helps us to disentangle the effects of accreted debris and in situ star formation, thus providing a clearer picture of the Milky Way’s formation history and its subsequent evolution. However, we note that this separation, while useful, does not lead to a purely in situ population.

In conclusion, our study provides strong evidence for the gradual spin-up of the metal-poor high-α disc, suggesting a continuous and progressive evolution from the proto-Galaxy to a rotationally supported high-α disc over increasing metallicities, which has been recently reported in Horta & Schiavon (2025) using APOGEE-Gaia data’s detailed chemical abundances and velocity distributions. The sharp transition observed in the low-α sample underscores the impact of accreted material on the Galaxy’s dynamical structure. These findings contribute to our understanding of the early evolutionary processes of the Milky Way and highlight the need for detailed chemical and kinematic studies to unravel the complex history of our Galaxy. With the advent of more sophisticated distribution functions and increased availability of abundances for a larger number of stars, we be better equipped to perform detailed analyses of in situ stars, further unraveling the intricate history of the Milky Way and refining our understanding of its earliest formation and evolution.

Acknowledgements

We express our gratitude to the reviewer for dedicating their valuable time and providing insightful comments that significantly enhanced the quality of our manuscript. This work was developed for the CCA Pre-Doctoral Program in 2023 at the Flatiron Institute. We thank them for their generous support. AV thanks Vasily Belokurov for suggesting the test using high-resolution spectroscopic surveys. AV also thanks Ewoud Wempe for the helpful discussions on Gaussian processes, Tom Callingham for the helpful discussions on orbital circularity, Amina Helmi for the useful comments on this work, that helped improve our model, and Alis Deason for the insightful conversations on this topic. AV gratefully acknowledges support from the Canadian Institute for Theoretical Astrophysics (CITA) through a CITA National Fellowship and the International Astronomical Union (IAU) and the Gruber Foundation through a IAU Gruber Fellowship. DH was supported by the UKRI Science and Technology Facilities Council under project 101148371 as a Marie Curie Research Fellowship. ES acknowledges funding through VIDI grant “Pushing Galactic Archaeology to its limits” (with project number VI.Vidi.193.093) which is funded by the Dutch Research Council (NWO). This work has been partially supported by a Spinoza Prize from NWO. 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. AV also thanks the availability of the following packages and tools that made this work possible: vaex (Breddels & Veljanoski 2018), pandas (Reback et al. 2022), astropy (Astropy Collaboration 2022), NumPy (Oliphant 2006; Van Der Walt et al. 2011), SciPy (Jones et al. 2001), matplotlib (Hunter 2007), seaborn (Waskom et al. 2016), agama (Vasiliev 2019), gala (Price-Whelan 2017), galpy (Kluyver et al. 2016), healpy (Zonca et al. 2019), gaiadr3-zeropoint (Lindegren et al. 2021), jax (Schoenholz & Cubuk 2019), numpyro (Phan et al. 2019), scikit-learn (Pedregosa et al. 2011), pyGMMis (Melchior & Goulding 2016), JupyterLab (Kluyver et al. 2016), and topcat (Taylor 2018).

Appendix A [M/H]−v_ϕ plane towards and away from the inner Galaxy

Examining the [M/H]−v_ϕ plane for low-α stars in the APOGEE (Figure 4) data, we can see that the halo-like (isotropic) population at low metallicity is almost fully disconnected from the higher metallicity disc-like component. However, this feature is not as immediately clear when inspecting the same stellar populations in the Gaia XP (Figure 3) data (there appears to be more of a connection between the halo-like and disc-like population for low-α stars). In this Appendix, we set out to investigate if this difference is due to the sample selection difference between APOGEE and Gaia, that could lead to more high-α stars contaminating our low-α star sample, by looking at the difference between the [M/H]−v_ϕ plane towards and away from the Galactic centre.

One possible reason for the discrepancy is that proto-Galactic fragments tend to be more centrally concentrated (Horta et al. 2024). As APOGEE is a near-infrared survey that can penetrate through prevalent dust extinction in the central regions of the Galaxy more easily, it is likely that we are probing the proto-Galaxy better with APOGEE than with the (optical) Gaia survey. Thus, if our α cut was designed well, the low-α star sample should not be as centrally concentrated as the high-α population, that hosts both more centrally concentrated (high-α) disc stars (Hayden et al. 2015; Imig et al. 2023) and proto-Galactic populations (Horta et al. 2024).

Figure A.1 shows the [M/H]−v_ϕ plane for all, low-α, and high-α stars towards Galactic longitudes in the direction of the Milky Way Centre (i.e. |ℓ|<30^°, top panels), and towards Galactic longitudes away from the Galactic Centre (i.e. the Galactic anti-centre, |ℓ|>30^°, bottom panels) for our Gaia XP sample. For low-α stars, we can see that the halo-like population at lower metallicities is fully disjoint from the disc-like population at higher metallicities for stars towards the Galactic anti-centre (similar to what is seen in APOGEE, Figure 4), whereas the transition in v_ϕ across different metallicities is smoother towards the Galactic centre, similar to high-α median tracks. We reason that this is because, when looking towards the Galactic anti-centre, the two dominant populations contributing to the low-α regime are: 1) the (outer) low-α disc that is highly rotating; 2) the debris from the GES merger, which are highly radial. This also affects the all stars panel for stars towards the Galactic centre, as seen in the top panels of Figure A.1.

This leads us to infer that there exists a slowly rotating in situ proto-Galaxy-like population contaminating our low-α sample that is centrally concentrated, which is likely not present in the APOGEE stars. This arises mainly due to the higher abundance precision in APOGEE (high-resolution spectroscopic survey) compared to Gaia-XP inferred α-abundances.

In Figure A.2, we show the on-sky distribution of stars in Galactic coordinates in Gaia XP sample with a healpix level of 7 (top panels) and APOGEE sample with a healpix level of 4 (bottom panels) for both high- and low-α selection between the metallicities of −0.7 and −1.3. We choose this metallicity range because this is the range between which we see the difference in the v_ϕ median tracks between low-α stars towards and away from the Galactic centre as shown in Figure A.1. In both APOGEE and Gaia XP high-α stars, we see centrally concentrated distribution of metal-poor stars, reminiscent of the ‘poor old heart’ (Rix et al. 2022) also known as the proto-Galactic in situ population. However, we also see a higher density of stars towards the inner Galaxy in the Gaia XP low-α panel, which is not as clear an overdensity in the APOGEE low-α panel. Within the inner Galaxy (|ℓ|<30^° and |b|<30^°), the high-α stars are overdense compared to the low-α stars by a factor of 8 for APOGEE sample while the high-α stars are overdense compared to the low-α stars by a factor of 4 for the Gaia XP sample. Therefore, this mismatch between APOGEE and Gaia XP is not just due to lower number statistics in APOGEE. This difference could be due to unreliable α estimates and simple definition of α-separation. Therefore, in situ versus accreted separation using Gaia XP α estimates is not as reliable and the low-α selection using Gaia-XP sample is more contaminated than the APOGEE sample, which is most likely the reason for a shallower step function in the low-α v_ϕ median tracks in Gaia XP (bottom left panel of Figure 3) compared to the steeper one in APOGEE (middle panel of Figure 4). However, we do not expect our conclusions to be affected by this contamination as our high-α (mostly in situ) selection is still pure.

Appendix B Comparison of orbital circularity versus metallicity with Chandra et al. (2024) results

In this section, we compare the evolution of orbital circularity versus metallicity in our work with the results from Chandra et al. (2024). Both the works use the same input catalogue based on Gaia XP spectra. The main differences between our approaches are the difference in the α-separation, and the way the column normalisation is performed. In Figure B.1 top panels, we compare the column-normalised by amplitude 2D histograms of orbital circularity versus metallicity and their corresponding η, and v_ϕ 1D histograms for Chandra et al. (2024) α-selection and the α-selection described in this work (equations 1 and 2). In all the panels, the 16th, 50th, and 84th percentile tracks are shown as black lines, mean tracks are shown as black dashed lines and mode tracks are shown as grey dashed lines. It is important to note that the mode tracks trace the underlying distribution the best when the column normalisation is calculated by the amplitude of the distribution. In Figure B.1 bottom panels, we compare the column-normalised by sum 2D histograms of orbital circularity versus metallicity and their corresponding η, and v_ϕ 1D histograms for Chandra et al. (2024) α-selection and the α-selection described in this work (equations 1 and 2). This is equivalent to what is presented in the rest of this paper (tracing the PDF of each distribution, equivalent to normalising such that the area under the curve is equal to 1). In all the panels, the 16th, 50th, and 84th percentile tracks are shown as black lines.

There are two main inferences that can be made from this comparison, as discussed below:

• The α-selection described in this work is more efficient in removing the last major merger (accreted GES) from the high-α selection (in situ equivalent), better than the α-selection described by Chandra et al. (2024). Therefore, our high-α stars are cleaner than the high-α stars from Chandra et al. (2024).

• Column normalisation of 2D histograms can be done in many ways. From our comparison, we conclude that the column normalisation is scaled by the amplitude of the distribution in Chandra et al. (2024), whereas in this work, we column normalise the histograms by the sum of the distribution. Column normalisation by amplitude scales the peak of each 1D histogram to 1.0 whereas column normalisation by sum scales the sum of the histogram (equivalent to integral under the histogram curve) to be equal to 1.0 thereby tracing the probability distribution function of each 1D histogram. Column normalisation by amplitude traces the peak of each 1D histogram with no easily interpretable connection between each 1D histograms, and therefore produces a noisy view of the underlying data. This can be seen by the noisy streaks in the 2D histograms and the noisy mode tracks in the top panels of Figure B.1. Column normalisation by amplitude traces the mode of the distribution, which is also tricky in case of bimodal distributions. The metal-poor end of our underlying azimuthal velocity is bimodal due to a slowly rotating proto-Galactic population and the high-α remnants of accreted stellar systems (mostly isotropic). The mode simply traces the peak of the distribution and therefore looks like a step function when one Gaussian dominates over the other in relative weights. This affects the high-α population the most as the low-α stars are already two almost disjoint distributions in metallicities (thin disc at higher metallicities and accreted halo at lower metallicities). Therefore, we emphasise that column normalisation by sum is more appropriate to understand the underlying distribution of orbital circularity or azimuthal velocity as it traces the PDF of the distribution and not just the peak/the mode.

Fig. A.1 Column-normalised (by sum) 2D histogram of stars in the [M/H]−vϕ plane for all the stars (right), high-α selection (left), and low-α selection (middle) towards the Galactic centre (top panels) and away from the Galactic centre (anti-centre, bottom panels). The running median track is shown as dashed black line and the 16th and 84th percentile tracks are shown as black lines in all panels. We can see that the low-α stars are not fully two separate populations towards the inner Galaxy due to contamination from the proto-Galaxy in the low-α end, that creates a small connection between the thin disc and halo populations in the overall [M/H]−vϕ plane for low-α stars as seen in Figure 3. Low- α stars towards the Galactic anti-centre show a cleaner, two separate population of thin disc and halo stars as expected. This effect is also reflected in the all stars panel. This effect is minimal in the APOGEE α-selection, most likely due to the high-resolution α-measurements.

These two reasons together explain why one would interpret the spin up to be rapid and drastic, whereas the underlying distribution is slowly gaining rotation over a wide range of metallicities.

In order to better understand the difference in the α-separation between our work and Chandra et al. (2024), we also show row normalised 2D histograms of [M/H]–η space for high-α stars in this work and high-α stars from Chandra et al. (2024) in the top and bottom panels of Figure B.2. Because of row-normalisation, higher metallicity stars are more highlighted in these figures. In both these panels, we can see the highly rotating high-α disc dominating higher metallicities ([M/H] > –0.4). Below these metallicities, we see high-α stars with a broad range of circularities, even down to retrograde orbits, but have metallicities that are representative of high-α disc. The most plausible origin for these stars are that they we born in the old high-α disc and got kicked-up into halo-like orbits by the last major merger, GES (Bonaca et al. 2017; Helmi et al. 2018; Belokurov et al. 2020). The most interesting difference between our high-α stars and that of Chandra et al. (2024) is that their high-α stars have larger number of stars in isotropic orbits around [M/H] ∼− 1.2, reminiscent of the GES merger. This is almost absent in our high-α stars. This leads us to believe that our high-α star selection is purer than the one in Chandra et al. (2024). It is easier for the eye to trace the excess of retrograde low-metallicity stars in both the panels, but these are only highlighted due to the row-normalisation (because retrograde stars are almost only present due to halo accretion events in the lower-metallicity end) and in reality, there are far fewer of these stars. However, we still see a population of accreted stars in our high-α selection (to a much lower extent than using Chandra et al. (2024) α-separation), which can be attributed to the evolution of any stellar system that has a high-α low-metallicity tail which cannot simply be removed using a simple α-separation. We model this population along with the evolution of high-α disc in Sections 4.2 and 5.

Fig. A.2 Logarithmic Galactic map of our sample of stars with high-α (left) and low-α (right) selections for the all-sky Gaia XP sample (top) and APOGEE DR17 sample (bottom), in the metallicity regime where the gradual spin-up in high-α vs the step function due to two different stellar populations in higher and lower metallicities in low-α is reflected the most: −1.3<[M/H]<−0.7. In the Gaia XP sample, we can see that the poor old heart or the proto-Galaxy is still present (contamination) in low-α towards the inner Galaxy, however, to a smaller extent than in high-α. This is not the case for APOGEE stars which have a much cleaner low-α population towards the inner Galaxy with the bulk of proto-Galaxy towards the inner Galaxy in the high-α Galactic map.

Appendix C Velocity evolution versus metallicity for high-α stars used as a cosmic clock

In this section, we show the evolution of the most prograde GMM component (based on the GMM runs explained in Section 6.3) in its different velocity components for high-α stars. In Figure C.1, we see the evolution of (V_r, V_ϕ, V_z) in the top panel, the evolution in their velocity dispersion (σ_r, σ_ϕ, σ_z) in the middle panel and their rotational support (V_ϕ/σ_z or V_ϕ/σ_tot) in the bottom panel across different [M/H] bins with the respective component names next to each curve. It is important to note that these velocities and their dispersions trace the present dynamical evolution and not those at formation. Using FIRE-2 simulations, McCluskey et al. (2024) shows that the rotational velocities can increase now compared to at formation in the pre-disc era due to stars being torqued into rotational orbits as the disc settles, and the rotational velocities can decrease now compared to at formation in the late-disc era due to dynamical heating. In case of rotational velocity dispersion, it also does not directly reflect the formation history, as it monotonically increases due to post-formation dynamical heating, adding to the velocity dispersion at formation. Therefore, even though metal-poor stars trace old stellar populations, the velocities do not simply trace the formation velocities. However, this simple analysis of velocity evolution can give us an idea of the evolution of the high-α disc as we see it now.

In the top panel of Figure C.1, we see that the radial velocity and vertical velocity is almost close to zero, with the rotational velocity (and the total velocity) increasing slowly with increasing metallicities. This is reminiscent of a slowly rotating proto-Galactic population that settles into a high-α disc at higher metallicities. We see that this spinning-up phase is more gradual over increasing metallicities than previously reported. However, if this gradual spin-up comes from the time of its formation or due to post-formation heating is an open question. In the middle panel of Figure C.1, we see the trends for velocity dispersions in all three directions decreasing with increasing metallicities, as the high-α disc begin to settle. In the bottom panel of Figure C.1, we see the evolution of rotational support as a function of metallicity. McCluskey et al. (2024) show that in late- and early-disc era, both v_ϕ/σ_z and v_ϕ/σ_tot decrease by a factor of 2 between formation and present state, due to post-formation dynamical heating that increases their velocity dispersions. In the bottom panel of Figure C.1, we see the evolution of rotational support gradually increasing with increasing metallicity, also reminiscent of a proto-Galactic population slowly gaining rotation and settling into the high-α disc. The metallicity at which this population becomes more rotation supported (‘discy’) is difficult to infer, as it differ between σ_z and σ_tot, and also that the overall V_ϕ/σ_tot reduces at present when compared to what it was at formation. Therefore, our main conclusion from these velocity evolution curves is that the spin-up phase of a slowly rotating proto-Galaxy is gradual across a large range of metallicities and not as rapid as previously reported.

Fig. B.1 Orbital circularity vs metallicity for high-α, low-α, and all stars, and 1D histogram of azimuthal velocity and orbital circularity at different metallicity bins for column-normalisation by amplitude (top panels) and by sum (bottom panels) using Chandra et al. (2024) α-selection (left panels) vs the α-selection implemented in this work (right panels). Running median (black line), mean (black dashed line), mode (grey dashed line), 16th and 84th percentile (black lines) tracks are overlaid for the 2D histograms with column normalisation by amplitude. We can see that the mean and median tracks follow each other closely except at higher metallicities where the non-Gaussianity of thin and thick discs dominates, whereas the mode tracks follow the peaks of the background and is very noisy and resembles the step function behaviour seen by Chandra et al. (2024), even for high-α stars where we see a gradual spin-up (across increasing [M/H]) using column normalisation by sum. Running median, 16th and 84th percentile tracks are overlaid as black lines for the 2D histograms with column normalisation by sum. We also see that our α-selection has a cleaner high-α sample and isolates the bulk of GES into the low-α sample.

Fig. B.2 Row-normalised conditional metallicity distribution for high-α stars across orbital circularities for our α-selection (top) and Chandra et al. (2024) α-selection (bottom). The high-α disc, in situ halo (hot high-α disc), and accreted components can be seen in both panels. However, the contamination from the accreted GES merger is stronger in the Chandra et al. (2024) α-selection than the α-selection proposed in this work.

Fig. C.1 Mean velocities (top), velocity dispersions (middle) and Vϕ/σz ratio (bottom) showing how rotationally supported the stars are, for the most prograde Gaussian component in each metallicity bin from our GMM runs for high-α stars.

References

All Tables

All Figures

	Fig. 1 Logarithmic density of [α/M] vs [M/H]. The purple band represents the high- and low-α sequence separation defined in this work (see text for details). Stars in the purple band are excluded. The bulk of accreted last major merger (GES) is primarily restricted to the low-α population with our selection.
In the text

Fig. 2 Distribution of stars in cylindrical Galactic coordinates (R–z plane) colour-coded by their mean metallicities for all the stars in our sample (top), high-α selection (bottom left), and low-α selection (bottom right). We can see that dust near the Milky Way’s midplane significantly impacts our survey selection. We see a sharp negative metallicity gradient with respect to height above the disc plane in all stars and low-α stars. High- α stars have a shallower negative metallicity gradient.

In the text

Fig. 3 Column-normalised (by sum) 2D histogram of stars in the [M/H]–vϕ plane (azimuthal velocity vs metallicity) for all the stars (top), high-α selection (bottom left), and low-α selection (bottom right). The running median track is shown as dashed black line and the 16th and 84th percentile tracks are shown as black lines in all panels. The running median tracks for the low-α panel look more like a step-function, while high-α tracks is shallower, supporting the gradual monotonic increase, that can be interpreted as a gradual spin-up from old proto-Galactic populations to the present day high-α disc with respect to metallicities.

In the text

Fig. 4 Column-normalised (by sum) 2D histogram of stars in the [M/H]−vϕ plane for APOGEE DR 17 giants (with the same high- and low-α separation as shown in Figure 1) for all the stars (right), high-α selection (left), and low-α selection (middle). The running median track is shown as dashed black line and the 16th and 84th percentile tracks are shown as black lines in all panels. The median, 16th and 84th percentile tracks for our Gaia XP sample are shown in grey for comparison. The tracks between Gaia XP and APOGEE samples show remarkable resemblance in high-α and small differences for low-α and all stars due to lower contamination in low-α selection in APOGEE (see Appendix A).

In the text

Fig. 5 Graphical model representation of the evolving means model, modelling the conditional distribution P(vϕ |[M/H]) with a two-component GMM, where the means and standard deviations are fixed (left) and varying (right) with respect to the underlying metallicities. The evolving or varying means model is also implemented for circularity as the conditional distribution or P(η |[M/H]). The subscript 1 and 2 stands for the disc-like and halo-like components, respectively, with μ and σ as the Gaussian mean and standard deviations in vϕ or η, conditioned on the value of [M/H]. w stands for the relative contribution of the disc-like component with the relative contribution of the halo-like component defined as (1−w). A list of the model parameters and their priors and functional forms are given in Table 1.

In the text

Fig. 6 Fractional contribution of disc-like (purple) and halo-like (orange) GMM components as a function of metallicity for high- and low-α stars. The scatter points show the location of the knots chosen to run our model. The bands are 1 σ uncertainties on the converged weights. These fractions are computed for a fixed Gaussian for halo and disc-like stars when the MCMC is converged. The converged azimuthal velocity and velocity dispersion is shown in the legend, shown in units of km s−1. The low-α model clearly shows two separate populations (thin disc and accreted) shown by the step function that could be misclassified as a rapid spin-up, while the high-α model shows a shallower monotonically increasing(decreasing) profile for the disc(halo) with [M/H].

In the text

Fig. 7 High- α model with evolving means and standard deviation for azimuthal velocity of the disc-like component. Left: Evolution of azimuthal velocity and azimuthal velocity dispersion as a function of metallicity, capturing a monotonic transition of a chaotic dispersion-dominated state (predominantly made of proto-Galactic populations) to a rotationally supported state with disc-like motion in purple; the velocity dispersion evolution of halo is shown in orange. Centre: Evolution of azimuthal velocity over azimuthal velocity dispersion, Vϕ/σϕ as a function of metallicity for the disc phase and halo component is shown in purple and orange, respectively. Right: Fractional contribution of disc phase (proto-Galaxy to high-α disc) (purple) and halo-like (orange) GMM components as a function of metallicity. In all the panels, the scatter points show the location of the knots chosen to run our model. The bands are 1 σ errors on the converged weights, velocities, and velocity dispersions. The fractions are computed for an evolving Gaussian for the disc phase and fixed Gaussian for the halo phase when the MCMC is converged.

In the text

Fig. 8 Distribution of sampled data column normalised by sum (left), model colour-coded by the probability density function (centre), and the normalised residual (right) in the [M/H]−vϕ plane (azimuthal velocity vs metallicity) for high-α evolving means model. The flat distribution of red stars at higher metallicities in the residuals show that the model is not able to fit the non-Gaussianity (long tail towards lower velocities) of the high-α disc. The systematic scatter in the periphery is due to random error in the data. However, for [M/H]<−0.7, the lack of strong systematic patterns, the low-amplitude, and the small scatter of the residuals in the central regions (with the major bulk of the data) verifies the validity of the model. The median, 16th and 84th percentile tracks for our data are shown in black lines to guide the eye towards the bulk of the data.

In the text

Fig. 9 Top: column-normalised (by sum) 2D histogram of stars in the [M/H]−η plane (orbital circularity vs metallicity, top) and [M/H]−vϕ plane (azimuthal velocity vs metallicity, middle bottom) for all the stars (right), high-α selection (left), and low-α selection (midde). The running median track is shown as dashed black line and the 16th and 84th percentile tracks are shown as black lines in all panels. The running median tracks for the all stars and low-α panels look more like a step-function, while the high-α tracks are shallower, supporting the more gradual spin-up phase with respect to metallicities. (middle top) Probability density functions (PDF) of η in bins of [M/H] for all the stars (right), high-α selection (left), and low-α selection (middle). (bottom) Probability density functions (PDF) of vϕ in bins of [M/H] for all the stars (right), high-α selection (left), and low-α selection (middle). We note the second peak with low net spin in high-α panel slowly gains rotation with increasing [M/H] and spins-up into the old high-α disc in both middle top and bottom panels.

In the text

Fig. 10 High- α model with evolving means and standard deviation for orbital circularity of the disc-like component. Left: Evolution of orbital circularity, and orbital circularity dispersion as a function of metallicity, capturing a gradual transition of a chaotic state (predominantly made of proto-Galactic populations) to a rotationally supported disc-like motion in purple and the orbital circularity dispersion evolution of halo is shown in orange. Right: Fractional contribution of disc phase (proto-Galaxy to high-α disc) (purple) and halo-like (orange) GMM components as a function of metallicity. In all the panels, the scatter points show the location of the knots chosen to run our model. The bands are 1 σ uncertainties on the converged weights, circularity, and circularity dispersions. The fractions are computed for an evolving Gaussian for the disc phase and fixed Gaussian for the halo phase when the MCMC is converged.

In the text

Fig. 11 Vϕ/σz ratio (azimuthal velocity to vertical velocity dispersion) for the most prograde Gaussian component in each metallicity bin (spaced evenly on a log scale) from our GMM runs in Section 6.3. The square size indicates the relative weight of the most prograde Gaussian in each metallicity bin. This ratio is a unitless quantity indicating the level of rotational support within a disc, essentially quantifying how dynamical cold it is at z=0. The dashed line at Vϕ/σz=1 marks the boundary between pressure-supported (halo-like) vs rotation-dominated kinematics (disc-like).

In the text

Fig. A.1 Column-normalised (by sum) 2D histogram of stars in the [M/H]−vϕ plane for all the stars (right), high-α selection (left), and low-α selection (middle) towards the Galactic centre (top panels) and away from the Galactic centre (anti-centre, bottom panels). The running median track is shown as dashed black line and the 16th and 84th percentile tracks are shown as black lines in all panels. We can see that the low-α stars are not fully two separate populations towards the inner Galaxy due to contamination from the proto-Galaxy in the low-α end, that creates a small connection between the thin disc and halo populations in the overall [M/H]−vϕ plane for low-α stars as seen in Figure 3. Low- α stars towards the Galactic anti-centre show a cleaner, two separate population of thin disc and halo stars as expected. This effect is also reflected in the all stars panel. This effect is minimal in the APOGEE α-selection, most likely due to the high-resolution α-measurements.

In the text

Fig. A.2 Logarithmic Galactic map of our sample of stars with high-α (left) and low-α (right) selections for the all-sky Gaia XP sample (top) and APOGEE DR17 sample (bottom), in the metallicity regime where the gradual spin-up in high-α vs the step function due to two different stellar populations in higher and lower metallicities in low-α is reflected the most: −1.3<[M/H]<−0.7. In the Gaia XP sample, we can see that the poor old heart or the proto-Galaxy is still present (contamination) in low-α towards the inner Galaxy, however, to a smaller extent than in high-α. This is not the case for APOGEE stars which have a much cleaner low-α population towards the inner Galaxy with the bulk of proto-Galaxy towards the inner Galaxy in the high-α Galactic map.

In the text

Fig. B.1 Orbital circularity vs metallicity for high-α, low-α, and all stars, and 1D histogram of azimuthal velocity and orbital circularity at different metallicity bins for column-normalisation by amplitude (top panels) and by sum (bottom panels) using Chandra et al. (2024) α-selection (left panels) vs the α-selection implemented in this work (right panels). Running median (black line), mean (black dashed line), mode (grey dashed line), 16th and 84th percentile (black lines) tracks are overlaid for the 2D histograms with column normalisation by amplitude. We can see that the mean and median tracks follow each other closely except at higher metallicities where the non-Gaussianity of thin and thick discs dominates, whereas the mode tracks follow the peaks of the background and is very noisy and resembles the step function behaviour seen by Chandra et al. (2024), even for high-α stars where we see a gradual spin-up (across increasing [M/H]) using column normalisation by sum. Running median, 16th and 84th percentile tracks are overlaid as black lines for the 2D histograms with column normalisation by sum. We also see that our α-selection has a cleaner high-α sample and isolates the bulk of GES into the low-α sample.

In the text

Fig. B.2 Row-normalised conditional metallicity distribution for high-α stars across orbital circularities for our α-selection (top) and Chandra et al. (2024) α-selection (bottom). The high-α disc, in situ halo (hot high-α disc), and accreted components can be seen in both panels. However, the contamination from the accreted GES merger is stronger in the Chandra et al. (2024) α-selection than the α-selection proposed in this work.

In the text

	Fig. C.1 Mean velocities (top), velocity dispersions (middle) and Vϕ/σz ratio (bottom) showing how rotationally supported the stars are, for the most prograde Gaussian component in each metallicity bin from our GMM runs for high-α stars.
In the text