© The Authors 2025

1 Introduction

The Galactic centre region and its gravitational potential have been of great interest for several decades. Before the orbits of the stars around the supermassive black hole Sgr A^⋆ could be measured, which ultimately constrained its mass with unmatched accuracy and precision (e.g. GRAVITY Collaboration 2022), other methods were applied. Extended stellar and gas kinematics were used to constrain the mass distribution, including the mass of Sgr A^⋆, M_•. Similar methods are still used today to derive the mass distributions of star clusters and external galaxies.

The Galactic centre contains two rotating stellar structures: The nuclear star cluster (NSC), with an effective radius R_e ~ 5 pc and a flattening (minor-to-major axis ratio) q ~ 0.7 (Schödel et al. 2014; Fritz et al. 2016; Gallego-Cano et al. 2020), and the surrounding nuclear stellar disc (NSD), with scale length of ~90 pc, a scale height of 28–45 pc, and q ~ 0.35 (Launhardt et al. 2002; Nishiyama et al. 2013; Schödel et al. 2014; Gallego-Cano et al. 2020; Sormani et al. 2022). The NSC dominates the stellar surface density in the projected inner r ~ 7 pc of the Milky Way (Feldmeier-Krause et al. 2025), and the NSD becomes dominant farther out.

Stellar dynamical models of the NSC often made assumptions on its shape and velocity structure. Many studies have assumed a spherical mass distribution and an isotropic velocity distribution (e.g. McGinn et al. 1989; Genzel et al. 1996; Haller et al. 1996; Chakrabarty & Saha 2001; Trippe et al. 2008; Oh et al. 2009; Fritz et al. 2016). Some studies relaxed the assumption and allowed anisotropy (e.g. Genzel et al. 2000; Schödel et al. 2009; Do et al. 2013b; Magorrian 2019) and also a flattening of the cluster (Feldmeier et al. 2014; Chatzopoulos et al. 2015; Feldmeier-Krause et al. 2017b). Many of these studies focussed on data in the inner 1 pc, and some had stellar kinematics data out to a few parsec of the Milky Way. Their goal often was to constrain the mass of Sgr A^⋆, which they were not always able to do, and many underestimated it. A few works focussed on regions on the scale of ~100 pc (Lindqvist et al. 1992; Deguchi et al. 2004; Sormani et al. 2020, 2022; Sanders et al. 2024), where the flatter NSD dominates the gravitational potential. The intermediate region (5–40 pc) lacked extended data so far, and the gravitational potential was inferred via interpolation.

One limitation of many stellar dynamical models is that they require the user to spatially bin the data to a mean velocity and velocity dispersion. By binning data, however, we lose information because we have to combine a sufficiently large number of velocities to obtain a robust estimate of the mean velocity and velocity dispersion. Only some of the aforementioned works on the Galactic centre used discrete data, but it has become more common in recent years (e.g. Chakrabarty & Saha 2001; Do et al. 2013a; Magorrian 2019; Sormani et al. 2022; Sanders et al. 2024). Another disadvantage when the data are binned is that we have to assume that the stars that are considered for binning all belong to the system that is being modelled. This may not be true because the samples can be contaminated with foreground or background stars. Using discrete dynamical models, we can include such a contaminant component and assign each star a probability of being a member star or a contaminant (Watkins et al. 2013).

It is also possible to include knowledge of the stellar populations in discrete models. For example, the stellar metallicity [M/H] can be used to differentiate multiple stellar populations, and their velocity structures can then be constrained separately. This has been done for Galactic globular clusters or dwarf spheroidal galaxies (Zhu et al. 2016a,b; Kamann et al. 2020; Kacharov et al. 2022). This approach is interesting for the Galactic centre because the NSC and NSD have a broad [M/H] distribution, and this indicates a mixture of stellar populations (Do et al. 2015; Feldmeier-Krause et al. 2017a, 2020; Fritz et al. 2021). Most stars are metal rich, but a low-metallicity tail also extends to sub-solar metallicity ([M/H] ≳−1.5 dex). Previous studies suggested that the spatial distribution and/or kinematic properties of the sub-solar [M/H] stars are different from those of the dominating super-solar [M/H] stars (Feldmeier-Krause et al. 2020; Do et al. 2020; Schultheis et al. 2021; Feldmeier-Krause et al. 2025). While they appear to rotate as fast as or even faster than the metal-rich stars in the NSC and inner NSD (Do et al. 2020; Feldmeier-Krause et al. 2025; longitude l ≲ 30 pc), they rotate much slower at larger scales in the NSD (l ≲ 200 pc] and resemble the motion of stars in the Galactic bulge (Schultheis et al. 2021). Some sub-solar [M/H] stars may be remnants from past infall events of star clusters or a dwarf galaxy (Arca Sedda et al. 2020), or they trace close passages of globular clusters (Ishchenko et al. 2023).

We combine recent developments of discrete chemodynamical axisymmetric Jeans models with spectroscopic data presented in recent works (Feldmeier-Krause et al. 2025; Xu et al., in prep.). The data cover the outer part of the NSC and extend out to longitude l=33 pc along the Galactic plane. These data add new information that has been missing in previous models. In addition, we use spectroscopic data from the literature to cover the inner NSC (Feldmeier-Krause et al. 2017a, 2020; Feldmeier-Krause 2022) and match the data with proper motion catalogues (Fritz et al. 2016; Libralato et al. 2021; Smith et al. 2025).

This paper is organised as follows: we describe the data in Sect. 2 and the modelling procedure in Sect. 3. We summarise the results for one-population Jeans models in Sect. 4 and for two-population models in Sect. 5. We discuss our result in Sect. 6 and conclude in Sect. 7.

Fig. 1 Spatial distribution of all stars we used for the model. The black crosses denote stars with VLOS and [M/H], red diamonds show stars with additional proper motions from Fritz et al. (2016), blue circles show stars from Libralato et al. (2020), and orange crosses show stars from Smith et al. (2025).

2 Data set

Our dynamical models require a stellar surface density profile and stellar kinematic data. The stellar surface density profile is used to predict the part of the gravitational potential that is generated by the stars. In addition, we included contributions of the gravitational potential to the supermassive black hole Sgr A^⋆, and we tested how important it is to include a contribution from dark matter. Based on the adopted total gravitational potential, the dynamical model predicts the stellar kinematics. The latter are compared to the stellar kinematic data, which thus constrain the total gravitational potential. For the stellar kinematics, we used a combination of line-of-sight velocities V_LOS and proper motions. To convert the proper motions from mas yr⁻¹ into km s⁻¹, we assumed a galactocentric distance of 8.3 kpc (GRAVITY Collaboration 2022) for all stars. Thus, 1″ corresponds to 0.04 pc. In addition, our data set contained measurements of the overall metallicity [M/H] and photometric colour H − K_S for each star.

2.1 Discrete line-of-sight velocity and metallicity catalogue

We combined several data sets for the discrete stellar line-of-sight (LOS) velocity V_LOS and metallicity [M/H]. We show the spatial distribution of the combined catalogue in Fig. 1 (black pluses). The density is highest in the centre because it contains the most stars. Density variations are also caused by the combination of different data sets that sample stars that are located in various regions and at varying depths.

Feldmeier-Krause et al. (2017a) and Feldmeier-Krause et al. (2020) published measurements of V_LOS and [M/H] of ~700 stars each. The stars are located within a projected distance of 127″ (~5 pc) to Sgr A*. Feldmeier-Krause (2022) published measurements for another ~400 stars, half of which are located at a projected distance of ~20 pc (0.° 14) to the Galactic east and the other half lie west of Sgr A^⋆. Xu et al. (in prep.) provide a catalogue of ~900 stars that are located within a circle with a radius of 247″(~0.° 067 ~9.6 pc) around Sgr A^⋆. All these data were observed with KMOS¹ (Sharples et al. 2013) at the Very Large Telescope (VLT) in the near-infrared K band with a spectral resolution of R~4000. While the depth of the catalogues varies, the applied methods for deriving V_LOS and [M/H] are identical between the various studies. In short, they used a full spectral fitting with the code STARKIT (Do et al. 2015; Kerzendorf & Do 2015) and synthetic model spectra from the PHOENIX spectral library (Husser et al. 2013).

For an even more extended data set, we used the catalogue of Feldmeier-Krause et al. (2025). These data were observed with Flamingos-2 at Gemini-South (Eikenberry et al. 2004), also in the K band, and at a slightly lower spectral resolution (R~3 400). The data extend from Sgr A^⋆ to 33 pc along the Galactic plane, both towards the Galactic east and west. The covered latitude is rather small, only ~1 pc to the Galactic south and 1 pc (2 pc in a central subfield) to the Galactic north. They applied the same method and models to derive V_LOS and [M/H].

