EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 695 (March 2025) A&A, 695 (2025) A220 Full HTML
Open Access
Issue
A&A
Volume 695, March 2025
Article Number A220
Number of page(s) 17
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202453304
Published online 25 March 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

Our knowledge about the Milky Way (MW) structure results from observations of individual stars. In 1785, William and Caroline Herschel made an extremely important discovery about the structure of our Galaxy. By counting stars in the sky, they found that most of the stars they could see lie in a flattened structure, and the number of stars was about the same in any direction. They concluded that the star system to which the Sun belongs is in the shape of a disc or wheel, and the Sun must be near the middle of this structure (Herschel 1785). Since then, the systematic counting of stars has been continuously increasing the number of known stars and improved our understanding of their host – the MW galaxy (e.g. Shapley 1918; Kapteyn & van Rhijn 1920; Aitken & Doolittle 1932; Gliese & Jahreiß 1979; van Altena et al. 1995; Hauck & Mermilliod 1998).

Two centuries later, in 1989, the European Space Agency (ESA) launched the Hipparcos spacecraft, which operated in orbit for 37 months, performing photometric and astrometric measurements and delivering positions, proper motions, and parallaxes of the MW stars at the millisecond level of accuracy. The observations resulted in the Hipparcos catalogue of approximately 100 thousand stars (Perryman et al. 1997), which enabled the precise determination of the distance to the Galactic center (Paczyński & Stanek 1998, the most accurate determination was provided later by GRAVITY Collaboration 2021) and the constraining of key kinematic parameters of the solar neighbourhood, revealing the complex motions of the nearby stellar populations (Dehnen & Binney 1998; Dehnen 1998). Auxiliary Tycho-1 and Tycho-2 catalogues (Høg et al. 1997, 2000) have made it possible to study millions of stars in the MW with more precise proper motions and, thus, have influenced many areas of modern astronomy (Perryman 2009). Nowadays, nearly a decade post-launch, it is hard to overstate the success of the still ongoing ESA space mission Gaia (Gaia Collaboration 2016, 2021). Perhaps the primary and the most fundamental outcome of this survey is our newfound ability to perceive the ensemble of the MW stars as a galaxy.

While the coverage by various stellar surveys continues to expand, thus reducing the ’Galaxia incognita’1 within the MW, the number of stars with known positions, stellar parameters, and kinematic information still peaks near the Sun, failing to reflect the underlying stellar density distribution. This discrepancy mostly stems from a fundamental limitation of our observations originating from the solar position within the MW disc midplane. Our location limits the volume of potentially observable stars, dependent on their apparent brightness, which is further influenced by interstellar medium (ISM) extinction and source crowding in the dense inner regions of the Galaxy, among other phenomena. Although spectroscopic observations across different wavelength ranges, such as infrared (e.g. APOGEE (Majewski et al. 2017) and soon MOONS (Gonzalez et al. 2020)), partially mitigate the impact of dust, they do not fundamentally alter the distribution of observed stars as a function of distance from the Sun. In this case, the observations cover only a fraction of the selected area and thus might not be representative of the stellar population of this area in the sky. To probe the underlying stellar populations, one needs to know what fraction of stars was observed relative to the total number of stars in the same area and correct for possible selection biases or demonstrate their absence (Chen et al. 2018; Nandakumar et al. 2017; Wojno et al. 2017).

Different techniques have been introduced to correct for the observational biases, assuming certain completeness of the observed sample as a function of some parameters (see, e.g. Bovy et al. 2014; Nidever et al. 2014; Sharma et al. 2014; Mints & Hekker 2019; Boubert et al. 2021; Boubert & Everall 2022; Castro-Ginard et al. 2023; Wang et al. 2024). For instance, Chen et al. (2018) used weighting schemes to correct for the metallicity bias introduced by the target selection in their sample of SEGUE main-sequence turn-off stars and G and K dwarf stars (see also Schlesinger et al. 2012). Nandakumar et al. (2017) found no large selection function effect on the metallicity distribution function (MDF) and the vertical metallicity gradient for APOGEE (DR13, Albareti et al. 2017). However, their sample is restricted from 7 to 9 kpc and within 2 kpc from the midplane. Mints & Hekker (2019) showed that the effect of the selection function is more prominent for surveys that have a complex target allocation strategy and whose footprint is more patchy (e.g. APOGEE and Gaia-ESO), while for surveys with contiguous footprints, such as GALAH, RAVE, and LAMOST, the selection function has almost no effect in the extended solar vicinity (Wojno et al. 2017; Chen et al. 2018).

Generally speaking, various sorts of analyses of the MW stellar populations rely on relatively small sub-samples of stars drawn from the spectroscopic catalogues by applying user-defined cuts motivated by the scientific problem and the desired precision of stellar parameters. Hence, the selection function of that sample of stars is always unique and thus requires particular treatment; otherwise, the recovered kinematics and stellar abundance variations may not be representative of particular regions (or components) in the MW. Therefore, there is still a strong need for a selection function independent analysis of data from ongoing and future spectroscopic surveys (Rix & Bovy 2013).

Over the past years, the availability of large amounts of high-precision observational data, particularly from recent and ongoing ground-based spectroscopic surveys (e.g., Gaia-ESO (Randich et al. 2022), RAVE (Steinmetz et al. 2020), LAMOST (Zhao et al. 2012), APOGEE (Majewski et al. 2017), GALAH (Buder et al. 2021)), and space mission Gaia (ESA Gaia Collaboration 2016, 2023b), has revolutionised our understanding of the Galactic stellar populations. Various datasets offer a wealth of information that can be leveraged to reconstruct the MW present-day structure (Poggio et al. 2018; Antoja et al. 2018; Khoperskov et al. 2020; Cautun et al. 2020; Anders et al. 2022) and uncover some episodes of its assembly histories (Belokurov et al. 2018; Haywood et al. 2018; Helmi et al. 2018; Kruijssen et al. 2019; Helmi 2020). In parallel, using this information, multiple models of the MW have been developed, which one can divide into ones tailored to constrain the formation and enrichment history of the Galaxy (see, e.g. Schönrich & Binney 2009; Sharma et al. 2021; Frankel et al. 2019; Prantzos et al. 2023) and models aimed to recover its present-day structure (see, e.g. Drimmel & Spergel 2001; Robin et al. 2003; Binney 2012; Czekaj et al. 2014; Trick et al. 2016; Price-Whelan et al. 2021; Binney & Vasiliev 2024, 2023). However, understanding the MW as a galaxy remains a formidable challenge with its complex interplay of stars, gas, and dark matter, encoding the underlying structure and kinematics of the Galaxy.

Outside the MW, significant progress has been made in developing sophisticated techniques to probe the dynamics and structure of external galaxies. One such method that has attracted significant attention is the Schwarzschild orbit superposition technique (Schwarzschild 1979, see also Pfenniger 1984; Richstone & Tremaine 1988; Rix et al. 1997; van der Marel et al. 1998; Thomas et al. 2004; Valluri et al. 2004; Cappellari et al. 2006; van den Bosch et al. 2008; Vasiliev & Valluri 2020a; Jethwa et al. 2020; Zhu et al. 2020; Neureiter et al. 2021). This powerful approach offers a unique way to model the gravitational potential and kinematics of galaxies by superposing a library of stellar orbits. The application of the Schwarzschild method has seen remarkable success in modelling the dynamics of external galaxies, providing insights into their dark matter content (see, e.g. Rix et al. 1997; van de Ven et al. 2010; Thomas et al. 2011; Comerón et al. 2023; Lu et al. 2024), bar structures (see, e.g. Wang et al. 2013; Vasiliev & Athanassoula 2015; Vasiliev & Valluri 2020b; Petersen et al. 2021; Tahmasebzadeh et al. 2022, 2024; Kacharov et al. 2024), central black holes (see, e.g. van der Marel & van den Bosch 1998; Verolme et al. 2002; Gebhardt et al. 2003; Houghton et al. 2006; Krajnović et al. 2009; Kormendy & Ho 2013; Saglia et al. 2016; van den Bosch 2016), assembly history (see, e.g. Ding et al. 2023; Poci et al. 2021; Zhu et al. 2022), and orbital structure (Krajnović et al. 2005; van den Bosch et al. 2008; Zhu et al. 2018b; den Brok et al. 2021; Tahmasebzadeh et al. 2024).

In this paper, we present a novel orbit superposition approach based on the APOGEE-like mock observations of a simulated MW-like barred galaxy. Such an approach has not been adapted in the context of the MW so far, as it presents specific challenges, but also new opportunities for a deeper exploration of observational data very much needed in light of the coming new generation of spectroscopic surveys, such as 4MOST (de Jong et al. 2019), SDSS-V (Almeida et al. 2023), MOONS (Gonzalez et al. 2020), and WEAVE (Jin et al. 2024). Our method follows the original orbit superposition approach first introduced by Schwarzschild (1979), where a library of stellar orbits was used to make a model of a triaxial elliptical galaxy without constraining the solution by kinematic information. We demonstrate that under reasonable assumptions regarding the galactic potential, this method makes it possible to correct the APOGEE-like mock data for the selection function and, even more, to go beyond a specific survey spatial footprint and recover the unbiased kinematics and abundance information across the entire galaxy. In the follow-up works, we apply this technique to reconstruct the three-dimensional stellar density distribution, kinematics, and chemical abundance composition of different components of the MW galaxy (Khoperskov et al. 2024a,b).

The structure of the paper is as follows. In Section 2 we explain in detail the methodology of the orbit superposition method in the context of the MW-like data and construction of the mock APOGEE-like data. In Section 3, we present the results of the orbit superposition using two mock catalogues based on different selections of star particles, and we test the reconstructed kinematics of stars and the ability of the method to reproduce the stellar metallicity distribution. We discuss the limitations and future prospects of our approach in Section 4. The summary of the paper is given in Section 5.

2 Methodology

2.1 Orbit superposition approach

The objective of a classic Schwarzschild orbit superposition method, in the context of external galaxies, is to find the parameters of gravitational potential (dark matter (DM) and baryons contributions) and to recover the phase-space structure of stellar components constrained by the observed distribution of stars (imaging) and line-of-sight kinematic information. Solving this problem requires finding a consistent solution to both the collisionless Boltzmann equation and the Poisson equation. This is usually accomplished by combining some basis orbits (or library of orbits) with corresponding weights to reproduce a density (or light) distribution that aligns with the potential used to generate the orbits. The total galactic potential can be represented by a sum of contributions of different components (dark matter, stars, gas, and central black hole): Ψtot(r)=ΨDM(r)+Ψs(r)+Ψg(r)+ΨBH(r).$\Psi_{tot}(\mathbf{r})=\Psi_{D M}(\mathbf{r})+\Psi_{s}(\mathbf{r})+\Psi_{g}(\mathbf{r})+\Psi_{B H}(\mathbf{r}).$(1)

While in external galaxies the total potential Ψtot(r) is unknown a priori, thanks to the amount and various sources of data, the total MW mass (Erkal et al. 2019, 2020; Cautun et al. 2020; Vasiliev et al. 2021), the contribution of different components (Bovy 2015; McMillan 2017; Pouliasis et al. 2017; Eilers et al. 2019), halo shape (Vera-Ciro & Helmi 2013; Posti & Helmi 2019; Vasiliev et al. 2021), and even the 3D structure of the stellar disc, including the complex bar/bulge inner region (Wegg & Gerhard 2013; Wegg et al. 2015; Portail et al. 2017b,a; Sormani et al. 2022), are constrained quite well. Therefore, if we assume that the total potential of the galaxy and, more specifically, its stellar component are known or well constrained by certain independent measurements, then the 3D stellar density distribution ρs(r) can be easily recovered from the Poisson equation: ∇2Ψs(r)=4πGρs(r),$\nabla^{2} \Psi_{s}(\mathbf{r})=4 \pi G \rho_{s}(\mathbf{r}),$(2)

