EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 700 (August 2025) A&A, 700 (2025) A90 Full HTML
Open Access
Issue
A&A
Volume 700, August 2025
Article Number A90
Number of page(s) 23
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202453306
Published online 08 August 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

Although the exact mechanism is still debated (Sellwood & Gerhard 2020), nowadays, there is a general consensus that the bulge of the Milky Way (MW) is the product of the thickening of the bar (see, e.g. Combes & Sanders 1981; Raha et al. 1991; Athanassoula 2008; Debattista et al. 2005, 2004, 2006; Shen et al. 2010; Martinez-Valpuesta & Gerhard 2011). Since galactic bars are a natural outcome of disc instability (Ostriker & Peebles 1973; Combes & Sanders 1981; Sellwood 1981; Efstathiou et al. 1982; Sellwood & Wilkinson 1993), the MW bulge is predominantly made of the disc stars, and thus inherits their age and abundance composition ‘scrambled’ during formation and secular evolution of the bar (Zhao et al. 1994; Zhao 1996; Zoccali et al. 2006; Soto et al. 2007; Bensby et al. 2010; Babusiaux et al. 2010; Johnson et al. 2012; Ness et al. 2013b, 2014; Wegg et al. 2015; McWilliam 2016; Haywood et al. 2016a; Di Matteo 2016; Nataf 2017; Fragkoudi et al. 2017b, 2020; Johnson et al. 2012, 2022). There is no reason to believe that the boxy/peanut structures observed in external galaxies have experienced different evolutionary histories (Bureau & Freeman 1999; Erwin & Debattista 2017). Hence, the MW bulge serves as a probe for understanding general concepts of the secular evolution of barred galaxies (Kormendy & Kennicutt 2004; Gadotti 2011; Sellwood 2014).

Large-scale near-infrared observations first unveiled the peculiar morphology of the MW bulge (Okuda et al. 1977; Dwek et al. 1995), a feature that was further refined with mid-infrared WISE data, revealing its distinct X-shaped residual structure (Ness & Lang 2016). The extensive analysis of the bulge using VVV red clump stars (Minniti et al. 2010; Smith et al. 2018) has made possible the reconstruction of a complete 3D mass (baryons and dark) distribution (Wegg & Gerhard 2013; Portail et al. 2015, 2017) and kinematics (Clarke et al. 2019; Sanders et al. 2019b) in great detail. In particular, the density distribution of the Galactic bulge is characteristic of a strongly boxy/peanut-shaped bulge within a barred galaxy qualitatively very similar to that predicted by N-body simulations of the buckling instability of the disc (Raha et al. 1991; Athanassoula & Misiriotis 2002; Debattista et al. 2005; Martinez-Valpuesta et al. 2006; Quillen et al. 2014). It turned out that the chemical abundance patterns of the bulge stars can be found in populations of stars located at the solar radius (Alves-Brito et al. 2010; Bensby et al. 2010, 2017; Ryde et al. 2010; Gonzalez et al. 2011, 2015b; Ness & Freeman 2016; Haywood et al. 2016a; Bensby et al. 2017, 2020, 2021; Nandakumar et al. 2024), suggesting that the bulge is composed of disc populations that extend to the Sun. Although the disc of the MW near the solar radius is a superposition of geometric thick and thin components or α-rich and α-poor populations, respectively (Bensby et al. 2003; Reddy et al. 2003; Nordström et al. 2004; Holmberg et al. 2009; Bovy et al. 2012; Haywood et al. 2013; Minchev et al. 2013), the distinction between geometrically and chemically defined disc components has no trivial solution within the inner MW region and the individual contribution of different (chemically and geometrically defined) disc components, cannot easily be constrained, partially because of the lack of panoramic coverage of the bulge by spectroscopic data.

Observational data and models tailored to match the MW bar and bulge show that stellar populations with different kinematics are mapped differently into these components (Rich 1990; Zhao et al. 1995; Minniti et al. 1995; Soto et al. 2007; Babusiaux et al. 2010; Hill et al. 2011; Bekki & Tsujimoto 2011a; Kunder et al. 2012; Ness et al. 2013b,a; Di Matteo 2016; Zoccali et al. 2017; Wylie et al. 2021). Debattista et al. (2017) introduced this effect as kinematic fractionation, which was also observed in isolated (Fragkoudi et al. 2017b; Di Matteo et al. 2019a) and self-consistent galaxy-formation simulations (Buck et al. 2018; Fragkoudi et al. 2020). According to this mechanism, stars with hot kinematics, before the emergence of the bulge, show a less prominent boxy-peanut structure, while the coldest populations reveal the sharpest boxy-peanut morphology. This phenomenon is seen in the MW, for example, in the magnitude distribution of red giant clump stars in the MW (Nataf et al. 2010; McWilliam & Zoccali 2010; Ness et al. 2012; Rojas-Arriagada et al. 2014), where, in particular, for kinematically hot populations, the two lobes of the peanut appear at higher latitudes compared to the colder populations. This, however, requires the association of kinematics with the chemistry (or geometry) of stellar populations, which is not so straightforward because some simulations (see, e.g. Fragkoudi et al. 2017a) suggest that in the outer parts of the boxy/peanut, the hot population tends to have a higher line-of-sight velocity, while in the inner regions, this trend can be reversed, depending also on the orientation of the bar with respect to the line of sight. Such a picture seems to have been confirmed observationally by Ness et al. (2012) and Rojas-Arriagada et al. (2014) using ARGOS and GES data, respectively, but is not so certain in Ness et al. (2016) with the APOGEE data.

Since the early studies, there has been no consensus on whether the bulge is dominated by metal-rich (Rich 1988; Whitford & Rich 1983; Rich et al. 2007; Schultheis et al. 2015) or old and metal-poor stellar populations (Ortolani et al. 1995; Feltzing & Gilmore 2000; Zoccali et al. 2003). With the arrival of data from massive spectroscopic surveys, the detailed metallicity distribution function (MDF) exploration of the bulge became possible (Zoccali et al. 2003, 2008; Bensby et al. 2013; Ness et al. 2013a; Bensby et al. 2017; Schultheis et al. 2017). Traditionally, the bulge MDF has been described as a mixture of several Gaussian components with varying contributions depending on the field location (Ness et al. 2013a; Rojas-Arriagada et al. 2014, 2017). Using data from the Gaia-ESO survey, Zoccali et al. (2017) proposed a bi-modal bulge MDF (see, also Babusiaux et al. 2010; Hill et al. 2011; Gonzalez et al. 2013; Rojas-Arriagada et al. 2014, 2017; Zoccali et al. 2017; Schultheis et al. 2017; Zoccali et al. 2018; Kunder et al. 2020). Rojas-Arriagada et al. (2020) identified three MDF components in the bulge using the APOGEE DR16 data, while Wylie et al. (2021) found that the bulge is best fit by four Gaussians, where the most metal-poor one is not very significant, thus agreeing with the three dominant components. Using datasets that differ completely in size and spatial coverage, Ness et al. (2013a) (ARGOS) and Bensby et al. (2013, 2017) (microlensed dwarf and subgiant stars) identified five bulge MDF components whose locations match each other reasonably well. These components were tentatively associated with different stellar populations spanning the metallicity range from ≈−2 to ≈0.5: a thin boxy/peanut-shaped bulge, a thicker boxy peanut-shaped bulge, the pre-instability thick disc, the metal-weak thick disc and the stellar halo. Nevertheless, all the above-mentioned studies agree that the MDF varies significantly depending on the position in the bulge region. Also, the resulting MDF may be affected by the intrinsic properties of a given survey, such as the metallicity scale, the selection function, and the survey footprint (Wylie et al. 2021). The latter is extremely patchy for spectroscopic surveys as the bulge stars are very hard to observe due to high extinction towards the Galactic centre, especially close to the mid-plane. Therefore, obtaining unbiased spectroscopic data for the bulge region is challenging and remains the focus of several future high-resolution MW disc surveys like 4MOST (Bensby et al. 2019), SDSSV (Almeida et al. 2023) and especially MOONS (Gonzalez et al. 2020), which will deliver high-signal-to-noise spectra in the near-infrared (Cirasuolo et al. 2020). Therefore, so far, the question of whether the bulge is a metal-rich or metal-poor component of the MW galaxy has remained disputable (Rojas-Arriagada et al. 2020; Fragkoudi et al. 2020).

It is believed that the bulge contains a (very) metal-poor component associated with a spheroidal or classical bulge, characterised by its lack of rotation (Zoccali et al. 2008; Dékány et al. 2013; Kunder et al. 2016; Savino et al. 2020; Arentsen et al. 2020; Barbuy et al. 2023). However, the mass contribution of this component appears to be quite small, accounting for only a few percent of the total bulge mass (Shen et al. 2010; Kunder et al. 2012; Ness et al. 2013a; Kunder et al. 2016; Di Matteo et al. 2014; Debattista et al. 2017). The presence of a spheroidal component can be confused with a contamination of the stellar halo. Outside the Galactic centre, the stellar halo is made of a mixture of accreted (Gaia-Sausage-Enceladus merger debris Belokurov et al. 2018; Haywood et al. 2018b; Helmi et al. 2018) and heated disc stars (Splash/Plume in-situ populations Belokurov et al. 2020; Di Matteo et al. 2019b) that, as was recently demonstrated, have different kinematics and chemical abundances; hence, they can be disentangled from the in-situ spheroid showing pre-disc chemical composition and a lack of rotation (Aurora, pre-disc in-situ populations Belokurov & Kravtsov 2022, 2023), unless it was spun up into cylindrical rotation (Fux 1997; Saha et al. 2016; Saha & Gerhard 2013; Saha et al. 2012; Di Matteo et al. 2015; Saha et al. 2012; Di Matteo et al. 2014). However, we need to underline that the discoveries regarding accreted and in-situ stellar halo components are mostly based on the samples of stars heavily weighted towards a solar radius with very little knowledge or predictions on how these components should look like in the inner Galaxy, where the chemo-kinematic MW composition is governed by an extremely complex boxy/peanut bulge and bar interplay. Therefore, understanding the MW bulge is crucial not only for constraining its formation and the present-day structure but also for helping us to isolate its contamination on the structure of accreted systems buried in the heart of our Galaxy.

This paper, part of a series, aims to describe the present-day composition of the MW bulge using the orbit superposition approach and the APOGEE data. In previous papers, we have introduced the method using mock MW data (Khoperskov et al. 2025b) (hereafter Khoperskov et al. 2025a) and demonstrated its advantage in the reconstruction of the survey footprint-independent chemo-kinematic structure of the MW disc (Khoperskov et al. 2025a) (hereafter Khoperskov et al. 2025a). We also showed that the APOGEE DR 17 giant stars sample is already sufficient for obtaining the unbiased stellar mass-weighted MDF together with morphology and stellar kinematics across the entire MW. The paper is structured as follows. In Section 2 we briefly describe the methodology and the APOGEE dataset, which were deeply covered in previous papers in the series. Section 3 focuses on the spatial structure of the bulge recovered using our orbit superposition approach. The kinematics and orbital composition of the bulge populations are presented in Section 4. Section 5 covers the elemental abundance structure, gradient spatial variations, MDF, and decomposition of the bulge population using Gaussian mixture models (GMMs). We summarise our main findings in Section 6.

2 Data and method

2.1 APOGEE data

In this work, we use the giant star sample from the APOGEE DR 17, selecting stars without problematic flags and according to surface gravity. The details of the selection can be found in the second paper of the series (Khoperskov et al. 2025a) from which the input sample of stars is adopted. We use radial velocities, atmospheric parameters and stellar abundances ([Fe/H] and [Mg/Fe]) from the APOGEE DR17 (Abdurro’uf et al. 2022), which were complemented by the proper motions from the Gaia DR3 catalogue (Gaia Collaboration 2023b). We use only stars with radial velocity uncertainties better than 2 km s−1, distance errors <20% and proper motion errors better than 10%, as this is more critical to calculate stellar orbits. In order to cover a larger area across the MW disc, we select giant stars with logg < 2.2 and flagged them as ASPCAPFLAG = 0 and EXTRATARG = 0. Our final selection includes approximately 80 000 stars spanning the MW disc, with around 13 000 of these located in the bulge region, which is comparable to the sample sizes used in recent bulge studies with APOGEE (Rojas-Arriagada et al. 2020; Boin et al. 2024). Although the raw number of stars in the bulge is modest, the use of an orbit superposition technique effectively amplifies the sampling by nearly two orders of magnitude, resulting in a significantly enhanced coverage of the bulge region.

The selection of giant stars from the APOGEE catalogues, despite their stellar parameters being less precise compared to those of dwarfs, is driven by the necessity to cover a larger area of the MW (see the comparison in Imig et al. 2023). As has been demonstrated in our previous work (Paper I), this broader coverage is crucial for the effective application of the orbit superposition method. Here, we also do not remove stars with relatively high abundance uncertainties from our sample, as they are being used to propagate the chemical information along the orbits of the stars.

In the present analysis, we do not include stellar age information, as it may lack sufficient precision for smaller samples in the inner disc, particularly for low-metallicity and high-α populations. Additionally, ages do not provide significant new insights due to the well-established tight age-metallicity relation observed in high-resolution spectroscopic datasets of the bulge (Bensby et al. 2013; Haywood et al. 2016a; Bensby et al. 2017; Haywood et al. 2018a), hence the metallicity can be safely used to trace the structure populations of different age. In our case, the age catalogue from Stone-Martinez et al. (2024) excludes recommended ages for stars with [Fe/H] < −0.7 dex, which represents about 5% of the total sample and ≈13% of the inner (<3.5 kpc) MW mass (see Fig. A.1). Although the broad age range (≈6–12 Gyr) and its distribution are consistent with the above-mentioned studies, we remain cautious about delving into the finer details of the bulge’s age structure. Furthermore, we have examined several large age catalogues in our orbit superposition analysis (Leung & Bovy 2019; Leung et al. 2023; Queiroz et al. 2023) and found quality issues affecting a substantial number of stars, limiting their utility in exploring the MW bulge.

2.2 Orbit superposition method

In this work, we use the output of the orbit superposition results obtained in the previous paper of the series (Khoperskov et al. 2025a), while the method is described in detail in Paper I. Hence, we briefly mention the key steps of the approach. We adopt the 3D mass distribution of the MW, including its stellar component, from Sormani et al. (2022), which is an updated analytic model of the potential constructed by Portail et al. (2017). This analytic potential, available in AGAMA (Vasiliev 2019), shows the correct behaviour of the mass distribution outside the bar region and reproduces well the 3D density of the bar, including the X-shape structure of the bulge (Wegg & Gerhard 2013; Wegg et al. 2015).

We integrated orbits of the APOGEE stars, assuming a constant bar pattern speed of 37 km s−1 kpc−1, in agreement with various different studies (Portail et al. 2017; Bovy et al. 2019; Sanders et al. 2019a; Khoperskov et al. 2020; Clarke & Gerhard 2022; Khoperskov & Gerhard 2022). The weights of the orbits in the rotating rest frame were calculated by adjusting their total 3D density to the analytic solution for the stellar component from Sormani et al. (2022). Each orbit was transformed to 500 phase-space coordinates1, and chemical abundances for these data points were assigned using normal distributions using the uncertainties as a width of the distribution. In other words, each orbit can be considered as a sample of stars following each other along the orbit with similar stellar parameters.