The data sets overlap to some extent. The Flamingos-2 data cover the same regions as the various KMOS catalogues, with an overlap of ~60–300 stars, depending on the catalogue. We compared the [M/H] measurements and found that they are mostly consistent. The median [M/H] difference (Δ[M/H]) of all stars with two [M/H] measurements is −0.09 to +0.04 dex (depending on the KMOS data set; see also Appendix B of Feldmeier-Krause et al. 2025). The median statistical uncertainty of all [M/H] is 0.18 dex. The two [M/H] measurements are usually consistent within their statistical uncertainties (i.e. $\mathrm{\Delta }[\mathrm{M}/\mathrm{H}]/\sqrt{(}{\sigma }_1^2+{\sigma }_2^2)\lesssim 1$, where σ₁ and σ₂ denote the [M/H] uncertainties from the two independent measurements). For our final catalogue, we only combined the multiple measurements of a star when the catalogue coordinates agreed within $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}2$, the measurements of V_LOS lay within 30 km s⁻¹, and K_S lay within 0.6 mag. For these stars, we combined the multiple measurements of V_LOS and [M/H] with a simple mean. As uncertainties, we used the propagated error of the mean (i.e. $0.5\cdot \sqrt{(}{\sigma }_1^2+{\sigma }_2^2)$) or the standard deviation of the two measurements, whichever was larger. The V_LOS was corrected for perspective rotation that is caused by the motion of the Sun around the Galactic centre. We computed this correction using the equations of van de Ven et al. (2006), and we found that it varied spatially by ±1 km s⁻¹ at most.

Our final catalogue contains 4791 stars with V_LOS and [M/H] measurements. All the stars are cool, with typical absorption lines that were used for the measurements of V_LOS and [M/H]. The colour-magnitude diagram (CMD) in Fig. 2 shows the stars are on the red giant branch, but brighter than the red clump, which is at K_S ~ 15 mag in the Galactic centre. The stellar photometry comes from the GALACTICNUCLEUS (GNS) survey catalogue (Nogueras-Lara et al. 2019).

The CMD also shows that most stars have H − K_S>1.3 mag, indicating that they are near the Galactic centre. Stars with bluer colours are likely foreground stars. The catalogue contains only a few such stars, and we excluded stars with either H − K_S<1.3 mag or H − K_S > 3.5 mag as foreground and background stars, respectively. Our sample thus contains 4600 unique stars.

Fig. 2 Colour-magnitude diagram H − KS vs. KS of the stars with VLOS and [M/H]. The vertical green lines denote the colour cuts we used to remove likely foreground (H − KS ≤ 1.3 mag) and background stars (H − KS ≥3.5 mag). The colour represents the density of stars.

2.2 Proper motion catalogues

We complemented the LOS data with the proper motion (PM) catalogues of Fritz et al. (2016) for the inner r ≲ 89″(~3.6 pc) and Libralato et al. (2020) for outer regions r ≥187≳(~7.5 pc). The spatial distribution of the PM catalogues where they overlap with the LOS data is illustrated in Fig. 1. The proper motion catalogue of Smith et al. (2025, VIRAC² 2) extends over the entire region of the LOS data, except for the most crowded inner ~1 pc region.

We only used stars of the PM catalogues that also had entries in the LOS data to ensure that the stars were red giants and not hot young stars, because these last can be found across the Galactic centre (e.g. Cotera et al. 1999; Mauerhan et al. 2010b,a; Dong et al. 2011; Clark et al. 2021; Feldmeier-Krause 2022). The spatial distribution (e.g. Do et al. 2013a; Støstad et al. 2015; Feldmeier-Krause et al. 2015; Clark et al. 2021) and kinematics (e.g. von Fellenberg et al. 2022; Hosek et al. 2022) of hot young stars are different from those of the red giant stars we used as tracers of the gravitational potential. While the mass contribution of the young stars is negligible, their kinematics might still be used as independent tracers of the underlying total gravitational potential. The spatial distribution of young stars cannot be measured robustly, however, so that the resulting tracer density might bias the inferred gravitational potential. We therefore chose not to use young stars in our analysis. The Galactic centre contains three young clusters, but only one of them is located in the field of view of our data, in the central parsec. It has a mass of ~10⁴ M_⊙ (e.g. Paumard et al. 2006). The stellar surface density of young stars and red giants in the innermost parsec is comparable (Do et al. 2013a), but the gravitational potential in this region is dominated by Sgr A^⋆. Farther out, the stellar number density of young stars decreases steeply, and they contribute only little overall to the mass at larger radii.

From the Smith et al. (2025) PM catalogue, we excluded stars with a PM exceeding 500 km s⁻¹ at a distance of 8.3 kpc and stars with variable K-band photometry (σ_KS ≥ 0.55 mag) because their PMs are uncertain. Based on the resulting matches with the LOS catalogue (within $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}2$), we excluded stars with an H- or K_S-band discrepancy of more than 0.5 mag as likely mismatches. This left 2688 stars with Smith et al. (2025) PMs. We found 875 matches of the LOS data with Fritz et al. (2016) and 493 matches with Libralato et al. (2020). For the stars with PM measurements in two data sets (≲500 stars), we computed the mean values. We obtained a PM measurement from at least one of these catalogues for 3567 unique stars.

The PM catalogues list the PMs in equatorial coordinates, and we converted them into Galactic coordinates using the equations in Poleski (2013). The Libralato et al. (2020) and Smith et al. (2025) data are in the absolute reference system of the Gaia DR2 catalogue. We corrected their PMs for the motion of Sgr A^⋆ relative to the Sun using the measurements of Reid & Brunthaler (2020) of −6.411 mas yr⁻¹ along the Galactic plane and −0.219 mas yr⁻¹ towards the Galactic north pole. A correction for perspective rotation, as made for the V_LOS data, was unnecessary because the Sun does not move significantly towards or away from Sgr A^⋆.

2.3 Stellar density distribution and multi-Gaussian expansion fit

We used the stellar density map of Gallego-Cano et al. (2020), produced from a combination of HAWK-I³ /VLT (Kissler-Patig et al. 2008) and NACO⁴/VLT (Rousset et al. 2003) K_S-band data. The HAWK-I data were obtained as part of the GNS survey (Nogueras-Lara et al. 2018, 2019) and cover a region of 84.4 pc × 21 pc with a point spread function of $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}2$ full width at half maximum. The central region, where crowding is severe, was complemented with NACO data of ~3.5 pc × 3.5 pc (details on the data reduction in Gallego-Cano et al. 2018). The NACO data have the benefit of a higher spatial resolution ($0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}05$) due to adaptive optics.

The data were cleaned from spectroscopically classified early-type stars (i.e. hot young OB-type stars) using the catalogue of Do et al. (2013a), and the density map was mostly created from red giant stars, but it may contain a few early-type stars that were not yet spectroscopically classified. The fraction of early-type stars drops rapidly beyond the central ~1 pc region of the Galactic centre, however (Støstad et al. 2015; Feldmeier-Krause et al. 2015). More extended spectroscopic studies, including those that provide our discrete catalogue, found only a few isolated early-type stars (Feldmeier-Krause 2022; Feldmeier-Krause et al. 2025, Xu et al., in prep.). We conclude that the contamination of early-type stars in the density map is ≲2% and is therefore negligible.

The stellar photometry was corrected for extinction before only stars in the range 9.0≤K_S,ext≤14.0 mag were selected. This corresponds to a magnitude range ~11.0 ≤ K_S ≤ 16.0 mag. Hence, there is significant overlap with the kinematic tracer stars, which mainly lie at 9.0≤ K_S ≤ 15.0 mag (see Fig. 2). The photometric data are deeper and have a higher completeness than the kinematic tracer stars.

Regions with high extinction can severely affect the density of detected sources, and Gallego-Cano et al. (2020) masked them in the stellar density map. The final map had a pixel scale of 5 arcsec pixel⁻¹ and covered 84.4 pc along the Galactic longitude and 21 pc along the Galactic latitude. It is currently the largest and most complete stellar number density map of the Galactic centre.

The discrete Jeans models use the stellar surface brightness distribution of the tracer stars in the form of a multi-Gaussian expansion (MGE, Emsellem et al. 1994). We derived it using the PYTHON package MGEFIT provided by Cappellari (2002). We measured the counts of the density map in 5-degree-wide sectors centred on Sgr A^⋆. We used our knowledge of the position angle of the major axis (along the Galactic plane), and the flattening (minor- to major-axis ratio) of the NSC (q ~ 0.71, Gallego-Cano et al. 2020). In the MGE fit, we constrained the value q to >0.3 because the flattening of the NSD is larger than ~0.3 (Gallego-Cano et al. 2020; Sormani et al. 2022). We further set the OUTER_SLOPE keyword, which forces the model to a profile at least as steep as R^{−outer_slope}, to 0.5, to allow for a flat profile at large radii. The default range of 1–4 in MGEFIT is too steep for the Galactic centre.

Finally, we converted the output values of the MGE fit into the values required for the Jeans model, following the equations provided in the MGEFIT README. To convert the central surface brightness into units of L_⊙ pc⁻², we assumed an exposure time = 1 s, a zero point = 0 mag, and an extinction A=0 mag, and we adjusted the value of the absolute magnitude for the Sun until our density distribution matched the surface brightness distribution of Feldmeier-Krause et al. (2017b) in the 4.5 μm Spitzer band in the centre. We list our final stellar surface density MGE in Table 1 and show the surface density profiles in the top panel of Fig. 3. The residuals are shown in the middle panel. They are highest in the central ≲2 pc because the density map has some fluctuations caused by the varying extinction. Farther out, the photometry counts are measured in larger regions and combine many pixels, which produces less scatter and lower residuals.