and only the complete phase-space distribution function remains unknown. In order to obtain the complete distribution function, one can use the orbits of stars calculated in the total potential of the galaxy, Ψtot(r). This can be done using separable triaxial potentials, where all orbits are regular, conserve three integrals of motion and can be calculated analytically (de Zeeuw 1985; Schwarzschild 1993). In practice, the orbital library is usually limited by several families of orbits (Statler 1987; Zhao 1996; van den Bosch et al. 2008): short-axis tubes, outer and inner longaxis tubes, box orbits and so on. In the case of the MW, once a certain potential is adopted, one can integrate orbits of real stars using positions and velocities available from Gaia and groundbased spectroscopic surveys as the initial conditions. Therefore, the library of orbits for the orbit superposition modelling can be based on the orbits of real objects and include all sorts of orbital families governed by the potential and the input phasespace distribution function.

In order to fill in the unknown (or unobserved by a given Galactic survey) parts of the distribution function, following the canonical Schwarzschild modelling, one can solve the linear equation: ρs(r)=∑iNwiρ(r)i,$\rho_{s}(\mathbf{r})=\sum_{i}^{N} w_{i} \rho(\mathbf{r})_{i},$(3)

where wi are the orbit weights, and each orbit has a 3D density distribution ρ(r)i, representing how much time a star spends in a given position. The solution of Eq. (3) provides us with the weights wi for the orbit of each star in the observational sample whose superposition results in the adopted stellar density.

It is crucial to emphasise the key distinctions between the orbit superposition approach we propose here for the MW and traditionally employed Schwarzschild models used to derive parameters of external galaxies. The primary divergence lies in the assumptions made regarding the knowledge of the gravitational potential, including individual contributions from dark matter, stars, and gas. For the MW, it is reasonable to assume that the total potential (Ψtot(r)), as well as the contributions from stellar component (Ψs(r) or ρs(r)), are known. Conversely, in the case of external galaxies, the gravitational potential is ab initio unknown, and its parameters are typically constrained using available kinematic and/or luminosity data. This can be achieved through an iterative fitting procedure, allowing for adjustments to the parameters of the potential (disc and DM halo scales) and its functional form.

Therefore, the question we are seeking an answer to is not the parameters of the MW potential but rather its unbiased phasespace distribution function, which is normally obscured by the selection function of different surveys and their spatial incompleteness. In particular, the resulting orbit superposition solution makes it possible to map stellar kinematics beyond the survey footprint, essentially encompassing the entire galaxy. The resulting velocity distribution will remain unbiased against the survey selection function. If stellar parameters such as abundances and ages are available for stars in the sample used to construct the orbital library, it becomes possible to transfer these ‘labels’ along the orbits and ‘paint’ their distribution across the galaxy, which again will take into account the weights of populations with different stellar parameters, as they are evolutionarily linked to the orbits of stars.

2.2 Simulation and potential reconstruction

For the above-described method verification, we use a single snapshot from a simulation of an isolated MW-like barred galaxy from an N-body/hydrodynamical simulation, presented in Vislosky et al. (2024) and used to match the present-day MW kinematics. In this work, we do not re-scale the simulation to match the MW, as we do not aim to reproduce any MW parameters but rather to validate the orbit superposition approach. The simulation also includes star formation and the basic chemical abundance cycle (Khoperskov et al. 2021; Haywood et al. 2024) from which we use only metallicity information to test the goodness of our approach in the reconstruction of stellar parameters evolutionary-coupled with the orbits of stars. We therefore select a snapshot when the bar and the X-shaped bulge are already well developed, and the simulated galaxy looks visually similar to the MW (see Fig. 1 in Vislosky et al. 2024). However, in order to highlight the limitations of our approach, we have deliberately chosen a snapshot with relatively strong spiral arms and not relaxed vertical motions of stars, which by construction cannot be reproduced by Schwarzschild orbit superposition methods assuming a dynamical equilibrium. For the snapshot analysis, we orientated the bar by 27∘ relative to the Sun-galactic centre line.

The distinct feature of the MW stellar component is the bar, and in this work, we aim to include such a non-axisymmetric component in the orbit superposition modelling. In such a case, we will construct the orbital library in a rest frame rotating with the rotational frequency of the bar. This allows us to consider the positions of stars along the orbits time-independent because the bar position is always fixed to the present-day orientation relative to the Sun. We illustrate this coordinates transformation in Fig. 1 where we show two orbits of random disc particles in a barred potential in the inertial (left) and rotating (right) frames. From the perspective of Galactic stellar populations, replacing individual stars with orbits implies that for each observed star, numerous other stars exist that are unobserved for various reasons yet follow identical orbits. That is why using a rotating rest frame is vital for ‘observing’ these hypothetical stars at the same moment in time. It also allows us to avoid the phase-mixing of the distribution function naturally happening in an evolving potential, for instance, in a rotating bar.

As the first step, we took particle positions and masses and create a smooth approximation of the galactic potential using the AGAMA code (Vasiliev 2019). Since the simulation includes star particles of different ages and, thus, different kinematics, and we aim to reconstruct the potential of the bar and the X-shaped structure in as much detail as possible for the 3D stellar density (potential) approximation, we use the spherical-harmonic expansion with lmax = 20 and mmax = 20. We recall that in such a case, we are unable to take into account spiral arms in the approximated version of the gravitational potential, which will result in minor but expected deviations of stellar kinematics and abundance distribution obtained in the orbital space compared to the original snapshot data. We approximate the potential of the simulated galaxy using three components: DM, stars and gas. Therefore, we can obtain the approximated 3D density distributions for each component separately. In Fig. 2 we compare the face-on stellar surface density in the snapshot (left) against its approximation with the spherical-harmonic expansion (right). The residual map demonstrates a very good approximation of the stellar disc mass distribution, including the 3D structure of the inner region and the central bar.

thumbnail Fig. 1

Example of two orbits (red and blue) of disc stars in a barred potential in inertial (left) and rotating with the bar (right) reference frames. In the bottom panels, the orbits are colour-coded to indicate the progression of time. The ovals in the figure represent the bar orientation, with three different configurations colour-coded to denote changes in its orientation over time in the inertial frame. Notably, the orientation of the bar varies with time in the inertial reference frame while, by definition, it remains fixed in the rotating frame. Therefore, the orbits in the rotating frame allow us to exclude the time dependence of particle positions along the orbits and prevent phase-mixing, enabling us to use them in the orbit superposition approach.

2.3 Selection of APOGEE-like samples from simulation

So far, we do not have access to the 6D phase-space information based on all stars across the MW, and the spatial selection function (or distribution function) of large-scale disc surveys is quite complex. Therefore, it is vital for us to understand how much the results obtained using orbits of a limited sample of stars provide meaningful information about a larger area across the Galaxy. In the following, we use the snapshot simulations to extract particles, which will be used to make mock observations. The spatial distribution of these stars and their metallicities will be used as inputs for the orbit superposition models. The main purpose is to test how a selection bias (spatial footprint) affects the results of the modelling. We, therefore, mimic the APOGEE-like selection footprint by sampling the N-body simulation particles according to the distance distribution in the APOGEE catalogue in each HEALpix, assuming a HEALpix level of 4, which is a fairly low sky resolution, but it gives a reasonably good mock catalogue by design reproducing the spatial footprint of the MW survey.

In order to test the effect of the survey footprint and, thus, different samples of orbits, we created two mock catalogues. The first one (Mock 1), is based on the APOGEE giant stars selection, where we adopt log g < 2.2. Limiting to smaller log g values, the selection minimises potential systematic uncertainties in abundance measurements, but the higher luminosity of these giant stars allows them to cover a larger area across the MW (see, e.g. Eilers et al. 2022; Cerqui et al. 2023; Imig et al. 2023; Haywood et al. 2024). The second sample (Mock 2) includes stars with 2.5 < log g < 3.6, making it relatively local but still covering a substantial area around the Sun. This sample of the MW stars is often characterised by a relatively lower age determination uncertainty (Anders et al. 2023; Leung et al. 2023). The number of star particles in the APOGEE-like mock selections includes 40.354 (Mock 1) and 39.596 (Mock 2), while the total number of star particles in the snapshot is 7.343.898. Both mock samples are shown in Fig. 3, where the density peaks at the solar position and shows typical pencil-beam-like features. Below, we refer to the results of the orbit superposition modelling as Model 1 and Model 2, based on Mock 1 and Mock 2, respectively.

3 Results

3.1 Orbit weight calculation

In Section 2 we described a general idea behind the orbit superposition modelling in the context of the MW. In this section, we focus on more practical steps that we developed using the simulated MW-like galaxy.

We used the approximated 3D stellar density distribution of the galactic stellar disc component, whose calculation is detailed in Sect. 2.2. The stellar density distribution (ρs(r)) was then binned onto a 3D cube of 30 kpc in size (from −15 kpc to 15 kpc), with 50 bins along each direction. This Cartesian grid was also employed to determine the contribution of individual orbits (ρ(r)i).

The orbits of stars with initial conditions from both Mock 1 and Mock 2 were integrated over 5 Gyr in a reference frame rotating at a constant angular speed corresponding to the measured bar pattern speed of 22 km s−1 kpc−1. We used an instantaneous value of the bar pattern speed measured directly using the bar orientation in the previous and subsequent snapshots. Next, each orbit was discretised into 500 data points (3 coordinates and 3 velocities) equally cadenced in time, which were then translated into density values on the 3D Cartesian grid mentioned above. Furthermore, to ensure the axial symmetry of the reconstructed density, each orbit was mirrored along every axis. We have experimented with other possible symmetries and did not find any substantial differences (see the detailed exploration of the orbit mirroring on the Schwarzschild modelling solution in Thater et al. 2022); thus, we stick to the simple mirroring of the orbits which also reduces computational costs. We note that the mirroring of each orbit does not imply that we assume all orbits to be symmetric relative to the principal axes. In fact, it means that for each asymmetric orbit, there are corresponding mirrored counterparts in our library.

We therefore have constructed the discrete versions of ρs(r) and ρ(r)i where r = {xj, yj, zj}j = 1, …, 503. Generally speaking, a solution of Eq. (3) is unique if the number of variables wi matches the number of the discrete density values ρs(r), and the basis functions (library of orbits) represent an orthogonal basis (see discussion in Schwarzschild 1979). However, our approach does not guarantee these conditions, potentially leading to the degeneracy of the solution (Zhao 1996). To address this issue, for a given orbital library, we calculate non-negative weights wi using the least square minimisation for 10% randomly (without repetitions) selected subset of orbits (we experimented in 2–30% range and find 10% to be optimal in terms of solution stability and computational efforts), ensuring that every subset fits well the entire 3D stellar density ρs(r). This procedure has been repeated ten times for the whole sample of orbits, so we obtained ten realizations of weights for each orbit. Then, in order to minimise a possible degeneracy of the solution, we use the mean values of weights across these different realizations for further analysis. In Fig. 4, we present the variation of weights across different realizations for Model 1, displaying the weight distribution for five out of ten realizations (top row) and the mean weights (top right). The weight distributions exhibit considerable width, spanning nearly four orders of magnitude, but the distributions closely resemble each other across the various realizations. Upon comparing the mean weights with those obtained in individual realizations, we observe a close one-to-one relation (middle row of Fig. 4); however, some scatter remains noticeable. Analysis of the ratio between the mean weights and the standard deviation across ten realizations reveals that individual measurements typically do not deviate by more than a factor of 2–3 (middle row, right panel). Nevertheless, the stellar density reconstructed using the individual weights is identical to that obtained using the mean weights (bottom panels). Therefore, we conclude that employing the averaged weights while helping to mitigate the degeneracy of the orbit superposition solution results in the desired solution for the disc stellar density.

thumbnail Fig. 2

Approximation of the stellar disc density from the simulation. The surface density of star particles in a snapshot is shown in the left panel. In the middle panel, we present the surface density of the stellar disk obtained from the potential approximated using multipole expansion (see Section 2.2 for details). The rightmost panel shows the difference between the two maps, divided by the snapshot stellar density. Notably, the approximated density (and potential) reproduces well the large-scale structure of the simulated galaxy, including the bar 3D structure; however, by construction, small-scale asymmetries such as spiral arms are not captured. The black circle shows the solar radius of 8.2 kpc. The bar orientation is chosen to be 27∘ relative to the horizontal axis, reproducing the MW bar positioning.
thumbnail Fig. 3

