  Next-generation planetary tracking methods, such as interplanetary laser
ranging (ILR) and same-beam interferometry (SBI) promise an orders-of-magnitude
increase in the accuracy of measurements of solar system dynamics. This
requires a reconsideration of modelling strategies for the translational and
rotational dynamics of natural bodies, to ensure that model errors are well
below the measurement uncertainties.
  The influence of the gravitational interaction of the full mass distributions
of celestial bodies, the so-called figure-figure effects, will need to be
included for selected future missions. The mathematical formulation of this
problem to arbitrary degree is often provided in an elegant and compact manner
that is not trivially relatable to the formulation used in space geodesy and
ephemeris generation. This complicates the robust implementation of such a
model in operational software packages. We formulate the problem in a manner
that is directly compatible with the implementation used in typical dynamical
modelling codes: in terms of spherical harmonic coefficients and Legendre
polynomials. An analytical formulation for the associated variational equations
for both translational and rotational motion is derived.
  We apply our methodology to both Phobos and the KW4 binary asteroid system,
to analyze the influence of figure-figure effects during estimation from
next-generation tracking data. For the case of Phobos, omitting these effects
during estimation results in relative errors of $0.42\%$ and $0.065\%$ for the
$\bar{C}_{20}$ and $\bar{C}_{22}$ spherical harmonic gravity field
coefficients, respectively. These values are below current uncertainties, but
orders of magnitude larger than those obtained from past simulations for
accurate tracking of a future Phobos lander.
