The formation and survival of the Milky Way’s oldest stellar disk

The formation of our Galaxy’s disk proceeded in two major phases that are almost disjoint in time^1,2. The earlier phase, which occurred in the first ~5 Gyr, resulted in the formation of a rapidly enriched (or high-α) old disk. Its stars are separated from those of the younger disk formed in the later phase by chemistry and kinematics^3,4,5,6,7, although the detailed picture of transition and interplay between these disks is still under debate^8,9,10,11,12.

Recent studies with precise stellar age, kinematics and abundances suggest that the onset of the Galactic old stellar disk was in the first billion years of the universe (age τ ≈ 13 Gyr, or redshift z ≈ 7), when gases were in a relatively metal-poor ([Fe/H] ≲ −1), chaotic state^2,13,14. The early formation of this old disk apparently preceded the epoch at which much of the stellar halo population was established by means of accretion of satellite galaxies². Early galactic disks have also been revealed in the distant Universe at z > 7 by recent James Webb Space Telescope (JWST) observations¹⁵. On the other hand, galaxy formation scenarios and simulations suggest that galaxy assembly in the earliest epoch was dynamically violent as merger events happened frequently^16,17, which was detrimental to the creation and maintenance of galactic stellar disks. It thus remains a mystery when our Milky Way first formed stars as disks, and what the oldest disk component is that did not get destroyed by subsequent mergers.

Although stellar chemistry and kinematics provide central information on the Milky Way’s old disk formation, the most intuitive way to understand the disk’s temporal evolution is to characterize how its structure depends on the age of the stellar population^18,19,20. For this, precise stellar ages for a large sample of stars with a well-defined selection function are needed. Such stellar samples have become available only in recent years thanks to large surveys. Such attempts have been made using spectroscopic data sets from the Apache Point Observatory Galactic Evolution Experiment (APOGEE)¹⁹ or from the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST) survey²⁰. However, due to relatively large uncertainties (20–40%) in the spectroscopic age estimates, it is difficult for those studies to unravel details about the earliest formation history of the high-α, old disk. Fortunately, stellar ages can now be estimated with much higher precision by combining the spectroscopic stellar parameters with precise astrometric data from the Gaia mission^21,22,23. In particular, in ref. ², ages were derived, with a median precision of ~8%, for a quarter of a million subgiant stars using Gaia astrometric parallax and photometry, as well as spectroscopic parameters from the LAMOST Galactic surveys^24,25,26.

Here we present a reconstruction of the stellar mass density profiles for different mono-age and mono-abundance subpopulations of the Galactic old, high-α disk, utilizing an updated version of the subgiant star sample of ref. ². The spatial distribution of the updated subgiant sample in the R–Z plane of the Galactic coordinate is presented in Extended Data Fig. 1. The selection process of the high-α, rotating disk stellar populations used in the current work is illustrated in Extended Data Figs. 2 and 3. We describe the stellar mass density distribution as \(\rho (\mathbf{r}| \tau ,[{\rm{Fe}}/{\rm{H}}];\{\,{\rho }_{0},{H}_{R},{H}_{Z}\})\), where the model parameter ρ₀ refers to the disk mass density at the solar position, H_R is the disk scale length and H_Z is the scale height. These parameters are determined through forward modelling that rigorously incorporates the selection function of the sample, and global parameter fitting is carried out using the Markov Chain Monte Carlo (MCMC) method (see Methods for details). An example of the model fits to the observed distance distributions for a few mono-age and mono-abundance subpopulations are shown in Extended Data Fig. 4.

The best-fit parameters for each mono-age and mono-abundance populations are shown in Fig. 1. Figure 1a reveals a high-α disk sequence with high local density along a tight age–[Fe/H] relationship from [Fe/H] ≈ −1 at 13 Gyr ago to [Fe/H] ≈ 0.5 at 8 Gyr ago. The metal-poor tail of this sequence, in the range −2.5 ≲ [Fe/H] ≲ −1.0, exhibits dramatically lower local stellar density. A more intuitive presentation of the local mass density trend is plotted in Fig. 1d, which shows a dramatic jump of ρ₀ at [Fe/H] ≈ −1 at 13 Gyr ago, as ρ₀ increases by about an order of magnitude from [Fe/H] ≲ −1.2 to [Fe/H] = −1.0. This is an independent verification of the previous finding through the [Fe/H]–[α/Fe] morphology that the star formation efficiency experienced a sudden increase at the early epoch, for [Fe/H] ≈ −1.3 (ref. ¹⁴).