Initial selection of star particles from the simulation mimicking the APOGEE footprint. In the left panel, the surface density of all star particles in the snapshot is shown. The middle and right panels show two mock catalogues which qualitatively reproduce the APOGEE DR17 giants (Mock 1) and red clump (Mock 2) stars footprints, respectively (see Sect. 2.3 for details). The number of star particles in the snapshot is 7.343.898, while the APOGEE-like selections include 40.354 (Mock 1) and 39.596 (Mock 2). In all the panels, the bar orientation is chosen to be 27∘ relative to the horizontal axis, reproducing the MW bar positioning.
thumbnail Fig. 4

Analysis of weights of orbits across various orbit superposition model realizations using Mock 1. The five left panels of the top row illustrate the weight distribution in five realizations of the orbit superposition using different 10% randomly sampled orbits, while the rightmost column shows the distribution of mean weights averaged across these realizations. In the second row, the five left panels reveal the relationship between mean weights and those obtained in different realizations of orbit superposition. The rightmost panel in this row presents the ratio between the standard deviation of weights and their mean. The bottom panels depict the relative difference between the stellar surface density obtained in various realizations and the surface density obtained using the mean orbital weights. Although a certain level of stochasticity in the distribution of the weights is seen, all the realizations successfully recover the same density distribution, which is nearly identical to the density based on the mean weights. The observed discrepancy in the orbital weight distribution across different realizations is likely attributed to a non-orthogonal orbital library. We, therefore, recommend using weights averaged across several orbit superposition realizations for more stable results.

3.2 Stellar density reconstruction

The results of the stellar density reconstruction are depicted in Fig. 5. In the left panel, we present a stacked distribution of all orbits obtained in our models (∑iNρ(r)i)$(\sum_{i}^{N} \rho(\mathbf{r})_{i})$. Notably, a striking difference emerges between our mock APOGEE-like selections. Model 1 yields an orbital density peaking at the galactic centre, indicating that our sample includes a significant number of stars with apocenters within our initial selection. Consequently, when mapping the orbits of these stars, they naturally contribute to the inner galaxy. Just beyond the bar, there appears to be an artificial (as it does not appear in the snapshot data) ring-like overdensity, likely populated by stars on nearly circular orbits and dominated at the solar radius. Although the model based on Mock 2 exhibits a similar orbital density, the bar-like overdensity is less pronounced. This difference is evidently due to the mock selection being more tightly confined to the solar radius, with less extension towards the centre of the galaxy. This result agrees with the finding of Wylie et al. (2022), who used the unweighted orbit superposition and found a similar inner ring-like structure in the MW between the planar bar and corotation.

The middle panels of Fig. 5 display the best solutions for the recovered stellar density, or, in other words, the weighted density distribution of all orbits in each model. As seen, the unweighted orbits based on Mock 1 and 2 are not identical (see left panels), as their sample is controlled not only by the gravitational potential but also by the initial distribution function defining the orbital library in each model. Here, the distinction between the orbit superposition solutions for Model 1 and Model 2 is quite evident. Model 1 exhibits a striking resemblance to the approximated stellar density, a similarity further emphasised in the residual density map depicted in the rightmost panel where the relative difference between the two densities is less than 5%. For Model 2, although the overall density structure looks similar to a barred galaxy, the recovered density distribution is far from the one desired. In particular, the bar overdensity is less prominent than the approximated density, while the surrounding inner region is too dense and shows a prominent spheroid-like component in the centre. The quality of the stellar density reconstruction in Model 1 and Model 2 is somewhat reflected by the distribution of the mean weights, shown in Fig. 6. The best solution based on Mock 2 exhibits an extremely broad distribution of weights, spanning approximately nine orders of magnitude, compared to the more compact distribution observed in Model 1. This likely serves as an indicator of the degeneracy of the orbit superposition model, which can be helpful if the target density is weakly constrained. Therefore, we suggest that such a broad distribution of weights may signal a problematic solution alongside poor convergence of the stellar density solution derived (see Fig. 5) from the orbit superposition.

These best solutions for the stellar surface density in Fig. 5 demonstrate a rather obvious result that a superposition of not every random sample of stars is sufficient to recover the density structure of the galaxy. One can consider an extreme scenario where, for various reasons, the selection of stars includes only cold-orbit stars, and none of them passes through the Galactic centre. In terms of Schwarzschild modelling, this implies that the orbital library does not include the full range of orbital families in a given potential. However, it appears that present-day spectroscopic data from APOGEE DR 17 (Mock 1) alone are adequate to explore the entire MW using orbit superposition methods. Obviously, one can leverage datasets such as Gaia DR3, where the spatial coverage of stars with relatively precise 6D phasespace information extends even beyond the Galactic centre. If the high-precision chemical abundance information from larger samples of stars is needed, a promising approach might be to use the data from several spectroscopic surveys using homogenised spectroscopic datasets (Thomas et al. 2024).

At this point, one might consider disregarding the analysis of Model 2, given its failure to recover the stellar density precisely. However, since the ground truth about the real MW is unknown, it remains crucial to illustrate potential biases and their origin that could come out in the orbit superposition models for a large but obscured initial phase-space selection. Also, it is important to showcase situations in which the results obtained can still offer valuable information.

thumbnail Fig. 5

Results of the orbit superposition method. The left panels display the stellar density of orbits with constant weights, essentially presenting the stacked library of orbits of all star particles from Mock 1 (top) and Mock 2 (bottom). The middle panels depict the stellar density reconstructed using weighted orbit superposition where the weights are calculated as the mean across ten orbit superposition realizations (see details in Sect. 3.1). The rightmost panels illustrate the relative residuals between the orbit superposition reconstruction and gravitational potential-based approximation (see Fig. 2 and Sect. 2.2 for more details).
thumbnail Fig. 6

Mean weights distribution in two orbit superposition models obtained in ten realizations for the APOGEE-like Mock 1 (blue) and Mock 2 (blue; see details in Sect. 2.3). The weight values are given in relative units. The broad distribution of weights in Model 2 indicates a worse convergence of the solution compared to Model 1.
thumbnail Fig. 7

Distribution of velocity components (VR in left, Vϕ in the middle, and Vz in the right) in the galactic rest frame within three 2 kpc-wide annuli centred at 1 (top), 8 (middle), and 11 kpc (bottom). The snapshot data are represented by the grey-filled areas, while the results of the orbit superposition based on Mock 1 and Mock 2 are depicted in blue and red colours, respectively. The blue and red-filled areas correspond to the initial mock data. For each panel, the numbers correspond to the mean and standard deviation of the velocity distribution in km s−1. We emphasise that the orbit superposition method we employ aims to recover the approximated stellar density distribution; the figure shows that Model 1 reproduces the kinematic information very well. At the same time, this approach is unable to capture spiral arms or the non-equilibrium kinematic features of the disc and bar. Furthermore, the effectiveness of the orbit superposition technique is contingent upon the accuracy of the potential approximation.

3.3 Kinematics

In the following section, we analyze the stellar kinematics of the simulated galaxy recovered using the orbit superposition approach. It is important to note that, unlike traditional Schwarzschild methods, no kinematic information from our sample of stars has been used to calculate the weights of the orbits.

First, we examined the distribution of the orbit massweighted velocity components at different galactocentric radii and compare them to the ground truth (snapshot data). In Fig. 7, we display the distributions of radial (VR), azimuthal (Vϕ), and vertical velocities (Vz) in three circular two kpc-wide annuli centred at 2, 8, and 11 kpc. We recall that our orbit superposition approach aimed to recover the approximated stellar density distribution rather than the snapshot density distribution itself. Thus, despite a very good density approximation (see Fig. 2), some differences between the snapshot data and the reconstructed stellar kinematics may be expected.

The orbit superposition model based on Mock 1 reproduces the snapshot kinematics with great fidelity, and we observe only a few minor discrepancies. In particular, the radial velocity component distribution (VR, left column) almost perfectly recovered at all radii highlighted by the match between the blue lines (Model 1) and grey shaded area (snapshot data). At the same time, in the inner galaxy, the azimuthal velocity distributions (Vϕ, top middle panel) do not perfectly align; it appears that the reconstructed velocity distribution is shifted towards lower absolute values, as evidenced by the difference in the mean velocity (−83.5 km s−1 versus −88.9 km s−1 in the snapshot data). Since we do not observe such a mismatch at larger radii (middle and bottom panels of the middle row), it is likely that the approximated density, which we recover using orbit superposition, does not perfectly reproduce some small details (lopsidedness or a weak misalignment relative principal axes) of the inner bar region. Additionally, there are slight discrepancies in the vertical velocity distribution, where the central peaks of the reconstructed velocity distributions are slightly lower. Nevertheless, these differences are hardly noticeable when comparing the mean and dispersion values, as indicated by the numbers in each panel.

While the mean and velocity dispersion values in Model 2 (red colour) remain very close to those of the snapshot data, as anticipated, it fails to capture multiple features of the corresponding distribution functions. Although the overall shapes of the radial and vertical velocity distributions (left and right columns) appear to be in qualitative agreement with the snapshot data, they do not reproduce the tails and central values of these distributions. A more striking artefact emerges in the azimuthal velocity distribution, which reveals a secondary component at low absolute velocity values. This kinematic feature, which either rotates slowly in the centre or remains non-rotating at larger radii, is naturally associated with a spheroid-like component, seen already in the stellar density solution in Fig. 5 (see Model 2 in the bottom middle panel). The emergence of such a feature in the orbit superposition solution is readily interpretable, as this method leverages existing star orbits to populate the 3D stellar density within a specified volume. To determine the stellar density solution at the centre, significant weight must be assigned to orbits that pass through the mock footprint. However, Mock 2 does not seem to capture a substantial number of stars near the galactic centre, relying only on hot radial orbits for the possible solution. These orbits possess minimal angular momentum, and when substantial mass is attributed to them, their cumulative effect manifests as an artificial spheroid with low net rotation. This underscores the critical importance of the initial star sample selection for the orbit superposition method, as it should allow the construction of a representative distribution function. Consequently, the failure of Model 2 is anticipated, yet it effectively demonstrates the machinery behind the orbit superposition approach.

As one can notice in Fig. 7, the azimuthal velocity distributions exhibit a more complex behaviour compared to the radial and vertical velocities. Therefore, it is vital to assess how well the azimuthal velocity distribution as a function of distance is recovered. This distribution is also quite interesting because it gives an idea of how well the rotation curve of the galaxy can be reconstructed. Figure 8 shows the density maps in the r − Vϕ plane for the snapshot data (left) and reconstructed densities in Model 1 (middle) and Model 2 (right). The mean azimuthal velocity and ± velocity dispersion from the snapshot data are shown with cyan lines across the panels. Model 2 shows a number of features non-observed in the snapshot data, where the low-Vϕ secondary overdensity in Fig. 7 (right) represents the artificial spheroid whose origin is discussed above. On the contrary, Model 1 shows a very good agreement with the rotational velocity distribution in the snapshot. Some small-scale features can be different; however, this can be explained by the lack of spiral arms in the reconstructed velocity distribution, which, for instance, are responsible for some of the features, such as diagonal ridges (Antoja et al. 2018; Hunt et al. 2019; Khoperskov & Gerhard 2022). Another source of discrepancy stems from the fundamental assumption underlying the orbit superposition method: the reliance on the galaxy being in dynamic equilibrium. This assumption is, of course, not exactly correct. For instance, for the orbit integration, we adopted an instantaneous value of the bar pattern speed; however, bars in the MW-mass galaxies are expected to slowdown (Debattista & Sellwood 1998, 2000; Athanassoula 2002). On top of this effect, the parameters of bars (length, strength, pattern speed) oscillate on a short time scale during a connection-disconnection cycle with slowlyrotating spiral arms (Hilmi et al. 2020; Vislosky et al. 2024; Marques et al. 2025). Nevertheless, both models capture well the bar OLR-related ridge whose high-Vϕ tail is seen on top of the distribution at ≈8 kpc (Fragkoudi et al. 2019, 2020) suggesting that the instantaneous bar pattern speed value works reasonably well in terms of the large-scale disc kinematics. At the same time, it allows us to study the second-order non-equilibrium effects by subtracting the kinematic orbit superposition solution from the original dataset (snapshot data or observations).

We conclude this section by presenting the in-plane distribution of the mean and velocity dispersion components. As demonstrated earlier, orbit superposition Model 2 fails to capture several important kinematic characteristics of the simulated galaxy, but it introduces some artificial features. Therefore, in Fig. 9, we solely compare the snapshot data with velocity distributions from Model 1.

