A slow spin to win - the gradual evolution of the proto-Galaxy to the old disc

The purple band in Figure 1 shows the selection cut used to separate high-/low-$\alpha$ populations. Stars in the purple band are cut out to ensure we have purer samples of high-/low-$\alpha$ stars. This is especially important for high-$\alpha$ stars since we use the high-$\alpha$ sample as the sample of stars tracing the evolution of the proto-Galaxy to high-$\alpha$ disc over metallicity and any contamination from the denser low-$\alpha$ disc can affect the purity of the high-$\alpha$ 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 $\alpha$-separation, the consequences of which are discussed in section 6.

The selection is defined using the following equations:

This selection is somewhat different from the one implemented by Chandra et al. (2023) who use the same sample and abundances from Gaia XP spectra. The main difference between Chandra et al. (2023) and our selection are: (i) we use a more stringent selection for high-$\alpha$ to make sure that the bulk of accreted Gaia-Enceladus-Sausage (GES) merger (Belokurov et al., 2018; Helmi et al., 2018) does not contaminate the high-$\alpha$ population (see Figure 1), and (ii) the purple band is also quantitatively larger to ensure a sample with highest purity. We test the validity of this separation by using the $\alpha$-selection on APOGEE DR17 data in the [M/H]-[$\alpha$/M] space and using [Al/Fe]-[Mg/Mn] space to see if they occupy accreted, low-$\alpha$ or high-$\alpha$ 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-$\alpha$ stars respectively. The high-$\alpha$ selection used here is stricter than the selection used in the literature, in order to remove as much accreted low-$\alpha$ 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.

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 using reliable astrometry, 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 center (GRAVITY Collaboration et al., 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 for every grouping of stars in each 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-/low-$\alpha$ samples, respectively. It is worth noting that the high-$\alpha$ stars probe a smaller range in R-z compared to low-$\alpha$ stars; most of the low-$\alpha$ stars at larger scale heights are expected to come from accretion events, and this is likely why the low-$\alpha$ sample covers a wider scale height. We see a strong and steep negative metallicity gradient with increasing z for the low-$\alpha$ sample; we also see this in the full sample, likely because it is dominated by stars in our low-$\alpha$ selection, while the high-$\alpha$ sample has a shallower negative metallicity gradient with increasing z. The high-$\alpha$ disc has a larger scale height compared to the low-$\alpha$ disc and we also see a clear disc flaring for the low-$\alpha$ disc (as seen for the Galactic thin disc) with a metallicity gradient towards larger R, reminiscent of radial migration (see for e.g., Haywood et al., 2013; Hayden et al., 2015; Ratcliffe et al., 2024).

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-$\alpha$ 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²²2In this work, 2D column-normalised histogram automatically means column normalisation by sum, such that the sum under the histogram curve equals 1. This is directly proportional to the probability density function of azimuthal velocity at each metallicity bin. for all, high-, and low-$\alpha$ 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 acreted stars), and low-$\alpha$ 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-$\alpha$ sequence have a relatively consistent and monotonic AMR. They also show that the high-$\alpha$ disc reached solar metallicities at around $\sim 8$ Gyr in lookback time. Therefore, the high-$\alpha$ panel should trace the first $\sim$ 5-6 Gyr of the Milky Way’s evolution, which is also proposed by Chandra et al. (2023) in their orbital circularity versus metallicity plane. All the 2D histogram in Figure 3 have 16^th, 50^th, and 84^th 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 will discuss the different panels separately in the following subsections.

Figure 4 shows 2D column-normalised histograms of [M/H]-v_ϕ plane for high-$\alpha$ (left), low-$\alpha$ (center), and all stars (right) for APOGEE DR17 stars with the same $\alpha$-selection described by equation 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 [$\alpha$/M]. The running 16^th, 50^th, and 84^th percentile tracks for the APOGEE data are shown as black lines and the tracks from our final sample in Figure 3 are shown in gray in all three panels. We find very good resemblance between the track in the APOGEE and Gaia XP data, especially for high-$\alpha$ stars stars. It shows that a scenario of a the gradual spin-up of proto-Galaxy to the high-$\alpha$ disc population is not due to large errors, as it is supported by a much more reliable high-resolution chemical abundance data. The low-$\alpha$ 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-$\alpha$ sample, while APOGEE has a much cleaner and purer low-$\alpha$ selection, due to higher quality of the chemical abundances. This can also be seen by the thin disc and halo populations being completely disjoint in the low-$\alpha$ 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-$\alpha$ 2D histograms. 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-$\alpha$ disc and kinematically colder low-$\alpha$ 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 will 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.