3 Discrete dynamical modelling

We used the discrete Jeans anistotropic MGE (JAM) modelling code by Watkins et al. (2013), CJAM, which is based on equations and the formalism presented by Cappellari (2008). The basic assumptions are that (i) the distribution function of the stellar system satisfies the collisionless Boltzmann equation, (ii) the model is axisymmetrical, and (iii) the velocity ellipsoid is aligned with cylindrical polar coordinates. The model parameters were fitted using EMCEE (Foreman-Mackey et al. 2013), which performed a Markov chain Monte Carlo (MCMC) sampling of the parameter space.

Fig. 3 Upper panel: surface density profile derived from the stellar density maps of Gallego-Cano et al. (2020, blue diamonds). The uncertainties are measured using the corresponding uncertainty map and are only shown for a subset of data points to improve visibility. We matched the profile to the centre value of Feldmeier-Krause et al. (2017b) to convert it into units of L⊙ pc−2. The solid black line denotes the profile along the major axis, the dashed red line shows the profile along the minor axis, and the dot-dashed orange line shows the profile along a position angle of 45°. Middle panel: residuals of the star count data after subtracting the MGE fit. Lower panel: projected axis ratio qproj.

3.1 Gravitational potential

The gravitational potential of the Galactic centre region was modelled by a combination of the supermassive black hole Sgr A^⋆ and the stellar mass distribution. In a subset of models, we also introduced a dark matter (DM) mass distribution.

For Sgr A^⋆, we used the well-measured mass of 4.3 × 10⁶ M_⊙ (GRAVITY Collaboration 2022) and modelled it as a single Gaussian with a width of $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}01$. This spatial scale is much smaller than the projected distance of the stars in our data set to the position of Sgr A^⋆, with the shortest being $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}9$. Even for the closest star, Sgr A^⋆ therefore produced a Keplerian potential in our model. In one set of models, we did not fix M_•, but fitted its value.

For the stellar mass distribution, we used the MGE of Table 1, multiplied by a mass conversion factor. This is equivalent to the mass-to-light ratio 𝚼 used in Jeans models that rely on the stellar surface brightness and not on the stellar number density. Because we scaled our MGE to the surface brightness-based MGE of Feldmeier-Krause et al. (2017b), however, the mass conversion factor is equivalent to a mass-to-light ratio in the 4.5 μm band. We expect it to have a value of ~0.6–1.0 (Norris et al. 2014; Kettlety et al. 2018).

In most models, we assumed a constant value for 𝚼, but we also tried a radially varying value for 𝚼, parameterised as a modified Osipkov-Merritt profile (Osipkov 1979; Merritt 1985) by $$a(r)=a_0+\frac{a_{\mathrm{\infty }}-a_0}{1+{\left(\frac{R_a}{r}\right)}^2},$$(1)

where a₀ is the central value of the quantity (in this case, 𝚼), a_∞ is the asymptotic (outer) value (for r → ∞), and R_a is the turn-over radius.

When 𝚼 (or β_z or κ; see Sect. 3.2) varied in the model, we needed seven different values for 𝚼, that is, a separate value for each of the seven MGE components. As in Kamann et al. (2020), we evaluated Eq. (1) at the radius r at which a given MGE component contributed most to the combined profile. To deproject the surface number density profile, we also required the inclination i. We assumed i=90°, that is, we see the Galactic centre edge-on.

In addition, we included a DM component in a subset of models. We assumed a spherical dark matter distribution (Navarro et al. 1996), $${\rho }_{DM}(r)={\rho }_S{\left(\frac{r}{r_{DM}}\right)}^{-\gamma }{[1+\left(\frac{r}{r_{DM}}\right)]}^{-(3-\gamma )},$$(2)

where ρ_S is the DM scale density, r_DM is the DM scale radius, which we fixed to 8 kpc (i.e. a value far beyond the extent of our data), and γ is the inner DM density slope. When the profile was cored, γ ≈ 0, and when the profile was cuspy, γ ≈ 1. We ran models for the two cases, γ=0 and γ=1, and left only ρ_S as a free parameter. The DM profile was then approximated with a one-dimensional MGE fit.

3.2 Anisotropy β_z and rotation κ

The velocity anisotropy β_z and the rotation parameter κ describe the velocity structure of the stellar system. The anisotropy β_z in our models is a measure of the imbalance of radial motion and vertical motion. The models assume a cylindrically aligned velocity ellipsoid (i.e. $\overline{v_R{v}_z}=0^{}$⁵), and β_z is defined as ${\beta }_z=1-\overline{v_z^2}/\overline{v_R^2}$ in the meridional plane⁶. When β_z=0, there is no imbalance, and the velocity structure of the model is isotropic. For β_z<0, there is a tangential anisotropy, and for β_z>0, there is a radial anisotropy, and the stellar orbits are more radial.

It is possible to assume a radially constant value for β_z or to let it vary for each stellar MGE component. In some models, we parameterised β_z(r) via Eq. (1) (i.e. three parameters: β₀,β_∞, and R_β).

The rotation parameter κ quantifies the amount of rotation of the system. κ is defined as $\kappa =\overline{v_{\varphi }}/(\sqrt{}(\overline{v_{\varphi }^2}-\overline{v_R^2})$, and it is a measure of the ratio of ordered to disordered motion⁷. When κ=0, there is no rotation. An isotropic rotator would be described by β_z=0 and κ = ±1. The sign of κ denotes the direction of the rotation. In the Galactic centre, we expect significant rotation and a negative value for κ (based on the definition of the coordinate system by Watkins et al. 2013). A positive value would mean that the sign of the angular momentum vector was reversed.

As with β, it is possible to implement different κ parametrisations: a constant value of κ (i.e. the same κ for each stellar MGE component, one additional parameter), or κ as a function of radius, parametrised via Eq. (1) (i.e. three parameters: κ₀, κ_∞, and R_κ). This three-parameter model allowed us to change the rotation from the NSC to the NSD.

3.3 Population components

It is possible to build multi-population models with discrete posterior distribution functions, as shown for example by Zhu et al. (2016a,b); Kamann et al. (2020); Kacharov et al. (2022), and we followed a similar procedure. We either included one or two population components and also a background component that represented the Galactic bar. In principle, it is possible to include even more components, but this would require more precise [M/H] data or additional data, such as elemental abundances or reliable colours, which are not available for a sufficiently large sample in the Galactic centre.

We solved the Jeans equations for both components in the same gravitational potential, but with a different tracer density I, β_z, and κ for the two components. We estimated the probability of each star in our kinematic sample to belong to one of the modelled stellar populations (k) or a background component (bg) given its on-sky position, velocity (V_LOS and, if available, PMs), and chemical population tracer ([M/H]). In summary, there were up to three probabilities for each star: the dynamical probability, the spatial probability, and (when k>1) the chemical probability.

3.4 Dynamical probability

Given the gravitational potential, β_z, and κ, CJAM solves the Jeans equations and predicts the first and second velocity moments in three directions j and at the location of each star i in our sample. The first moment is the mean velocity μ, and the second moment is the quadratic sum of the mean velocity μ and the velocity dispersion σ.

Based on the measured velocity V_LOS,i ± δV_LOS,i, and for some stars, also the proper motions V_l,i ± δV_l,i and V_b,i ± δV_b,i, we computed the dynamical probability of star i to belong to population k in one dimension j, $$P_{\mathrm{d}\mathrm{y}\mathrm{n},j,i}^k=\frac{1}{\sqrt{{\left({\sigma }_{j,i}^k\right)}^2+{\left(\delta V_{j,i}\right)}^2}}\cdot \exp [-\frac{1}{2}\frac{{(V_{j,i}-{\mu }_{j,i}^k)}^2}{{\left({\sigma }_{j,i}^k\right)}^2+{\left(\delta V_{j,i}\right)}^2}],$$(3)

where ${\mu }_{j,i}^k$ and ${\sigma }_{j,i}^k$ are the velocity and velocity dispersion predicted by the dynamical model of population k at the sky position of star i, either along the line of sight, parallel to the Galactic plane (l), or perpendicular to it (b). When a star only had $V_{\text{LOS}},P_{\text{dyn},i}^k=P_{\text{dyn},j,i}^k$. When star i also had proper motions, the dynamical probability was $P_{\mathrm{d}\mathrm{y}\mathrm{n},i}^k={\prod }_{j=1}^3{P}_{\mathrm{d}\mathrm{y}\mathrm{n},j,i}^k$.

For the background component, we did not solve the Jeans equations, but assumed a Gaussian velocity distribution with a fixed mean velocity μ^bg=0 km s⁻¹ at the location of each star i and in each direction j and a velocity dispersion of σ^bg=130 km s⁻¹ (Portail et al. 2017) that is consistent with the Galactic bar.

3.5 Spatial probability

We assumed that our MGE (Table 1) represented the total stellar surface density distribution I in the Galactic centre. For two populations, their densities I^k add up to I. The fraction h^k contributed by population k was able to vary as a function of radius, and we parametrised the change in h^k with radius with Eq. (1), that is, with a monotonous three-parameter function. Because we had only two populations, we had I¹(x_i, y_i) = h(r) · I(x_i, y_i) and I²(x_i, y_i) = (1 − h(r)) · I(x_i, y_i). We used the respective I^k to solve the Jeans equations for the two populations separately.

We assumed that the background surface number density is uniform, which is a fair assumption given the small extent of the modelled region in comparison to the size of the bar. We introduced the parameter ε, which denotes the fraction of contaminant stars with respect to the central surface density I(0, 0).

The spatial probability for star i at the position (x_i, y_i) to belong to population k is $$P_{\mathrm{s}\mathrm{p}\mathrm{a},i}^k=\frac{I^k(x_i,y_i)}{I(x_i,y_i)+I_{bg}}=\frac{I^k(x_i,y_i)}{I(x_i,y_i)+ϵ\cdot I(0,0)},$$(4)

and the spatial probability for star i to belong to the contaminant population is simply $$P_{\mathrm{s}\mathrm{p}\mathrm{a},i}^{bg}=1-\sum _k{P}_{\mathrm{s}\mathrm{p}\mathrm{a},i}^k.$$(5)

3.6 Chemical probability

The NSC and NSD may contain stars that belong to different stellar populations with different dynamical properties (β_z, κ). This means that stars may have formed at different times and from different materials. This may be reflected in the stellar metallicity [M/H], and we can use it as a stellar population tracer.

We assumed that each population k has a Gaussian metallicity distribution with a mean metallicity $Z_0^k$ and a metallicity dispersion ${\sigma }_Z^k$. For star i with a measured metallicity Z_i ± δZ_i, the chemical probability for population k becomes $$P_{\mathrm{c}\mathrm{h}\mathrm{m},i}^k=\frac{1}{\sqrt{2\pi [{\left({\sigma }_Z^k\right)}^2+{\left(\delta Z_i\right)}^2)]}}\cdot \exp [-\frac{1}{2}\frac{{(Z_i-Z_0^k)}^2}{{\left({\sigma }_Z^k\right)}^2+{\left(\delta Z_i\right)}^2}].$$(6)

