© The Authors 2025

1 Introduction

It is now widely accepted that the Milky Way (MW) experienced a major merger with the Gaia-Sausage/Enceladus (GSE) progenitor approximately 8–10 billion years ago, depositing a substantial population of ex situ stars into the halo (Belokurov et al. 2018; Helmi et al. 2018). In addition, observations of galaxies at cosmic noon (z ∼ 2) show that stellar discs were already in place by this time (Förster Schreiber et al. 2009; Stott et al. 2016; Lian & Luo 2024; Liu et al. 2024; Tsukui et al. 2025; Xiang et al. 2025), with similar evidence emerging from studies of the oldest stars in the MW (e.g. Conroy et al. 2022; Xiang et al. 2025). These findings suggest that the merger scattered a portion of in situ stars from the primordial disc into the halo under dynamically intense and turbulent conditions. This population of dynamically heated in situ stars in the MW is commonly recognised as the Splash (Bonaca et al. 2017; Di Matteo et al. 2019; Belokurov et al. 2020).

Yet, our understanding of the full consequences of this likely gas-rich merger for the MW’s evolutionary history is far from complete. Rapid gas accretion during the merger, along with tidal torques, may have triggered vigorous star formation by compressing gas to high densities and driving shocks (e.g. Mihos & Hernquist 1996; Barnes & Hernquist 1996). If so, stars formed in this environment may retain distinct chemical and kinematic signatures that distinguish them from those accreted or dynamically heated. Identifying these newly formed stars and decoding their signatures in today’s stellar populations are therefore crucial steps towards reconstructing the MW’s early assembly history and clarifying the significance of the GSE merger in the context of galaxy formation and evolution.

In support of this picture, An et al. (2023) present preliminary evidence for such merger-driven star formation, based on phase-space diagrams constructed using photometric metallicity estimates and Gaia astrometry (Gaia Collaboration 2023). According to the analysis, stars within the local volume (< 3 kpc) exhibit a scale-height distribution in the [Fe/H] versus rotational velocity (v_ϕ)¹ plane that appears to be sculpted by a broad, valley-like feature. This distinctive pattern – effectively traced by a coherent assembly of high proper-motion stars – stays close to v_ϕ ∼ 0 km s⁻¹ at low metallicities ([Fe/H] ≲−1.0), but rises sharply to ∼ 180 km s⁻¹ over the narrow interval −1.0 ≲ [Fe/H] ≲−0.5. The reduced scale heights relative to the surrounding volume are a telltale sign of the GSE merger, which occurred through a low-inclination, radial collision (e.g. Naidu et al. 2021). This implies that some of the stars along the valley likely originated from the same merger event, as suggested by their phase-space overlap with both the accreted GSE and dynamically heated in situ stars.

Motivated by the identification of this high proper-motion sequence (HPMS; originally termed the ‘Galactic starburst sequence’), An et al. (2023) propose that a portion of its metal-rich segment resulted from enhanced star formation triggered by a gas-rich merger. This interpretation is partly supported by recent numerical simulations of MW-like galaxies by Grand et al. (2020), who demonstrate that stars formed during merger-driven starbursts are relatively metal-rich (−1.0 ≲ [Fe/H] ≲ −0.5) and span a broad range of rotational velocities (0 ≲ v_ϕ ≲ 200 km s⁻¹), significantly overlapping with stars scattered from the primordial disc. Notably, this region of phase space aligns with the location of the metal-rich segment of the HPMS, suggesting that some of these stars may have formed during a burst of star formation triggered by the collision between the GSE progenitor and the young MW.

Additional independent support for merger-driven star formation comes from Ciucă et al. (2024), who analyse spectroscopic data from the Apache Point Observatory Galactic Evolution Experiment (APOGEE; Majewski et al. 2017) in combination with asteroseismic age estimates. They identified a clump of stars associated with the so-called ‘great Galactic starburst’ phase, characterised by rapid chemical enrichment, with [Fe/H] increasing from approximately −0.5 to −0.3 over a period of 12–13 billion years (on their asteroseismic age scale, where ages of the oldest stars exceed 18 billion years). This enrichment is accompanied by a rise in [Mg/Fe], consistent with the inflow and mixing of fresh gas during a gas-rich merger event. Together, these features provide further support for the idea that such a merger played a critical role in triggering a major starburst episode in the early MW.

However, it remains unclear how effectively and to what extent these processes enhanced shocks and/or gravitational instabilities sufficient to form a new generation of stars. A key uncertainty concerns the origin of the star-forming gas clouds. If new stars formed primarily from well-mixed primordial disc material (i.e. an in situ origin), their chemical and kinematic properties would likely resemble those of the thick-disc population. Conversely, if the star-forming material originated from the progenitor dwarf galaxy (i.e. an ex situ origin), or even from inflowing gas clouds of intergalactic origin (Renaud et al. 2021), the resulting stars would be expected to exhibit distinct chemical signatures – particularly in α-element and light-element abundances – relative to the pre-existing MW populations. Beyond the origin of the gas, it remains uncertain where these bursts of star formation occurred within the early MW. The gas-rich disc of the primordial MW may have experienced a global enhancement in star formation, or, alternatively, star formation may have been concentrated in localised regions where turbulence and gas compression were strongest (Sparre et al. 2017; Yu et al. 2021; Renaud et al. 2022).

To address these outstanding questions about the GSE merger and its role in the MW’s evolution, this study investigated the HPMS through detailed spectroscopic analyses and refined photometric metallicity estimates. Specifically, we hypothesised that the HPMS contains a significant contribution from stars formed in merger-driven star formation, alongside accreted GSE debris and dynamically heated disc stars. To test this hypothesis, we searched for stars of merger origin that exhibit distinct chemical enrichment signatures – including abundances of α-elements, sodium (Na), and aluminium (Al) – by leveraging high-resolution spectroscopic data. Through this hypothesis-driven approach, we aimed to characterise the phase-space distribution of these stars and quantify their relative fraction.

This paper is organised as follows. In Sect. 2, we introduce the HPMS and examine its presence in spectroscopic datasets. In Sect. 3, we analyse the metallicity and angular momentum distributions of HPMS stars and identify a subset of stars likely formed through merger-driven star formation, based on distinctive Na and Al enhancements. We further investigate their spatial and dynamical properties and compare them with results from cosmological zoom-in simulations. In Sect. 4, we estimate the relative fraction of this population by decomposing the metallicity distributions in spectroscopic samples, while mitigating selection biases using photometric data. In Sect. 5, we summarise our findings and discuss their implications for the MW’s evolutionary history.

2 The high proper-motion sequence (HPMS)

Figure 1 presents the distribution of high proper-motion stars in [Fe/H] versus specific angular momentum perpendicular to the Galactic plane (L_Z), which forms the basis of our analysis in this study. Panel a displays the distribution of stars based on photometry from the Sloan Digital Sky Survey (SDSS) Data Release (DR) 14 (Abolfathi et al. 2018) and the SkyMapper Southern Survey (SMSS) DR4 (Onken et al. 2024), while the remaining panels show the [Fe/H]-L_Z distributions of high proper-motion stars in large spectroscopic datasets from SDSS (Abdurro’uf et al. 2022), the Large-sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST) DR6 (Cui et al. 2012; Luo et al. 2015), the GALactic Archaeology with HERMES (GALAH) DR4 (Buder et al. 2025), and APOGEE DR17 (Majewski et al. 2017; Abdurro’uf et al. 2022). All datasets were cross-matched with the Gaia DR3 catalogue (Gaia Collaboration 2023) to obtain astrometric parameters. We included stars within 10 kpc and applied the same astrometric quality cuts across all samples, keeping only those with parallax uncertainties below 20% and proper-motion uncertainties below 30%. The Gaia DR3 parallaxes were corrected for the zero-point bias following the prescription of Lindegren et al. (2021). Further details on sample selection and the estimation of L_Z are described below.

As is illustrated in panel b of Fig. 1, the HPMS emerges as a coherent structure in phase space, consisting of stars with high proper motions that span a broad range in [Fe/H] and L_Z. Rather than representing a single stellar population, the HPMS forms a chain of previously identified components such as GSE debris (Belokurov et al. 2018; Helmi et al. 2018; Naidu et al. 2021), the Splash (Bonaca et al. 2017; Di Matteo et al. 2019; Belokurov et al. 2020), and disc stars, which together form a continuous sequence in this space. In the metal-poor regime, GSE stars naturally stand out in our high proper-motion samples, as their highly radial orbits result in large proper motions when observed locally, in contrast to the more isotropic motions of in situ halo stars. In the metal-rich regime, heated disc stars are also readily apparent, reflecting their more eccentric orbits compared to the nearly circular motions of disc stars.

Fig. 1 High proper-motion sequence (HPMS) in the [Fe/H]-LZ plane from various samples. Only stars with Zmax > 2 kpc are shown. (a) Photometric sample with proper motions exceeding 60 mas yr−1. Metallicities were estimated via Bayesian inference using SDSS and SMSS photometry, and LZ values were computed along the Galactic prime meridian without radial velocity information. (b) Schematic locations and extents of major stellar populations for reference. (c, d) Medium-resolution spectroscopic samples from SDSS and LAMOST, respectively, showing metallicities against LZ derived from full space motions. Stars with proper motions exceeding 20 mas yr−1 in SDSS and 40 mas yr−1 in LAMOST were selected. (e, f) High-resolution spectroscopic samples from GALAH and APOGEE, respectively, limited to stars with proper motions above 20 mas yr−1. For these samples, LZ values were also computed using complete kinematic information.

2.1 Photometric sample

