Propagation and Estimation of the Dynamical Behaviour of Gravitationally Interacting Rigid Bodies

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.


Introduction
For the robust analysis of tracking data from planetary missions, the dynamics of solar system bodies under investigation should ideally be modelled to well below the observational accuracy and precision. Several exceptionally accurate tracking-data types are emerging for planetary missions, such as multi-wavelength radiometric range and Doppler measurements (Dehant et al. 2017), same-beam interferometry (SBI) (Kikuchi et al. 2009;Gregnanin et al. 2012), and interplanetary laser ranging (ILR) (Degnan 2002;Turyshev et al. 2010;Dirkx 2015). For the analysis of these data, dynamical models for natural bodies need to be developed and implemented to beyond the current state-of-the-art of typical state propagation and estimation software.
Examples of such software tools are GEODYN (Genova et al. 2016), GINS (Marty et al. 2009), GMAT (Hughes et al. 2014), NOE (Lainey et al. 2004), OREKIT (Maisonobe and Pommier-Maurussane 2010) and Tudat (which we use in this manuscript, see Appendix C). We stress that the full functionality of several of these codes (GMAT, OREKIT and Tudat being the exceptions) cannot be transparently determined, as up-to-date source code and documentation is not openly available for them. In this article we discuss, and present models to mitigate, one of the common challenges that these tools face for the analysis of future planetary tracking data.
The specific physical effects that must be incorporated for future missions depend strongly on the object under consideration, and the available tracking data types. For both SBI and ILR, there is a need for sub-mm accurate dynamical models over the time span of the mission. For Doppler data, variations in range at the level of 1-10 µm/s need to be accounted for. To meet the requirements that result from these tracking accuracies, various dynamical models may need to be improved, depending on the situation under consideration. Examples of such models include: a fully consistent tidal-rotational-translational dynamics model, realistic models for frequency-dependent tidal dissipation, detailed nonconservative force models for small bodies, and figure-figure gravitational interactions between massive bodies. The development of a model of the latter, for the purpose of state and parameter estimation, is the topic of the present paper.
Modelling barycentric motion to the mm-level is presently limited by the knowledge of the properties of small solar system bodies, in particular main-belt asteroids (e.g., Kuchynka et al. 2010). However, the relative motion of solar system bodies in close proximity (e.g., planetary satellite system or multiple asteroid system) is dominated by the gravity fields of these bodies themselves. Their relative motion is only weakly influenced by the gravitational fields of the other solar system bodies. As a result, the uncertainties in the local dynamics of such systems stem largely from uncertainties in, and mismodelling of the effects of, the gravitational interactions in the system itself. Measurements of the local dynamics can be instrumental in improving the estimates of the properties of the bodies in the system (e.g., Lainey et al. 2007;Folkner et al. 2014;Dirkx et al. 2016), provided that the dynamical model can be set up and parameterized to sufficient accuracy.
Currently, the dynamical models of planetary/asteroid systems that are used in typical state propagation and estimation software cannot robustly capture the motion to the measurement accuracies of ILR and SBI (e.g., Dirkx 2015;Dehant et al. 2017), which would prevent the data from being optimally exploited. Among others, characterizing such systems' dynamics will require a new level of detail for the models used to describe the gravitational interaction between extended bodies. Specifically, the coupling between higher-order terms in the gravity field expansions (so-called figure-figure effects; Bois et al. 1992) will need to be included when propagating and estimating the translational and rotational dynamics of such bodies. Although such models are incorporated in LLR data analysis software (although not necessarily for arbitrary degree and order), the underlying models are not clearly described in literature, nor are these software frameworks openly available. The resulting mathematical problem is also termed the full two-body problem (F2BP). The level to which these effects need to be included will depend strongly on the system under consid-eration. However, the a priori assumption that such figurefigure terms can be neglected (e.g., Lainey et al. 2007) will no longer be a given for many cases with high-accuracy, next-generation tracking systems. At the very least, an evaluation of the magnitude of the influence of these terms should be made before performing the actual data processing.
The influence of low-order figure-figure terms on the translational and/or rotational dynamics of solar system bodies has been analyzed for a variety of cases, such as the Moon (Bois et al. 1992;Müller et al. 2014), Phobos (Borderies and Yoder 1990;Rambaux et al. 2012), and binary asteroids (Fahnestock and Scheeres 2008;. Their analyses show that including figure-figure interactions is important for accurate dynamical modelling of selected systems of interacting bodies. For multiple asteroid systems, the higher-order gravitational interactions are especially strong, as a result of their highly irregular shapes and close orbits. As discussed by Batygin and Morbidelli (2015), understanding the spin interaction of these bodies is crucial in building a complete picture of the dynamical evolution of the solar system. A body's rotational state is a key parameter in determining the influence of dissipative effects, which in turn play an important role in a body's long-term evolution.
A general formulation of mutual gravitational interaction potential of two extended bodies, which can be used to fully model such effects, was developed by Sidlichovsky (1978) and Borderies (1978). Subsequently, Maciejewski (1995) used these results to set up general translational and rotational equations of motion, later formally derived by Lee et al. (2007), and extended to N bodies by Jiang et al. (2016), including the static electric and magnetic potential. This method is described and applied further by Mathis and Le Poncin-Lafitte (2009), and Compère and Lemaître (2014), using symmetric trace-free (STF) tensors (Hartmann et al. 1994) and mass multipole moments. Recently, an efficient representation of this problem was introduced by Boué (2017) by applying angular momentum theory. An equivalent formulation of the problem, in terms of inertia integrals instead of mass multipole moments, was developed by Paul (1988), Tricarico (2008), , with a highly efficient implementation presented by Hou (2018). In an alternative approach, a formulation of the mutual interaction of homogeneous bodies is derived by Werner and Scheeres (2005), Fahnestock and Scheeres (2006), Hirabayashi and Scheeres (2013) based on polyhedron shape models, which is highly valuable for the simulation of small bodies, such as binary asteroids (Fahnestock and Scheeres 2008).
Explicit expansions of the mutual two-body interaction to low order have been derived by Giacaglia and Jefferys (1971), Schutz (1981), Ashenberg (2007), Boué and Laskar (2009), and Dobrovolskis and Korycansky (2013), using a variety of approaches. For the analysis of future trackingdata types, the inclusion of higher-order interactions effects will be relevant, especially for highly non-spherical bodies in close orbits, such as binary asteroids . The need for figure-figure interactions in lunar rotational dynamics when analyzing LLR data is well known (Eckhardt 1981). Recent re-analysis by Hofmann (2017) has shown the need to use the figure-figure interactions up to degree 3 in both the rotational and translational dynamics of the Earth-Moon system, for the analysis of modern LLR data. With the exception of LLR, figurefigure interactions have not been applied in tracking data analysis, and full algorithms to do so are not available in literature, nor is software to perform these analyses. Errors in dynamical modelling during data analysis can lead to biased estimates, and a true estimation error that is many times larger than the formal estimation error.
The formulation of equations of motion in the F2BP does not trivially lend itself to the direct implementation in typical state propagation and estimation software tools. In such codes, the gravitational potential is described by the (normalized) spherical harmonic coefficients and Legendre polynomials, as opposed to the multipole moments/STF tensors and inertia integrals used in the F2BP. This gap between theoretical description and practical implementation must be closed before the figure-figure effects can be routinely included up to arbitrary degree in tracking-data analysis for (future) missions. Moreover, transparently and consistently including the figure-figure effects in a general manner in orbit determination and ephemeris generation algorithms has not yet been explored in detail.
In this article, the main goal is to derive a direct link between the theoretical model for the F2BP and the implementation of the one-body potential, for both the propagation and estimation of the translational and rotational dynamics of the system. This will bridge the existing gap between theory and implementation in the context of spacecraft tracking and planetary geodesy. We start our development in Sect. 2 from the formulation of Mathis andLe Poncin-Lafitte (2009), Compère andLemaître (2014) and Boué (2017), and derive a direct and explicit link with typical one-body implementations. In Sect. 3, we present the equations of motion and derive the associated variational equations, allowing the models to be used in orbit determination and ephemeris generation. A consistent formulation of the variational equations is crucial for the extraction of physical signatures from the (coupled) translational and rotational dynamics, from tracking data. In Sect. 4, we illustrate the impact of our method, by analyzing how estimation errors of the gravity field of Phobos, and the two bodies in the KW4 binary asteroid system, are affected if figure-figure interactions are omitted during the estimation. Section 5 summarizes the main results and findings.
Our focus is on the development of an explicit link between the F2BP formulation and the formulations used in orbit determination/ephemeris generation, while ensuring a computational efficiency not prohibitive from a practical point of view. Our goal is not to improve the current state of the art in terms of computational performance (Boué 2017;Hou 2018).

Gravitational potential
We start by reviewing the formulation of the one-body potential in Sect. 2.1, followed by a discussion of the full twobody potential in Sect. 2.2. We discuss the transformation of the spherical harmonic coefficients between two reference frames in Sect. 2.3. Finally, we provide explicit expressions relating the computation of terms from the one-body potential to that of the full two-body potential in Sect. 2.4.

Single-body potential and notation
Applications in planetary geodesy typically represent the gravitational field of a single extended body by means of a spherical harmonic expansion of its gravitational potential (e.g., Montenbruck and Gill 2000;Lainey et al. 2004): where r denotes the position at which the potential is evaluated and s denotes the position inside the body of the mass element dM. The spherical coordinates (radius, longitude, latitude) in a body-fixed frame are denoted by r, ϑ and ϕ. The reference radius of the body is denoted by R and μ is the body's gravitational parameter. P lm and Y lm denote the unnormalized Legendre polynomials and spherical harmonic basis functions, respectively (both at degree l and order m). The term U lm is the full contribution from a single degree l and order m to the total potential. M lm represents the unnormalized mass multipole moments (typically used in the F2BP), and C lm and S lm are the spherical harmonic coefficients (typically used in spacecraft tracking analysis). The terms M lm are related to C lm and S lm as: where * indicates the complex conjugate. The spherical harmonic basis functions in Eq.
(3) can be expressed as: Often the mass multipoles and basis functions are represented in a normalized manner. In this manuscript, we will apply two normalizations: 4π -normalized (for which quantities are represented with an overbar), and Schmidt semi-normalized (for which quantities are represented with a tilde). The Schmidt semi-normalized formulation is obtained from: which are used in the formulations of Compère and Lemaître (2014) and Boué (2017). In planetary geodesy, the 4πnormalized coefficients are typically used, for which: For the 4π -normalized coefficients (which we shall simply refer to as 'normalized' from now on): =C lm − iS lm (13) where we have introduced theC lm ,S lm notation (which we stress are distinct from C lm and S lm ) to avoid awkward expressions in later derivations. None of the final quantities that are needed in the computation are complex (Sect. 3). However, the complex number notation is more concise, so we retain it in some sections to keep the derivation and analytical formulation tractable.

Two-body interaction potential
For the interaction between extended bodies, we use the mutual force potential introduced by Borderies (1978). It is ob-tained from the following: where d 12 denotes the distance between the mass elements dM 1 and dM 2 and r 1 and r 2 denote the inertial positions of the centers of mass of bodies 1 and 2, respectively. B 1 and B 2 denote full volume of bodies 1 and 2, respectively. The rotation from a frame A to a frame B is denoted as R B /A , while F i denotes the frame fixed to body i and I denotes a given inertial frame (such as J2000).
The double integral in Eq. (16) can be expanded in terms of the mass multipoles and spherical harmonics (Borderies 1978). We use a slightly modified 1 form of the notation used by Mathis and Le Poncin-Lafitte (2009), Compère and Lemaître (2014): where the distance between the centers of mass r 21 = r 2 −r 1 is written as r. The F 1 superscript denotes that a vector is represented in the body-fixed frame of body 1. The angles ϑ and ϕ denote the latitude and longitude of body 2, expressed in frame F 1 . The termγ l 1 ,m 1 l 2 ,m 2 is a scaling term (Mathis and Le Poncin-Lafitte 2009), which can be written as: from which it follows thatγ l 1 ,m 1 l 2 ,m 2 = 1, if l 1 = 0 or l 2 = 0. As discussed in Sect. 1, our goal is to find an explicit expression relating the implementation of Eq. (17) to that of Eq. (15), which uses 4π -normalized mass multipoles. Using Eqs. (6)-(14), the terms V l 1 ,m 1 l 2 ,m 2 can be rewritten explicitly as follows: 1,2;F 1 l 1 ,m 1 ,l 2 ,m 2 cos |m 1 + m 2 |ϑ 1 we use a ∼ to denote the semi-normalized parameters used by Compère and Lemaître (2014), Boué (2017) to distinguish the terms from our (un)normalized formulations. + i s m 1 +m 2 sin |m 1 + m 2 |ϑ P l 1 +l 2 ,|m 1 +m 2 | (sin ϕ) r (20) Here, we have introduced effective two-body multipole mo-mentsM 1,2;F 1 l 1 ,m 1 ,l 2 ,m 2 , defined by: The moments are expressed in the frame of body 1, and are therefore dependent on R F 1/F 2 if l 2 > 0.
The real part of the formulation for V l 1 ,m 1 l 2 ,m 2 in Eq. (20) is similar to a single term U lm of the one-body potential in Eq. (2). Consequently, this formulation lends itself to the implementation in typical state propagation and estimation software (see Sect. 1) by the correct change of variables, as we will discuss in detail in Sect. 2.4. In later sections, the following decomposition for V l 1 ,m 1 l 2 ,m 2 will ease some derivations: l 1 +l 2 ,m 1 +m 2 (ϑ, ϕ) r l 1 +l 2 +1 (26) which explicitly separates the dependency on R F 1/F 2 and r F 1 . We assume that the mass multipolesM i,F i are timeindependent (in their local frames F i ). In principle, the inclusion of tidal effects (Mathis and Le Poncin-Lafitte 2009) does not fundamentally change the formulation of the mutual force potential. However, it does make the M i,F i lm terms dependent on the relative positions and orientations of the bodies, substantially complicating the analytical formulation of the derivatives of these terms w.r.t. position and orientation (Sect. 3). Therefore, we limit ourselves to static gravity fields in this article, focussing on the relation between the one-body and two-body potential.

Transformation of the gravity field coefficients
The main complication of using the mutual force potential in Eq. (17) lies in the orientation dependency ofM 2,F 1 lm . Determining these values requires a transformation of multipole momentsM 2 lm from F 2 (in which they are typically defined) to F 1 . A transformation from the semi-normalized multipolesM 2,F 2 lm toM 2,F 1 lm is given by Boué (2017), based on the methods from Wigner and Griffin (1959), discussed in detail by Varshalovich et al. (1988): where D l mk represents the Wigner D-matrix of degree l (with −l ≤ m, k ≤ l). Expressions for D l mk can be found in literature in terms of Euler angles and Cayley-Klein parameters (among others). Here, we choose to express it in terms of the non-singular Cayley-Klein parameters, defined by two complex parameters a and b, which are closely related to the unit quaternion more typically used in celestial mechanics (Appendix A.2). We denote the vector containing the four elements of a and b as c: We follow the same computational scheme as Boué (2017) to determine the Wigner D matrices, which is a recursive formulation based on Gimbutas and Greengard (2009). Analytical formulations for D 0 mk and D 1 mk are given in terms of a and b by Varshalovich et al. (1988) and Boué (2017) as: The following recursive formulation is then applied for m ≥ 0 and l ≥ 2: where D l mk = 0 if |m| > l or |k| > l. For m < 0: The transformation in Eq. (28) can be rewritten in terms of fully normalized multipole momentsM by using Eqs. (8)-(11) to obtain: Now, by virtue of Eq. (12) the transformation only needs to be performed for positive m, and we can rewrite Eq. (36) as: which allows for a direct and transparent relation to spherical harmonic coefficients to be made (see Sect. 2.4).

Explicit formulation in terms of one-body potential
Here, we explicitly provide equations linking the F2BP to the governing equations of the one-body potential in terms of spherical harmonics. From Eqs. (12) and (38), the explicit equations for the transformed spherical harmonic coefficients become: A single term V l 1 ,m 1 l 2 ,m 2 of the mutual force potential coefficients can be computed using the exact same routines as those for computing a single term of the one body potential U lm by a correct substitution of variables. Comparing Eq. (15) for U lm with the real part of Eq. (20) for V l 1 ,m 1 l 2 ,m 2 , we obtain, with Eq. (25): with theC,S terms defined by Eq. (13). In the above, we have introduced the (C,S) l 1,2 ;m 1,2 notation to denote the effective one-body spherical harmonic coefficients that are to be used to evaluate a single term V l 1 ,m 1 l 2 ,m 2 . By substituting Eq. (12) into Eqs. (45) and (46) we obtain (omitting the F 1 superscripts): where the s m and σ m functions are defined in Eq. (22). These equations provide the direct formulation of the effective two-body spherical harmonic coefficients in terms of the respective (transformed) one-body spherical harmonic coefficients: × P lm (sin ϕ)(C l 1,2 ;m 1,2 cos mϑ +S l 1,2 ;m 1,2 sin mϑ) with m = m 1 + m 2 and l = l 1 + l 2 .

Degree-two interactions-circular equatorial orbit
To gain preliminary insight into the influence of figurefigure interactions, we perform a simplified analytical anal-ysis of their effects. As a test case, we take two bodies in mutual circular equatorial orbits, with both bodies' rotations tidally locked to their orbit. In this configuration, the tidal bulges always lie along the same line, with a constant distance between the two bodies, and the frames F 1 and F 2 are equal, so that (C,S) 2,F 1 l,m = (C,S) 2,F 2 l,m . This represents a highly simplified model for, for instance, a binary asteroid system.
Consequently, the figure-figure interactions of degree 2 present themselves in a similar manner as one-body interactions with (l, m) = (4, 0), (4, 2) and (4, 4). As a result, the process of determining the relevance of figure-figure interactions is analogous, in our simplified situation, to determining relevance of one-body interactions of degree 4 terms, under the suitable change of variables forC 4,0 ,C 4,2 andC 4,4 . For more realistic cases (non-zero eccentricity and/or inclination), the effective spherical harmonic coefficients will be time-dependent. However, for low-eccentricity/inclination situations, the deviations will be relatively low, as long as bodies remain tidally locked, and consequently the frames F 1 and F 2 remain close.
The above analysis is also applicable in the case where the rotation of body 1 is not tidally locked to body 2, so long as we setC 1 22 to zero (no tidal bulge). Such a situation would be representative of a planetary satellite orbiting its host planet.

Dynamical equations
Here, we set up our equations of translational and rotational motion in Sect. 3.1, following the approach of Maciejewski (1995) and Compère and Lemaître (2014). We use the algorithm by Boué (2017) for the calculation of the torques, as it provides an efficient, elegant and non-singular implementation. Subsequently, we derive an analytical formulation for the variational equations in Sect. 3.2. We define the governing equations of the dynamics of the bodies in an inertial frame, as opposed to the mutual position and orientation of the bodies that are used by Maciejewski (1995) and Compère and Lemaître (2014). Such an approach is more in line with typical implementation of few-body codes (e.g., solar system simulators; see Sect. 1). Also, it enables an easier implementation for the interaction of N extended bodies.
We provide the formulation in which a quaternion defines each body's orientation, instead of the full rotation matrices or Euler angles. Quaternions are singularity-free, and their use was found by Fukushima (2008) to be most efficient in terms of numerical error. We provide some basic aspects of quaternions in Appendix A.1, and present the relations with Cayley-Klein parameters (see Sect. 2.3) in Appendix A.2.

Equations of motion
To describe the complete two-body dynamics, we propagate the translational and rotational state of both bodies. Our state vector x is defined as follows: where r i and v i denote the inertial position and velocity of body i. The term ω F i i denotes the angular velocity vector of body i w.r.t. the I frame, expressed in frame F i . The quaternion q I /F i describes the quaternion rotation operator from frame F i to frame I . The variable q I /F i i represents the vector containing the four entries of the quaternion (see Appendix A.1 for more details).
In the remainder, we omit the superscripts for ω F i i and q I /F i i (writing them as ω i and q i ), reintroducing an explicit frame notation only if it differs from the standard one in Eq. (53). We denote the quaternion that defines the full rotation as from F 2 to F 1 as q. We stress that neither the quaternion vector q nor the Cayley-Klein vector c (containing the entries of the complex numbers a and b, see Eq. (29)) represents a vector in the typical use of the term (quantity with magnitude and direction). In this context we use the more general definition of vector used in computer science, i.e., a container of numerical values.
Due to the symmetry of Eq. (17) in r 1 and r 2 , the translational equations of motion of the two bodies in an inertial frame are obtained immediately from Maciejewski (1995) as: The mutual force potential depends on r i through the spherical relative coordinates r, ϑ and ϕ, as shown in Eq. (20), which is identical to the case for the single-body potential.
Consequently, the calculation of the potential gradient can be done using standard techniques in space geodesy (Montenbruck and Gill 2000), facilitated by our relation between the one-body and two-body potentials in Sect. 2.4. The rotational dynamics is described by (e.g., Fukushima 2008): where we take the inertia tensor I i of body i in coordinates fixed to body i. For our application, we setİ = 0 (see Sect. 2.2). Theq = Qω formulation is used for the numerical propagation. For the computation of gravitational torques in the F2BP, Boué (2017) has very recently introduced a novel approach to compute the terms M I i , based on Varshalovich et al. (1988), which is computationally more efficient, and allows the torque to be evaluated without resorting to Euler angles. We express a single term V l 1 ,m 1 l 2 ,m 2 as in Eq. (26), which allows for expressing M F 1 2 as follows using the angular momentum operatorĴ : The termsĴ (D l 2 m 2 ,k 2 ) can be evaluated directly in Cartesian coordinates, from the expressions given by Boué (2017): where the terms D l m,k are computed as described in Sect. 2.3. The formulation of Eq. (61) can be evaluated using the approach of Sect. 2.4. The difference that needs to be introduced when evaluatingĴ subsequently replacing ( , ) l mk in Eq. (40) and (41) with ( , ) l mk providesĴ (C 2,F 1 lm ) andĴ (S 2,F 1 lm ), respectively. Continuing in the same manner, Eqs. (47) and (48) are adapted to obtainĴ (C l 1,2 ;m 1,2 ) andĴ (S l 1,2 ;m 1,2 ), respectively. Finally, we evaluate, analogously to Eq. (49): × Ĵ (C l 1,2 ;m 1,2 ) cos mϑ +Ĵ (S l 1,2 ;m 1,2 ) sin mϑ (67) To complete the rotational equations of motion, the torques M 1 and M 2 are related as: as a consequence of conservation of angular momentum,

Variational equations
To estimate the rotational and translational behavior in the F2BP from planetary tracking data, we need a formulation of the variational equations (Montenbruck and Gill 2000;Tapley et al. 2004;Milani and Gronchi 2010) for the dynamical model defined in Sect. 3.1. This requires the computation of the partial derivatives of the accelerations and torques w.r.t. the current state, as well as a parameter vector. These partial derivatives can be computed numerically, but an analytical formulation is often computationally more efficient and less prone to numerical error. This is especially true in the case of gravitational accelerations, for which the computation can be performed in a recursive manner based on the acceleration components. Typically, only the translational motion is estimated dynamically (i.e., represented in the state vector x) from the tracking data. The rotational behavior of the bodies is most often parameterized as a mean rotational axis and rate (possibly with slow time variations) in addition to a spectrum of libration amplitudes. These spectra may be obtained from fitting observations of solar system bodies (Archinal et al. 2018), or a numerical integration of the rotational equations of motion, as done by Rambaux et al. (2012) for Phobos, based on current models for the physical properties of the system. This approach partly decouples the translational and rotational dynamics in the estimation. A coupled determination of the initial rotational and translational state has been performed for only a limited number of cases, most notably in the case of the Moon using LLR data (e.g., Newhall and Williams 1996;Folkner et al. 2014;Hofmann 2017). Dynamical estimation of rotational motion of asteroid Bennu was performed by Mazarico et al. (2017) during a simulation study in preparation for the OSIRIS-REx mission. There, the possibly significant wobble of the asteroid is shown to require a dynamical approach to rotation characterization when performing accurate proximity operations. We propose a similar approach here, in which the full state vector x is dynamically determined, and all couplings are automatically included in the estimation model.
We stress that the approach of estimating libration amplitudes, rotation rates, etc., is the preferred choice for analyses of data from most current solar system missions. The signatures of the full rotational behavior are typically too weak in these data sets to warrant the full dynamical estimation (with the clear exception of lunar rotation from LLR data). However, such a full dynamical approach was shown to be crucial by  for the realistic data analysis of a Phobos Laser Ranging mission, both to ensure full consistency between all estimated parameters, and to prevent the estimation of an excessive number of correlated libration parameters.
To determine the full state behavior as a function of time, the state x at some time t 0 (denoted x 0 ), as well as a set of parameters represented here as p, are estimated. The influence of initial state and parameter errors are mapped to any later time by means of the state transition and sensitivity ma-trices, Φ(t, t 0 ) and S(t), defined as: The differential equations governing the time-behavior of these matrices are given by: where the state derivative modelẋ is defined by Eqs. (55)-(58). In Sects. 3.2.1 and 3.2.2 we present the detailed formulation of the terms in Eqs. (72) and (73), respectively.

State transition matrix
Writing out the partial derivatives in Eq. (72), we obtain four matrix blocks of the following structure (with i, j = 1...2): In this section, we will explicitly derive equations that enable an analytical evaluation of the partial derivatives. Table 1 provides an overview of the result of the derivation for the terms in the above matrix equations. Note that this approach is directly applicable to i, j > 2 as (in the absence of tides) any two-body interaction is independent of additional bodies in the system under consideration.
Since the equations of motion only contain the relative position r, not the absolute position r i , we have the following symmetry relation: Also, the symmetry expressed by Eq. (56), allows for the computation of the partial derivatives of v 2 from the associated partials ofv 1 as: No such symmetry exists for the partial derivatives w.r.t. q j , since the potential is explicitly dependent on both R F 1/F 2 and R F 1/I , but not R F 2/I . This is due to the fact that the angles ϑ and ϕ, representing the relative position of body 2 w.r.t. body 1, are expressed in a frame fixed to body 1.
Having defined the relevant symmetry relations, we derive expressions for the partial derivatives in Eq. (74). The partial derivatives ofv 1 w.r.t. positions r i are obtained from Eqs. (55): requiring the computation of the second derivatives of the potential components w.r.t. r F 1 , which we discuss later in this section. Note the post-multiplication with R F 1/I (in both Eq. (77) and Eq. (85)), which is due to the potential Hessian being computed w.r.t. the body-fixed position r F 1 , whereas the required partial derivative for Eq. (74) is required w.r.t. the inertial position r.
The expression for the partial derivatives ofv 1 w.r.t. the orientations q j are obtained from: where the derivatives of R I /F 1 are only non-zero for j = 1 (i.e., for the partial w.r.t. the orientation of body 1). Since the potential is written in terms of r F 1 , not the state variable r, the term ∂ ∂q j ( ∂V 1-2 ∂r F 1 ) is computed using: where * = ∂V 1-2 ∂r F 1 when evaluating Eq. (78). The crossderivative of V 1-2 is discussed later in this section. The partial derivatives ∂c/∂q i are computed from ∂c/∂q and ∂q/∂q i , which are explicitly given in Appendices A.1 and A.2, respectively.
The derivative ofω i w.r.t. ω i , which follows directly from Eq. (57), withİ i = 0: where the term [c] × is an anti-symmetric matrix used to represent the cross-product (as this will aid later derivations), and is defined as: For the partial derivatives ofω i w.r.t. the r j and q j , we obtain from Eq. (58): The partial derivative of M F 1 1 w.r.t. r are obtained from Eq. (68): The final term can be computed using the same procedure as For the final partial derivatives ∂M F i i /∂q j , none of the preceding symmetry relations can be used. From Eq. (68), we obtain the following for i = 1.
where again the derivatives of R F 1/I are non-zero only for j = 1. Note that the final term in Eq. (87) is computed using Eq. (79). Finally, for the partial derivatives of M F 2 2 , we have from Eq. (69): From Eqs. (75)-(89), we can compute the terms in Eq. (74), as summarized in Table 1. The following partial derivatives of the mutual potential are needed in the formulation: where each of the partial derivatives can be evaluated in a component-wise manner on V l 1 ,m 1 l 2 ,m 2 andĴ (V l 1 ,m 1 l 2 ,m 2 ). The calculation of the ∂ 2 V ∂(r F 1 ) 2 can be performed using well-known techniques in space geodesy, e.g. Montenbruck and Gill (2000). We obtain the derivatives of V l 1 ,m 1 l 2 ,m 2 w.r.t. the Cayley-Klein parameters c from Eqs. (28) and (31): where we have introduced the real column vectorM 2,F 1 l 2 ,m 2 = [ (M 2,F 1 l 2 ,m 2 ), (M 2,F 1 l 2 ,m 2 )] T . The derivatives ∂D l mk /∂c for l = 0, 1 are obtained directly from Eqs. (29) and (30).
The derivatives of the angular momentum operators are obtained similarly: where ∂D l m,k ∂c follows directly from Eqs. (63) and (94).
In addition to the partial derivatives of the potential terms, partial derivatives between the angle representations are required for Eqs. (77) finalizing the completely analytical model for the evaluation ofΦ of the full two-body gravitational interaction.

Sensitivity matrix
A formulation for ∂ẋ/∂p is required for the evaluation of Eq. (73). The parameter vector p can contain any physical parameter of the environment/system/observable. In the context of the F2BP, key parameters are the static gravity field coefficients of both bodies. Here, we provide a general formulation of ∂ẋ/∂p, for the case where entries of p directly influence (C,S) i,F i lm . By using such a formulation, the influence of model parameters on the full dynamics are accurately and consistently represented. This will prevent errors in the state transition/sensitivity matrices from propagating into biased estimates of the gravity field parameters. This is illustrated with test cases of Phobos and KW4 in Sect. 4.
The following partial derivatives, along with the symmetry relation in Eq. (76), allows for an analytical evaluation of ∂ẋ/∂p (again omitting the contribution byİ): The derivatives of the mutual potential and angular momentum operator require the computation of the terms ∂ ∂p ( ∂V l 1 ,m 1 l 2 ,m 2 ∂r F 1 ) and ∂ ∂p (Ĵ (V l 1 ,m 1 l 2 ,m 2 )), which are obtained from: The partial derivative of the inertia tensor is computed from: where I is related to the unnormalized gravity field coefficients by: HereĪ denotes the body's mean moment of inertia. Note thatĪ could be one of the entries in the vector p.

Numerical implementation
In this section we summarize the implementation of the algorithm. Our implementation has been done in an extended version (Dirkx 2015) of the Tudat 2 software toolkit, see Appendix C, a generic and modular astrodynamics toolbox written in C++. We focus on the implementation of the inner loops of the function evaluation (i.e., terms inside one or more summations). Summarizing, the main steps in the evaluation of the governing equations are (per time step): 1. Computation of Wigner D-matrices D l mk from Eqs. (30)-(35), from the current relative orientation of the two bodies, expressed by the Cayley-Klein parameters a and b. These are obtained directly from the components of q 1 and q 2 in the state vector. The relation with Cayley-Klein parameters is given in Appendix A.2. 2. Computation ofC 2,F 1 l 2 m 2 andS 2,F 1 l 2 m 2 (for all l 2 , m 2 ), by inserting D k lm into Eqs. (40) and (41). 3. Evaluation of effective spherical harmonic coefficients C l 1,2 ;m 1,2 andS l 1,2 ;m 1,2 for each combination of l 1 , m 1 , l 2 , m 2 from Eqs. (47) and (48).

Model results
We consider two test cases to illustrate our methodology, and to analyze the need for the use of figure-figure interactions in accurate ephemeris generation. First, we analyze the dynamics of Phobos, motivated by various recent analyses of high-accuracy tracking to future Phobos landers (e.g. Turyshev et al. 2010;Le Maistre et al. 2013;). Second, we apply our method to the dynamics of the KW4 double asteroid. This double asteroid has been the ubiquitous example in analyses of the figure-figure interactions. Our goal in this section is to ascertain the consequence of neglecting figure-figure effects during these bodies' state estimation, if high-accuracy tracking data were available. In our simulations, we introduce a difference between the truth model (used to simulate observations of the dynamics) and the estimation model (which is used to fit the observations). Our truth model includes the figure-figure interactions as in Sect. 3.1, while our estimation model does not. We consider the estimation of only the translational dynamics, and gravity field coefficients of the bodies, as is done in data analysis of current observations (see Sect. 3.2).
We use simulated observations of the full three-dimensional Cartesian state of the body under consideration (w.r.t. its primary). These observations cannot be realized in practice, but are most valuable in determining the sensitivity of the dynamics to various physical effects (in this case the figure-figure interactions). Specifically, it allows us to determine how much the figure-figure effects are absorbed into the estimation of other parameters when omitting these effects during the estimation. Mathematical details of this estimation approach are given by Dirkx et al. (2016). We use Fig. 1 Full two-body acceleration components acting on Phobos (index 1) due to its interaction with Mars (index 2), up to degree and order 2 noise-free observations, to ensure that any estimation error is due to the dynamical effects, not the observation uncertainty. The dynamical model during the estimation (e.g. without figure-figure effects) reduces to that used by, e.g., Lainey et al. (2004). This dynamical model is equivalent to our formulation, without the terms where both l 1 > 0 and l 2 > 0.

Phobos
We take the Phobos gravity field coefficients from Jacobson and Lainey (2014), who obtainC P 20 = −0.0473 ± 0.003 andC P 22 = 0.0229 ± 0.0006. The other Phobos gravity field coefficients are set to zero. We use the rotational model by Rambaux et al. (2012), which was obtained from numerical integration of rotational equations of motion, including the influence of figure-figure interactions. In their rotation model, the Phobos orbit is fixed to that produced by Lainey et al. (2007).
In our model, the dynamics of Phobos is numerically propagated including the figure-figure interactions between Mars and Phobos up to l 1 = m 1 = l 2 = m 2 = 2. In Fig. 1, the magnitude of the separate terms of the series expansion of the gravitational acceleration in Eq. (55) are shown over several orbits. The strongest figure-figure interactions (l 1 = l 2 = m 1 = 2, m 2 = 0) have a magnitude that is about 0.5% that of the weakest point-mass interaction (l 1 = 2, l 2 = m 1 = m 2 = 0). We use the Mars gravity field by Genova et al. (2016), and the Mars rotation model by Konopliv et al. (2006).
For geophysical analysis of Phobos, theC P 20 andC P 22 coefficients are of prime interest. Results of a consider covariance analysis by  have shown thatC P 22 could be determined to at least 10 −9 using a laser ranging system on a Phobos lander, operating over a period of 5 years (and close to 10 −7 after only 1 year). ForC 20 , they obtain errors of 10 −4 (relative error 0.2%) and 10 −7 (relative error 0.0002%) after 1 and 5 years, respectively. Note that they only consideredC 20 andC 22 for the Phobos gravity field. Here, the influence of omitting the figure-figure interactions during the estimation of these coefficients is quantified. Ideal observations of Phobos' position are simulated, with all interactions up to l 1 = l 2 = m 1 = m 2 = 2. Then, these simulated data are used to recover the initial state of Phobos, as well as its full degree two gravity field, using a dynamical model without any figure-figure interactions.
The results of the estimation are shown in Fig. 2, where the errors inC P 2m andS P 2m are shown as a function of the duration of the simulations. The errors inC P 20 andC P 22 due to neglecting figure-figure interactions in the estimation model are, at 2 × 10 −4 (relative error 0.42%) and 1.5 × 10 −5 (relative error 0.065%) respectively, much larger than the uncertainties obtained by . This unambiguously shows that precision tracking of a Phobos lander will require the use of the figure-figure interactions in both the dynamical modelling and estimation of orbit/physical parameters. Also, it shows that for the analysis of existing data of Phobos, neglecting figure-figure interactions is an acceptable assumption, as the errors we obtain are smaller by more than an order of magnitude than the formal uncertainties of these parameters reported by Jacobson and Lainey (2014) (shown in Fig. 2 as dashed lines).   Jacobson and Lainey (2014) our simulations, estimating the gravity field of Phobos up to degree and order 4, instead of 2.
Unfortunately, the resulting estimation problem becomes ill-conditioned, preventing results from being obtained. This indicates that the dynamics of Phobos does not contain the required information to independently estimate its full gravity field up to degree four. This is a general property for satellites in a (near-)circular and (near-)equatorial orbit, as the contribution of the degree two and degree four coefficients cannot be distinguished in the dynamics: both show the same temporal signature (with different magnitudes) in the observations, leading to ill-posedness in the estimation. The approach taken by e.g. Yoder et al. (2003) is to estimate 'lumped' gravity field coefficients, essentially acknowledging that an estimation of theC 2,0 term includes contributions fromC 4,0 ,C 6,0 , etc. As a consequence of this degeneracy, estimating higher-order gravity field coefficients will require tracking of a spacecraft around/near Phobos, as it cannot be fully estimated from Phobos' orbital dynamics alone.
This section shows that future precision tracking of Phobos, in particular in the case of landers, must take into account the figure-figure interactions up to at least the 2nd degree of Phobos' and Mars' gravity field. This article provides the required algorithms to achieve this. The values of the gravity field coefficients encode information concerning the mass distribution inside Phobos, provide insight into the presence or absence of voids and ices in the interior, and possible lateral density variations. Failure to include these effects into the estimation model can lead to erroneous interpretation of Phobos' interior structure (Le Maistre et al. 2019).

KW4
The KW4 binary asteroid system, which consists of two bodies termed Alfa and Beta, has been the 'standard' test case for the influence of figure-figure interactions. The shapes of the two asteroids were determined by Ostro et al. (2006) using radar imaging. Combined with determination of the mass of the two bodies, and the assumption of homogeneous interior mass distribution, this provides a model for both bodies' gravity field, facilitating the analysis of their dynamics. The full two-body dynamics of the system was studied in detail by Fahnestock and Scheeres (2008). The dynamics of this system was later used as a test case for models of the F2BP by (e.g., Boué and Laskar 2009;Compère and Lemaître 2014;Boué 2017;Hou 2018).
We use the spherical harmonic coefficients of the two bodies up to degree 4, as provided by Compère and Lemaître (2014). The body Beta (secondary) is propagated w.r.t. Alfa (primary), using the unexcited initial conditions given by Fahnestock and Scheeres (2008). We plot the resulting accelerations during several orbits in Fig. 3. In this plot, accelerations have been summed over both m 1 and m 2 in Eq. (55). This figure shows the significant influence that the figurefigure interactions have, with an acceleration magnitude of the strongest figure-figure interaction (total l 1 = l 2 = 2 terms) at about 4 × 10 −5 that of the primary (l 1 = l 2 = 0 term). In fact, the combined l 1 = l 2 = 2 terms are only about one order of magnitude weaker than the combined (l 1 = 0, l 2 = 2) terms. This exceptionally strong influence of figure-figure interactions is a consequence of the close orbit of the two bodies (R 1 /a ≈ 0.26, R 2 /a ≈ 0.09, with a the semi-major axis of the mutual orbit), and the highly irregular gravity fields of the two bodies. Using the same methodology as for Phobos, we perform two analyses for KW4, in which we estimate only the gravity field coefficients of either Alfa or Beta, respectively. Estimating the coefficients of both bodies simultaneously results in an ill-conditioned problem, indicating that the relative motion of the bodies alone cannot be used to constrain both bodies' gravity field coefficients. The underlying reason is similar to that for the ill-posedness encountered in the Phobos analysis: the contributions of, for instance, both bodies'C 2,0 coefficients cannot be independently determined from observations of the bodies' relative dynamics. Determining both bodies' gravity fields simultaneously will require spacecraft orbit/flyby observations, which would improve the conditioning of the problem. Such a mission design is outside of the scope of this article, and we limit ourselves to the relative motion of the two bodies.
Therefore, the results obtained here are not a direct indication of the attainable gravity field estimation accuracy during data analysis when omitting figure-figure effects, unlike the results for Phobos in Sect. 4.1. Instead, our results for KW4 are indicative of the relevance of figure-figure interactions, when attempting to infer the interior structure of the bodies from their relative motion.
With l max = 2 for both the simulation and estimation model, the resulting errors in the gravity field coefficients of Alfa and Beta (which are obtained in separate estimations) are shown in Fig. 4 with solid lines. As was the case for Phobos, significant errors in the estimated degree-2 gravity-field coefficients are obtained when omitting the figure-figure effects during estimation. Unlike the Phobos simulations, however, the magnitude of the errors does not always converge to a constant value for sufficiently long observation times.
The analysis where simulation of observations was done up to l max = 2 (with figure-figure interactions), with an esti-mation model up to l max = 4 (without figure-figure interactions) has also been performed. Unlike the case of Phobos, where the estimation problem failed to converge in this situation, the KW4 simulations consistently produce results. The errors in the degree two coefficients of both bodies are shown in Fig. 4 with dashed lines. Comparing to the case where l max = 2 during both observation simulation and estimation, there is a decrease of several orders of magnitude in the errors of the degree 2 gravity field coefficients. Moreover, the values of the estimation error show a much more stable behavior for long simulation times. Errors in the estimated degree-3 gravity-field coefficients are limited, at 10 −5 at most, and <10 −10 for many coefficients. Errors in degree 4 coefficients are significant, however, in particular for Beta, where the absolute errors inC 40 andC 42 reach values of 0.075 and 0.012, respectively. For Alfa, these errors are at the level of 10 −3 and 2 × 10 −4 . Note that Fig. 4 only shows the errors in degree 2 gravity field coefficients, not in the coefficients of degree 3 or 4.
These results indicate that the signature of figure-figure interactions at degree 2 can be partly absorbed by the estimation of gravity field coefficients of degree 4, when omitting figure-figure effects during estimation. In this situation, the degree 4 coefficients are used for a similar purpose as empirical accelerations, which are typically used in spacecraft orbit determination (Montenbruck and Gill 2000): to absorb part of the dynamical mismodelling, and reduce the degree to which these model errors impact the estimation quality.
The impact that the estimation of higher-degree coefficients has on the attainable estimation accuracy will be strongly dependent on the system under consideration, though. Also, applying this method requires the estimation of a substantial number of extra parameters, potentially weakening the stability of the numerical solution. We stress that, even in cases where this approach could be robustly used to estimate degree-2 gravity-field coefficients, the resulting degree-4 coefficients will not have a physical meaning, as they have absorbed part of the signature of the figurefigure effects. When using our methodology outlined in this article for the estimation, the coefficients will not absorb dynamical model errors, and the estimated degree-4 coefficients would be representative of the bodies' actual mass distribution.
In summary, our results indicate that the determination of the gravity fields of Alfa and Beta based (in part) on the observed mutual dynamics of the two bodies should take into account figure-figure interactions. Failure to do so results in substantial errors in the degree-2 or degree-4 coefficients, depending on the degree to which the gravity fields are expanded. The attainable accuracy of gravity field coefficient estimation is strongly dependent on tracking data types/uncertainties, lander/orbiter mission geometry, etc. For a specific mission design, our methodology can be used to quickly assess which effects should and should not be considered during the data analysis.

Conclusions
We have derived in detail the equations that relate the full two-body gravitational interaction, including all figurefigure effects, to the typical one-body gravity field representation in terms of (fully-normalized) spherical harmonic coefficients. Our derivation provides the explicit link between the elegant representations found in literature and the need for a practical implementation in existing state propagation and estimation software, in particular for applications in solar system dynamics and orbit determination.
Our framework enables robust modelling of the full gravitational interactions in both simulation studies and data analysis. The approach will be crucial for future tracking techniques, such as ILR and SBI, in which (sub-)mm relative model accuracy between the positions of celestial bodies is required over a period of years. In addition to capturing the full dynamics, the formulation allows the rotational dynamics of a celestial body to be estimated in a manner that prevents the estimation of a broad libration spectrum, and mitigates problems associated with decoupled rotationaltranslational dynamics in the estimation model. We provide explicit formulations for transformed spherical harmonic coefficients, given by Eqs. (39)-(41), the effective coefficients for use in the one-body potential formulation in Eqs. (47)-(48), and the explicit transformation from two-body to one-body potential terms in Eqs. (42)- (46) and (49). Associated equations for unnormalized spherical harmonic coefficients are given in Appendix B. These formulations are crucial to transparently implement the F2BP in existing tracking data analysis suites.
To be able to apply our model to state estimation of natural bodies, we have derived an analytical formulation of the variational equations in Sect. 3.2. The resulting differential equations for the state transition matrix Φ(t, t 0 ) and sensitivity matrix S(t) allow the direct use of the models in existing state propagation and estimation software, and automatically capture any dependency of the mutual dynamics on either body's state, gravity field coefficient, or other physical parameter. The use of analytical partial derivatives is crucial to prevent an excessive computation load. The required equations are summarized in Table 1, and have not been found comprehensively in literature for any representation of the F2BP.
The relevance of our method is predicated on the assumption that figure-figure interactions will be relevant for data analysis of (future) planetary missions. We have analyzed the impact of ignoring figure-figure interactions during state and gravity-field estimation (as is typical standard practice) for both Phobos and the KW4 binary asteroid. For Phobos, relative errors in the estimated values of theC 20 andC 22 coefficients are at 0.42% and 0.065%, respectively. These values are below the current uncertainty in Phobos' gravity field, validating the omission of figure-figure effects in previous studies. It was shown by  that accurate laser tracking of a Phobos lander will allow these coefficients to be estimated to much greater accuracy, requiring the figure-figure effects to be included in the estimation model. The same will be true for Doppler tracking of a Phobos lander (Le Maistre et al. 2013). A similar analysis of the mutual dynamics of the KW4 binary asteroid system was performed, again indicating significant errors in the gravity field coefficient estimation when omitting figure-figure effects. The errors in the degree two coefficients are reduced significantly, to the range 10 −9 -10 −6 when increasing the degree to which gravity field coefficients are estimated from 2 to 4, while still omitting figure-figure interactions in the estimation. However, this is at the expense of significant errors in the degree 4 gravity field coefficients.
Each of the existing formulations of the F2BP, for instance in terms of mass multipoles (Compère and Lemaître 2014), inertia integrals ) and polyhedrons (Fahnestock and Scheeres 2008), has its specific (dis)advantages. Our formulation for the F2BP in terms of spherical harmonic coefficients is purposefully designed for the use in space-mission tracking data analysis, maximizing compatibility with currently typical approaches to state and gravity field parameter estimation.
In this appendix, we discuss the relation between the angle representations we use: 3-1-3 Euler angles, rotation matrices and quaternions, discussed in greater detail by Diebel (2006). Additionally, we derive explicit formulations for the partial derivatives between these representations that are needed in the evaluation of the full two-body gravitational problem.