We remind the reader that our approach heavily relies on the assumption about the dynamical equilibrium of the MW disc. While there is evidence that this is not entirely true near the Sun (see, e.g. Siebert et al. 2011; Williams et al. 2013; Gaia Collaboration 2018; Antoja et al. 2018) and in the outer disc of the MW (see, e.g. Poggio et al. 2018; Gaia Collaboration 2021; Poggio et al. 2025), so far, we have no clues of ongoing disequilibrium associated with the ongoing X-shaped structure evolution or the bar growth. The apparent evolution of the bar may include spiral-bar interconnection (Hilmi et al. 2020; Vislosky et al. 2024) resulting in a short timescale of the bar parameters, as they are quasi-periodic, depending on the beat frequency of spirals and the bar; however, the time-averaged properties of the bar remain unchanged. A possible ongoing deceleration of the MW bar caused by the interaction with the DM halo (Debattista & Sellwood 1998, 2000) can be more problematic for our model. For instance, Chiba et al. (2021) estimated the current bar slowing rate as of ≈ −4.5 km s−1 kpc−1 Gyr−1, potentially resulting in ≈10% (≈ 300 − 450 pc) increase in the bar size over the last Gyr. Note, however, that for the ≈9 Gyr old MW bar (Sanders et al. 2024; Haywood et al. 2024), such a rapid recent transformation is rather unlikely, as N-body models clearly show much slower bar parameters change 1–2 Gyr after its formation (Debattista & Sellwood 1998, 2000; Athanassoula 2003; Łokas 2019; Fragkoudi et al. 2021). Finally, the deceleration rate was obtained by Chiba et al. (2021), who used a pure m = 2 bar to explain the existence of the Hercules stream via bar deceleration; the latter, however, is not needed to shape the solar neighbourhood kinematics if the bar density distribution also accounts for an m = 4 component (Hunt & Bovy 2018). Considering the high level of uncertainties about the time-dependence or ongoing evolution of the MW bar parameters, we can safely assume that the inner MW is in dynamic equilibrium, obviously apart from the rigid rotation of the bar itself, which we consider in the present work.

thumbnail Fig. 1

On sky distribution of the APOGEE stars in our initial sample (left) and stellar mass-weighted density projection obtained by superposition of orbits of these stars (right). The top panels cover the entire disc of the MW, while the bottom ones zoom in on the MW bulge region. The orbit superposition technique not only reconstructs the 3D density distribution of the MW disc and bulge but also maps the abundance patterns of APOGEE stars across these components, effectively eliminating the limitations caused by the survey’s footprint and correcting for selection function biases. The reconstructed density shows a prominent boxy/X-shaped bulge structure, which is asymmetric along the longitude due to the bar’s orientation.

3 Bulge structure and orbital composition

3.1 3D density distribution in the inner MW

In order to demonstrate the advantages of the orbit superposition approach in Fig. 1, we compare the obtained projected stellar density map with the footprint of the APOGEE DR 17 in the Galactic coordinates (l, b). The top panels show the entire data sample and total MW density, while the bottom ones are zoomed into the bulge region, within <3.5 kpc from the Galactic centre. Since APOGEE is a near-IR survey, it covers well the inner MW region but still does not capture the regions very close to the midplane (|b| < 2–3°); moreover, the spatial footprint is very patchy and represented by several isolated observed fields. On the contrary, in our approach, we are able to fill all these gaps by ‘painting’ the abundance data from APOGEE along the orbits in the bulge, providing the full 3D chemo-kinematic picture of the inner MW, and effectively correcting the spatial footprint of the survey. The importance of the complete MW bulge coverage is vital as its chemo-kinematic characteristics change a lot as a function of the sky coordinates (see, e.g. Di Matteo 2016; Barbuy et al. 2018; Gerhard 2020, and reference in Introduction), not surprisingly, as stellar populations experience very complex dynamical evolution (see, e.g. Quillen et al. 2014; Di Matteo et al. 2014, 2019a; Fragkoudi et al. 2017a; Fragkoudi et al. 2018, 2019, 2020; Buck et al. 2019; Debattista et al. 2020, 2023; Ciambur et al. 2021; Beraldo e Silva et al. 2023; Anderson et al. 2024; Boin et al. 2024).

The bottom-right panel of Fig. 1 illustrates the structure of the MW bulge, as reconstructed from the orbits of real stars in the APOGEE survey. This reveals the characteristic boxy and the X-shaped bulge structure, which has been previously detected through the split in the magnitude distribution of red clump stars along specific lines of sight (see, e.g., McWilliam & Zoccali 2010; Nataf et al. 2010; Saito et al. 2011) and mapped in full 3D by Wegg & Gerhard (2013). At the same time, the density contours highlight the boxy structure and the X-shaped component, which is very prominent further away from the midplane, b > 5°. Two lobes of the bulge are asymmetric in longitude relative to the Galactic centre because of the bar orientation relative to the Sun-Galactic centre line, making the nearby lobe more vertically extended than the farther one.

As the bar and bulge density structure varies a lot as a function of the position in the disc, in Fig. 2, we present stellar column density maps in (l, b) coordinates at varying distances from the Sun. From top to bottom of Fig. 2, the panels display density slices that capture the asymmetry in the disc caused by the Galactic bar and bulge. The maximum density shifts from left to right as we slice the bar at different distances from the Sun. At around −4 kpc from the Galactic centre, the near side of the bar becomes visible as an overdense region with a moderate thickening, shifted by approximately l ≈ −20° relative to the Galactic centre. As we approach −2.7 kpc from the centre, the bulge reveals a distinctive diamond-like structure aligned with the bar, which appears most prominent at around −1.3 kpc from the Galactic centre. The middle panel demonstrates the symmetric, boxy component of the inner bulge. On the far side (three bottom panels), the bar and bulge structure mirrors itself across the longitude axis, appearing less extended in the latitude direction due to geometric effects.

In Fig. 3 to highlight the ability of our orbit superposition method in the reconstruction of the details of the density distribution, we present the one-dimensional stellar density distributions along various lines of sight towards the Galactic bulge. The top left panel illustrates the positions of 5 × 11 fields (l = −15° to 15°) selected at five latitudes (b = 0°, 3°, 7°, 10°, 15°), overlaid on the stellar density map. Notably, two density peaks are already visible at b = 15°, emphasising the vertical extent of the X-shaped bulge structure. While shifting towards each other, these peaks become more pronounced at lower latitudes and are predominantly observed for l > 0° (represented by shades of red lines). The increase in separation between the peaks of density for greater heights above the plane has been used as a manifestation of the peanut structure of the MW bulge (Ness et al. 2012; Rojas-Arriagada et al. 2014; Gonzalez et al. 2015b; Di Matteo 2016; Gómez et al. 2016) as it is shown in various MW-like barred galaxy simulations (Shen et al. 2010; Debattista et al. 2017; Portail et al. 2017; Fragkoudi et al. 2017b; Buck et al. 2018). The picture we observe here is consistent with other observations finding a single peak density distribution near the midplane (Babusiaux & Gilmore 2005; Rattenbury et al. 2007; Pietrukowicz et al. 2012; Gonzalez et al. 2013) and becoming bimodal at higher latitudes (McWilliam & Zoccali 2010; Saito et al. 2011; Wegg & Gerhard 2013; Ness et al. 2013a).

We conclude that our method successfully recovers the full 3D density distribution of the MW bulge, which is expected given the proven strength of orbit superposition techniques in capturing fine details of bars and bulges (Li & Shen 2012; Wang et al. 2012, 2013; Vasiliev & Valluri 2020; Tahmasebzadeh et al. 2024). What sets our approach apart is its use of parameters from resolved stellar populations, enabling us to explore the bulge’s structure based on the chemical composition of its subpopulations (Poci et al. 2019; Zhu et al. 2020, 2022).

thumbnail Fig. 2

Variation of the inner MW density structure along the line-of-sight reconstructed using orbit superposition and APOGEE data. From top to bottom, the panels show the stellar density in (l,b) coordinates in 0.6 kpc-wide slices with increasing distances from the Sun, as marked in the bottom left of each panel. The middle panel shows the boxy bulge structure at the Galactic centre, while the upper and lower panels show the near and far sides of the bar-bulge, respectively. As we slice the MW bar inclined relative to the Sun – Galactic centre line, the off-centred X-shaped structure transitions from negative to positive longitudes with increasing distance.

3.2 Orbits of the MW bulge

The presented bulge density structure provides us with a very important conclusion that the orbit superposition approach allows us to recover the details of the MW bulge where stars from the APOGEE DR 17 giants sample precisely reproduce the orbital backbone of the X-shaped structure. A detailed description of the orbital families of the X-shaped structures in simulated systems is given in Parul et al. (2020), while here we focus on the dominant components that we find in the reconstructed bulge of the MW. Using the standard techniques (Pfenniger & Friedli 1991; Athanassoula 2003; Parul et al. 2020), we calculated frequencies along the Cartesian coordinates (fx and fz) and cylindrical radius (fr) for all orbits passing through the inner 3.5 kpc. In Fig. 4, we show the in-plane frequency ratio distribution for the APOGEE stars whose orbits have masses obtained in the orbit superposition modelling. The most prominent peak, around fr/fx ≈ 2, highlights the orbits supporting the bar (see, e.g. Patsis & Katsanikas 2014; Portail et al. 2015; Patsis & Athanassoula 2019), while the rest can be considered as the disc as stars on these orbits do not follow the elongated density distribution along the bar’s major axis (see more details in Parul et al. 2020). Based on the orbital frequencies, the definition of the bar yields ≈11% of the total MW stellar mass, or ≈5.86 × 109 M⊙.

Next, we decomposed the bar-like orbits according to the vertical to in-plane frequency ratio, fz/fx, shown in the bottom panel of Fig. 4. The stellar mass-weighted distribution roughly follows the ones presented in Portail et al. (2015); however, in our case, the peak at fz/fx ≈ 2, corresponding to the banana-like orbits (Pfenniger & Friedli 1991; Patsis et al. 2002; Skokos et al. 2002) is more prominent. We find that about 36% of the bar mass is represented by orbits of this class. Unlike Portail et al. (2015), who analysed only the bar orbits (fr/fx ∈ 2 ± 0.1) in a narrow range of vertical-to-in-plane orbital frequencies (1.4 < fz/fx < 2.1), we present the full orbital composition of the entire orbital frequency ratio range of the bar. In Fig. 5 we show the face-on and side-on density structure of the orbital families in different ranges of the fz/fx ratio. Similarly to Portail et al. (2015) one can see that the X-shaped bulge component is dominated by the orbits with 1.4 < fz/fx < 1.9 (first four columns) corresponding to the pretzel-like orbits. In this case, the X-shaped structure has very sharp edges (see also Patsis & Katsanikas 2014; Patsis & Harsoula 2018) with a total mass of ≈1.7 × 109 M⊙. Higher-order orbital members (1.9 < fz/fx < 3) represent the peanut component, where the banana-like (fz/fx ≈ 2) orbits dominate. The highest-order orbits correspond to the edges of the bar major axis. In our case, the fraction of the X-shaped orbits is somewhat lower, ≈30% versus 44% found by Portail et al. (2015). This discrepancy comes from a higher mass of the peanut and the long bar components, which was increased in Portail et al. (2017), providing a better agreement for the bar pattern speed, and propagated to the analytic potential we adapt in our modelling. We remind the reader that the masses of the bulge components correspond to the orbital-based definition of the bar (fr/fx ∈ 2 ± 0.1) and do not include the non-resonant disc component or higher-order periodic orbits embedded into the inner Galaxy, which in fact are far less significant building blocks of bars (Patsis & Athanassoula 2019).

thumbnail Fig. 3

Line-of-sight density structure of the MW bulge reconstructed using orbit superposition and APOGEE data. The top left panel shows the selection of the bulge fields marked by circles of different colours with the total column stellar density on the background. Other panels show the stellar density distribution along the line-of-sight at different latitudes, as is marked in the top left panel. The colour of the lines matches the colour of the circles in the top left panel. Only fields with b ≥ 0 are shown as the data are symmetric relative to the midplane. An animated version is available online.

4 Bulge line-of-sight kinematics

As was discussed above, the intricate 3D density distribution and orbital structure of the MW bulge naturally lead to complex kinematic behaviour among its stellar populations. The known kinematic trends in the MW bulge region include a cylindrical rotation with the decrease and flattening of the velocity dispersions with vertical distance from the disc plane (Howard et al. 2009; Kunder et al. 2012; Ness et al. 2013b; Zoccali et al. 2014). However, the details of the stellar kinematics vary with the metallicity of stellar populations (Fragkoudi et al. 2017a; Debattista et al. 2017; Di Matteo et al. 2019a; Debattista et al. 2019; Boin et al. 2024), allowing various scenarios of MW bulge formation to be tested.

Before exploring the metallicity-dependent variations in bulge kinematics, we first present the overall trends. In Fig. 6, we display the projected, averaged line-of-sight velocity and velocity dispersion for stars within a heliocentric distance of 8.12 ± 3.5 kpc, covering the entire bulge-bar region. The kinematic maps reveal strong rotational motion, with the strength of rotation decreasing at higher latitudes. At any given latitude, populations located further out in the disc have higher line-of-sight velocities. The kinematic disc component is clearly present, as is emphasised by the narrowing of iso-velocity contours towards the midplane. In the innermost region (l < 20°), the velocity isocontours gradually broaden towards higher latitudes. The velocity dispersion map highlights the presence of colder disc populations near the mid-plane, showing a trend with latitude. This trend is in the sense that the velocity dispersion decreases with increasing latitude. The most striking feature in the bulge region is the peak velocity dispersion, reaching up to ≈150 km s−1 near the Galactic centre, rapidly decreasing and saturating around 100 km s−1 at latitudes |b| > 5°. The picture we see here is in agreement with data from the BRAVA, ARGOS and APOGEE surveys, which have shown that the bulge is cylindrically rotating (Howard et al. 2008; Ness et al. 2013b, 2016), which is the typical velocity field produced by bars. In this case, the velocity dispersion peak in the very centre does not imply the presence of the classical bulge, as it can be explained by a pure buckled bar models (Shen et al. 2010; Bekki & Tsujimoto 2011a; Li & Shen 2012; Saha & Gerhard 2013; Di Matteo et al. 2015).

While maps in the (l, b) perspective provide us with the projected kinematics typically available in observations, the structure of the bulge as a function of distance can help us understand the 3D motions of stars in the inner MW. Figure 7 shows the mean line-of-sight velocity and velocity dispersion in several latitude bins (0–5, 5–10, 10–15 and 15–20°). Such information is hard to obtain in the existing observational datasets, as it requires a large-scale spatial coverage, but more importantly, precise distance and radial velocity determination far away from the Sun. In the case of our approach, these quantities are known from the orbits of APOGEE stars and do not blur the observed picture. The kinematic maps presented in the figure show, first of all, large-scale disc rotation, which in the centre is affected by the presence of the bar. Since the bar is inclined relative to the Sun–Galactic centre line, we can see a prominent tilt of the radial velocity associated with the streaming motion along its major axis, often observed in extragalactic IFU data (Aguerri et al. 2015; Gonzalez et al. 2016; Molaeinezhad et al. 2016; Guo et al. 2019; Kacharov et al. 2024) and simulations (Debattista et al. 2015, 2019). As we explored in Khoperskov et al. 2025a, such a velocity tilt is difficult to obtain using the raw Gaia or APOGEE data, as the precision of the distance determination is not sufficient in this region (Gaia Collaboration 2023a). The radial velocity tilt in the bar region is the most prominent at latitudes close to the midplane, and the pattern fades away with increasing distance from the midplane.

In Fig. 7, the line-of-sight velocity dispersion face-on maps also reveal several peculiarities in the inner bulge region. The velocity dispersion has a peak near the Galactic centre (<0.5 kpc), but it disappears rapidly at latitudes higher than b ≈ 5°. Such a feature is easy to explain as this region is dominated by the shortest pretzel-like orbits, which have very complex shapes (see Fig. 5). Once stars with different orbital phases are mixed, they result in a very high line-of-sight velocity dispersion. Interestingly, the velocity dispersion decreases along the bar major axis with local minima on the edges; the latter can be understood as there is not much diversity in the orbital composition of the bar’s edge, which is dominated by orbits with the highest fz/fx ratio (see Fig. 5). The radial velocity dispersion plateau, observed along the minor bar axis, corresponds to the X-shaped component.