Panel a of Fig. 1 presents the chemo-dynamical distribution of high proper-motion (> 60 mas yr⁻¹) stars in the [Fe/H]–L_Z plane, based on SDSS DR14 and SMSS DR4 photometry. As is detailed in An & Beers (2020, 2021a, b) and An et al. (2023), photometric metallicity estimates were obtained by matching individual stars to an empirically calibrated set of stellar isochrones for main-sequence stars (5.0 < M_G < 7.5), which relate photometric colours to stellar effective temperature, luminosity, and metallicity. As is shown in An et al. (2023), metallicity maps derived from SDSS and SMSS photometry offer all-sky coverage in both the northern and southern hemispheres, excluding regions along the Galactic plane, and cover a local volume extending to approximately 3 kpc.

In this work, we adopted revised photometric metallicity estimates derived from a Bayesian approach. These estimates incorporate Gaia DR3 parallaxes and a three-dimensional extinction map (An et al. 2024) ² as priors. The methodology and resulting data products will be detailed in a forthcoming publication (An et al., in prep.). Compared to our previous estimates (An et al. 2023), this approach provides a more accurate treatment of unresolved binary populations within photometric catalogues, while preserving the same metallicity scale as the previous approach and maintaining consistency with the spectroscopic scale. We selected reliable [Fe/H] estimates with uncertainties of less than 0.4 dex, focusing on stars that are likely single or are unresolved binaries with low mass ratios (< 0.5).

2.2 Low-resolution spectroscopic samples

The SDSS data were compiled from the original SDSS (York et al. 2000), the Sloan Extension for Galactic Understanding and Exploration (SEGUE; Yanny et al. 2009; Rockosi et al. 2022), the Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al. 2013), and the extended Baryon Oscillation Spectroscopic Survey (eBOSS; Blanton 2017). Low-resolution stellar spectra (R ∼ 1800) from these surveys were analysed using an updated version of the SEGUE Stellar Parameter Pipeline (SSPP; Allende Prieto et al. 2008; Lee et al. 2008a,b, 2011; Smolinski et al. 2011; Lee et al. 2013). The SSPP was used to determine stellar atmospheric parameters, including effective temperature (T_eff), surface gravity (log g), metallicity ([Fe/H]), and abundance ratios such as [α/Fe] and [C/Fe], with typical 1σ uncertainties of 180 K, 0.24 dex, 0.23 dex, less than 0.1 dex, and less than 0.3 dex, respectively.

As is demonstrated by Lee et al. (2015), the SSPP has been successfully applied to LAMOST DR6 ³, owing to its similar wavelength coverage (3800–9000 Å) and spectral resolution (R ∼ 1800) to those of SDSS. In this study, we used stellar atmospheric parameters and elemental abundance ratios – including [Fe/H] and [α/Fe]− derived with the SSPP from approximately 6 million spectra in LAMOST DR6. Among ∼44 000 stars common to both SDSS and LAMOST, the systematic differences are small: 5 K in T_eff, 0.04 dex in log g, 0.1 dex in [Fe/H], and less than 0.02 dex in [α/Fe]. As these offsets are well within the typical measurement uncertainties, we did not apply any corrections in our analysis. We refer to the approach outlined by Lee et al. (2023), who present the combined SDSS and LAMOST dataset, although their sample is restricted to stars on highly eccentric orbits.

Panels c and d of Fig. 1 present high proper-motion stars from SDSS and LAMOST, respectively. These stars, selected within 4400 < T_eff < 7000 K, have well-measured [Fe/H] values, with metallicity and [α/Fe] uncertainties of less than 0.3 and 0.1 dex, respectively, and a minimum signal-to-noise ratio (S/N) of 10. For the LAMOST sample, we adopted a propermotion threshold of > 40 mas yr⁻¹, which is less stringent than the > 60 mas yr⁻¹ cut used for the photometric sample, as the latter would lead to a substantial loss of stars. As is noted in An et al. (2023), lowering the proper-motion threshold increases contamination from disc stars, which in turn diminishes the contrast and coherence of the HPMS structure. This trade-off is further shaped by survey-specific selection effects. Accordingly, the adopted proper-motion cut represents a compromise between maximising sample size and preserving the distinctiveness of the HPMS in phase space. For the SDSS sample, an even more inclusive threshold of > 20 mas yr⁻¹ was applied, owing to the smaller number of stars tracing the vertical segment of the HPMS.

2.3 High-resolution spectroscopic samples

The GALAH sample in panel e of Fig. 1 includes [Fe/H] estimates derived from high-resolution spectra (R ∼ 28 000) from DR4 (Buder et al. 2025). The L_Z of individual stars was computed using full astrometric and radial velocity information from Gaia and GALAH. The α-element abundances ([α/Fe]) were computed as a weighted mean of [Mg/Fe], [Si/Fe], [Ca/Fe], and [Ti/Fe], including only measurements with quality flags not set (flag_X_fe). We selected stars within the temperature range 4800 < T_eff < 7000 K that have reliable [Fe/H] and [α/Fe] measurements (uncertainties below 0.3 dex and 0.1 dex, respectively, and no flag set in flag_fe_h). Additionally, we required an S/N (snr_px_ccd3) greater than 10 and ensure that measurements are not flagged in the GALAH pipeline (flag_sp, flag_sp_fit, and flag_red).

As is shown in panel (f), we also utilised accurate stellar parameter estimates, including [Fe/H] and [α/M]⁴, from APOGEE DR17 (Majewski et al. 2017; Abdurro’uf et al. 2022), selecting stars with spectra (R ∼ 22 500) having S/N > 10, and requiring that the pipeline processing flag ASPCAPFLAG is not set. To ensure reliable parameter estimates, we further restricted the sample to stars with 4400 < T_eff < 7000 K, and with uncertainties of less than 0.3 dex in [Fe/H] and less than 0.1 dex in [α/M]. Given the relatively small number of high proper-motion stars along the HPMS in these high-resolution samples, we applied the same proper-motion threshold as in the SDSS sample (>20 mas yr⁻¹).

2.4 Orbital integration

The L_Z of individual stars was computed using parallaxes and rotational velocities (v_ϕ) in the Galactocentric cylindrical coordinate system. We adopted a Galactic centre distance of 8.34 kpc (Reid et al. 2014) and the Sun’s vertical displacement of 20.8 pc (Bennett & Bovy 2019). The Sun’s motion was taken as (U, V, W) = (11.1, 12.24, 7.25) km s⁻¹ relative to the Local Standard of Rest, which has a circular velocity of v_ϕ = 238 km s⁻¹ (Schönrich et al. 2010; Schönrich 2012).

Additional orbital parameters of individual stars – including the maximum vertical excursion from the Galactic plane (Z_max), orbital eccentricity (e) and apogalacticon (r_ap) – were computed using the galpy package (Bovy 2015)⁵, adopting the MW gravitational potential from McMillan (2017). In the case of the photometric sample, orbital parameters were derived using Gaia’s proper motions along the Galactic prime meridian (within ± 30° of l = 0° and l = 180°), without relying on radial velocities, as such data were unavailable for the full photometric dataset (see An & Beers 2020). The uncertainties in the derived orbital parameters were estimated via a Monte Carlo simulation, where each orbit was integrated 500 times with input values perturbed according to the astrometric measurement uncertainties. For the stars of interest near L_Z ∼ 0 kpc km s⁻¹, the median uncertainties are σ(r_ap) = 0.1 kpc, σ(Z_max) = 0.2 kpc, σ(e) = 0.02, and σ(L_Z) = 45 kpc km s⁻¹.

For all HPMS samples in this study, we required Z_max to be greater than 2 kpc in order to reduce contamination from disc stars. The main data sets used in the following analysis comprise 9012, 9902, 2778, and 1333 stars from SDSS, LAMOST, GALAH, and APOGEE, respectively, together with 11506 stars from the photometric sample.

2.5 Biases and systematics

Various surveys often adopt complex and distinct target-selection strategies, making direct comparisons between them nontrivial. Nonetheless, the appearance of the HPMS across datasets, as shown in Fig. 1, provides a practical diagnostic for assessing the relative biases of each survey. Compared to the HPMS identified in the photometric sample, the SDSS sample is strongly biased towards GSE stars with lower metallicity and smaller L_Z. In contrast, both GALAH and APOGEE include more metal-rich stars, which better delineate the vertical segment of the HPMS, while showing a significantly weaker presence of the GSE component.

Among the spectroscopic datasets, the LAMOST sample exhibits a shape and distribution most similar to the photometric HPMS (panel a). Moreover, its overall extent and peak density in phase space are comparable to those of the photometric sample. While photometric catalogues generally contain significantly more stars than spectroscopic ones, this apparent similarity arises from specific limitations of our photometric sample. First, we restricted the sample to stars located near the Galactic prime meridian in order to convert proper motion into v_ϕ in the absence of radial velocity information. Second, our photometric metallicity estimator was calibrated only for mainsequence stars, thereby confining the sample to the local volume (< 3 kpc). In contrast, spectroscopic surveys such as LAMOST include giant stars at greater distances, enabling them to probe a more extended region of the MW with significantly larger samples.

Despite inherent sample biases, the spectroscopic samples shown in Fig. 1 confirm that the HPMS is a genuine structure, rather than an artefact arising from systematic errors in photometric parameter estimates. Furthermore, unlike the photometric sample, the spectroscopic data enable a detailed examination of the elemental-abundance distributions of stars along the sequence. In this regard, our datasets contain a sufficient number of HPMS stars to provide a robust foundation for subsequent chemical analyses. Although selection effects in the spectroscopic samples may still influence the inferred chemical trends, the photometric data offer a valuable means of evaluating and correcting for these biases.

