Modelling Rigid Body Potential of Small Celestial Bodies for Analyzing Orbit–Attitude Coupled Motions of Spacecraft

Abstract

The present study aims to propose a general framework of modeling rigid body potentials (RBPs) suitable for analyzing the orbit–attitude coupled motion of a spacecraft (S/C) near small celestial bodies, regardless of gravity estimation models. Here, ‘rigid body potential’ refers to the potential of a small celestial body integrated across the finite volume of an S/C, assuming that the mass of the S/C has no influence on the motion of the small celestial body. First proposed is a comprehensive formulation for modeling the RBP including its associated force, torque, and Hessian matrix, which is then applied to three gravity estimation models. The Hessian of potential plays a crucial role in calculating the RBP. This study assesses the RBP via numerical simulations for the purpose of determining proper gravity estimation models and seeking modeling conditions. The gravity estimation models and the associated RBP are tested for eight small celestial bodies. In this study, we utilize distance units (DUs) instead of SI units, where the DU is defined as the mean radius of the given small celestial body. For a given specific distance in Dus, the relative error of the gravity estimation model at this distance has a similar value regardless of the small celestial body. However, the difference value between the potential and RBP depends on the DU; in other words, it depends on the size of the small celestial body. This implies that accurate gravity estimation models are imperative for conducting RBP analysis. The overall results can help develop a propagation system for orbit–attitude coupled motions of an S/C in the vicinity of small celestial bodies.

1. Introduction

After the initial success of the Hayabusa mission, most deep space missions exploring small celestial bodies have opted to include proximity operations [1,2,3,4,5,6,7,8]. Given the distinctive dynamical characteristics of small celestial bodies due to their light masses, irregular shapes, and potentially variable spin axes/rates, proximity operations necessitate a meticulous analysis of the dynamical environment near small celestial bodies for successful missions.

Among a variety of dynamic characteristics, particularly noteworthy are the orbit–attitude interactions pertinent to gravity. The interaction caused by gravity becomes more influential as the orbital radius becomes smaller and the size of the spacecraft (S/C) becomes larger [9]. Given that more than 95% of these asteroids have diameters less than 1 km [10,11], these interactions are not only conspicuous but also substantially significant. Hence, this study primarily addresses the interaction between the orbit and attitude motions of an S/C induced by gravity.

When modeling the orbit–attitude coupled motion of an S/C, many studies have separately modeled orbit and attitude dynamics and combined them for simplicity of analysis, or they have employed mutual potentials. Most common is the former approach, which does not require any attitude information for analyzing S/C orbital motion but does require S/C positions to calculate gravitational torques [12]. This approach is applicable when the primary body significantly outweighs the S/C, as the contribution of the S/C volume to the gravitational force is negligible. Given the diminutive size of small celestial bodies, the contribution of S/C volume often becomes non-negligible. Though some recent studies discussed the influences of S/C attitude on orbital motions [13,14], they do not ensure that the total energy of the orbit–attitude coupled motions remain constant, since orbit and attitude motions are developed separately. As this study also aims to formulate the orbit–attitude coupled dynamics that naturally conserve total energy and momentum. The former approach is not suitable for our purposes. In contrast to the former approach, the latter employs the full two-body problem with mutual potential [15,16,17]. It addresses the orbit–attitude motions of both the primary body and the S/C, enabling the conservation of total energy within the entire system. Unnecessary computations arise, however, when the motion of small celestial bodies is not of primary concern. The present study focuses on the orbit–attitude coupled motions of an ‘S/C’ in proximity to small celestial bodies. In this context, it is important to differentiate between the ‘rigid body potential (RBP)’ and the ‘mutual potential’. Any motions of small celestial bodies themselves are not of our concern. The S/C is assumed to be an extended rigid body with finite volume; the ‘restricted’ assumption holds, i.e., the S/C does not exert any gravitational force on the small celestial body. With these assumptions, the interplay between the orbit and attitude of an S/C can be modeled through the process of volume integration using the gravity estimation model.

A multitude of research has been carried out with regard to RBP estimations in the vicinity of small celestial bodies. These studies are differentiated based on how to approximate the body of an S/C and how to estimate gravity. With regard to the methods for approximating the S/C body, some studies have developed higher-order integrations for the S/C body [18], while others rely on using the moment of inertia (MOI) for extended bodies [19]. In this study, the latter approach is selected, wherein the S/C is treated as an extended body. This choice is made to leverage its extensive applicability within dynamic and control systems.

The main contributions of our present study can be summarized as follows: (1) We present a comprehensive formulation of the RBP, including its associated force and torque components. The proposed framework of the RBP is applicable regardless of the gravity estimation models. (2) We apply this comprehensive formulation to three gravity estimation models, which can be naturally extended into other gravity estimation models. The estimation process of the RBP and its resultant force/torque requires the Hessian matrix of potential to be derived through each gravity estimation model. (3) This study assesses the RBP calculated using the three gravity estimation models. The primary objectives of this assessment are to determine the most suitable model for proximity operations in the vicinity of small celestial bodies and to determine conditions of applicability such as small body size. The numerical simulation reveals that the efficacy and precision of the RBP are contingent upon the chosen gravity estimation models, in addition to the accuracy of the Hessian calculations.

The remaining discussion is composed of five main sections. Section 2 introduces the background of this study and three gravity estimation models. Section 3 presents two propositions, encompassing a comprehensive formulation for the RBP, along with the related force and torque, and the second partial derivatives of the RBP concerning position. By amalgamating the contents from Section 2 and Section 3, Section 4 furnishes the RBP, accompanied by its corresponding force and torque, with three distinct gravity estimation models. Section 5 evaluates the associated potentials and RBPs. Section 6 summarizes and concludes the overall discussion.

2. Background and Gravity Estimation Models

2.1. Background

This study focuses on analyzing the orbit and attitude motions of an S/C orbiting small celestial bodies. If the masses/sizes of both the S/C and a small celestial body are similar to each other, the motions of both of them should be analyzed. In this case, it is convenient to select their barycenter as the coordinate origin. However, if the S/C is significantly smaller and lighter than a small celestial body, the gravitational influences of the S/C on the small celestial body is negligible. In such instances, selecting the center of the small celestial body as the coordinate origin becomes more convenient. This study considers the ‘in-between’ of two aforementioned cases; it is assumed that the mass of the S/C is considerably light, but its size is not negligibly small, compared to the small celestial body. There is typically an increase in the masses of small celestial bodies as their sizes grow, while the exact values vary based on their density [20]. In this section, it is confirmed that the ‘restricted’ assumption still holds, i.e., the S/C does not exert any gravitational force on the small celestial body based on existing exploration missions.

Table 1. The ratio of the masses of the S/C to the asteroid, the distance between the barycenter and the center of small celestial bodies, and the ratio of the acceleration caused by the S/C to that caused by the Sun [21,22,23,24,25,26].

Index #	Name	$\mathit{a}$ (AU)	$\mathit{e}$	${\mathit{R}}_{\mathit{S}}$ (m)	${\mathit{M}}_{\mathit{S}/\mathit{C}}$ (kg)	${\mathit{M}}_{\mathit{S}/\mathit{C}}/{\mathit{M}}_{\mathit{S}}$	$‖{\mathit{x}}_{\mathit{B}}-{\mathit{x}}_{\mathit{S}}‖$ (m)	$‖{\mathit{a}}_{\mathit{S}/\mathit{C}}‖/‖{\mathit{a}}_{⨀}‖$
4	Vesta	$2.36$	$0.089$	$2.88\times {10}^5$	$1108$	$4.28\times {10}^{-18}$	$3.69\times {10}^{-12}$	$9.95\times {10}^{-16}$
243	Ida	$2.86$	$0.043$	$2.76\times {10}^4$	$2717$	$2.72\times {10}^{-14}$	$2.25\times {10}^{-9}$	$3.57\times {10}^{-13}$
433	Eros	$1.46$	$0.22$	$1.60\times {10}^4$	$805$	$1.20\times {10}^{-13}$	$5.78\times {10}^{-9}$	$1.12\times {10}^{-13}$
951	Gaspra	$2.21$	$0.17$	$1.08\times {10}^4$	$2717$	$7.60\times {10}^{-13}$	$2.45\times {10}^{-8}$	$1.76\times {10}^{-12}$
1036	Ganymed	$2.67$	$0.53$	$3.43\times {10}^4$	$3000$*	$1.80\times {10}^{-14}$	$1.85\times {10}^{-9}$	$4.81\times {10}^{-13}$
1620	Geographos	$1.24$	$0.34$	$3.13\times {10}^3$	$805$	$2.01\times {10}^{-10}$	$1.89\times {10}^{-6}$	$2.55\times {10}^{-12}$
4179	Toutatis	$2.54$	$0.63$	$2.54\times {10}^3$	$1750$	$3.50\times {10}^{-11}$	$2.67\times {10}^{-7}$	$5.23\times {10}^{-11}$
4769	Castalia	$1.06$	$0.48$	$9.85\times {10}^2$	$805$	$5.75\times {10}^{-10}$	$1.70\times {10}^{-6}$	$2.30\times {10}^{-11}$
25143	Itokawa	$1.32$	$0.28$	$4.37\times {10}^2$	$510$	$1.45\times {10}^{-8}$	$1.91\times {10}^{-5}$	$8.57\times {10}^{-11}$
99942	Apophis	$0.92$	$0.19$	$9.96\times {10}^2$	$2110$	$3.46\times {10}^{-8}$	$1.03\times 10^{-4}$	$2.88\times {10}^{-11}$

* The selection of this value is arbitrary.