Next, we proceed with a more conventional analysis of the MW bulge kinematics that is widely used in the literature. In Fig. 8 we show the line-of-sight velocity and velocity dispersion longitude profiles for different latitude slices (different colour) and also separated for populations of different metallicity. For reference, in each panel, the total velocity and velocity dispersion profiles for all stars located inside the bulge (<3.5 kpc) are shown with black lines. Although the MW bulge MDF is known to possess several components, the metallicity bins were chosen arbitrarily for the current analysis, while the complete MDF decomposition is investigated in Section 5.3. The following trends can be observed in Fig. 8:

  • In all the panels, the rotation is cylindrical within |l| < 10°, with very little difference in rotation speed at different latitudes, except for the highest latitude (b = 10–15°, dark blue lines) showing a shallower velocity gradient in this longitude range. This behaviour is observed in all metallicity bins, except for the most metal-poor one [Fe/H] < −1, showing almost negligible difference in rotation as a function of latitude. Such uniform cylindrical rotation is observed in many models where the bulge is produced via vertical instability of the bar (Shen et al. 2010; Fragkoudi et al. 2017a; Debattista et al. 2019). We note that almost no rotation is found for the very metal-poor bulge stars (below [Fe/H] ≈ −1.5; Arentsen et al. 2020), which the APOGEE sample does not account for; however, we estimate that the total mass of this population is less than 1%.

  • For a given metallicity range (except [Fe/H] < −1), the velocity profiles diverge at larger longitudes, l > 10, showing the decrease in the velocity (or shallower velocity gradient) with increasing latitude, as expected due to asymmetric drift. The velocity profiles are nearly symmetric around l = 0, with a velocity being slightly higher for l > 0 compared to l < 0 (see also Cabrera-Lavers et al. 2007; Rangwala et al. 2009; Vásquez et al. 2013; Zoccali et al. 2014). The latter is likely the result of the projection effect, where the far side of the bulge is seen at lower latitudes (see Figs. 2 and 7). A weak longitudinal asymmetry of the velocity dispersion can also be seen in the bottom panels.

    thumbnail Fig. 5

    Orbital decomposition of the MW bulge in face-on (top) and side-on (bottom) projections. The panels show the projected stellar density of orbits classified by their vertical-to-in-plane frequency ratio, fz/fx (see Fig. 4). The bar’s major axis is aligned along the X-axis in both projections. Each orbital class represents a distinct dynamical population, and the stellar mass associated with each class is indicated in the top left corner of the corresponding panels.
    thumbnail Fig. 6

    Kinematics of the MW bulge recovered using orbit superposition approach. The panels show the mean line-of-sight velocity (top) and the line-of-sight velocity dispersion (bottom) in the MW bulge region, limited by a heliocentric distance of 8.12 ± 3.5 kpc, around the Galactic centre. The grey contours reflect the shape of the stellar density distribution, while the white ones correspond to the constant velocity and velocity dispersion values. The bulge region shows a prominent rotation due to its disc origin, while the asymmetry of the radial velocity profile is due to projection effects. The velocity dispersion map shows a prominent peak extended perpendicular to the midplane.

  • There is a trend in increasing velocity outside l ≃ 10° with increasing metallicity from [Fe/H] = −1.5 to ≈0 dex, which then is weakly reversed, as is seen for the higher latitudes. In a scenario where the low and high-metallicity bulge populations correspond to the thick and thin discs, respectively, this is explained if thin disc stars lose more angular momentum once trapped by the boxy/peanut bulge than thick disc populations (Fragkoudi et al. 2017a) (see more details and comprehensive analysis in Boin et al. 2024).

  • The velocity dispersion gradually decreases with increasing metallicity, as one can expect, assuming a reasonable age-velocity dispersion relation and an increase in metallicity with age. This, however, does not mean that the observed velocity dispersion remained unaffected by the bulge emergence, as suggested by Debattista et al. (2017) the differential impact of the vertical bar instability on pre-existing disc populations results in their kinematic fractionation, preserving but not completely erasing the ‘initial’ kinematics (see also Fragkoudi et al. 2017a; Di Matteo et al. 2019a).

  • For a given metallicity range, at high latitudes (b > 5–6°) the velocity dispersion profiles are almost flat. At latitudes < 5 we can see that the velocity dispersion increases towards l = 0°; however, thanks to the panoramic coverage of the reconstructed MW bulge data, inaccessible for the sparse raw spectroscopic datasets, we can see a prominent flattening of the profile around l ≈ 10°, seen in simulations of MW-like bulge galaxies (Debattista et al. 2017). This flattening of the radial velocity dispersion profile reflects the plateau we discussed in Fig. 7, which we attribute to the X-shaped component.

  • The central velocity dispersion peak, as was discussed above, is caused by the inner boxy bulge component. Interestingly, its peak value slightly increases from solar towards very metal-rich populations, suggesting that they are trapped more efficiently in the boxy/peanut instability, which can be the case if they were initially colder compared to less metal-rich populations (Debattista et al. 2017; Fragkoudi et al. 2017a). These findings provide an interpretation of the recently rediscovered metal-rich knot in the innermost MW (Rix et al. 2024).

thumbnail Fig. 4

Orbital frequency distribution of the MW bulge populations. The top panel shows the stellar mass-weighted distribution of the in-plane orbital frequency ratio, fr/fx. The elongated orbits, fr/fx ≈ 2, correspond to the bar, while the rest do not support the bar and represent the disc. The bottom panel shows the distribution of the vertical-to-radial orbital frequencies, fz/fx. The frequency ratio fz/fx ≈ 2 corresponds to the banana-like orbits. The lower values (fz/fx ∈ 1.4–2) are pretzel-like populations, dominating in the X-shaped bulge structure, whose extent along the bar major axis decreases with decreasing frequency ratio. Orbits with fz/fx > 2 correspond to the longer bar component (>2.5–3 kpc) outside the bulge.

5 Chemical abundance trends

Various studies suggest that the MW bulge does not consist of a single component but is a superposition of multiple components (see Introduction) or perhaps a continuum of stellar populations with different chemical and kinematic properties shaped by both varying physical conditions of the corresponding star formation and redistribution during dynamical evolution of the bar and bulge. In this section, using the stellar mass distribution obtained in the orbit superposition reconstruction, we dissect the MW bulge in terms of its chemical abundance profiles, gradients, and distribution functions.

thumbnail Fig. 7

Kinematics of the inner MW and bulge in the face-on projection. The panels show the mean heliocentric line-of-sight velocity (top) and velocity dispersion (bottom) maps at different latitudes. The grey contours show the total stellar isodensity levels, marking the extent and orientation of the bar. The dashed lines highlight the longitude values of l ± 20°. The yellow asterisks mark the position of the Sun.
thumbnail Fig. 8

Line-of-sight kinematics of the MW bulge. The panels show the mean line-of-sight radial velocity (top) and the line-of-sight velocity dispersion (bottom) versus longitude measured at different latitudes and in [Fe/H] bins, as marked in the header of the top panels. The line colour indicates the latitude range. The panels are based on the data within <3.5 kpc around the Galactic centre. In each panel, the black lines show the averaged velocity and velocity dispersion profiles for all stellar populations of the entire region.

5.0.1 Bulge MDF

In Figure 9, we compare the raw APOGEE data with the reconstructed stellar mass-weighted distribution in the [Fe/H]-[Mg/Fe] plane for the region within <3.5 kpc. The white line divides the distribution into high- and low-α populations, following the method used in the analysis of the MW disc (Khoperskov et al. 2025a). The relative mass distributions of the high- and low-α populations along the [Fe/H] and [Mg/Fe] axes are represented by the blue and red lines, respectively, while the total distribution functions for [Mg/Fe] and [Fe/H] are shown in yellow. The MW bulge closely resembles the chemical abundance pattern of the MW disc, which is expected, as the bulge is composed of thin (mostly low-α) and thick (mostly high-α) disc stars, as was discussed earlier (see, e.g. Bensby et al. 2010; Gonzalez et al. 2011; Johnson et al. 2014; Bensby et al. 2014).

The fractional contribution of the high- and low-α components in the bulge region is reversed compared to what is observed across the entire MW disc. In particular, around 40% of the bulge consists of low-α populations, whereas approximately 44% of the total MW mass is high-α, based on the same definition (Khoperskov et al. 2025a). It is important to note that this reversal is not due to differences in radial density profiles, as the monoabundance populations in the bulge have similar scalelengths of ≈2 kpc (Khoperskov et al. 2025a). Additionally, the orbit-superposition-based mass fraction for the bulge shows a higher proportion of high-α stars compared to the ratio derived from the raw APOGEE data, where about 32% is low-α. This discrepancy arises because our approach allows us to account for regions of the bulge near the midplane, which are underrepresented in the APOGEE footprint (see Fig. 1). Therefore, we conclude that the [Mg/Fe] distribution is bimodal with relatively close fractions of low- and high-α populations (see also Rojas-Arriagada et al. 2019).

Although the stellar mass-weighed MDF (right panel in Fig. 9) is very wide, it is characterised by two main components, centred at −0.7 and +0.3 dex, corresponding to the centres of distributions of the high- and low-α disc populations. Several weak modulations are observed between these metallicity values, potentially indicating the presence of several components claimed by Ness et al. (2013a); Bensby et al. (2010, 2011, 2013, 2017). Before delving into the details of the bulge MDF, we provide more comparisons with the raw APOGEE data. In Fig. 10, we show the MDF of the bulge stars inside R < 3.5 kpc, limiting our spatial selections by different latitude values.

Overall, we notice that the global shape of the MDFs for the raw APOGEE and stellar mass-weighted datasets very closely resemble each other and can be broadly described as bimodal (Hill et al. 2011; Uttenthaler et al. 2012; Bensby et al. 2011; Gonzalez et al. 2015b,a; Zoccali et al. 2017; Rojas-Arriagada et al. 2014, 2017; Schultheis et al. 2017). However, there are several small but significant differences. In all the panels, the metal-rich peak is shifted to slightly lower values in the orbit superposition results than the raw APOGEE. Near the midplane (b < 1.7–5°), the APOGEE MDF shows 3–4 peaks (Ness et al. 2013a), while the mass-weighted method fills the gaps between them. The best match we observe is for the middle right panel (|b| < 7.5°). Once we consider populations at higher latitudes, the APOGEE MDF lacks a number of metal-poor stars, which is seen in the bottom panels. We conclude that the correction of the MDF provided by the orbit superposition is rather modest, highlighting a very good APOGEE selection function in the bulge region, weakly affecting the metallicity distribution relative to the underlying mass distribution (Rojas-Arriagada et al. 2017; Zasowski et al. 2017; Nandakumar et al. 2017).

Using the ARGOS spectroscopic survey of multiple fields in the bulge region, Ness et al. (2013a) identified five components in the bulge MDF and associated them with different stellar populations, with decreasing metallicity: A – a thin, boxy/peanut-shaped ([Fe/H] ≈ +0.15); B – a thicker, boxy, peanut-shaped ([Fe/H] ≈ −0.25); C – the inner thick disc ([Fe/H] ≈ −0.7); D – the metal-weak thick disc ([Fe/H] ≈ −1.2); and E – the inner stellar halo ([Fe/H] ≈ −1.7). The five bulge MDF components were identified by Bensby et al. (2013) using microlensed dwarfs with a good agreement with the ones from Ness et al. (2013a) (see also Bensby et al. 2017, for the most recent results and comparison). The agreement of these studies, based on different datasets with different spatial coverage, supports the existence of multiple distinct bulge populations with different spatial distributions and kinematics. In the case of the orbit superposition, while we notice some small amplitude modulation between the main high- and low-α components, it is not obvious whether these MDF features can be associated with the components known in the literature.

Since we have identified stars in the APOGEE dataset whose orbits maintain the MW bulge boxy/peanut structure, it is natural to associate corresponding orbital families with the MDF features. In different panels of Fig. 11 we present the MDF of the orbits in several orbital frequency ratio ranges (see Fig. 5), representing the inner boxy orbits (fz/fx = 1.4–1.6), X-shaped orbits (fz/fx = 1.6–1.9), banana orbits (fz/fx = 1.9–2.2) and long bar orbits (fz/fx > 2.2). The red vertical lines highlight the location of the MDF components from Bensby et al. (2017). The grey histogram shows the MDF of the disc stars, i.e. stars with orbital frequencies which do not support the bar, |fz/fr − 2| > 0.1.

We find that all orbital families have broad metallicity distributions with multiple peaks whose location differs from one orbital family to another. Some of the peaks match the position of components from Bensby et al. (2017). We need to keep in mind that the metallicity scale of the APOGEE can be different from the one in Ness et al. (2013a) (ARGOS) or Bensby et al. (2017) samples (Rojas-Arriagada et al. 2017). Nevertheless, we discuss the features of individual orbital families and the differences between them.

The inner boxy orbits are the only ones whose MDF matches reasonably well almost all peaks from Bensby et al. (2017), and it is the only orbital family which shows a small peak near [Fe/H] ≈ −1. Ness et al. (2013a) associated this peak (D) with the thick metal-weak disc present in the MW before the bulge formation. The next metallicity peak, around [Fe/H] ≈ −0.7, is seen for all orbital families supporting the bulge. However, its significance decreases from inner boxy to X-shaped orbits, but it is still somewhat present in the banana orbits and has a very small contribution in the long bar orbital families. This MDF peak is the most dominant one, corresponding to the centre of the high-α populations, or inner thick disc (component C in Ness et al. 2013a). However, it is also seen for stars not supporting the bar orbits (see grey histogram). The two following MDF peaks (≈−0.2 and ≈+0.1 dex) are also seen among all bulge-supporting orbital families. The most metal-rich MDF peak matches the best by the X-shaped and banana orbits; interestingly, these orbits capture more metal-rich stars compared to the background disc whose metal-rich MDF peak is visibly shifted towards lower metallicities by ≈0.15 dex (see, e.g. Buck et al. 2018).

To summarise, we demonstrated that the MDFs of different orbital families supporting the MW bulge indeed reveal multiple peaks, some of which can be associated with the MDF components found in the literature. On the contrary, the background disc orbits MDF is dominated by the two main rather broad components associated with high- and low-α populations of the bulge.

