EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 711 (July 2026) A&A, 711 (2026) A98 Full HTML
Open Access
Issue
A&A
Volume 711, July 2026
Article Number A98
Number of page(s) 16
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202659834
Published online 03 July 2026

© The Authors 2026

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. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1 Introduction

Galaxy centres play an important role in the formation and evolution of galaxies. Inflow and outflow of gas to and from the galaxy centre is common, and found in about 50% of nearby low-ionisation nuclear emission-line regions (LINERs, Hermosa Muñoz et al. 2022). Although galaxy centres only extend over a small spatial scale of the galaxy (up to a few hundred par-secs), tight scaling relations between properties of the galaxy centre and the more extended galaxy have been found. These scaling relations imply that the large-scale galaxy assembly and the build-up of the galaxy centre are connected. For example, there is a correlation between the mass of the supermassive black hole (M•) and the host galaxy stellar velocity dispersion (e.g. Ferrarese & Merritt 2000; Gebhardt et al. 2000; Kormendy & Ho 2013), and a correlation between the mass of the nuclear star cluster (NSC) and the galaxy stellar mass (e.g. Ferrarese et al. 2006; Scott & Graham 2013; Georgiev et al. 2016). The slope of these relations, their scatter, and behaviour at different mass regimes and environments are still areas of active research.

The most reliable way to measure an NSC mass is via stellar dynamical modelling (Neumayer et al. 2020), using stellar kinematic data. In most NSCs, we can so far only use the spec-troscopic integrated stellar light for dynamical modelling (Barth et al. 2009; Nguyen et al. 2018). The Milky Way (MW)'s centre is one of the few galaxy centres where we can resolve single stars and use those to constrain the mass distribution. This makes the Galactic centre (GC) an excellent benchmark object where we can test the methods and techniques used to infer the mass distributions in other galaxies.

The GC contains several massive structures, each of which dominates the total gravitational potential at a certain scale: (1) the mildly flattened NSC with an effective radius, Re, of about 5 pc (Fritz et al. 2016; Gallego-Cano et al. 2020) and a total mass of a few 107 M⊙ (e.g. Feldmeier et al. 2014; Chatzopoulos et al. 2015; Feldmeier-Krause et al. 2017b); (2) the more flattened surrounding nuclear stellar disc (NSD) with a scale radius of about 100 pc, a scale height of ~30 pc (Philipp et al. 1999; Launhardt et al. 2002; Sormani et al. 2022), and a mass of ~109 M⊙ (Sormani et al. 2022); and (3) the central supermassive black hole Sgr A⊙. The mass of Sgr A⊙ is known with high precision, M• =(4.30± 0.012)·•106 M⊙, (GRAVITY Collaboration 2022), thanks to the long-time monitoring of stellar orbits (e.g. Boehle et al. 2016; Gillessen et al. 2017; Do et al. 2019). The mass distributions of the MW's NSC and NSD have been studied with various dynamical modelling approaches, some of which assume spherical symmetry (e.g. Magorrian 2019), axisymmetry (e.g. Chatzopoulos et al. 2015; Sormani et al. 2022; Feldmeier-Krause et al. 2025b; Vasiliev et al. 2026), or triaxiality (Feldmeier-Krause et al. 2017b). More recent works used discrete stellar velocities rather than binning the data, and combined the line-of-sight velocities, VLOS, obtained from spectroscopy with proper motions. This is usually not possible for extragalactic systems, where the kinematic information is often limited to line-of-sight data from the integrated light, and thus spatially binned rather than discrete.

Nonetheless, besides a map of the integrated VLOS, more information can be extracted from integrated light data, such as maps of the velocity dispersion, σLOS, or higher-order moments of the line-of-sight velocity distribution (LOSVD, e.g. Cappellari & Emsellem 2004). Even non-parametric LOSVDs have been extracted from optical spectroscopic data (Mehrgan et al. 2019; Falcón-Barroso & Martig 2021). In combination with the stellar light distribution, such kinematic data can be used to infer the total mass distribution. Stellar dynamical modelling approaches make assumptions about the shape of the stellar systems (i.e. spherical, axisymmetric, or triaxial) and assume that the gravitational potential is static. Models based on the Jeans equations (Jeans 1922) require additional assumptions about the velocity structure. One of the commonly used approaches that does not make any such assumptions is based on Schwarzschild (1979); for example, in van den Bosch et al. (2008), Vasiliev (2013) and Neureiter et al. (2021). The codes have been improved and further developed in recent years (e.g. Poci et al. 2019; Jethwa et al. 2020; Thater et al. 2022a), and can even include barred structures (Vasiliev & Valluri 2020; Tahmasebzadeh et al. 2024; Jin et al. 2025). This so-called orbit-based modelling has been validated on simulations several times (e.g. Jin et al. 2019; Zhu et al. 2020; Neureiter et al. 2023; Jin et al. 2025). Unbarred triaxial Schwarzschild models have been used to model the nuclear regions of galaxies and infer the central mass distribution and stellar orbital structure (Lyubenova et al. 2013; Feldmeier-Krause et al. 2017b; Fahrion et al. 2019; den Brok et al. 2021; Thater et al. 2023, 2026; Lamprecht et al. 2026).

The GC is an interesting object with which to test stellar dynamical models. However, several studies based on Jeans models underestimated M• (e.g. Schödel et al. 2009; Feldmeier et al. 2014; Fritz et al. 2016). Also, triaxial (Feldmeier-Krause et al. 2017b) and spherical (Magorrian 2019) orbit-based models could barely recover the mass of Sgr A⋆. Their kinematic data covered only the inner region of the NSC. More recently, discrete axisymmetric Jeans models of Feldmeier-Krause et al. (2025b), using more extended kinematic data and a precise stellar density distribution, recovered M• very accurately as (4.35−0.23+0.24)⋅106M⊙Mathematical equation: $\left( {4.35_{ - 0.23}^{ + 0.24}} \right) \cdot {10^6}{M_ \odot }$. In this work, we test if we can also obtain the correct value for M• using triaxial orbit-based models and integrated light stellar kinematics, utilising the same density distribution and kinematic data that extend over a similar area as in Feldmeier-Krause et al. (2025b). Our models also deliver the mass distribution and orbital structure of the NSC and inner NSD. We decompose the stellar orbits into dynamically cold (high angular momentum) and hot (low angular momentum) components. These shed light on the orbit distribution of the GC stellar structures.

This paper is organised as follows. We describe the spec-troscopic data and kinematics in Sect. 2, and the orbit-based modelling approach in Sect. 3. We present our results in Sect. 4, discuss them in Sect. 5 and conclude in Sect. 6.

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Spatial coverage of our F2 spectroscopic data. The data extend ~66 pc along the Galactic longitude l, centred on Sgr A⋆ (marked as a red plus symbol), and ~1 pc to the Galactic north and south, except for the centre region, which extends further to the Galactic north (~2pc). The image was constructed from the spectroscopic scans. We show the Galactic co-ordinate grid as dashed lines.

2 Spectroscopic data

2.1 Observations and data cube construction

The observations were taken on five nights (June 24, 25, 26, 27, and 29, 2015) with Flamingos-2 (F2, Eikenberry et al. 2004) at the Gemini-South telescope. We observed in the near-infrared K band, which is less affected by interstellar extinction (AKS≈2.5 magMathematical equation: ${A_{{K_S}}} \approx 2.5{\rm{ mag}}$ mag, e.g. Nogueras-Lara et al. 2018) than optical wavelengths (AV≈30 mag, Fritz et al. 2011). The entire coverage of the spectroscopic data, meaning all fields combined, extends over ~66 pc along the Galactic longitude, l, with Sgr A⋆ in the centre, and ~2 pc along the Galactic latitude, b (~3 pc in the inner |l| ≲8pc); see Fig. 1.

The observations and basic data reduction are described in detail in Feldmeier-Krause et al. (2025a). Here we just give a summary: We observed five regions of the GC, which we labelled outer west, inner west, central, inner east, and outer east. These names are for Galactic co-ordinates and relative to Sgr A⋆. In each region, we observed 50 exposures (87 exposures in the central region) with a 6′ long slit mask, and the telescope drifted by 1″ per 300 s exposure perpendicular to the slit mask. After 6–22 subsequent exposures, which we call a sequence, during which the telescope continuously drifted towards the Galactic South, we interrupted the series for sky and telluric calibration exposures. After these, we continued to scan the region to obtain 50 exposures in total. Due to these telescope offsets, in some cases, we have overlapping exposures or small gaps (see the horizontal white line in the inner west of Fig. 1) between the observation sequences. We observed 20 sequences in total with the 6′ long slit. The slit mask consists of six approximately 1′ long slits or slitlets, aligned in a single row, with five small regions that were not cut to stabilise the mask. These regions cause small gaps in our data every 1′ along l (see vertical white lines in Fig. 1).