The radial velocity component of the simulated galaxy exhibits a well-known quadrupole or so-called butterfly pattern in the bar region (Ceverino & Klypin 2007; Fragkoudi et al. 2020), a phenomenon that has recently been directly observed in the MW as well (Bovy et al. 2019; Gaia Collaboration 2023a). Outside the bar region, we see a large-scale radial velocity pattern associated with the presence of a spiral structure in the simulation. The reconstructed in Model 1 radial velocity effectively captures the butterfly pattern. It is worth noting that Mock 1 includes only a small fraction of stars from this inner region (see Fig. 3); however, the superposition of orbits enables us to access the entire pattern. Interestingly, in the recovered velocity distribution outside the bar region, the mean radial velocity is not zero but exhibits weak yet systematic asymmetric variations still caused by the bar (Monari et al. 2016b). As we have shown already, the spiral arms behaviour is not captured by our orbit superposition, which is very well seen in the residual map (top right panel), depicting the radial velocity pattern outside the bar region. We note that in the case of multi-arm spirals, the radial velocity variations do not trace exactly the density waves.

The main features of the azimuthal velocity distribution (second row in Fig. 9) are also precisely reproduced by our model. The bar influences the oval-like shape of the rising inner part of the rotational velocity. Once again, we observe some discrepancies outside the bar region, attributable to the presence of spiral arms, as discussed in the context of the radial velocity distributions. The residual map displays certain features, which we attribute to the precision of the stellar density distribution approximation in the bar region. Although we have demonstrated that the approximated density effectively recovers the snapshot data, minor deviations (see Fig. 5) may lead to the observed small-scale discrepancies in the kinematic data, which, however, for the azimuthal velocity do not exceed 10 km s−1.

Unlike the in-plane velocity components, the reconstructed vertical velocity map exhibits significant discrepancies with the snapshot data (third row in Fig. 9). The velocity distribution based on Mock 1 appears essentially featureless, as expected from the orbits superposition modelling, which provides an equilibrium solution. However, the simulation snapshot displays some spiral-like vertical velocity corrugations. Since these variations trace the radial velocity pattern associated with spirals, we assume that the vertical velocity features mainly result from the presence of spiral arms, consistent with some previous theoretical works (Debattista 2014; Faure et al. 2014; Monari et al. 2016a); however the effect of the X-shaped bulge formation (Łokas 2019; Khoperskov et al. 2019) and vertical perturbations excited by the halo particles (Chequers et al. 2018) might be non-negligible.

Overall, the velocity dispersion components (three bottom rows in Fig. 9) are recovered quite precisely by our Model 1, exhibiting only minor inconsistencies compared to the snapshot data, typically less than 10 km s−1. As discussed earlier, these discrepancies likely stem from a slight mismatch of the density (gravitational potential) approximation and the presence of non-equilibrium kinematic features in the disc, as discussed above.

To summarise the results of this section. We have demonstrated that our orbit superposition approach allows not only the recovery of the adopted 3D stellar density distribution but also the reproduction of the kinematics of stars in great detail across the entire galaxy. The quality of the reconstruction of the kinematic parameters severely depends on the adopted star sample used to calculate the orbital library. However, using mock catalogues reproducing the APOGEE footprint, we demonstrated that the present-day data are already sufficient for such purposes.

thumbnail Fig. 8

Reconstruction of the rotational velocity distribution as a function of the galactocentric distance. The left panel displays the stellar density distribution (normalized by the total stellar mass) in the simulation, while the middle and right panels depict the stellar density reconstructed using orbit superposition based on Mock 1 and Mock 2, respectively. In all panels, the cyan lines represent the mean (solid) and mean ± standard deviation (dashed) of the snapshot data. Correspondingly, the white lines denote the same quantities derived from the orbit superposition. As evident from the middle panes, Model 1 is able to recover not only the mean rotational velocity but also the velocity dispersion as a function of galactocentric distance, including some small details associated with the presence of the bar.
thumbnail Fig. 9

Comparison of the face-on kinematics in the simulation snapshot and orbit superposition Model 1. From top to bottom: the mean velocities VR, Vϕ, Vz, and velocity dispersion components σR, σϕ, and σz. The left column corresponds to the snapshot data, the middle one depicts the orbit superposition, and the right column displays their residuals. The orbit superposition model reproduces kinematic features typical for MW-like barred galaxies, including the quadrupole VR pattern, the axisymmetric rise of the rotational velocity with distance from the centre, Vϕ, and the 2D velocity dispersion profiles with ≈10 km s−1 precision. At the same time, the model is unable to recover the spiral arms-induced velocity perturbations and a weak disequilibrium seen in the Vz maps.
thumbnail Fig. 10

Reconstruction of energy-angular momentum space. The top row displays the initial distribution of stars from the simulation in the E – Lz space in two APOGEE-like mock spatial selections adopted for analysis. The bottom row presents a comparison between the snapshot data (left) and the results of the orbit superposition based on Mock 1 (middle) and Mock 2 (right). Across all panels, the grey-filled area represents the generalised distributions of E and Lz in the snapshot data, contrasted with Mock 1 (blue) and Mock 2 (red) in the middle, the results of the models are at the bottom. The density distribution maps are normalized by the maximum value in a given panel. Both initial mocks are biased towards rotationally supported disc components, as evidenced by their density distributions peaking around the solar radius (see Fig. 3), thereby diminishing contributions from the innermost galaxy. However, the superposition models in the bottom row effectively correct such selection function biases, accurately reproducing the snapshot data.

3.4 Energy, angular momentum, and circularity reconstruction

Although we have already explored the quality of the velocity and velocity dispersion components reconstruction by the orbit superposition method, in this section, we aim to address the behaviour of several other parameters often used to explore the MW and external galaxies. In particular, we look at the energy-angular momentum space (E – Lz) and circularity distributions. We emphasise that, unlike axisymmetric potentials, these quantities are not conserved along orbits in a steadily rotating barred potential (Ceverino & Klypin 2007; Binney & Tremaine 2008). Hence, a single orbit can contribute with different energy values, angular momentum and circularity. However, it is important to recall that the weight attached to different values of kinematic parameters along a single orbit remains constant throughout.

In Fig. 10 we show the stellar density distribution in the E − Lz coordinates. The top panels present the distributions of the input Mock 1 and Mock 2 star particles selection, which do not cover the whole disc (see the spatial footprints in Fig. 3) but are somewhat weighted to the solar radius, which is evident from the panels, showing the dominance of kinematically cold discy component. These APOGEE-like footprints capture relatively few low-energy and low-angular momentum stars, which can be found throughout the entire disc but may not necessarily represent the genuine innermost region of the galaxy. Interestingly, both mock catalogues (in Mock 2, it is more prominent, top right panel) reveal a blob at high energy values and weak-to-no rotation separated from the low-energy tail by a gap. As it is known in the literature (Lane et al. 2022; Belokurov & Kravtsov 2022), and highlighted by Fig. 10, this E – Lz distribution structure is solely caused by the APOGEE footprint. We recall that our simulated galaxy evolved in isolation and, thus, does not include accreted populations and is obviously not affected by external interactions, highlighting the importance of critical analysis of the spatial selection function impact on the E – Lz space (and other kinematic or actions coordinates) often used for the kinematic substructures identification and merger debris analysis in the MW halo (Jean-Baptiste et al. 2017; Khoperskov et al. 2023a; Pagnini et al. 2023; Mori et al. 2024).

The snapshot data in the E – Lz plane (bottom left panel in Fig. 10) reveals at least one striking feature associated with the presence of the bar. Several horizontally aligned overdensities are visible, representing the resonances of the bar in the stellar disc (Dillamore et al. 2023, 2024). The signal observed in our disc-galaxy simulation is much stronger compared to the one found by Dillamore et al. (2023); however, one can expect some contamination from disc stars in the halo. This is partly because these two components have no sharp boundaries, especially if the inner MW halo is composed predominantly of heated disc stars (Di Matteo et al. 2019; Belokurov et al. 2020).

Once we compared the reconstructed stellar density distribution with the snapshot data, we observed that Model 1 reproduces very well the distribution of E, Lz separately, as well as the detailed structure of the E – Lz space. The model captures nearly constant energy overdensities, highlighting their bar-related origin. The Model 2 also closely resembles the true distributions, recovering the Lz distribution well but exhibiting less accuracy in the E distribution. The latter displays several peaks and multiple artificial features in E – Lz coordinates, likely resulting from a limited and possibly not fully representative orbital library. In this situation, the orbit superposition adds too much weight to certain existing orbit families while many important ones are missing, as we already discussed in this and previous sections. Therefore, caution should be exercised when interpreting features in the E – Lz analysis based on the orbit superposition models, as the outcome of such models should be thoughtfully tested depending on the input sample of stars reflecting the survey selection function.

Another parameter widely used to explore both structural and kinematic structure of galaxies is orbital circularity, the ratio of stellar specific angular momentum and the angular momentum of a perfectly circular orbit with the same energy (Abadi et al. 2003; Brook et al. 2004; Tissera et al. 2012), which also plays a vital role in the interpretation of orbit superposition models of external galaxies (Jin et al. 2020; Zhu et al. 2018b,a). In Fig. 11, we depict the circularity distribution as a function of galactocentric distance for the mock selections of star particles (top middle and right), snapshot data (bottom left), and reconstructed distributions in Model 1 and Model 2 (bottom middle and right). The top left panel displays the circularity distribution in the snapshot data (grey) and the distributions obtained using orbit superposition in Model 1 and Model 2, depicted in blue and red, respectively. As we showed earlier, here, Model 1 works perfectly and recovers not only the overall distribution of circularity but also its variation with the galactocentric distance. The fraction of hot orbits dominating the inner galaxy and the transition to colder disc-like populations are very much identical to the snapshot data. Conversely, Model 2 exhibits an artificial overabundance of hot orbit stars with circularity values near zero. This discrepancy was previously observed and discussed in detail above (see Figs. 8 and 10).

thumbnail Fig. 11

Reconstruction of circularity distribution as a function of the galactocentric distance. Similar to Fig. 10, the panels illustrate the circularity distribution as a function of galactocentric distance in the mock selections (top) and the results of orbit superposition (bottom), compared to the snapshot data (bottom left). The top left panel displays the circularity distribution obtained in the orbit superposition (blue and red for Model 1 and 2, respectively) compared to the snapshot data (grey).
thumbnail Fig. 12

Left panel: [Fe/H] distribution for stars in the initial model selections (blue and red). Right panel: results of the orbital decomposition. In both panels, the MDF of stars in the simulation snapshot is represented by the grey-filled histogram. The MDFs of the initial selections are biased against relatively metal-rich stars, which is corrected by the orbit superposition modelling. We emphasise that our approach enables the recovery of the complex shape of the MDF, characterized by multiple peaks, including a small fraction of extremely metal-rich stars originating from the centre of the simulated galaxy. This underscores the capability of our method to reconstruct any characteristics (abundances and ages) of stellar populations that are evolutionarily linked to the orbital parameters of stars.

3.5 Reconstruction of metallicity distribution

The ultimate stress test for our orbit superposition method, based on the MW-like data, lies in understanding the reconstruction of stellar parameter variations across the galaxy. It is important to note that no stellar parameters (chemical abundances and stellar ages) were used in the modelling process or generation of the mock catalogues. However, given the evolutionary coupling between chemical abundances and orbits of stars (Bensby et al. 2003; Haywood 2008; Holmberg et al. 2009; Bovy et al. 2012; Haywood et al. 2013; Minchev et al. 2014; Ness et al. 2019), we anticipate that our model can also reasonably capture the metallicity structure of the simulated galaxy. We remind the reader that here, we use an isolated galaxy model where the enrichment history and, thus, the metallicity scale may not be representative of the MW.