By construction, the 3D density distribution of the MW bulge is defined by the adopted potential and the precision of its recovery using the orbit superposition. However, our approach does not control the dependence of the bulge density structure (and kinematics) for populations with different abundances. Therefore, it is very interesting to see how different stellar populations are mapped into the bulge. We conclude this section by presenting the MW bulge density structure as a function of metallicity. In Figure 12, we show the projected stellar density distribution in (l, b) coordinates for six metallicity bins for stars inside 3.5 kpc from the Galactic centre. For consistency, we use the same metallicity bins for which we analysed the bulge kinematics (see Fig. 8).

  • The most metal-poor bin (−1.5 < [Fe/H] < −0.8) shows an oval-like distribution elongated along the midplane but with no clear signature of the discy or the X-shaped component. These populations are the most metal-poor stars in our APOGEE sample, and they are likely contaminated by the halo (in-situ and accreted) stars (Di Matteo et al. 2019b; Belokurov et al. 2020) whose ages are typically >10–11 Gyr. Moreover, their initial hot orbits make it hard for them to be captured by the X-shaped structure, which should appear not earlier than the MW bar itself, which is ≈ 8–9 Gyr old (Sanders et al. 2024; Haywood et al. 2024). Nevertheless, as they show little rotation (see Fig. 8), they are still dominated by the metal-poor, early MW disc potentially corresponding to the spin-up phase (Belokurov & Kravtsov 2022; Semenov et al. 2024; Zhang et al. 2024).

    thumbnail Fig. 12

    Shape of the MW bulge stellar density across different metallicity bins. The corresponding metallicity range and stellar mass are specified at the top of each panel. The top three panels represent the high-α populations, while the bottom panels depict the low-α populations. The mapping of these populations into the bulge reveals distinct trends where the contrast of the X-shaped/boxy structure increases with metallicity. High-α stars, corresponding to a kinematically hot thick disc component, are more widely dispersed in the bulge. In contrast, low-α stars are part of a colder population that is more closely confined to the midplane but still efficiently mapped into the bulge. The most (extremely) metal-rich populations exhibit a vertically extended, highly concentrated distribution towards the centre of the MW without displaying a prominent disc-like morphology. An animated version is available online.

  • The following two metallicity bins (−0.8 < [Fe/H] < −0.5 and −0.5 < [Fe/H] < −0) correspond to the bulk of the high-α sequence or the inner thick disc. This is reflected by their large vertical extension. Both populations demonstrate the boxy component and the X-shaped structure is more prominent for the more metal-rich bin.

  • Two supersolar metallicity bins (metal-rich 0 < [Fe/H] < 0.25 and very metal-rich 0.25 < [Fe/H] < 0.5) represent the bulk of the low-α stars and, as expected, the initially thin disc populations are more tightly confined to the midplane. However, towards the Galactic centre, the thickness increases sharply, making rise to the bulge, where the X-shaped structure is the most prominent among all populations. Interestingly, the very metal-rich populations appear to be thicker than the metal-rich ones. This difference could be attributed to the vertical bar instability, which impacts the colder (very metal-rich) populations a bit more efficiently, pushing them further from the midplane. In contrast, the metal-rich stars were already kinematically hotter at the time of bulge formation, making them less susceptible to this effect.

  • The last, extremely metal-rich populations (EMR, [Fe/H] > 0.5) reveal a very compact, boxy structure without a strong signature of the X-shaped component. These stars show a very thin and compact disc component outside the bulge region (Rix et al. 2024) and, while being the most metal-rich population of the MW formed about 5–9 Gyr ago (Khoperskov et al. 2025a), can be formed only in a disc-like configuration with little radial redistribution over time. The appearance of these compact populations as a vertically thick component simply reflects their influence by the vertical bar instability, even if they formed after the emergence of the main X-shaped/boxy component (Debattista et al. 2019). Therefore, this apparent spheroid-like shape of the EMR stellar populations distribution does not suggest the presence of a classical bulge.

thumbnail Fig. 9

[Mg/Fe]-[Fe/H] relation for the MW bulge. The left panel shows the distribution of APOGEE stars, while the right panel depicts the mass-weighted distribution obtained using the orbit superposition method within <3.5 kpc from the Galactic centre. The yellow lines show the distribution functions of [Mg/Fe] and [Fe/H] separately. The solid white line shows the border used to separate high- and low-α populations. The number of stars and stellar mass corresponding to the high- and low-α populations are marked in the panels. While the orbit superposition method assigns more mass to the metal-rich populations, both APOGEE data and the orbit superposition-based density reveal rather similar abundance distributions.
thumbnail Fig. 10

Comparison of the MW bulge MDF limited by different distances in the vertical direction from the midplane as marked in each panel. The blue line is the MDF weighed using the orbit superposition method versus the raw APOGEE sample counts in grey. The orbit superposition method and APOGEE show similar MDFs close to the midplane; however, some systematic differences become prominent once the input sample covers regions further out.
thumbnail Fig. 11

MDFs of different orbital families of the MW bulge. The orbital families are classified according to the in-plane to vertical frequencies ratio: inner boxy orbits (fz/fx = 1.4–1.6), X-shaped orbits (fz/fx = 1.6–1.9), banana orbits (fz/fx = 1.9–2.2) and long bar orbits (fz/fx > 2.2) orbits. In each panel, the distributions are normalised to the maximum value. For reference, the vertical red dashed lines mark the locations of the peaks in the MDF based on the microlensed dwarf and subgiant stars in the MW bulge from Bensby et al. (2017). The same grey histogram in the background shows the MDF of the disc orbits with |fz/fr − 2| > 0.1.

5.1 Abundance variations and gradients in bulge

The determination of the abundance gradients and spatial variations in the bulge is fundamental because its exact shape can be related to several different bulge formation channels. The presence of the vertical metallicity gradient within the inner 2 kpc was demonstrated by Minniti et al. (1995). Later, Rich et al. (2007) found no evidence of major iron abundance or abundance ratio gradient between the inner field and Baade window, suggesting the lack of the vertical metallicity gradient. On the contrary, Zoccali et al. (2008) demonstrated the existence of a vertical metallicity gradient along the bulge minor axis and the photometric metallicity map of Gonzalez et al. (2013) indicates a strong vertical gradient of ≈ −0.28 dex kpc−1. The secular MW-like bulge formation models suggest that the metallicity gradients are inevitable because of the pre-existing radial and vertical metallicity gradients in the MW disc(s) stellar populations, which are then modified, depending on the kinematics of corresponding populations, by the X-shaped structure but cannot be erased (Debattista et al. 2017).

The comparison of the exact values of the bulge metallicity gradient is challenging as different surveys cover different regions of the bulge with varying depths along the line of sight. With the orbit superposition data, we are no longer affected by the survey spatial footprint and can analyse the abundance variations across the whole bulge region. In Fig. 13, we present the maps of the mean [Fe/H], [Mg/Fe] and low-α mass fraction. The latter is useful for understanding the contribution of the main bulge (disc) components in different regions. The mean metallicity and [Mg/Fe] maps show very complex distributions, which, in fact, follow the density distribution (grey background contours). The metallicity is higher in the densest regions close to the midplane and decreases away from it, implying a negative metallicity gradient. The central region of the bulge seems to have a flat metallicity distribution with the mean metallicity around ≈−0.1, which is somewhat lower compared to the regions outside the bulge (|l| > 20°). This might give the impression that the bulge is more metal-poor compared to the surrounding disc populations or that there is a positive radial metallicity gradient in the inner MW. This behaviour has been demonstrated in several spectroscopic studies demonstrating that the MW bulge has a significant metal-poor component at all latitudes (|l, b| < 10°), which leads to the MW bulge being on average metal-poor (Ness et al. 2013a; Ness & Freeman 2016; Rojas-Arriagada et al. 2017; Zoccali et al. 2017). We argue, however, that this is a projection effect because, although the radial metallicity profile is flattened in the inner 3–4 kpc, the radial gradient is still negative along the bar major axis (see Fig. 21 in Khoperskov et al. 2025a). This is supported by the Auriga simulations (Fragkoudi et al. 2020) where boxy/peanut bulges are predominantly metal-rich, in agreement with IFU spectroscopic studies of the inner regions of external galaxies (Gonzalez et al. 2017; Neumann et al. 2020) and, hence, the MW is a typical barred galaxy.

At higher latitudes (|b| > 2.5°), the metallicity variations along longitude are not monotonic (see Fig. 13). The metallicity drops once we move inwards (along longitude) from the disc to the bulge region, and it rises again once we cross the bulge. This is the result of the vertical thickening of the metal-rich low-α populations whose fraction is about 50% at these latitudes in the bulge. Notably, the low-α stars are also expected to thicken but at a higher rate compared to the colder metal-rich populations, as predicted in the secular bulge formation scenario. On top of the non-monotonic metallicity variations, we observe an asymmetric behaviour relative to l = 0. This is caused by the bar orientation where the nearby lobe is seen as a more vertically extended metallicity distribution at l < 0 compared to the outer lobe (l > 0).

In Figure 13, the [Mg/Fe] ratio variations show an inverted pattern compared to the [Fe/H] profiles, which is a trivial result coming from the [Mg/Fe]-[Fe/H] relation where metal-poor stars have higher [Mg/Fe] values and vice versa (see Fig. 9). To assess the metallicity gradient variations across the bulge in Fig. 14 we show the 2D maps of the radial (along longitude) and vertical (along latitude) slopes of the metallicity profiles measured in each pixel within ±1 deg. One can see that the vertical or longitudinal metallicity gradients in the MW bulge cannot be characterised by a single value. Thus, it is unsurprising that surveys covering various parts of the bulge arrive at different solutions. The vertical gradient is flat in the very centre (l < 5 and |b| < 2). Further out, we observe the increase in the vertical negative metallicity gradient value with both latitude and longitude. The maximum values of the vertical gradient are observed near the edge of the X-shaped structure, between the regions dominated by low- and high-α populations (see rightmost panels in Fig. 13). The longitudinal metallicity gradients also change with the position, tracing the X-shaped bulge component well.

The existence of the metallicity gradients in the bulge is predicted by N-body and hydrodynamical simulations of secular bulge formation (Martinez-Valpuesta & Gerhard 2013; Fragkoudi et al. 2017b, 2020). The question is whether this is a pure effect of the mapping in the vertical direction of horizontal (Martinez-Valpuesta & Gerhard 2011) or vertical (Bekki & Tsujimoto 2011b; Fragkoudi et al. 2017b) metallicity gradients initially present in the disc, or a combination of the two (Di Matteo et al. 2015, 2019a). Answering this question will require reproducing the abundance maps presented in this section in simulations. However, we cannot exclude certain degeneracy between different scenarios, especially if both thin and thick MW discs existed and both had radial and vertical gradients before the X-shaped structure formation. The question of how young stars populated the boxy peanut structure after its emergence is poorly explored in simulations (but see Debattista et al. 2019), but it is exciting and may help to provide additional constraints on the origin of the MW bulge and conditions in the MW at the time of its emergence.

5.2 MW bulge: Metal-rich or metal-poor

In this section, we address a simple but rather debatable issue in the literature: whether the MW bulge is metal-poor or metal-rich. Since the earliest photometric studies, the MW bulge was traditionally considered a predominantly metal-rich component (see, e.g. Rich 1988; Feltzing & Gilmore 2000; Ortolani et al. 1995; Fulbright et al. 2007). However, with the advent of large-scale spectroscopic surveys, it has become evident that the bulge also harbours a substantial metal-poor population, present across all latitudes within l < 10° (see, e.g. Ness et al. 2013a; Ness & Freeman 2016; Rojas-Arriagada et al. 2017; Zoccali et al. 2017; Fragkoudi et al. 2018). This metal-poor component has been highlighted, for instance, in Bovy et al. (2019) and Queiroz et al. (2021), where the central region of the MW is shown to be relatively metal-poor. In this context, the bulge metallicity should be considered relative to the surrounding disc populations. If it is metal-poor, the MW would be pretty unusual among the nearby barred galaxies with pseudo-bulges (Neumann et al. 2020, 2024) and contrary to what is suggested in modern galaxy-formation simulations (Buck et al. 2019; Debattista et al. 2020; Fragkoudi et al. 2020), suggesting that centres of bars are metal-rich. The metal-poor MW bulge would also suggest that the radial metallicity distribution of the MW disc is not typical and has a broken profile, as reported in Lian et al. (2023); Bovy et al. (2019). However, in the previous work of the series, we demonstrated that this is not entirely true once the selection function and masses of mono-abundance populations are taken into account (Khoperskov et al. 2025a). We found that the radial metallicity gradient of the MW, although flattened in the inner region, is still negative, and the MW bar appears to be metal-rich.

The MW bulge metallicity structure is hard to analyse in full 3D, as the existing spectroscopic surveys have limited coverage in the inner MW; however, it is not an issue for the orbit superposition reconstructed MW. In Fig. 15, we show the mean stellar metallicity in the face-on disc projection (XY), considering different distances from the midplane. In other words, the metallicity is averaged in horizontal slabs with increasing (from left to right) thickness from 0.1 kpc to a given value indicated in each panel. Close to the midplane (|z| < 0.1–0.5 kpc, two leftmost panels), the mean metallicity is low near the Galactic centre compared to the surrounding disc region, as expected. However, when we consider populations farther away from the midplane, an inverted trend emerges, revealing the bulge as a metal-rich component. This behaviour can be explained by assuming that the bulge consists of both a thin, metal-rich population and a thicker, metal-poor one. The deficit of metal-rich stars near the midplane in the bulge region results from vertical bar instability, which causes these populations to be ‘puffed up’ above the plane. In contrast, the more metal-poor stars, which are kinematically hotter, remain less affected by this instability during bulge formation and dominate at smaller distances from the midplane while experiencing less vertical thickening. As we move away from the plane, we capture the metal-rich stars that were displaced from the midplane, increasing the average metallicity in the bulge region compared to the surrounding disc. This inversion in the metallicity profiles is seen along a longitude in Fig. 13. The explanation we present here clarifies the findings of Queiroz et al. (2021), who claimed that the MW bulge is metal-poor (see their Fig. 7). This conclusion is naturally influenced by the fact that their APOGEE sample is biased toward the midplane and does not account for the APOGEE footprint or selection function properly. As a result, Queiroz et al. (2021) obtained a metallicity distribution similar to the leftmost panels in Fig. 15, which, as we now understand, does not present a complete picture.

While the MW bulge has slightly subsolar metallicity overall, it still appears as a metal-rich component, similar to (pseudo) bulges observed in external galaxies (Moorthy & Holtzman 2006; Sánchez-Blázquez et al. 2014; Breda & Papaderos 2018; Neumann et al. 2024). However, in certain projections – such as in Galactic coordinates (l, b) shown in Fig. 13 – the bulge can appear metal-poor due to these projection effects and also because it is being compared to the metal-rich thin disc stars in the midplane. In Fig. 15 we also see that kinematically colder, metal-rich stars are more concentrated along the bar’s major axis, producing the metal-rich bar. The latter, however, does not imply the metal-rich stars formed along the bar major axis as the star formation is usually suppressed in the bar region (Khoperskov et al. 2018b; George et al. 2019; Géron et al. 2021; Fraser-McKelvie et al. 2020). The latter is seen as the lack of recent significant star formation in the inner (<4–5 kpc, but excluding the central molecular zone) MW as well (see, e.g. Elia et al. 2022; Khoperskov et al. 2025a). Instead, the ability of younger, more metal-rich populations to be trapped by the bar more efficiently compared to the older metal-poor counterparts results in the more metal-rich MW bar, as is suggested by the kinematic fractionation mechanism (Debattista et al. 2017; Fragkoudi et al. 2017a; Khoperskov et al. 2018a) and in agreement with extragalactic observations (Neumann et al. 2024).

thumbnail Fig. 13

Chemical abundance profiles across the bulge region. The left and middle panels in the top row show the mean stellar mass-weighted [Fe/H] and [Mg/Fe], respectively, while the right one depicts the stellar mass fraction of the low-α populations (see Fig. 9 for the definition). The bottom panels show one-dimensional profiles of the same quantities as a function of longitude and at different latitudes, marked with different colours. The bottom group of panels shows the same, but for the vertical profiles, along latitude. The black lines at the bottom show the latitude- or longitude-averaged profiles. The grey background lines in the 2D abundance maps depict the stellar isodensity contours.
thumbnail Fig. 14

Maps of the radial (top) and vertical (bottom) metallicity gradients in the MW bulge region. The maps show the slopes of the metallicity profiles measured in each pixel within ±1 deg in the corresponding direction. The radial and vertical gradient variations trace well the X-shaped and boxy density components, respectively. The gradients are steeper near the transition regions between high- and low-α populations. The grey background lines reflect the projected stellar density distribution.

5.3 Chemo-kinematics of bulge populations: 2D GMM decomposition

In this section, we investigate the chemo-kinematic composition of the MW bulge populations, which are defined using the GMMs of the MDF as is widely done in the literature. In our case, thanks to the complete coverage of the bulge and the knowledge of the masses of different mono-abundance populations, instead of using the MDF only, we use the GMM approach to define the bulge populations in the two-dimensional [Mg/Fe]–[Fe/H] plane. This allows us to link the chemically defined populations of the bulge with the conventional chemical abundance structure of the MW disc(s).

To identify the centres of the GMM components in the [Mg/Fe]–[Fe/H] plane, we used the entire bulge sample (<3.5 kpc from the Galactic centre). Following previous studies (Bensby et al. 2017; Rojas-Arriagada et al. 2020; Wylie et al. 2021), we ran the 2D GMM search for randomly sampled 10% of the bulge populations. We repeated this procedure 50 times, and we adopted a full covariance GMM mode for each sample. We fitted the [Mg/Fe]–[Fe/H] density distributions with an increasing number of Gaussian components, from two to six, and analysed the output using the Bayesian information criteria (BIC, Schwarz 1978), which is minimal for five components, precisely the same number as in Ness et al. (2013a) or in Bensby et al. (2017). However, we underline that our GMM decomposition is different from the results of these studies. At this stage, we record only the positions ([Mg/Fe] and [Fe/H]) of the Gaussian components. The final centres of the 2D GMM components were obtained by averaging the positions obtained across all 50 random subsamples. Next, we ran the GMM decomposition with the fixed centres of the components in different latitude slabs, allowing the model to find the independent orientation of each component at different positions relative to the midplane.