A.1 Quaternions
We nominally use quaternions to represent the orientation of celestial bodies. Quaternions provide a singularity-free representation, while ensuring that the minimum number of required variables are propagated (Hughes 2004). Moreover, propagating the quaternions was shown by Fukushima (2008) to be the optimal choice (of the broad range of options that they considered) in terms of numerical integration errors. The quaternion has four entries, denoted as q 0 , q 1 , q 2 and q 3 , which can be directly related to an angle-axis transformation of an angle θ about a unit-axisn through: where the first term denotes the q 0 entry, and the subsequent vector the q 1 , q 2 and q 3 terms. We again use q for the vector representation of the quaternion and q for the operator representation. For a quaternion representing a rotation, the norm is exactly 1, so: Consequently, the first and second derivative of a quaternion w.r.t. to any parameter α must satisfy: Defining a new quaternion operator as: and denoting a quaternion as per Eq. (107), we obtain the following: using the notation of Eqs. (81) and (107). This allows partials of a composite quaternion to be determined from its constituents' partials.

A.2 Cayley-Klein parameters
Cayley-Klein parameters are directly related to quaternions by the following: where we use the same sign-convention as Boué (2017). From here, the derivatives w.r.t. q are directly obtained from:

A.3 Rotation matrix partial derivatives
The evaluation of the rotational equations of motion require the determination of the terms ∂R ij /∂α for a number of quantities α (where R ij denotes the entry of the rotation matrix R). However, as was discussed for quaternions in Appendix A.1, the entries R ij are not independently chosen.
In particular, the columns R i must satisfy the following six relations (e.g., Maciejewski 1995): Consequently, the derivatives of the rotation matrices w.r.t. some quantity α must obey: Here we choose the terms R 13 , R 23 and R 31 to be independent, 3 which we collectively term P. The vector of the remaining six dependent entries is denoted R . The vector of dependent partial derivatives, denoted ∂R ∂α , is now obtained from the solution of: where the matrix A is expressed in terms of P and R , and b is a function of R and ∂P ∂α . The entries of the matrix A represent and the vector b are directly obtained from Eq. (118) for each of the six equations in Eq. (118). The above can be used to obtain ∂R ij /∂q k from the corresponding derivatives of P.

