  A tracer sample in a gravitational potential, starting from a generic initial
condition, phase-mixes towards a stationary state. This evolution is
accompanied by an entropy increase, and the final state is characterized by a
distribution function (DF) that depends only on integrals of motion (Jeans
theorem). In terms of angle-action variables, the final state is uniform in
angles (high entropy) and maximally clustered in actions (low entropy). We
present a method exploring this fact to constrain a gravitational potential
using a sample that is stationary in it. We estimate the entropy in the action
space of trial potentials and recover the true potential by minimizing this
entropy. This method avoids assuming a known DF, and may be applicable to other
sets of integrals. We provide expressions for the entropy of DFs depending on
energy, $f(E)$, energy and angular momentum, $f(E,L)$, or three actions,
$f(\vec{J})$, and investigate the bias and fluctuations in their estimates. We
show that the method correctly recovers the potential parameters for spherical
and axisymmetric potentials. We also present a methodology to characterize the
posterior probability distribution of the parameters with an Approximate
Bayesian Computation, and indicate a pathway for application to observational
data. Using $N=10^4$ tracers with $20\%$-uncertainties in the 6D coordinates,
we recover the flattening parameter $q$ of an axisymmetric potential with
$\sigma_q/q\sim 10\%$.