The results of the 2D GMM decomposition of the [Mg/Fe]–[Fe/H] relation for the bulge are presented in Fig. 16. The lines of different colours show the individual GMM components in [Fe/H]–[Mg/Fe] plane (left) and stacked [Fe/H] (middle) and [Mg/Fe] (right) distributions. The fractional contributions of each component at different heights are marked with the same colour in the middle panels.

The most natural explanation of the most metal-poor component ([Fe/H] ≈ −0.8 dex, dark blue) is directly related to its location in the [Mg/Fe]-[Fe/H] plane, where chemically defined halo and its ex-situ populations are located (Helmi et al. 2018; Hayes et al. 2018; Di Matteo et al. 2019b; Mackereth et al. 2019). This is reinforced by its slightly lower [Mg/Fe] values compared to the high-α in-situ populations. The fractional contribution of the accreted populations increases from 5% in the centre to 12% at high latitudes. Interestingly, these values match very well various estimations of the mass of the stellar halo or spheroidal component in the MW bulge (Shen et al. 2010; Kunder et al. 2012; Di Matteo et al. 2014; Kunder et al. 2016; Debattista et al. 2017).

The following GMM components essentially decompose the [Fe/H]–[Mg/Fe] in-situ stars distribution into four populations. The light blue and red correspond to the bulk of high- and low-α populations, respectively, while the intermediate (cyan and green) two represent the tails of the bridge connecting the main components. Therefore, the physical meaning of the MW bulge GMM components is rather obvious. As demonstrated in the previous section, individual peaks of the MDF seen in our data do not correspond to specific orbital families, which have broad MDFs. Moreover, once we apply GMM decomposition, we do not recover these peaks, which seem to exist on top of the smooth distribution governed by the chemical evolution track(s) along metallicity. Essentially, the strength of the high- and low-α GMM components is defined by the masses of the thick and thin disc populations in the bulge region. The transition between these two populations, or the α-knee, shapes the other two intermediate GMM components. The fact that we observe a gap in the [Mg/Fe]–[Fe/H] relation around metallcity ≈−0.05 dex (see Fig. 9, right panel) may imply the star formation quenching episode (Haywood et al. 2016b, 2018a). The manifestation of the quenching in the [Mg/Fe]–[Fe/H] plane allows us to break down the bridge connecting the high- and low-α populations into two GMM components.

Although we do not suggest that the intermediate components (cyan and green in Fig. 16) should be associated with distinct bulge populations, they can have very different spatial and kinematic characteristics, especially if the MW experienced rapid transformation at the corresponding epoch. To illustrate this in Fig. 17, we present the face-on and (l, b) stellar density projections of populations identified by the 2D GMM, which are obtained by multiplying the stellar density by the GMM component probability distributions. The maps show the trends presented in Fig. 12. The most metal-poor component (first column) appears as an elongated structure along the bar’s major axis, with a weak flattening towards the midplane but lacks a prominent boxy or X-shaped morphology. It is worth noting that even when using the GMM probabilities as weights, some contamination from the high-α bulge population cannot be entirely ruled out because the membership for each component is based on the distance from the component centre. Hence, overlapping 2D distributions cannot be fully separated unless more dimensions are used, which is the case of accreted/in-situ population distinction (Feuillet et al. 2021; Buder et al. 2022; Myeong et al. 2022; Khoperskov et al. 2023b). Nevertheless, the structure of this component is distinct from the next one (second column), which represents the bulk of the high-α populations. While the X-shaped structure is noticeable, it is not particularly sharp, giving the bar an overall oval-like appearance. In the third GMM component, the bar increases in size. Surprisingly, the longest bar is observed in the fourth GMM component, representing the low-[Fe/H] tail of the low-α populations. This likely corresponds to the coldest stellar populations most tightly trapped by the bar. The bar remains relatively long for the most metal-rich component, though it is somewhat less extended than the previous component because of the rapidly decreasing scalelength of mono-abundance populations at [Fe/H] > +0.25 (Khoperskov et al. 2025a, Rix et al. 2024).

The kinematics of the four individual in-situ components follow the trend we discussed above (see Section. 4): increasing rotational support (and decreasing velocity dispersion) with increasing metallicity (see a relevant figure in Appendix A.2). The most metal-poor, potentially accreted GMM component still shows very weak rotation, which, as we already mentioned, can be contaminated by the metal-poor tail of the high-α populations.

Since the release of the second Gaia catalogue, which provided phase-space data for millions of stars in the MW, numerous studies have sought to identify whether various stellar populations have in-situ or accreted origin using kinematic spaces (see Malhan et al. 2022, for compilation of the results). One of the most commonly used methods for identifying accreted systems is through the energy-angular momentum space (E − Lz), as these quantities are generally assumed to be conserved after a smaller galactic system falls into the MW (Helmi & de Zeeuw 2000). While this assumption is applicable for low-mass accreted systems, it is not entirely correct for massive galaxies, whose remnants are affected by dynamical friction (Jean-Baptiste et al. 2017; Pagnini et al. 2023; Mori et al. 2024), and because the MW mass distribution (stellar and DM) is non-axisymmetric and evolves over time (Panithanpaisal et al. 2021; Khoperskov et al. 2023a). Although our orbit superposition approach is primarily focused on the main disc stellar populations and was not explicitly designed to accommodate accreted systems, it is worth understanding the structure of the in-situ component in E − Lz space. This provides a reference framework, helping to distinguish accreted systems from genuine MW stars in future, especially in the innermost region of the bulge, where the remnants of anciently accreted systems should be buried.

In Figure 18, we present the E − Lz space for populations identified using a 2D GMM. The top row illustrates the GMM probability-weighted mass distributions in E − Lz space, revealing a double-component structure across all populations, regardless of their metallicity. However, the specific details of these distributions vary slightly. To better understand the origin of this double-peak distribution, we separated the populations into two groups: stars located along the bar (within ±1 kpc of the bar major axis) and those outside the bar (within <3.5 kpc but not in the bar region). The bottom two rows in the figure demonstrate that bimodal E − Lz distribution clearly arises from the superposition of bar and disc components. As expected, the bar component exhibits minimal rotation, as stars in this region tend to follow radial orbits, while the disc populations display weak but prograde rotation. Here, we underline that the non-rotating component, Lz ≈ 0, is not a signature of a spheroidal component but the bar.

This bimodal angular momentum distribution is not surprising and was also observed in an N-body simulation of a MW-like galaxy by Fragkoudi et al. (2018). In their findings, thin disc stars were shown to lose more angular momentum due to bar/bulge formation, resulting in a more pronounced low-angular-momentum peak, whereas thick disc stars were less efficiently trapped by the bar and thus exhibited less pronounced bimodality in Lz. This trend is consistent with what we observe in Fig. 18, where the fraction of low-Lz stars increases with metallicity for each GMM-identified population, as indicated by the values displayed in each panel.

In most metal-poor E − Lz bin, we observe several distinct features or clumps, which could potentially be remnants of yet unknown accreted galaxies, bulge/bar resonances (Dillamore et al. 2024b,a), or stars stripped from globular clusters in the bulge region (Ferrone et al. 2023). Determining the true origin of these features requires further detailed investigation, and we do not assert their definitive nature at this stage. However, the fact that these features appear exclusively among the most metalpoor populations is in favour of the hypothesis that they are of ex-situ origin or related to old (not necessarily accreted) globular clusters.

thumbnail Fig. 15

Face-on maps of the mean stellar metallicity [Fe/H] distribution across the MW disc. From left to right, the maps represent stellar populations within horizontal slabs at increasing distances from the midplane. Near the midplane (two panels on left), the bulge region appears to be metal-poor because of the deficit of the metal-rich (low-α) populations expelled to higher distances perpendicular to the disc plane via bar/bulge instability. However, the actual metallicity of the MW bulge is higher compared to the surrounding disc once its full vertical extension is taken into account. The bar appears to be metal-rich along the major axis because kinematically colder low-α (metal-rich) populations are more efficiently captured by the bar’s gravitational potential compared to the hotter high-α (metal-poor) stars dominating along the bar minor axis.
thumbnail Fig. 16

Two-dimensional decomposition of the MW bulge in the [Mg/Fe]–[Fe/H] plane as a function of latitude based on the data inside 3.5 kpc from the Galactic centre. The left panels show the normalised stellar density maps, where the contours of different colours show the components identified using the 2D GMM approach. The grey histograms in the middle and right columns show the generalised [Fe/H]– and [Mg/Fe]–distribution functions, respectively. The coloured lines show the one-dimensional projections of the 2D-GMM components, while the black lines represent the sum of these components. In the middle row, the fractional contribution of each component is marked by the numbers of the same colour. Our analysis suggests that the [Mg/Fe]–[Fe/H] bulge composition is best represented by the five GMM components, where only two are the main ones and correspond to the bulk of high- and low-α populations. Two intermediate components allow for the decomposition of the transition between the main ones by fitting the α-knee and likely do not correspond to distinct bulge populations. The most metal-poor GMM component likely corresponds to the accreted populations.
thumbnail Fig. 17

Density structure of the MW bulge components identified using 2D GMM in the [Fe/H]–[Mg/Fe] plane. The top panels show the stellar mass density of each component obtained by multiplying the total stellar mass in the plane (background contours) by the 2D GMM probability for a given component. The middle and bottom rows show the mass×GMM-probability weighed density distribution in the XY and (l, b) coordinates, respectively. The density is scaled by the maximum value in each panel.

6 Summary

Nowadays, the MW bulge is seen as the inner structure of the bar and is certainly a complex environment. It is a composite system with a mixture of different stellar populations, similar to other nearby extragalactic systems. Studying the bulge, which is a major component of our Galaxy, is an important step in clarifying how the MW formed and evolved. In this work, we have applied our novel orbit superposition approach combined with the APOGEE DR 17 data to provide a unified present-day picture of the correlations between the spatial structure, chemical abundances, kinematics, and orbital composition of stellar populations in the MW bulge. The main results of the work are summarised as follows:

  • We have shown that our orbit superposition method is precise enough, and the APOGEE giant stars sample provides sufficient data to successfully reconstruct the 3D stellar density structure of the MW bulge and bar region. This includes capturing the distinct X-shaped/boxy structure and the bimodal peanut density distribution observed along different lines of sight towards the bulge, further validating the effectiveness of our approach in recovering the complex structure of the inner MW;

  • The fractional contributions of various orbital families that support the X-shaped boxy/peanut structure of the MW bulge highlight its intricate orbital dynamics, as predicted in secular formation scenarios. Our analysis reveals that banana-like orbits play a substantial role, contributing approximately 30% to the overall bulge and bar structure. Meanwhile, the inner X-shaped orbits account for a more prominent 45%, underscoring their critical importance in shaping the distinctive morphology of the inner bulge;

  • In the absence of any effects from the APOGEE footprint, we analysed the complete kinematics of the bulge. We are able to detect an increase (decrease) in the cylindrical rotation (velocity dispersion) rate with increasing metallicity in the range from −1.5 to +0.65. Thanks to the high spatial coverage of the bulge, we detect small-scale asymmetries introduced by the bar, the major axis of which lies at an angle to the Sun–Galactic centre line. Additionally, we observed flattening in the 1D velocity dispersion profiles, characterised by shoulder-like features at various latitudes, which correspond to the X-shaped component of the bulge;

    thumbnail Fig. 18

    Energy-angular momentum structure of the MW bulge components identified using 2D GMM in the [Fe/H]–[Mg/Fe] plane. The top panels are the same as in Fig. 17 and are given for reference. The bottom group of panels show the normalised stellar mass for the entire component (top) and the region inside (middle) and outside (bottom) the bar. The numbers in the middle and bottom rows correspond to the mass fraction of a given component inside and outside the bar. The white dashed lines mark the minimum possible energy for a given value of the angular momentum.

  • Chemical abundance profiles vary greatly in the bulge region along longitude and latitude. The radial metallicity profiles across the bulge region are not monotonic, while the vertical ones are negative but change with longitude. Close to the midplane (|b| < 2°), the mean metallicity decreases towards the Galactic centre, suggesting the lower metallicity of the MW bulge compared to the surrounding disc (the metal-poor bulge in the (l, b) projection). This picture is already inverted at latitudes b > 2°. Once we measure the mean metallicity, including high distances from the midplane, we find that the inner bulge region and the bar major axis are more metal-rich than the surrounding disc. Therefore, while the mean metallicity of the bulge is slightly subsolar, we argue that it is metal-rich, which is in perfect agreement with nearby boxy/peanut barred galaxies and simulations;

  • Overall, there is no universal metallicity gradient value that can characterise the MW bulge. The radial gradients very closely trace the X-shaped bulge density structure, reaching the maximum of 0.01 dex deg−1 near the tips of the X-shaped component. The spatial variations in the vertical gradient better trace the boxy component, reaching the maximum values (−0.07 dex deg−1) at high latitudes, on the edge of the bulge, where the transition from a low- to high-α dominated region is very prominent. The vertical metallicity gradient reaches the smallest values in the innermost region (|l| < 3° and |b| < 2°), where the metallicity profile is flat;

  • We have shown that the MW bulge is composed of two primary components: high- and low-α populations. These correspond to the inner thick disc, which is metal-poor with subsolar metallicity, and the thin disc, which is metal-rich with a supersolar metallicity, respectively. The mass ratio between these chemically distinct populations is 4:3, indicating the slight predominance of the metal-poor component over the metal-rich one;

  • The MDF of the MW bulge is broad, featuring two dominant peaks at approximately −0.7 dex and +0.3 dex, which correspond to the maximum density of the high- and low-α populations. Between these two main components, we observe multiple narrow peaks, which are present in the orbital families that support the bulge but are absent in the MDF of the background disc populations. However, we do not propose a direct one-to-one connection between individual MDF peaks and specific orbital families, as the latter exhibit a range of metallicities and have relatively broad MDFs themselves. The latter illustrates that stars formed at different epochs can be trapped by the X-shaped/boxy structure with varying efficiency, depending on the initial distribution function of these populations and the evolution of the bar/bulge parameters;

  • Although the MW bulge is primarily dominated by two main components, its broad MDF cannot be fully explained by these components alone. To delve into the complexity of the bulge chemical abundance composition, we decomposed it using the 2D GMM approach in the stellar-mass weighted [Fe/H]-[Mg/Fe] plane. We find that the optimal number of 2D GMM components is five, where the most metal-poor one is dominated by stars with a chemical composition typical of accreted populations. Although our orbit superposition approach is not tailored to include distinct halo populations, we estimate that the mass of stars with an exsitu chemical composition – those accreted from outside the MW – amounts to approximately 8 × 108M⊙, for stars with a metallicity above ≈−1.2 dex. Between the primary high- and low-α populations, we identified two additional 2D GMM components. These metal-intermediate components provide a more detailed breakdown of the α-knee and help to characterise the ‘bridge’ that connects the high- and low-α populations. We therefore suggest that these intermediate MDF components do not correspond to distinct bulge populations but rather describe the transition between the two main chemical components.

In this study, we provide a comprehensive view of the present-day MW bulge. By analysing its chemo-kinematic properties, we confirm that the bulge is a product of vertical bar instability, which captured and reshaped the inner populations of both the thin and thick discs over time. Aside from a relatively small fraction of stars from accreted populations, we find no significant evidence of a spheroidal component in the bulge. This suggests that the MW bulge shares many characteristics with secularly formed bulges in other barred galaxies, reinforcing the idea that the MW is a typical barred galaxy.