a–c, The distribution of best-fit parameters in the age–[Fe/H] plane. Panels show results for the local stellar mass density (a), scale length (b) and scale height (c). The dashed lines delineate the parameter window in age–[Fe/H] where the high-α disk is expected to dominate. We retained only these subsamples for the subsequent analysis, as contamination by the low-α stellar populations (due to measurement errors) may dominate beyond. d–f, The best-fit parameters as a function of age for stellar populations in the selected age windows of the upper panels, for local stellar mass density (d), scale length (e) and scale height (f). The error bars shown in the figure represent the typical range of parameter estimates (that is, 84th–16th percentile of the integrated likelihood function) for several populations of different metallicity.

Source data

Full size image

Figure 1b,e suggests that the scale length of this old, high-α sequence is 1–3 kpc, which evolves only marginally in the full range of metallicity and age, except for the extremely metal-rich ([Fe/H] ≈ 0.4, τ < 10 Gyr) populations. However, the scale height decreases by an order of magnitude, from 2 kpc for [Fe/H] < −1 to 0.2 kpc for [Fe/H] ≳ 0 (Fig. 1c,e). For a constant [Fe/H], the scale height changes only moderately with age. These results suggest that the high-α disk has experienced a much more dramatic structure evolution in the vertical direction than in the radial direction.

We note that Fig. 1 also displays notable populations of relatively young (≲8 Gyr), metal-poor ([Fe/H] ≲ −0.5) stars for which our model fitting results in a negative scale height (Fig. 1b). These populations are likely to be contamination of the low-α disk as a consequence of imperfect population distinction due to measurement errors in the [α/Fe]. Therefore, in this work, we confined our discussion about the high-α disk for populations in an age window along the distinguished, tight age–[Fe/H] sequence, as marked by the solid lines in the upper panels of Fig. 1. Also note that here we have ignored the disk flaring effect in the modelling. Similar results are found after taking the flaring effect into consideration (see Extended Data Fig. 5 and more detailed discussions in the Methods and the Supplementary Information).

Figure 2 compares the scale length and scale height of the high-α disk for the mono-age and mono-abundance subpopulations. It illustrates that the scale height H_Z correlates tightly with age. The vast majority of subpopulations of ages τ ≲ 13 Gyr, or metallicity [Fe/H] ≳ −1, exhibit a disk-like geometry, with H_Z smaller than H_R. Such a disk-like geometry even presents for older populations (τ ≈ 13.5 Gyr), although the majority of the oldest populations (τ ≳ 14 Gyr) have a scale height comparable to the scale length. The scale height drops below half the scale length for subpopulations with τ ≲12.5 Gyr.

The left and right panel are colour-coded by metallicity and age of the stellar subpopulations, respectively. The open circles without error bars in the figure refer to results where the MCMC fit failed to yield a valid parameter constraint. The scale height is smaller than the scale length for nearly all components with τ ≲ 13 Gyr, and also for some subpopulations with even older ages. The shaded regions mark the parameter space where the scale height is smaller than half the scale length. Some of the most metal-poor ([Fe/H] ≲ −1.5) subpopulations exhibit smaller scale height than scale length, but many of them have ages younger than 13 Gyr.

Source data

Full size image

However, the metallicity dependence for populations of [Fe/H] ≲ −1 is less clear. Some of these metal-poor subpopulations exhibit scale heights comparable to the scale length, but a notable number of them, even for those of [Fe/H] ≲ −1.5, can exhibit disk-like geometry. On the other hand, whereas the majority of the subpopulations with [Fe/H] ≈ −1 exhibit disk-like geometry, some of them may have a scale height comparable to the scale length. This complex, non-monotonic relationship between the geometry and metallicity might indicate a complex early assembly and enrichment history of our Galaxy.

Interestingly, many of the most metal-poor subpopulations (\(\left[{\rm{Fe}}/{\rm{H}}\right]\lesssim -1.5\)) that exhibit disk-like geometry tend to have younger ages (τ < 13 Gyr) than the populations of similar metallicity but with H_Z ≈ H_R. This implies that these metal-poor but younger disk-like populations may have experienced a different origin to the older populations of similar metallicity. It is possible that they are a part of accreted populations with coplanar orbits.