First, we access the reconstruction of the MDF of the entire galaxy. In Fig. 12 we show the MDF of the simulated snapshot (grey) considering all stars formed in the simulation and located within 15 kpc. The Mock 1 and Mock 2 selections are shown with blue and red histograms in the left panels, revealing a bias towards lower metallicity stars, as is expected from their spatial footprints, which do not capture many relatively metal-rich stars in the inner galaxy. Nevertheless, both initial selections still cover the full range of metallicity, yet the fractional contributions of stars in both samples do not reproduce the true MDF. The full range of metallicity in the initial sample is crucial to recovering the true MDF, again highlighting the importance of the initial sample of stars and the need for the complete coverage of the parameter space. The right panel shows the orbit superposition-based MDFs (Model 1 and Model 2) against the simulation snapshot data. We notice that both models recover the true MDF quite well. They recover all three peaks of the MDF, even the most metal-rich one, which was poorly populated in the corresponding mocks. The figure demonstrates the advantage of the orbit superposition methods, as it allows to rebalance the metallicity (or any other stellar parameters) distribution in a given spectroscopic dataset, thus providing selection-function corrected distributions.

Finally, we compare the 2D face-on maps of the mean stellar metallicity obtained in our orbit superposition models. Figure 13 shows the comparison between the true (simulation snapshot data) metallicity distribution (top left), reconstructed maps based on orbit superposition (top middle and right) and corresponding residuals in the bottom panels. We notice that the metallicity distribution in the snapshot data is quite complex. Again, it might not be typical for the disc galaxies and the MW in particular; however, it allows us to test the ability of our approach to reconstruct the complex abundance distributions. Nevertheless, some elements are consistent with our expectations about the metallicity distribution across discs of barred galaxies. The metallicity peaks in the very compact region representing a metal-rich core likely associated with efficient enrichment in the region fuelled by the gas inflow along the bar (Ratcliffe et al. 2024) and recently explored in the MW (Rix et al. 2024). Inside the bar region, the metallicity is lower with a modest enhancement along the bar (Di Matteo et al. 2013; Fragkoudi et al. 2018; Neumann et al. 2024). Outside the bar region, the metallicity is higher, and it decreases towards the disc outskirts. The spiral arms are seen as having slightly higher metallicity compared to the surrounding area, as expected in various scenarios (Sánchez-Menguiano et al. 2016; Spitoni et al. 2019; Khoperskov et al. 2018, 2023b).

Model 1 captures most of the elements of the metallicity map we have in the simulation: the central metallicity peak, lower metallicity region inside the bar region and a decline of the mean stellar metallicity outside the bar. We also note that there is a prominent large-scale azimuthal variation of the mean metallicity outside the bar, which is not associated with spiral arms and seems to be observed in the MW (Hawkins 2023), as we suggest can be caused by the bar. The spiral arms-related increase of the metallicity is not seen in the reconstructed metallicity map based on Model 1; however, it is prominent in the residuals. Model 2 surprisingly well recovers these features. However, the solution is quite noisy, the metal-rich core is completely missing, and the residuals show more peculiar features. Nevertheless, the overall metallicity map in Model 2 is quite adequate, considering how poorly it recovers the density and kinematic information.

4 Discussion

4.1 Limitations of the orbit superposition approach

While we have demonstrated the robustness of our approach in mapping the galactic disc using the orbit superposition, we need to acknowledge the limitations inherent in the method. The method we describe aims to reconstruct the present-day complete distribution function of the galaxy using patchy information only about the distribution of its resolved stellar populations, and formally, it does not provide or use any information about the galactic structure in the past. The orbital library we construct is obtained in a rotating reference frame, where the time variable is essentially excluded. Although we have created a galaxy model implementing the rotating bar, axisymmetric models can be useful for various applications, such as the halo structure reconstruction, where the impact of the bar and other disc asymmetries are not dominant.

As discussed in the relevant sections, we underline that our model cannot capture the non-equilibrium kinematic features naturally present in galactic discs. Spiral arms, in principle, can be included in the potential reconstruction and orbital library calculation; however, this would imply that spirals are steady structures that maintain their strength and pattern speed on a long time scale. This is quite a very strong requirement and is not supported by the modern literature (Dobbs & Baba 2014; Sellwood & Masters 2022), suggesting that spirals are transient, recurrent and initiated by swing amplified instabilities. A constant speed of the bar is another assumption which might not perfectly represent the situation in barred galaxies. As we discussed above, non-equilibrium potential, such as a rapidly evolving (buckling) bulge or bar slowing down or accelerating, cannot be adequately captured. The latter, however, can be implemented but requires a creative approach for projecting the orbits of stars into the present-day stellar populations allowing to avoid phase-mixing.

It is well-known that the outer disc of the MW is warped, which is seen not only from the star count (Drimmel & Spergel 2001; López-Corredoira et al. 2002) but also in the kinematics of stars (Poggio et al. 2018; Schönrich & Dehnen 2018). Our model is not able to reproduce such a phenomenon unless the approximated gravitational potential includes such a feature. This can be implemented; however, this would require a steady rotation of the warp or an assumption about its evolution in the past, which might not be so trivial taking into account its precession (Poggio et al. 2020; Jónsson & McMillan 2024).

The results of the orbit superposition are governed by the adopted potential of the galaxy, which is trivial to obtain using a simulated galaxy; however, in practice, the true potential of the MW is not known. Therefore, although several models of the MW potential are available (Bovy 2015; McMillan 2017; Pouliasis et al. 2017; Portail et al. 2017a; Sormani et al. 2022), great caution should be addressed to the potential adopted in the context of the real MW data. In this case, the pattern speed and the present-day orientation of the MW bar can play a vital role, too, which is not a major problem in the simulated data. At the same time, a non-reliable potential is likely to result in odd kinematic features compared to the real data, and we suggest that at least the rotational velocity should be well-recovered by the orbit superposition.

In this work, we showed that two different selections of star particles can result in different solutions. This aligns with the results from Schwarzschild modelling of external galaxies, which underscore the significance of large-area coverage and high spatial resolution constraints offered by IFU data compared to the long-slit observations (Cappellari & McDermid 2005). The enhanced capabilities of IFU data make it possible to resolve degeneracies in determining the galaxy’s gravitational potential, including more precise constraints on the parameters of the central SMBHs. To obtain more robust results in the context of the MW, the sample of stars (providing orbits used for superposition) should at least partially include a contribution of stars originating from the inner galaxy and whose apocenters lie inside the initial sample area (mock catalogue). This is particularly important for the reconstruction of the galactic nuclei and the bulge region. Similarly, the orbits of those stars in the outer disc whose pericenters fall outside the selection footprint will not be accounted for in the model, which will be unrealistic at the outskirts.

The selected sample of stars also should cover the whole range of stellar parameters if one is interested in mapping metallicities, elemental abundances and ages. The orbit superposition models developed here do not provide additional information about the stellar populations but rather allow us to extrapolate or project their parameters along their orbital space. Therefore, the model may not adequately represent reality if the abundance or age information in the input sample of stars is incomplete or heavily biased. In any case, the ground truth about the MW stellar populations is unknown as different methods can deliver different abundances from the same stellar spectra (see, e.g. Hayes et al. 2022) and various age determination methods, although they generally agree with each other, may provide different ages for the same stars (see, e.g. Anders et al. 2023; Queiroz et al. 2023). Therefore, the orbit superposition models should be considered as a projection or extension of a given set of parameters to the whole galaxy under certain assumptions regarding its gravitational potential and orbital library.

thumbnail Fig. 13

Reconstruction of the face-on stellar disc metallicity maps. The top panels display the mean [Fe/H] maps for the snapshot data (left) and the orbit superposition modelling based on Model 1 (middle) and Model 2 (right). The bottom panels illustrate the residuals between the snapshot data and the reconstructed metallicity maps. While the models effectively reproduce the overall 2D behaviour of metallicity, they fail to capture the enhancement of the mean metallicity associated with spiral arms.

4.2 Future prospects

It is important to recognize that the limitations discussed in the previous section allow us to investigate processes influenced by factors not included in our model. For instance, we have shown the metallicity and velocity patterns associated with spiral arms are quite prominent in the residual maps (see Figs. 13 and 9). Thus, such systematic variation of stellar parameters is seen once the data are compared with the outcome of the orbit superposition models. This allows us to subtract the background information (orbit superposition results) from the data and focus on specific phenomena, such as spirals and the warp.

The origin of the galactic spiral structure remains largely elusive; hence, testing what we can learn from models with a steady or evolving spiral structure can be very promising. By experimenting with the strength and the pattern speed of spirals, one can find to what extent the steady spiral structure is not applicable for recovering the observed data. Since its discovery (Blitz & Spergel 1991; Nakada et al. 1991; Weinberg 1992), a lot of efforts has been paid to pin down the parameters of the MW bar (Bland-Hawthorn & Gerhard 2016). A similar strategy we propose for spiral arms study can be pursued for finding the bar parameters if one can vary the structure and kinematics of the bar and use some known parameters of stellar populations, kinematics or abundance distribution, to find the most appropriate solution. For instance, one can try to constrain the bar parameters by matching its resonance ‘location’ in various kinematic spaces. In some sense, such an approach would be similar to the made-to-measure (M2M) models (de Lorenzi et al. 2007; Dehnen 2009; Hunt & Kawata 2013; Syer & Tremaine 1996; Long & Mao 2010). However, the main advantage compared to the existing M2M models is that in the orbit superposition approach, one can also use the stellar parameters information, thus putting more constraints on the model and learning more about the present-day structure of the MW stellar populations.

Finally, the implemented orbit superposition model is relatively simple, as it assumes that the total and stellar disc potentials (or density distributions) are known, allowing us to bypass the use of kinematic data. The latter statement, however, may not be entirely precise if the method is applied to the MW because the parameters of its main components (DM, stellar disc and gas density distributions) are, in fact, constrained using the kinematics of various tracers which are not part of the orbit superposition modelling. Nevertheless, despite progress in constraining the 3D density distribution of the MW, as discussed in the paper, some uncertainties remain. Therefore, we propose that future developments of orbit superposition models for the MW should relax the requirement for a known galactic potential. Instead, the potential could be constrained by incorporating kinematic information (higher-order moments of the distribution function) from the initial sample of stars used to create the orbital library. This step is not trivial, as it would require a convolution of the orbit superposition model results with observational errors (uncertainties in radial velocity, distance, and proper motions) and their ‘re-observation’ using the same spatial footprint as the input dataset. This formalism then will be identical to the modern Schwarzschild orbit superposition models currently developed only for external galaxies.

5 Summary

In this paper, we described a novel orbit superposition method in the context of the MW data about its resolved stellar populations. The method was extensively tested on a simulated MW-like barred galaxy and can be outlined in the following steps:

  1. Define the total galactic gravitational potential and the threedimensional stellar mass density distribution component;

  2. Integrate the orbits of stars from the observational sample within the total potential;

  3. Project the integrated orbits and the 3D stellar density onto a grid;

  4. Calculate non-negative weights for each orbit to ensure their superposition in each grid cell corresponds to the adopted stellar mass.

The resulting orbital weights represent the stellar mass-weighed complete distribution function.

The main results of the paper are as follows:

  • To test the impact of the initial sample of stars on the orbit superposition model results, we used two mock catalogues generated from a simulated MW-like barred galaxy. These catalogues were designed to mimic the spatial distribution of giant stars (Mock 1) and red clump stars (Mock 2) in the APOGEE DR17 dataset (see Fig. 3);

  • We demonstrated that orbit superposition Model 1 (based on Mock 1) shows remarkable agreement with the simulated galaxy in terms of the mean velocity and velocity dispersion components, as well as circularity, energy, and angular momentum distributions (see Sections 3.2–3.4);

  • Model 2 fails to reproduce certain details of the velocity distributions accurately, but it introduces several artificial density and velocity features, particularly in the inner region of the galaxy. As a result, the corresponding spatial footprint (Mock 2) is less suitable for applications related to the MW. The divergence in results between the considered models stems from the incomplete orbital families in Model 2. This model input mock data lacks a sufficient number of stars with apocentres in the innermost region of the galaxy, instead including only hot orbit stars that pass through the initial footprint and, as a result, the orbits of these stars receive too much weight and produce an artificial spheroid-like stellar component, not present in the simulation;

  • We showed that the orbit superposition approach allows the recovery of the MDF of the entire galaxy and the face-on mean stellar metallicity distributions. This showcases the ability of our approach to correct observed stellar parameters distributions and kinematic structure for selection functions of Galactic surveys and deliver reliable information about the chemo-kinematic properties of stellar populations even outside the survey spatial footprint;

  • We have shown that the present-day data, in particular, APOGEE DR 17, is already sufficient to produce a highquality orbit superposition analysis of the MW galaxy. However, the number of stars with precise kinematic and stellar parameter measurements is poised for a significant increase with the upcoming release of data from surveys such as 4MOST, SDSS-V, MOONS, and WEAVE. Therefore, the proposed orbit superposition approach and its further extensions can then be used to project these data and reveal the MW stellar populations in great detail.