A few pieces still remain unresolved in the MW bugle story, particularly regarding the precise mechanism behind the bulge formation. It is unclear whether it was shaped by vertically asymmetric buckling (Pfenniger & Friedli 1991; Merritt & Sellwood 1994) or a slower, more symmetric resonant heating process relative to the midplane (Quillen 2002; Quillen et al. 2014). Additionally, the exact timing of the X-shaped structure’s formation remains unconstrained. Although N-body simulations indicate that vertical instability often occurs relatively quickly after the bar’s development, no clear chemo-kinematic indicators that allow us to constrain the epoch of bulge formation in the MW have yet been identified in models (see the discussion in Debattista et al. 2019). Further theoretical efforts are required to address these open questions. On the observational side, more data will be provided in the inner parts of the MW by forthcoming spectroscopic surveys like 4MOST (de Jong et al. 2019), SDSS-V (Almeida et al. 2023), and MOONS (Gonzalez et al. 2020). Additionally, the continued application of the orbit super-position approach and, more generally, Schwarzschild methods will be crucial in advancing our understanding of the MW bulge and its formation history.

Data availability

Movies associated to Figs. 3 & 12 are available at https://www.aanda.org

Acknowledgements

We thank the anonymous referee for their constructive comments. SK expresses gratitude to Ortwin Gerhard for an insightful discussion of relevant topics in April 2024. This research was supported by the Munich Institute for Astro-, Particle and BioPhysics (MIAPbP) which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany´s Excellence Strategy-EXC-2094-390783311 The authors thank the organizers of the “A new dawn of dwarf galaxy research” workshop and the Lorentz Center (Leiden) for the organizational and financial support of the meeting, which provided an insightful atmosphere which stimulated the progress of this work. This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia Multi-Lateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia Archive website is http://archives.esac.esa.int/gaia. 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 acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss4.org. SDSS 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 (CfA), 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.

Appendix A Extra figures

thumbnail Fig. A.1

Stellar mass-weighted age distribution (top) and age-metallicity relation (bottom) in the MW bulge region (< 3.5 kpc from the Galactic centre). The metallicity range ([Fe/H] < −0.65 dex, Stone-Martinez et al. (2024)) with uncertain ages is shaded in the bottom panel and not included in the top distribution, accounting for about 13% of the bulge mass.
thumbnail Fig. A.2

Kinematics of the MW bulge populations identified using 2D GMM in the [Fe/H]-[Mg/Fe] plane. The top panels are the same as in Fig. 17 and are given for reference. The middle and bottom panels show each component’s line-of-sight velocity and dispersion maps, weighted by the mass multiplied by the GMM probabilities. The grey contours show the corresponding density distribution, while the white ones highlight the velocity/dispersion structure.