In order to quantify the spin-up of the high-$\alpha$ 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 kms^-1 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$-\alpha$ disc 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 kms^-1 and standard deviation of 50 kms^-1, with low/high limits restricted between 50 and 200 kms^-1.

It is important to note that the high-$\alpha$ 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-$\alpha$ 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-$\alpha$ disc. The metallicity range used for high-$\alpha$ 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 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⁵⁵5This is because Gaussian processes can be seen as an infinite-dimensional generalization of multivariate normal distributions. described below:

where $\sigma^{2}$ is the overall variance, $\ell$ is the characteristic length scale of the covariance function, that describes the range in [M/H] to which the covariance function typically holds, and $\alpha$ is the positive scale-mixture parameter of the covariance function ($\alpha$ ¿ 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 centered at 120 kms^-1 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 centered at 150 kms^-1. 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 number of samples to generate from the Markov Chain for the high-$\alpha$ 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-$\alpha$ 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-$\alpha$ sample.

Figure 7 shows the converged parameters as a function of metallicity. As in Fig 6, the metallicity knots are shown as scatter points and $1\leavevmode\nobreak\ \sigma$ 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-$\alpha$ disc-like component. Here, the trend gradually increases from v${}_{\phi}\sim 50$ kms^-1 at [M/H] $\sim-2$ to v${}_{\phi}\sim 200$ kms^-1 at [M/H] $\sim-0$. This result can be interpreted as the high-$\alpha$ 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., 2023; Zhang et al., 2023), this is the first time that high-$\alpha$ 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.5 dex (between –2 and –1.5 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] $\sim-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${}_{\phi}=0$ kms^-1) 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 kms^-1) compared to the rest of the halo resting at about 100 kms^-1 velocity dispersion (see Figure 3). The rise in standard deviation at higher metallicities (when the high-$\alpha$ 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-$\alpha$ disc. This is one of the disadvantages of using a Gaussian distribution.

The middle panel of Figure 7 shows the ratio of v_ϕ to $\sigma_{\phi}$, 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${}_{\phi}/\sigma_{\phi}$, reaching up to a ratio of 6 at solar metallicity. Moreover, the right panel of Figure 7 shows the relative contribution of halo-like (orange) and disc-like (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]$\lesssim-1.5$. In the very metal-poor end ([M/H]¡–2) of our high-$\alpha$ 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 modeling purpose. This component captures the proto-Galaxy to spin-up phase to high-$\alpha$ disc. The disc-like component’s fractional contribution increases slowly from [M/H]$\sim-2$, with an approximately constant gradient, up to [M/H]$\sim$-0.5. Upon reaching this point, the disc-like component (i.e., high-$\alpha$ disc) dominates the sample. All these results are highly in favour of a gradual spin-up for the disc-like component. 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-$\alpha$ disc.

In Figure 8, we show the [M/H]-v_ϕ plane for high-$\alpha$ 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 16^th, 50^th, and 84^th 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 kms^-1 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 less stars than what is present in reality. We conjecture this is because the v_ϕ distribution of high-$\alpha$ 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-$\alpha$ disc than the low-$\alpha$ 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-$\alpha$ 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-$\alpha$ 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-$\alpha$ 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-$\alpha$ 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 in v_ϕ across 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-$\alpha$ (or proto-Galaxy population) gradually spins-up to a rotation-supported high-$\alpha$ 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-$\alpha$ regime, capturing the evolution of the first 5-6 Gyr of the of proto-Milky Way populations to the high-alpha 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.8).

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 circulairity (as it abruptly cuts at -1 and 1 for a perfectly circular retrograde and prograde orbit 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-$\alpha$ stars because they trace the transition of the proto-Galaxy to rotation-supported high-$\alpha$ disc more clearly as seen in Figure 9, while the low-$\alpha$ 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 low/high limits restricted between 0.2 and 0.7 across the metallicity knots. We use the same metallicity range as in subsection 4.1, 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 covaraince matrix to enable the information of between each components across the metallicity is unchanged in this circularity model, with the standard deviation on the mean of the circularity has 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 centered at 0.3 and the standard deviation is described by a half normal prior centered 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] $\sim$ –1.7 to –1 and is still steadily increasing up to –0.5. However, the interpretation of metal-rich stars ([M/H]¿–1.0) is trickier given that the long asymmteric 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] $\sim$ –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, as opposed to the rapid spin up shown in the literature (see for e.g., Chandra et al., 2023; 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] $\sim$ –1.2) compared to velocity ([M/H] $\sim$ –1.6), and (ii) the mean orbital circularity is more circular ($\eta$ $\sim$ 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_ϕ $\sim$ 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-$\alpha$ sample. The latter difference could also arise from the fact that in the literature, the metal-poor high-$\alpha$ stars are not modeled 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-$\alpha$ 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-$\alpha$ disc is more gradual and not as rapid at [M/H] $\sim$ –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.

In this work, we have set out to model the azimuthal velocity and orbital circularity evolution of high-/low-$\alpha$ 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-$\alpha$ disc. At first, we model the conditional distribution P(v${}_{\phi}\mid$[M/H]) using a 2-component (disc-like and halo-like) Gaussian Mixture Model with frozen means and standard deviations for the azimuthal velocities for the high- and low-$\alpha$ stars. We find the inferred posterior matches better with the data for low-$\alpha$ stars better than high-$\alpha$ stars. The fractional contribution from disc-like and halo-like components look like a step function at [M/H]$\sim$–1 for low-$\alpha$ stars, while the high-$\alpha$ stars have a relatively gradual rise over increasing metallicities in the fractional contribution (Figure 6). Secondly, given that the high-$\alpha$ stars are not modeled well by the frozen means model, we model the conditional distribution P(v${}_{\phi}\mid$[M/H]) using a 2-component (disc-like and halo-like) Gaussian Mixture Model with evolving means and standard deviations for the azimuthal velocities for the high-$\alpha$ 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_ϕ$\sim$50 km/s over increasing metallicities (–1.7¡[M/H]¡–1.0). Thirdly, we perform the same exercise for high-$\alpha$ stars in orbital circularity versus metallicity plane, [M/H]-$\eta$, given that the orbital circularity is position independent, as opposed to azimuthal velocity. From this exercise, we also see a gradual rise in oribital circularity starting from an $\eta\sim$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-$\alpha$ 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-$\alpha$ disc (likely the proto-Galaxy) to a rotationally-supported high-$\alpha$ disc at higher [M/H]. On the contrary, due to the superposition of the GES debris and the low-$\alpha$ disc in the low-$\alpha$ 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]$\sim-1$. Due to the GES debris dominating the metal-poor sample for all stars in our data set, this yields a similar profile when inspecting the Gaia XP sample without any $\alpha$ selection. Thus, our results highlight the importance of the [$\alpha$/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, eventually settling into a high-$\alpha$ disc and the disc formation was not as rapid or dramatic across a narrow range of metallicities as previously reported.

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 w.r.t. the Galactic disc (v${}_{\phi}\sim 50$ km/s), should primarily be concentrated towards the innermost Galactic regions ($r\lesssim 10$ kpc), and should host stars with high [$\alpha$/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. 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 data set of [M/H] (Andrae et al., 2023a) and now also [$\alpha$/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-$\alpha$ disc (e.g., Chandra et al., 2023; Zhang et al., 2023). However, most of these studies either: ($i$) examine the profile of azimuthal velocity, v_ϕ, and/or orbital circularity, $\eta$, with [M/H] agglomerating stars with different [$\alpha$/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 rotatinally dominated populations in the velocities-[M/H] plane using disconnected Gaussian Mixture Models (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.

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 v_ϕ-metallicity plane; these authors found that the Milky Way “spun up” rapidly between $-1.3<\mathrm{[Fe/H]}<-0.9$ (see also Zhang et al. (2023) for a similar result using Gaia XP metallicities). Along those lines, Chandra et al. (2023) used the available [$\alpha$/Fe] for the Gaia XP sample to study the profile of [M/H]-$\eta$, finding that the high-$\alpha$ disc emerges rapidly at $\mathrm{[M/H]}\sim-1$.

In this work, we have found that when performing a similar [$\alpha$/M] cut to Chandra et al. (2023), our transition from proto-Milky Way populations to the high-$\alpha$ disc (Figure 9) shows a much more gradual increase when compared to these previous works. In more detail, both when examining and modelling the v_ϕ-[M/H] and [M/H]-$\eta$ 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<\mathrm{[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 Belokurov & Kravtsov (2022) already shows preliminary hints of a gradual spin up (see their Figure 5); ($ii$) Zhang et al. (2023) perform a traditional GMM without any $\alpha$-selection, resulting in the GES debris driving the rapid spin up inference at a narrow range in metallicities; ($iii$) (Chandra et al., 2023) shows a qualitative view of [M/H]-$\eta$ plane, while their high-$\alpha$ 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-$\alpha$ disc transition has likely been gradual over increasing metallicities and not as dramatic over a narrow range of metallicities as previously reported.

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]$\approx-0.8$. We have tested to ensure 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-$\alpha$ 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\lesssim r\lesssim 30$ kpc: Belokurov et al., 2018; Helmi et al., 2018; Koppelman et al., 2018, e.g.,). This merger possibly catalyzed the formation of the in-situ halo, or ”hot high-$\alpha$ disc” comprising stars on halo-like orbits with chemistry similar to the high-$\alpha$ 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 Milky Way disc.

In this subsection, we perform the traditional $N$-component Gaussian Mixture Model (GMM) decomposition of high-$\alpha$, low-$\alpha$, 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.

Typically, within any GMM framework, each sub-population is characterized 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. (2023) for the GMM decomposition, to be able to compare the effects of $\alpha$-separation in this type of GMM decomposition. To perform the GMM analysis, we utilized 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-$\alpha$ and high-$\alpha$ 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-)$\alpha$ samples. The all stars sample is fit using the same six bins as in the low-$\alpha$ sample⁶⁶6For the all stars and low-$\alpha$ stars samples, we neglect the most metal-rich bin as it is affected by the disc’s asymmetry and suggests up to 8 components as the optimal fit, which is un-physical.. Zhang et al. (2023) 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. (2023). 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 $\sim 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-$\alpha$ and 2 bins in low-$\alpha$ and all stars sub-samples. We do not compute the uncertainty in the GMM parameters as Zhang et al. (2023) found them to be generally around $\sim 0.1$ kms^-1. 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_ϕ/$\sigma_{{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 is capturing well the halo and disc components, we verify that the other GMM components the model finds fit for a halo-like substructure (V_ϕ/$\sigma_{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-$\alpha$ stars gradually increase in their V_ϕ/$\sigma_{{z}}$ (in purple) towards higher metallicity bins while low-$\alpha$ (in orange) and all stars (in black) have a more rapid increase in V_ϕ/$\sigma_{z}$ between two [M/H] bins, $\sim-1.3$ and $\sim-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 $\alpha$-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-$\alpha$ populations, we find that the high-$\alpha$ prograde GMMs favour a gradual spin-up of a metal-poor high-$\alpha$ population (likely the proto-Galaxy) to a metal-rich high-$\alpha$ disc. The GMMs from Zhang et al. (2023) 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 $\alpha$ estimate from Li et al. (2024), whereas Zhang et al. (2023) only used the metallicity estimates and therefore, we end up having much less VMP stars than theirs. Zhang et al. (2023)’s results also suggest a more rapid growth in V$\phi$/$\sigma_{\phi}$ 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-$\alpha$ 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-$\alpha$ stars tracing the evolution of the proto-galaxy is shown in Appendix C.

In summary, our results reveal that the high-$\alpha$ subsample shows a gradual rise in V_ϕ/$\sigma_{{z}}$, favouring a gradual spin-up phase for the Milky Way’s high-$\alpha$ disc.

In this work, we use an $\alpha$-separation on Gaia XP based [$\alpha$/M] and [M/H] abundances to separate high-/low-$\alpha$ stars and use the high-$\alpha$ stars to trace the azimuthal velocity (v_ϕ) evolution of the old Milky Way disc over a range of metallicities ($-2.5<\rm[M/H]<0.1$). The high-$\alpha$ 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 mergers other populations (e.g., Heracles and the high-alpha tail of the GES and possibly other merger remnants) overlap, so a simple $\alpha$-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 section 4.1 and 5.1. However, this model assumes a Gaussian distribution of azimuthal velocity at any given metallicity, which is not strictly true for the high-$\alpha$ 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-$\alpha$ 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-$\alpha$ 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-$\alpha$ disc state. 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-$\alpha$ disc state as a function of increasing metallicity. We reserve this exploration for 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 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 [$\alpha$/Fe] with our high-$\alpha$ sample, given our $\alpha$ separation. Recent studies such as Belokurov & Kravtsov (2022); Conroy et al. (2022) used APOGEE and H3 surveys, respectively, to show that the [$\alpha$/Fe] ratio gradually declines between $-3>$ [Fe/H] $>-1.3$, and then rises to a higher value to meet the high-$\alpha$ disc population. It is very likely that we are not catching this dip in [$\alpha$/Fe] with our simple $\alpha$-separation. However, this should not affect the conclusions of this work.