At [Fe/H] > −1, the HPMS appears as a nearly vertical structure in all samples, as shown in Fig. 1, although its exact location varies due to differences in the metallicity scales of each dataset. For the SDSS and LAMOST samples, the median [Fe/H] of stars with 1500 < L_Z < 1750 kpc km s⁻¹, corresponding to the most prominent part of the vertical segment, is approximately −0.70, reflecting the use of the same analysis tool for both survey datasets. In contrast, the GALAH and APOGEE samples show the vertical arm shifted towards slightly higher metallicities, with [Fe/H] values of approximately −0.50 and −0.44, respectively. Notably, these values are consistent with that of the photometric HPMS (−0.49), lending support to the reliability of the photometric metallicity scale. Furthermore, all samples exhibit a nearly zero mean L_Z for GSE stars at [Fe/H] < −1. This consistency particularly reinforces the reliability of our L_Z estimates in the photometric sample, which were derived using proper motions alone.

3 Tracing stars formed in a Galactic merger

3.1 Metallicity distributions as a function of L_z

Figure 2 presents the metallicity distributions of our samples from photometry, SEGUE, LAMOST, GALAH, and APOGEE, respectively. Each dataset was binned by L_Z to highlight the varying contributions of different stellar populations across the HPMS. In the spectroscopic samples, each bin was further divided into low-α (red histograms) and high-α (blue histograms) sequences (see Appendix A). For the GALAH and APOGEE samples, we adopted an uncertainty of Δ[α/Fe] = ± 0.02 dex for the division between the two α sequences, while a slightly larger uncertainty of Δ[α/Fe] = ± 0.03 dex was adopted for the medium-resolution spectroscopic samples. Green histograms indicate the level of uncertainty in the population fractions arising from variations in the α-based division.

At high metallicity ([Fe/H] > −1.0), the spectroscopic samples show that stars on the high-α sequence (hereafter, high-α stars) clearly dominate over those on the low-α sequence (low-α stars) across most L_Z bins. The peak metallicity, corresponding to the vertical part of the HPMS in Fig. 1, remains nearly constant at [Fe/H] ≈ −0.5 to −0.7 over the range −500 < L_Z < 1500 kpc km s⁻¹, above which it tends to increase mildly. These high-α stars are associated either with the canonical thick disc or with stars formed through dynamical heating of the primordial disc. The photometric sample in the left-most column also shows a peak at [Fe/H] ∼−0.5; however, this metallicity peak becomes less prominent at low L_Z(< 500 kpc km s⁻¹), indicating that the metal-rich component has a weaker influence than in the spectroscopic samples. Such systematic differences likely arise from sample biases.

In contrast, the low-metallicity peak at [Fe/H] = −1.3 to −1.4, seen in all samples, is dominated by low-α stars, as is evident from the spectroscopic samples. These low-α stars are predominantly found in the low-L_Z regime (−500 ≲ L_Z ≲ 500 kpc km s⁻¹), consistent with the interpretation that they are mostly debris from the GSE progenitor galaxy. The contribution of low-α stars declines rapidly at L_Z > 1000 kpc km s⁻¹, but they reappear at larger L_Z, forming a rotationally supported component. These stars, however, exhibit significantly higher metallicities of [Fe/H] > −0.7 and are more closely related to the canonical thin disc, despite the Z_max > 2 kpc cut applied to the sample.

The bimodal distribution of L_Z among low-α stars is also evident in Fig. 3, which displays the distribution of [α/Fe] as a function of L_Z along the vertical arm of the HPMS (−1.0 < [Fe/H] < −0.5). For stars whose orbits are primarily confined to the inner region (r_ap < 10 kpc; top panel), the low-α stars exhibit two distinct peaks at L_Z ∼ 0 and 1800 kpc km s⁻¹, corresponding to the GSE and disc populations, respectively. A similar bimodal pattern is observed among low-α stars in the outer region (r_ap > 12 kpc; bottom panel), indicating that these populations extend far into the outskirts of the MW.

On the contrary, high-α stars with small r_ap (blue circles in the top panel of Fig. 3) maintain a nearly constant [α/Fe] across a wide range of L_Z. Notably, this sequence extends to negative rotation (L_Z ≈ −500 kpc km s⁻¹), indicating the presence of a dynamically heated population originating from a chemically well-mixed primordial disc. Meanwhile, the outer region (bottom panel), which consists of stars with large r_ap, does not prominently exhibit the Splash population. This suggests that the heated-disc population is primarily confined to the inner MW, consistent with a scenario in which the primordial disc was relatively compact when it was disrupted by a merger (e.g. Belokurov et al. 2020; An & Beers 2021b).

Based on an inspection of the [Fe/H], [α/Fe], and L_Z distributions in Figs. 2 and 3, it is evident that stars along the HPMS comprise at least three distinct and well-characterised populations:

• GSE debris: These metal-poor stars (primarily with [Fe/H] ≲ −1.0) follow the low-α sequence and are characterised by low L_Z and negligible net rotation.

• High-α population: These metal-rich stars (−1.0 ≲ [Fe/H] ≲ −0.4) span L_Z values from approximately −500 to 2000 kpc km s⁻¹. This group includes thick-disc stars as well as a dynamically heated population (often referred to as Splash), both likely sharing a common chemical origin.

• Metal-rich, low-α stars: With [Fe/H] ≳ −0.7 and high angular momentum (L_Z ≳ 1500 kpc km s⁻¹), these low-α stars may be part of the vertically extended tail of the canonical thin-disc population.

In addition to these populations, the spectroscopic metallicity distributions from GALAH (Col. (d) in Fig. 2) clearly reveal the presence of a fourth population. In the lower two panels (−500 < L_Z < 500 kpc km s⁻¹), the metallicity distribution of low-α stars (red histograms) displays an excess of metal-rich stars, with peaks at nearly the same metallicity as the heated-disc population ([Fe/H] ≈ −0.65). Their low α-element abundances, combined with relatively high metallicities, suggest a possible connection to thin-disc stars. However, their low L_Z values rule out the possibility that they are low-α counterparts of the heated-disc population. Moreover, while high-α stars consistently exhibit elevated α-element abundances across the full L_Z range, the low-α stars display a bimodal L_Z distribution (Fig. 3). Rather than being associated with the heated disc, the combination of low α-element abundances and low L_Z instead points to a closer connection with the GSE population. To further investigate their origins, we examined high-resolution spectroscopic abundances from GALAH and APOGEE in the following section, which reveal distinct chemical signatures that clearly distinguish them from the accreted GSE population.

Fig. 2 Metallicity distributions of the HPMS samples (Fig. 1) binned by LZ. In all columns except Col. (a), stars belonging to the low-α and high-α sequences are represented by red and blue histograms, respectively, while the green histogram indicates the uncertainty range due to variations in the α-based division (see text for details). To guide the eye, vertical dashed lines mark the median [Fe/H] of all stars within 1500 < LZ < 1750 kpc km s−1 in each sample. Cols. (d) and (e) refer to the high-resolution samples from GALAH and APOGEE, respectively.

Fig. 3 Distribution of [α/Fe] for GALAH stars in the vertical segment of the HPMS (−1.0 < [Fe/H] < −0.5). The top and bottom panels show [α/Fe] as a function of LZ for stars in the inner (rap < 10 kpc) and outer (rap > 12 kpc) regions, respectively. High-α stars are shown in blue and low-α stars in red. The high-α population dominates the inner region and exhibits a relatively constant [α/Fe] over a wide range of LZ. In contrast, the low-α stars display a bimodal LZ distribution: the low-LZ stars are primarily associated with GSE, while the high-LZ stars may trace the high-Zmax tail of the thin-disc population.

3.2 Detailed chemical signatures

The first line of evidence suggesting that the excess of high [Fe/H], low-α stars has an origin distinct from the directly accreted GSE population comes from detailed elemental abundance ratios in the GALAH and APOGEE samples. This is demonstrated clearly in Fig. 4, which presents [α/Fe], [Na/Fe], and [Al/Fe] as functions of [Fe/H] from GALAH. In this analysis, we included stars with upper measurement limits (shown as triangles), but excluded those with flagged abundance values (flag_X_fe). All panels show stars with |L_Z| < 600 kpc km s⁻¹, where the excess of low-α stars is most prominent in the metallicity distributions (Fig. 2).

The top panels of Fig. 4 divide the sample into low-α (left) and high-α (right) stars. Each sequence follows a relatively tight [α/Fe]–[Fe/H] relation and remains largely distinct from the other. The GSE debris is also well represented in this L_Z range, as indicated by the prominent low-α sequence. However, as is shown in the middle left panel, the low-α stars branch into two groups in the [Fe/H]–[Na/Fe] plane at [Fe/H] ≳−0.9. The main body of GSE stars has [Na/Fe] ≈ −0.2 at [Fe/H] ∼−1.3 and trends towards lower [Na/Fe] with increasing metallicity, reaching [Na/Fe] ≈ −0.3 at [Fe/H] ∼−0.7. In contrast, at −0.9 < [Fe/H] < −0.4, where the excess of low-α stars appears in the metallicity distributions, more than half of them (48 out of 69) exhibit elevated Na abundances around [Na/Fe] ∼ 0.

The red squares in the left panels of Fig. 4 highlight these low-α, high-Na (LAHN) stars, which are defined as those lying within the dashed blue box in the middle left panel. While some stars at higher [Fe/H], or even some high-α stars with systematically lower Na abundances, may be included in this selection, the [Fe/H]–[Na/Fe] cut is designed to isolate the LAHN stars and enable comparisons with the low-α, Na-poor GSE stars in the same metallicity range. Under this selection, the LAHN stars in the top left panel exhibit a subtle but systematic shift towards higher [α/Fe] relative to their Na-poor counterparts. A similar trend is seen in [Al/Fe] in the bottom left panel: although Al detections in GALAH are more challenging, with many measurements reported as upper limits, the available data reveal a population of low-α stars with enhanced Al abundances compared to other low-α stars at similar [Fe/H].