The newly developed approach potentially allows for the recovery of an unbiased present-day age, chemical abundance, and kinematic structure of the MW galaxy, which can then be used to reconstruct the formation of its main components and their evolutionary connection. The application of the developed method to the real MW data will be presented in follow-up papers.

Acknowledgements

SK thanks Vasily Belokurov for a discussion about the impact of the APOGEE footprint on the distribution function and James Binney for some useful suggestions regarding the application of the method. 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 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 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. 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.

References

  1. Abadi, M. G., Navarro, J. F., Steinmetz, M., & Eke, V. R. 2003, ApJ, 597, 21 [NASA ADS] [CrossRef] [Google Scholar]
  2. Aitken, R. G., & Doolittle, E. 1932, New General Catalogue of Double Stars within 120∘ of the North Pole [Google Scholar]
  3. Albareti, F. D., Allende Prieto, C., Almeida, A., et al. 2017, ApJS, 233, 25 [Google Scholar]
  4. Almeida, A., Anderson, S. F., Argudo-Fernández, M., et al. 2023, ApJS, 267, 44 [NASA ADS] [CrossRef] [Google Scholar]
  5. Anders, F., Khalatyan, A., Queiroz, A. B. A., et al. 2022, A&A, 658, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Anders, F., Gispert, P., Ratcliffe, B., et al. 2023, A&A, 678, A158 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Antoja, T., Helmi, A., Romero-Gómez, M., et al. 2018, Nature, 561, 360 [Google Scholar]
  8. Athanassoula, E. 2002, ApJ, 569, L83 [NASA ADS] [CrossRef] [Google Scholar]
  9. Belokurov, V., & Kravtsov, A. 2022, MNRAS, 514, 689 [NASA ADS] [CrossRef] [Google Scholar]
  10. Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, MNRAS, 478, 611 [Google Scholar]
  11. Belokurov, V., Sanders, J. L., Fattahi, A., et al. 2020, MNRAS, 494, 3880 [Google Scholar]
  12. Bensby, T., Feltzing, S., & Lundström, I. 2003, A&A, 410, 527 [CrossRef] [EDP Sciences] [Google Scholar]
  13. Binney, J. 2012, MNRAS, 426, 1328 [CrossRef] [Google Scholar]
  14. Binney, J., & Tremaine, S. 2008, Galactic Dynamics, 2nd edn. [Google Scholar]
  15. Binney, J., & Vasiliev, E. 2023, MNRAS, 520, 1832 [NASA ADS] [CrossRef] [Google Scholar]
  16. Binney, J., & Vasiliev, E. 2024, MNRAS, 527, 1915 [Google Scholar]
  17. Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529 [Google Scholar]
  18. Blitz, L., & Spergel, D. N. 1991, ApJ, 379, 631 [Google Scholar]
  19. Boubert, D., & Everall, A. 2022, MNRAS, 510, 4626 [NASA ADS] [CrossRef] [Google Scholar]
  20. Boubert, D., Everall, A., Fraser, J., Gration, A., & Holl, B. 2021, MNRAS, 501, 2954 [NASA ADS] [CrossRef] [Google Scholar]
  21. Bovy, J. 2015, ApJS, 216, 29 [NASA ADS] [CrossRef] [Google Scholar]
  22. Bovy, J., Nidever, D. L., Rix, H.-W., et al. 2014, ApJ, 790, 127 [CrossRef] [Google Scholar]
  23. Brook, C. B., Kawata, D., Gibson, B. K., & Freeman, K. C. 2004, ApJ, 612, 894 [NASA ADS] [CrossRef] [Google Scholar]
  24. Bovy, J., Rix, H.-W., Hogg, D. W., et al. 2012, ApJ, 755, 115 [Google Scholar]
  25. Bovy, J., Leung, H. W., Hunt, J. A. S., et al. 2019, MNRAS, 490, 4740 [Google Scholar]
  26. Buder, S., Sharma, S., Kos, J., et al. 2021, MNRAS, 506, 150 [NASA ADS] [CrossRef] [Google Scholar]
  27. Cappellari, M., & McDermid, R. M. 2005, Class. Quant. Grav., 22, S347 [NASA ADS] [Google Scholar]
  28. Cappellari, M., Bacon, R., Bureau, M., et al. 2006, MNRAS, 366, 1126 [Google Scholar]
  29. Castro-Ginard, A., Brown, A. G. A., Kostrzewa-Rutkowska, Z., et al. 2023, A&A, 677, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Cautun, M., Benítez-Llambay, A., Deason, A. J., et al. 2020, MNRAS, 494, 4291 [Google Scholar]
  31. Cerqui, V., Haywood, M., Di Matteo, P., Katz, D., & Royer, F. 2023, A&A, 676, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Ceverino, D., & Klypin, A. 2007, MNRAS, 379, 1155 [NASA ADS] [CrossRef] [Google Scholar]
  33. Chen, B. Q., Liu, X. W., Yuan, H. B., et al. 2018, MNRAS, 476, 3278 [NASA ADS] [CrossRef] [Google Scholar]
  34. Chequers, M. H., Widrow, L. M., & Darling, K. 2018, MNRAS, 480, 4244 [CrossRef] [Google Scholar]
  35. Comerón, S., Trujillo, I., Cappellari, M., et al. 2023, A&A, 675, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Czekaj, M. A., Robin, A. C., Figueras, F., Luri, X., & Haywood, M. 2014, A&A, 564, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. de Jong, R. S., Agertz, O., Berbel, A. A., et al. 2019, The Messenger, 175, 3 [NASA ADS] [Google Scholar]
  38. de Lorenzi, F., Debattista, V. P., Gerhard, O., & Sambhus, N. 2007, MNRAS, 376, 71 [Google Scholar]
  39. de Zeeuw, T. 1985, MNRAS, 216, 273 [NASA ADS] [CrossRef] [Google Scholar]
  40. Debattista, V. P. 2014, MNRAS, 443, L1 [NASA ADS] [CrossRef] [Google Scholar]
  41. Debattista, V. P., & Sellwood, J. A. 1998, ApJ, 493, L5 [Google Scholar]
  42. Debattista, V. P., & Sellwood, J. A. 2000, ApJ, 543, 704 [Google Scholar]
  43. Dehnen, W. 1998, AJ, 115, 2384 [NASA ADS] [CrossRef] [Google Scholar]
  44. Dehnen, W. 2009, MNRAS, 395, 1079 [NASA ADS] [CrossRef] [Google Scholar]
  45. Dehnen, W., & Binney, J. J. 1998, MNRAS, 298, 387 [NASA ADS] [CrossRef] [Google Scholar]
  46. den Brok, M., Krajnović, D., Emsellem, E., Brinchmann, J., & Maseda, M. 2021, MNRAS, 508, 4786 [NASA ADS] [CrossRef] [Google Scholar]
  47. Di Matteo, P., Haywood, M., Combes, F., Semelin, B., & Snaith, O. N. 2013, A&A, 553, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Di Matteo, P., Haywood, M., Lehnert, M. D., et al. 2019, A&A, 632, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Dillamore, A. M., Belokurov, V., Evans, N. W., & Davies, E. Y. 2023, MNRAS, 524, 3596 [NASA ADS] [CrossRef] [Google Scholar]
  50. Dillamore, A. M., Belokurov, V., & Evans, N. W. 2024, MNRAS, 532, 4389 [NASA ADS] [CrossRef] [Google Scholar]
  51. Ding, Y., Zhu, L., van de Ven, G., et al. 2023, A&A, 672, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Dobbs, C., & Baba, J. 2014, PASA, 31, e035 [NASA ADS] [CrossRef] [Google Scholar]
  53. Drimmel, R., & Spergel, D. N. 2001, ApJ, 556, 181 [Google Scholar]
  54. Eilers, A.-C., Hogg, D. W., Rix, H.-W., & Ness, M. K. 2019, ApJ, 871, 120 [Google Scholar]
  55. Eilers, A.-C., Hogg, D. W., Rix, H.-W., et al. 2022, ApJ, 928, 23 [NASA ADS] [CrossRef] [Google Scholar]
  56. Erkal, D., Belokurov, V., Laporte, C. F. P., et al. 2019, MNRAS, 487, 2685 [Google Scholar]
  57. Erkal, D., Belokurov, V. A., & Parkin, D. L. 2020, MNRAS, 498, 5574 [Google Scholar]
  58. Faure, C., Siebert, A., & Famaey, B. 2014, MNRAS, 440, 2564 [CrossRef] [Google Scholar]
  59. Fragkoudi, F., Di Matteo, P., Haywood, M., et al. 2018, A&A, 616, A180 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Fragkoudi, F., Katz, D., Trick, W., et al. 2019, MNRAS, 488, 3324 [NASA ADS] [Google Scholar]
  61. Fragkoudi, F., Grand, R. J. J., Pakmor, R., et al. 2020, MNRAS, 494, 5936 [Google Scholar]
  62. Frankel, N., Sanders, J., Rix, H.-W., Ting, Y.-S., & Ness, M. 2019, ApJ, 884, 99 [Google Scholar]
  63. Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Gaia Collaboration (Drimmel, R., et al.) 2023a, A&A, 674, A37 [CrossRef] [EDP Sciences] [Google Scholar]
  66. Gaia Collaboration (Vallenari, A., et al.) 2023b, A&A, 674, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Gebhardt, K., Richstone, D., Tremaine, S., et al. 2003, ApJ, 583, 92 [Google Scholar]
  68. Gliese, W., & Jahreiß, H. 1979, A&AS, 38, 423 [Google Scholar]
  69. Gonzalez, O. A., Mucciarelli, A., Origlia, L., et al. 2020, The Messenger, 180, 18 [NASA ADS] [Google Scholar]
  70. GRAVITY Collaboration (Abuter, R., et al.) 2021, A&A, 647, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Hauck, B., & Mermilliod, M. 1998, A&AS, 129, 431 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Hawkins, K. 2023, MNRAS, 525, 3318 [NASA ADS] [CrossRef] [Google Scholar]
  73. Hayes, C. R., Masseron, T., Sobeck, J., et al. 2022, ApJS, 262, 34 [NASA ADS] [CrossRef] [Google Scholar]
  74. Haywood, M. 2008, MNRAS, 388, 1175 [NASA ADS] [CrossRef] [Google Scholar]
  75. 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]
  76. Haywood, M., Di Matteo, P., Lehnert, M. D., et al. 2018, ApJ, 863, 113 [Google Scholar]
  77. Haywood, M., Khoperskov, S., Cerqui, V., et al. 2024, A&A, 690, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Helmi, A. 2020, ARA&A, 58, 205 [Google Scholar]
  79. Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Nature, 563, 85 [Google Scholar]
  80. Herschel, W. 1785, Philos. Trans. Roy. Soc. Lond. Ser. I, 75, 213 [Google Scholar]
  81. Hilmi, T., Minchev, I., Buck, T., et al. 2020, MNRAS, 497, 933 [Google Scholar]
  82. Høg, E., Bässgen, G., Bastian, U., et al. 1997, A&A, 323, L57 [Google Scholar]
  83. Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27 [Google Scholar]
  84. Holmberg, J., Nordström, B., & Andersen, J. 2009, A&A, 501, 941 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Houghton, R. C. W., Magorrian, J., Sarzi, M., et al. 2006, MNRAS, 367, 2 [NASA ADS] [CrossRef] [Google Scholar]
  86. Hunt, J. A. S., & Kawata, D. 2013, MNRAS, 430, 1928 [NASA ADS] [CrossRef] [Google Scholar]
  87. Hunt, J. A. S., Bub, M. W., Bovy, J., et al. 2019, MNRAS, 490, 1026 [NASA ADS] [CrossRef] [Google Scholar]
  88. Imig, J., Price, C., Holtzman, J. A., et al. 2023, ApJ, 954, 124 [CrossRef] [Google Scholar]
  89. Jean-Baptiste, I., Di Matteo, P., Haywood, M., et al. 2017, A&A, 604, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Jethwa, P., Thater, S., Maindl, T., & Van de Ven, G. 2020, DYNAMITE: DYnamics, Age and Metallicity Indicators Tracing Evolution, Astrophysics Source Code Library [record ascl:2011.007] [Google Scholar]
  91. Jin, Y., Zhu, L., Long, R. J., et al. 2020, MNRAS, 491, 1690 [NASA ADS] [Google Scholar]
  92. Jin, S., Trager, S. C., Dalton, G. B., et al. 2024, MNRAS, 530, 2688 [NASA ADS] [CrossRef] [Google Scholar]
  93. Jónsson, V. H., & McMillan, P. J. 2024, A&A, 688, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Kacharov, N., Tahmasebzadeh, B., Cioni, M.-R. L., et al. 2024, A&A, 692, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Kapteyn, J. C., & van Rhijn, P. J. 1920, ApJ, 52, 23 [NASA ADS] [Google Scholar]
  96. Khoperskov, S., & Gerhard, O. 2022, A&A, 663, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Khoperskov, S., Di Matteo, P., Haywood, M., & Combes, F. 2018, A&A, 611, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Khoperskov, S., Di Matteo, P., Gerhard, O., et al. 2019, A&A, 622, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  99. Khoperskov, S., Gerhard, O., Di Matteo, P., et al. 2020, A&A, 634, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  100. Khoperskov, S., Haywood, M., Snaith, O., et al. 2021, MNRAS, 501, 5176 [NASA ADS] [CrossRef] [Google Scholar]
  101. Khoperskov, S., Minchev, I., Libeskind, N., et al. 2023a, A&A, 677, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  102. Khoperskov, S., Sivkova, E., Saburova, A., et al. 2023b, A&A, 671, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  103. Khoperskov, S., Di Matteo, P., Steinmetz, M., et al. 2024a, A&A, submitted [arXiv:2411.18182] [Google Scholar]
  104. Khoperskov, S., Steinmetz, M., Haywood, M., et al. 2024b, A&A, submitted [arXiv:2411.16866] [Google Scholar]
  105. Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511 [Google Scholar]
  106. Krajnović, D., Cappellari, M., Emsellem, E., McDermid, R. M., & de Zeeuw, P. T. 2005, MNRAS, 357, 1113 [CrossRef] [Google Scholar]
  107. Krajnović, D., McDermid, R. M., Cappellari, M., & Davies, R. L. 2009, MNRAS, 399, 1839 [Google Scholar]
  108. Kruijssen, J. M. D., Pfeffer, J. L., Reina-Campos, M., Crain, R. A., & Bastian, N. 2019, MNRAS, 486, 3180 [Google Scholar]
  109. Lane, J. M. M., Bovy, J., & Mackereth, J. T. 2022, MNRAS, 510, 5119 [NASA ADS] [CrossRef] [Google Scholar]
  110. Leung, H. W., Bovy, J., Mackereth, J. T., & Miglio, A. 2023, MNRAS, 522, 4577 [NASA ADS] [CrossRef] [Google Scholar]
  111. Łokas, E. L. 2019, A&A, 629, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  112. Long, R. J., & Mao, S. 2010, MNRAS, 405, 301 [NASA ADS] [Google Scholar]
  113. López-Corredoira, M., Cabrera-Lavers, A., Garzón, F., & Hammersley, P. L. 2002, A&A, 394, 883 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  114. Lu, S., Zhu, K., Cappellari, M., et al. 2024, MNRAS, 530, 4474 [NASA ADS] [CrossRef] [Google Scholar]
  115. Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, AJ, 154, 94 [NASA ADS] [CrossRef] [Google Scholar]
  116. Marques, L., Minchev, I., Ratcliffe, B., et al. 2025, A&A, submitted [arXiv:2502.02651] [Google Scholar]
  117. McMillan, P. J. 2017, MNRAS, 465, 76 [NASA ADS] [CrossRef] [Google Scholar]
  118. Minchev, I., Chiappini, C., Martig, M., et al. 2014, ApJ, 781, L20 [Google Scholar]
  119. Mints, A., & Hekker, S. 2019, A&A, 621, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  120. Monari, G., Famaey, B., & Siebert, A. 2016a, MNRAS, 457, 2569 [Google Scholar]
  121. Monari, G., Famaey, B., Siebert, A., et al. 2016b, MNRAS, 461, 3835 [Google Scholar]
  122. Mori, A., Di Matteo, P., Salvadori, S., et al. 2024, A&A, 690, A136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  123. Nakada, Y., Onaka, T., Yamamura, I., et al. 1991, Nature, 353, 140 [NASA ADS] [CrossRef] [Google Scholar]
  124. Nandakumar, G., Schultheis, M., Hayden, M., et al. 2017, A&A, 606, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  125. Ness, M. K., Johnston, K. V., Blancato, K., et al. 2019, ApJ, 883, 177 [CrossRef] [Google Scholar]
  126. Neumann, J., Thomas, D., Maraston, C., et al. 2024, MNRAS, 534, 2438 [NASA ADS] [Google Scholar]
  127. Neureiter, B., Thomas, J., Saglia, R., et al. 2021, MNRAS, 500, 1437 [Google Scholar]
  128. Nidever, D. L., Bovy, J., Bird, J. C., et al. 2014, ApJ, 796, 38 [Google Scholar]
  129. Paczyński, B., & Stanek, K. Z. 1998, ApJ, 494, L219 [Google Scholar]
  130. Pagnini, G., Di Matteo, P., Khoperskov, S., et al. 2023, A&A, 673, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  131. Perryman, M. 2009, Astronomical Applications of Astrometry: Ten Years of Exploitation of the Hipparcos Satellite Data [Google Scholar]
  132. Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, A&A, 323, L49 [Google Scholar]
  133. Petersen, M. S., Weinberg, M. D., & Katz, N. 2021, MNRAS, 500, 838 [Google Scholar]
  134. Pfenniger, D. 1984, A&A, 141, 171 [Google Scholar]
  135. Poci, A., McDermid, R. M., Lyubenova, M., et al. 2021, A&A, 647, A145 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  136. Poggio, E., Drimmel, R., Lattanzi, M. G., et al. 2018, MNRAS, 481, L21 [Google Scholar]
  137. Poggio, E., Drimmel, R., Andrae, R., et al. 2020, Nat. Astron., 4, 590 [NASA ADS] [CrossRef] [Google Scholar]
  138. Portail, M., Gerhard, O., Wegg, C., & Ness, M. 2017a, MNRAS, 465, 1621 [NASA ADS] [CrossRef] [Google Scholar]
  139. Portail, M., Wegg, C., Gerhard, O., & Ness, M. 2017b, MNRAS, 470, 1233 [NASA ADS] [CrossRef] [Google Scholar]
  140. Posti, L., & Helmi, A. 2019, A&A, 621, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  141. Pouliasis, E., Di Matteo, P., & Haywood, M. 2017, A&A, 598, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  142. Prantzos, N., Abia, C., Chen, T., et al. 2023, MNRAS, 523, 2126 [NASA ADS] [CrossRef] [Google Scholar]
  143. Price-Whelan, A. M., Hogg, D. W., Johnston, K. V., et al. 2021, ApJ, 910, 17 [NASA ADS] [CrossRef] [Google Scholar]
  144. Queiroz, A. B. A., Anders, F., Chiappini, C., et al. 2023, A&A, 673, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  145. Randich, S., Gilmore, G., Magrini, L., et al. 2022, A&A, 666, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  146. Ratcliffe, B., Khoperskov, S., Minchev, I., et al. 2024, A&A, 690, A352 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  147. Richstone, D. O., & Tremaine, S. 1988, ApJ, 327, 82 [NASA ADS] [Google Scholar]
  148. Rix, H.-W., & Bovy, J. 2013, A&A Rev., 21, 61 [NASA ADS] [CrossRef] [Google Scholar]
  149. Rix, H.-W., de Zeeuw, P. T., Cretton, N., van der Marel, R. P., & Carollo, C. M. 1997, ApJ, 488, 702 [NASA ADS] [CrossRef] [Google Scholar]
  150. Rix, H.-W., Chandra, V., Zasowski, G., et al. 2024, ApJ, 975, 293 [NASA ADS] [CrossRef] [Google Scholar]
  151. Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  152. Saglia, R. P., Opitsch, M., Erwin, P., et al. 2016, ApJ, 818, 47 [Google Scholar]
  153. Sánchez-Menguiano, L., Sánchez, S. F., Kawata, D., et al. 2016, ApJ, 830, L40 [CrossRef] [Google Scholar]
  154. Schlesinger, K. J., Johnson, J. A., Rockosi, C. M., et al. 2012, ApJ, 761, 160 [Google Scholar]
  155. Schönrich, R., & Binney, J. 2009, MNRAS, 396, 203 [Google Scholar]
  156. Schönrich, R., & Dehnen, W. 2018, MNRAS, 478, 3809 [Google Scholar]
  157. Schwarzschild, M. 1979, ApJ, 232, 236 [NASA ADS] [CrossRef] [Google Scholar]
  158. Schwarzschild, M. 1993, ApJ, 409, 563 [NASA ADS] [CrossRef] [Google Scholar]
  159. Sellwood, J. A., & Masters, K. L. 2022, ARA&A, 60, 73 [NASA ADS] [CrossRef] [Google Scholar]
  160. Shapley, H. 1918, ApJ, 48, 154 [NASA ADS] [CrossRef] [Google Scholar]
  161. Sharma, S., Bland-Hawthorn, J., Binney, J., et al. 2014, ApJ, 793, 51 [NASA ADS] [CrossRef] [Google Scholar]
  162. Sharma, S., Hayden, M. R., & Bland-Hawthorn, J. 2021, MNRAS, 507, 5882 [NASA ADS] [CrossRef] [Google Scholar]
  163. Sormani, M. C., Gerhard, O., Portail, M., Vasiliev, E., & Clarke, J. 2022, MNRAS, 514, L1 [NASA ADS] [CrossRef] [Google Scholar]
  164. Spitoni, E., Cescutti, G., Minchev, I., et al. 2019, A&A, 628, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  165. Statler, T. S. 1987, ApJ, 321, 113 [NASA ADS] [Google Scholar]
  166. Steinmetz, M., Guiglion, G., McMillan, P. J., et al. 2020, AJ, 160, 83 [NASA ADS] [CrossRef] [Google Scholar]
  167. Syer, D., & Tremaine, S. 1996, MNRAS, 282, 223 [Google Scholar]
  168. Tahmasebzadeh, B., Zhu, L., Shen, J., Gerhard, O., & van de Ven, G. 2022, ApJ, 941, 109 [NASA ADS] [CrossRef] [Google Scholar]
  169. Tahmasebzadeh, B., Zhu, L., Shen, J., et al. 2024, MNRAS, 534, 861 [NASA ADS] [CrossRef] [Google Scholar]
  170. Thater, S., Jethwa, P., Tahmasebzadeh, B., et al. 2022, A&A, 667, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  171. Thomas, J., Saglia, R. P., Bender, R., et al. 2004, MNRAS, 353, 391 [NASA ADS] [CrossRef] [Google Scholar]
  172. Thomas, J., Saglia, R. P., Bender, R., et al. 2011, MNRAS, 415, 545 [Google Scholar]
  173. Thomas, G. F., Battaglia, G., Gran, F., et al. 2024, A&A, 690, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  174. Tissera, P. B., White, S. D. M., & Scannapieco, C. 2012, MNRAS, 420, 255 [NASA ADS] [CrossRef] [Google Scholar]
  175. Trick, W. H., Bovy, J., & Rix, H.-W. 2016, ApJ, 830, 97 [Google Scholar]
  176. Valluri, M., Merritt, D., & Emsellem, E. 2004, ApJ, 602, 66 [Google Scholar]
  177. van Altena, W. F., Lee, J. T., & Hoffleit, E. D. 1995, The general catalogue of trigonometric [stellar] parallaxes [Google Scholar]
  178. van de Ven, G., Falcón-Barroso, J., McDermid, R. M., et al. 2010, ApJ, 719, 1481 [Google Scholar]
  179. van den Bosch, R. C. E. 2016, ApJ, 831, 134 [NASA ADS] [CrossRef] [Google Scholar]
  180. van den Bosch, R. C. E., van de Ven, G., Verolme, E. K., Cappellari, M., & de Zeeuw, P. T. 2008, MNRAS, 385, 647 [Google Scholar]
  181. van der Marel, R. P., & van den Bosch, F. C. 1998, AJ, 116, 2220 [NASA ADS] [CrossRef] [Google Scholar]
  182. van der Marel, R. P., Cretton, N., de Zeeuw, P. T., & Rix, H.-W. 1998, ApJ, 493, 613 [Google Scholar]
  183. Vasiliev, E. 2019, MNRAS, 482, 1525 [Google Scholar]
  184. Vasiliev, E., & Athanassoula, E. 2015, MNRAS, 450, 2842 [Google Scholar]
  185. Vasiliev, E., & Valluri, M. 2020a, ApJ, 889, 39 [Google Scholar]
  186. Vasiliev, E., & Valluri, M. 2020b, in Galactic Dynamics in the Era of Large Surveys, 353, eds. M. Valluri & J. A. Sellwood, 176 [NASA ADS] [Google Scholar]
  187. Vasiliev, E., Belokurov, V., & Erkal, D. 2021, MNRAS, 501, 2279 [NASA ADS] [CrossRef] [Google Scholar]
  188. Vera-Ciro, C., & Helmi, A. 2013, ApJ, 773, L4 [Google Scholar]
  189. Verolme, E. K., Cappellari, M., Copin, Y., et al. 2002, MNRAS, 335, 517 [NASA ADS] [CrossRef] [Google Scholar]
  190. Vislosky, E., Minchev, I., Khoperskov, S., et al. 2024, MNRAS, 528, 3576 [NASA ADS] [CrossRef] [Google Scholar]
  191. Wang, Y., Mao, S., Long, R. J., & Shen, J. 2013, MNRAS, 435, 3437 [CrossRef] [Google Scholar]
  192. Wang, Z., Sharma, S., Hayden, M. R., et al. 2024, MNRAS, 534, 1175 [NASA ADS] [CrossRef] [Google Scholar]
  193. Wegg, C., & Gerhard, O. 2013, MNRAS, 435, 1874 [Google Scholar]
  194. Wegg, C., Gerhard, O., & Portail, M. 2015, MNRAS, 450, 4050 [NASA ADS] [CrossRef] [Google Scholar]
  195. Weinberg, M. D. 1992, ApJ, 384, 81 [NASA ADS] [CrossRef] [Google Scholar]
  196. Wojno, J., Kordopatis, G., Piffl, T., et al. 2017, MNRAS, 468, 3368 [NASA ADS] [CrossRef] [Google Scholar]
  197. Wylie, S. M., Clarke, J. P., & Gerhard, O. E. 2022, A&A, 659, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  198. Zhao, H. 1996, MNRAS, 283, 149 [NASA ADS] [CrossRef] [Google Scholar]
  199. Zhao, G., Zhao, Y.-H., Chu, Y.-Q., Jing, Y.-P., & Deng, L.-C. 2012, Res. Astron. Astrophys., 12, 723 [CrossRef] [Google Scholar]
  200. Zhu, L., van de Ven, G., van den Bosch, R., et al. 2018a, Nat. Astron., 2, 233 [Google Scholar]
  201. Zhu, L., van den Bosch, R., van de Ven, G., et al. 2018b, MNRAS, 473, 3000 [Google Scholar]
  202. Zhu, L., van de Ven, G., Leaman, R., et al. 2020, MNRAS, 496, 1579 [Google Scholar]
  203. Zhu, L., Pillepich, A., van de Ven, G., et al. 2022, A&A, 660, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