In Fig. 3, we further analyse the age dependence of the high-α disk structure, looking at the scale height-to-length ratio H_Z/H_R as a function of stellar age. We compare the observational results with analogous quantities for 31 Milky Way-like galaxies from the TGN50-1 simulation^27,28. These simulated galaxies were selected from the sample of ref. ²⁹, which contains 198 analogues of M31 or the Milky Way. In each of them, we chose in situ stars with angular momentum L_z > 500 kpc km s⁻¹ and galactic distances from 6 to 12 kpc. When calculating the present H_Z/H_R in each age bin, we incorporated a normal age error for all simulated stars. Here we only selected the 31 Milky Way analogues (31% among a total of 99) whose present-day H_Z/H_R-age trend resembles our Milky Way: that is, had a smoothly decreasing H_Z/H_R with age (solid grey line). For these 31 simulations, we compared the present-day structure with the same stars’ structure at birth (blue dashed line). Figure 3 illustrates that the simulated stellar populations can be born disk-like at the beginning, within the first gigayear of the Universe, but they are considerably thicker now than at birth due to dynamical heating. These results imply that even the oldest stars in our observation, which have H_Z/H_R ≈ 1, might have initially formed in a disk-like structure but were subsequently destroyed by frequent galaxy mergers at the early epochs. They might now form the high-angular-momentum tail of the “poor old heart of the Milky Way”^13,30.

The dots are our measurements, colour-coded by [Fe/H] of the stellar populations. Only populations with relative uncertainties in H_z/H_R < 50% are shown. The solid red line is the mean observed H_z/H_R(τ) relationship marginalized over abundances. The solid and dashed lines in black show the average height-to-length ratio for Milky Way analogues in the TNG50 simulation, with the shaded regions indicating the variance among individual galaxies. The solid black line shows the stellar distribution morphology at present, after adding an age error of 8% to the simulation output to mimic the observation. The dashed line shows the morphology at birth of the corresponding population. Populations with formal ages older than the Universe result from the age uncertainties.

Source data

Full size image

Nevertheless, there is a slight difference in the present-day height-to-length ratio between the simulation and our measurements. In the simulations, most components older than ~10.5 Gyr have H_Z/H_R ≳ 0.5, whereas in the Milky Way this is only the case for stars of τ > 12.5 Gyr. The Milky Way’s H_Z/H_R is below 1σ border of the simulated galaxies, which means that only ≲1/6 of the simulated Milky Way-like galaxies have surviving disk populations with H_Z/H_R comparable to the Milky Way. The reason for such a difference remains to be investigated. It may point towards a problem in the current simulations, but it may also be an indication that our Milky Way is extremely quiet so that the disk heating is not efficient compared with its analogues. We know that our Milky Way has experienced at least one prominent merger beyond the first 2 Gyr, the Gaia–Encelaus–Sausage^31,32,33, whose final phase was completed less than 12.5 Gyr ago. This suggests that the impact of this early merger to our Milky Way’s disk has been less traumatic than seen in the simulations.

Our modelling also enables a rigorous estimate of the total stellar birth mass of each mono-age and mono-abundance disk component by integrating the density profile over all radii and heights. This yields the star formation history of the high-α disk as a whole. Figure 4 shows this as the star formation rate. It started out at 2 M_⊙ per year in the first few hundred million years (τ > 13 Gyr). This earliest high-α disk has built a total stellar mass of 3.7 × 10⁹ M_⊙, as of 13 Gyr ago. Of this initial mass, 2.2 × 10⁹ M_⊙ remain as of now, after considering the mass loss by stellar explosion, which reduces 40% of stellar mass according to the stellar evolution models. Most of these stellar masses, about 2.0 × 10⁹ M_⊙, survived as a disk that has H_Z/H_R < 0.75. The star formation rate increased to 6 M_⊙ per year 12.5 Gyr ago, and further to 11 M_⊙ per year 11 Gyr ago. Such a peak star formation rate is consistent very well with both the chemical modelling³⁴ and the TNG50 hydrodynamic simulations, in which the Milky Way-like galaxies exhibit a peak star formation rate of ~10 M_⊙ per year at the early epoch²⁹. However, we find that the star formation rate then dropped quickly about 10 Gyr ago, presumably because the gas was mostly depleted as a consequence of the high star formation rate before that epoch (however, see ref. ³⁵). Ultimately, the high-α disk assembled a total of 2 × 10¹⁰ M_⊙ stellar mass (including remnants) that remains at present.

The star formation rate (SFR) is derived by integrating the structural parameters in the spatial range of 0 < R < 20 kpc, −10 < Z < 10 kpc. At a given age, individual metallicity bins are co-added to generate the star formation rate shown in the figure. The error bars are estimated by propagating the fitting errors of the scale parameters. The cumulative stellar mass refers to the total stellar mass that the high-α disk formed before a given epoch. In fact, about 40% of this mass has been lost by means of stellar explosion, which turns the stars into gas and dust, providing the materials of chemical enrichment. Thus, one needs to multiply the number on the right-hand vertical axis by ~0.6 to obtain the present-day disk stellar mass.

Source data

Full size image