A weaker but qualitatively similar chemical trend is seen in the APOGEE sample (Fig. 5). Since APOGEE DR17 provides [Al/Fe] measurements but lacks [Na/Fe], we utilised [α/M] and the Al abundances. To guide our interpretation, we used the Al abundance range of the LAHN stars in the GALAH sample, as indicated by the dashed blue box in the bottom left panel. This selection is approximate and intended only as a reference. The stars falling within this selection box are highlighted by red squares. Although there is substantial overlap between LAHN stars and GSE debris, there is a modest but consistent indication that the LAHN stars exhibit systematically higher Al and α abundances compared to those likely associated with GSE.

Light elements such as Na and Al can be produced in asymptotic giant branch stars, but their primary synthesis occurs in core-collapse supernovae (SNe), where they are generated in greater quantities during the explosive nucleosynthesis of massive stars (see Johnson 2019, and references therein). Because core-collapse SNe are closely linked to active star formation, Na, Al, and α-elements serve as valuable tracers of star-formation history and chemical evolution. In this context, the presence of chemical bimodality in Na and Al among low-α stars, along with a weaker bimodality in the α elements, provides compelling evidence for two distinct chemical-enrichment pathways. The low-α, low-Na group corresponds to accreted stars from GSE that formed within its progenitor dwarf galaxy, where chemical enrichment was slower and less efficient. In contrast, the LAHN group likely formed during episodes of intense star formation, as indicated by their systematically elevated Na and Al abundances, along with slightly enhanced α-element levels compared to the GSE stars. These chemical signatures point to a stronger contribution from core-collapse SNe to their enrichment history.

The LAHN population shows a close connection to the intermediate [α/Fe] stars identified by Myeong et al. (2022, termed ‘Eos’). They selected stars with highly eccentric orbits (e > 0.85) from APOGEE DR17 and GALAH DR3, and identified a group with a mean [Fe/H] ≈ −0.7 and [α/Fe] ≈ +0.1 using a Gaussian mixture model. Despite considerable scatter, these stars occupy a region between the GSE and heated-disc populations in the [Fe/H]−[α/Fe] plane, consistent with our findings. Moreover, their mean elemental abundances of [Na/Fe] ≈ −0.1 and [Al/Fe] ≈ 0.0 strongly resemble ours. Furthermore, Eos stars fall within a similar angular momentum range (|L_Z| ≲ 500 kpc km s⁻¹; Matsuno et al. 2024), indicating that a fair fraction of our LAHN population belongs to the Eos group. Despite this notable overlap, our identification of the LAHN population builds on the characterisation of the HPMS (An et al. 2023), a previously unrecognised stellar sequence linked to the GSE merger. A more detailed discussion of the relationship between the LAHN and Eos populations will be presented in Sect. 5, following an examination of the spatial and orbital properties of the LAHN stars.

Fig. 4 Chemical abundances ([α/Fe], [Na/Fe], and [Al/Fe]) of the GALAH sample with |LZ|< 600 kpc km s−1. In the top panels, the dashed blue line indicates the division based on [α/Fe], separating stars into low −α (left) and high-α (right) sequences. Red squares and triangles in the left column mark low-α, high-Na (LAHN) stars. The dashed blue box in the middle left panel outlines this selection. A representative median error is shown.

3.3 Spatial and orbital characteristics

To further investigate the nature and origin of the LAHN stars, we examined their orbital properties within the MW. Figure 6 shows the distributions of L_Z, Z_max, orbital eccentricity, and apogalacticon for low-α stars in the GALAH sample, selected along the vertical arm of the HPMS (−0.9 < [Fe/H] < −0.4). The sample is divided into two groups based on their [Na/Fe] abundances: the LAHN stars (−0.2 < [Na/Fe] < +0.2; red circles) and the GSE population with [Na/Fe] < −0.2 (blue triangles). Histograms along each axis show the distributions of these two groups separately. As is discussed above (Fig. 3), low-α stars exhibit a bimodal L_Z distribution, with peaks at L_Z ∼ 0 and 1800 kpc km s⁻¹, separated by a sparsely populated region. To distinguish stars with GSE-like kinematics, we adopted a division at L_Z = 1000 kpc km s⁻¹. Stars that meet the same selection criteria as the LAHN group but lie above this L_Z threshold are shown as grey crosses. Additional plots for the high-α population are provided in Appendix B.

Three key features of the LAHN stars are evident in Fig. 6:

• They have high orbital eccentricities, comparable to those of accreted GSE stars (i.e. low-Na stars).

• They are more tightly confined to the inner MW compared to the accreted component. This is evident in their maximum vertical excursions: about 75% of the LAHN stars have Z_max < 6 kpc, whereas the corresponding fraction for the accreted stars is only 50%. Most LAHN stars also have apogalacticons within r_ap < 12 kpc, while the low-Na group is more extended, with a possible accumulation of GSE stars near r_ap ∼ 20 kpc (see Deason et al. 2018; Naidu et al. 2021).

• Among the LAHN stars, those with higher orbital eccentricities tend to reach higher Z_max. This coupling between Z_max and eccentricity reflects a shared dynamical origin of the LAHN stars and the dynamically heated population, likely tracing back to the same merger event (see also Khoperskov et al. 2023).

Since the orbital characteristics (eccentricities and L_Z) of the LAHN stars closely resemble those of stars accreted from GSE, it follows that the clouds in which the LAHN stars formed originally belonged to the progenitor galaxy of GSE rather than the primordial MW. Their closer alignment with the low-α sequence further supports this origin, indicating that these clouds were part of the dwarf galaxy. As star-forming gas clouds from the GSE progenitor lost angular momentum during the merger, they may have been funnelled into the inner MW, leading to the formation of LAHN stars with more compact orbits and tighter confinement to the disc plane than stars directly accreted from GSE. Therefore, these spatial and orbital signatures provide a second, independent line of evidence supporting their merger-driven origin.

In addition, Fig. 6 hints at a possible connection between the LAHN population and the canonical thin disc. The metallicity distributions of low-α, high-L_Z disc stars (Fig. 2) reveal a low-metallicity cut-off at [Fe/H] ≈ −0.7, coinciding with the metallicity range where LAHN stars are predominantly found. Furthermore, the LAHN stars exhibit near-solar Na and Al abundances – higher than those of GSE but lower than those of the canonical thick disc or heated-disc population. As is shown by the grey crosses in Fig. 6, these stars have orbits that are more confined to the disc, yet are more radially extended than those of the LAHN stars, consistent with an upside-down, inside-out formation of the Galactic thin disc (e.g. Bird et al. 2013). If thin-disc stars had formed exclusively from gas clouds in the primordial MW, such a coincidence in chemical and dynamical properties would be less likely.

Fig. 5 Same as Fig. 4, but displaying [α/M] and [Al/Fe] from the APOGEE dataset for stars with |LZ|< 600 kpc km s−1. The dashed blue box indicates the region in the [Fe/H]–[Al/Fe] plane where LAHN stars were identified in the GALAH sample (bottom left panel of Fig. 4).

3.4 Evidence from numerical simulations

The final piece of evidence supporting the interpretation of LAHN stars as a distinct population formed during an episode of intense star formation comes from a cosmological zoom-in simulation of a MW-like galaxy (Hirai et al. 2022), conducted using the N-body/smoothed particle hydrodynamics code ASURA (Saitoh et al. 2008, 2009). The simulation follows the chemo-dynamical evolution of satellite galaxies within a MW-mass host, with a virial mass of 1.2 × 10¹² M_⊙, making it well suited for studying the origin of stellar populations contributed by accreted dwarf galaxies. It includes key physical processes, such as metallicity-dependent radiative cooling and heating (Ferland et al. 2013), star formation based on the Schmidt (1959) law with realistic thresholds (Hirai et al. 2021), and chemical enrichment from core-collapse SNe, Type Ia SNe, and neutron star mergers (Saitoh 2017). Feedback mechanisms, such as Lyman-α heating from massive stars (Fujii et al. 2021) and turbulent metal diffusion (Hirai & Saitoh 2017), are also incorporated, allowing for realistic self-regulation of star formation. Star particles were modelled as simple stellar populations with the Chabrier (2003) initial mass function.

Figure 7 presents the spatial (left) and [Fe/H]–L_Z (right) distributions of star particles in the simulation at z = 0. Since our primary goal is to identify stars formed during mergerdriven star formation, we show only those that formed within the main halo (i.e. in situ stars). To generate mock results comparable to our observational data, we applied a spatial cut of 6 < R < 20 kpc and 2 < |Z| < 10 kpc, selecting stars within a region that encompasses the solar annulus. Unlike our observational samples, no additional selection based on proper motions was applied. In Fig. 7, the panels are divided into four age bins based on major transition periods in the galaxy’s evolution: (1) stars that formed prior to the last major merger (z ∼ 2; top row), (2) stars that formed during the bursty and episodic star formation between z ≈ 1.7 and z ≈ 1.1 (middle two rows), and (3) stars that formed after the period of intense star formation (bottom row).

By the time the last major merger occurred, approximately half of the halo mass at z = 0 had been assembled. As is shown in the top row of Fig. 7, most stars formed prior to the significant merger events originated in the primordial disc and were subsequently scattered into the halo, marking them as a heated-disc population. Their quasi-spherical spatial distribution is consistent with that observed in the MW (Belokurov et al. 2020; An & Beers 2021b). Their metallicity distribution (−1.2 ≲[Fe/H] ≲ −0.2), together with their broad range in L_Z(−500 ≲ L_Z ≲ 1500 kpc km s⁻¹), further aligns with the chemo-dynamical signatures of dynamically heated stars (see Fig. 1).

