Distribution function-based modelling of discrete kinematic datasets, in application to the Milky Way nuclear star cluster

Abstract: We present a method for constructing dynamical models of stellar systems described by distribution functions and constrained by discrete-kinematic data. We implement various improvements compared to earlier applications of this approach, demonstrating with several examples that it can deliver meaningful constraints on the mass distribution even in situations when the density profile of tracers and the selection function of the kinematic catalogue are unknown. We then apply this method to the Milky Way nuclear star cluster, using kinematic data (line-of-sight velocities and proper motions) for a few thousand stars within 10 pc from the central black hole, accounting for the contributions of the nuclear stellar disc and the Galactic bar. We measure the mass of the black hole to be 4x10^6 Msun with a 10% uncertainty, which agrees with the more precise value obtained by the GRAVITY instrument. The inferred stellar mass profile depends on the choice of kinematic data, but the total mass within 10 pc is well constrained in all models to be (2.0-2.3)x10^7 Msun. We make our models publicly available as part of the Agama software framework for galactic dynamics.

• The paper introduces a likelihood-based, action-Df model that accurately recovers both SMBH and NSC mass profiles from discrete kinematic observations.
• It demonstrates methodological advancements in handling unknown selection functions, improving upon traditional Jeans and orbit-superposition techniques.
• The framework integrates line-of-sight velocities and proper motions of ~3200 stars using self-consistent, iterative modelling with the Agama software suite.

Distribution Function-Based Modelling of Discrete Kinematic Data in the Milky Way Nuclear Star Cluster

Introduction and Motivation

The analysis focuses on constructing equilibrium dynamical models for stellar systems using parametrized distribution functions (DFs), constrained solely by discrete kinematic measurements rather than full phase-space information—including spatial density profiles. The application is the Milky Way's nuclear star cluster (NSC), a unique system with abundant kinematic data within 10 pc of the central supermassive black hole (SMBH), Sgr A*. Existing dynamical mass measurements in the Galactic center, both for the SMBH and for the extended NSC mass profile, often rely on Jeans models or orbit-superposition methods, which present notable limitations—especially for discrete datasets and situations where the tracer spatial profile or the selection function is uncertain. This work introduces improvements to the DF-based modelling workflow, enabling robust recovery of the mass profile even under unknown selection functions or spatial incompleteness.

Key aims include:

1. Developing a practical, likelihood-based method for fitting equilibrium, action-based DFs to discrete kinematic data, accommodating uncharacterized selection effects.
2. Applying the method to the Milky Way NSC, testing SMBH mass recovery, and producing DF-based models using several thousand stars with line-of-sight velocities (V_LOS) and proper motions (PMs).

Distribution Function Formalism and Modelling Framework

The core approach utilizes self-consistent, action-based DFs for both the NSC and the nuclear stellar disc (NSD). Specifically, the NSC is modelled with a DoublePowerLaw DF parameterized by 11 shape, mass, and dynamical parameters, including mass ($M_{\mathrm{NSC}}$), SMBH mass, inner/outer power-law slopes ($\Gamma$, B), anisotropy parameters, and rotation controls. The NSD is included as a component with fixed parameters from prior large-scale fits.

A critical aspect is the self-consistent construction of the potential: the DF determines density, which via the Poisson equation yields the potential, thereby closing the loop. The Agama software framework is employed for the iterative solution and likelihood evaluation.

Handling Discrete and Incomplete Data: Selection Functions and Likelihood Construction

A major technical emphasis is on the likelihood function in the presence of unknown or non-uniform selection. While previous works often assume knowledge of the spatial density of tracers or treat the selection function as uniform, in reality observational samples have complex incompleteness—especially in the Galactic center with variable extinction. The likelihood function formulation distinguishes among cases:

• Complete uniform selection: Joint probability in full phase space (cf. Figure 1).
• Known but non-uniform spatial selection: Requires normalization over the selection-weighted DF.
• Unknown/Unformalizable selection: Only conditional velocities at known positions are used, effectively fitting the model VDFs to the observed kinematic sample.

  Figure 1: Schematic showing likelihood construction under uniform selection, spatially-restricted selection, and selection with unformalizable bias.