B.2 Gravity field coefficient transformations
The explicit formation of the transformed unnormalized coefficients are obtained from Eqs. (8) and (28). The result for normalized coefficients was given in Eqs. (40) and (41). The corresponding equations for unnormalized coefficients are obtained byC → C,S → S andν lmk → ν lmk , with: The expressions for the normalized equivalent of Eqs. (45) and (46) are then obtained directly from Eqs. (15) and (20). The resulting equations are obtained by replacing both the coefficients C, S, C, S and the coefficient γ by their normalized coefficients (denoted by the straight overbar). These equations enable the computation of the mutual force potential terms using subroutines for the one-body normalized spherical harmonic potential in Eq. (14). This is achieved by applying the scaling rules in Eqs. (42)-(46), and replacing the effective one-body coefficients with the normalized counterparts from Eqs. (47) and (48).

Appendix C: Tudat software
The Tudat software, which was used to generate the results presented in this article, is a multi-purpose, modular, astrodynamics software suite written in C++, developed by staff and students of the Asrodynamics & Space Missions group of Delft University of Technology. The project was started in 2010, with the current architecture having been set up in 2015. Aspects of the software are discussed by Dirkx (2015), with detailed feature documentation (including installation guide and tutorials) at tudat.tudelft.nl, and code documentation at doxygen.tudat.tudelft.nl. The code is licensed under the BSD 3-Clause "New" or "Revised" License, and freely available from github.com/tudat. All Tudat functionality is tested in an associated (unit) test. All tests can be rerun at will by users to verify the integrity of the code after installation and/or modification.
In this appendix, we recall the main features of the Tudat software suite, limiting ourselves to orbit propagation and estimation functionality (omitting various other features related to, e.g., mission design).