The onset of the Galactic disk has been suggested to happen at an early epoch with [Fe/H] ≈ −1.3, which is inferred either by modelling stellar distribution in the [Fe/H]–[Mg/Fe] plane¹⁴, or by characterizing the rise of stellar tangential velocities¹³. Combining precise isochrone ages with our new geometry measurements shows that the (still surviving) Milky Way’s high-α disk formation dates back to more than 13.5 Gyr ago. For many of the very oldest subpopulations in the sample (τ ≈ 14 Gyr), we do not find a clear disk-like morphology at present. It is, however, possible that these very oldest populations were once born as a disk but were destroyed by early mergers, as also suggested by simulations. Therefore, we suggest that our Galaxy may have started as a disk. We dub the earliest disk formed in the first few hundred million years (τ > 13 Gyr) that has survived from early mergers and remains disk-like today as PanGu(盘古), the God who created heaven in Chinese mythology. Coincidentally, the two Chinese characters of the name are translated literally as old disk, as Pan(盘) means disk and Gu(古) means old. The Galactic archaeology results resonate with the substantive incidence of disk galaxies seen in z > 6 JWST observations¹⁵.

Our results further suggest that PanGu has a present-day stellar mass of M_* ≈ 2 × 10⁹ M_⊙. Even more stellar mass may have been born as a disk along with PanGu but was destroyed by subsequent mergers. These destroyed disk populations may have partially contributed to the Galactic metal-poor ([M/H] < −1.5) heart³⁰ (see also ref. ¹³, for Aurora). In ref. ³⁰, a mass of M_* ≳ 10⁸M_⊙ was revealed for the (metal) “poor old heart of the Milky Way”. This is smaller than the mass of the PanGu, which implies that PanGu was likely to have been a dominant component of our Galaxy at the earliest epoch.

Our results also confirm that the vast majority of stars in the high-α disk formed at an extended burst that lasted from ~12 to 10 Gyr ago, which has been well established in previous work^2,8,36,37. The present work now can add an actual rate estimate at the peak, of 11 M_⊙ per year, and also a decent estimate of the total stellar mass assembled in the high-α disk, which is 2 × 10¹⁰M_⊙ at present. This estimate of star formation rate, by direct age-dependent star counting, is remarkably consistent with chemical modelling³⁴ and hydrodynamical simulations^17,29.

Our results have revealed that the disk scale height decreased by nearly an order of magnitude with decreasing age over the first 5 Gyr. At the same time, the variation of scale length was only moderate. This is likely to be a consequence of two effects: early on stars were born in a thicker disk (upside-down formation; refs. ^38,39,40,41), or they may have been heated subsequently (for example, see refs. ^38,39,40,42). Both effects are seen in the simulations. Nevertheless, a final scenario for the assembly of the high-α, old disk formation should explain, in addition to the scale height-to-length ratio evolution, other observations, for example, the negative vertical metallicity gradients^43,44,45. Future work could carry out a comprehensive comparison between simulations or scenarios and observations.

As mentioned above, our results also hint that the assembly and chemical enrichment history in the first billion years are complex, rather than the simple, closed-box gas enrichment that results in a single, monolithic age-metallicity relationship. The most metal-poor stellar populations can exhibit disk-like geometry, which confirms, from a direct view, the existence of a metal-weak (poor) thick disk^{46,47,48,49,50}. The results further show evidence that some of these metal-poor disk populations are probably results of early accretion events, as they tend to be younger than in situ stellar populations of similar metallicity.

Finally, a comparison of age evolution of the disk morphology with the TNG50 simulations suggests that only ≲1/6 of the simulated Milky Way-like galaxies have surviving disk populations with H_Z/H_R comparable with the Milky Way, which indicates that our Milky Way has experienced a very dynamical quiescent history. The impact of early mergers, such as the Gaia–Encelaus–Sausage^31,32, with our Milky Way’s disk may have been less traumatic than seen in the simulations.

An updated LAMOST–Gaia subgiant star sample