Toy models—single and two-population Maxwellian vertical discs; spherical double-power-law systems—quantitatively demonstrate the ability to constrain the potential and mass profile from kinematics alone, provided the DF exhibits sufficient complexity (non-Gaussian features or multiple kinematic populations).

Figure 2: Illustration of how the velocity distribution function (VDF) shape, varying with height or radius, informs constraints on the gravitational potential even in the absence of explicit positional density information.

Figure 3: Confidence intervals on mass density recovery from toy model spherical systems under various completeness and self-consistency assumptions.

Application to the Milky Way Nuclear Star Cluster

The observational sample (~3200 stars within 10 pc) combines V_LOS and proper motions from several sources. Notably, the model fitting explicitly marginalizes over unknown phase-space coordinates (e.g., line-of-sight position, missing velocity components) and integrates over uncertainties. The likelihood for each star is a conditional probability of the observed velocities given position, marginalized over unknowns and convolved with measurement errors.

Robustness checks include a mock catalogue with injected contamination and uncertainties, successfully recovering NSC and SMBH mass to within 1σ.

Figure 4: Parameter recovery for mock catalogues, showing posterior constraints on NSC and SMBH mass and 3D density profiles in ensemble fits.

Empirical Constraints: Mass, Structure, and SMBH Recovery

Multiple model series are fitted:

• Using only high-precision (AO-based) proper motions (central 2 pc)
• Extended with VIRAC2 proper motions (out to 10 pc)

Key results:

• SMBH mass consistently recovered: Within uncertainties, the model recovers $$M_{\mathrm{BH}} = (3.9 \pm 0.4)\text{--}(4.1 \pm 0.3) \times 10^6\,M_\odot$$, in agreement with GRAVITY's orbital-fitting measurements.
• Total mass within 10 pc robustly constrained: Both data-sets yield $M(<10\,\text{pc})$ in the $(2.0$–$2.3)\times 10^7\,M_\odot$ range, with differences in the NSC mass profile at larger radii due to the fixed NSD model.

The marginalized posterior for the radial density and mass profiles, as well as kinematic structure (dispersion, rotation) as a function of radius, are reported.

Figure 5: Mass profile, surface, and 3D density constraints, comparing different model series and literature reference fits.

Differences in the outer cluster slopes primarily reflect model assumptions (fixed NSD vs. mass assigned to NSC in the fits) rather than firm constraints. The inner NSC structure is inferred to be nearly isotropic in the central parsec, with axis ratio $q \sim 0.8$, aligning with photometric estimates.

Figure 6: Radial profiles of 1D velocity dispersions and anisotropy measures for various model series.

Treatment of Foreground and Bar Contamination

Model likelihoods incorporate explicitly a contamination model for foreground stars from the Galactic bar, using empirically determined velocity distribution functions and fixed normalizations—found to be essential for unbiased SMBH mass recovery.

Additional test fits demonstrate that improper treatment of bar contamination can systematically bias mass estimates downward by up to a third, underscoring the necessity for physically motivated mixture components.

Theoretical and Practical Implications

This study demonstrates that, under realistic conditions of incomplete or ill-characterized selection, action-based DF models can deliver accurate constraints on the gravitational potential and the SMBH mass. This robustness is particularly significant given that photometrically determined tracer densities are often unreliable due to extinction or confusion, especially in high-density, central Galactic environments.

The framework is flexible for generalization, including:

• Extension to non-axisymmetric geometries
• Inclusion of multiple chemical or age subpopulations (with "extended DFs" incorporating chemical abundance dependence)
• Joint inference on DF and selection function parameters

Potential future developments include simultaneous inference on both NSC and NSD DFs, inclusion of more extended outer kinematic datasets, and application to extragalactic NSCs or integrated light kinematic data.

Conclusion

The methodology described provides a statistically rigorous, computationally tractable approach for equilibrium dynamical modelling in Galactic nuclei and related systems, robust to uncertainties in selection and spatial distribution. The explicit demonstration that SMBH mass and extended mass profiles can be recovered from the kinematics of discrete tracers—even without direct use of the spatial density profile—clarifies the power and limitations of DF-based mass modelling. The implementation is released in the Agama suite, enabling broader adoption.

---

• "Distribution function-based modelling of discrete kinematic datasets, in application to the Milky Way nuclear star cluster" (2603.29502)