C.1 State propagation
A strong design driver of the Tudat software is modularity: the software architecture makes no a priori assumptions on bodies being considered, or on models being used for physical properties of bodies, accelerations, etc. Instead, users are free to set up the simulation according to their own needs, choosing from a broad variety of models. To make implementation straightforward, various default models are provided (which can be overridden at will), for instance for solar system ephemerides, gravity fields and rotational models. Tudat makes no distinction between natural or artificial bodies, this distinction is introduced purely by the properties assigned to a given body.
Tudat is applicable to a broad range of topics, and has been used for, among others, simulations of solar system dynamics (Dirkx et al. 2018), interplanetary trajectory design (Musegaas 2013), planetary system dynamics (Kumar et al. 2015;Dirkx et al. 2016), space debris impact predictions (Ronse and Mooij 2014) and atmospheric re-entry (Dirkx and Mooij 2014).
Setting up a numerical state propagation in Tudat consists of defining the following models: • Environment models: all properties of bodies (ephemerides, gravity fields, shapes, atmospheres, etc.) are stored in the environment. In the Tudat architecture, this includes properties of artificial bodies, such as engine models. • Dynamical model type(s) to be propagated. Tudat is currently capable of numerically propagating a body's translational state, rotational state and mass (a custom state derivative model is also provided to provide flexibility for other types of dynamics). A combination of any or all of these types of dynamics, for any number of bodies may be provided. The type of dynamics need not be the same for each body, so that the translational state and mass of a spacecraft, and translational state and rotational state of a planet, may be propagated. The dynamics may be propagated hierarchically (e.g. spacecraft w.r.t. Moon, Moon w.r.t. Earth and Earth w.r.t. Barycenter) in a concurrent fashion. • Dynamical model formulation. Depending on the dynamics type, multiple formulations for the governing equations may be used, such as Cartesian, Keplerian and Modified Equinoctial Elements for translational dynamics, and quaternions or modified Rodrigues parameters, in combination with angular velocities, for rotational dynamics. States may be propagated in a single-arc or multi-arc fashion, as well as a combination of the two, as described by (Dirkx et al. 2018). • State derivative models. Depending on the types of dynamics being propagated, models for, e.g., torques and accelerations must be provided. In Tudat, these models are defined by their type (e.g. aerodynamic, spherical harmonic gravity, third-body point-mass gravity, radiation pressure), the body exerting and body undergoing the torque/acceleration and, if needed, additional information (such as maximum degree/order for a spherical harmonic acceleration). • Numerical integrator settings: Tudat provides implementations of fixed-and variable time-step multi-stage methods (Runge-Kutta-Fehlberg), multi-step methods (Adams-Bashforth-Moulton) and extrapolation methods (Bulirsch-Stoer). • Output settings. By default, Tudat outputs only the numerically propagated states, but a wide array of additional variables may be saved in addition. These include (but are not limited to), acceleration/torque terms, local atmospheric conditions, relative orientations of bodies, etc. Such settings may also be used to define termination conditions for a numerical propagation, for instance when two bodies get within a certain user-defined distance.