lists the physical properties of exploration/flyby missions targeting small celestial bodies [21,22,23,24,25,26]. The third, fourth, and fifth columns indicate the semi-major axis $a$, the eccentricity $e$, and the mean radius $R_S$ of a small celestial body, respectively. $M_{S/C}$ and $M_S$ denote the masses of the S/C and small celestial body, respectively. The third column from the right side of Table 1 shows the ratio of the S/C mass to that of the asteroid. The second column from the right side indicates the distance between the barycenter ${\mathit{x}}_B$ and the center ${\mathit{x}}_S$ of the small celestial bodies when the distance between the S/C and the small celestial body is $3R_S$. The maximum and minimum distances are $1.03\times {10}^{-4} \mathrm{m}$ and $3.69\times {10}^{-12 }\mathrm{m}$, respectively. These maximum and minimum values are $1.01\times 10^{-7}$ times and $1.28\times 10^{-17}$ times the average radius of the target bodies $R_S$, respectively. Even the maximum value is negligibly small when considering the difference in orders of magnitude between this value and the position values of the S/C. The first column from the right indicates the magnitude ratio of the acceleration ${\mathit{a}}_{S/C}$ caused by the S/C to the acceleration ${\mathit{a}}_{⨀}$ caused by the Sun when the small celestial bodies are located at their apoapsis. The maximum ratio is $2.30\times {10}^{-11}$, indicating that the gravitational force of the Sun predominates over the motions of small celestial bodies, rather than the gravitational force of the S/C.

2.2. Gravity Estimation Models

This section introduces three gravity estimation models used in this study: the Point Mass (PM) model, Extended Body (EB) model [27], and Triaxial Ellipsoid (TE) model [28]. All gravity estimation models are built upon the shape model of a small celestial body. Although these three models offer rather simple representations of the dynamical environments surrounding small celestial bodies, they still serve as valuable examples for elucidating the application of the RBP.

2.2.1. Point Mass Model

The PM model characterizes small celestial bodies as point masses, effectively representing restricted two-body motion. Let $\overrightarrow{\mathit{x}}\in {\mathbb{R}}^3$ be a position vector of a S/C. The potential $V^{PM}\in \mathbb{R}$ of a small body, its gradient ${\overrightarrow{\mathit{V}}}_x^{PM}\in {\mathbb{R}}^3$, and Hessian matrix ${\mathit{V}}_{xx}^{PM}\in {\mathbb{R}}^{3\times 3}$ are expressed as

$$V^{PM}=\frac{{\mu }_s}{‖\overrightarrow{\mathit{x}}‖},$$

$${\overrightarrow{\mathit{V}}}_x^{PM}=-\frac{{\mu }_s}{{‖\overrightarrow{\mathit{x}}‖}^3}\overrightarrow{\mathit{x}}={\overrightarrow{\mathit{f}}}_V^{PM},$$

(1)

$${\mathit{V}}_{xx}^{PM}=\frac{{\mu }_s}{{‖\overrightarrow{\mathit{x}}‖}^3}\left(3\widehat{\mathit{x}}{\widehat{\mathit{x}}}^T-{\mathbb{I}}_3\right)$$

(2)

where ${\mu }_s\in \mathbb{R}$ and ${\mathbb{I}}_3\in {\mathbb{R}}^{3\times 3}$ are a standard gravitational parameter and the identity matrix, respectively. $\widehat{\mathit{x}}$ and $‖\overrightarrow{\mathit{x}}‖$ denote a unit vector and the Euclidian norm for a given vector $\overrightarrow{\mathit{x}}$, respectively. The gravitational torque ${\overrightarrow{\mathit{\tau }}}_V^{PM}$ is given by

$${\overrightarrow{\mathit{\tau }}}_V^{PM}=3\frac{{\mu }_s}{{‖\overrightarrow{\mathit{x}}‖}^3}{\left({\mathit{R}}^T\widehat{\mathit{x}}\right)}^{\times }\mathcal{J}{\mathit{R}}^T\widehat{\mathit{x}}$$

where $\mathcal{J}\in {\mathbb{R}}^{3\times 3}$ and $\mathit{R}\in SO\left(3\right)$ are the MOI and attitude of S/C, respectively. ${\left(\cdot \right)}^{\times }:{\mathbb{R}}^3\to {\mathbb{R}}^{3\times 3}$ is a skew-symmetry operator and ${\left(\cdot \right)}^T:{\mathbb{R}}^{3\times 3}\to {\mathbb{R}}^{3\times 3}$ is a transpose operator.

2.2.2. Extended Body Model

The EB model utilizes the MOI associated with a small celestial body to depict its geometric characteristics [27]. It is not universally embraced for estimating the gravity of small bodies. Let ${\mathcal{J}}_s\in {\mathbb{R}}^{3\times 3}$ and $G\in \mathbb{R}$ be the MOI of a small celestial body and the gravitational constant, respectively. The potential $V^{EB}\in \mathbb{R}$ pertaining to the small celestial body is expressed as

$$V^{EB}=\frac{{\mu }_s}{‖\overrightarrow{\mathit{x}}‖}+\frac{G}{2{‖\overrightarrow{\mathit{x}}‖}^3}\left(tr\left[{\mathcal{J}}_s\right]-3{\widehat{\mathit{x}}}^T{\mathcal{J}}_s\widehat{\mathit{x}}\right),$$

where $tr\left[\cdot \right]:{\mathbb{R}}^{3\times 3}\to \mathbb{R}$ is a trace operator. The successive derivatives of $V^{EB}$ with respect to $\overrightarrow{\mathit{x}}$ yield both its gradient ${\overrightarrow{\mathit{V}}}_x^{EB}\in {\mathbb{R}}^3,$ and its Hessian matrix ${\mathit{V}}_{xx}^{EB}\in {\mathbb{R}}^{3\times 3}$ are formulated as

$${\overrightarrow{\mathit{V}}}_x^{EB}=-\frac{{\mu }_s}{{‖\overrightarrow{\mathit{x}}‖}^3}\overrightarrow{\mathit{x}}-\frac{3}{2}\frac{G}{{‖\overrightarrow{\mathit{x}}‖}^5}\left\{tr\left[{\mathcal{J}}_s\right]{\mathbb{I}}_3+2{\mathcal{J}}_s-5{\widehat{\mathit{x}}}^T{\mathcal{J}}_s\widehat{\mathit{x}}\right\}\overrightarrow{\mathit{x}}={\overrightarrow{\mathit{f}}}_V^{EB},$$

$${\mathit{V}}_{xx}^{EB}=3\frac{{\mu }_s}{{‖\overrightarrow{\mathit{x}}‖}^7}\left\{{‖\overrightarrow{\mathit{x}}‖}^2{\mathbb{I}}_3+\frac{5}{2}\frac{G}{M_s}\left(tr\left[{\mathcal{J}}_s\right]{\mathbb{I}}_3+4{\mathcal{J}}_s-7{\widehat{\mathit{x}}}^T{\mathcal{J}}_s\widehat{\mathit{x}}{\mathbb{I}}_3\right)\right\}\overrightarrow{\mathit{x}}{\overrightarrow{\mathit{x}}}^T$$

$$-\frac{{\mu }_s}{{‖\overrightarrow{\mathit{x}}‖}^5}\left\{{‖\overrightarrow{\mathit{x}}‖}^2{\mathbb{I}}_3-\frac{3}{2}\frac{G}{M_s}\left(tr\left[{\mathcal{J}}_s\right]{\mathbb{I}}_3+2{\mathcal{J}}_s-5{\widehat{\mathit{x}}}^T{\mathcal{J}}_s\widehat{\mathit{x}}{\mathbb{I}}_3\right)\right\}$$

where $M_s={\mu }_s/G$ is the mass of the small body.

2.2.3. Triaxial Ellipsoid Model

The TE model approximates a small celestial body as a triaxial ellipsoid, characterized by dimensions of $2\alpha \times 2\beta \times 2\gamma$ ($\alpha \ge \beta \ge \gamma$). Here, the dimensional parameters $\alpha$, $\beta$, and $\gamma$ are the semi-major, intermediate, and semi-minor axes of the approximated ellipsoid, respectively. The potential is constructed to account for a triaxial ellipsoidal body with uniform density [28]. Given specific values for $\alpha$, $\beta$, $\gamma$, and $M_s$, the potential $V^{TE}\in \mathbb{R}$, its gradient ${\overrightarrow{\mathit{V}}}_x^{TE}\in {\mathbb{R}}^3$, and its Hessian matrix ${\mathit{V}}_{xx}^{TE}\in {\mathbb{R}}^{3\times 3}$ are formulated as

$$V^{TE}=\frac{3{\mu }_s}{4}{\int }_{u^{\prime }\left(\overrightarrow{\mathit{x}}\right)}^{\infty }\varphi \left(\overrightarrow{\mathit{x}};u\right)\frac{du}{\mathsf{\Delta }\left(u\right)},$$

$${\overrightarrow{\mathit{V}}}_x^{TE}=\frac{3{\mu }_s}{4}{\int }_{u^{\prime }\left(\overrightarrow{\mathit{x}}\right)}^{\infty }{\overrightarrow{\mathit{\varphi }}}_x\left(\overrightarrow{\mathit{x}};u\right)\frac{du}{\mathsf{\Delta }\left(u\right)}={\overrightarrow{\mathit{f}}}_V^{TB},$$

(3)

$${\mathit{V}}_{xx}^{TE}=\frac{3{\mu }_s}{4}{\int }_{u^{\prime }\left(\overrightarrow{\mathit{x}}\right)}^{\infty }{\mathit{\varphi }}_{xx}\left(\overrightarrow{\mathit{x}};u\right)\frac{du}{\mathsf{\Delta }\left(u\right)}-\psi \left(\overrightarrow{\mathit{x}};u^{\prime }\right){{\overrightarrow{\mathit{\varphi }}}^{\prime }}_x\left(\overrightarrow{\mathit{x}};u^{\prime }\right){{{\overrightarrow{\mathit{\varphi }}}^{\prime }}_x}^T\left(\overrightarrow{\mathit{x}};u^{\prime }\right)$$

where

$$\varphi \left(\overrightarrow{\mathit{x}};u\right)=\frac{x_1^2}{{\alpha }^2+u}+\frac{x_2^2}{{\beta }^2+u}+\frac{x_3^2}{{\gamma }^2+u}-1,$$

(4)