The middle two rows of Fig. 7 present the phase-space distributions of stars with ages of 7.9–9.9 billion years, subdivided into two age bins. These stars formed during bursty episodes of intense star formation triggered by the merger and supernova feedback, with successive peaks separated by a few hundred million years, highlighting the imprint of the GSE-like merger event. They are confined to a relatively narrow metallicity range (−0.6 < [Fe/H] < −0.1) and exhibit low L_Z values (|L_Z| ≲ 1000 kpc km s⁻¹), consistent with our findings for the LAHN population. Moreover, their spatial distribution at z = 0 is more compact than that of older stars and becomes increasingly flattened towards the Galactic plane for the slightly younger population. This spatial behaviour is also characteristic of the LAHN stars, which show a more concentrated distribution compared to accreted GSE stars. Interestingly, the phase-space distribution in [Fe/H]–L_Z remains nearly unchanged between the two age bins. This suggests that the star-formation process was highly chaotic, rather than one in which star formation occurred as gas clouds gradually settled into the disc.

Unlike in the observational samples, we did not classify star particles in the simulation based on their α-element abundances. The primary purpose of our α-based division was to distinguish heated-disc stars; however, such a division is unnecessary in the simulation, where heated-disc populations can be directly identified based on stellar ages (top panels in Fig. 7). Nonetheless, elemental abundances in numerical simulations provide valuable constraints on the properties of starburst populations, and will be explored in greater detail in a forthcoming study (Hirai et al., in prep.). Briefly, our simulation, which includes a detailed chemical evolution model, shows that stars formed in a mergerdriven starburst exhibit significant chemical inhomogeneities in both α–and light elements. This suggests that some stars formed directly from the ejecta of core-collapse SNe, while others were primarily enriched by Type Ia SNe from earlier generations.

Following the period of intense star formation, the simulated galaxy experienced more continuous star formation, building up the disc in an inside-out fashion, as evident in the bottom-left panel of Fig. 7. Our simulation indicates that young disc stars formed from gas enriched by inflowing clouds associated with merger events, supporting the idea that the material that gave rise to the LAHN stars may have also contributed to the formation of the MW’s thin disc. Interestingly, these young disc stars exhibit a negative correlation between L_Z and [Fe/H], mirroring the anti-correlation observed in the Galactic thin disc (e.g. Lee et al. 2011). Nonetheless, our observed HPMS sample does not show such a relation, but appears clustered in phase space at the high-v_ϕ end. This is due to the stringent proper-motion constraints imposed on the sample, which excluded high-metallicity thin-disc stars with v_ϕ values similar to that of the Sun.

We also note that our simulation does not produce a well-defined thick-disc structure, which would otherwise show a positive correlation between [Fe/H] and L_Z. This may be due to the selection criteria used for the simulated galaxy, which were based on mass and assembly history, rather than specific properties such as the presence of a thick disc. Conversely, it could be that stars in the primordial-disc component were fully dispersed into the halo by dynamical heating, leaving behind no coherent thick-disc population distinguishable in the chemo-dynamical space.

Fig. 6 Spatial and orbital properties of low-α stars with −0.9 < [Fe/H] < −0.4 in the GALAH sample. The stars were divided into two groups based on their [Na/Fe] abundances (see the middle left panel of Fig. 4): the LAHN stars are displayed as red circles, while other low-α stars with lower [Na/Fe], likely associated with the accreted GSE debris, are indicated by blue triangles. Error bars are shown only for the LAHN group. Histograms display the distributions of each group using the same colour scheme. In addition to these two groups, high-LZ(> 1000 kpc km s−1.; disc-like) stars that meet the same abundance criteria as the LAHN population are marked with grey crosses. Corresponding plots for high-α stars are provided in Appendix B.

Fig. 7 Cosmological zoom-in simulation of a MW-like galaxy. Only stars that formed within the main halo are shown, excluding those accreted from dwarf galaxies. The left column displays the spatial distributions in an edge-on X-Z projection at z = 0, divided into four age bins. The right column shows the corresponding phase-space distributions of star particles in the [Fe/H]–LZ space, selected from a region near the solar circle (6 < R < 20 kpc and 2 < |Z|< 10 kpc). Star particles with ages between 7.9 and 9.9 billion years (middle panels) originated during a burst of star formation following the major merger at z ∼ 2. Stars formed before this period (top panels) exhibit characteristics similar to the heated-disc population in the MW, whereas younger stars (bottom panels) contributed to the disc’s inside-out growth.

4 Relative fraction of the LAHN population

Having established that the LAHN stars likely formed as a result of star formation triggered by the GSE merger, we estimated their relative fraction with respect to both the accreted GSE population and the heated-disc stars. This was achieved by decomposing the observed metallicity distributions from the spectroscopic samples for the high-α and low-α groups, respectively (see Fig. 2). We focused on the lower portion of the HPMS (L_Z ∼ 0 kpc km s⁻¹), where the accreted, heated, and LAHN populations overlap, and adopted a relatively large bin size (|L_Z|< 600 kpc km s⁻¹) to ensure sufficient sample size for this analysis.

Figures 8–9 show the decomposition results for each of the four spectroscopic samples. Error bars indicate the uncertainties arising from the α-based division (Appendix A). For the high-α stars (Fig. 8), which primarily belong to the heateddisc component, we modelled the metallicity distribution over the range −1.0 < [Fe/H] < −0.1 with a single Gaussian function, shown as a solid line. While an extended metal-poor tail is evident, a significant portion of it may be due to uncertainties in the α-based classification. The peak and width of the distributions are similar across the samples, with (mean, standard deviation) values of (−0.60, 0.18), (−0.59, 0.15), (−0.69, 0.17), and (−0.67, 0.16) from panel a to d, respectively. The small differences in [Fe/H] may reflect slight systematic variations in the metallicity calibration (see Sect. 2.5).

For the low-α stars (Fig. 9), we decomposed the [Fe/H] distribution into two Gaussian components over the range −1.7 < [Fe/H] < −0.4, motivated by the presence of the LAHN stars identified in the high-resolution spectroscopic data sets (Sect. 3.2). Even without this assumption, the GALAH and APOGEE samples clearly reveal a secondary peak at [Fe/H] ≈ −0.7, associated with the LAHN stars, superimposed on a broader distribution centred at [Fe/H] ≈ −1.3, corresponding to GSE. Reassuringly, the ∼ 1: 2 population ratio between the primary (accreted) and secondary (LAHN) components over −0.9 < [Fe/H] < −0.4 aligns with our chemical identification of the LAHN stars in the [Na/Fe]−[Fe/H] plane (Fig. 4).

Following this approach, Fig. 10 presents the decomposition of metallicity distributions derived from photometric data, ranging from L_Z = −600 kpc km s⁻¹ to 1800 kpc km s⁻¹. The error bars in the histograms were computed using bootstrap resampling. Unlike the spectroscopic samples, photometric data are limited by the absence of α-element abundance measurements. However, since the LAHN population exhibits a [Fe/H] distribution nearly identical to that of the heated-disc stars, the overall distribution can still be effectively modelled with two Gaussian components: one representing accreted stars with a low-metallicity peak, and the other corresponding to the high-metallicity component that includes both heated-disc stars and the LAHN population. The fit to the photometric data using these two components is shown by the solid lines in Fig. 10. The fitting was carried out independently within each L_Z bin. Compared to the spectroscopic data, the best-fitting Gaussian components display broader distributions, primarily due to the larger [Fe/H] uncertainties in the photometric sample. At [Fe/H] > −1, the fraction of GSE stars decreases with increasing L_Z, while the combined fraction of the LAHN and heated-disc populations increases accordingly. This trend provides further support for the transition from accreted to in situ stellar populations along the HPMS.

The estimated fractions of the accreted (GSE), heated-disc, and LAHN populations in the HPMS samples are provided in the upper part of Table 1 (labelled ‘raw fractions’) and displayed in the top panel of Fig. 11. The uncertainties in these fractions were derived by propagating the fitting errors and accounting for covariances in the decomposition process. Notably, the raw fraction of the LAHN population remains relatively constant, ranging from 7 to 19% across all samples. For the LAMOST sample, we took advantage of the larger dataset to subdivide it into two L_Z bins (see Appendix C); however, no clear trend was detected in the population fractions across these bins.

The estimated contribution of the LAHN stars relative to the accreted component is significantly lower than that reported by Lee et al. (2023), who analysed the same SDSS and LAMOST datasets used in this study. They identified roughly equal proportions of the starburst (i.e. LAHN-like) and GSE populations among low-α stars with e > 0.7, but found only a negligible contribution of starburst stars among high-α stars, which were dominated by the heated-disc population. However, the population divisions in Lee et al. (2023) were primarily phenomenological in nature, based on a statistical decomposition of observed trends in orbital inclination and radial velocity dispersion with metallicity. Moreover, their analysis did not incorporate light-element abundances, and thus the stars they classified as starburst-origin may not correspond to the chemically defined LAHN population in this work. Consequently, the assignment of stars to specific populations was inherently model-dependent, leaving room for re-interpretation of the resulting population fractions.

As is shown in Table 1, the different spectroscopic samples yield significantly varying population fractions, largely due to sample selection biases. For instance, the SDSS sample shows the highest contribution from the GSE population, whereas the APOGEE sample shows the lowest. This trend is reversed for the heated-disc population. To correct for the sampling bias inherent in each spectroscopic sample, we adopted the total fraction of the heated and LAHN populations inferred from the photometric sample (25% ± 5%) as a reference for all spectroscopic datasets. The resulting bias-corrected values are presented in the lower part of Table 1 (labelled ‘normalised fractions’) and illustrated in the bottom panel of Fig. 11. Following this correction, the relative contributions of the heated and LAHN populations show good agreement across the various spectroscopic datasets. The resulting weighted mean fractions for stars within |L_Z| < 600 kpc km s⁻¹ and Z_max > 2 kpc are 19% ± 3% for the heated population and 6% ± 3% for the LAHN population in our HPMS samples.