References

  1. Abdurro’uf, Accetta, K., Aerts, C., et al. 2022, ApJS, 259, 35 [NASA ADS] [CrossRef] [Google Scholar]
  2. Aguerri, J. A. L., Méndez-Abreu, J., Falcón-Barroso, J., et al. 2015, A&A, 576, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Almeida, A., Anderson, S. F., Argudo-Fernández, M., et al. 2023, ApJS, 267, 44 [NASA ADS] [CrossRef] [Google Scholar]
  4. Alves-Brito, A., Meléndez, J., Asplund, M., Ramírez, I., & Yong, D. 2010, A&A, 513, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Anderson, S. R., Gough-Kelly, S., Debattista, V. P., et al. 2024, MNRAS, 527, 2919 [Google Scholar]
  6. Antoja, T., Helmi, A., Romero-Gómez, M., et al. 2018, Nature, 561, 360 [Google Scholar]
  7. Arentsen, A., Starkenburg, E., Martin, N. F., et al. 2020, MNRAS, 491, L11 [Google Scholar]
  8. Athanassoula, E. 2003, MNRAS, 341, 1179 [Google Scholar]
  9. Athanassoula, E. 2008, in IAU Symposium, 245, Formation and Evolution of Galaxy Bulges, eds. M. Bureau, E. Athanassoula, & B. Barbuy, 93 [Google Scholar]
  10. Athanassoula, E., & Misiriotis, A. 2002, MNRAS, 330, 35 [Google Scholar]
  11. Babusiaux, C., & Gilmore, G. 2005, MNRAS, 358, 1309 [NASA ADS] [CrossRef] [Google Scholar]
  12. Babusiaux, C., Gómez, A., Hill, V., et al. 2010, A&A, 519, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Barbuy, B., Chiappini, C., & Gerhard, O. 2018, ARA&A, 56, 223 [Google Scholar]
  14. Barbuy, B., Friaça, A. C. S., Ernandes, H., et al. 2023, MNRAS, 526, 2365 [NASA ADS] [CrossRef] [Google Scholar]
  15. Bekki, K., & Tsujimoto, T. 2011a, MNRAS, 416, L60 [NASA ADS] [CrossRef] [Google Scholar]
  16. Bekki, K., & Tsujimoto, T. 2011b, ApJ, 738, 4 [Google Scholar]
  17. Belokurov, V., & Kravtsov, A. 2022, MNRAS, 514, 689 [NASA ADS] [CrossRef] [Google Scholar]
  18. Belokurov, V., & Kravtsov, A. 2023, MNRAS, 525, 4456 [NASA ADS] [CrossRef] [Google Scholar]
  19. Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, MNRAS, 478, 611 [Google Scholar]
  20. Belokurov, V., Sanders, J. L., Fattahi, A., et al. 2020, MNRAS, 494, 3880 [Google Scholar]
  21. Bensby, T., Feltzing, S., & Lundström, I. 2003, A&A, 410, 527 [CrossRef] [EDP Sciences] [Google Scholar]
  22. Bensby, T., Feltzing, S., Johnson, J. A., et al. 2010, A&A, 512, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Bensby, T., Adén, D., Meléndez, J., et al. 2011, A&A, 533, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Bensby, T., Yee, J. C., Feltzing, S., et al. 2013, A&A, 549, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Bensby, T., Feltzing, S., & Oey, M. S. 2014, A&A, 562, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Bensby, T., Feltzing, S., Gould, A., et al. 2017, A&A, 605, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Bensby, T., Bergemann, M., Rybizki, J., et al. 2019, The Messenger, 175, 35 [NASA ADS] [Google Scholar]
  28. Bensby, T., Feltzing, S., Yee, J. C., et al. 2020, A&A, 634, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Bensby, T., Gould, A., Asplund, M., et al. 2021, A&A, 655, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Beraldo e Silva, L., Debattista, V. P., Anderson, S. R., et al. 2023, ApJ, 955, 38 [NASA ADS] [CrossRef] [Google Scholar]
  31. Boin, T., Di Matteo, P., Khoperskov, S., et al. 2024, A&A, 692, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Bovy, J., Rix, H.-W., Liu, C., et al. 2012, ApJ, 753, 148 [Google Scholar]
  33. Bovy, J., Leung, H. W., Hunt, J. A. S., et al. 2019, MNRAS, 490, 4740 [Google Scholar]
  34. Breda, I., & Papaderos, P. 2018, A&A, 614, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Buck, T., Ness, M. K., Macciò, A. V., Obreja, A., & Dutton, A. A. 2018, ApJ, 861, 88 [Google Scholar]
  36. Buck, T., Ness, M., Obreja, A., Macciò, A. V., & Dutton, A. A. 2019, ApJ, 874, 67 [NASA ADS] [CrossRef] [Google Scholar]
  37. Buder, S., Lind, K., Ness, M. K., et al. 2022, MNRAS, 510, 2407 [NASA ADS] [CrossRef] [Google Scholar]
  38. Bureau, M., & Freeman, K. C. 1999, AJ, 118, 126 [NASA ADS] [CrossRef] [Google Scholar]
  39. Cabrera-Lavers, A., Hammersley, P. L., González-Fernández, C., et al. 2007, A&A, 465, 825 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Chiba, R., Friske, J. K. S., & Schönrich, R. 2021, MNRAS, 500, 4710 [Google Scholar]
  41. Ciambur, B. C., Fragkoudi, F., Khoperskov, S., Di Matteo, P., & Combes, F. 2021, MNRAS, 503, 2203 [NASA ADS] [CrossRef] [Google Scholar]
  42. Cirasuolo, M., Fairley, A., Rees, P., et al. 2020, The Messenger, 180, 10 [NASA ADS] [Google Scholar]
  43. Clarke, J. P., & Gerhard, O. 2022, MNRAS, 512, 2171 [NASA ADS] [CrossRef] [Google Scholar]
  44. Clarke, J. P., Wegg, C., Gerhard, O., et al. 2019, MNRAS, 489, 3519 [Google Scholar]
  45. Combes, F., & Sanders, R. H. 1981, A&A, 96, 164 [NASA ADS] [Google Scholar]
  46. de Jong, R. S., Agertz, O., Berbel, A. A., et al. 2019, The Messenger, 175, 3 [NASA ADS] [Google Scholar]
  47. Debattista, V. P., & Sellwood, J. A. 1998, ApJ, 493, L5 [Google Scholar]
  48. Debattista, V. P., & Sellwood, J. A. 2000, ApJ, 543, 704 [Google Scholar]
  49. Debattista, V. P., Carollo, C. M., Mayer, L., & Moore, B. 2004, ApJ, 604, L93 [NASA ADS] [CrossRef] [Google Scholar]
  50. Debattista, V. P., Carollo, C. M., Mayer, L., & Moore, B. 2005, ApJ, 628, 678 [Google Scholar]
  51. Debattista, V. P., Mayer, L., Carollo, C. M., et al. 2006, ApJ, 645, 209 [Google Scholar]
  52. Debattista, V. P., Ness, M., Earp, S. W. F., & Cole, D. R. 2015, ApJ, 812, L16 [Google Scholar]
  53. Debattista, V. P., Ness, M., Gonzalez, O. A., et al. 2017, MNRAS, 469, 1587 [Google Scholar]
  54. Debattista, V. P., Gonzalez, O. A., Sanderson, R. E., et al. 2019, MNRAS, 485, 5073 [Google Scholar]
  55. Debattista, V. P., Liddicott, D. J., Khachaturyants, T., & Beraldo e Silva, L. 2020, MNRAS, 498, 3334 [Google Scholar]
  56. Debattista, V. P., Liddicott, D. J., Gonzalez, O. A., et al. 2023, ApJ, 946, 118 [NASA ADS] [CrossRef] [Google Scholar]
  57. Dékány, I., Minniti, D., Catelan, M., et al. 2013, ApJ, 776, L19 [Google Scholar]
  58. Di Matteo, P. 2016, PASA, 33, e027 [CrossRef] [Google Scholar]
  59. Di Matteo, P., Haywood, M., Gómez, A., et al. 2014, A&A, 567, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Di Matteo, P., Gómez, A., Haywood, M., et al. 2015, A&A, 577, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Di Matteo, P., Fragkoudi, F., Khoperskov, S., et al. 2019a, A&A, 628, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Di Matteo, P., Haywood, M., Lehnert, M. D., et al. 2019b, A&A, 632, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Dillamore, A. M., Belokurov, V., & Evans, N. W. 2024a, MNRAS, 532, 4389 [NASA ADS] [CrossRef] [Google Scholar]
  64. Dillamore, A. M., Monty, S., Belokurov, V., & Evans, N. W. 2024b, ApJ, 971, L4 [Google Scholar]
  65. Dwek, E., Arendt, R. G., Hauser, M. G., et al. 1995, ApJ, 445, 716 [CrossRef] [Google Scholar]
  66. Efstathiou, G., Lake, G., & Negroponte, J. 1982, MNRAS, 199, 1069 [NASA ADS] [CrossRef] [Google Scholar]
  67. Elia, D., Molinari, S., Schisano, E., et al. 2022, ApJ, 941, 162 [NASA ADS] [CrossRef] [Google Scholar]
  68. Erwin, P., & Debattista, V. P. 2017, MNRAS, 468, 2058 [Google Scholar]
  69. Feltzing, S., & Gilmore, G. 2000, A&A, 355, 949 [NASA ADS] [Google Scholar]
  70. Ferrone, S., Di Matteo, P., Mastrobuono-Battisti, A., et al. 2023, A&A, 673, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Feuillet, D. K., Sahlholdt, C. L., Feltzing, S., & Casagrande, L. 2021, MNRAS, 508, 1489 [NASA ADS] [CrossRef] [Google Scholar]
  72. Fragkoudi, F., Di Matteo, P., Haywood, M., et al. 2017a, A&A, 606, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Fragkoudi, F., Di Matteo, P., Haywood, M., et al. 2017b, A&A, 607, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Fragkoudi, F., Di Matteo, P., Haywood, M., et al. 2018, A&A, 616, A180 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Fragkoudi, F., Katz, D., Trick, W., et al. 2019, MNRAS, 488, 3324 [NASA ADS] [Google Scholar]
  76. Fragkoudi, F., Grand, R. J. J., Pakmor, R., et al. 2020, MNRAS, 494, 5936 [Google Scholar]
  77. Fragkoudi, F., Grand, R. J. J., Pakmor, R., et al. 2021, A&A, 650, L16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Fraser-McKelvie, A., Merrifield, M., Aragón-Salamanca, A., et al. 2020, MNRAS, 499, 1116 [CrossRef] [Google Scholar]
  79. Fulbright, J. P., McWilliam, A., & Rich, R. M. 2007, ApJ, 661, 1152 [NASA ADS] [CrossRef] [Google Scholar]
  80. Fux, R. 1997, A&A, 327, 983 [NASA ADS] [Google Scholar]
  81. Gadotti, D. A. 2011, MNRAS, 415, 3308 [Google Scholar]
  82. Gaia Collaboration (Katz, D., et al.) 2018, A&A, 616, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  83. Gaia Collaboration (Antoja, T., et al.) 2021, A&A, 649, A8 [EDP Sciences] [Google Scholar]
  84. Gaia Collaboration (Drimmel, R., et al.) 2023a, A&A, 674, A37 [CrossRef] [EDP Sciences] [Google Scholar]
  85. Gaia Collaboration (Vallenari, A., et al.) 2023b, A&A, 674, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. George, K., Subramanian, S., & Paul, K. T. 2019, A&A, 628, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Gerhard, O. 2020, in IAU Symposium, 353, Galactic Dynamics in the Era of Large Surveys, eds. M. Valluri, & J. A. Sellwood, 26 [Google Scholar]
  88. Géron, T., Smethurst, R. J., Lintott, C., et al. 2021, MNRAS, 507, 4389 [CrossRef] [Google Scholar]
  89. Gómez, A., Di Matteo, P., Stefanovitch, N., et al. 2016, A&A, 589, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Gonzalez, O. A., Rejkuba, M., Zoccali, M., et al. 2011, A&A, 530, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Gonzalez, O. A., Rejkuba, M., Zoccali, M., et al. 2013, A&A, 552, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Gonzalez, O. A., Zoccali, M., Debattista, V. P., et al. 2015a, A&A, 583, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Gonzalez, O. A., Zoccali, M., Vasquez, S., et al. 2015b, A&A, 584, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Gonzalez, O. A., Gadotti, D. A., Debattista, V. P., et al. 2016, A&A, 591, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Gonzalez, O. A., Debattista, V. P., Ness, M., Erwin, P., & Gadotti, D. A. 2017, MNRAS, 466, L93 [Google Scholar]
  96. Gonzalez, O. A., Mucciarelli, A., Origlia, L., et al. 2020, The Messenger, 180, 18 [NASA ADS] [Google Scholar]
  97. Guo, R., Mao, S., Athanassoula, E., et al. 2019, MNRAS, 482, 1733 [Google Scholar]
  98. Hayes, C. R., Majewski, S. R., Shetrone, M., et al. 2018, ApJ, 852, 49 [NASA ADS] [CrossRef] [Google Scholar]
  99. Haywood, M., Di Matteo, P., Lehnert, M. D., Katz, D., & Gómez, A. 2013, A&A, 560, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  100. Haywood, M., Di Matteo, P., Snaith, O., & Calamida, A. 2016a, A&A, 593, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  101. Haywood, M., Lehnert, M. D., Di Matteo, P., et al. 2016b, A&A, 589, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  102. Haywood, M., Di Matteo, P., Lehnert, M., et al. 2018a, A&A, 618, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  103. Haywood, M., Di Matteo, P., Lehnert, M. D., et al. 2018b, ApJ, 863, 113 [Google Scholar]
  104. Haywood, M., Khoperskov, S., Cerqui, V., et al. 2024, A&A, 690, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  105. Helmi, A., & de Zeeuw, P. T. 2000, MNRAS, 319, 657 [Google Scholar]
  106. Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Nature, 563, 85 [Google Scholar]
  107. Hill, V., Lecureur, A., Gómez, A., et al. 2011, A&A, 534, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  108. Hilmi, T., Minchev, I., Buck, T., et al. 2020, MNRAS, 497, 933 [Google Scholar]
  109. Holmberg, J., Nordström, B., & Andersen, J. 2009, A&A, 501, 941 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  110. Howard, C. D., Rich, R. M., Clarkson, W., et al. 2009, ApJ, 702, L153 [NASA ADS] [CrossRef] [Google Scholar]
  111. Howard, C. D., Rich, R. M., Reitzel, D. B., et al. 2008, ApJ, 688, 1060 [NASA ADS] [CrossRef] [Google Scholar]
  112. Hunt, J. A. S., & Bovy, J. 2018, MNRAS, 477, 3945 [NASA ADS] [CrossRef] [Google Scholar]
  113. Imig, J., Price, C., Holtzman, J. A., et al. 2023, ApJ, 954, 124 [CrossRef] [Google Scholar]
  114. Jean-Baptiste, I., Di Matteo, P., Haywood, M., et al. 2017, A&A, 604, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  115. Johnson, C. I., Rich, R. M., Kobayashi, C., & Fulbright, J. P. 2012, ApJ, 749, 175 [NASA ADS] [CrossRef] [Google Scholar]
  116. Johnson, C. I., Rich, R. M., Kobayashi, C., Kunder, A., & Koch, A. 2014, AJ, 148, 67 [Google Scholar]
  117. Johnson, C. I., Rich, R. M., Simion, I. T., et al. 2022, MNRAS, 515, 1469 [NASA ADS] [CrossRef] [Google Scholar]
  118. Kacharov, N., Tahmasebzadeh, B., Cioni, M.-R. L., et al. 2024, A&A, 692, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  119. Khoperskov, S., & Gerhard, O. 2022, A&A, 663, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  120. Khoperskov, S., Di Matteo, P., Haywood, M., & Combes, F. 2018a, A&A, 611, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  121. Khoperskov, S., Haywood, M., Di Matteo, P., Lehnert, M. D., & Combes, F. 2018b, A&A, 609, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  122. Khoperskov, S., Gerhard, O., Di Matteo, P., et al. 2020, A&A, 634, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  123. Khoperskov, S., Minchev, I., Libeskind, N., et al. 2023a, A&A, 677, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  124. Khoperskov, S., Minchev, I., Steinmetz, M., et al. 2023b, arXiv e-prints, [arXiv:2310.05287] [Google Scholar]
  125. Khoperskov, S., Steinmetz, M., Haywood, M., et al. 2025a, A&A, 700, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  126. Khoperskov, S., van de Ven, G., Steinmetz, M., et al. 2025b, A&A, 695, A220 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  127. Kormendy, J., & Kennicutt, Robert C., J. 2004, ARA&A, 42, 603 [NASA ADS] [CrossRef] [Google Scholar]
  128. Kunder, A., Koch, A., Rich, R. M., et al. 2012, AJ, 143, 57 [NASA ADS] [CrossRef] [Google Scholar]
  129. Kunder, A., Rich, R. M., Koch, A., et al. 2016, ApJ, 821, L25 [CrossRef] [Google Scholar]
  130. Kunder, A., Pérez-Villegas, A., Rich, R. M., et al. 2020, AJ, 159, 270 [NASA ADS] [CrossRef] [Google Scholar]
  131. Leung, H. W., & Bovy, J. 2019, MNRAS, 489, 2079 [CrossRef] [Google Scholar]
  132. Leung, H. W., Bovy, J., Mackereth, J. T., & Miglio, A. 2023, MNRAS, 522, 4577 [NASA ADS] [CrossRef] [Google Scholar]
  133. Li, Z.-Y., & Shen, J. 2012, ApJ, 757, L7 [NASA ADS] [CrossRef] [Google Scholar]
  134. Lian, J., Bergemann, M., Pillepich, A., Zasowski, G., & Lane, R. R. 2023, Nat. Astron., 7, 951 [CrossRef] [Google Scholar]
  135. Łokas, E. L. 2019, A&A, 629, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  136. Mackereth, J. T., Schiavon, R. P., Pfeffer, J., et al. 2019, MNRAS, 482, 3426 [Google Scholar]
  137. Malhan, K., Ibata, R. A., Sharma, S., et al. 2022, ApJ, 926, 107 [NASA ADS] [CrossRef] [Google Scholar]
  138. Martinez-Valpuesta, I., & Gerhard, O. 2011, ApJ, 734, L20 [NASA ADS] [CrossRef] [Google Scholar]
  139. Martinez-Valpuesta, I., & Gerhard, O. 2013, ApJ, 766, L3 [NASA ADS] [CrossRef] [Google Scholar]
  140. Martinez-Valpuesta, I., Shlosman, I., & Heller, C. 2006, ApJ, 637, 214 [NASA ADS] [CrossRef] [Google Scholar]
  141. McWilliam, A. 2016, PASA, 33, e040 [NASA ADS] [CrossRef] [Google Scholar]
  142. McWilliam, A., & Zoccali, M. 2010, ApJ, 724, 1491 [NASA ADS] [CrossRef] [Google Scholar]
  143. Merritt, D., & Sellwood, J. A. 1994, ApJ, 425, 551 [NASA ADS] [CrossRef] [Google Scholar]
  144. Minchev, I., Chiappini, C., & Martig, M. 2013, A&A, 558, A9 [CrossRef] [EDP Sciences] [Google Scholar]
  145. Minniti, D., Olszewski, E. W., Liebert, J., et al. 1995, MNRAS, 277, 1293 [NASA ADS] [CrossRef] [Google Scholar]
  146. Minniti, D., Lucas, P. W., Emerson, J. P., et al. 2010, New A, 15, 433 [Google Scholar]
  147. Molaeinezhad, A., Falcón-Barroso, J., Martínez-Valpuesta, I., et al. 2016, MNRAS, 456, 692 [NASA ADS] [CrossRef] [Google Scholar]
  148. Moorthy, B. K., & Holtzman, J. A. 2006, MNRAS, 371, 583 [Google Scholar]
  149. Mori, A., Di Matteo, P., Salvadori, S., et al. 2024, A&A, 690, A136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  150. Myeong, G. C., Belokurov, V., Aguado, D. S., et al. 2022, ApJ, 938, 21 [NASA ADS] [CrossRef] [Google Scholar]
  151. Nandakumar, G., Schultheis, M., Hayden, M., et al. 2017, A&A, 606, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  152. Nandakumar, G., Ryde, N., Mace, G., et al. 2024, ApJ, 964, 96 [NASA ADS] [CrossRef] [Google Scholar]
  153. Nataf, D. M. 2017, PASA, 34, e041 [NASA ADS] [CrossRef] [Google Scholar]
  154. Nataf, D. M., Udalski, A., Gould, A., Fouqué, P., & Stanek, K. Z. 2010, ApJ, 721, L28 [NASA ADS] [CrossRef] [Google Scholar]
  155. Ness, M., & Lang, D. 2016, AJ, 152, 14 [NASA ADS] [CrossRef] [Google Scholar]
  156. Ness, M., & Freeman, K. 2016, PASA, 33, e022 [Google Scholar]
  157. Ness, M., Freeman, K., Athanassoula, E., et al. 2012, ApJ, 756, 22 [NASA ADS] [CrossRef] [Google Scholar]
  158. Ness, M., Freeman, K., Athanassoula, E., et al. 2013a, MNRAS, 430, 836 [NASA ADS] [CrossRef] [Google Scholar]
  159. Ness, M., Debattista, V. P., Bensby, T., et al. 2014, ApJ, 787, L19 [NASA ADS] [CrossRef] [Google Scholar]
  160. Ness, M., Freeman, K., Athanassoula, E., et al. 2013b, MNRAS, 432, 2092 [NASA ADS] [CrossRef] [Google Scholar]
  161. Ness, M., Zasowski, G., Johnson, J. A., et al. 2016, ApJ, 819, 2 [NASA ADS] [CrossRef] [Google Scholar]
  162. Neumann, J., Fragkoudi, F., Pérez, I., et al. 2020, A&A, 637, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  163. Neumann, J., Thomas, D., Maraston, C., et al. 2024, MNRAS, 534, 2438 [NASA ADS] [Google Scholar]
  164. Nordström, B., Mayor, M., Andersen, J., et al. 2004, A&A, 418, 989 [Google Scholar]
  165. Okuda, H., Maihara, T., Oda, N., & Sugiyama, T. 1977, Nature, 265, 515 [NASA ADS] [CrossRef] [Google Scholar]
  166. Ortolani, S., Renzini, A., Gilmozzi, R., et al. 1995, Nature, 377, 701 [CrossRef] [Google Scholar]
  167. Ostriker, J. P., & Peebles, P. J. E. 1973, ApJ, 186, 467 [Google Scholar]
  168. Pagnini, G., Di Matteo, P., Khoperskov, S., et al. 2023, A&A, 673, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  169. Panithanpaisal, N., Sanderson, R. E., Wetzel, A., et al. 2021, ApJ, 920, 10 [NASA ADS] [CrossRef] [Google Scholar]
  170. Parul, H. D., Smirnov, A. A., & Sotnikova, N. Y. 2020, ApJ, 895, 12 [Google Scholar]
  171. Patsis, P. A., & Katsanikas, M. 2014, MNRAS, 445, 3525 [NASA ADS] [CrossRef] [Google Scholar]
  172. Patsis, P. A., & Harsoula, M. 2018, A&A, 612, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  173. Patsis, P. A., & Athanassoula, E. 2019, MNRAS, 490, 2740 [Google Scholar]
  174. Patsis, P. A., Skokos, C., & Athanassoula, E. 2002, MNRAS, 337, 578 [NASA ADS] [CrossRef] [Google Scholar]
  175. Pfenniger, D., & Friedli, D. 1991, A&A, 252, 75 [NASA ADS] [Google Scholar]
  176. Pietrukowicz, P., Udalski, A., Soszyński, I., et al. 2012, ApJ, 750, 169 [NASA ADS] [CrossRef] [Google Scholar]
  177. Poci, A., McDermid, R. M., Zhu, L., & van de Ven, G. 2019, MNRAS, 487, 3776 [Google Scholar]
  178. Poggio, E., Drimmel, R., Lattanzi, M. G., et al. 2018, MNRAS, 481, L21 [Google Scholar]
  179. Poggio, E., Khanna, S., Drimmel, R., et al. 2025, A&A, 699, A199 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  180. Portail, M., Wegg, C., Gerhard, O., & Martinez-Valpuesta, I. 2015, MNRAS, 448, 713 [NASA ADS] [CrossRef] [Google Scholar]
  181. Portail, M., Gerhard, O., Wegg, C., & Ness, M. 2017, MNRAS, 465, 1621 [NASA ADS] [CrossRef] [Google Scholar]
  182. Queiroz, A. B. A., Chiappini, C., Perez-Villegas, A., et al. 2021, A&A, 656, A156 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  183. Queiroz, A. B. A., Anders, F., Chiappini, C., et al. 2023, A&A, 673, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  184. Quillen, A. C. 2002, AJ, 124, 722 [NASA ADS] [CrossRef] [Google Scholar]
  185. Quillen, A. C., Minchev, I., Sharma, S., Qin, Y.-J., & Di Matteo, P. 2014, MNRAS, 437, 1284 [Google Scholar]
  186. Raha, N., Sellwood, J. A., James, R. A., & Kahn, F. D. 1991, Nature, 352, 411 [Google Scholar]
  187. Rangwala, N., Williams, T. B., & Stanek, K. Z. 2009, ApJ, 691, 1387 [Google Scholar]
  188. Rattenbury, N. J., Mao, S., Sumi, T., & Smith, M. C. 2007, MNRAS, 378, 1064 [NASA ADS] [CrossRef] [Google Scholar]
  189. Reddy, B. E., Tomkin, J., Lambert, D. L., & Allende Prieto, C. 2003, MNRAS, 340, 304 [Google Scholar]
  190. Rich, R. M. 1988, AJ, 95, 828 [NASA ADS] [CrossRef] [Google Scholar]
  191. Rich, R. M. 1990, ApJ, 362, 604 [NASA ADS] [CrossRef] [Google Scholar]
  192. Rich, R. M., Origlia, L., & Valenti, E. 2007, ApJ, 665, L119 [NASA ADS] [CrossRef] [Google Scholar]
  193. Rix, H.-W., Chandra, V., Zasowski, G., et al. 2024, ApJ, 975, 293 [NASA ADS] [CrossRef] [Google Scholar]
  194. Rojas-Arriagada, A., Recio-Blanco, A., Hill, V., et al. 2014, A&A, 569, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  195. Rojas-Arriagada, A., Recio-Blanco, A., de Laverny, P., et al. 2017, A&A, 601, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  196. Rojas-Arriagada, A., Zoccali, M., Schultheis, M., et al. 2019, A&A, 626, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  197. Rojas-Arriagada, A., Zasowski, G., Schultheis, M., et al. 2020, MNRAS, 499, 1037 [Google Scholar]
  198. Ryde, N., Gustafsson, B., Edvardsson, B., et al. 2010, A&A, 509, A20 [CrossRef] [EDP Sciences] [Google Scholar]
  199. Saha, K., & Gerhard, O. 2013, MNRAS, 430, 2039 [NASA ADS] [CrossRef] [Google Scholar]
  200. Saha, K., Martinez-Valpuesta, I., & Gerhard, O. 2012, MNRAS, 421, 333 [Google Scholar]
  201. Saha, K., Gerhard, O., & Martinez-Valpuesta, I. 2016, A&A, 588, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  202. Saito, R. K., Zoccali, M., McWilliam, A., et al. 2011, AJ, 142, 76 [NASA ADS] [CrossRef] [Google Scholar]
  203. Sánchez-Blázquez, P., Rosales-Ortega, F. F., Méndez-Abreu, J., et al. 2014, A&A, 570, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  204. Sanders, J. L., Kawata, D., Matsunaga, N., et al. 2024, MNRAS, 530, 2972 [NASA ADS] [CrossRef] [Google Scholar]
  205. Sanders, J. L., Smith, L., & Evans, N. W. 2019a, MNRAS, 488, 4552 [NASA ADS] [CrossRef] [Google Scholar]
  206. Sanders, J. L., Smith, L., Evans, N. W., & Lucas, P. 2019b, MNRAS, 487, 5188 [CrossRef] [Google Scholar]
  207. Savino, A., Koch, A., Prudil, Z., Kunder, A., & Smolec, R. 2020, A&A, 641, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  208. Schultheis, M., Cunha, K., Zasowski, G., et al. 2015, A&A, 584, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  209. Schultheis, M., Rojas-Arriagada, A., García Pérez, A. E., et al. 2017, A&A, 600, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  210. Schwarz, G. 1978, Ann. Statist., 6, 461 [Google Scholar]
  211. Sellwood, J. A. 1981, A&A, 99, 362 [NASA ADS] [Google Scholar]
  212. Sellwood, J. A. 2014, Rev. Mod. Phys., 86, 1 [Google Scholar]
  213. Sellwood, J. A., & Gerhard, O. 2020, MNRAS, 495, 3175 [Google Scholar]
  214. Sellwood, J. A., & Wilkinson, A. 1993, Rep. Progr. Phys., 56, 173 [CrossRef] [Google Scholar]
  215. Semenov, V. A., Conroy, C., Chandra, V., Hernquist, L., & Nelson, D. 2024, ApJ, 962, 84 [NASA ADS] [CrossRef] [Google Scholar]
  216. Shen, J., Rich, R. M., Kormendy, J., et al. 2010, ApJ, 720, L72 [NASA ADS] [CrossRef] [Google Scholar]
  217. Siebert, A., Famaey, B., Minchev, I., et al. 2011, MNRAS, 412, 2026 [Google Scholar]
  218. Skokos, C., Patsis, P. A., & Athanassoula, E. 2002, MNRAS, 333, 847 [Google Scholar]
  219. Smith, L. C., Lucas, P. W., Kurtev, R., et al. 2018, MNRAS, 474, 1826 [Google Scholar]
  220. Sormani, M. C., Gerhard, O., Portail, M., Vasiliev, E., & Clarke, J. 2022, MNRAS, 514, L1 [NASA ADS] [CrossRef] [Google Scholar]
  221. Soto, M., Rich, R. M., & Kuijken, K. 2007, ApJ, 665, L31 [NASA ADS] [CrossRef] [Google Scholar]
  222. Stone-Martinez, A., Holtzman, J. A., Imig, J., et al. 2024, AJ, 167, 73 [NASA ADS] [CrossRef] [Google Scholar]
  223. Tahmasebzadeh, B., Zhu, L., Shen, J., et al. 2024, MNRAS, 534, 861 [NASA ADS] [CrossRef] [Google Scholar]
  224. Uttenthaler, S., Schultheis, M., Nataf, D. M., et al. 2012, A&A, 546, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  225. Vasiliev, E. 2019, MNRAS, 482, 1525 [Google Scholar]
  226. Vasiliev, E., & Valluri, M. 2020, ApJ, 889, 39 [Google Scholar]
  227. Vásquez, S., Zoccali, M., Hill, V., et al. 2013, A&A, 555, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  228. Vislosky, E., Minchev, I., Khoperskov, S., et al. 2024, MNRAS, 528, 3576 [NASA ADS] [CrossRef] [Google Scholar]
  229. Wang, Y., Mao, S., Long, R. J., & Shen, J. 2013, MNRAS, 435, 3437 [CrossRef] [Google Scholar]
  230. Wang, Y., Zhao, H., Mao, S., & Rich, R. M. 2012, MNRAS, 427, 1429 [NASA ADS] [CrossRef] [Google Scholar]
  231. Wegg, C., & Gerhard, O. 2013, MNRAS, 435, 1874 [Google Scholar]
  232. Wegg, C., Gerhard, O., & Portail, M. 2015, MNRAS, 450, 4050 [NASA ADS] [CrossRef] [Google Scholar]
  233. Whitford, A. E., & Rich, R. M. 1983, ApJ, 274, 723 [Google Scholar]
  234. Williams, M. E. K., Steinmetz, M., Binney, J., et al. 2013, MNRAS, 436, 101 [Google Scholar]
  235. Wylie, S. M., Gerhard, O. E., Ness, M. K., et al. 2021, A&A, 653, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  236. Zasowski, G., Cohen, R. E., Chojnowski, S. D., et al. 2017, AJ, 154, 198 [Google Scholar]
  237. Zhang, H., Belokurov, V., Evans, N. W., et al. 2024, MNRAS, 535, 2873 [Google Scholar]
  238. Zhao, H. 1996, MNRAS, 283, 149 [NASA ADS] [CrossRef] [Google Scholar]
  239. Zhao, H., Spergel, D. N., & Rich, R. M. 1994, AJ, 108, 2154 [NASA ADS] [CrossRef] [Google Scholar]
  240. Zhao, H., Spergel, D. N., & Rich, R. M. 1995, ApJ, 440, L13 [Google Scholar]
  241. Zhu, L., van de Ven, G., Leaman, R., et al. 2020, MNRAS, 496, 1579 [Google Scholar]
  242. Zhu, L., van de Ven, G., Leaman, R., et al. 2022, A&A, 664, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  243. Zoccali, M., Renzini, A., Ortolani, S., et al. 2003, A&A, 399, 931 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  244. Zoccali, M., Lecureur, A., Barbuy, B., et al. 2006, A&A, 457, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  245. Zoccali, M., Hill, V., Lecureur, A., et al. 2008, A&A, 486, 177 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  246. Zoccali, M., Gonzalez, O. A., Vasquez, S., et al. 2014, A&A, 562, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  247. Zoccali, M., Vasquez, S., Gonzalez, O. A., et al. 2017, A&A, 599, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  248. Zoccali, M., Valenti, E., & Gonzalez, O. A. 2018, A&A, 618, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