$${\overrightarrow{\mathit{\varphi }}}_x\left(\overrightarrow{\mathit{x}};u\right)=2{\left[\frac{x_1}{{\alpha }^2+u}, \frac{x_2}{{\beta }^2+u},\frac{x_3}{{\gamma }^2+u}\right]}^T\equiv 2{{\overrightarrow{\mathit{\varphi }}}^{\prime }}_x\left(\overrightarrow{\mathit{x}};u\right),$$

$${\mathit{\varphi }}_{xx}\left(\overrightarrow{\mathit{x}};u\right)=2\left[\begin{array}{ccc}{\left({\alpha }^2+u\right)}^{-1}& 0& 0\\ 0& {\left({\beta }^2+u\right)}^{-1}& 0\\ 0& 0& {\left({\gamma }^2+u\right)}^{-1}\end{array}\right]\equiv 2{{\mathit{\varphi }}^{\prime }}_{xx}\left(\overrightarrow{\mathit{x}};u\right),$$

$$\psi \left(\overrightarrow{\mathit{x}};u^{\prime }\right)=\frac{3{\mu }_s}{\mathsf{\Delta }\left(u^{\prime }\right)}{\left\{\frac{x_1^2}{{\left({\alpha }^2+u^{\prime }\right)}^2}+\frac{x_2^2}{{\left({\beta }^2+u^{\prime }\right)}^2}+\frac{x_3^2}{{\left({\gamma }^2+u^{\prime }\right)}^2}\right\}}^{-1}=\frac{3{\mu }_s}{\mathsf{\Delta }\left(u^{\prime }\right){‖{\overrightarrow{\mathit{\varphi }}}_x^{\prime }\left(\overrightarrow{\mathit{x}};u^{\prime }\right)‖}^2},$$

$$\mathsf{\Delta }\left(u\right)=\sqrt{\left({\alpha }^2+u\right)\left({\beta }^2+u\right)\left({\gamma }^2+u\right)}.$$

3. Rigid Body Potential

‘Rigid body potential’ is defined as the potential of a small celestial body integrated across the finite volume of the S/C, with the assumption that the mass of the S/C has no influence on the motion of the small celestial body. The RBP is distinguishable from the ‘mutual potential’ in that it does not consider the motion of the small celestial body. Small celestial bodies are approximated as polyhedrons with uniform density, and the S/C is considered as an extended body with finite volume. This enables us to utilize a variety of gravity estimation models for comparison. Figure 1 visualizes the coordinate system for the RBP. In this section, the RBP is denoted as U and the subscript “x” specifically refers to the partial derivatives with respect to the variable x.

Proposition 1.

Let the potential, its gradient, and its Hessian be expressed as $V\in \mathbb{R}$, ${\overrightarrow{\mathit{V}}}_x\in {\mathbb{R}}^3$, and ${\mathit{V}}_{xx}\in {\mathbb{R}}^{3\times 3}$ for a given pair of position vector $\overrightarrow{\mathit{x}}\in {\mathbb{R}}^3$ and attitude $\mathit{R}\in SO\left(3\right)$. Then, the RBP $U\in \mathbb{R}$, its force $\overrightarrow{\mathit{f}}\in {\mathbb{R}}^3$, its Hessian matrix ${\mathit{U}}_{xx}$, and its torque $\overrightarrow{\mathit{\tau }}\in {\mathbb{R}}^3$ are generalized as follows:

$$U=VM-\frac{1}{2}{\sum }_{n=1}^3{\lambda }_n{\overrightarrow{\mathit{v}}}_n^T\mathit{R}\mathcal{J}{\mathit{R}}^T{\overrightarrow{\mathit{v}}}_n+\frac{1}{4}tr\left[\mathcal{J}\right]tr\left[{\mathit{V}}_{xx}\right],$$

(5)

$$M\overrightarrow{\mathit{f}}=M{\overrightarrow{\mathit{V}}}_x+\frac{tr\left[\mathcal{J}\right]}{2}{\sum }_{n=1}^3{\mathsf{\Lambda }}_n{\overrightarrow{\mathit{v}}}_n-{\sum }_{n=1}^3{\mathsf{\Lambda }}_n\mathit{R}\mathcal{J}{\mathit{R}}^T{\overrightarrow{\mathit{v}}}_n,$$

(6)

$${\mathit{U}}_{xx}=M{\mathit{V}}_{xx}+\frac{tr\left[\mathcal{J}\right]}{2}{\sum }_{n=1}^3\left({\mathit{u}}_{x1}^n+\frac{\partial {\overrightarrow{\mathit{v}}}_n}{\partial \overrightarrow{\mathit{x}}}{\mathsf{\Lambda }}_n^T\right)-{\sum }_{n=1}^3\left({\mathit{u}}_{x2}^n+\frac{\partial {\overrightarrow{\mathit{v}}}_n}{\partial \overrightarrow{\mathit{x}}}\mathit{R}\mathcal{J}{\mathit{R}}^T{\mathsf{\Lambda }}_n^T\right),$$

(7)

$$\mathcal{J}\overrightarrow{\mathit{\tau }}={\sum }_{n=1}^3{\lambda }_n{\left({\mathit{R}}^T{\overrightarrow{\mathit{v}}}_n\right)}^{\times }\mathcal{J}{\mathit{R}}^T{\overrightarrow{\mathit{v}}}_n.$$

(8)

${\mathit{\Lambda }}_n={\lambda }_n\left(\frac{\partial {\overrightarrow{\mathit{v}}}_n}{\partial \overrightarrow{\mathit{x}}}\right)+\frac{1}{2}\left(\frac{\partial {\lambda }_n}{\partial \overrightarrow{\mathit{x}}}\right){\overrightarrow{\mathit{v}}}_n^T$ where ${\lambda }_n\in \mathbb{R}$ and ${\overrightarrow{\mathit{v}}}_n\in {\mathbb{R}}^3$ are the n^th eigenvalue of ${\mathit{V}}_{xx}$ and its associated eigenvector, respectively. ${\mathit{u}}_{xm}^n={\left[{\overrightarrow{\mathit{u}}}_{xm}^{1n}, {\overrightarrow{\mathit{u}}}_{xm}^{2n}, {\overrightarrow{\mathit{u}}}_{xm}^{3n}\right]}^T\in {\mathbb{R}}^{3\times 3}$ for $m=1,2$. Here, ${\overrightarrow{\mathit{u}}}_{x1}^{in}=\frac{\partial {\mathsf{\Lambda }}_n}{\partial x_i}{\overrightarrow{\mathit{v}}}_n\in {\mathbb{R}}^3$ and ${\overrightarrow{\mathit{u}}}_{x2}^{in}=\frac{\partial {\mathsf{\Lambda }}_n}{\partial x_i}\mathit{R}\mathcal{J}{\mathit{R}}^T{\overrightarrow{\mathit{v}}}_n\in {\mathbb{R}}^3$ where

$$\frac{\partial {\mathsf{\Lambda }}_n}{\partial x_i}=\frac{3}{2}\frac{\partial {\lambda }_n}{\partial x_i}\frac{\partial {\overrightarrow{\mathit{v}}}_n}{\partial \overrightarrow{\mathit{x}}}+{\lambda }_n\frac{{\partial }^2{\overrightarrow{\mathit{v}}}_n}{\partial x_i\partial \overrightarrow{\mathit{x}}}+\frac{1}{2}\frac{{\partial }^2{\lambda }_n}{\partial x_i\partial \overrightarrow{\mathit{x}}}{\overrightarrow{\mathit{v}}}_n^T$$

(9)

for $i=1,2,3$. $M\in {\mathbb{R}}^{+}$ is the total mass of the S/C.

Proof of Proposition 1.

Volume integration of the S/C body $\mathcal{B}$ yields $U$ as

$$U={\int }_{\mathcal{B}}^{ }V\left(\overrightarrow{\mathit{x}}+\mathit{R}\overrightarrow{\mathit{\rho }}\right)dm\left(\overrightarrow{\mathit{\rho }}\right)$$

where $\overrightarrow{\mathit{\rho }}\in {\mathbb{R}}^3$ indicates the position of the mass element $dm$ in the body-fixed frame of the S/C as in Figure 1. It is assumed that $‖\overrightarrow{\mathit{\rho }}‖\ll ‖\overrightarrow{\mathit{x}}‖$. Let $ϵ$ and $\kappa$ be defined such that $ϵ=‖\overrightarrow{\mathit{\rho }}‖/‖\overrightarrow{\mathit{x}}‖$ and $\cos \kappa =\widehat{\mathit{x}}\cdot \mathit{R}\widehat{\mathit{\rho }}$, respectively. $‖\overrightarrow{\mathit{x}}+\mathit{R}\overrightarrow{\mathit{\rho }}‖$ is rewritten as

$$‖\overrightarrow{\mathit{x}}+\mathit{R}\overrightarrow{\mathit{\rho }}‖=‖\overrightarrow{\mathit{x}}‖{\left(1+2ϵ\cos \kappa +{ϵ}^2\right)}^{\frac{1}{2}}.$$

Then, it is possible to approximate the potential $V$ at $\overrightarrow{\mathit{x}}$ up to the second order as

$$V\left(x_{ϵ}\right)=V\left(\overrightarrow{\mathit{x}}\right)+ϵ{\left[\frac{\partial V}{\partial ϵ}\right]}_{ϵ=0}+\frac{{ϵ}^2}{2}{\left[\frac{{\partial }^{2}V}{\partial {ϵ}^2}\right]}_{ϵ=0}+\mathcal{O}\left(ϵ\right)$$

$$=V\left(\overrightarrow{\mathit{x}}\right)+{\left(\mathit{R}\overrightarrow{\mathit{\rho }}\right)}^T{\overrightarrow{\mathit{V}}}_x+\frac{1}{2}{\left(\mathit{R}\overrightarrow{\mathit{\rho }}\right)}^T{\mathit{V}}_{xx}\mathit{R}\overrightarrow{\mathit{\rho }}+\mathcal{O}\left(ϵ\right)$$