Photometric metallicity estimates are also subject to sample biases. In the most restricted case, with volume- and mass-limited cuts (e.g. An et al. 2013), there is a known bias against metal-rich main-sequence stars, which are intrinsically fainter than metal-poor stars of the same colour and are thus underrepresented near the survey’s brightness limit. However, the sample shown in Fig. 10 is less susceptible to such biases owing to our selection of nearby, relatively bright stars with reliable astrometric measurements. Furthermore, the ratio between the heated and LAHN populations (∼ 3: 1) is expected to remain robust, as both populations exhibit similar metallicities.

Fig. 8 Metallicity distributions of high-α stars with |LZ|< 600 kpc km s−1, representing a predominantly heated population. The [Fe/H] distribution was modelled with a single Gaussian (solid line). Error bars reflect uncertainties arising from variations in the α-based division (see text).

Fig. 9 Metallicity distributions of low-α stars with |LZ|< 600 kpc km s−1. Two Gaussian components were used to model the observed distribution: a broad one representing accreted stars from GSE (dashed line), and another capturing the secondary peak linked to the LAHN population. The sum of both components is shown as a solid line. Error bars reflect uncertainties due to variations in the α-based classification (see text).

5 Summary and discussion

In this study, we investigated the origin and characteristics of the HPMS (An et al. 2023), identifying a distinct population of LAHN stars that likely formed during the GSE merger event. Using high-resolution spectroscopic data, we demonstrated that this population exhibits unique chemical signatures – particularly enhanced Na and Al abundances – distinguishing them from both the accreted stars of GSE and the heated-disc population. Their spatial and orbital properties further suggest that they formed from gas clouds originating in the GSE progenitor galaxy, which underwent rapid star formation before being mixed with clouds in the MW. Our findings are supported by numerical simulations of a MW-like galaxy, which reveal a burst of in situ star formation triggered by the merger. The simulated starburst populations exhibit spatial, dynamical, and chemical properties consistent with our observational results, reinforcing the hypothesis that the LAHN stars are direct descendants of such an event. Through decomposition of the metallicity distribution in the HPMS samples, we estimated that the LAHN population constitutes approximately 6% of the local high proper-motion stars with GSE-like kinematics, with the remaining fractions attributed to the GSE debris (∼ 75%) and dynamically heated-disc stars (∼ 19%). These results provide strong evidence that the gas-rich merger not only deposited ex situ stars into the MW, but also triggered a significant episode of in situ star formation, playing a crucial role in shaping the early MW.

Fig. 10 Decomposition of the metallicity distribution from the photometric sample. The distribution was modelled using two Gaussian components: one peaking at [Fe/H] ∼−1.3, representing the GSE debris, and a more metal-rich component reflecting a mixture of high-α (heated-disc) stars and the LAHN population. Error bars indicate uncertainties estimated via bootstrap resampling.

Fig. 11 Fractions of individual stellar populations identified in the HPMS samples. Top: raw fractions of the accreted, heated-disc, and LAHN populations, derived from Gaussian decomposition (Figs. 8–10). For the photometric sample, only the accreted component is shown, because the heated and LAHN populations cannot be separated based on photometric metallicity estimates alone. Bottom: normalised fractions of the heated and LAHN populations, scaled so that their combined contribution matches that derived from the photometric sample.

5.1 Hypothesis-driven perspective on Eos

We demonstrated that the HPMS, traced by high proper-motion stars, includes a population of stars formed during an episode of intense star formation. However, not all stars within the HPMS are direct products of this process. Rather, the HPMS traces a broader sequence of stars that were accreted, dynamically modified, or formed during the interaction between the GSE and the MW. It highlights the intricate interplay of accretion, heating, and star formation triggered by the GSE merger, offering deeper insights into how major Galactic collisions influenced the chemical and kinematic properties of stellar populations.

In this sense, our study offers a complementary perspective on the Eos group identified by Myeong et al. (2022). Using Gaussian mixture models with elemental abundances from APOGEE DR17 and GALAH DR3, combined with orbital energy, they found that Eos exhibits higher [Al/Fe] than GSE stars and is more tightly bound to the MW. Notably, they proposed that Eos bridges the GSE population at low metallicity and the low-α thin disc at higher metallicity. While the LAHN stars share many chemical similarities with Eos, our study adopted an independent approach, identifying a distinct stellar population likely linked to merger-driven star formation. The key differences relative to Eos can be summarised as follows:

• i)

  The primary distinction lies in the methodological approach adopted. Myeong et al. (2022) employed an unsupervised machine learning method to identify substructure without invoking a specific formation scenario. In contrast, our study was hypothesis-driven: we searched for stars that could plausibly have formed through merger-driven star formation by focusing on the HPMS and identifying those with chemical and dynamical properties expected for a population formed during the GSE merger event.

• ii)

  Within our hypothesis-driven framework, we identified a broader eccentricity distribution among LAHN stars, spanning e ∼ 0.55 to 1.0, with one possible outlier at e ∼ 0.35. In contrast, Myeong et al. (2022) applied a selection cut of e > 0.85, effectively restricting their sample to highly radial orbits. Our results imply that stars with similar chemical characteristics to Eos are not limited to extreme orbits, and that imposing strict eccentricity thresholds may overlook a substantial fraction of this chemically distinct population.

• iii)

  Our analysis revealed that the LAHN population occupies more spatially confined orbits than GSE, with smaller apogalacticons and vertical excursions. While Eos is also known to exhibit tightly bound orbits (Matsuno et al. 2024), our results revealed an additional layer of dynamical structure: a systematic correlation between orbital eccentricity and Z_max (see below).

• iv)

  Our study uncovered a clear bimodality in the L_Z distribution of low-α, high-Na stars. One group, referred to as the LAHN stars, is concentrated around L_Z ∼ 0 kpc km s⁻¹, while the other, with LAHN-like chemical abundances, overlaps with the high-L_Z regime typical of disc stars. The notable absence of such stars between these two L_Z groups suggests discrete episodes of star formation, likely reflecting a non-continuous inflow of star-forming gas.

5.2 The precursor to the merger-driven starburst

Our analysis of the LAHN population provided new insights into the aftermath of the GSE merger and its role in the early evolution of the MW. Although modest in number, these stars exhibit notable chemical homogeneity and intermediate abundance patterns, offering critical clues about the origin, timescale, and environment of their formation.

For the LAHN stars in the GALAH sample, the 1 σ dispersion in [Fe/H] is 0.1 dex, while [Na/Fe] and [Al/Fe] show dispersions as small as ∼ 0.1 dex and 0.2 dex, respectively (Figs. 4 and 9). The dispersion in [Fe/H] is about half that observed in the heated-disc component. Such tight chemical abundance patterns are difficult to achieve in systems with extended, multi-phase star formation, where progressive chemical enrichment or incomplete mixing of the interstellar medium could produce significant abundance spreads (e.g. Escala et al. 2018). Star formation in a turbulent disc environment is also plausible, but in that case, these stars would be expected to have more circularised orbits (e.g. Brook et al. 2004; Sales et al. 2009) rather than GSE-like kinematics with extreme eccentricities.

Instead, the observed chemical homogeneity points towards a scenario in which rapid star formation occurred from a relatively compact, well-mixed gas reservoir. This may align with the fact that the LAHN stars comprise only a modest fraction of the MW’s stellar populations. Cosmological zoom-in simulations show that stars can form in short, intense bursts during merger events (e.g. Renaud et al. 2014; Moreno et al. 2015; Ma et al. 2017; Sparre et al. 2017; Yu et al. 2021; Hirai et al. 2022; Renaud et al. 2022; Hirai et al. 2024). These bursts typically last less than a few hundred million years. Based on our estimated population fraction – roughly 1 LAHN star for every 10 GSE stars – we inferred a total stellar mass of ∼ 5 × 10⁷ M_⊙ for the LAHN population, assuming a stellar mass of 5 × 10⁸ M_⊙ for the GSE progenitor (Naidu et al. 2021). This implies a star-formation rate of the order of 0.1 M_⊙ yr⁻¹ for the gas clouds that gave rise to the LAHN stars.

The star formation rate associated with the LAHN population appears modest, likely comparable to that of the Central Molecular Zone in the present-day MW (An et al. 2011), which harbours a dense gas reservoir of ∼ 5 × 10⁷ M_⊙ (Pierce-Price et al. 2000). Still, the inferred star-formation rate remains well below the threshold typically associated with classical starbursts (∼ 10− 100 M_⊙ yr⁻¹; Kennicutt 1998). It is also significantly lower than the peak rates of MW analogues at cosmic noon, estimated to be on the order of a few tens of M_⊙ yr⁻¹ (Papovich et al. 2015). These considerations suggest that, although the star formation responsible for this population was locally significant, it did not dominate the global star-forming activity in the early MW.

Indeed, compact, high-intensity star-forming regions have also been observed in high-redshift galaxies. For example, Liu et al. (2024) reported a dusty star-forming clump in a z ∼ 1.5 galaxy using JWST and ALMA, with properties consistent with a short-duration, high-efficiency burst. Likewise, the gravitationally lensed images of a galaxy at z ∼ 1.4 (Mowla et al. 2022; Claeyssens et al. 2023) reveal compact clumps within its main galaxy that may represent either young globular clusters or sites of localised starbursts. The chemically coherent LAHN population we identified may represent a fossil record of such bursty, clump-scale star formation common at high redshift.