We used the overall metallicity [M/H] and its uncertainty derived from full spectral fitting. For the background component, we used a Gaussian distribution centred on $Z_0^{bg}=0.07$ dex and with a width of ${\sigma }_Z^{bg}=0.3\mathrm{d}\mathrm{e}\mathrm{x}$, as found by Fritz et al. (2021) on a bulge control field observed with KMOS.

3.7 Total probability

The total probability for star i in a given model is then $$L_i=(\sum _{k\ne bg}{P}_{\mathrm{s}\mathrm{p}\mathrm{a},i}^k\cdot P_{\mathrm{c}\mathrm{h}\mathrm{m},i}^k\cdot P_{\mathrm{d}\mathrm{y}\mathrm{n},i}^k)+P_{\mathrm{s}\mathrm{p}\mathrm{a},i}^{bg}\cdot P_{\mathrm{c}\mathrm{h}\mathrm{m},i}^{bg}\cdot P_{\mathrm{d}\mathrm{y}\mathrm{n},i}^{bg},$$(7)

where k denotes the two populations that are not the background. The total likelihood L of a model is given by $$L=\prod _{i=1}^N{L}_i,$$(8)

where i runs over all N stars.

For a given star, its likelihood of belonging to population k is $L_i^k=P_{\mathrm{s}\mathrm{p}\mathrm{a},i}^k\cdot P_{\mathrm{c}\mathrm{h}\mathrm{m},i}^k\cdot P_{\mathrm{d}\mathrm{y}\mathrm{n},i}^k$, and the probability of belonging to a given population is $P_i^k=L_i^k/(L_i^{bg}+{\sum }_k{L}_i^k)$. When we fit only one population, we only have the probability for a single population and the background, and we still lack the chemical probability $(P_{\mathrm{c}\mathrm{h}\mathrm{m},i}^k=P_{\mathrm{c}\mathrm{h}\mathrm{m},i}^{bg}=1)$.

3.8 One-population parameter summary

We ran several sets of models with only one population and the background component, where we only considered the dynamical and spatial probability. These models had the parameters listed below.

• (i)

  The background percentage parameter ϵ. This is one parameter, given in percent, and has a uniform prior [0, 10].

• (ii)

  The anisotropy β_z. This is three parameters, and the uniform priors β_0,∞ are [−1.5, 0.9], log(R_β)=[1.079, 2.91], and R_β in arcseconds.

• (iii)

  The rotation κ. This is three parameters with uniform priors κ_0,∞ of [−1.5, 0.4], log(R_κ)=[1.079, 2.91], and R_κ in arcseconds.

Each run focused on a different part of the gravitational potential, as listed below.

• (i)

  On the mass of the supermassive black hole M_•. This is one parameter, and the uniform prior is [2, 6] × 10⁶ M_⊙.

• (ii)

  On the mass-to-light conversion 𝚼. This is three parameters, and the uniform priors 𝚼_0,∞ are [0.1, 3.0], log(R_𝚼)=[1.079, 2.91], and R_𝚼 in arcseconds.

• (iii)

  On the dark matter scale density ρ_S. This is one parameter, and it is in the range log(ρ_S)=[−5, 5].

We always included a single value for ϵ and a modified Osipkov-Merrit profile (Eq. (1)) for β_z and κ. We also used 𝚼 in each model, but implemented different parametrisations. In Sections 4.1 and 4.3, we use a constant value for 𝚼, that is, we used a single parameter that did not vary with the radius, and we either included M_• (nine free parameters in total) or ρ_S (nine free parameters in total) in the fit. In Sect. 4.2 we fix M_•, neglect DM, and fit a radially varying 𝚼 parametrised with a modified Osipkov-Merrit profile (ten free parameters in total; see also Table 2). The bounds of the scale radii (R_β, R_K, R_𝚼) were chosen to be larger than the innermost MGE width σ and not larger than the outermost discrete kinematic measurement.

We started the chains at β_0,∞=0 (i.e. isotropy), with moderate rotation (κ₀=−0.9, κ_∞=−0.5) with a transition at R_β,κ,𝚼 = 45 arcsec, 𝚼=1, ϵ=5%, M_•=4.3 × 10⁶ M_⊙, and log(ρ_S)=0. For each model set-up, we ran EMCEE (Foreman-Mackey et al. 2013) with 100 walkers and at least 1900 steps, and we discarded the first 20% of the steps before the chains converged. We also tested longer chains for the non-DM models with more than 6000 steps, but our best-fit parameters did not change significantly.

3.9 Two-population parameter summary

In the two-population model, we considered the total probability and a background component. The parameters are listed below.

• (i)

  The background percentage parameter ϵ. This is one parameter, given in percent, and has a uniform prior [0, 10].

• (ii)

  The anisotropy ${\beta }_z^k$. These are two parameters, each with a uniform prior range [−1.5, 0.9] and constant with radius, for k=2 populations.

• (iii)

  The rotation κ^k. These are two parameters, each with a uniform prior range [−1.5, 0.4] and constant with radius, for k=2 populations.

• (iv)

  The mass-to-light conversion 𝚼. This is one parameter with a uniform prior range [0.1, 3.0] and constant with radius. We use the same 𝚼 for both populations k.

• (v)

  The population fraction h. This is three parameters with uniform priors h_0,∞ of [0.5, 1.0], log(R_h)=[1.079, 2.91], and R_h in arcseconds.

• (vi)

  The Gaussian mean metallicity $Z_0^k$ and dispersion ${\sigma }_Z^k$. These are four parameters, and the uniform priors are $Z_0^1=[0.0,0.8]\mathrm{d}\mathrm{e}\mathrm{x},{\sigma }_Z^1=[0.2,1.5]\mathrm{d}\mathrm{e}\mathrm{x},Z_0^2=[-1.6,0.2]\mathrm{d}\mathrm{e}\mathrm{x}$, and ${\sigma }_Z^2=[0.2,1.5]$ dex.

We fit a constant value for ϵ, while the population fraction h was parametrised with Eq. (1). We started the chains at ${\beta }_z^k=0$ (i.e. isotropy), moderate rotation (κ^k=−0.8), $Z_0^1=0.35\mathrm{d}\mathrm{e}\mathrm{x},Z_0^2=-0.6\mathrm{d}\mathrm{e}\mathrm{x},{\sigma }_Z^k=0.4\mathrm{d}\mathrm{e}\mathrm{x},h_{0,\mathrm{\infty }}=0.75,R_h=100\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c},ϵ=2.5\mathrm{\%}$, and 𝚼=0.71. We had 13 parameters in total. For this model setup, we ran EMCEE with 100 walkers, 4700 steps, and a burn-in fraction of 0.3.