where $x_{ϵ}=‖\overrightarrow{\mathit{x}}‖{\left(1+2ϵ\cos \kappa +{ϵ}^2\right)}^{\frac{1}{2}}$ and $\mathcal{O}\left(ϵ\right)$ contain all higher-order terms. The RBP then becomes

$$U\left(\overrightarrow{\mathit{x}}\right)\approx {\int }_{\mathcal{B}}\left(V\left(\overrightarrow{\mathit{x}}\right)+{\overrightarrow{\mathit{\rho }}}^T{\mathit{R}}^T{\overrightarrow{\mathit{V}}}_x\left(\overrightarrow{\mathit{x}}\right)+\frac{1}{2}{\overrightarrow{\mathit{\rho }}}^T{\mathit{R}}^T{\mathit{V}}_{xx}\left(\overrightarrow{\mathit{x}}\right)\mathit{R}\overrightarrow{\mathit{\rho }}\right)dm\left(\overrightarrow{\mathit{\rho }}\right)$$

$$=V\left(\overrightarrow{\mathit{x}}\right)M+\frac{1}{2}{\int }_{\mathcal{B}}{\overrightarrow{\mathit{\rho }}}^T{\mathit{R}}^T{\mathit{V}}_{xx}\left(\overrightarrow{\mathit{x}}\right)\mathit{R}\overrightarrow{\mathit{\rho }}dm\left(\overrightarrow{\mathit{\rho }}\right).$$

(10)

Now, the Hessian ${\mathit{V}}_{xx}$ is decomposed by using scalar and vector pairs in order to rewrite the second term in Equation (10). As the eigenvalues of a real symmetric matrix are always real and the associated eigenvectors are orthogonal to each other [29], the Hessian ${\mathit{V}}_{xx}$ can be represented using three real eigenvalues and their associated eigenvectors. Let ${\lambda }_n$ and ${\overrightarrow{\mathit{v}}}_n$ be eigenvalues and the corresponding eigenvectors of ${\mathit{V}}_{xx}$ for $n=1, 2, 3$, respectively. Then, ${\mathit{V}}_{xx}$ is represented as

$${\mathit{V}}_{xx}={\sum }_{n=1}^3{\lambda }_n{\overrightarrow{\mathit{v}}}_n{\overrightarrow{\mathit{v}}}_n^T.$$

(11)

Note that the other pair $\left(a_n,{\overrightarrow{x}}_n\right)$ can be taken into account if ${\mathit{V}}_{xx}$ can be written as ${\sum }_{n=1}^N{a}_n{\overrightarrow{\mathit{x}}}_n{\overrightarrow{\mathit{x}}}_n^T$ for $n=1, 2,\cdots ,N$. In this study, the eigenvalues and eigenvectors are utilized to derive the RBP without loss of generality. Introducing Equation (10) into (11) yields the RBP as

$$U\left(\overrightarrow{\mathit{x}}\right)\approx V\left(\overrightarrow{\mathit{x}}\right)M+\frac{1}{2}{\sum }_{n=1}^3{\lambda }_n{\int }_{\mathcal{B}}{\overrightarrow{\mathit{\rho }}}^T{\mathit{R}}^T{\overrightarrow{\mathit{v}}}_n{\overrightarrow{\mathit{v}}}_n^T\mathit{R}\overrightarrow{\mathit{\rho }}dm\left(\overrightarrow{\mathit{\rho }}\right)$$

$$=V\left(\overrightarrow{\mathit{x}}\right)M+\frac{1}{2}{\sum }_{n=1}^3{\lambda }_n{\overrightarrow{\mathit{v}}}_n^T\mathit{R}{\int }_{\mathcal{B}}\overrightarrow{\mathit{\rho }}{\overrightarrow{\mathit{\rho }}}^{T}dm\left(\overrightarrow{\mathit{\rho }}\right){\mathit{R}}^T{\overrightarrow{\mathit{v}}}_n.$$

Since ${\int }_{\mathcal{B}}\overrightarrow{\mathit{\rho }}{\overrightarrow{\mathit{\rho }}}^{T}dm\left(\overrightarrow{\mathit{\rho }}\right)=\frac{1}{2}tr\left[\mathcal{J}\right]{\mathbb{I}}_3-\mathcal{J}$, the RBP is finally expressed as

$$U\left(\overrightarrow{\mathit{x}}\right)=V\left(\overrightarrow{\mathit{x}}\right)M-\frac{1}{2}{\sum }_{n=1}^3{\lambda }_n{\overrightarrow{\mathit{v}}}_n^T\mathit{R}\mathcal{J}{\mathit{R}}^T{\overrightarrow{\mathit{v}}}_n+\frac{1}{4}tr\left[\mathcal{J}\right]tr\left[{\mathit{V}}_{xx}\right].$$

(12)

The gravitational force $\overrightarrow{\mathit{f}}$ is obtained by taking partial derivative of Equation (12):

$$M\overrightarrow{\mathit{f}}=M{\overrightarrow{\mathit{V}}}_x+\frac{tr\left[\mathcal{J}\right]}{2}{\sum }_{n=1}^3\left\{{\lambda }_n\left(\frac{\partial {\overrightarrow{\mathit{v}}}_n}{\partial \overrightarrow{\mathit{x}}}\right)+\frac{1}{2}\left(\frac{\partial {\lambda }_n}{\partial \overrightarrow{\mathit{x}}}\right){\overrightarrow{\mathit{v}}}_n^T\right\}{\overrightarrow{\mathit{v}}}_n$$

$$-{\sum }_{n=1}^3\left\{{\lambda }_n\left(\frac{\partial {\overrightarrow{\mathit{v}}}_n}{\partial \overrightarrow{\mathit{x}}}\right)+\frac{1}{2}\left(\frac{\partial {\lambda }_n}{\partial \overrightarrow{\mathit{x}}}\right){\overrightarrow{\mathit{v}}}_n^T\right\}\mathit{R}\mathcal{J}{\mathit{R}}^T{\overrightarrow{\mathit{v}}}_n$$

$$=M{\overrightarrow{\mathit{V}}}_x+\frac{tr\left[\mathcal{J}\right]}{2}{\sum }_{n=1}^3{\mathsf{\Lambda }}_n{\overrightarrow{\mathit{v}}}_n-{\sum }_{n=1}^3{\mathsf{\Lambda }}_n\mathit{R}\mathcal{J}{\mathit{R}}^T{\overrightarrow{\mathit{v}}}_n.$$

The gravitational torque $\overrightarrow{\mathit{\tau }}$ can be derived by taking the partial derivative of Equation (12) with respect to $\mathit{R}$:

$$\frac{\partial U}{\partial \mathit{R}}=\frac{\partial }{\partial \mathit{R}}\left\{VM+\frac{tr\left[\mathcal{J}\right]}{4}{\sum }_{n=1}^3{\lambda }_n{\overrightarrow{\mathit{v}}}_n^T{\overrightarrow{\mathit{v}}}_n-\frac{1}{2}{\sum }_{n=1}^3{\lambda }_n{\overrightarrow{\mathit{v}}}_n^T\mathit{R}\mathcal{J}{\mathit{R}}^T{\overrightarrow{\mathit{v}}}_n\right\}$$

$$=-{\sum }_{n=1}^3{\lambda }_n{\overrightarrow{\mathit{v}}}_n{\overrightarrow{\mathit{v}}}_n^T\mathit{R}\mathcal{J}=-V_{xx}\mathit{R}\mathcal{J}.$$

(13)

Substituting Equation (13) into ${\left(\mathcal{J}\overrightarrow{\mathit{\tau }}\right)}^{\times }={\left(\frac{\partial U}{\partial \mathit{R}}\right)}^T\mathit{R}-{\mathit{R}}^T\frac{\partial U}{\partial \mathit{R}}$ provides

$${\left(\mathcal{J}\overrightarrow{\mathit{\tau }}\right)}^{\times }={\sum }_{n=1}^3{\lambda }_n\left\{-\mathcal{J}{\mathit{R}}^T{\overrightarrow{\mathit{v}}}_n{\overrightarrow{\mathit{v}}}_n^T\mathit{R}+{\mathit{R}}^T{\overrightarrow{\mathit{v}}}_n{\overrightarrow{\mathit{v}}}_n^T\mathit{R}\mathcal{J}\right\}.$$

Using $\mathit{y}{\mathit{x}}^T-\mathit{x}{\mathit{y}}^T={\left({\mathit{x}}^{\times }\mathit{y}\right)}^{\times }$ for $\mathit{x}, \mathit{y}\in {\mathbb{R}}^3$, the torque due to the RBP can be given as Equation (8).

In order to derive the Hessian matrix ${\mathit{U}}_{xx}$, start by taking partial derivatives of ${\overrightarrow{\mathit{U}}}_x$ with respect to $x_i$ as

$$\frac{\partial }{\partial x_i}\left(\frac{\partial U}{\partial \overrightarrow{\mathit{x}}}\right)=\frac{\partial }{\partial x_i}\left(M{\overrightarrow{\mathit{V}}}_x+\frac{tr\left[\mathcal{J}\right]}{2}{\sum }_{n=1}^3{\mathsf{\Lambda }}_n{\overrightarrow{\mathit{v}}}_n-{\sum }_{n=1}^3{\mathsf{\Lambda }}_n\mathit{R}\mathcal{J}{\mathit{R}}^T{\overrightarrow{\mathit{v}}}_n\right)$$

$$=M\frac{\partial {\overrightarrow{\mathit{V}}}_x}{\partial x_i}+\frac{tr\left[\mathcal{J}\right]}{2}\sum _{n=1}^3\left({\overrightarrow{\mathit{u}}}_{x1}^{in}+{\mathsf{\Lambda }}_n\frac{\partial {\overrightarrow{\mathit{v}}}_n}{\partial x_i}\right)-\sum _{n=1}^3\left({\overrightarrow{\mathit{u}}}_{x2}^{in}+{\mathsf{\Lambda }}_n\mathit{R}\mathcal{J}{\mathit{R}}^T\frac{\partial {\overrightarrow{\mathit{v}}}_n}{\partial x_i}\right)$$