1

We increased the resolution to 2000 data points for the orbital frequency analysis in Section 4.

All Figures

thumbnail Fig. 1

On sky distribution of the APOGEE stars in our initial sample (left) and stellar mass-weighted density projection obtained by superposition of orbits of these stars (right). The top panels cover the entire disc of the MW, while the bottom ones zoom in on the MW bulge region. The orbit superposition technique not only reconstructs the 3D density distribution of the MW disc and bulge but also maps the abundance patterns of APOGEE stars across these components, effectively eliminating the limitations caused by the survey’s footprint and correcting for selection function biases. The reconstructed density shows a prominent boxy/X-shaped bulge structure, which is asymmetric along the longitude due to the bar’s orientation.
In the text
thumbnail Fig. 2

Variation of the inner MW density structure along the line-of-sight reconstructed using orbit superposition and APOGEE data. From top to bottom, the panels show the stellar density in (l,b) coordinates in 0.6 kpc-wide slices with increasing distances from the Sun, as marked in the bottom left of each panel. The middle panel shows the boxy bulge structure at the Galactic centre, while the upper and lower panels show the near and far sides of the bar-bulge, respectively. As we slice the MW bar inclined relative to the Sun – Galactic centre line, the off-centred X-shaped structure transitions from negative to positive longitudes with increasing distance.
In the text
thumbnail Fig. 3

Line-of-sight density structure of the MW bulge reconstructed using orbit superposition and APOGEE data. The top left panel shows the selection of the bulge fields marked by circles of different colours with the total column stellar density on the background. Other panels show the stellar density distribution along the line-of-sight at different latitudes, as is marked in the top left panel. The colour of the lines matches the colour of the circles in the top left panel. Only fields with b ≥ 0 are shown as the data are symmetric relative to the midplane. An animated version is available online.
In the text
thumbnail Fig. 4

Orbital frequency distribution of the MW bulge populations. The top panel shows the stellar mass-weighted distribution of the in-plane orbital frequency ratio, fr/fx. The elongated orbits, fr/fx ≈ 2, correspond to the bar, while the rest do not support the bar and represent the disc. The bottom panel shows the distribution of the vertical-to-radial orbital frequencies, fz/fx. The frequency ratio fz/fx ≈ 2 corresponds to the banana-like orbits. The lower values (fz/fx ∈ 1.4–2) are pretzel-like populations, dominating in the X-shaped bulge structure, whose extent along the bar major axis decreases with decreasing frequency ratio. Orbits with fz/fx > 2 correspond to the longer bar component (>2.5–3 kpc) outside the bulge.
In the text
thumbnail Fig. 5

Orbital decomposition of the MW bulge in face-on (top) and side-on (bottom) projections. The panels show the projected stellar density of orbits classified by their vertical-to-in-plane frequency ratio, fz/fx (see Fig. 4). The bar’s major axis is aligned along the X-axis in both projections. Each orbital class represents a distinct dynamical population, and the stellar mass associated with each class is indicated in the top left corner of the corresponding panels.
In the text
thumbnail Fig. 6

Kinematics of the MW bulge recovered using orbit superposition approach. The panels show the mean line-of-sight velocity (top) and the line-of-sight velocity dispersion (bottom) in the MW bulge region, limited by a heliocentric distance of 8.12 ± 3.5 kpc, around the Galactic centre. The grey contours reflect the shape of the stellar density distribution, while the white ones correspond to the constant velocity and velocity dispersion values. The bulge region shows a prominent rotation due to its disc origin, while the asymmetry of the radial velocity profile is due to projection effects. The velocity dispersion map shows a prominent peak extended perpendicular to the midplane.
In the text
thumbnail Fig. 7

Kinematics of the inner MW and bulge in the face-on projection. The panels show the mean heliocentric line-of-sight velocity (top) and velocity dispersion (bottom) maps at different latitudes. The grey contours show the total stellar isodensity levels, marking the extent and orientation of the bar. The dashed lines highlight the longitude values of l ± 20°. The yellow asterisks mark the position of the Sun.
In the text
thumbnail Fig. 8

Line-of-sight kinematics of the MW bulge. The panels show the mean line-of-sight radial velocity (top) and the line-of-sight velocity dispersion (bottom) versus longitude measured at different latitudes and in [Fe/H] bins, as marked in the header of the top panels. The line colour indicates the latitude range. The panels are based on the data within <3.5 kpc around the Galactic centre. In each panel, the black lines show the averaged velocity and velocity dispersion profiles for all stellar populations of the entire region.
In the text
thumbnail Fig. 9

[Mg/Fe]-[Fe/H] relation for the MW bulge. The left panel shows the distribution of APOGEE stars, while the right panel depicts the mass-weighted distribution obtained using the orbit superposition method within <3.5 kpc from the Galactic centre. The yellow lines show the distribution functions of [Mg/Fe] and [Fe/H] separately. The solid white line shows the border used to separate high- and low-α populations. The number of stars and stellar mass corresponding to the high- and low-α populations are marked in the panels. While the orbit superposition method assigns more mass to the metal-rich populations, both APOGEE data and the orbit superposition-based density reveal rather similar abundance distributions.
In the text
thumbnail Fig. 10

Comparison of the MW bulge MDF limited by different distances in the vertical direction from the midplane as marked in each panel. The blue line is the MDF weighed using the orbit superposition method versus the raw APOGEE sample counts in grey. The orbit superposition method and APOGEE show similar MDFs close to the midplane; however, some systematic differences become prominent once the input sample covers regions further out.
In the text
thumbnail Fig. 11

MDFs of different orbital families of the MW bulge. The orbital families are classified according to the in-plane to vertical frequencies ratio: inner boxy orbits (fz/fx = 1.4–1.6), X-shaped orbits (fz/fx = 1.6–1.9), banana orbits (fz/fx = 1.9–2.2) and long bar orbits (fz/fx > 2.2) orbits. In each panel, the distributions are normalised to the maximum value. For reference, the vertical red dashed lines mark the locations of the peaks in the MDF based on the microlensed dwarf and subgiant stars in the MW bulge from Bensby et al. (2017). The same grey histogram in the background shows the MDF of the disc orbits with |fz/fr − 2| > 0.1.
In the text
thumbnail Fig. 12

Shape of the MW bulge stellar density across different metallicity bins. The corresponding metallicity range and stellar mass are specified at the top of each panel. The top three panels represent the high-α populations, while the bottom panels depict the low-α populations. The mapping of these populations into the bulge reveals distinct trends where the contrast of the X-shaped/boxy structure increases with metallicity. High-α stars, corresponding to a kinematically hot thick disc component, are more widely dispersed in the bulge. In contrast, low-α stars are part of a colder population that is more closely confined to the midplane but still efficiently mapped into the bulge. The most (extremely) metal-rich populations exhibit a vertically extended, highly concentrated distribution towards the centre of the MW without displaying a prominent disc-like morphology. An animated version is available online.
In the text
thumbnail Fig. 13

Chemical abundance profiles across the bulge region. The left and middle panels in the top row show the mean stellar mass-weighted [Fe/H] and [Mg/Fe], respectively, while the right one depicts the stellar mass fraction of the low-α populations (see Fig. 9 for the definition). The bottom panels show one-dimensional profiles of the same quantities as a function of longitude and at different latitudes, marked with different colours. The bottom group of panels shows the same, but for the vertical profiles, along latitude. The black lines at the bottom show the latitude- or longitude-averaged profiles. The grey background lines in the 2D abundance maps depict the stellar isodensity contours.
In the text
thumbnail Fig. 14

Maps of the radial (top) and vertical (bottom) metallicity gradients in the MW bulge region. The maps show the slopes of the metallicity profiles measured in each pixel within ±1 deg in the corresponding direction. The radial and vertical gradient variations trace well the X-shaped and boxy density components, respectively. The gradients are steeper near the transition regions between high- and low-α populations. The grey background lines reflect the projected stellar density distribution.
In the text
thumbnail Fig. 15

Face-on maps of the mean stellar metallicity [Fe/H] distribution across the MW disc. From left to right, the maps represent stellar populations within horizontal slabs at increasing distances from the midplane. Near the midplane (two panels on left), the bulge region appears to be metal-poor because of the deficit of the metal-rich (low-α) populations expelled to higher distances perpendicular to the disc plane via bar/bulge instability. However, the actual metallicity of the MW bulge is higher compared to the surrounding disc once its full vertical extension is taken into account. The bar appears to be metal-rich along the major axis because kinematically colder low-α (metal-rich) populations are more efficiently captured by the bar’s gravitational potential compared to the hotter high-α (metal-poor) stars dominating along the bar minor axis.
In the text
thumbnail Fig. 16

Two-dimensional decomposition of the MW bulge in the [Mg/Fe]–[Fe/H] plane as a function of latitude based on the data inside 3.5 kpc from the Galactic centre. The left panels show the normalised stellar density maps, where the contours of different colours show the components identified using the 2D GMM approach. The grey histograms in the middle and right columns show the generalised [Fe/H]– and [Mg/Fe]–distribution functions, respectively. The coloured lines show the one-dimensional projections of the 2D-GMM components, while the black lines represent the sum of these components. In the middle row, the fractional contribution of each component is marked by the numbers of the same colour. Our analysis suggests that the [Mg/Fe]–[Fe/H] bulge composition is best represented by the five GMM components, where only two are the main ones and correspond to the bulk of high- and low-α populations. Two intermediate components allow for the decomposition of the transition between the main ones by fitting the α-knee and likely do not correspond to distinct bulge populations. The most metal-poor GMM component likely corresponds to the accreted populations.
In the text
thumbnail Fig. 17

Density structure of the MW bulge components identified using 2D GMM in the [Fe/H]–[Mg/Fe] plane. The top panels show the stellar mass density of each component obtained by multiplying the total stellar mass in the plane (background contours) by the 2D GMM probability for a given component. The middle and bottom rows show the mass×GMM-probability weighed density distribution in the XY and (l, b) coordinates, respectively. The density is scaled by the maximum value in each panel.
In the text
thumbnail Fig. 18

Energy-angular momentum structure of the MW bulge components identified using 2D GMM in the [Fe/H]–[Mg/Fe] plane. The top panels are the same as in Fig. 17 and are given for reference. The bottom group of panels show the normalised stellar mass for the entire component (top) and the region inside (middle) and outside (bottom) the bar. The numbers in the middle and bottom rows correspond to the mass fraction of a given component inside and outside the bar. The white dashed lines mark the minimum possible energy for a given value of the angular momentum.
In the text
thumbnail Fig. A.1

Stellar mass-weighted age distribution (top) and age-metallicity relation (bottom) in the MW bulge region (< 3.5 kpc from the Galactic centre). The metallicity range ([Fe/H] < −0.65 dex, Stone-Martinez et al. (2024)) with uncertain ages is shaded in the bottom panel and not included in the top distribution, accounting for about 13% of the bulge mass.
In the text
thumbnail Fig. A.2

Kinematics of the MW bulge populations identified using 2D GMM in the [Fe/H]-[Mg/Fe] plane. The top panels are the same as in Fig. 17 and are given for reference. The middle and bottom panels show each component’s line-of-sight velocity and dispersion maps, weighted by the mass multiplied by the GMM probabilities. The grey contours show the corresponding density distribution, while the white ones highlight the velocity/dispersion structure.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.