4 One-population dynamical models

In this section, we describe the results we obtained when we only fitted the stellar kinematic data and did not try to separate the data into two populations. We ran several sets of models that tested whether our results were consistent with the known M_• mass (Sect. 4.1), if 𝚼 varied with radius (Sect. 4.2), and if a DM component was necessary (Sect. 4.3).

The results of the runs are summarised in Table 2. The percentage of background stars ϵ, the velocity anisotropy β_z, and the rotation parameter κ did not alter significantly when we changed the components that contribute to the gravitational potential. The results are consistent because the DM contribution is low, and our M_•. result agrees with the reference value. We use this knowledge for the two-population dynamical models in Sect. 5.

4.1 Fitting M_•

The mass of Sgr A^⋆, M_•, is well known to be (4.30±0.012) × 10⁶ M_⊙ (GRAVITY Collaboration 2022). We tested whether our data were consistent with this value by fitting the value of M_•.

In this fit with a radially constant 𝚼, our value for $M_{•}=\left({4.35}_{-0.23}^{+0.24}\right)\times {10}^6{M}_{\odot }$ is remarkably close to the stellar orbit measurements reported for instance by Boehle et al. (2016); GRAVITY Collaboration (2022). The mass of the supermassive black hole, Sgr A^⋆, has been attempted to be determined numerous times via dynamical models. The results from the stellar orbits often do not match, and the mass of Sgr A^⋆ is often under-estimated with dynamical models (see e.g. Magorrian 2019, and references therein). Our result is close to the reference value. This shows us that our data do not disagree with M_• =4.3 × 10⁶ M_⊙, and hence, they cause no biases in the other parameter fits. We therefore fixed M_• to the literature value in all future runs. We show the corner plot of the posterior distribution in Fig. A.1.

4.2 Radially varying 𝚼

We tested whether the conversion of the stellar number density (which traces the stellar light) into the stellar mass is constant as a function of radius. We parametrised 𝚼 with an Osipkov-Merritt profile (Eq. (1)), which means that there is an inner value 𝚼₀, an outer value 𝚼_∞, and a transition radius R_𝚼.

The inner value of 𝚼₀ is higher than the outer value 𝚼_∞, but because of the uncertainties, the difference is small (smaller than 1.5 times the sum of the uncertainties). The value of 𝚼_∞ is very well constrained, with much smaller uncertainties than for 𝚼₀. The radius at which the transition occurs also has a large upper uncertainty, and its 1σ upper limit is at ~13 pc. We show the 𝚼 profile in the bottom panel of Fig. 4. The radially constant 𝚼 that we fitted in the previous section lies in between our 𝚼₀, and 𝚼_∞, and it is consistent with 𝚼_∞ within the uncertainties. This suggests that the firm constraints on 𝚼 come from data in the outer region. The corner plot of the posterior distribution is shown in Fig. A.2.

The mass-to-light ratio in the mid-infrared is rather constant compared to other bands, and variations with age or metallicity are small, with only a factor of two for ages between 3 and 10 Gyr (Meidt et al. 2014). The values of 0.6–0.9 are consistent with ages of 5–10 Gyr and with a sub-solar to supersolar metallicity (Norris et al. 2014). Older stellar populations and populations with lower metallicity tend to have a higher 𝚼. This means that the higher 𝚼₀ agrees with either an older population or with a lower [Fe/H] compared to 𝚼_∞. A high concentration of dark remnants (stellar mass black holes or neutron stars) around Sgr A^⋆ may also increase 𝚼₀ compared to 𝚼_∞, however.

In summary, the profile 𝚼 only shows a slight radial dependence that only affects the innermost r ≲ 2 pc, and all other parameters are robust. We therefore assumed that 𝚼 is constant in our DM and two-population models.

Fig. 4 Radial profile of the velocity anisotropy βz (top), the rotation κ (middle), and the mass-to-number density conversion 𝚼 (bottom) as derived in Sect. 4.2 using a three-parameter function. The red lines denote the median of the posterior distribution, and the grey lines are 50 randomly drawn realisations from the posterior distribution. The vertical solid lines denote 1 Re = 5 pc, the dotted lines show the outer limit of the kinematic data, and the dot-dashed lines show the median of Rβ, Rκ, and R𝚼, respectively.

4.3 Dark matter scale density

We tested whether we needed to include DM in our models by fitting Jeans models with either a DM cusp (γ = 1) or a DM core (γ = 0). We only fitted the DM scale density ρ_S. We assumed that 𝚼 is constant with radius, and we fixed M_• to the literature value.

For DM cusp and DM core, the upper value of ρ_S is more tightly constrained than the lower value, especially for the cored DM profile. For both DM profiles, the enclosed mass assigned to DM within 1 R_e of the NSC is significantly below the stellar mass. In particular, the enclosed DM mass for a DM core profile at deprojected radius r=5 pc=1 R_e is below 0.5% (84th percentile upper limit) and lower than ~12% of the total mass at 33 pc (the outer extent of the kinematic data). For a DM cusp model, the enclosed mass of DM is higher, at ~1–4% of the total enclosed mass at 5 pc (16th and 84th percentile limits) and ~6–25% at 33 pc. The DM scale radius in these models was fixed to be much farther out than our data (r_DM=8 kpc). We tried a higher value (r_DM=17.2 kpc), and the DM mass fractions were in the same ranges. We show the corner plots of the posterior distributions in Figs. A.3–A.4.

Simulations of Milky Way analogues appear to favour an inner logarithmic density slope of 0.5–1.3 and exclude a cored DM profile (Hussein et al. 2025). Our assumed DM profiles of γ=1 may still not reflect the true distribution. Close to Sgr A^⋆, gravitational encounters, self-annihilation, or scattering may alter the density profile within r<1 pc (Merritt 2010). These processes are not likely to affect the outer regions, however. We tested two extreme cases for the DM profile, and both resulted in a small contribution to the total mass in the region of our data. We note that several model parameters (ϵ, β₀, β_∞, R_β, κ₀, κ_∞, R_κ) changed only slightly compared to the models without DM and within their uncertainties. The only exception was 𝚼, which changed by more than ~1σ for the cusp profile. Because most parameters were barely affected by the inclusion of a DM profile, we did not include a DM component for the two-population models because the values of ρ_S and γ are poorly constrained and the computing time increased significantly when we included DM.

Fig. 5 Map of the stellar positions (top) and position-velocity plots along the Galactic longitude l for the PM along l (second panel), along b (third panel), and along the line of sight (fourth panel). The colour-coding is from the median realisation of the one-population model in Sect. 4.2. Stars with a probability to be a member star higher than 0.5 are shown as an orange circle, and background stars are shown as a black pluses.

4.4 Model properties

Our one-population models give fully consistent results for the fraction of background stars ϵ, the velocity anisotropy β_z, and the rotation parameter κ. ϵ is consistently at 2.4%, which corresponds to ~450 stars with a probability $P_i^{bg}>0.5$. We show their spatial distribution and location in a position-velocity diagram in Fig. 5 using the 50th percentile realisation of Sect. 4.2 with a radially varying 𝚼 as an example. Other models gave similar results. The background stars were identified via high velocities in the PM or V_LOS component. Because only ~75% of the stars in our catalogue have PM measurements, the value of ϵ may be underestimated.

The models found a mild tangential anisotropy, with β_z ~ −0.1 to −0.2. They did not find a strong variation in the anisotropy with radius. β₀ and β_∞ differed by less than their uncertainties, and R_β was poorly constrained. We conclude that using a constant value for β_z is a fair assumption. An even more complicated β_z profile such as the profile used by Do et al. (2013b), with a fourth parameter to adjust the sharpness of the transition from β₀ to β_∞, is not required to describe our kinematic data. With a larger data set or more precise kinematic data, a more complicated functional form for β_z would be worth investigating. We show the β_z profiles drawn from the posterior distribution of Sect. 4.2 in the top panel of Fig. 4 as an example. The profiles were very similar in the other runs.

The rotation parameter κ is also consistent among the different models, and we show an example from Sect. 4.2 in the middle panel of Fig. 4. The radial change from κ₀ to κ_∞ indicates a lower ratio of the ordered to disordered motion in the centre than farther out. This change occurs within the central ≲1.2 pc, where Sgr A^⋆ dominates the gravitational potential. The value of κ₀ is even consistent with zero in the centre within the 84th percentile. This only holds for the inner ~1 pc. There is strong and significant rotational support at larger radii. The innermost 1 pc may have different dynamical properties, which may be caused by the influence of Sgr A^⋆. The velocity anisotropy or rotation parameters in the region of the inner few parsecs and the outer regions do not differ much, however, and we assumed that they are constant with radius in the two-population models.

The velocity dispersion maps in three dimensions (along the Galactic longitude l, along the latitude b, and along the line of sight), the V_LOS map, and V_LOS/σ_LOS maps of the same model in Fig. 6 show us that within the region of the NSC, there are distinct peaks in V_LOS(~ ±41 km s⁻¹) and V_LOS/σ_LOS that are followed by a drop at l ~ 10 pc. |V_LOS| and V_LOS/σ_LOS increase again at larger l. The peak value within the NSC is at V_LOS/σ_LOS = 0.68, the minimum is at 0.47, and it then continues to rise further to V_LOS/σ_LOS > 1.2 and V_LOS ~ ±58 km s⁻¹ at l ~30 pc.