(14)

where ${\overrightarrow{\mathit{u}}}_{x1}^{in}=\frac{\partial {\mathsf{\Lambda }}_n}{\partial x_i}{\overrightarrow{\mathit{v}}}_n$ and ${\overrightarrow{\mathit{u}}}_{x2}^{in}=\frac{\partial {\mathsf{\Lambda }}_n}{\partial x_i}\mathit{R}\mathcal{J}{\mathit{R}}^T{\overrightarrow{\mathit{v}}}_n$. The first term of Equation (14) leads to

$$\frac{\partial {\overrightarrow{\mathit{V}}}_x}{\partial \overrightarrow{\mathit{x}}}={\left[\frac{\partial {\overrightarrow{\mathit{V}}}_x}{\partial x_1},\frac{\partial {\overrightarrow{\mathit{V}}}_x}{\partial x_2},\frac{\partial {\overrightarrow{\mathit{V}}}_x}{\partial x_3}\right]}^T={\mathit{V}}_{xx}.$$

(15)

The common term $\frac{\partial {\mathsf{\Lambda }}_n}{\partial x_i}$ of ${\overrightarrow{\mathit{u}}}_{x1}^{in}$ and ${\overrightarrow{\mathit{u}}}_{x2}^{in}$ yields

$$\frac{\partial {\mathsf{\Lambda }}_n}{\partial x_i}=\frac{3}{2}\frac{\partial {\lambda }_n}{\partial x_i}\frac{\partial {\overrightarrow{\mathit{v}}}_n}{\partial \overrightarrow{\mathit{x}}}+{\lambda }_n\frac{{\partial }^2{\overrightarrow{\mathit{v}}}_n}{\partial x_i\partial \overrightarrow{\mathit{x}}}+\frac{1}{2}\frac{{\partial }^2{\lambda }_n}{\partial x_i\partial \overrightarrow{\mathit{x}}}{\overrightarrow{\mathit{v}}}_n^T.$$

The second term of Equation (14) provides

$$\frac{\partial }{\partial \overrightarrow{\mathit{x}}}\left({\sum }_{n=1}^3{\mathsf{\Lambda }}_n{\overrightarrow{\mathit{v}}}_n\right)={\left[{\sum }_{n=1}^3\left({\overrightarrow{\mathit{u}}}_{x1}^{1n}+{\mathsf{\Lambda }}_n\frac{\partial {\overrightarrow{\mathit{v}}}_n}{\partial x_1}\right),{\sum }_{n=1}^3\left({\overrightarrow{\mathit{u}}}_{x1}^{2n}+{\mathsf{\Lambda }}_n\frac{\partial {\overrightarrow{\mathit{v}}}_n}{\partial x_2}\right),{\sum }_{n=1}^3\left({\overrightarrow{\mathit{u}}}_{x1}^{3n}+{\mathsf{\Lambda }}_n\frac{\partial {\overrightarrow{\mathit{v}}}_n}{\partial x_3}\right)\right]}^T$$

$$={{\sum }_{n=1}^3\left\{\left[{\overrightarrow{\mathit{u}}}_{x1}^{1n},{\overrightarrow{\mathit{u}}}_{x1}^{2n},{\overrightarrow{\mathit{u}}}_{x1}^{3n}\right]+{\mathsf{\Lambda }}_n\left[\frac{\partial {\overrightarrow{\mathit{v}}}_n}{\partial x_1},\frac{\partial {\overrightarrow{\mathit{v}}}_n}{\partial x_2},\frac{\partial {\overrightarrow{\mathit{v}}}_n}{\partial x_3}\right]\right\}}^T$$

$$={\sum }_{n=1}^3\left({\mathit{u}}_{x1}^n+\frac{\partial {\overrightarrow{\mathit{v}}}_n}{\partial \overrightarrow{\mathit{x}}}{\mathsf{\Lambda }}_n^T\right).$$

(16)

Similarly, the third term of Equation (14) becomes

$$\frac{\partial }{\partial \overrightarrow{\mathit{x}}}\left\{{\sum }_{n=1}^3\left({\overrightarrow{\mathit{u}}}_{x2}^{in}+{\mathsf{\Lambda }}_n\mathit{R}\mathcal{J}{\mathit{R}}^T\frac{\partial {\overrightarrow{\mathit{v}}}_n}{\partial x_i}\right)\right\}={\sum }_{n=1}^3\left({\mathit{u}}_{x2}^n+\frac{\partial {\overrightarrow{\mathit{v}}}_n}{\partial \overrightarrow{\mathit{x}}}\mathit{R}\mathcal{J}{\mathit{R}}^T{\mathsf{\Lambda }}_n^T\right).$$

(17)

Finally, the Hessian of the RBP can be obtained by combining Equations (15)–(17). □

Proposition 1 can be extended to cases where the Hessian matrix of potential ${\mathit{V}}_{xx}$ is decomposed into two distinct vectors. Refer to Appendix A for details.

4. Application of Rigid Body Potential

Section 4 presents the RBP estimated through three gravity estimation models, namely, the PM, EB, and TE models introduced in Section 3. As stipulated by Proposition 1, the derivation of the $U$ and its corresponding torque necessitates both the potential $V$ itself and its associated Hessian ${\mathit{V}}_{xx}$. Furthermore, the calculation of the force stemming from the RBP mandates the partial derivatives of certain scalar quantities and vectors linked to ${\mathit{V}}_{xx}$ and requires ${\overrightarrow{\mathit{V}}}_x$. In this section, the superscript “A” indicates its computation through the utilization of the A gravity estimation model.

4.1. Point Mass Model

This section presents the RBP, along with its gradient ${\overrightarrow{\mathit{U}}}_x^{PM}$ and its Hessian ${\mathit{U}}_{xx}^{PM}$ when a small celestial body is conceptualized as a PM. Although it is challenging to assert that the PM model is accurate enough, its usefulness lies in understanding how Proposition 1 can be applied to a variety of gravity estimation models.

Example 1.

The RBP $U^{PM}\in \mathbb{R}$ with the PM model, along with its associated force ${\overrightarrow{\mathit{f}}}^{PM}\in {\mathbb{R}}^3$, Hessian ${\mathit{U}}_{xx}^{PM}\in {\mathbb{R}}^{3\times 3}$, and torque ${\overrightarrow{\mathit{\tau }}}^{PM}\in {\mathbb{R}}^3$, is represented as

$$U^{PM}=\frac{{\mu }_s}{‖\overrightarrow{\mathit{x}}‖}\left\{M-\frac{3}{2}\frac{{\widehat{\mathit{x}}}^T\mathit{R}\mathcal{J}{\mathit{R}}^T\widehat{\mathit{x}}}{{‖\overrightarrow{\mathit{x}}‖}^2}+\frac{1}{2}\frac{tr\left[\mathcal{J}\right]}{{‖\overrightarrow{\mathit{x}}‖}^2}\right\},$$

(18)

$$M{\overrightarrow{\mathit{f}}}^{PM}=-\frac{{\mu }_s}{{‖\overrightarrow{\mathit{x}}‖}^3}\left\{M{\mathbb{I}}_3-\frac{15}{2}\frac{{\widehat{\mathit{x}}}^{\mathit{T}}\mathit{R}\mathcal{J}{\mathit{R}}^{\mathit{T}}\widehat{\mathit{x}}}{{‖\overrightarrow{\mathit{x}}‖}^2}{\mathbb{I}}_3+3\frac{\mathit{R}\mathcal{J}{\mathit{R}}^{\mathit{T}}}{{‖\overrightarrow{\mathit{x}}‖}^2}+\frac{3}{2}\frac{tr\left[\mathcal{J}\right]}{{‖\overrightarrow{\mathit{x}}‖}^2}{\mathbb{I}}_3\right\}\overrightarrow{\mathit{x}}={\overrightarrow{\mathit{U}}}_x^{PM},$$

$${\mathit{U}}_{xx}^{PM}=\frac{{\mu }_s}{{‖\overrightarrow{\mathit{x}}‖}^3}\left\{M\left(3\widehat{\mathit{x}}{\widehat{\mathit{x}}}^T-{\mathbb{I}}_3\right)-\frac{15}{2}\frac{{\widehat{\mathit{x}}}^T\mathit{R}\mathcal{J}{\mathit{R}}^{\mathit{T}}\widehat{\mathit{x}}}{{‖\overrightarrow{\mathit{x}}‖}^2}\left(7\widehat{\mathit{x}}{\widehat{\mathit{x}}}^T-{\mathbb{I}}_3\right)+3\frac{\mathit{R}\mathcal{J}{\mathit{R}}^{\mathit{T}}}{{‖\overrightarrow{\mathit{x}}‖}^2}\left(10\widehat{\mathit{x}}{\widehat{\mathit{x}}}^T-{\mathbb{I}}_3\right)+\frac{3}{2}\frac{tr\left[\mathcal{J}\right]}{{‖\overrightarrow{\mathit{x}}‖}^2}\left(5\widehat{\mathit{x}}{\widehat{\mathit{x}}}^T-{\mathbb{I}}_3\right)\right\},$$

(19)

$$\mathcal{J}{\overrightarrow{\mathit{\tau }}}^{PM}=3\frac{{\mu }_s}{{‖\overrightarrow{\mathit{x}}‖}^3}{\left({\mathit{R}}^T\widehat{\mathit{x}}\right)}^{\times }\mathcal{J}{\mathit{R}}^T\widehat{\mathit{x}}.$$

Proof of Example 1.

The initial step involves identifying scalar–vector pairs that facilitate the representation of ${\mathit{V}}_{xx}^{PM},$ as demonstrated in Equation (11). As evident from Equation (2), the first term of ${\mathit{V}}_{xx}^{PM}$ already resembles the desired form of Equation (11), while the other term is a diagonal matrix. Any diagonal matrix can be transformed into the desired form by employing the standard basis and considering its diagonal elements as