1

The term ‘terra incognita’ employed in ancient cartography denotes regions that have not been mapped or documented. It is believed to have first appeared in Ptolemy’s Geography around 150 AD.

All Figures

thumbnail Fig. 1

Example of two orbits (red and blue) of disc stars in a barred potential in inertial (left) and rotating with the bar (right) reference frames. In the bottom panels, the orbits are colour-coded to indicate the progression of time. The ovals in the figure represent the bar orientation, with three different configurations colour-coded to denote changes in its orientation over time in the inertial frame. Notably, the orientation of the bar varies with time in the inertial reference frame while, by definition, it remains fixed in the rotating frame. Therefore, the orbits in the rotating frame allow us to exclude the time dependence of particle positions along the orbits and prevent phase-mixing, enabling us to use them in the orbit superposition approach.
In the text
thumbnail Fig. 2

Approximation of the stellar disc density from the simulation. The surface density of star particles in a snapshot is shown in the left panel. In the middle panel, we present the surface density of the stellar disk obtained from the potential approximated using multipole expansion (see Section 2.2 for details). The rightmost panel shows the difference between the two maps, divided by the snapshot stellar density. Notably, the approximated density (and potential) reproduces well the large-scale structure of the simulated galaxy, including the bar 3D structure; however, by construction, small-scale asymmetries such as spiral arms are not captured. The black circle shows the solar radius of 8.2 kpc. The bar orientation is chosen to be 27∘ relative to the horizontal axis, reproducing the MW bar positioning.
In the text
thumbnail Fig. 3