The data reduction includes dark subtraction, persistence subtraction, flat field division, cosmic ray removal, rectification, sky subtraction, and telluric correction. For each of the 20 sequences, divided into six slitlet regions, we constructed a stitched image. These 120 images were each cross-correlated with a Vista Variables in the Via Lâctea KS image (Saito et al. 2012) to obtain their astrometric calibration. Figure 1 is a mosaic of these 120 images after astrometric calibration.

In Feldmeier-Krause et al. (2025a), we extracted and analysed the spectra of the brightest stars, but there is also valuable information in the light of the fainter stars in the data. For this reason, we here use these data to construct data cubes of the unresolved light. In the data cubes, we masked all foreground and bright stars, as a few bright stars can dominate the light and outshine the unresolved stars (Feldmeier et al. 2014; Davidge 2020). We created masks as follows. We used the JHKS band photometry of the GALACTICNUCLEUS (GNS) survey catalogue (Nogueras-Lara et al. 2019) to identify the approximate positions of stars on the slit. Starting with the brightest star per observing sequence, we fitted a Gaussian function to determine its exact location along the slit and then masked out a region of 7σ. For stars fainter than KS = 14 mag, we considered only the primary exposure where the slit covers the star and no adjacent exposures. But brighter stars (KS < 14 mag) can also contribute significantly to the exposure taken immediately before and after the one where the star is centred due to the seeing, and we also masked the flux of such stars in those exposures. We created separate masks for foreground stars, which we identified using their H – KS colour. Foreground stars are less reddened than GC stars, and thus bluer. Our foreground star maps include all the stars with H – KS ≤ 1.3 mag (Feldmeier-Krause 2022), and for stars with unknown colour, where either the H or KS band photometry are missing. This can be the case for bright stars that are saturated in the GNS survey.

As done for the stitched images, we combined subsequently taken exposures to stitched data cubes. We used the previously created masks to remove the light from foreground stars, stars with unknown colour (as they may also be foreground stars), and bright stars (KS,0=KS−AKS<9.0 magMathematical equation: ${K_{S,0}} = {K_S} - {A_{{K_S}}} < 9.0{\rm{ mag}}$). The brightness cut was chosen to be close to KS,cut = 11.5 mag, to match the cut chosen by Feldmeier et al. (2014). However, instead of using a fixed value of KS,cut for the observed KS band photometry, we fixed the value of the extinction corrected KS band photometry, KS,0, to account for spatial variations in the extinction, AKSMathematical equation: ${A_{{K_S}}}$. We used the extinction map of Feldmeier-Krause et al. (2025a) and the respective mean extinction AKs in the 120 image regions. For confirmed and possible foreground stars, we used an even stricter magnitude cut and masked them down to KS < 15 mag.

We applied a wavelength calibration correction determined on the sky lines (see Feldmeier-Krause et al. 2025a) for each exposure. We combined the 120 stitched data cubes such that the 50 exposures (86 in the centre) along latitude b are combined to a single data cube, but keeping the separation into six 1′ wide slitlets, using the code IFSR_MOSAIC.pro1. This results in 30 (5 regions × 6 slitlets) data cubes. We then rebinned the data cubes with IFSR_REBIN.pro to ~1″ pixel−1.

The median number of masked foreground stars per exposure and slitlet (covering ~1′x 1″) is 0.8. The median number of masked bright stars is one. The final data cube contains only stars fainter than KS ~ 11.5 mag. According to the GNS catalogue, the median number of remaining stars with 11.5 < KS < 15 mag is 16 per exposure and slitlet. The minimum number is ~8 stars in this magnitude range in the outer regions, and in the denser central region, we have a maximum of ~70 stars.

2.2 Spatial binning

We binned our 30 data cubes spatially before we measured the stellar kinematics. To ensure that we have a sufficiently large number of stars per bin, we used the S/N of the data cubes and applied the Voronoi binning code of Cappellari & Copin (2003). We slightly modified the PYTHON code such that it does not start at the pixel with the highest flux but rather in the centre of the respective data cube. We required a target S/N of 75 per bin, leading to up to 45 bins in the central field data cubes but only 1 bin in some outer data cubes, where the stellar density is lower. In total, we have 197 Voronoi bins distributed over the entire region. For each Voronoi bin, we have a spectrum that contains the integrated flux of the unresolved stars in the region.

Using a higher value for the target S/N per bin leads to coarser spatial binning in the centre of our field, averaging out spatial information. A lower target S/N produces a finer spatial binning, but at the cost of noisier kinematic maps. We tried different target S/N values (ranging from 45–95) and found that a target S/N of 75 leads to robust kinematic results without compromising valuable spatial information.

2.3 Stellar kinematics