$${\mathit{V}}_{xx}^{PM}=\sum _{n=1}^4{\lambda }_n^{PM}{\overrightarrow{\mathit{v}}}_n^{PM}{{\overrightarrow{\mathit{v}}}_n^{PM}}^T$$

(20)

where ${\overrightarrow{\mathit{v}}}_n^{PM}$, for $n=1, 2, 3$, represents the n^th standard basis vector ${\widehat{\mathit{e}}}_n$, and ${\lambda }_n^{PM}$ corresponds to its associated scalar $-\frac{{\mu }_S}{{‖\overrightarrow{\mathit{x}}‖}^3}$. The scalar–vector pair that remains is $\left({\lambda }_4, {\overrightarrow{\mathit{v}}}_4\right)=\left(3\frac{{\mu }_S}{{‖\overrightarrow{\mathit{x}}‖}^5},\overrightarrow{\mathit{x}}\right)$. The remaining proof involves the application of Corollary A1 from Appendix A with Equation (20), and the details are provided in Appendix B.1. □

4.2. Extended Body Model

This section presents the RBP and its associated parameters by employing the EB model where both the S/C and small celestial body are considered as extended bodies with finite volumes.

Example 2.

Let $M_s\in {\mathbb{R}}^{+}$ and ${\mathcal{J}}_s\in {\mathbb{R}}^{3\times 3}$ be the mass and the MOI of a small celestial body, respectively. For a given $V^{EB}\in \mathbb{R}$, accompanied by its gradient ${\overrightarrow{\mathit{V}}}_x^{EB}\in {\mathbb{R}}^3$, its Hessian ${\mathit{V}}_{xx}^{EB}\in {\mathbb{R}}^{3\times 3}$, the RBP $U^{EB}\in \mathbb{R}$ with the EB model, along with its corresponding ${\overrightarrow{\mathit{f}}}^{EB}\in {\mathbb{R}}^3$, ${\mathit{U}}_{xx}^{EB}\in {\mathbb{R}}^{3\times 3},$ and ${\overrightarrow{\mathit{\tau }}}^{EB}\in {\mathbb{R}}^3$, can be expressed as

$$U^{EB}=V^{EB}M-\frac{1}{2}{\sum }_{n=1}^4{{\overrightarrow{\mathit{v}}}_{2n}^{EB}}^T\mathit{R}\mathcal{J}{\mathit{R}}^T{\overrightarrow{\mathit{v}}}_{1n}^{EB}+\frac{1}{4}tr\left[\mathcal{J}\right]tr\left[{\mathit{V}}_{xx}^{EB}\right],$$

$$M{\overrightarrow{\mathit{f}}}^{EB}=M{\overrightarrow{\mathit{V}}}_x^{EB}-\frac{1}{4}{\sum }_{\begin{array}{c}i,j\in \left\{1,2\right\}\\ i\ne j\end{array}}{\sum }_{n=1}^4\left(\frac{\partial {\overrightarrow{\mathit{v}}}_{in}^{EB}}{\partial \overrightarrow{\mathit{x}}}\right)\left(2\mathcal{J}{\mathit{R}}^T+tr\left[\mathcal{J}\right]{\mathbb{I}}_3\right){\overrightarrow{\mathit{v}}}_{jn}^{EB}={\overrightarrow{\mathit{U}}}_x^{EB},$$

$${\mathit{U}}_{xx}^{EB}=M{\mathit{V}}_{xx}^{EB}-{\mathit{u}}_x^4-{\mathit{u}}_x^{123}-{\mathit{H}}^{EB},$$

(21)

$$\mathcal{J}{\overrightarrow{\mathit{\tau }}}^{EB}=\sum _{n=1}^4{\left({\mathit{R}}^T{\overrightarrow{\mathit{v}}}_{1n}^{EB}\right)}^{\times }\mathcal{J}{\mathit{R}}^T{\overrightarrow{\mathit{v}}}_{2n}^{EB}$$

where

