  We present a method to calculate gravitational potential gradients within
regions containing few tens of thousands stars with known phase space
coordinates. The central idea of the method is to calculate orbital arcs for
each star within a given region for a certain parametrised potential
(gravitational acceleration) and to assume that position of each star on its
orbital arc is a random variable with a uniform probability density in time.
Thereafter, by combining individual probability densities of stars it is
possible to calculate the overall probability density distribution and
likelihood for a given region as a function of gravitational acceleration
parameters. The likelihood has a maximum if the calculated probability
distribution and the observed distribution of stars in phase space are
consistent. This allows us to constrain gravitational accelerations and
potential gradient values. The method assumes that phases of stars are mixed
within the regions where stellar orbits are calculated. We tested the method
for 12 small rectangular regions within simulated disc galaxy from Gaia Wiki.
Tests show that even with a rather simple acceleration form the calculated
accelerations in galactic plane coincide with their true values from simulation
about 5 per cent, misalignment between the calculated and true acceleration
vector directions is less than 1 degree (median values). The model can be used
with the Milky Way Gaia complete solution data.