C.2 State estimation
In addition to the numerical propagation modules of Tudat, state estimation (using a batch least square filter) can also be performed. This functionality was used by Bauer et al. (2016) for the orbit determination for LRO using laser ranging data. Simulation models for various astrometric and radiometric data types are available, and applied by e.g. Dirkx et al. ( , 2016Dirkx et al. ( , 2017Dirkx et al. ( , 2018, with recent publications focussing on the JUICE-PRIDE radio tracking experiment (Gurvits et al. 2013). Models for the detailed analysis of real radio and optical data are currently under development.
In addition to the settings listed above for the numerical propagation, Tudat requires the following information to perform a state estimation: • A list of parameters to estimate, which may include initial rotational and/or translational state (single-arc, multi-arc, or a combination of the two), gravity field coefficients, Parametric Post-Newtonian (PPN) parameters, etc. A user may, but need not, provide an a priori covariance matrix for these parameters. • Observation models and settings. Various types of observation models (such as range, range-rate, angular position) can be modelled by Tudat, where a user has the freedom to add bias models, atmospheric, relativistic correc-tions, and other settings that may be needed, such as integration time for closed-loop Doppler observables. • Observations and associated weights. The input for the estimation is a list of realizations of observations, with associated time tag and weights. These observations may be simulated by Tudat, or obtained from data archived (e.g. PDS) or other simulation tools.
Using the above information, acceleration and observation partial derivative models are automatically set up, and the associated variational equations for Φ(t, t 0 ) and S(t) are solved numerically.

C.3 External libraries
In addition to its own codebase, Tudat links to a number of external software libraries, primarily: • Boost: The Boost libraries 4 are a collection of C++ libraries, which is used in support of various Tudat functions, primarily unit testing, multi-dimensional arrays, file reading/writing and file system access. • CSPICE: The SPICE toolbox 5 (Acton 1996) is used in Tudat primarily to retrieve ephemerides, rotational states, and other physical quantities, such as gravitational parameters and radii of solar system bodies. • SOFA: The Standards of Fundamental Astronomy (SOFA) library 6 is used in Tudat primarily for Earth orientation and time conversion functionality • Eigen: Tudat uses the Eigen library 7 for all its linear algebra operations • Pagmo2: Tudat provides an optional interface to the Pagmo2 software library, 8 a global optimization toolbox developed by ESA's Advanced Concept Team (ACT).
When cloning our "Tudat bundle" repository, 9 Tudat and all other required libraries are automatically downloaded.