$${\overrightarrow{\mathit{v}}}_{1n}^{EB}=\left\{\begin{array}{cc}3\frac{{\mu }_s}{{‖\overrightarrow{\mathit{x}}‖}^7}{B}_{nn}^{EB}{\widehat{\mathit{e}}}_n,& n=1,2,3\\ 3\frac{{\mu }_s}{{‖\overrightarrow{\mathit{x}}‖}^9}{\mathit{A}}^{EB}\overrightarrow{\mathit{x}},& n=4\end{array}\right.,$$

(22)

$${\overrightarrow{\mathit{v}}}_{2n}^{EB}=\left\{\begin{array}{cc}{\widehat{\mathit{e}}}_n,& n=1,2,3\\ \overrightarrow{\mathit{x}},& n=4\end{array}\right.,$$

(23)

$$\frac{\partial {\overrightarrow{\mathit{v}}}_{1n}^{EB}}{\partial \overrightarrow{\mathit{x}}}=\left\{\begin{array}{cc}-21\frac{{\mu }_s}{{‖\overrightarrow{\mathit{x}}‖}^9}\overrightarrow{\mathit{x}}{\widehat{\mathit{e}}}_n^T{B}_{nn}^{EB}+3\frac{{\mu }_s}{{‖\overrightarrow{\mathit{x}}‖}^7}{\mathit{D}}^{EB}\overrightarrow{\mathit{x}}{\widehat{\mathit{e}}}_n^T,& n=1,2,3\\ 3\frac{{\mu }_s}{{‖\overrightarrow{\mathit{x}}‖}^9}{\mathit{A}}^{EB}+3\frac{{\mu }_s}{{‖\overrightarrow{\mathit{x}}‖}^9}{\mathit{C}}^{EB}\overrightarrow{\mathit{x}}{\overrightarrow{\mathit{x}}}^T,& n=4\end{array}\right.,$$

(24)

$$\frac{\partial {\overrightarrow{\mathit{v}}}_{2n}^{EB}}{\partial \overrightarrow{\mathit{x}}}=\left\{\begin{array}{cc}{\mathbb{O}}_3,& n=1,2,3\\ {\mathbb{I}}_3,& n=4\end{array} \right.,$$

(25)

$${\mathit{u}}_x^4=\frac{1}{4}\left({\mathit{E}}^{EB}\overrightarrow{\mathit{x}}{\overrightarrow{\mathit{x}}}^T+\overrightarrow{\mathit{x}}{\overrightarrow{\mathit{x}}}^T{\mathit{E}}^{EB}\right),$$

$${\mathit{u}}_x^{123}=\frac{1}{4}\left({\mathit{F}}^{EB}\overrightarrow{\mathit{x}}{\overrightarrow{\mathit{x}}}^T+\overrightarrow{\mathit{x}}{\overrightarrow{\mathit{x}}}^T{\mathit{F}}^{EB}+{\mathit{G}}^{EB}\right),$$

$${\mathit{A}}^{EB}={‖\overrightarrow{\mathit{x}}‖}^2{\mathbb{I}}_3+\frac{5}{2}\frac{G}{M_s}\left(tr\left[{\mathcal{J}}_s\right]{\mathbb{I}}_3+4{\mathcal{J}}_s-7\frac{{\overrightarrow{\mathit{x}}}^T{\mathcal{J}}_s\overrightarrow{\mathit{x}}}{{‖\overrightarrow{\mathit{x}}‖}^2}{\mathbb{I}}_3\right),$$

$${\mathit{B}}^{EB}=\frac{3}{4}\frac{G}{M_s}\left(tr\left[{\mathcal{J}}_s\right]{\mathbb{I}}_3+2{\mathcal{J}}_s-5\frac{{\overrightarrow{\mathit{x}}}^T{\mathcal{J}}_s\overrightarrow{\mathit{x}}}{{‖\overrightarrow{\mathit{x}}‖}^2}{\mathbb{I}}_3\right)-\frac{1}{2}{‖\overrightarrow{\mathit{x}}‖}^2{\mathbb{I}}_3,$$

$${\mathit{C}}^{EB}=2{\mathbb{I}}_3+35\frac{G}{M_s}\frac{{\overrightarrow{\mathit{x}}}^T{\mathcal{J}}_s\overrightarrow{\mathit{x}}}{{‖\overrightarrow{\mathit{x}}‖}^4}{\mathbb{I}}_3-35\frac{G}{M_s}\frac{{\mathcal{J}}_s}{{‖\overrightarrow{\mathit{x}}‖}^2}-\frac{9}{{‖\overrightarrow{\mathit{x}}‖}^2}{\mathit{A}}^{EB},$$

$${\mathit{D}}^{EB}=-\frac{1}{2}\left(15\frac{G}{M_s}\frac{{\mathcal{J}}_s}{{‖\overrightarrow{\mathit{x}}‖}^2}-15\frac{G}{M_s}\frac{{\overrightarrow{\mathit{x}}}^T{\mathcal{J}}_s\overrightarrow{\mathit{x}}}{{‖\overrightarrow{\mathit{x}}‖}^4}{\mathbb{I}}_3+{\mathbb{I}}_3\right),$$

$${\mathit{E}}^{EB}=\frac{3}{2}\frac{{\mu }_s}{{‖\overrightarrow{\mathit{x}}‖}^9}\left\{8{\mathbb{I}}_3-10{\mathit{C}}^{EB}+{‖\overrightarrow{\mathit{x}}‖}^2{\mathit{C}}^{EB}+9{\mathit{A}}^{EB}-70\frac{G}{M_s}\frac{{\mathcal{J}}_s}{{‖\overrightarrow{\mathit{x}}‖}^2}\right\},$$

$${\mathit{F}}^{EB}=\frac{9}{2}\frac{{\mu }_s}{{‖\overrightarrow{\mathit{x}}‖}^{10}}\left\{\frac{21}{{‖\overrightarrow{\mathit{x}}‖}^2}tr\left[{\mathit{B}}^{EB}\right]{\mathbb{I}}_3-14{\mathit{D}}^{EB}+\frac{15}{2}\frac{G}{M_s}\frac{1}{{‖\overrightarrow{\mathit{x}}‖}^2}\left(2{\mathcal{J}}_s-{\widehat{\mathit{x}}}^T{\mathcal{J}}_s\widehat{\mathit{x}}{\mathbb{I}}_3\right)\right\},$$

$${\mathit{G}}^{EB}=9\frac{{\mu }_s}{{‖\overrightarrow{\mathit{x}}‖}^7}\left({\mathit{D}}^{EB}-\frac{7}{{‖\overrightarrow{\mathit{x}}‖}^2}tr\left[{\mathit{B}}^{EB}\right]{\mathbb{I}}_3\right),$$

$${\mathit{H}}^{EB}=\frac{1}{4}\left\{\left(\frac{\partial {\overrightarrow{\mathit{v}}}_{14}^{EB}}{\partial \overrightarrow{\mathit{x}}}\right)\left(2\mathcal{J}{\mathit{R}}^T+tr\left[\mathcal{J}\right]{\mathbb{I}}_3\right)+\left(2\mathcal{J}{\mathit{R}}^T+tr\left[\mathcal{J}\right]{\mathbb{I}}_3\right){\left(\frac{\partial {\overrightarrow{\mathit{v}}}_{14}^{EB}}{\partial \overrightarrow{\mathit{x}}}\right)}^T\right\}.$$

Proof of Example 2.

Refer to Appendix B.2. □

4.3. Triaxial Ellipsoidal Model

This section presents the RBP $U^{TE}$ obtained using the TE model, together with its associated force ${\overrightarrow{\mathit{f}}}^{TE}$, torque ${\overrightarrow{\mathit{\tau }}}^{TE}$, and Hessian ${\mathit{U}}_{xx}^{TE}$.

Example 3.

Let the potential with the TE model, its corresponding gradient, and its Hessian matrix be denoted as $V^{TE}\in \mathbb{R}$,  ${\overrightarrow{\mathit{V}}}_x^{TE}\in {\mathbb{R}}^3$, and  ${\mathit{V}}_{xx}^{TE}\in {\mathbb{R}}^{3\times 3}$, respectively. Then, the RBP  $U^{TE}\in \mathbb{R}$, along with  ${\overrightarrow{\mathit{f}}}^{TE}\in {\mathbb{R}}^3$,  ${\mathit{U}}_{xx}^{EB}\in {\mathbb{R}}^{3\times 3}$, and  ${\overrightarrow{\mathit{\tau }}}^{TE}\in {\mathbb{R}}^3$, can be expressed as

$$U^{TE}=M{V}^{TE}-\frac{3\mu }{4}{\int }_{u^{\prime }}^{\infty }\left[tr\left[\mathcal{J}{\mathit{\varphi }}_{xx}^{\prime }\right]-\frac{1}{4}tr\left[\mathcal{J}\right]tr\left[{\mathit{\varphi }}_{xx}^{\prime }\right]\right]\frac{du}{\mathsf{\Delta }\left(u\right)}+\psi {{\overrightarrow{\mathit{\varphi }}}_x^{\prime }}^T\left({\mathbb{I}}_3-\frac{1}{2}\mathit{R}\mathcal{J}{\mathit{R}}^T\right){\overrightarrow{\mathit{\varphi }}}_x^{\prime },$$

$$M{\overrightarrow{\mathit{f}}}^{TE}=M{\overrightarrow{\mathit{V}}}_x^{TE}-\frac{\psi }{4}tr\left[\mathcal{J}\right]\left({\mathit{\varphi }}_{xx}^{\prime }+tr\left[{\mathit{\varphi }}_{xx}^{\prime }\right]{\mathbb{I}}_3\right){\overrightarrow{\mathit{\varphi }}}_x^{\prime }+{\mathit{A}}^{TE}={\overrightarrow{\mathit{U}}}_x^{TE},$$

(26)

$${\mathit{U}}_{xx}^{TE}=M{\mathit{V}}_{xx}^{TE}+\left(\frac{tr\left[\mathcal{J}\right]}{2}-\mathit{R}\mathcal{J}{\mathit{R}}^T\right){\mathit{B}}^{TE}+\frac{tr\left[\mathcal{J}\right]}{2}{\mathit{u}}_{x1}^{TE}-{\mathit{u}}_{x2}^{TE},$$

(27)

$$\mathcal{J}{\overrightarrow{\mathit{\tau }}}^{TE}=\psi {\left({\mathit{R}}^T{\overrightarrow{\mathit{\varphi }}}_x^{\prime }\right)}^{\times }\mathcal{J}{\mathit{R}}^T{{\overrightarrow{\mathit{\varphi }}}^{\prime }}_x$$

where

$${\mathit{A}}^{TE}=\frac{\psi }{2}\left\{2{\mathit{\varphi }}_{xx}^{\prime }\left({\mathbb{I}}_3-{{\widehat{\mathit{\varphi }}}^{\prime }}_x{{{\widehat{\mathit{\varphi }}}^{\prime }}_x}^T\right)-tr\left[{{\mathit{\varphi }}^{\prime }}_{xx}\right]{{\widehat{\mathit{\varphi }}}^{\prime }}_x{{{\widehat{\mathit{\varphi }}}^{\prime }}_x}^T\right\}\mathit{R}\mathcal{J}{\mathit{R}}^T{{\overrightarrow{\mathit{\varphi }}}^{\prime }}_x,$$

$${\mathit{B}}^{TE}={\mathit{\varphi }}_{xx}^{\prime }\left\{\psi \left({\mathbb{I}}_3-{{\widehat{\mathit{\varphi }}}^{\prime }}_x{{{\widehat{\mathit{\varphi }}}^{\prime }}_x}^T\right){\mathit{\varphi }}_{xx}^{\prime }-\frac{1}{2}\psi tr\left[{{\mathit{\varphi }}^{\prime }}_{xx}\right]{{\widehat{\mathit{\varphi }}}^{\prime }}_x{{{\widehat{\mathit{\varphi }}}^{\prime }}_x}^T\right\},$$

$${\mathit{u}}_{x1}^{TE}=\frac{1}{2}\psi \left(tr\left[{{\mathit{\varphi }}^{\prime }}_{xx}\right]{\mathbb{I}}_3+2{{\mathit{\varphi }}^{\prime }}_{xx}\right)\left(tr\left[{{\mathit{\varphi }}^{\prime }}_{xx}\right]{\mathbb{I}}_3+3{{\mathit{\varphi }}^{\prime }}_{xx}\right){\widehat{\mathit{\varphi }}\prime }_x{\widehat{\mathit{\varphi }}\prime }_x^T-\frac{1}{2}\psi {{\mathit{\varphi }}^{\prime }}_{xx}^2,$$

$${\mathit{u}}_{x2}^{TE}=\frac{1}{2}\psi \left(tr\left[{\mathit{\varphi }\prime }_{xx}\right]{\mathbb{I}}_3+2{\mathit{\varphi }\prime }_{xx}\right)\left(tr\left[{{\mathit{\varphi }}^{\prime }}_{xx}\right]\mathit{R}\mathcal{J}{\mathit{R}}^T+tr\left[\mathcal{J}\right]{\mathit{\varphi }}_{xx}^{\prime }\right){\widehat{\mathit{\varphi }}\prime }_x{\widehat{\mathit{\varphi }}\prime }_x^T-\frac{1}{2}\psi tr\left[\mathcal{J}\right]{{\mathit{\varphi }}^{\prime }}_{xx}^2.$$

Proof of Example 3.

Refer to Appendix B.3. □

5. Numerical Simulations

Section 5 presents numerical analyses of the RBP $V$ and its associated Hessian $U$ with three aforementioned gravity estimation models. The direct integration model is chosen as a reference for comparison, since it is the most accurate for a given polyhedral body [30,31]. Refer to Appendix C for details on the direct integration model. Eight small celestial bodies, i.e., 4 Vesta, 243 Ida, 951 Gaspra, 1036 Ganymed, 2063 Bacchus, 4769 Castalia, 25143 Itokawa, and 99942 Apophis, have been chosen for evaluating/testing the gravity estimation models. Their masses range from ${10}^{10} \mathrm{k}\mathrm{g}$ to ${10}^{20} \mathrm{k}\mathrm{g}$, while the number of faces is distributed across a spectrum from 2024 to 32040. The outcomes related to three pairs of small bodies characterized by similar masses, i.e., (253 Ida, 1036 Ganymed), (2163 Bacchus, 4769 Castalia), and (26,143 Itokawa, 99942 Apophis), serve to demonstrate the impact of both the number of faces and the shape of these celestial bodies.

Table 2. Mass and 1 DU values for small celestial bodies.

Small Body	$\mathbf{Mass} (10^{15} \mathbf{k}\mathbf{g}$)	$1 \mathbf{DU} (\mathbf{k}\mathbf{m}$)	Number of Faces
4 Vesta	$2.59\times {10}^5$	$287.7538$	$5040$
243 Ida	$100$	$27.5884$	$\mathrm{32,040}$
951 Gaspra	$3.57$	$10.7637$	$\mathrm{32,040}$
1036 Ganymed	$167$	$34.28$	$2040$
2063 Bacchus	$0.0033$	$0.6617$	$4092$
4769 Castalia	$0.0014$	$0.9849$	$4092$
25143 Itokawa	$3.5\times {10}^{-5}$	$0.4373$	$\mathrm{12,192}$
99942 Apophis	$2.70\times {10}^{-5}$	$0.9957$	$2024$

shows the physical/modelling properties of eight small celestial bodies. $1R_S$ is selected as 1 distance unit (DU). The MOI and mass of the S/C are selected as $diag\left(\left[200,500,300\right]\right)\mathrm{k}\mathrm{g}\cdot {\mathrm{m}}^2$ and $600 \mathrm{k}\mathrm{g}$, respectively.

5.1. Analysis on Gravity Estimation Models

The aforementioned gravity estimation models are numerically analyzed first, followed by analysis of the RBP. Each gravity estimation model is assessed based on criteria that include potential accuracy. Accuracy evaluations are performed with respect to the radial distance along the x-, y-, and z-axes. Figure 2 shows the relative errors of the potential attributed to each gravity estimation model when applied to the eight small celestial bodies. The left and right illustrate the relative errors of the potential at $1DU$ and at $3DU$, respectively. The outcomes pertaining to the x-, y-, and z-axes are denoted as circle, square, and triangle markers, respectively. Blue, red, and yellow markers indicate results calculated using the PM, EB, and TE models, respectively. In the left figure, the relative errors tend to be larger compared to those in the right figure, and the accuracy order for each small celestial body is usually consistent regardless of the radial distance. The TE model generally yields the smallest errors in most cases, except for the results at $3DU$ for 99942 Apophis. The PM model gives relatively accurate results when applied to spherical bodies (i.e., 4 Vesta and 1036 Ganymed), similar to the EB model. Although the results are ordered based on the mass of small celestial bodies, it is unclear to identify any clear tendency with respect to mass. The relative errors associated with heavier bodies decrease at a slower rate compared to those of lighter bodies, with respect to radial distance, if $1DU$ is considered to be a significant magnitude for heavier bodies. These results will be further discussed in terms of the conditions suitable for the utilization of RBP, in conjunction with the analysis of RBP. Figure 3 presents the relative errors, mean, and standard deviation for each gravity estimation models in a single plot, without distinguishing between the small celestial bodies. Blue, red, and yellow markers denote the relative errors at $1DU$, $2DU$, and $3DU$, respectively. The relative errors of all gravity estimation models decrease in relation to radial distance, regardless of the small celestial bodies and evaluation directions. The TE model exhibits the best performance among the three gravity estimation models, as evidenced by its smallest mean and standard deviation. The observed trend, where the PM model exhibits the largest relative errors, suggests that its accuracy is mainly reliant on the shape of the small celestial bodies. This is because the PM model considers a small celestial body as a sphere, regardless of its actual shape. Thus, the PM model provides the most accurate results when a sphere with uniform density is employed as a small celestial body.

5.2. Analysis on Rigid Body Potential

The RBP is composed of the potential $V$ and additional terms $U^a$. Here, we focus on evaluating $U^a$. Seven different types of attitudes are considered for evaluating the RBP, as outlined in

Table 3. Seven types of the fixed attitude of S/C.

Case	1	2	3	4	5	6	7
$\widehat{\mathit{q}}$	$\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right]$	$\left[\begin{array}{c}0\\ 1\\ 0\\ 0\end{array}\right]$	$\left[\begin{array}{c}0\\ 0\\ 1\\ 0\end{array}\right]$	$\left[\begin{array}{c}0\\ 0\\ 0\\ 1\end{array}\right]$	$\frac{1}{\sqrt{2}}\left[\begin{array}{c}1\\ 1\\ 0\\ 0\end{array}\right]$	$\frac{1}{\sqrt{2}}\left[\begin{array}{c}1\\ 0\\ 1\\ 0\end{array}\right]$	$\frac{1}{\sqrt{2}}\left[\begin{array}{c}1\\ 0\\ 0\\ 1\end{array}\right]$

, where $\widehat{\mathit{q}}$ denotes the quaternion. Different colors distinguish the gravity estimation models from each other. Blue-based, red-based, and yellow-based colors indicate the results calculating using the PM, EB, and TE models, respectively. Outcomes related to the x-, y-, and z-axes are represented with circle, square, and triangle markers, respectively.

It can be seen from Proposition 1 that $U$ is calculated by combining eigenvalues and their corresponding eigenvectors, with the MOI of the S/C also being one of the parameters. Figure 4 shows the almost proportional relationship between $U^{\mathit{a}}$ and the magnitude of ${\mathit{V}}_{xx}$. All results are represented, irrespective of small celestial bodies and the attitudes of S/C. The maximum eigenvalues are used as the magnitude of ${\mathit{V}}_{xx}$. The left figure depicts $U^a=a{\left|{\mathit{V}}_{xx}\right|}^b$, where $a$ and $b>0$ are real numbers. Both the x- and y-axes are displayed on a logarithmic scale to clearly visualize the tendency. On the right, it is evident that $b$ is nearly 1, indicating an almost proportional tendency between $U^a$ and the magnitude of ${\mathit{V}}_{xx}$. Figure 5 presents the relative $U^a$ and relative errors of $V$ for 25143 Itokawa, in a single plot. The relative $U^a$ represents the magnitude of $U^a$ over $V$ for each gravity estimation model, and it is denoted as a lighter color. The relative error of $V$ is the absolute error of the given gravity estimation model with respect to its potential value. The standard value for the absolute error is the potential calculated using the direct integration model. The relative error of $V$ is represented in a darker color. The attitude of the S/C is specified as $\widehat{\mathit{q}}={\left[1, 1, 0, 0\right]}^T/\sqrt{2}$. Unfortunately, the relative errors of the three gravity estimation models surpass the magnitude of $U^a$. This indicates that inaccurate gravity estimation models are unsuitable not only for potential estimation but also for estimating the RBP. Figure 6 illustrates the ratio of $U^a$ to $V$ (left) and the relative error of $V$ (right), with respect to radial distance in kilometers. All results from the three gravity estimation models are depicted in the same colors, with color variations denoting different small celestial bodies. $U^a$ is derived from ${\mathit{V}}_{xx}$, which decreases as the radial distance increases. Consequently, the $U^a$ ratio diminishes with increasing radial distance, irrespective of any specific gravity estimation models, as evident in Figure 6. On the other hand, as depicted on the right, the relative error values of $V$ estimated using each gravity estimation model are similar. Although the three gravity estimation models utilized in this study may not be adequate for the analysis on RBP, other models exhibiting relative errors below ${10}^{-5}$ are deemed suitable for analyzing the RBP near relatively small bodies, such as 99942 Apophis.

The propagated trajectories of the S/C with three different types of potentials are compared with each other to validate the proposed RBP. As all gravity estimation models considered in this study are rather limited in accuracy, a uniform sphere is chosen as an artificial small celestial body. This allows us to observe the effects of the additional terms in the RBP on the motion of the S/C, even when the PM model is adopted, as discussed in Section 5.1. The mass, radius, spin rate, and spin axis of the artificial/spherical small celestial body are selected as $4.19\times {10}^9 \mathrm{k}\mathrm{g}$, $1.2 \mathrm{k}\mathrm{m}$, $1.13\times {10}^{-5} \mathrm{r}\mathrm{a}\mathrm{d},$ and the z-axis, respectively. The initial position and velocity are chosen as $\left[1.30, 0, 0\right] \mathrm{k}\mathrm{m}$ and $\left[0, 0, 0\right] \mathrm{k}\mathrm{m}/\mathrm{s}$ in the principal axis frame, respectively, which correspond to a geosynchronous orbit in the restricted two-body context. The initial attitude of the S/C is case 7 in Table 3, and the initial angular velocity is set to be the same as the angular velocity of the artificial small celestial body, which ensures that the S/C maintains a consistent orientation relative to the artificial small celestial body. These initial conditions are suitable for observing the influence of additional terms in the RBP, i.e., the second and third terms in Equation (18). Now compared are three dynamic models: the orbital motion in the restricted two-body problem, the orbit–attitude coupled motion with the RBP, and the orbit–attitude coupled motion with the mutual potential. Figure 7 illustrates the deviations from the prescribed initial positions of the S/C in the principal axis frame and displays trajectories for $500 \mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$ in the xy-plane. The blue solid, red dashed, and yellow solid lines show the propagations with the RBP, the mutual potential, and the potential $V$, respectively. Since the initial condition is set to form a geosynchronous orbit, the propagated position with the potential $V$ remains in its initial position after $500 \mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$. The propagation with the RBP shows a bigger difference between the initial and final positions than the propagation with the mutual potential. This is because the propagation with the RBP does not consider the gravitational influence of the S/C on the motion of the artificial small celestial body. Additionally, both images in Figure 7 illustrate that the propagation with the mutual potential (rather than the propagation with the potential) shows a similar tendency to the propagation with the RBP. This implies that the analysis of orbit–attitude coupled motion with the RBP is a reasonable approximation and efficient alternative, when the mass of the S/C is considerably light but its size is not negligibly small, compared with the small celestial body.

6. Conclusions

This study has established a comprehensive framework for calculating the RBP, along with its first and second derivatives. The overall analysis does not take into account the motion of the small celestial body itself. The terminology ‘rigid body potential’ is defined to distinguish it from the mutual potential. S/Cs are assumed to be extended rigid bodies with finite volume to ensure their applicability within dynamic and control systems. The mass ratio between a S/C and a small celestial body ensures that their barycenter can be located at the center of the small celestial body. Additionally, the ‘restricted’ assumption, which stipulates that the S/C does not exert any gravitational force on the small celestial body, has been numerically validated to be reasonable by comparing the propagated trajectories with three different types of potentials (i.e., potential V, rigid body potential (RBP) U, and mutual potential) with each other. Three gravity estimation models based on the shape model of the small celestial body have been introduced: the PM, EB, and TE models. The formulations of the RBP, along with its first and second order derivatives are proposed and are implemented in conjunction with the three gravity estimation models. The Hessian matrix of the potential plays a crucial role in constructing the RBP. Analyses of the gravity estimation models and the RBP were conducted and numerically tested for eight small celestial bodies, chosen arbitrarily. The relative error of the potential decreases as the radial distance increases, irrespective of the gravity estimation model. However, when the radial distance in DU is similar, the relative error of the same gravity estimation model is also similar, regardless of the choice of small celestial body. On the other hand, the additional terms of the RBP decrease as the DU increases, which indicates dependence on the size of the small celestial body. This implies that accurate gravity estimation models are imperative for conducting RBP analysis. Nevertheless, the comparative analysis of the propagated trajectories with a uniform sphere as a small celestial body suggests that the RBP should be a reasonable approximation and an efficient alternative to the mutual potential for analyzing the orbit–attitude coupled motion of S/C ‘only’. This observation motivates us to further apply the proposed framework of RBP to other gravity estimation models, such as the mass concentration model, spherical harmonics model, and direct integration model.