We show the deprojected cumulative total mass profiles with a spherical radius in Fig. 7 with a constant 𝚼 and a free M_• (top left panel), or with varying 𝚼 and fixed M_• (top right panel), and for comparison, we show several models from the literature. Our mass profiles are lower than the literature profiles in the region ~5–20 pc, for some, even out to 60 pc. This is probably caused by the lack of data in this region in the literature. They had sparse data at best in this region (e.g. ~120 stars at r ≳ 6 pc in Feldmeier et al. 2014, ~200 stars at r ≳ 3.5 pc in Chatzopoulos et al. 2015, and fewer than 140 stars at r<20 pc in Sormani et al. 2022). A low number of stars that may be contaminated by foreground stars can cause an overestimation of the velocity dispersion, and hence, of the total dynamical mass. The cumulative mass of the model with a DM cusp is shown in the bottom left panel. We also list the enclosed total mass at four different spherical deprojected radii in Table 3 and the mass for projected radii in Table 4. Our resulting masses mostly agree with each other within the uncertainties.

Fig. 6 Map of the velocity dispersion in three directions (σl, σb, and σLOS), VLOS, and VLOS/σLOS for the median model with a radially varying 𝚼.

Fig. 7 Total enclosed mass as a function of spherical deprojected radius for the one-population fits with free M• (top left), radially varying 𝚼 (top right), a DM cusp (bottom left), and for the two-population fits (bottom right). The vertical solid lines denote 1 Re of the NSC, the dotted lines show the outer limit of the kinematic data, the dot-dashed line in the top right panel represents the median of R𝚼, and in the bottom right panel the median of Rh.

5 Two-population dynamical models

In this section, we describe the results we obtained when we fitted two stellar populations and an additional background component. As a chemical tracer, we used the overall metallicity [M/H].

For each population, we fitted a separate radially constant ${\beta }_z^k$ and κ^k. The two populations shared the same radially constant 𝚼. The black hole mass was fixed to the literature value, and any DM contribution was neglected. The population fraction h was left free to vary as a function of radius, parametrised with Eq. (1). Our results are summarised in Table 5.

The fraction of high-[M/H] stars is high, more than ~91%, which corresponds to ~3950 stars with a probability higher than 50% of belonging to this group. Only ~250 stars have a probability higher than 50% of belonging to the low-[M/H] population, and ~330 stars may belong to the background. About 70 stars cannot be attributed to either group and have ${\mathrm{P}}_i^k<50\mathrm{\%}$.

The two populations have a similar [M/H] width ${\sigma }_Z^k$. The [M/H] histogram is shown in Fig. 8. Both populations co-rotate, but the high-[M/H] population has a higher ratio of ordered-to-disordered motion (see also the fourth row of Fig. 10). Further, the velocity anisotropy is significantly different: the anisotropy is mildly tangential for the high-[M/H] population (slightly negative ${\beta }_z^1$), and it is significantly radial for the low-[M/H] population (positive ${\beta }_z^2$).

𝚼 is slightly higher than for the one-population models. The background fraction ϵ is significantly lower, only 1.71%. Some of the stars that are attributed to the background in the one-component models, including those on apparently counter-rotating orbits, are attributed to the radial anisotropic low-[M/H] population in the two-population models (compare Figs. 5 and 9). The high-[M/H] population is dominant at all radii, but it is slightly more dominant in the centre, with $h_0^1=0.98$, after which it decreases to $h_{\mathrm{\infty }}^1=0.91$. The transition radius of h lies within the NSC, with R_h ≈ 1.8 pc.

The high-[M/H] population resembles the one-population model in terms of the strong rotation signal and the shape of the V_LOS/σ_LOS map. The highest value of V_LOS/σ_LOS is again at the large |l|, but it is now lower and only reaches ~1.1 and not 1.25, as it does in the one-population fit. On the other hand, the low-[M/H] population has a higher σ_LOS (~80–90 km s⁻¹) than the high-[M/H] stars (50–60 km s⁻¹; see Fig. 10, third row), except in the innermost parsec region.

Fig. 8 [M/H] histogram. The grey histogram denotes all stars, the vertical dotted lines show the 50th percentile of Zk, and the dashed lines represent the resulting distributions of Zk and ${\sigma }_Z^k$. Pink denotes the high[M/H] population, cyan shows the low-[M/H] population, and orange shows the background population.

6 Discussion

Our discrete axisymmetric Jeans models constrained the central Galactic gravitational potential and the velocity structure of the stars very well. Our enclosed mass profile in the radial range of 5–30 pc is significantly lower than in previous studies, up to ~3.5 × 10⁷ M_⊙ below the Sormani et al. (2020) profile at ~30 pc (using the two-population mass profile as reference), although the largest discrepancy in terms of percent is at 15 pc, where the enclosed mass we found is lower by almost 75%. The agreement is better with a discrepancy of only ~1.5 × 10⁷ M_⊙ (at ~17 pc) with the Sormani et al. (2022) mass profile, and the mass is ≲60% lower (highest discrepancy at 12 pc) when we added a black hole mass of 4.3 × 10⁶ M_⊙ to their model. The triaxial orbit-based model of Feldmeier-Krause et al. (2017b) agrees better with our results, with a discrepancy of ~7 × 10⁶ M_⊙ within the range of the data at most, and ≲40% (highest discrepancy at 13 pc). Beyond the extent of our data, the different models diverge further in terms of cumulative mass, but their relative differences are ≲25% at 60 pc. As noted above, these studies were based on a smaller sample size than ours in the range of 5–30 pc. Their samples might also be contaminated with high-velocity stars, which are attributed to the background component that represents the Milky Way bar in our model. Both factors can bias the velocity dispersion to higher values, which then results in an overestimation of the enclosed mass.

Our result for M_• perfectly fits the value from stellar orbit fitting, which is remarkable considering the discrepant results of other stellar dynamical models. For example, Feldmeier et al. (2014) found $\left({1.7}_{-1.1}^{+1.4}\right)\times {10}^6{M}_{\odot }$, Chatzopoulos et al. (2015) (3.86±0.14) × 10⁶ M_⊙, Feldmeier-Krause et al. (2017b) $\left({3.0}_{-1.3}^{+1.0}\right)$ × 10⁶ M_⊙, and Magorrian (2019) (3.76±0.22)×10⁶ M_⊙. Our mass-to-light conversion factor 𝚼 also perfectly agrees with expectations for the 4.5 μm band (to which we scaled our number density distribution) from stellar population studies (Meidt et al. 2014; Norris et al. 2014; Kettlety et al. 2018).

Our fits with DM suggest that a core profile contributes a few percent to the total mass in the inner ~30 pc at best. Simulations excluded a core profile (Hussein et al. 2025) and favoured a steeper inner profile (γ=0.5 to 1.3), however, and our cusp profile is in this range. With a cusp profile, the DM contribution in our models is higher and reaches up to 6–25% at 33 pc. The DM volume density of simulations is not resolved in the region of our data, and Hussein et al. (2025) only showed it at r≳500 pc. We extrapolated our cusp DM profile and found that our DM volume density is higher by about a factor of ten than that of Hussein et al. (2025), and it is also higher than the DM density of the dynamical bar models by Portail et al. (2017). Our value of ρ_S has a firm upper limit, but a broader tail to lower values. This suggests that our DM contribution may be overestimated. Our other model parameters are robust and barely affected by the inclusion or exclusion of a DM component. This means that it was justified to neglect the DM contribution in the two-population models.

The velocity anisotropy β_z and κ vary little as a function of radius. Most of the stars are attributed to a high-[M/H] population, which has very mild tangential anisotropy. Moreover, the models reported by Sormani et al. (2022) indicated a tangential anisotropy in the inner 30 pc of the NSD.