We built a sample of subgiant stars from the seventh data release (DR7) LAMOST Galactic survey (https://dr7.lamost.org)⁵¹ and third data release of the Gaia mission²³, following ref. ². Subgiants are stars in a brief evolutionary phase where their core contracts and releases gravitational potential energy after the end of core hydrogen burning. Stars in this phase allow the most precise age determination, as their luminosity is a direct and sensitive indicator of their ages⁵².

In ref. ², 247,104 subgiant stars were identified from the T_eff–M_K diagram utilizing the LAMOST spectroscopic and Gaia DR3 astrometric data sets, where T_eff refers to effective temperature and M_K refers to absolute magnitude in the infrared K band. These data afforded a median age uncertainty of only ~8%. Getting such a clean sample involved a number of steps: first, they discarded possible binary stars with a Gaia renormalized unit weight error (RUWE)^53,54 value larger than 1.2. They also further identified binary stars based on the excess of their parallax-based ‘geometric’ M_K with respect to their spectroscopic M_K, and discarded them from the sample. The ‘geometric’ M_K refers to K band absolute magnitude inferred from distance modulus using Gaia parallax, and indicates the total luminosity of the binary system, whereas the spectroscopic M_K refers to that deduced from the normalized LAMOST spectra with a data-driven approach, which provides a good luminosity indicator of the primary star, in particular, for equal-mass binary systems for which the spectral lines of the individual components are indistinguishable^2,55. This method effectively eliminates binaries with high mass ratio (≳0.7), which may lead to erroneous age estimates if treating as single stars. This is different from, and thus complementary to, the Gaia RUWE method, which may recognize binaries most effectively for a lower mass ratio. Reference ² also discarded stars with M_K < 0.5 from the sample to eliminate potential contamination of core helium burning blue horizontal branch (BHB) stars. Although these criteria leave a clean sample, they also result in incompleteness for stellar mass density reconstruction that needs to be characterized.

To consider the effect of α element enhancement on the age estimation, ref. ² derived the stellar ages with the Yonsei–Yale (YY) isochrones^56,57,58 of different [α/Fe], namely, [α/Fe] = 0, 0.2, 0.4. To compute the best age, they combined the age estimates from these isochrones, weighted by Gaussian whose mean and σ were taken from the LAMOST [α/Fe] measurements. It was subsequently realized that this strategy was imperfect because the typical measurement error in the LAMOST [α/Fe] (<0.05 dex; ref. ⁵⁹) is much smaller than the step of the isochrone grid (0.2 dex).

Therefore, in the current work, we update the Gaia and LAMOST subgiant star sample of ref. ², by inheriting their overall approach, but with several updates. First, we included binary stars in our sample by discarding both the Gaia RUWE and M_K criteria for binary elimination of ref. ², but with a careful treatment of their age and distance estimates. For single stars, both geometric M_K and spectroscopic M_K were adopted for age estimation. For binary stars, as their geometric M_K are overluminous due to the contribution from the secondary, only the spectroscopic M_K were used for a proper age estimation for the primary stars (thus the binary system). On the contrary, for distance estimates, it turns out that the geometric distance is a better choice for binaries, but the spectro-photometric distance can be wrong. So we used only the geometric distance for binaries. The inclusion of binaries added 70,228 stars to the sample in ref. ², among which 41,747 were recovered by removing the M_K criterion, and the others were recovered by removing the Gaia RUWE criterion.

Second, we included subgiant stars with M_K < 0.5 in our sample but eliminated BHB star contamination by setting a metallicity cut of [Fe/H] > −0.8. BHB stars are metal-poor halo populations and can be effectively eliminated by the metallicity cut. Meanwhile, subgiant stars within this magnitude range are expected to be young and metal rich, so that they can well survive the metallicity cut. This update added 6,076 stars to the sample. However, we note that these updates did not make a significant impact on the present work as we focused only on the old disk populations (τ ≳ 8 Gyr).

Third, we computed three ages for each star using the YY isochrones of [α/Fe] = 0, 0.2, 0.4, respectively, and then delivered the final age estimate by linear interpolation to fit the measured [α/Fe]. This treatment slightly improved the age of stars with [α/Fe] ≈ 0.1 or [α/Fe] ≈ 0.3, the middle intervals of the isochrone grids, although, on the whole, the impact on the age estimates was minor.

With these updates, we obtained a sample of 320,028 subgiant stars that had precise age, abundances and orbital parameters⁶⁰. The spatial distribution of these stars in the R–Z plane can be found in Extended Data Fig. 1.

Selection of high-α stars

We selected high-α disk stars through their abundances and kinematics, using the sa.me criteria as in ref. ² (see also Extended Data Fig. 2),

These criteria led to 86,001 high-α sample stars. In the subsequent analysis, we further discarded stars younger than 5 Gyr, as they may be dominated by low-α young disk stars that are misclassified due to uncertainties in abundance determinations. Also, a smaller fraction of them might be blue stragglers of the high-α disk population whose ages are incorrectly inferred by assuming single-star evolution isochrones⁶¹.

It was found that the high-α stellar populations for [Fe/H] > −1.2 were dominated by a rotating component with a peak at L_Z ≈ 1,200 kpc km s⁻¹, whereas disk stars that were splashed into halo orbits by ancient mergers such as the Gaia–Sausage–Enceladus^13,35,62 only contributed a small low-angular-momentum (L_Z ≲ 500 kpc km s⁻¹) tail (Extended Data Fig. 3). For [Fe/H] ≲ −1.2, the angular momentum shows a dominant peak at L_Z ≈ −200 kpc km s⁻¹, which is mainly contributed by the accreted Gaia–Sausage–Enceladus stars. These accreted stars survived the [α/Fe] selection because they have comparable [α/Fe] values to the in situ stars at such a low metallicity. Nevertheless, Extended Data Fig. 3 clearly shows the existence of a rotating stellar population for populations of −1.4 < [Fe/H] < −1.2 and −1.6 < [Fe/H] < −1.4. Such a low-metallicity, rotating stellar population has been revealed extensively in previous studies, and has been explained as the metal-weak tail of the thick disk populations^{46,47,48,49,50,63,64,65,66,67}. For the more metal-poor case of −1.8 < [Fe/H] < −1.6, it is harder to identify such a rotating population as the low-angular-momentum population becomes too dominant. In this case, a large portion of stars selected by the kinematic criterion L_Z > 500 kpc km s⁻¹ might be just a part of the large-radius tail of the ‘poor old heart’ of the Milky Way³⁰ or the Aurora¹³.

Modelling the density distribution of subgiant stars

We now consider different subsets of subgiant stars, each in a narrow age and abundance range, mono-age and mono-abundance population. For each of these, we reconstructed their mass density distribution at birth, \(\rho (\mathbf{r}| \tau ,[{\rm{Fe}}/{\rm{H}}];\mathbf{\uptheta })\), characterized with a set of structural parameters θ. We note that stars will never stay at their birth place. When specifying stellar mass density at birth, we need to distinguish from the stellar mass density at present, as the former reflects the underlying star formation history while the latter does not, because of mass loss due to stellar explosion.

The connection between the underlying stellar mass density and the number density found in the subgiant sample is complex and needs to be addressed carefully.

where \(\mathbf{r}:= \{l,b,d\}\), with l referring to the Galactic longitude, b the Galactic latitude and d the distance. The m_* and \({m}_{* }^{{\prime} }\) refer to the present and initial stellar mass, respectively. \(\xi ({m}_{* }^{{\prime} })\) is the stellar initial mass function (IMF), \(I({m}_{* }| {m}_{* }^{{\prime} },\tau ,[{\rm{Fe}}/{\rm{H}}])\) is the initial and present mass relationship, which can be delivered from stellar evolution models, and \(S(\mathbf{r},{m}_{* })\) is the selection function of the sample. The effective mean mass of the subgiant stars is \({\bar{m}}_{* }\), which is derived as

This expression shows that the effective mean mass of the sample’s subgiant stars, \({\bar{m}}_{* }(\mathbf{r}| \tau ,[{\rm{Fe}}/{\rm{H}}])\), depends on their distance or position \(\mathbf{r}\). This is because the survey’s apparent magnitude limit includes bright subgiants (high mass) within a far larger volume than a faint subgiant (low mass). For subpopulations selected to be mono-age and mono-abundance, this effect is greatly reduced, as the absolute magnitude M_K of subgiant stars of given age and metallicity is nearly constant.

We parametrize the spatial structure the stellar mass density distribution of each subpopulation as \(\rho (\mathbf{r}| \tau ,[{\rm{Fe}}/{\rm{H}}];\mathbf{\uptheta })\), presuming that they follow an exponential profile in both the radial and vertical directions.

where α ≔ 1/H_R is the inverse of the disk scale length. It was shown in refs. ^19,68,69 that the mono-age or mono-abundance stellar populations may show a complex radial profile, with some of them showing a density profile increasing outwards. Using α ≔ 1/H_R avoids divergences, but requires that α be positive or negative. We adopted a Galactocentric radius of the Sun's position R_⊙ of 8.178 kpc from the Gravity measurement⁷⁰. For simplicity, we also fixed the value of Z₀ to be 0. Therefore, the model contains three free parameters, \(\mathbf{\uptheta }=\{\,{\rho }_{{R}_{\odot },0},{H}_{R},{H}_{Z}\}\).

To link the stellar mass density with the stellar number density, we adopted the Kroupa IMF⁷¹ that we presumed to be universal for all ages and abundances. Recent studies have presented evidence for varying IMF with age and abundance⁷², but not in the relevant mass range 0.7 M_⊙ < M_*< 1.5 M_⊙. In our case, ignoring the variation of IMF might cause small uncertainties in the derived stellar mass density \({\rho }_{{R}_{\odot },0}\) and star formation rate. However, the scale parameters {H_R, H_Z} are not affected by the shape of the IMF.

As in ref. ², we adopted the YY isochrones to the initial-to-present mass relationship, and converted the fundamental stellar parameters {m_*, τ, [Fe/H]} to photometric and spectroscopic observables. We note that, because systematic differences exist among different isochrones, it is essential to use the same isochrones as that for the age determination to ensure precise mass estimation for the current purpose.

The sample selection function

The sample’s selection procedures are implemented in the survey’s footprint (l, b) and photometric colour–magnitude diagram (c, m) as well as subsequent cut-off in the spectroscopic T_eff-M_K diagram. The selection function thus need be further modelled as

where S_L(l, b) is the selection function accounting for the survey’s footprint, which is defined as

where S_L(c, m∣l, b) is the LAMOST survey’s selection function in the colour-magnitude diagram, and is defined as

Here, N_LAMOST refers to the number of stars that have been targeted by LAMOST, have S/N > 20, and have spectroscopic stellar labels derived with the DD-Payne⁵⁹. We adopt the Gaia eDR3 (BP − RP, G) diagram as the parent photometric colour-magnitude diagram. The (BP − RP, G) diagram is split into cells of 0.2 × 0.5 mag to build a grid of selection function for each sky field (l, b) of 2° × 2°. The selection function of each star is then interpolated from the grid.

S_sg(T_eff, M_K∣τ, [Fe/H]) is the part of the selection function that accounts for our choice of ‘subgiant’ stars, defined as a region in the (T_eff–M_K) diagram; its value is unity for stars within our adopted T_eff and M_K border, and is zero for stars outside the border.

Finally, ∣ J(c, m, T_eff, M_K; d, m_*, τ, [Fe/H], A)∣ is the Jacobian matrix that converts the fundamental stellar physical parameters to observables. It also consists of several parts

where J_s(T_eff, M_K; m_*, τ, [Fe/H]) converts the fundamental stellar parameters to the spectroscopic parameters with the YY stellar evolution models and J_p(c, m; m_*, τ, [Fe/H], d, A) converts the fundamental stellar parameters, distance, and extinction to the colour and magnitude using distance modulus

where the vectors refer to the corresponding magnitudes and extinction in the G, BP, RP bands.

Combining all these elements, we define an effective selection function,

which links the incidence of subgiant stars in our sample to the stellar mass density distribution at birth,

Also,

where ∣ J(R, Z; l, b, d)∣ is the Jacobian function for the coordinate transformation from (l, b, d) to (R, Z).

The 3D extinction map

To model the expected colour and magnitude for any envisioned subgiant stars we needed a three-dimensional (3D) extinction map of the Milky Way, F(A∣l, b, d). For the present purpose, we built such a map using the LAMOST stars. As introduced in ref. ², we derived the distance of the LAMOST DR7 stars by combining the spectro-photometric distance with the geometric distance from the Gaia parallax. We estimated the intrinsic stellar colours of the individual LAMOST stars using their spectroscopic parameters, and then derived the reddening and extinction with colour excesses (see detailed introduction in ref. ². Given the high precision of the spectroscopic parameters, the internal uncertainty in the reddening E(B − V) estimate for individual stars was typically only ~0.01 mag.

For each of the 2° × 2° field, we built an extinction-distance relationship by fitting the LAMOST stars with a series of functions of the form

where p_i are coefficients determined by the fitting. For each field, we chose the set of functions that best fit the LAMOST data as the adopted extinction-distance relationship. The E(B − V) values are converted to extinction by multiplying the total-to-selective extinction coefficients, as described in ref. ².

In Supplementary Fig. 1, we show a comparison of our 3D extinction with the extinction map of Schlegel et al. (SFD; ref. ⁷³), and with the 3D extinction map of Green et al. (Bayestar19; ref. ⁷⁴). We found remarkably good agreement with these previous maps, with a mean difference and dispersion of only −0.002 and 0.008 mag, respectively, in the E(B − V) difference between ours and the SFD map for stars with ∣b∣ > 20°, and a mean difference and dispersion of 0.013 and 0.063 mag, respectively, in the E(B − V) difference between ours and the Bayestar19 map for stars with ∣b∣ < 10°. The SFD extinction map is a two-dimensional integrated map. The good agreement with our 3D map is likely to be a reflection of the fact that most of the extinction is contributed by gas and dust in a local, thin disk. The relatively large dispersion with the Bayestar19 values for stars at low Galactic latitudes might be partially contributed by spatial inhomogeneity of the E(B − V) distribution in a sky area of 2° × 2°, that is, our 3D map may fail to accurately describe the spatially inhomogeneous extinction in these regions. However, as an uncertainty of 0.06 mag in E(B − V) estimates may only induce ~10% uncertainty in the distance modelling, it will not make a serious impact on our conclusions.

Model optimization with MCMC

The disk structural parameters θ can be estimated by taking the modelling of n_sg in equation (11) as a Poisson point process (see, for example, ref. ⁷⁵). The likelihood to be maximized for each mono-age and mono-abundance population (τ, [Fe/H]) can be expressed as

where δ refers to the Dirac function. In the first term, only the density \(\rho (l,b,d| \mathbf{\uptheta })\) depends on θ, whereas the summations (in logarithmic scale) of the other quantities are constant, and thus the equation can be further simplified^9,68. Given this likelihood function, the model parameters θ were optimized by means of an MCMC method in the Interactive Data Language (IDL) environment⁷⁶ (Supplementary Fig. 2). Only populations with N > 20 were considered here to ensure robust structure fits.

The set of parameters that yielded the maximal likelihood were adopted as the best-fitting parameters, and the parameters at 16th and 84th percentile of the integrated likelihood distribution function were adopted as the lower and upper error parameters, respectively. To have an intuitive knowledge on how good the fits are, Extended Data Fig. 4 shows the distribution of stellar distance modulus for three mono-age and mono-abundance subpopulations as examples. On the whole, the fits with exponential disk models are satisfactory.

Validating the assumption of a disk profile

Although a single exponential disk profile that we adopted (equation (4)) has been shown to be a good description for relatively metal-rich ([Fe/H] ≳ −1), high-α stellar populations^9,68, it is less clear whether such a profile provides an optimal description for populations in the low-metallicity tail of our high-α, rotating stellar sample. Therefore, we implemented a test by fitting the stellar distribution with a two-axial power-law ellipsoid halo profile

where q is the ellipticity of the halo and n is the power index for density drop from the Galactic centre. A spheroidal halo has q = 1, whereas q < 1 refers to an oblate halo and a smaller q means a profile that is closer to a disk. For the Milky Way stellar halo population, canonical values for q and n are 0.62 and 2.77, respectively⁷⁷.

We obtained a much smaller q value (0.1–0.5) but a much larger n value (~6) for all the mono-age populations than the canonical value for the stellar halo (Supplementary Fig. 3). The oldest population (15.5 Gyr) has a q value comparable to the canonical stellar halo, but the n value (~4.3) is still significantly larger. We also directly compared the radial and vertical variations of stellar density for different profiles (Supplementary Figs. 4 and 5). We found that the observed distribution for all the stellar subpopulations in our sample can be reasonably well described by the exponential disk profile, but are very different from the canonical stellar halo profile.

It has also been suggested that disk flaring phenomenon presents even for the high-α mono-abundance populations^19,69,78 or for the oldest stellar population²⁰, which may be related to its formation history⁷⁹. In the above analysis, we opted to omit the flare, as the inclusion of it would induce more complexity, considering that the H_Z/H_R is no longer an intuitive description of the disk geometry. However, to validate the effectiveness of such a simplification, we examined results for model fits that include the flaring effect,

where an exponential function f(R) is adopted to describe the flaring effect, that is, the increase of scale height as a function of R. The strength of the flaring effect is described by the parameter β.

The best-fitting parameters from this flared disk model for mono-age and mono-abundance populations are presented in Extended Data Fig. 5. The results for \({\rho }_{{R}_{\odot },0}\) (left panel), H_R (middle-left panel), and H_Z (middle-right panel) are quite similar to that presented in Fig. 1, with the exception that the value of 1/H_R now becomes larger: that is, typical scale length H_R decreased from ~2 kpc in Fig. 1 to ~1.5 kpc in Extended Data Fig. 5. This is because of the existence of the flaring effect, which exhibits a strength varying from β ≈ 0.05 for the young (τ ≈ 8 Gyr), metal-rich subpopulations to β ≈ 0.2 for the old (τ ≈ 13 Gyr), metal-poor subpopulations.

In the presence of flare, the H_Z/H_R varies with Galactocentric radius. In Supplementary Fig. 6, we present detailed comparisons between H_Z and H_R for three different radii. We found that the scale height is smaller than the scale length for most of the subpopulations at the inner part of the disk. In particular, for R ≲ 6 kpc, even the oldest stellar populations of τ > 13.5 Gyr may exhibit a clear disk geometry. As we expect that most stellar mass of the old, high-α disk populations are confined in a relatively small Galactic radius (R ≲ 6 kpc), we conclude that imposing the disk flare in the modelling does not affect the conclusions in the main text.