We measured the stellar kinematics using the PYTHON code PPXF (Cappellari 2023) on the Voronoi binned spectra. This code requires template spectra, and we used the high-resolution spectral library of late-type stars by Wallace & Hinkle (1996), convolved to the respective spectral resolution of the data, in the wavelength region of 2.21–2.365 μm. This wavelength region does not include the Na I doublet, as we found systematic residuals in this region in some of the spectra, but includes the Ca I triplet, and four CO transitions (12CO (v=2−0) 2.2935 μm, 12CO (3−1) 2.3227 μm, 12CO (4−2) 2.3525 μm and 13CO (2−0) 2.3448 μm. We fitted four Gauss-Hermite moments of the LOSVD. The first and second moments correspond to the velocity VLOS and velocity dispersion σLOS in the limit where the Gauss-Hermite coefficients h3 and h4 are equal to zero. We used additive Legendre polynomials with degree=4 to correct the template continuum shape during the fit, as recommended for kinematic fits.

The spectral resolution was measured in Feldmeier-Krause et al. (2025a) by fitting Gaussians to the sky emission lines. They found that the spectral resolution varies both as a function of wavelength and spatially for individual slitlets, with a maximum spectral resolution R=3400. For each slitlet, they derived a spectral resolution function, which is a second-degree polynomial function of the wavelength, and we used these functions for the spectral analysis.

It is necessary to include H2 gas emission lines to fit the inner regions of our data. The gas emission originates from the circum-nuclear disc (CND) of molecular gas in the inner r ≲ 3 pc around Sgr A⋆. We included the transitions at 2.2235 μm (H2 1–0 S(0)) and 2.2477 μm (H2 2–1 S(1)), constrained to the same gas velocity and velocity dispersion. There are other H2 transitions in the fitted wavelength region (e.g. 2.3556 μm H2 2−1 S(0)), but they are very weak, and we did not include them in the fit.

To reach our final kinematic results, we started with an initial fit, where we masked the wavelength region 2.3145–2.3175 μm, as this region can have residuals from the sky subtraction procedure. In a second fit, we masked in addition all pixels where the absolute value of the residual spectrum (= data – best-fit) from the first fit exceeds 0.15. In a third fit, we set in addition the PPXF keyword CLEAN, which does iterative sigma-clipping to remove unmasked bad pixels (Cappellari et al. 2002). We usually used the kinematic results from the third fit, unless more than 20 percent of all pixels were masked in this fit, then we used the result obtained by the second fit. The results are usually close, differences are only a few kilometres per second and within the measurement uncertainties. We show the example of a typical spectrum and the best-fit stellar and gas model in Fig. 2. To obtain the kinematic errors, we performed Monte Carlo simulations. This was done by adding noise to the spectrum in 500 iterations and fitting the kinematics (with CLEAN=False). We use the standard deviation of these fits as our kinematic uncertainties.

We tested fits where we include the Na I doublet, or to limit the fitting range to the first two CO lines (12CO (v=2–0) 2.2935 μm, 12CO (3–1) 2.3227 μm). While these generally achieve similar results, our chosen wavelength range is better at dealing with spectra that are affected by sky residuals in the outer spatial regions of our data, leading to more robust results.

We also tried different template spectra for the fit. The X-SHOOTER spectral library (XSL) single stellar population (SSP) models (Verro et al. 2022) have a sufficiently high spectral resolution and long wavelength range, but these templates had larger residuals at the Ca I feature in comparison to the Wallace & Hinkle (1996) templates. This may be explained in the context of spectral index measurements of resolved stars. Stars in the GC have stronger Ca I and Na I EW values compared to Galactic disc stars (e.g. Blum et al. 1996; Feldmeier-Krause et al. 2017a). This trend also holds for integrated light spectra (Davidge 2020). The Wallace & Hinkle (1996) stellar templates may have more flexibility in finding the best template than the XSL SSP models. The widely used SSP templates of Vazdekis et al. (2016, EMILES) and Conroy et al. (2018) have a lower spectral resolution than our data; hence, we chose not to use them.

We applied point-symmetrisation on our kinematic maps, using the SYMMETRIZE_VELFIELD.py code in the PLOT-bin python package2. Symmetrising data reduces noise and removes systematic effects such as the systemic velocity in an extragalactic system (van den Bosch & de Zeeuw 2010), without causing significant biases on the modelling results (Walsh et al. 2012; Thater et al. 2022b). We do not alter the uncertainties of our maps. We show the point-symmetrised kinematic maps (left) and their uncertainties (right) in Fig. 3. The uncertainties are lower in the central region of the maps. Note the rotation seen in the Vlos map and the anti-correlation between VLos and h3.

Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Example of a stellar kinematic fit with PPXF. The black line shows the data, the red line the best-fit stellar model, and the orange and pink lines the best-fit gas emission. Green symbols denote residuals; the shaded grey regions were masked in the fit due to sky residuals or bad pixels.

3 Orbit-based modelling

We used orbit-based modelling based on van den Bosch et al. (2008) with the tool DYNAMITE (DYnamics, Age and Metallicity Indicators Tracing Evolution, Jethwa et al. 2020; Thater et al. 2022a). The code finds the best combination of orbits in a given potential, and the best set of hyperparameters to describe the gravitational potential. The models are constrained by the four stellar kinematic maps shown in Fig. 3. In this section, we describe how we model the gravitational potential and compute the orbit library.

3.1 Gravitational potential

We assumed a gravitational potential that consists of a supermassive black hole (SMBH) and the stellar distribution. In a subset of models, we also included a spherical dark matter (DM) component. In our models, we assumed a galactocentric distance of 8.3 kpc (GRAVITY Collaboration 2022).

3.1.1 Supermassive black hole

The SMBH was modelled as a Plummer sphere, with a scale radius of 0.″01. The scale radius is much smaller than the spatial resolution and pixel scale of our data; hence, the potential is a good approximation of a point mass. The SMBH mass, M•, is a free hyperparameter and was fit on a logarithmic scale. We restricted M• to 106.201–106.878 Μ⊙ = (1.6–7.6) • 106 M⊙. The well-constrained value of M• = (4.3 ± 0.012) • 106 M⊙ (GRAVITY Collaboration 2022) was used as starting value.

3.1.2 Stellar mass

We used the stellar number density map of Gallego-Cano et al. (2020), which traces red giant stars in the GC (~84.4pc × 21 pc). The stars are in the extinction corrected magnitude range of 9.0 ≤ KS,0 ≤ 14.0 mag. Our kinematic maps also trace red giant stars KS,0 ≤ 9.0mag (see Sect. 2.1). Feldmeier-Krause et al. (2025b) fitted a two-dimensional multi-Gaussian expansion (MGE, Emsellem et al. 1994; Cappellari 2002) to the Gallego-Cano et al. (2020) surface density map, and we used this MGE to model the surface stellar density (see Table 1 in Feldmeier-Krause et al. 2025b). The MGE was scaled to the 4.5 μm band surface light distribution of Feldmeier-Krause et al. (2017b) in the centre.

In our DYNAMITE models, the deprojection of the surface stellar density MGE is determined by the space orientation angles (θ, ϕ, ψ), which can be converted analytically to three intrinsic shape parameters (q, p, u). These are ratios of the long, intermediate, and short axes a, b, and c of a triaxial system, with q = c/a, p = b/a, and u = a′/a, where a′ denotes the length of the longest axis a as projected on the sky. The flattest MGE component puts the strongest constraint on the deprojection, which means we have three free parameters (qmin, pmin, umin) in our models (see also Zhu et al. 2018c). We limited qmin to the range 0.1–0.2999 (the upper limit is given by the MGE), pmia to 0.6–0.9999, and umin to 0.8–1.0, with starting points at qmin = 0.29, pmin= 0.94, and umin = 0.99.

The deprojected stellar density distribution is multiplied by the dynamical mass-to-light ratio, ϒ. We assumed that this factor is spatially constant (see Feldmeier-Krause et al. 2025b, who found that varying ϒ is unnecessary), and used it to convert the stellar light density to a stellar mass density. We restricted its value to 0.1–2.0, with a starting value at ϒ = 1.0. The hyper-parameters qmin, pmin, umin, and ϒ thus describe the stellar gravitational potential.

Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

Stellar kinematic maps (left) after symmetrisation and their respective uncertainties (right). The plots show, from top to bottom, VLOS, σLOS, h3, and h4.

3.1.3 Dark matter

In a subset of models, we included a spherical DM component. The DM distribution is parameterised with a Navarro-Frenk-White (NFW Navarro et al. 1996) profile with a mass-concentration from Dutton & Macciò (2014). Thus, we are left with only one parameter, the dark matter fraction, f = M200/M* where M200 is the mass enclosed within a radius R200, the radius at which the average density is 200 times the critical density, and M* is the total stellar mass. The value of the hyperparameter log10 f is in the limits of −7 to +7 on a logarithmic scale, with a starting value at −3.

3.2 Constructing orbit solutions

For each gravitational potential, constrained by the hyperparameters qmin, pmin, umin, M• (and in a subset f), we numerically integrated a sample of orbits to form an orbit library, following the orbit-sampling scheme of van den Bosch et al. (2008). The hyperparameter Τ only scales the potential and does not require a separate orbit integration. As also M• is scaled by ϒ, the range of M• being modelled is ~(1.3–10) • 106 M⊙.

We have three orbit libraries per modelled potential: tube orbits, counter-rotating (CR) tube orbits, and box orbits. Each orbit library samples a combination of the three integrals of motion (E, I2, and I3, see van den Bosch et al. 2008), with a grid of nE×nI2×nI3=35×11×11Mathematical equation: ${n_E} \times {n_{{I_2}}} \times {n_{{I_3}}} = 35 \times 11 \times 11$ orbits. The energy grid is sampled using the relation of energy to a circular orbit, and the radii of the circular orbits are in the range of 100.5–103.86 arc-sec, corresponding to 0.127–291.5 pc. We note that the inner and outer radii correspond to ~0.3 × σmin and ~5 × σmax, respectively, where σmin and σmax denote the σ of the innermost and outermost MGE components in Sect. 3.1.2. We applied dithering with 33, such that we actually have orbit bundles, and in total 3×35×11×11×33= 343 035 orbits. The orbits were integrated for 200 periods, with a sampling of 50 000 points per orbit.

The orbit library weights were then fitted to reproduce the stellar LOSVD maps (V, σ, h3, h4) with the additional constraints that the deprojected 3D MGE stellar density distribution is reproduced within 1%, and the observed 2D MGE surface density distribution within 2% using a non-negative least squares (NNLS) minimisation. This is done for each set of hyperparameters. We computed the χ2 using the LOSVD maps, and found the best-fit model that minimises χ2. We refer to van den Bosch et al. (2008) and Zhu et al. (2018c) for further details on the modelling and fitting procedure.

We sampled the parameter space with the DYNAMITE option LEGACYGRIDSEARCH, which starts with a provided initial guess, then adds new models in later iterations. At each iteration, new models are seeded from existing models deemed acceptable based on a χ2 criterion. New models are seeded with gravitational potential parameters adjusted by some given step-size. The step-size decreases with later iterations to refine measurements. In total, we have computed >14 000 models without DM in 16 iterations, and >3 000 models with DM, in a narrower region of M• =(3.7–4.8) • 106 M⊙. We illustrate the sampled grid in Figs. A.1 and A.2.

For the 1σ confidence level, we select all models within a tolerance Δχ2 of the best-fitting model. To determine Δχ2, one can use Δχ2=2⋅nGH⋅NkinMathematical equation: $\Delta {\chi ^2} = \sqrt {2 \cdot {n_{GH}} \cdot {N_{kin}}} $, where nGH = 4 LOSVD parameters, and Nkin = 197 Voronoi bins (see Sect. 2.2). However, the validity of this level can be tested in a bootstrapping process (as in Tahmasebzadeh et al. 2024; Jin et al. 2025), to account for the model's numerical noise and non-uniqueness of the orbit weight distribution. We perturbed our kinematic data 1000 times, assuming Gaussian error distributions given by the uncertainty of each of the 4 × 197 measurements, and re-symmetrised the kinematic maps. Then we re-fit the orbit weights for these perturbed kinematic maps using the gravitational potential and orbit library of the best-fit model without DM. The resulting χ2 distribution has a Gaussian σ that is ~1.6×2⋅nGH⋅NkinMathematical equation: $\~1.6 \times \sqrt {2 \cdot {n_{GH}} \cdot {N_{kin}}} $, and we use this larger value as adjusted Δχ2 for the 1σ confidence level of the modelling.

4 Results

Our best-fit results are listed in Table 1. The models with and without DM obtain the same best-fit results for M•, ϒ, qmin, pmin, and umin. The value of the DM fraction, f, is only poorly constrained. The best-fit model (without DM) flux and kinematics are compared to the data in Fig. 4. The model and data agree very well, as indicated by the residual maps (right panel).

Table 1

Best-fit model parameters.

4.1 Shape and orientation

Our best-fit shape parameters qmin and umin prefer values near the edge of the model grid, at 0.2999 and 1.0, respectively. This indicates that the best-fit model is edge-on, and that the long axis is parallel to the Galactic plane. As the location of the Sun is relatively close to the Galactic mid-plane (~21 pc, Bennett & Bovy 2019), this finding suggests that the GC region we model is aligned with the larger MW, and has the same inclination. The short axis is parallel to Galactic latitude, b, and the intermediate axis is along the line of sight towards Sgr A⋆.

Our value of pmin =0.92 deviates from 1.0, and thus from perfect oblate axisymmetry. Given the orientation of the axes, p indicates compression along the line of sight. We show its radial behaviour in Fig. 5. However, pmin is closer to one than the result of Feldmeier-Krause et al. (2017b, pmin = 0.64), though the values agree within 1σ. The profile of q indicates increasing flattening towards larger radii, which is expected as the inner NSC is less flattened than the surrounding NSD.

The triaxiality parameter, T = (1 – p2)/(1 – q2), quantifies the deviation from axisymmetry. A value of T = 0 denotes oblate axisymmetry, T = 1 prolate axisymmetry. As shown in Fig. 5, we obtain T in the range of roughly 0.03–0.27, with the best-fit value at T ≈ 0.15. This brings the models in the range of being at most mildly triaxial (T = 0.1–0.3) according to the classification of Santucci et al. (2022).

4.2 Mass distribution

The value of M• matches the initial guess and true value of 4.3 • 106 M⊙, though we have large uncertainties. Nonetheless, this confirms the validity of our models. Our result is an improvement on the triaxial Schwarzschild models of Feldmeier-Krause et al. (2017b), who obtained M• = 3.0 • 106 M⊙, and the results agree within the uncertainties. Our value of ϒ =1.0 also matches the expectations. We scaled the central surface density to the surface brightness profile of Feldmeier-Krause et al. (2017b, ϒ = 0.9), and our results for ϒ are in agreement.

We show the mass distribution for the best-fit model without DM and the 1σ uncertainty in Fig 6, and list the enclosed mass at four different radii in Table 2. The best-fit model with DM results in a very similar total mass distribution as the model without DM. This is because the best-fit DM mass is more than five orders of magnitude lower than the total mass. For comparison, we also show the stellar mass profiles from other studies, and we will discuss the differences in Sect. 5.

Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Comparison of observations (left in rotated image), the best-fit model (middle, no DM), and the residuals (right). The rows show, from top to bottom: stellar flux, VLOS, σLOS, h3, and h4.
Thumbnail: Fig. 5 Refer to the following caption and surrounding text. Fig. 5

Intrinsic shape parameters and triaxiality as a function of radius r. The black lines denote q = c/a, blue lines p = b/a, and red lines triaxiality T = (1 – p2)/(1 – q2). Solid lines denote the best-fit model parameters, the dashed coloured lines their 1 σ uncertainties. The vertical solid line denotes 1 Re of the NSC, the dotted line the outer limit of the kinematic data. Horizontal lines mark steps of 0.1.
Table 2

Enclosed total mass for the model with no DM, total mass, and DM mass for the model with DM at different deprojected spherical radii.

4.3 Orbit circularity distribution

The circularity, λz=Lz¯/(r×Vc¯)Mathematical equation: ${\lambda _z} = \overline {{L_z}} /\left( {r \times \overline {{V_c}} } \right)$, indicates the orbit angular momentum, Lz¯Mathematical equation: $\overline {{L_z}} $, around the short axis, normalised by the angular momentum of a circular orbit of the same binding energy. Thus, λz = 1 presents a strongly rotating and dynamically cold short-axis tube orbit, and λz = −1 its retrograde counterpart, while λz = 0 is a dynamically hot box orbit. For each orbit, we computed the mean orbital radius, r¯Mathematical equation: ${\bar r}$, from the radii stored at each equal time step during orbit integration. The distribution of λz as a function of the mean orbital radius, r¯Mathematical equation: ${\bar r}$, of a model allows us to assess which region is dominated by dynamical cold, warm, or hot components.

We show our circularity distribution λz in Fig. 7. To create this plot, we used the orbit distributions of >800 models within the 1σ uncertainty limit and computed the entropy-regularised Wasserstein barycentre (Benamou et al. 2015; Flamary et al. 2021, 2024), which finds the transport-based middle of the 2-dimensional weight distributions. The weight distributions of the >800 models can differ by small shifts, and this method preserves the structure better, and produces sharper and physically more meaningful results than a simple average or median map.

We distinguish hot from warm orbits at λz = ±0.25, and warm from cold orbits at λz = ±0.8, using the same λz cuts as Santucci et al. (2022) and Thater et al. (2023). To better understand how the relative contribution of different orbit types changes with r¯Mathematical equation: ${\bar r}$, we show the relative weight of the cold, warm, hot, and CR orbits at a given mean orbital radius r¯Mathematical equation: ${\bar r}$ in Fig. 8, the coloured bands show the 1σ percentile ranges of >800 models within the 1σ uncertainty limit.

In the region dominated by the NSC (r¯≲7 pc)Mathematical equation: $\left( {\bar r \mathbin{\lower.3ex\hbox{$\buildrel<\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 7{\rm{ pc}}} \right)$, the models reveal that the largest contribution is from warm orbits, followed by hot, CR, and cold orbits. At a mean orbital radius of 5–15 pc, that is 1–3 Re of the NSC, the contribution of hot orbits increases from 0.3 to ~0.6, whereas the weight fraction of warm orbits decreases from 0.4 to ~0.25. These trends revert again at ~20 pc, when the warm orbits' weight fraction starts to increase, while the hot orbit weight fraction decreases. The CR orbits decrease from ~0.2 within the NSC to ~0.1 in the inner NSD. Cold rotating orbits contribute ~0.1 at these radii.

Our 1σ models have large fractions of cold orbits and CR orbits in the outer NSD, though their exact location varies between individual models (r¯~80−140 pcMathematical equation: $\bar r\~80 - 140{\rm{ pc}}$). We note that the field of view (FOV) of the kinematic data extends only to l = 33 pc. Orbits with r¯>33pcMathematical equation: $\bar r > 33{\rm{pc}}$ pc spend at least half of the time at intrinsic radii beyond the extent of the data. But stars on such orbits can cross our FOV if their projected radius is smaller. The presence of these cold orbits is constrained by our models, though the exact value of r¯Mathematical equation: ${\bar r}$ is not.

Another way to differentiate hot from cold orbits, but without distinction in co-rotating and CR orbits, is via Rmax vs zmax plots. We show examples of such plots using the best-fit model in Appendix B to enable a comparison with Nieuwmunster et al. (2024), who integrated orbits of NSD stars in a fixed gravitational potential.

Thumbnail: Fig. 6 Refer to the following caption and surrounding text. Fig. 6

Total enclosed mass as a function of spherical deprojected radius. Shaded regions show 1σ uncertainties. The vertical solid line denotes 1 Re of the NSC, the dotted line the outer limit of the kinematic data.
Thumbnail: Fig. 7 Refer to the following caption and surrounding text. Fig. 7

Orbit circularity distribution of λz as a function of mean orbital radius r¯Mathematical equation: ${\bar r}$, computed using over >800 models within 1σ. Colour indicates the orbit density in the phase space, horizontal dashed lines divide the orbits into cold (λz>0.8), warm (0.25<λz ≤0.8), hot (−0.25<λz<0.25), CR warm (−0.8 < λz ≤ −0.2), and CR cold (λz ≤ −0.8) orbits, vertical lines are as in Fig. 5.
Thumbnail: Fig. 8 Refer to the following caption and surrounding text. Fig. 8

Relative orbit weight profile as function of mean orbital radius r¯Mathematical equation: ${\bar r}$, computed from the λz distribution in Fig. 7. The different colours denote different orbit types, the shaded regions the 1σ percentiles of the 800 models. Blue denotes cold, orange warm, red hot, and cyan CR orbits, vertical lines are as in Fig. 5.

4.4 Best-fit model orbit decomposition

Here, we decompose the stellar orbits of our best-fit model. Using the orbit weight distribution and λz, we show the spatial and kinematic distributions of different orbital types. This method has been applied to various galaxies in (Zhu et al. (2018c,a); Jin et al. (2020); Santucci et al. (2022); Breda et al. (2026)), but barely on NSCs (but see Lamprecht et al. 2026).

After grouping the orbits according to their we computed the surface brightness, VLOS, and σLOS maps that the groups contribute to the total model. These maps are shown in Fig. 9. They take into account the contribution of all orbits that pass our FOV. We see that CR warm orbits have a steep decline in the surface brightness; they contribute most in the inner ±2 pc (~10%), and little further out (≲5%). Cold orbits and hot orbits have the highest overall surface brightness contributions. Even though cold orbits appear sparse in the inner circularity plot (Figs. 7–8), they are important at large r¯Mathematical equation: ${\bar r}$. These orbits cross the FOV, so that they contribute significantly to the surface brightness of the model. Except for the inner ≲2pc, where warm and hot orbits dominate, cold orbits have the highest relative weight to the surface brightness maps.

The VLOS maps reveal that the cold and warm orbits have VLOS values up to approximately ±80 km s−1, but as there are also CR orbits and hot box orbits with lower absolute VLOS (≲ ± 10 km s−1 ), the VLOS value of all orbits combined adds up to only ±40 km s−1. Interestingly, the CR cold orbits have the highest VLOS peaks in the inner ~1 pc with about ±125 km s−1. The σLOS maps show that cold orbits have lower than warm orbits, and hot orbits have the highest σLOS of all orbit types, as expected.

5 Discussion

5.1 Comparison to dynamical models in the literature

Our best-fit model fits the data very well. Given our kinematic uncertainties and the flexibility of the Schwarzschild models, we found many models that fit similarly well. Subsequently, the parameter uncertainties are large, in particular, compared to Feldmeier-Krause et al. (2017b), who used an older version of the DYNAMITE code and constrained the models with different kinematic data and a different stellar MGE distribution. However, this is because we use a more conservative Δχ2 cut to define the parameter uncertainties than Feldmeier-Krause et al. (2017b) did, who used Δχ2=5.9 and 18.2 as 1σ and 3σ limits. With these bounds (see Table 1), our uncertainties are much lower and similar to the uncertainties of Feldmeier-Krause et al. (2017b).

Nonetheless, we can constrain the shape and orbital distributions. The intrinsic shape parameters, in particular q, have only a rather small uncertainty (see Fig. 5), and agree very well with the flattening of the stellar density measured by Gallego-Cano et al. (2020), that is 0.71±0.10 for the NSC and 0.338±0.002 for the NSD. We find that the value of p, the intermediate-to-major axis ratio, is relatively constant with radius, and less compressed than found by Feldmeier-Krause et al. (2017b, pmin = 0.64) in the inner ~6 pc. Feldmeier-Krause et al. (2017b) speculate that interstellar dust within the GC dominantly affects the stars further away along the LOS, and that the integrated light is biased to the near side of the GC. A possible consequence is that the GC appears compressed along the LOS, and the value of pmin may be underestimated. Our FOV is larger and extends to less reddened regions, and our spectroscopic observations are deeper (we have integrated five times as long). These two factors may, to some extent, mitigate this LOS bias in our data. The anisotropy profiles (Fig. C.1) show a remarkable agreement, with a minimum value within r ≲ 3 pc, and close to isotropy at r ≳ 4 pc.

We compare the stellar mass profiles from Feldmeier-Krause et al. (2017b, binned triaxial orbit-based models), Sormani et al. (2022, discrete axisymmetric distribution function models), and Feldmeier-Krause et al. (2025b, discrete axisymmetric Jeans models) in Fig. 6. The stellar mass profiles agree mostly within the uncertainties, though they deviate beyond our uncertainties at r ≳ 45 pc, where we have no kinematic data. The data of Feldmeier-Krause et al. (2017b) are limited to l ≲ 6 pc. The models of Sormani et al. (2022) are focused on the more extended NSD. At l ≲ 20 pc, their sample is rather small, as their data lies mostly at larger distances. We used the same light distribution as Feldmeier-Krause et al. (2025b), but our mass distribution has larger uncertainties, by about a factor of 10. One reason for this is that Jeans models have less flexibility in comparison to the orbit-based modelling, and hence unrealistically small uncertainties. Feldmeier-Krause et al. (2025b) had discrete data rather than binned data, which was more extended (covering a circle of r ≲ 9.5pcinadditiontodatainour FOV),andpropermotionsfor ~75% of the 4 600 stars with Vlos, leading to tighter constraints on the enclosed mass and mass-to-light ratio (ϒ=0.75±0.02). Further, the models of Feldmeier-Krause et al. (2025b) have a separate component for a background contribution of bar stars, which we do not have. This may bias our velocity dispersion to higher values, especially at the outer region of our kinematic data, as the bar contribution becomes more important at larger radii (Feldmeier-Krause et al. 2025a). This may cause a higher ϒ and thus mass estimate.

Most of the previously mentioned works on the extended GC mass distribution neglected DM. Feldmeier-Krause et al. (2025b) included a DM component to a subset of their models. Like us, they obtained large uncertainties. An extrapolation of the DM volume density to ≳100 pc allows for a comparison with the simulations of Hussein et al. (2025) and the bar models of Portail et al. (2017). While the Feldmeier-Krause et al. (2025b) models lie mostly above the values of Portail et al. (2017) and Hussein et al. (2025), our DM volume density tends to be lower, yet both agree within the uncertainties. To be more quantitative, the Feldmeier-Krause et al. (2025b) DM mass contribution to the total enclosed mass lies at 6–25% at 33 pc, whereas we obtain only ≲0.8%. The true value lies probably in between.

Thumbnail: Fig. 9 Refer to the following caption and surrounding text. Fig. 9

Orbit decomposition of the best-fit model in cold, warm, counter-rotating (cr) cold, cr warm, hot, and all orbits. The rows show, from left to right: stellar surface brightness, Vlos, σLOS.

5.2 Circularity distribution and orbit decomposition

Our circularity distribution shows contributions from cold, warm, hot, and CR orbits. The circularity contains information on the origin of the stars. Boecker et al. (2023) studied the circularity of the inner 500 pc3 of TNG50 galaxies and found that, at the stellar mass range of the MW (a few times 1010 M⊙), stars that migrated to the galaxy centre have on average more rotational support (colder orbits, λz ~ 0.5) compared to in situ stars, which show more random motion (λz ~ 0.25). Ex situ stars, that is, stars accreted from galaxy mergers, are on average on hot orbits, but there is a large galaxy-to-galaxy variation. The circularity depends on the time of the merger and the orbit configuration between the host and merging satellite, and can result in CR orbits. Boecker et al. (2023) also found that stars on colder orbits are, on average, younger, and the orbits become hotter due to dynamical heating over time (see also Breda et al. 2024). The cold orbits we found at the outer NSD (r¯≳80 pcMathematical equation: $\bar r \mathbin{\lower.3ex\hbox{$\buildrel>\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 80{\rm{ pc}}$) may indeed be populated by younger stars, as the inside-out formation scenario for NSDs (Bittner et al. 2020; Nogueras-Lara et al. 2023) suggests. We note that most of the giant molecular gas structures of the central molecular zone are located within ~ 100 pc (except for the 1.3° cloud complex, Henshaw et al. 2016). However, as our data covers only a narrow range of latitude (≲3 pc), the circularity of cold orbits at r¯≳80 pcMathematical equation: $\bar r \mathbin{\lower.3ex\hbox{$\buildrel>\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 80{\rm{ pc}}$ may be overestimated by our models, and the orbits may be warmer. Nonetheless, we note that Nieuwmunster et al. (2024), who integrated the orbits of >1 100 stars at l≲210 pc in the NSD, also found a large fraction of tube orbits (≳65%).

In general, dynamical heating makes a dynamically cold, thin disc thicker and hotter. It drives the system towards isotropy, erasing signatures of the ex situ origin of stars. The heating mechanisms in the GC are two-body relaxation (among stars or with stellar remnants, Alexander 2005), and massive perturbers (e.g. giant molecular clouds or star clusters, Perets et al. 2007). While the former is dominant in the inner parsec, the latter is dominant at r ≳1.5 to ~100 pc (Perets et al. 2007). The Galactic bar further influences stellar orbits and redistributes angular momentum within the Galaxy (Athanassoula 1992). Our finding that the inner regions are dominated by warm and hot orbits rather than cold orbits is in agreement with these regions being older and having experienced more dynamical heating.

Counter-rotating orbits may trace infalling star clusters. Tsatsi et al. (2017) studied the consecutive infall of star clusters to the GC and found that these can produce kinematic substructures if they come from random orbital directions. Though such signatures will probably be washed out by dynamical relaxation in a few gigayears (Arca Sedda et al. 2020). Ishchenko et al. (2023) found that there may be up to 3–4 close passages of globular clusters with the GC per 1 Gyr. During these interactions, the star clusters can lose some of their stars, which then leave a dynamical imprint in the circularity distribution.

Some of the hot orbits in our model may be chaotic orbits. Penoyre et al. (2025) computed the fraction of chaotic orbits among box orbits (λz = 0) as a function of radius and NSC flattening, and found values up to 40% for q=0.7. Assuming this fraction remains constant at r > 10 pc, we would have an overall fraction of ~15% chaotic orbits in our FOV. Galactic bar orbits are not explicitly included in our models. In the inner 100 pc, we would expect some contribution from x1, x2, and x3 orbits, extended along and perpendicular to the bar length. Though the overall contribution of the bar to the stellar surface density in the region of our data is ≲10% at 30 pc (Feldmeier-Krause et al. 2025a), and even less further in Zhu et al. (2018b) showed that triaxial orbit-based models that do not explicitly include a bar structure still reproduce the orbit distributions of simulated barred galaxies without strong biases. We conclude that this omission has no strong effect on our orbit distribution.

There have been speculations of a nuclear bar inside the GC (Alard 2001; Gonzalez et al. 2011; Fiteni et al. 2026), but no clear detection yet. Our models also show no sign of a nuclear bar. Bars produce diagonal lines in the λζ distribution plots (Tahmasebzadeh et al. 2024), which we did not obtain. We will investigate including bar orbits in orbit-based modelling of the GC in the future.

6 Summary and conclusions

We have computed a large suite of triaxial orbit-based dynamical models of the GC, constrained by the integrated line-of-sight kinematics of red giant stars. The kinematic data extend out to 33 pc from Sgr A⋆ towards the Galactic east and west, and 1–2 pc towards the north and south, and thus from the NSC to the inner part of the NSD.

Our models recover the mass of the central black hole correctly, though with large uncertainties: M•=(4.3−2.9+4.1)×106M⊙Mathematical equation: ${M_ \bullet } = \left( {4.3_{ - 2.9}^{ + 4.1}} \right) \times {10^6}{M_ \odot }$.

Our stellar mass agrees with axisymmetric studies, though we tend to obtain higher masses at r>30pc. Including a dark matter component does not alter the other best-fit parameters, as we find only a very low DM contribution: ≲1% of the total mass at 33 pc. Overall, our results are in good agreement with studies in the literature that use resolved kinematic data. This validates that orbit-based dynamical models with integrated light data, which are commonly used in extragalactic studies, produce robust results.

We find only mild triaxiality across the modelled region. The orbit circularity as a function of mean orbital radius shows mostly hot and warm orbits in the NSC and inner NSD. We detect a cold component at the outer part of the NSD (r¯~80−160 pcMathematical equation: $\bar r\~80 - 160{\rm{ pc}}$), though its exact location varies among 1σ models. This component has experienced less dynamical heating and may possibly be younger. We also detect some CR orbits, which may correspond to infalling star clusters.

To make the best use of orbit-based modelling in the GC, it would be ideal to fit the models to discrete velocity measurements rather than the integrated light maps and to include proper motion data. Such models have been attempted for axisymmetric (Chanamé et al. 2008) and spherical systems (Magorrian 2019), but with only a few applications. In the future, the DYNAMITE code will be extended to include proper motion and discrete data, and the GC is an ideal testbed for this code.

Acknowledgements

We thank the anonymous referee for their review and constructive comments. A.F.K. acknowledges funding from the Austrian Science Fund (FWF) [grant DOI 10.55776/ESP542]. I.B. has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Grant agreement ID no. 101059532. This project was extended for 6 months by the Franziska Seidl Funding Program of the University of Vienna. I.B. was supported by Fundaçâo para a Ciência e a Tecnologia (FCT) through national funds under the research grant UID/04434/2025 (DOI 10.54499/UID/04434/2025). We thank the Gemini Observatory staff for their support during the planning and execution of the observations and for their advice on data reduction. We thank David Rupke for providing a general-purpose library for IFU data cubes. The computational results have been achieved using the Austrian Scientific Computing (ASC) infrastructure. Based on observations obtained at the international Gemini Observatory, a program of NSF NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the U.S. National Science Foundation on behalf of the Gemini Observatory partnership: the U.S. National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea). Data were processed using the Gemini IRAF package. This research made use of Montage. It is funded by the National Science Foundation under Grant Number ACI-1440620, and was previously funded by the National Aeronautics and Space Administration's Earth Science Technology Office, Computation Technologies Project, under Cooperative Agreement Number NCC5-626 between NASA and the California Institute of Technology. This research has made use of NASA's Astrophysics Data System; the IPython package (Pérez & Granger 2007); Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2018, 2013); SciPy (Virtanen et al. 2020); matplotlib, a Python library for publication quality graphics (Hunter 2007); NumPy (Harris et al. 2020); ds9, a tool for data visualization supported by the Chandra X-ray Science Center (CXC) and the High Energy Astrophysics Science Archive Center (HEASARC) with support from the JWST Mission office at the Space Telescope Science Institute for 3D visualization.

References

  1. Alard, C. 2001, A&A, 379, L44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Alexander, T. 2005, Phys. Rep., 419, 65 [CrossRef] [Google Scholar]
  3. Arca Sedda, M., Gualandris, A., Do, T., et al. 2020, ApJ, 901, L29 [Google Scholar]
  4. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  6. Athanassoula, E. 1992, MNRAS, 259, 328 [NASA ADS] [CrossRef] [Google Scholar]
  7. Barth, A. J., Strigari, L. E., Bentz, M. C., Greene, J. E., & Ho, L. C. 2009, ApJ, 690, 1031 [NASA ADS] [CrossRef] [Google Scholar]
  8. Benamou, J.-D., Carlier, G., Cuturi, M., Nenna, L., & Peyré, G. 2015, SIAM J. Sci. Comput., 37, A1111 [CrossRef] [Google Scholar]
  9. Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bittner, A., Sánchez-Blázquez, P., Gadotti, D. A., et al. 2020, A&A, 643, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Blum, R. D., Sellgren, K., & Depoy, D. L. 1996, AJ, 112, 1988 [Google Scholar]
  12. Boecker, A., Neumayer, N., Pillepich, A., et al. 2023, MNRAS, 519, 5202 [NASA ADS] [CrossRef] [Google Scholar]
  13. Boehle, A., Ghez, A. M., Schödel, R., et al. 2016, ApJ, 830, 17 [Google Scholar]
  14. Breda, I., van de Ven, G., Thater, S., et al. 2024, A&A, 692, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Breda, I., van de Ven, G., Thater, S., et al. 2026, 709, A262 [Google Scholar]
  16. Cappellari, M. 2002, MNRAS, 333, 400 [NASA ADS] [CrossRef] [Google Scholar]
  17. Cappellari, M. 2023, MNRAS, 526, 3273 [NASA ADS] [CrossRef] [Google Scholar]
  18. Cappellari, M., & Copin, Y. 2003, MNRAS, 342, 345 [Google Scholar]
  19. Cappellari, M., & Emsellem, E. 2004, PASP, 116, 138 [Google Scholar]
  20. Cappellari, M., Verolme, E. K., van der Marel, R. P., et al. 2002, ApJ, 578, 787 [NASA ADS] [CrossRef] [Google Scholar]
  21. Chanamé, J., Kleyna, J., & van der Marel, R. 2008, ApJ, 682, 841 [Google Scholar]
  22. Chatzopoulos, S., Fritz, T. K., Gerhard, O., et al. 2015, MNRAS, 447, 948 [Google Scholar]
  23. Conroy, C., Villaume, A., van Dokkum, P. G., & Lind, K. 2018, ApJ, 854, 139 [Google Scholar]
  24. Davidge, T. J. 2020, AJ, 160, 146 [Google Scholar]
  25. den Brok, M., Krajnović, D., Emsellem, E., Brinchmann, J., & Maseda, M. 2021, MNRAS, 508, 4786 [NASA ADS] [CrossRef] [Google Scholar]
  26. Do, T., Hees, A., Ghez, A., et al. 2019, Science, 365, 664 [Google Scholar]
  27. Dutton, A. A., & Macciò, A. V. 2014, MNRAS, 441, 3359 [Google Scholar]
  28. Eikenberry, S. S., Elston, R., Raines, S. N., et al. 2004, SPIE Conf. Ser., 5492, 1196 [NASA ADS] [Google Scholar]
  29. Emsellem, E., Monnet, G., & Bacon, R. 1994, A&A, 285, 723 [NASA ADS] [Google Scholar]
  30. Fahrion, K., Lyubenova, M., van de Ven, G., et al. 2019, A&A, 628, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Falcón-Barroso, J., & Martig, M. 2021, A&A, 646, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Feldmeier, A., Neumayer, N., Seth, A., et al. 2014, A&A, 570, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Feldmeier-Krause, A. 2022, MNRAS, 513, 5920 [Google Scholar]
  34. Feldmeier-Krause, A., Kerzendorf, W., Neumayer, N., et al. 2017a, MNRAS, 464, 194 [Google Scholar]
  35. Feldmeier-Krause, A., Zhu, L., Neumayer, N., et al. 2017b, MNRAS, 466, 4040 [NASA ADS] [Google Scholar]
  36. Feldmeier-Krause, A., Neumayer, N., Seth, A., et al. 2025a, A&A, 696, A213 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Feldmeier-Krause, A., Veršič, T., van de Ven, G., Gallego-Cano, E., & Neumayer, N. 2025b, A&A, 699, A239 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Ferrarese, L., & Merritt, D. 2000, ApJ, 539, L9 [Google Scholar]
  39. Ferrarese, L., Côté, P., Dalla Bontà, E., et al. 2006, ApJ, 644, L21 [Google Scholar]
  40. Fiteni, K., Li, X., Sormani, M. C., et al. 2026, A&A, 709, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Flamary, R., Courty, N., Gramfort, A., et al. 2021, J. Mach. Learn. Res., 22, 1 [Google Scholar]
  42. Flamary, R., Vincent-Cuaz, C., Courty, N., et al. 2024, POT Python Optimal Transport (version 0.9.5) [Google Scholar]
  43. Fritz, T. K., Gillessen, S., Dodds-Eden, K., et al. 2011, ApJ, 737, 73 [Google Scholar]
  44. Fritz, T. K., Chatzopoulos, S., Gerhard, O., et al. 2016, ApJ, 821, 44 [Google Scholar]
  45. Gallego-Cano, E., Schödel, R., Nogueras-Lara, F., et al. 2020, A&A, 634, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJ, 539, L13 [Google Scholar]
  47. Georgiev, I. Y., Böker, T., Leigh, N., Lützgendorf, N., & Neumayer, N. 2016, MNRAS, 457, 2122 [Google Scholar]
  48. Gillessen, S., Plewa, P. M., Eisenhauer, F., et al. 2017, ApJ, 837, 30 [Google Scholar]
  49. Gonzalez, O. A., Rejkuba, M., Minniti, D., et al. 2011, A&A, 534, L14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. GRAVITY Collaboration (Abuter, R., et al.) 2022, A&A, 657, L12 [NASA ADS] [CrossRef] [Google Scholar]
  51. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [NASA ADS] [CrossRef] [Google Scholar]
  52. Henshaw, J. D., Longmore, S. N., Kruijssen, J. M. D., et al. 2016, MNRAS, 457, 2675 [Google Scholar]
  53. Hermosa Muñoz, L., Márquez, I., Cazzoli, S., Masegosa, J., & Agís-González, B. 2022, A&A, 660, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  55. Hussein, A., Necib, L., Kaplinghat, M., et al. 2025, arXiv e-prints [arXiv:2581.14868] [Google Scholar]
  56. Ishchenko, M., Sobolenko, M., Kuvatova, D., Panamarev, T., & Berczik, P. 2023, A&A, 674, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Jeans, J. H. 1922, MNRAS, 82, 122 [NASA ADS] [CrossRef] [Google Scholar]
  58. 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]
  59. Jin, Y., Zhu, L., Long, R. J., et al. 2019, MNRAS, 486, 4753 [Google Scholar]
  60. Jin, Y., Zhu, L., Long, R. J., et al. 2020, MNRAS, 491, 1690 [NASA ADS] [Google Scholar]
  61. Jin, Y., Zhu, L., Tahmasebzadeh, B., et al. 2025, A&A, 700, A249 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511 [NASA ADS] [CrossRef] [Google Scholar]
  63. Lamprecht, J., Feldmeier-Krause, A., Lyubenova, M., et al. 2026, A&A, 706, A373 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Launhardt, R., Zylka, R., & Mezger, P. G. 2002, A&A, 384, 112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Lyubenova, M., van den Bosch, R. C. E., Côté, P., et al. 2013, MNRAS, 431, 3364 [Google Scholar]
  66. Magorrian, J. 2019, MNRAS, 484, 1166 [NASA ADS] [CrossRef] [Google Scholar]
  67. Mehrgan, K., Thomas, J., Saglia, R., et al. 2019, ApJ, 887, 195 [NASA ADS] [CrossRef] [Google Scholar]
  68. Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [Google Scholar]
  69. Neumayer, N., Seth, A., & Böker, T. 2020, A&A Rev., 28, 4 [Google Scholar]
  70. Neureiter, B., Thomas, J., Saglia, R., et al. 2021, MNRAS, 500, 1437 [Google Scholar]
  71. Neureiter, B., de Nicola, S., Thomas, J., et al. 2023, MNRAS, 519, 2004 [Google Scholar]
  72. Nguyen, D. D., Seth, A. C., Neumayer, N., et al. 2018, ApJ, 858, 118 [NASA ADS] [CrossRef] [Google Scholar]
  73. Nieuwmunster, N., Schultheis, M., Sormani, M., et al. 2024, A&A, 685, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Nogueras-Lara, F., Gallego-Calvente, A. T., Dong, H., et al. 2018, A&A, 610, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Nogueras-Lara, F., Schödel, R., Gallego-Calvente, A. T., et al. 2019, A&A, 631, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Nogueras-Lara, F., Schultheis, M., Najarro, F., et al. 2023, A&A, 671, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Penoyre, Z., Rossi, E. M., & Stone, N. C. 2025, MNRAS, 542, 322 [Google Scholar]
  78. Perets, H. B., Hopman, C., & Alexander, T. 2007, ApJ, 656, 709 [CrossRef] [Google Scholar]
  79. Pérez, F., & Granger, B. E. 2007, Comput. Sci. Eng., 9, 21 [Google Scholar]
  80. Philipp, S., Zylka, R., Mezger, P. G., et al. 1999, A&A, 348, 768 [NASA ADS] [Google Scholar]
  81. Poci, A., McDermid, R. M., Zhu, L., & van de Ven, G. 2019, MNRAS, 487, 3776 [Google Scholar]
  82. Portail, M., Gerhard, O., Wegg, C., & Ness, M. 2017, MNRAS, 465, 1621 [NASA ADS] [CrossRef] [Google Scholar]
  83. Saito, R. K., Hempel, M., Minniti, D., et al. 2012, A&A, 537, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Santucci, G., Brough, S., van de Sande, J., et al. 2022, ApJ, 930, 153 [NASA ADS] [CrossRef] [Google Scholar]
  85. Schödel, R., Merritt, D., & Eckart, A. 2009, A&A, 502, 91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Schwarzschild, M. 1979, ApJ, 232, 236 [NASA ADS] [CrossRef] [Google Scholar]
  87. Scott, N., & Graham, A. W. 2013, ApJ, 763, 76 [NASA ADS] [CrossRef] [Google Scholar]
  88. Sormani, M. C., Sanders, J. L., Fritz, T. K., et al. 2022, MNRAS, 512, 1857 [CrossRef] [Google Scholar]
  89. Tahmasebzadeh, B., Zhu, L., Shen, J., et al. 2024, MNRAS, 534, 861 [NASA ADS] [CrossRef] [Google Scholar]
  90. Thater, S., Jethwa, P., Tahmasebzadeh, B., et al. 2022a, A&A, 667, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Thater, S., Krajnović, D., Weilbacher, P. M., et al. 2022b, MNRAS, 509, 5416 [Google Scholar]
  92. Thater, S., Lyubenova, M., Fahrion, K., et al. 2023, A&A, 675, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Thater, S., Chaturvedi, A., Krajnović, D., et al. 2026, A&A, accepted [arXiv:2685.28959] [Google Scholar]
  94. Tsatsi, A., Mastrobuono-Battisti, A., van de Ven, G., et al. 2017, MNRAS, 464, 3720 [Google Scholar]
  95. van den Bosch, R. C. E., & de Zeeuw, P. T. 2010, MNRAS, 401, 1770 [NASA ADS] [CrossRef] [Google Scholar]
  96. 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]
  97. Vasiliev, E. 2013, MNRAS, 434, 3174 [Google Scholar]
  98. Vasiliev, E., & Valluri, M. 2020, ApJ, 889, 39 [Google Scholar]
  99. Vasiliev, E., Feldmeier-Krause, A., & Sormani, M. C. 2026, ApJ, 1002, 71 [Google Scholar]
  100. Vazdekis, A., Koleva, M., Ricciardelli, E., Röck, B., & Falcón-Barroso, J. 2016, MNRAS, 463, 3409 [Google Scholar]
  101. Verro, K., Trager, S. C., Peletier, R. F., et al. 2022, A&A, 661, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  102. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  103. Wallace, L., & Hinkle, K. 1996, ApJS, 107, 312 [NASA ADS] [CrossRef] [Google Scholar]
  104. Walsh, J. L., van den Bosch, R. C. E., Barth, A. J., & Sarzi, M. 2012, ApJ, 753, 79 [NASA ADS] [CrossRef] [Google Scholar]
  105. Zhu, L., van de Ven, G., Méndez-Abreu, J., & Obreja, A. 2018a, MNRAS, 479, 945 [Google Scholar]
  106. Zhu, L., van de Ven, G., van den Bosch, R., et al. 2018b, Nat. Astron., 2, 233 [Google Scholar]
  107. Zhu, L., van den Bosch, R., van de Ven, G., et al. 2018c, MNRAS, 473, 3000 [Google Scholar]
  108. Zhu, L., van de Ven, G., Leaman, R., et al. 2020, MNRAS, 496, 1579 [Google Scholar]
3

σt=σϕ2+σθ2Mathematical equation: ${\sigma _t} = \sqrt {\sigma _\phi ^2 + \sigma _\theta ^2} $

Appendix A Corner plots

We show the parameter space spanned by our models and the location of the best-fitting model in Fig. A.1 (14000 models without DM) and in Fig. A.2 (3 000 models with DM). The shape parameter qmin is near the edge of the grid, and umin is close to 1, indicating that near edge-on models are preferred, with the major axis being parallel to the Galactic plane.

Appendix B Rmax versus zmax plots

In Figs. B.1 and B.2, we show the distribution of orbit weights of the best-fit model in Rmax vs zmax space. Rmax denotes the maximum of the orbits' intrinsic radii, while zmax denotes the orbits' maximum distance to the major axis. Hot orbits have higher values of zmax than cold orbits. The orbits in the inner region are finer sampled, as indicated by the higher density of points in the inner ≲33 pc region.

Appendix C Velocity anisotropy

The intrinsic velocity anisotropy profile ßr is shown in Fig. C.1. Anisotropy is defined as βr=1−σt2/σr2Mathematical equation: ${\beta _r} = 1 - \sigma _t^2/\sigma _r^2$ (σr, σϕ,σθ are the radial, azimuthal angular and polar angular velocity dispersion in spherical co-ordinates3), plotted along intrinsic radius, r. A value of zero denotes an isotropic velocity distribution; a higher (lower) value than zero denotes radially (tangentially) anisotropic orbits. Radial anisotropy is a sign of dynamically hot orbits, and tangential anisotropy indicates more circular orbits. Our models indicate some tangential anisotropy within the inner ~3 pc ofthe NSC, but are in agreement with isotropy for >10pc.

Thumbnail: Fig. A.1 Refer to the following caption and surrounding text. Fig. A.1

Illustration of the sampled parameter space, each point denotes one of >14 000 models. The symbol size and colour represent the value of χ2, as indicated by the colour bar; the black x-symbol denotes the best-fitting model. The parameters are, from top to bottom: umin, pmin, qmin, M•, and ϒ.
Thumbnail: Fig. A.2 Refer to the following caption and surrounding text. Fig. A.2

Same as Fig. A.1, but for models with a spherical DM component and a more narrow range of M•. The parameters are, from top to bottom: umin, pmin, qmin, f, M•, and ϒ.
Thumbnail: Fig. B.1 Refer to the following caption and surrounding text. Fig. B.1

Rmax vs. zmax diagram of the best-fit model. Each data point represents an orbit with non-zero weight, the colour represents the relative weight of the orbit, the solid line denotes 1 Re of the NSC, and the dashed line denotes the outer extent of the kinematic data. The panels show, from top to bottom, regular cold and warm tube orbits (λz>0.25), CR orbits (λz ≤-0.25), hot orbits, and all orbits.
Thumbnail: Fig. B.2 Refer to the following caption and surrounding text. Fig. B.2

Same as Fig. B.1, but zoomed into the region dominated by the NSC. Note the different colour scale to improve visibility of a range of orbit weights.
Thumbnail: Fig. C.1 Refer to the following caption and surrounding text. Fig. C.1

Radial velocity anisotropy ßr as function of intrinsic radius. The black line denotes the best-fit model without DM, the shaded region the 1σ uncertainty, the vertical solid line 1 Re of the NSC, the vertical dotted line the outer limit of the kinematic data.

All Tables

Table 1

Best-fit model parameters.

In the text
Table 2

Enclosed total mass for the model with no DM, total mass, and DM mass for the model with DM at different deprojected spherical radii.

In the text

All Figures

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Spatial coverage of our F2 spectroscopic data. The data extend ~66 pc along the Galactic longitude l, centred on Sgr A⋆ (marked as a red plus symbol), and ~1 pc to the Galactic north and south, except for the centre region, which extends further to the Galactic north (~2pc). The image was constructed from the spectroscopic scans. We show the Galactic co-ordinate grid as dashed lines.
In the text
Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Example of a stellar kinematic fit with PPXF. The black line shows the data, the red line the best-fit stellar model, and the orange and pink lines the best-fit gas emission. Green symbols denote residuals; the shaded grey regions were masked in the fit due to sky residuals or bad pixels.
In the text
Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

Stellar kinematic maps (left) after symmetrisation and their respective uncertainties (right). The plots show, from top to bottom, VLOS, σLOS, h3, and h4.
In the text
Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Comparison of observations (left in rotated image), the best-fit model (middle, no DM), and the residuals (right). The rows show, from top to bottom: stellar flux, VLOS, σLOS, h3, and h4.
In the text
Thumbnail: Fig. 5 Refer to the following caption and surrounding text. Fig. 5

Intrinsic shape parameters and triaxiality as a function of radius r. The black lines denote q = c/a, blue lines p = b/a, and red lines triaxiality T = (1 – p2)/(1 – q2). Solid lines denote the best-fit model parameters, the dashed coloured lines their 1 σ uncertainties. The vertical solid line denotes 1 Re of the NSC, the dotted line the outer limit of the kinematic data. Horizontal lines mark steps of 0.1.
In the text
Thumbnail: Fig. 6 Refer to the following caption and surrounding text. Fig. 6

Total enclosed mass as a function of spherical deprojected radius. Shaded regions show 1σ uncertainties. The vertical solid line denotes 1 Re of the NSC, the dotted line the outer limit of the kinematic data.
In the text
Thumbnail: Fig. 7 Refer to the following caption and surrounding text. Fig. 7

Orbit circularity distribution of λz as a function of mean orbital radius r¯Mathematical equation: ${\bar r}$, computed using over >800 models within 1σ. Colour indicates the orbit density in the phase space, horizontal dashed lines divide the orbits into cold (λz>0.8), warm (0.25<λz ≤0.8), hot (−0.25<λz<0.25), CR warm (−0.8 < λz ≤ −0.2), and CR cold (λz ≤ −0.8) orbits, vertical lines are as in Fig. 5.
In the text
Thumbnail: Fig. 8 Refer to the following caption and surrounding text. Fig. 8

Relative orbit weight profile as function of mean orbital radius r¯Mathematical equation: ${\bar r}$, computed from the λz distribution in Fig. 7. The different colours denote different orbit types, the shaded regions the 1σ percentiles of the 800 models. Blue denotes cold, orange warm, red hot, and cyan CR orbits, vertical lines are as in Fig. 5.
In the text
Thumbnail: Fig. 9 Refer to the following caption and surrounding text. Fig. 9

Orbit decomposition of the best-fit model in cold, warm, counter-rotating (cr) cold, cr warm, hot, and all orbits. The rows show, from left to right: stellar surface brightness, Vlos, σLOS.
In the text
Thumbnail: Fig. A.1 Refer to the following caption and surrounding text. Fig. A.1

Illustration of the sampled parameter space, each point denotes one of >14 000 models. The symbol size and colour represent the value of χ2, as indicated by the colour bar; the black x-symbol denotes the best-fitting model. The parameters are, from top to bottom: umin, pmin, qmin, M•, and ϒ.
In the text
Thumbnail: Fig. A.2 Refer to the following caption and surrounding text. Fig. A.2

Same as Fig. A.1, but for models with a spherical DM component and a more narrow range of M•. The parameters are, from top to bottom: umin, pmin, qmin, f, M•, and ϒ.
In the text
Thumbnail: Fig. B.1 Refer to the following caption and surrounding text. Fig. B.1

Rmax vs. zmax diagram of the best-fit model. Each data point represents an orbit with non-zero weight, the colour represents the relative weight of the orbit, the solid line denotes 1 Re of the NSC, and the dashed line denotes the outer extent of the kinematic data. The panels show, from top to bottom, regular cold and warm tube orbits (λz>0.25), CR orbits (λz ≤-0.25), hot orbits, and all orbits.
In the text
Thumbnail: Fig. B.2 Refer to the following caption and surrounding text. Fig. B.2

Same as Fig. B.1, but zoomed into the region dominated by the NSC. Note the different colour scale to improve visibility of a range of orbit weights.
In the text
Thumbnail: Fig. C.1 Refer to the following caption and surrounding text. Fig. C.1

Radial velocity anisotropy ßr as function of intrinsic radius. The black line denotes the best-fit model without DM, the shaded region the 1σ uncertainty, the vertical solid line 1 Re of the NSC, the vertical dotted line the outer limit of the kinematic data.
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.