The high-[M/H] population has strong rotational support, and it is close to an isotropic rotator because ${\beta }_z^1=-0.06$ is close to zero and κ¹ = −1.07 is close to −1. For an oblate isotropic rotator, we can use the minor-to-major axis ratio q to estimate V_LOS/σ_LOS. For the NSC, we have q=0.7, and thus V_LOS/σ_LOS = $\sqrt{(1-q)/q}=0.65$, for the NSD, q=0.35, and thus V_LOS/σ_LOS = 1.36. The Jeans model predicts that V_LOS/σ_LOS has a maximum within 1 R_e of the NSC at ~0.65 (see also Fig. 10, bottom left panel). This agrees with extragalactic NSCs in late-type galaxies, which have integrated values of (V_LOS/σ_LOS)_e (i.e. within 1 R_e) of 0.2–0.55 (Pinna et al. 2021). These values are lower than our maximum of V_LOS/σ_LOS, as Pinna et al. (2021) reported the flux-weighted integrated V_LOS/σ_LOS within 1 R_e, and these NSCs are not seen perfectly edge-on, as is the case for the Milky Way. Farther out, after a drop to V_LOS/σ_LOS ~0.43 at l ~9.5 pc, we found that V_LOS/σ_LOS predicted by the model increases even further to more than ~1. There is likely a maximum of V_LOS/σ_LOS, located in the NSD, but this is beyond the range of our data. In extragalactic NSDs, observations have detected maxima of V_LOS/σ_LOS, with values in the range of 1–3 and at radii of 200–1000 pc (Gadotti et al. 2020), which agrees with hydrodynamical simulations of a bar-built NSD (Cole et al. 2014). Our best-fit model is consistent with the observations of extragalactic NSDs, and the high V_LOS/σ_LOS of the high-[M/H] population supports the scenario that these stars formed in situ after gas inflow from the Galactic disc. To determine where the V_LOS/σ_LOS of the NSD has its maximum and what the maximum is, we need to model a larger region of the NSD. Kinematic data are available in the literature (Fritz et al. 2021), and we plan to do this in the future.

Our low-[M/H] population, on the other hand, is radially anisotropic and has weaker rotational support (V_LOS/σ_LOS ≲ 0.4 at l ≳2 pc), with an overall mean of V_LOS/σ_LOS = 0.2. The reason is twofold: The value of V_LOS near the Galactic plane (|b|≲0.5 pc) is ~±16 km s⁻¹, but ~ ±24–55 km s⁻¹ for the high-[M/H] population, whereas σ_LOS is ~90 km s⁻¹, but 50–60 km s⁻¹ for the high-[M/H] population. The velocity dispersions σ_l and σ_b are not significantly different for the two populations.

The fraction of background stars in the two-population models is lower than in the one-population model, and some of the stars attributed to the low-[M/H] population were previously assigned to the background component. The question therefore arises whether some stars of the low-[M/H] population are instead bar interlopers, with less extreme velocities and rather low [M/H]. While this may be the case for some stars, we found that the distributions of the low-[M/H] (i.e. stars with P² ≥ 0.5) and the background populations (i.e. stars with P^bg > 0.5) significantly different. We performed k-sample Anderson-Darling tests (Scholz & Stephens 1987) for V_l, V_b, V_LOS, [M/H], and the distance r to Sgr A^⋆ for the low-[M/H] and the background populations. In each of these tests, we obtained low p values, and the null hypothesis that the low-[M/H] and the background population samples come from the same distribution can be rejected with a confidence of 99.9%. We conclude that it is unlikely that the entire low-[M/H] population can be attributed to the background bar component.

One possible explanation for the low-[M/H] population is an external origin. Dynamical friction can move massive star clusters towards the Galactic centre, where the stars mix with the in situ stellar population (e.g. Tremaine et al. 1975; Capuzzo-Dolcetta 1993; Antonini et al. 2012; Arca Sedda et al. 2020; van Donkelaar et al. 2024). Tsatsi et al. (2017) showed with N-body simulations that even an NSC built entirely from 12 consecutively infalling star clusters coming from random orbital directions can have rotational support (V_LOS/σ_LOS = 0.1–0.7). Milky Way globular clusters on eccentric orbits may also approach the Galactic centre and deposit some of their stars. Ishchenko et al. (2023) found that some of the known globular clusters of the Milky Way have a minimum distance to the Galactic centre of several dozen parsec. The ten clusters with the highest probability of such an encounter in Ishchenko et al. (2023) include several with [Fe/H] in the range of −1.5 to −1.0 dex (Carretta et al. 2009), which is close to our low-[M/H] sample (−1.5 dex<[M/H]<−0.5 dex). Their contribution to our data set may be low. Former already-dissolved clusters may also have deposited stars in the past. An external origin might explain why the low-[M/H] population has different values for κ and β_z. These star clusters may have been on radial orbits that brought them close to the Galactic centre, which explains the radial anisotropy of the low-[M/H] stars we measured. As the two-body relaxation time in the NSC is shorter than ~3 Gyr (Arca Sedda et al. 2020), this suggests that such infalls or passages occurred in the last 3 Gyr. To constrain the number of these infall events, more precise [M/H] measurements and chemical abundances are needed.

Our axisymmetric Jeans modelling approach has some limitations because it (1) assumes an axisymmetric shape and (2) makes assumptions on the alignment of the velocity ellipsoid. Although axisymmetric models fit the data very well also for more extended NSD data (e.g., Sanders et al. 2024), it is unclear whether the NSD is truly axisymmetric or if the Milky Way has an inner bar (Alard 2001; Gonzalez et al. 2011), as was found in several galaxies (Gadotti et al. 2019; Erwin 2024) at scales of 0.1–1 kpc. Sormani et al. (2022) suggested that the kinematic signature of a nuclear or inner bar may be seen in an asymmetry of the V_l distribution. A larger statistical sample and a better understanding of the expected signature of a nuclear bar from theoretical studies are needed, however. The second limitation (assumptions on the alignment of the velocity ellipsoid) can be overcome by using orbit-based triaxial models (Zhu et al. 2020), but these models are not yet available for discrete data.

Fig. 9 Map of the stellar positions (top), position-velocity plots along the Galactic longitude l for the PM along l (second panel), along b (third panel), along the line of sight (fourth panel), and along the stellar metallicity [M/H]. The colour-coding is from the 50th percentile realisation of the two-population model and is shown by the colour bar on the top. Pink denotes a high probability to be a high [M/H] star, cyan shows a low [M/H] star, and the black plus shows a background star. The orange circles denote stars with a population membership probability ${\mathrm{P}}_i^k<0.5$.

Fig. 10 Spatial distribution of stars, colour-coded from top to bottom by the velocity dispersion along l, along b, along the line of sight, VLOS, VLOS/σLOS, for the 50 th percentile two-population model, with high-[M/H] stars (left), and with low-[M/H] stars (right).

7 Conclusion

We have constructed the first discrete axisymmetric chemodynamical Jeans models of the Galactic centre. The chemodynamic data consisted of V_LOS and [M/H] for 4600 unique stars, 3567 of which also had proper motions. Foreground and background stars were excluded based on H − K_S colour cuts. We accounted for the remaining contribution of bar stars by including a background component in the models. Our data extend out to |l|=33 pc from Sgr A^⋆ along the Galactic longitude, into a region in which other studies only had several hundred stars.

We constructed purely dynamical models for which only the stellar density distribution and kinematic data were used. These models showed that the velocity anisotropy β_z ≈ −0.1, the rotation parameter κ ≈ −1, and the mass-to-light ratio 𝚼 ≈ 0.7 depend little on the radius and can be assumed to be constant. When we fitted the mass of Sgr A^⋆, we obtained $({4.35}_{-0.23}^{+0.24}\times {10}^6)M_{\odot }$, which perfectly agrees with literature values of stellar orbits. We tested a cored (γ=0) and a cusped (γ=1) DM profile, but the mass contribution in the NSC was only a few percent in either case, and we decided to neglect DM. We listed the enclosed total mass of various models for several deprojected and projected radii.

In chemo-dynamical models, we assumed two populations k of stars, each with a unique but radially constant value for ${\beta }_z^k$ and κ^k, and Gaussian [M/H] distributions with a unique centroid $Z_0^k$ and width ${\sigma }_Z^k$. The two populations shared the same 𝚼; their relative contribution was a function of radius. We obtained a dominant (~90%) [M/H]-rich population with a mild tangential anisotropy (${\beta }_z^1=-0.06$) and strong rotation (κ¹ = −1.06), which is consistent with being formed in situ after gas inflow from the outer Galactic disc transported in via the bar. The sub-dominant [M/H]-poor population has a significant radial anisotropy (${\beta }_z^2=+0.64$) and is less well supported by rotation (κ² = −0.6), which indicates a different origin.

For a better understanding of the Galactic centre mass distribution and orbital structure, it would be interesting to relax some assumptions of the axisymmetric Jeans models and use triaxial orbit-based models, for example, which make no assumptions on the intrinsic velocity distribution. Concerning the data, a larger sample with more precise velocity and [M/H] measurements and possibly even chemical abundances would allow us to obtain a more thorough analysis of the subpopulations and even to include more than two populations. This would help us to better understand the history of merger and accretion events, and thus, the formation of the nuclear region of the Milky Way.

Acknowledgements

We thank the anonymous referees for their thorough review and constructive comments that improved the manuscript. A.F.K. acknowledges funding from the Austrian Science Fund (FWF) [grant DOI 10.55776/ESP542]. This research made use of NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), matplotlib, a Python library for publication quality graphics (Hunter 2007), Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018, 2022). Calculations were performed using supercomputer resources provided by the Vienna Scientific Cluster (VSC).

Appendix A MCMC post burn-in distributions

In Figs. A.1-A.5 we show the post burn-in probability distributions of our various EMCEE fits.