In our localised, burst-driven formation scenario, the LAHN stars formed from metal-enriched gas directly accreted from the GSE progenitor. This interpretation is supported by their mean metallicities, which lie at the upper end of the metallicity distribution observed in GSE stars. In addition, the increased gas density, or shocks induced by the merger, likely triggered a brief yet intense episode of star formation within the cloud, enhancing α – and light-element (Na and Al) abundances through corecollapse SNe. However, given the modest overall star-formation rate, the total yields of these elements would have been lower than those produced in the high-α disc population. This naturally explains why the elemental abundances of the LAHN stars occupy an intermediate regime between those of the GSE stars and the high-α, heated-disc population. A more intricate scenario involving pristine gas inflows from the intergalactic medium or the outer disc remains plausible (e.g. Renaud et al. 2021), implying a distinct enrichment history and potentially diverse origins. Nevertheless, the similar kinematics of these stars to the GSE debris argue against a purely outer disc or intergalactic origin.

In addition, the timing of star formation can be inferred from the orbital properties of the LAHN stars. As is shown in Fig. 6, most of them exhibit a tight correlation between orbital eccentricity and Z_max, such that more eccentric orbits tend to reach higher vertical distances. As is further demonstrated in Appendix B, high-α stars also follow a similar relation, which likely arises from stronger vertical perturbations during mergers. The fact that the LAHN stars align with this trend suggests that they formed during the same episode of dynamical disturbance, when infalling gas clouds from the dwarf galaxy collided with the primordial disc of the MW.

In the context of a merger origin for the LAHN population, Omega Centauri (ω Cen) serves as a valuable analog. ω Cen exhibits a large internal spread in metallicity, extending up to [Fe/H] ∼−0.7, and is widely regarded as the remnant nucleus of a disrupted dwarf galaxy (Lee et al. 1999; Pancino et al. 2000; Bekki & Freeman 2003). These characteristics suggest a prolonged enrichment history, consistent with multiple episodes of star formation and self-enrichment within a deep gravitational potential well. Notably, the most metal-rich stars in ω Cen, with [Fe/H] ∼−0.6 (Johnson & Pilachowski 2010; An et al. 2017), exhibit a wide range of α-element abundances, raising the possibility that some of our LAHN stars could trace the nuclear component of a progenitor dwarf galaxy.

Indeed, ω Cen members display elevated [Na/Fe] values (∼ 0.9 dex) in the highest-metallicity group ([Fe/H] > −0.9), compared to the majority of the more metal-poor stars in the cluster (Johnson & Pilachowski 2010). However, a key distinction lies in the detailed trends of the abundance ratios. While our sample of low-α stars includes a group with distinctly elevated [Na/Fe] (i.e. the LAHN group) atop a decreasing [Na/Fe] – [Fe/H] trend seen in GSE stars (Fig. 4), the stars in ω Cen follow a monotonically increasing [Na/Fe] trend with metallicity (see their Fig. 10). The [Al/Fe] abundances in ω Cen also display a more complex pattern: the metal-poor stars split into two distinct branches, while the most metal-rich stars show intermediate [Al/Fe] values, features not observed in our low −α sample (Figs. 4 and 5). These contrasting chemical trends suggest that, despite superficial similarities, the LAHN population and ω Cen followed distinct chemical enrichment pathways. This likely reflects the effect of enhanced star-forming activity in the aftermath of the GSE merger, which led to greater contributions from core-collapse SN(e) and the resulting chemical enrichment.

Following the intense but short-lived formation episode of the LAHN stars, star formation may have temporarily ceased, as suggested by the absence of low-α stars at intermediate L_Z (Fig. 3). The residual gas left behind after the merger may have funnelled into the inner MW, where it mixed with preexisting clouds in the primordial disc and eventually fuelled the formation of subsequent stellar generations (e.g. Buck 2020). A possible signature of this process is found in stars along the HPMS that exhibit similar chemical abundances to LAHN stars but have larger L_Z values (grey crosses in Fig. 6). These stars likely formed from gas clouds that acquired significant angular momentum during the merger process, possibly through tidal torques (e.g. Barnes 2002; Brook et al. 2012).

Within this framework, the starburst episode identified by Ciucă et al. (2024) may represent a later stage in a merger-driven, multi-phase star-formation history in the early MW. While their starburst sample spans a metallicity range comparable to that of our LAHN stars, it occupies a distinct kinematic regime, with stars predominantly at L_Z > 1000 kpc km s⁻¹, significantly larger than those of the LAHN population. This distinction suggests that the discontinuity in metallicity and [α/Fe] trends with stellar age reported by Ciucă et al. (2024) likely reflects a starburst event that occurred after the merger-driven gas clouds had settled into the pre-existing disc. In contrast, the LAHN stars may mark an earlier, dynamically chaotic phase of star formation, initiated by the same merger event. Taken together, these merger-driven populations may signify the transition from bursty to steady star formation in the MW, setting the stage for the subsequent, more extended phases of inside-out disc growth (e.g. Ma et al. 2017; Yu et al. 2021; Hirai et al. 2022; McCluskey et al. 2024).

6 Conclusions

In this study, we analysed a set of high proper-motion stars from photometric and spectroscopic surveys, selected to trace the aftermath of the GSE merger. Within these datasets, we identified and characterised a chemically and kinematically distinct stellar population, which we interpreted as the product of moderately intense star formation triggered by the merger. Our analysis leads to the following conclusions:

• A subset of metal-rich ([Fe/H] ≈ −0.6) stars with GSE-like kinematics exhibit low-α and high-Na abundances (LAHN stars). Their chemically distinct abundance patterns distinguish them from both the accreted GSE debris and the heated-disc population (Splash).

• While their chemical properties resemble those of the Eos population, the LAHN stars span a significantly broader range of orbital eccentricities (0.5 ≲ e ≲ 1), suggesting that Eos may represent the high-eccentricity tail of a more extensive population.

• The similarity in orbital properties (L_Z and e) between the LAHN stars and the GSE debris suggests that the LAHN stars may have formed from gas accreted from the GSE progenitor dwarf galaxy. Their present-day spatial-orbital structure, specifically the coupling between Z_max and e, also implies that they may have been subsequently displaced into the halo by dynamical heating.

• Nonetheless, the LAHN stars are more concentrated in the inner MW than the GSE stars. This finding is consistent with theoretical predictions that gas-rich mergers at high redshift can lead to the formation of stars in the inner regions of a host galaxy. Their elevated Na and/or Al abundances additionally support an in situ origin, indicating that these stars likely formed during active star formation rather than being directly accreted from a dwarf galaxy.

• The homogeneous chemical abundances of LAHN stars, together with their modest population ratio (∼ 1:10) relative to the accreted GSE stars, suggest that they formed within spatially compact gas clouds, analogous to clumpy star-forming regions observed in high-redshift galaxies.

• The modest population fraction also suggests that the star formation immediately following the GSE merger did not play a dominant role in the early assembly of the MW. Instead, the LAHN population may represent an early and localised episode of star formation, providing observational insight into the conditions that triggered starburst activity during cosmic noon.

Although the contribution of gas accreted from the GSE progenitor to the formation of the low-α disc is not yet fully constrained, our findings point to a scenario in which merger-driven star formation influenced the early assembly of the thin disc. Future work combining precise stellar ages with observed chemical and dynamical properties of stars will further clarify the connection between this merger-driven star formation and the assembly of the Galactic disc system. In particular, upcoming surveys such as the Legacy Survey of Space and Time (LSST; Ivezić et al. 2019) will enable the detection of ancient merger-driven starburst populations across a wider spatial and dynamical range, thereby offering new insight into their role in the evolution of the MW.

Acknowledgements

D.A. acknowledges support provided by the National Research Foundation (NRF) of Korea grant funded by the Ministry of Science and ICT (RS-2021-NR058093). Y.S.L. acknowledges support from the NRF of Korea grant (RS-2024-00333766). Y.H. acknowledges support from the JSPS KAKENHI Grant Numbers JP22KJ0157, JP25H00664, and JP25K01046. Numerical computations and analysis were carried out on Cray XC50 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan. T.C.B. acknowledges partial support for this work from grant PHY 14-30152; Physics Frontier Center/JINA Center for the Evolution of the Elements (JINA-CEE), and OISE-1927130: The International Research Network for Nuclear Astrophysics (IReNA), awarded by the US National Science Foundation. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss4.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics I Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatório Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University. This work made use of the Fourth Data Release of the GALAH Survey (Buder et al. 2025). The GALAH Survey is based on data acquired through the Australian Astronomical Observatory, under programmes: A/2013B/13 (The GALAH pilot survey); A/2014A/25, A/2015A/19, A2017A/18 (The GALAH survey phase 1); A2018A/18 (Open clusters with HERMES); A2019A/1 (Hierarchical star formation in Ori OB1); A2019A/15, A/2020B/23, R/2022B/5, R/2023A/4, R2023B/5 (The GALAH survey phase 2); A/2015B/19, A/2016A/22, A/2016 B/10, A/2017 B/16, A/2018 B/15 (The HERMES-TESS programme); A/2015 A/3, A/2015 B/1, A/2015 B/19, A/2016 A/22, A/2016 B/12, A/2017 A/14, A/2020B/14 (The HERMES K2-follow-up programme); R/2022B/02 and A/2023A/09 (Combining asteroseismology and spectroscopy in K2); A/2023A/8 (Resolving the chemical fingerprints of Milky Way mergers); and A/2023B/4 (s-process variations in southern globular clusters). We acknowledge the traditional owners of the land on which the AAT stands, the Gamilaraay people, and pay our respects to elders past and present. This paper includes data that has been provided by AAO Data Central (datacentral.org.au). Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences.

