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.

Figures (13)

• Figure 1: Illustration of the fitting process in three scenarios, depending on the distribution of stars with kinematic measurements (black dots) and without (red). In panel A, the observed sample traces the entire population, i.e., the selection probability is uniform across the entire system. In panel B, only a https://galileo-unbound.blog/2019/06/16/vladimir-arnolds-cat-map/ is observed (the spatial selection function $S(\boldsymbol x)$ is zero outside the region, and needs not be uniform within it, but at least is known in advance). In panel C, stars are selected for kinematic observations based on some criterion that is difficult or impossible to formalise, making the selection function uncomputable. $\ln\mathcal{L} = \sum_{i=1}^{N_\text{stars}} \ln f(\boldsymbol x_i, \boldsymbol v_i)$.$\ln\mathcal{L} = \sum_{i=1}^{N_\text{stars}} \ln \frac{f(\boldsymbol x_i, \boldsymbol v_i)}{\mathscr N}$,$\ln\mathcal{L} = \sum_{i=1}^{N_\text{stars}} \ln \mathfrak{f}(\boldsymbol v_i \:|\: \boldsymbol x_i)$,$\mathscr N \equiv \iiint \text{d}^3 \boldsymbol x \iiint \text{d}^3\boldsymbol v\; f(\boldsymbol x, \boldsymbol v)\; S(\boldsymbol x)$.$\mathfrak{f}(\boldsymbol v \:|\: \boldsymbol x) \equiv \frac{f(\boldsymbol x, \boldsymbol v)}{\iiint \text{d}^3 \boldsymbol v\; f(\boldsymbol x, \boldsymbol v)} = \frac{f(\boldsymbol x, \boldsymbol v)}{\rho(\boldsymbol x)}$. $\blacktriangleleft$$\blacktriangleleft$Note that in the first two cases, $f(\boldsymbol x, \boldsymbol v)$ is the joint probability of observing a star with given phase-space coordinates, whereas in the last case, $\mathfrak{f}(\boldsymbol v \:|\: \boldsymbol x)$ is the conditional probability (i.e., the velocity distribution function at a given position). $\mathscr{N}$ is the normalisation factor, i.e., the integral of $f$ over the entire phase space, weighted by the selection function $S$ of the observed sample, which is usually purely spatial (independent of velocity). In the case of a complete or at least a uniform coverage ($S(\boldsymbol x)$=const), the second case is equivalent to the first one, but the third one is not: it entirely ignores the spatial distribution of stars, since one cannot assume that they faithfully trace the entire population. However, if this is actually the case, one can add another term in the likelihood function corresponding to the spatial probability, i.e., $\ln\mathcal{L} = \sum_{i=1}^{N_\text{stars}} \ln \rho(\boldsymbol x_i)$, and then the result becomes equivalent to the first case.
• Figure 2: Illustration of constraining the potential using the 1d toy model (Section \ref{['sec:toy_1d']}) of a two-component tracer population with different velocity dispersions. Left and centre column demonstrate the basic idea: the density profiles $\rho_\star(z)$ of both components depend on the potential $\Phi(z)$ through Equation \ref{['eq:df1d_rho']}, but this dependence is stronger for the colder component (blue) than for the warmer one (red). Centre column shows the velocity distribution functions (VDFs) at three values of $z$ (marked by vertical dashed lines in the left column), separately for the cold and warm components (blue and red Gaussians), as well as their sum; solid lines correspond to the true potential, and dashed -- to a different (wrong) potential. The amplitudes of the two Gaussians are the density values at the given $z$, so they change depending on the potential. At $z=0$ (top row) the potential is the same, and so are the VDFs (up to a small difference in the overall amplitude). At an intermediate value $z=1.5$ (middle row), the density of the cold component is much lower in the wrong potential than in the true one, so the shape of the VDF is sensitive to the choice of $\Phi$. At an even higher value $z=3$ (bottom row), the cold component makes negligible contribution in both potentials (the dashed blue line is not even visible on this plot), so the VDF is again nearly the same in both cases up to a difference in the overall amplitude. This suggests that one can constrain the potential from the shape of the VDF, as long as there is a mixture of populations with different temperatures. Right column shows the results of fitting the mass profile specified by Equation \ref{['eq:df1d_trialdensity']} to the sample of $10^4$ tracers, using both their positions and velocities in the likelihood function (Equation \ref{['eq:df1d_jointposvel']}, cyan contours) or only the velocities conditioned on the position (Equation \ref{['eq:df1d_onlyvel']}, yellow contours). The shaded regions correspond to 16/84% confidence intervals of the inferred potential (top) and density profiles of both components (bottom), while the dashed lines show the true profiles (same as solid lines in the left column). It demonstrates that one can get rather tight constraints in the region where the two kinematically distinct components overlap, but the uncertainties explode at $z\ge 2$ where only one (warmer) component is dominant.
• Figure 3: Illustration of the dynamical modelling precision in different scenarios using the toy spherical model (Section \ref{['sec:toy_spherical']}). Shown are constraints on the mass density (left and middle columns) and tracer density (right column), in the cases when the spatial selection function is known and uniform (scenario A in Figure \ref{['fig:selection_function']}; green and cyan) or unknown (scenario C, which uses only the kinematic information but not the spatial distribution of tracers; magenta and red curves). In the self-consistent case (left column), tracer and mass densities are identical, while in the opposite case, they are shown separately (middle and right columns). The top row shows the true profile as a black dashed line, and the 16/84% confidence intervals on $\rho(r)\,r^3$ as shaded regions, and bottom row shows the square root of the ratio of upper and lower bounds of these intervals (qualitatively, the relative uncertainty, with 1 corresponding to a precise constraint and $\gg1$ to a poor constraint). Naturally, the constraints are tighter at intermediate radii where the majority of tracers reside ($\rho r^3$ shows the mass of stars per logarithmic interval in radius). The radial range roughly encompasses all $10^3$ kinematic tracers used in the fit.
• Figure 4: Left half: Illustration of the full 3d velocity distribution shown as its 2d projections (bottom left corner) and 1d projections (diagonal panels) at a location $X=3$ pc, $Y=2$ pc from one of the models in the MCMC chain. Details to note are the obviously non-Gaussian shape of marginalised 1d distribution and non-elliptical contours in the 2d distributions, a slight tilt of the $V_\text{hor}$--$V_\text{ver}$ distribution (middle row, left panel), and more prominently, a banana-shaped peak region in the $V_\text{hor}$--$V_\textrm{LOS}$ distribution (bottom row, left panel). While the tilt is imprinted in the non-diagonal elements of the velocity dispersion tensor, the banana-shaped peak does not leave a specific signature in it. For comparison, the top right corner shows the approximation of the full velocity distribution (mirrored w.r.t. the diagonal from the bottom right corner, i.e., rotated by 90 degrees) as a product of 1d distributions, which lacks these asymmetries. The full 3d velocity distribution is used as the main method of likelihood computation (Section \ref{['sec:projected_df']}), while the product of marginalised 1d distributions is an alternative method used in past studies and described in Section \ref{['sec:vdf']}. Right half: Comparison of total likelihood values of the entire sample computed with these two methods for a range of models from the MCMC chain. $\ln\mathcal{L}_\text{3d VDF}$ uses the full 3d velocity distribution for each star (Equation \ref{['eq:conditional_projected_df']}), $\ln\mathcal{L}_\mathrm{3\!\times\!1d\, VDF}$ uses the product of marginalised 1d distributions (Equation \ref{['eq:vdf1d']}). Both are shifted by the same constant (the maximum value of the first quantity). The values of log-likelihoods for individual stars computed using both methods are fairly close, with a mean difference of only 0.02 and a scatter of 0.2; however, this nonzero mean offset accumulates over the entire sample to produce an overall shift in log-likelihoods of order a few dozen. There is still a good correlation between the total log-likelihoods computed using both methods, illustrated by the red dotted line with a unit slope, but the scatter of order a few in $\ln\mathcal{L}_\mathrm{3\times1d}-\ln\mathcal{L}_\mathrm{3d}$ may affect the relative odds of models in different parts of the parameter space.
• Figure 5: Test of modelling machinery on the mock data generated from the fiducial model with parameters given in Table \ref{['tab:fiducial_model']}. The left panels show the posterior distributions of the NSC total mass and the SMBH mass. The middle panel shows the 3d density profiles of both NSC and NSD measured along the major axis; the shaded region corresponds to 16th/84th percentiles (note that the NSD DF is kept fixed throughout the fit, and thus its density profile varies very little), and the short-dashed line shows the true profile of the fiducial model. The right panel shows the axis ratio of the NSC in the same way.