Fig. A.1 MCMC post burn-in distributions for the one-population models with M• as free parameter (Sect. 4.1). The scatter plots show the projected two-dimensional distributions, the red lines the respective 50th percentile, and the dashed grey lines the 16th and 84th percentiles. Their values are also written on top of each column, displaying the projected one-dimensional distributions. From top to bottom and left to right, the panels show the inner anisotropy β0, outer anisotropy β∞, the anisotropy transition radius log(Rβ), inner rotation parameter κ0, outer rotation parameter κ∞, the rotation transition radius log(Rκ), the mass-to-light ratio 𝚼, fraction of background stars ϵ in percent, and the black hole mass M• in 106 M⊙.

Fig. A.2 Same as Fig. A.1, but for the one-populations models with a varying three-parameter mass-to-light ratio (Sect. 4.2). From top to bottom and left to right, the panels show the inner anisotropy β0, outer anisotropy β∞, the anisotropy transition radius log(Rβ), inner rotation parameter κ0, outer rotation parameter κ∞, the rotation transition radius log(Rκ), the inner mass-to-light ratio 𝚼0, outer mass-to-light ratio 𝚼∞, mass-to-light ratio transition radius log(R𝚼), and the fraction of background stars ϵ in percent. The black hole mass is fixed to 4.3 × 106 M⊙, dark matter is not considered in the gravitational potential.

Fig. A.3 Same as Figs. A.1–A.2, but for the one-population models with a cuspy dark matter profile (Sect. 4.3). From top to bottom and left to right, the panels show the inner anisotropy β0, outer anisotropy β∞, the anisotropy transition radius log(Rβ), inner rotation parameter κ0, outer rotation parameter κ∞, the rotation transition radius log(Rκ), the mass-to-light ratio 𝚼, fraction of background stars ϵ in percent, and the dark matter scale density ρS. The black hole mass is fixed to 4.3×106 M⊙.

Fig. A.4 Same as Fig. A.3, but for the one-population models with a cored dark matter profile (Sect. 4.3).

Fig. A.5 Same as Figs. A.1-A.4, but for the two-population models (Sect. 5). From top to bottom and left to right, the panels show the mean metallicity $Z_0^1$, metallicity dispersion ${\sigma }_Z^1$, anisotropy ${\beta }_z^1$, rotation parameter κ1 for the high [M/H] population, $Z_0^2,{\sigma }_Z^2,{\beta }_z^2,{\kappa }^2$ for the low [M/H] population, inner population fraction h0, outer population fraction h∞, the transition radius log(Rh), fraction of background stars ϵ in percent, and the mass-to-light ratio 𝚼. The black hole mass is fixed to 4.3×106 M⊙, dark matter is not considered in the gravitational potential.

References

All Tables

All Figures

	Fig. 1 Spatial distribution of all stars we used for the model. The black crosses denote stars with VLOS and [M/H], red diamonds show stars with additional proper motions from Fritz et al. (2016), blue circles show stars from Libralato et al. (2020), and orange crosses show stars from Smith et al. (2025).
In the text

	Fig. 2 Colour-magnitude diagram H − KS vs. KS of the stars with VLOS and [M/H]. The vertical green lines denote the colour cuts we used to remove likely foreground (H − KS ≤ 1.3 mag) and background stars (H − KS ≥3.5 mag). The colour represents the density of stars.
In the text

Fig. 3 Upper panel: surface density profile derived from the stellar density maps of Gallego-Cano et al. (2020, blue diamonds). The uncertainties are measured using the corresponding uncertainty map and are only shown for a subset of data points to improve visibility. We matched the profile to the centre value of Feldmeier-Krause et al. (2017b) to convert it into units of L⊙ pc−2. The solid black line denotes the profile along the major axis, the dashed red line shows the profile along the minor axis, and the dot-dashed orange line shows the profile along a position angle of 45°. Middle panel: residuals of the star count data after subtracting the MGE fit. Lower panel: projected axis ratio qproj.

In the text

Fig. 4 Radial profile of the velocity anisotropy βz (top), the rotation κ (middle), and the mass-to-number density conversion 𝚼 (bottom) as derived in Sect. 4.2 using a three-parameter function. The red lines denote the median of the posterior distribution, and the grey lines are 50 randomly drawn realisations from the posterior distribution. The vertical solid lines denote 1 Re = 5 pc, the dotted lines show the outer limit of the kinematic data, and the dot-dashed lines show the median of Rβ, Rκ, and R𝚼, respectively.

In the text

Fig. 5 Map of the stellar positions (top) and position-velocity plots along the Galactic longitude l for the PM along l (second panel), along b (third panel), and along the line of sight (fourth panel). The colour-coding is from the median realisation of the one-population model in Sect. 4.2. Stars with a probability to be a member star higher than 0.5 are shown as an orange circle, and background stars are shown as a black pluses.

In the text

	Fig. 6 Map of the velocity dispersion in three directions (σl, σb, and σLOS), VLOS, and VLOS/σLOS for the median model with a radially varying 𝚼.
In the text

Fig. 7 Total enclosed mass as a function of spherical deprojected radius for the one-population fits with free M• (top left), radially varying 𝚼 (top right), a DM cusp (bottom left), and for the two-population fits (bottom right). The vertical solid lines denote 1 Re of the NSC, the dotted lines show the outer limit of the kinematic data, the dot-dashed line in the top right panel represents the median of R𝚼, and in the bottom right panel the median of Rh.

In the text

	Fig. 8 [M/H] histogram. The grey histogram denotes all stars, the vertical dotted lines show the 50th percentile of Zk, and the dashed lines represent the resulting distributions of Zk and ${\sigma }_Z^k$. Pink denotes the high[M/H] population, cyan shows the low-[M/H] population, and orange shows the background population.
In the text

Fig. 9 Map of the stellar positions (top), position-velocity plots along the Galactic longitude l for the PM along l (second panel), along b (third panel), along the line of sight (fourth panel), and along the stellar metallicity [M/H]. The colour-coding is from the 50th percentile realisation of the two-population model and is shown by the colour bar on the top. Pink denotes a high probability to be a high [M/H] star, cyan shows a low [M/H] star, and the black plus shows a background star. The orange circles denote stars with a population membership probability ${\mathrm{P}}_i^k<0.5$.

In the text

	Fig. 10 Spatial distribution of stars, colour-coded from top to bottom by the velocity dispersion along l, along b, along the line of sight, VLOS, VLOS/σLOS, for the 50 th percentile two-population model, with high-[M/H] stars (left), and with low-[M/H] stars (right).
In the text

Fig. A.1 MCMC post burn-in distributions for the one-population models with M• as free parameter (Sect. 4.1). The scatter plots show the projected two-dimensional distributions, the red lines the respective 50th percentile, and the dashed grey lines the 16th and 84th percentiles. Their values are also written on top of each column, displaying the projected one-dimensional distributions. From top to bottom and left to right, the panels show the inner anisotropy β0, outer anisotropy β∞, the anisotropy transition radius log(Rβ), inner rotation parameter κ0, outer rotation parameter κ∞, the rotation transition radius log(Rκ), the mass-to-light ratio 𝚼, fraction of background stars ϵ in percent, and the black hole mass M• in 106 M⊙.

In the text

Fig. A.2 Same as Fig. A.1, but for the one-populations models with a varying three-parameter mass-to-light ratio (Sect. 4.2). From top to bottom and left to right, the panels show the inner anisotropy β0, outer anisotropy β∞, the anisotropy transition radius log(Rβ), inner rotation parameter κ0, outer rotation parameter κ∞, the rotation transition radius log(Rκ), the inner mass-to-light ratio 𝚼0, outer mass-to-light ratio 𝚼∞, mass-to-light ratio transition radius log(R𝚼), and the fraction of background stars ϵ in percent. The black hole mass is fixed to 4.3 × 106 M⊙, dark matter is not considered in the gravitational potential.

In the text

Fig. A.3 Same as Figs. A.1–A.2, but for the one-population models with a cuspy dark matter profile (Sect. 4.3). From top to bottom and left to right, the panels show the inner anisotropy β0, outer anisotropy β∞, the anisotropy transition radius log(Rβ), inner rotation parameter κ0, outer rotation parameter κ∞, the rotation transition radius log(Rκ), the mass-to-light ratio 𝚼, fraction of background stars ϵ in percent, and the dark matter scale density ρS. The black hole mass is fixed to 4.3×106 M⊙.

In the text

	Fig. A.4 Same as Fig. A.3, but for the one-population models with a cored dark matter profile (Sect. 4.3).
In the text

Fig. A.5 Same as Figs. A.1-A.4, but for the two-population models (Sect. 5). From top to bottom and left to right, the panels show the mean metallicity $Z_0^1$, metallicity dispersion ${\sigma }_Z^1$, anisotropy ${\beta }_z^1$, rotation parameter κ1 for the high [M/H] population, $Z_0^2,{\sigma }_Z^2,{\beta }_z^2,{\kappa }^2$ for the low [M/H] population, inner population fraction h0, outer population fraction h∞, the transition radius log(Rh), fraction of background stars ϵ in percent, and the mass-to-light ratio 𝚼. The black hole mass is fixed to 4.3×106 M⊙, dark matter is not considered in the gravitational potential.

In the text