Initial selection of star particles from the simulation mimicking the APOGEE footprint. In the left panel, the surface density of all star particles in the snapshot is shown. The middle and right panels show two mock catalogues which qualitatively reproduce the APOGEE DR17 giants (Mock 1) and red clump (Mock 2) stars footprints, respectively (see Sect. 2.3 for details). The number of star particles in the snapshot is 7.343.898, while the APOGEE-like selections include 40.354 (Mock 1) and 39.596 (Mock 2). In all the panels, the bar orientation is chosen to be 27∘ relative to the horizontal axis, reproducing the MW bar positioning.
In the text
thumbnail Fig. 4

Analysis of weights of orbits across various orbit superposition model realizations using Mock 1. The five left panels of the top row illustrate the weight distribution in five realizations of the orbit superposition using different 10% randomly sampled orbits, while the rightmost column shows the distribution of mean weights averaged across these realizations. In the second row, the five left panels reveal the relationship between mean weights and those obtained in different realizations of orbit superposition. The rightmost panel in this row presents the ratio between the standard deviation of weights and their mean. The bottom panels depict the relative difference between the stellar surface density obtained in various realizations and the surface density obtained using the mean orbital weights. Although a certain level of stochasticity in the distribution of the weights is seen, all the realizations successfully recover the same density distribution, which is nearly identical to the density based on the mean weights. The observed discrepancy in the orbital weight distribution across different realizations is likely attributed to a non-orthogonal orbital library. We, therefore, recommend using weights averaged across several orbit superposition realizations for more stable results.
In the text
thumbnail Fig. 5

Results of the orbit superposition method. The left panels display the stellar density of orbits with constant weights, essentially presenting the stacked library of orbits of all star particles from Mock 1 (top) and Mock 2 (bottom). The middle panels depict the stellar density reconstructed using weighted orbit superposition where the weights are calculated as the mean across ten orbit superposition realizations (see details in Sect. 3.1). The rightmost panels illustrate the relative residuals between the orbit superposition reconstruction and gravitational potential-based approximation (see Fig. 2 and Sect. 2.2 for more details).
In the text
thumbnail Fig. 6

Mean weights distribution in two orbit superposition models obtained in ten realizations for the APOGEE-like Mock 1 (blue) and Mock 2 (blue; see details in Sect. 2.3). The weight values are given in relative units. The broad distribution of weights in Model 2 indicates a worse convergence of the solution compared to Model 1.
In the text
thumbnail Fig. 7

Distribution of velocity components (VR in left, Vϕ in the middle, and Vz in the right) in the galactic rest frame within three 2 kpc-wide annuli centred at 1 (top), 8 (middle), and 11 kpc (bottom). The snapshot data are represented by the grey-filled areas, while the results of the orbit superposition based on Mock 1 and Mock 2 are depicted in blue and red colours, respectively. The blue and red-filled areas correspond to the initial mock data. For each panel, the numbers correspond to the mean and standard deviation of the velocity distribution in km s−1. We emphasise that the orbit superposition method we employ aims to recover the approximated stellar density distribution; the figure shows that Model 1 reproduces the kinematic information very well. At the same time, this approach is unable to capture spiral arms or the non-equilibrium kinematic features of the disc and bar. Furthermore, the effectiveness of the orbit superposition technique is contingent upon the accuracy of the potential approximation.
In the text
thumbnail Fig. 8

Reconstruction of the rotational velocity distribution as a function of the galactocentric distance. The left panel displays the stellar density distribution (normalized by the total stellar mass) in the simulation, while the middle and right panels depict the stellar density reconstructed using orbit superposition based on Mock 1 and Mock 2, respectively. In all panels, the cyan lines represent the mean (solid) and mean ± standard deviation (dashed) of the snapshot data. Correspondingly, the white lines denote the same quantities derived from the orbit superposition. As evident from the middle panes, Model 1 is able to recover not only the mean rotational velocity but also the velocity dispersion as a function of galactocentric distance, including some small details associated with the presence of the bar.
In the text
thumbnail Fig. 9

Comparison of the face-on kinematics in the simulation snapshot and orbit superposition Model 1. From top to bottom: the mean velocities VR, Vϕ, Vz, and velocity dispersion components σR, σϕ, and σz. The left column corresponds to the snapshot data, the middle one depicts the orbit superposition, and the right column displays their residuals. The orbit superposition model reproduces kinematic features typical for MW-like barred galaxies, including the quadrupole VR pattern, the axisymmetric rise of the rotational velocity with distance from the centre, Vϕ, and the 2D velocity dispersion profiles with ≈10 km s−1 precision. At the same time, the model is unable to recover the spiral arms-induced velocity perturbations and a weak disequilibrium seen in the Vz maps.
In the text
thumbnail Fig. 10

Reconstruction of energy-angular momentum space. The top row displays the initial distribution of stars from the simulation in the E – Lz space in two APOGEE-like mock spatial selections adopted for analysis. The bottom row presents a comparison between the snapshot data (left) and the results of the orbit superposition based on Mock 1 (middle) and Mock 2 (right). Across all panels, the grey-filled area represents the generalised distributions of E and Lz in the snapshot data, contrasted with Mock 1 (blue) and Mock 2 (red) in the middle, the results of the models are at the bottom. The density distribution maps are normalized by the maximum value in a given panel. Both initial mocks are biased towards rotationally supported disc components, as evidenced by their density distributions peaking around the solar radius (see Fig. 3), thereby diminishing contributions from the innermost galaxy. However, the superposition models in the bottom row effectively correct such selection function biases, accurately reproducing the snapshot data.
In the text
thumbnail Fig. 11

Reconstruction of circularity distribution as a function of the galactocentric distance. Similar to Fig. 10, the panels illustrate the circularity distribution as a function of galactocentric distance in the mock selections (top) and the results of orbit superposition (bottom), compared to the snapshot data (bottom left). The top left panel displays the circularity distribution obtained in the orbit superposition (blue and red for Model 1 and 2, respectively) compared to the snapshot data (grey).
In the text
thumbnail Fig. 12

Left panel: [Fe/H] distribution for stars in the initial model selections (blue and red). Right panel: results of the orbital decomposition. In both panels, the MDF of stars in the simulation snapshot is represented by the grey-filled histogram. The MDFs of the initial selections are biased against relatively metal-rich stars, which is corrected by the orbit superposition modelling. We emphasise that our approach enables the recovery of the complex shape of the MDF, characterized by multiple peaks, including a small fraction of extremely metal-rich stars originating from the centre of the simulated galaxy. This underscores the capability of our method to reconstruct any characteristics (abundances and ages) of stellar populations that are evolutionarily linked to the orbital parameters of stars.
In the text
thumbnail Fig. 13

Reconstruction of the face-on stellar disc metallicity maps. The top panels display the mean [Fe/H] maps for the snapshot data (left) and the orbit superposition modelling based on Model 1 (middle) and Model 2 (right). The bottom panels illustrate the residuals between the snapshot data and the reconstructed metallicity maps. While the models effectively reproduce the overall 2D behaviour of metallicity, they fail to capture the enhancement of the mean metallicity associated with spiral arms.
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.