  One of the major goals of the field of Milky Way dynamics is to recover the
gravitational potential field. Mapping the potential would allow us to
determine the spatial distribution of matter - both baryonic and dark -
throughout the Galaxy. We present a novel method for determining the
gravitational field from a snapshot of the phase-space positions of stars,
based only on minimal physical assumptions. We first train a normalizing flow
on a sample of observed phase-space positions, obtaining a smooth,
differentiable approximation of the phase-space distribution function. Using
the collisionless Boltzmann equation, we then find the gravitational potential
- represented by a feed-forward neural network - that renders this distribution
function stationary. This method is far more flexible than previous parametric
methods, which fit narrow classes of analytic models to the data. This is a
promising approach to uncovering the density structure of the Milky Way, using
rich datasets of stellar kinematics that will soon become available.