Appendix A Empirical division of α sequences

Figure A.1 illustrates the separation between the low- and high-α sequences for the spectroscopic samples used in this study. The [Fe/H] vs. [α/Fe] distributions are shown across three L_Z bins, selected to highlight different stellar populations and their varying contributions. The high-α sequence spans the entire L_Z range, corresponding to thick-disc stars at high L_Z and heated-disc stars at low L_Z. In contrast, the low-α sequence consists mainly of accreted GSE stars at low L_Z and thin-disc stars at high L_Z. A sparsely populated region at intermediate L_Z provided an empirical basis for defining the division between the two sequences.

Fig. A.1 Distribution of [α/Fe] versus [Fe/H] for the spectroscopic samples used in this study. The dashed red lines, derived from these datasets by examining systematic variations across three LZ bins, were used to separate the high-α and low-α sequences within each sample.

Appendix B Spatial and orbital properties of high-α stars

Figure B.1 presents the spatial and orbital properties of the high-α population (blue crosses), which were not included in Fig. 6. For reference, the same set of LAHN stars from Fig. 6 is shown as red circles. The same selection criteria were applied, including Z_max > 2 kpc and −0.9 < [Fe/H] < −0.4. Additionally, only stars with L_Z < 1000 kpc km s⁻¹ were included, where the high-α and LAHN populations overlap. As is shown in the bottom two panels, both populations exhibit a tight correlation between orbital eccentricity and Z_max, while the correlation between eccentricity and r_ap is less pronounced.

Fig. B.1 Same as Fig. 6, but comparing spatial and orbital properties of high-α stars (blue crosses) and the LAHN population (red circles).

Appendix C Sub-division of the LAMOST sample

Given the sufficient number of stars in the LAMOST sample, we subdivided the data into two L_Z bins for the decomposition of its metallicity distribution: −600 < L_Z < 0 kpc km s⁻¹ and 0 < L_Z < 600 kpc km s⁻¹. For the high-α stars, we applied a single Gaussian fit, as in the main sample (Fig. 8). No significant variation was observed in the peak metallicity ([Fe/H] ≈ −0.7) between the two L_Z bins, consistent with the overall uniformity of the high-α sequence (Fig. 2). Likewise, the low-α sample was modelled with two Gaussian components, following the procedure used for the main sample (Fig. 9). The secondary component, corresponding to the LAHN population, consistently peaked at [Fe/H] ≈ −0.8. The relative fraction compared to the accreted component was estimated to be 0.23 ± 0.13 and 0.52 ± 0.24 in 0 < L_Z < 600 kpc km s⁻¹ and −600 < L_Z < 0 kpc km s⁻¹, respectively. Consequently, the large uncertainties prevent any firm conclusion regarding variation across the two L_Z bins.

Fig. C.1 Decomposition of the metallicity distribution of the LAMOST sample, subdivided into 0 < LZ < 600 kpc km s−1 (top) and −600 < LZ < 0 kpc km s−1 (bottom), respectively. The left panels show the case of high-α stars, and the right panels for low-α stars. The solid lines represent the best-fitting results from the same decomposition procedure as in the main samples (see Figs. 8 and 9).

References

All Tables

All Figures

Fig. 1 High proper-motion sequence (HPMS) in the [Fe/H]-LZ plane from various samples. Only stars with Zmax > 2 kpc are shown. (a) Photometric sample with proper motions exceeding 60 mas yr−1. Metallicities were estimated via Bayesian inference using SDSS and SMSS photometry, and LZ values were computed along the Galactic prime meridian without radial velocity information. (b) Schematic locations and extents of major stellar populations for reference. (c, d) Medium-resolution spectroscopic samples from SDSS and LAMOST, respectively, showing metallicities against LZ derived from full space motions. Stars with proper motions exceeding 20 mas yr−1 in SDSS and 40 mas yr−1 in LAMOST were selected. (e, f) High-resolution spectroscopic samples from GALAH and APOGEE, respectively, limited to stars with proper motions above 20 mas yr−1. For these samples, LZ values were also computed using complete kinematic information.

In the text

Fig. 2 Metallicity distributions of the HPMS samples (Fig. 1) binned by LZ. In all columns except Col. (a), stars belonging to the low-α and high-α sequences are represented by red and blue histograms, respectively, while the green histogram indicates the uncertainty range due to variations in the α-based division (see text for details). To guide the eye, vertical dashed lines mark the median [Fe/H] of all stars within 1500 < LZ < 1750 kpc km s−1 in each sample. Cols. (d) and (e) refer to the high-resolution samples from GALAH and APOGEE, respectively.

In the text

Fig. 3 Distribution of [α/Fe] for GALAH stars in the vertical segment of the HPMS (−1.0 < [Fe/H] < −0.5). The top and bottom panels show [α/Fe] as a function of LZ for stars in the inner (rap < 10 kpc) and outer (rap > 12 kpc) regions, respectively. High-α stars are shown in blue and low-α stars in red. The high-α population dominates the inner region and exhibits a relatively constant [α/Fe] over a wide range of LZ. In contrast, the low-α stars display a bimodal LZ distribution: the low-LZ stars are primarily associated with GSE, while the high-LZ stars may trace the high-Zmax tail of the thin-disc population.

In the text

Fig. 4 Chemical abundances ([α/Fe], [Na/Fe], and [Al/Fe]) of the GALAH sample with |LZ|< 600 kpc km s−1. In the top panels, the dashed blue line indicates the division based on [α/Fe], separating stars into low −α (left) and high-α (right) sequences. Red squares and triangles in the left column mark low-α, high-Na (LAHN) stars. The dashed blue box in the middle left panel outlines this selection. A representative median error is shown.

In the text

	Fig. 5 Same as Fig. 4, but displaying [α/M] and [Al/Fe] from the APOGEE dataset for stars with |LZ|< 600 kpc km s−1. The dashed blue box indicates the region in the [Fe/H]–[Al/Fe] plane where LAHN stars were identified in the GALAH sample (bottom left panel of Fig. 4).
In the text

Fig. 6 Spatial and orbital properties of low-α stars with −0.9 < [Fe/H] < −0.4 in the GALAH sample. The stars were divided into two groups based on their [Na/Fe] abundances (see the middle left panel of Fig. 4): the LAHN stars are displayed as red circles, while other low-α stars with lower [Na/Fe], likely associated with the accreted GSE debris, are indicated by blue triangles. Error bars are shown only for the LAHN group. Histograms display the distributions of each group using the same colour scheme. In addition to these two groups, high-LZ(> 1000 kpc km s−1.; disc-like) stars that meet the same abundance criteria as the LAHN population are marked with grey crosses. Corresponding plots for high-α stars are provided in Appendix B.

In the text

Fig. 7 Cosmological zoom-in simulation of a MW-like galaxy. Only stars that formed within the main halo are shown, excluding those accreted from dwarf galaxies. The left column displays the spatial distributions in an edge-on X-Z projection at z = 0, divided into four age bins. The right column shows the corresponding phase-space distributions of star particles in the [Fe/H]–LZ space, selected from a region near the solar circle (6 < R < 20 kpc and 2 < |Z|< 10 kpc). Star particles with ages between 7.9 and 9.9 billion years (middle panels) originated during a burst of star formation following the major merger at z ∼ 2. Stars formed before this period (top panels) exhibit characteristics similar to the heated-disc population in the MW, whereas younger stars (bottom panels) contributed to the disc’s inside-out growth.

In the text

	Fig. 8 Metallicity distributions of high-α stars with |LZ|< 600 kpc km s−1, representing a predominantly heated population. The [Fe/H] distribution was modelled with a single Gaussian (solid line). Error bars reflect uncertainties arising from variations in the α-based division (see text).
In the text

Fig. 9 Metallicity distributions of low-α stars with |LZ|< 600 kpc km s−1. Two Gaussian components were used to model the observed distribution: a broad one representing accreted stars from GSE (dashed line), and another capturing the secondary peak linked to the LAHN population. The sum of both components is shown as a solid line. Error bars reflect uncertainties due to variations in the α-based classification (see text).

In the text

	Fig. 10 Decomposition of the metallicity distribution from the photometric sample. The distribution was modelled using two Gaussian components: one peaking at [Fe/H] ∼−1.3, representing the GSE debris, and a more metal-rich component reflecting a mixture of high-α (heated-disc) stars and the LAHN population. Error bars indicate uncertainties estimated via bootstrap resampling.
In the text

Fig. 11 Fractions of individual stellar populations identified in the HPMS samples. Top: raw fractions of the accreted, heated-disc, and LAHN populations, derived from Gaussian decomposition (Figs. 8–10). For the photometric sample, only the accreted component is shown, because the heated and LAHN populations cannot be separated based on photometric metallicity estimates alone. Bottom: normalised fractions of the heated and LAHN populations, scaled so that their combined contribution matches that derived from the photometric sample.

In the text

	Fig. A.1 Distribution of [α/Fe] versus [Fe/H] for the spectroscopic samples used in this study. The dashed red lines, derived from these datasets by examining systematic variations across three LZ bins, were used to separate the high-α and low-α sequences within each sample.
In the text

	Fig. B.1 Same as Fig. 6, but comparing spatial and orbital properties of high-α stars (blue crosses) and the LAHN population (red circles).
In the text

	Fig. C.1 Decomposition of the metallicity distribution of the LAMOST sample, subdivided into 0 < LZ < 600 kpc km s−1 (top) and −600 < LZ < 0 kpc km s−1 (bottom), respectively. The left panels show the case of high-α stars, and the right panels for low-α stars. The solid lines represent the best-fitting results from the same decomposition procedure as in the main samples (see Figs. 8 and 9).
In the text