Determining the Milky Way gravitational potential without selection functions

TL;DR

The paper tackles the challenge of inferring the Milky Way's gravitational potential from stellar kinematics in the presence of spatial selection effects. It introduces a conditional Deep Potential variant that models the velocity distribution $p(\vec{v}\mid \vec{x})$ and the true tracer density $n(\vec{x})$ via the collisionless Boltzmann equation under stationarity, avoiding explicit selection-function modeling. By decomposing the distribution as $f(\vec{x},\vec{v})=n(\vec{x})\,p(\vec{v}\mid \vec{x})$ and employing continuous normalizing flows with Conditional Flow Matching, the method jointly learns the potential $Φ(\vec{x})$ and tracer densities. In a mock test with a complex 3D dust distribution, the conditional approach recovers the gravitational density with high fidelity and far fewer parameters than a selection-function-based variant, demonstrating scalability to Gaia-scale data and applicability to surveys with patchy sky coverage. This work thus enables robust MW potential inference in realistic observational conditions, reducing reliance on detailed modeling of intricate selection effects and broadening the method’s applicability to contemporary large-scale surveys.

Abstract

Selection effects, such as interstellar extinction and varying survey depth, complicate efforts to determine the gravitational potential - and thus the distribution of baryonic and dark matter - throughout the Milky Way galaxy using stellar kinematics. We present a new variant of the "Deep Potential" method of determining the gravitational potential from a snapshot of stellar positions and velocities that does not require any modeling of spatial selection functions. Instead of modeling the full six-dimensional phase-space distribution function $f\left(\vec{x},\vec{v}\right)$ of observed kinematic tracers, we model the conditional velocity distribution $p\left(\vec{v}\mid\vec{x}\right)$, which is unaffected by a purely spatial selection function. We simultaneously learn the gravitational potential $Φ\left(\vec{x}\right)$ and the underlying spatial density of the entire tracer population $n\left(\vec{x}\right)$ - including unobserved stars - using the collisionless Boltzmann equation under the stationarity assumption. The advantage of this method is that unlike the spatial selection function, all of the quantities we model, $p\left(\vec{v}\mid\vec{x}\right)$, $Φ\left(\vec{x}\right)$, and $n\left(\vec{x}\right)$, typically vary smoothly in both position and velocity. We demonstrate that this "conditional" Deep Potential method is able to accurately recover the gravitational potential in a mock dataset with a complex three-dimensional dust distribution that imprints fine angular structure on the selection function. Because we do not need to model the spatial selection function, our new method can effectively scale to large, complex datasets while using relatively few parameters, and is thus well suited to Gaia data.

Figures (5)

• Figure 1: Two-dimensional slices (the $z = 0$ plane) through the three-dimensional dust extinction density (left) and integrated extinction (as viewed from the origin; middle left), the resulting selection function (middle right), and the observed spatial distribution of stars (right). The observer is assumed to be at the origin. Dust extinction imprints fine angular structure on the selection function and the observed spatial distribution of stars, which can be difficult to capture with neural networks.
• Figure 2: The gravitational densities ($\rho = \nabla^2 \Phi / 4 \pi G$) and their residuals (vs. the true densities), as recovered by the three variants of Deep Potential. Top row:$\ln \rho$ as a function of radius $r$, with the solid curves showing the median (in spherical shells) and the shaded envelopes enclosing the 16th to 84th percentiles. The dotted curves show the true density. Second row: Residuals $\Delta \ln \rho = \ln \rho_{\mathrm{model}} - \ln \rho_{\mathrm{true}}$ as a function of $r$. Third row: Residuals as a function of the true selection function $S\left(\vec{x}\right)$, illustrating the deterioration in the recovered density as the fraction of stars observed decreases. Bottom row: Residuals in the selection function $S\left(\vec{x}\right)$ (as recovered by the "selection function" method) and true spatial density of tracers $n\left(\vec{x}\right)$ (as recovered by the "conditional" method), each as a function of the true $S\left(\vec{x}\right)$. In the top three rows, the vertical axes are the same for the "selection function" and "conditional" methods.
• Figure 3: Recovery of the selection function and true spatial density of tracers by the "selection function" and "conditional" methods, respectively, in a two-dimensional slice through the $x-y$ plane (at $z=0$). Top row: The recovered selection function (left), true selection function (middle), and residuals (right) for the "selection function" method. Bottom row: The recovered true spatial density of tracers (left), true spatial density of tracers (middle), and residuals (right) for the "conditional" method. Dust imprints fine angular structure on the selection function, necessitating much larger models ($\sim$1.4 million parameters in this case) than are needed to capture the much smoother true spatial density of tracers ($\sim$47 thousand parameters). As seen in the top-right panel, the selection-function residuals appear as long radial spurs, illustrating the difficulty of recovering the angular structures in $S\left(\vec{x}\right)$.
• Figure 4: Sky plots of the "selection function" method's prediction of the selection function $S(\vec{x})$, the true selection function, and their fractional residuals on the surface of a sphere of radius 4kpc, centered on the mock observer. Areas with $S_\mathrm{true} < 0.01$ are masked. The differences can be attributed to shot noise in areas of low completeness, and to the difficulty of capturing fine angular structure with a neural network. Nevertheless, with a mock selection function with structure down to $\sim$0.25 kpc scales, our selection-function implementation performs well, albeit at the cost of large model sizes (here, $\sim$1.4 million parameters).
• Figure 5: Top panel: Residuals between our normalizing flow model of the distribution of observed kinematic tracers and the ground truth: $\Delta \ln n_{\mathrm{obs}}\left(\vec{x}\right) \equiv \ln n_{\mathrm{obs}}\left(\vec{x}\right) - \ln \left[n\left(\vec{x}\right) S\left(\vec{x}\right)\right]$, where $n\left(\vec{x}\right)$ and $S\left(\vec{x}\right)$ are the true spatial density of tracers (including those that are not observed) and selection function, respectively. Bottom panel: Residuals between our conditional normalizing flow model of the velocity distribution and the ground truth. The conditional distribution $p\left(\vec{v}\mid\vec{x}\right)$ is recovered to within $\sim$2%, even for $S\left(\vec{x}\right)$ as low as $\sim$0.01. Due to the complex nature of the spatial selection function, the fractional errors in the recovered $n_{\mathrm{obs}}\left(\vec{x}\right)$ are larger, and the median residual diverges substantially from zero at low $S\left(\vec{x}\right)$.