  The advent of datasets of stars in the Milky Way with six-dimensional
phase-space information makes it possible to construct empirically the
distribution function (DF). Here, we show that the accelerations can be
uniquely determined from the DF using the collisionless Boltzmann equation,
providing the Hessian determinant of the DF with respect to the velocities is
non-vanishing. We illustrate this procedure and requirement with some analytic
examples. Methods to extract the potential from datasets of discrete positions
and velocities of stars are then discussed. Following Green & Ting
(arXiv:2011.04673), we advocate the use of normalizing flows on a sample of
observed phase-space positions to obtain a differentiable approximation of the
DF. To then derive gravitational accelerations, we outline a semi-analytic
method involving direct solutions of the over-constrained linear equations
provided by the collisionless Boltzmann equation. Testing our algorithm on mock
datasets derived from isotropic and anisotropic Hernquist models, we obtain
excellent accuracies even with added noise. Our method represents a new,
flexible and robust means of extracting the underlying gravitational
accelerations from snapshots of six-dimensional stellar kinematics of an
equilibrium system.
