Free-Vertex Tetrahedral Finite-Element Representation and Its Use for Estimating Density Distribution of Irregularly-Shaped Asteroids

Abstract

In this article, we present a free-vertex tetrahedral finite-element representation of irregularly shaped small bodies, which provides an alternative solution for estimating asteroid density distribution. We derived the transformations between gravitational potentials expressed by the free-vertex tetrahedral finite elements and the spherical harmonic functions. Inversely, the density of each free-vertex tetrahedral finite element can be estimated via the least-squares method, assuming a spherical harmonic gravitational function is present. The proposed solution is illustrated by modeling gravitational potential and estimating the density distribution of the simulated asteroid 216 Kleopatra.

1. Introduction

Asteroids, as remnants of the early stages of planet formation, have the potential to advance knowledge about the origin of the solar system. Knowledge of the internal structure of asteroids is of great significance for understanding the origin, evolution, and current state of asteroids [1]. Although the interior cannot be observed directly, other information can be used to assist in gaining an up-to-date understanding. For example, in the study of the internal structure of planetary moons, Durante et al. [2] used flyby data to determine the gravitational field to degree and order 5, which improved the estimated value of tidal Love number k2 and determined that Titan has a strongly differentiated interior. Cappuccio et al. [3] used Ganymede’s gravity field and tide response, determined by Doppler and range data from a 500 km circular orbit in the 3GM experiment simulation, to provide critical information for accurately modeling Ganymede’s internal structure. However, the size of most asteroids is smaller than these moons, and we have less prior knowledge of the internal structure, although many asteroid exploration missions have measured their gravity fields [4,5,6]. Additionally, many asteroids are considered to be rubble pile asteroids, and their internal composition and structure may be unconstrained [7]. Therefore, if it is possible to estimate the internal structure of asteroids with a tiny amount of information about them (only the gravity field and shape data are available); this will be of great help for research on asteroids. What factors are related to the estimation of the internal structure of asteroids? Since it is difficult to directly detect the internal structure, we can indirectly obtain the distribution of the internal structure from the density and gravity field information of asteroids. Density is indeed a fundamental property for understanding their composition and internal structure [1], which means that the density of asteroids may reflect their internal composition. Otherwise, the asteroid gravity field implies information about the density and internal structure of asteroids. An accurate gravity field will help asteroid exploration missions, such as in the design of frozen orbits [8]. Therefore, this paper aims to establish a new gravity field model of asteroids that can be used to estimate their density distribution. The performance of density estimation is investigated in two cases: (1) the internal structure with no prior information; and (2) the internal structure with prior information.

There are three fundamental solutions for modeling the gravitational field of an asteroid (or a comet): the spherical harmonic function [9], the polyhedron model [10], and the mascon method [11]. As each of these solutions has its pros and cons, a series of improved models have been proposed to either increase the computational efficiency or improve the fineness of computing irregularly shaped asteroid gravitation [12,13,14,15,16]. Among these improved models, the tetrahedral finite-element method (FEM) [17,18] attracts the most attention since it combines the advantages of the polyhedron (excellent asteroid irregular surface representation) and the mascon (non-uniform density distribution representation) approaches. The basic idea of the tetrahedral FEM is to use tetrahedral elements to fill up an asteroid’s spatial shape, and the faces of the outer-layer elements approximate the asteroid’s irregular surface. The volume of each tetrahedral element is determined by the specific meshing method with either coarse meshes or fine meshes. In the literature, there are two representative types of finite element meshes in the tetrahedral FEM, which are depicted in Figure 1: centroid-vertex meshes and free-vertex meshes. The former was first proposed by Scheeres et al. [16], in which at least one of the vertexes of all finite elements is constrained at the center of mass of the asteroid. However, this approach cannot directly express the internal density distribution of asteroids, such as the rubble piles model [19,20], the elongated bi-lobed model (4179 Toutatis [21] and comet 67P/Churyumov–Gerasimenko [22]) and the hypothetical model (the planar cut model, the surface layer model, the one core model, the two-core model, and the torus model, proposed by Takahashi et al. [23]). In contrast, the free-vertex meshes can remove the constraint assumed in the centroid-vertex meshes, providing a more flexible solution for generating tetrahedral finite elements. Therefore, this article focuses on presenting the free-vertex tetrahedral finite-element representation and its use for estimating the density distribution of irregularly shaped asteroids, which may be helpful for obtaining the internal structure of these small bodies.

The concept of the tetrahedral FEM for asteroid shape representation can be traced independently back to the concepts of the mascon method and the polyhedron model. To overcome the shortcomings of the mascon method in describing the shapes of asteroids, Pearl and Hitt proposed the tetrahedral FEM and higher-order polyhedral FEM meshes (which are called polyhedral-dual meshes, different from the polyhedral model) [24], and they found that the number of mascons derived from the FEM meshes can be reduced by 90% compared with the uniformly distributed mascons [14]. Pearl et al. [25] then compared the mascons derived from different FEM meshes with the polyhedron model and showed that the FEM-based mascons are more accurate for computing asteroid gravitation than the polyhedron model in some specific regions, such as a thin region that hugs the surface of the asteroid. Moreover, Rathinam et al. [26] suggested that the mascons based on the octree approach (a mesh generation method of FEM) are more computationally efficient than the polyhedron model. On the other hand, the tetrahedral FEM can also be obtained by improving the polyhedral model, which assumes a constant density distribution. Park et al. [27] used finite element spheres and finite element cubes to approximate the shape of asteroids, but there are gaps between the elements, and the surface is not smooth. Recently, Yu et al. [18] proposed an approach of computing the gravitational potentials of the binary asteroids using the free-vertex tetrahedral finite element and showed the effectiveness of the tetrahedral FEM. Yin et al. [17] found that the free-vertex tetrahedral finite element model can fit the structure of small bodies better than the polyhedron model.

The tetrahedral FEM with the variable meshing method is more suitable for the study of the asteroid density distribution than the polyhedron model and the mascon method, as it provides representations of the non-uniform density distribution of an asteroid with finer meshes. However, it is difficult to obtain the internal density distribution of the asteroid even if we can detect the gravity field and shape of the asteroid. A solution for estimating the density distribution of asteroids was proposed by Scheeres et al. [16]; the method uses the relationship between the spherical harmonic coefficients and the tetrahedral model to determine the density of each centroid-vertex tetrahedron by applying a least-squares method. Then, Park et al. [27] applied an alternative method using navigation measurements such as station-to-spacecraft tracking to estimate density distribution. However, the asteroid density distribution is quite unlikely to be accurately approximated by the centroid-vertex tetrahedral elements. Takahashi et al. [23,28] and Ledbetter et al. [29] divided the centroid-vertex tetrahedron element into layers, but the accuracy of density distribution estimation is still limited.

In summary, the motivation of this study is to investigate the use of the free-vertex tetrahedral FEM for representing an asteroid’s spatial shape and density distribution. It is deemed a follow-up advancement of density distribution estimation using the strip-shaped elements observed in the centroid-vertex tetrahedral FEM [16]. In Section 2, the free-vertex tetrahedral finite element gravitation model is introduced. In Section 3, we derive the transformations between gravitational potentials expressed by the free-vertex tetrahedral finite elements and the spherical harmonic functions, and a least-squares method is used to estimate the density of each free-vertex tetrahedral finite element. In Section 4, we consider the asteroid 216 Kleopatra [30] as an example to demonstrate the free-vertex tetrahedral FEM. Conclusions are presented in Section 5.

2. Free-Vertex Tetrahedral Finite Element Gravitation Model

The gravitational potential and force caused by each free-vertex tetrahedral finite element (which can also be called a tetrahedron) can be calculated by applying the polyhedron method [10], and then the values of each potential and force are added up to obtain the gravitational potential $U_{poly-tetra}$ and the gravitational force $\nabla U_{poly-tetra}$ for a specific asteroid, which are given below.

$$U_{poly-tetra}=\frac{1}{2}G\sum _{i=1}^n{\rho }_i\left(\sum _{e\in edges}{\mathbf{r}}_e^i·{\mathbf{E}}_e^i·{\mathbf{r}}_e^i·L_e^i-\sum _{f\in faces}{\mathbf{r}}_f^i·{\mathbf{F}}_f^i·{\mathbf{r}}_f^i·{\omega }_f^i\right)$$

(1)

$$\nabla U_{poly-tetra}=-G\sum _{i=1}^n{\rho }_i\left(\sum _{e\in edges}{\mathbf{E}}_e^i·{\mathbf{r}}_e^i·L_e^i-\sum _{f\in faces}{\mathbf{F}}_f^i·{\mathbf{r}}_f^i·{\omega }_f^i\right),$$

(2)

where G is the universal gravitational constant, i is the serial number of the tetrahedral elements, ${\rho }_i$ is the density of the numbered tetrahedral element, e represents the edge of the tetrahedron, and f represents the triangular face by three vertices of the tetrahedron.

$$\begin{array}{c}{\mathbf{F}}_f^i={\hat{n}}_f{\hat{n}}_f\hfill \\ {\mathbf{E}}_e^i={\hat{n}}_A{\hat{n}}_e^A+{\hat{n}}_B{\hat{n}}_e^B\hfill \\ L_e^i={\int }_e\frac{1}{r}ds={\int }_{P_l}^{P_j}\frac{1}{r}ds=ln\frac{{r^{\prime }}_l+{r^{\prime }}_j+e_{lj}}{{r^{\prime }}_l+{r^{\prime }}_j-e_{lj}}\hfill \\ {\omega }_f^i={\int \int }_{triangle}\frac{\Delta z}{r^3}dS=2·arctan\frac{{\mathbf{r}}_l·{\mathbf{r}}_j\times {\mathbf{r}}_k}{r_l{r}_j{r}_k+r_l\left({\mathbf{r}}_j·{\mathbf{r}}_k\right)+r_j\left({\mathbf{r}}_k·{\mathbf{r}}_l\right)+r_k\left({\mathbf{r}}_l·{\mathbf{r}}_j\right)},\hfill \end{array}$$

(3)

where ${\mathbf{F}}_f^i$ is the dyad related to the triangular face, ${\mathbf{E}}_e^i$ is the dyad related to the edge, $L_e^i$ is the dimensionless factor of each edge, ${\omega }_f^i$ is the dimensionless factor of each face, ${\hat{n}}_f$ represents the outward-pointing surface normal direction vector of the face, and ${\hat{n}}_e^f$ is the normal direction vector of the edge of the face lying in the face plane, which is orthogonal to both the face normal vector and the along-the-edge vector, and points outward. $\mathrm{A}$ and $\mathrm{B}$ are the symbols of the two faces that have common edges, respectively. ${\mathbf{r}}_l,{\mathbf{r}}_j$ and ${\mathbf{r}}_k$ represent the vectors from the field point to the three vertices of the triangle surface, respectively, and $r_l,r_j$, and $r_k$ are the modulus of the corresponding vectors. The edges and faces, as well as the normal direction vectors, are shown in Figure 2.

The free-vertex tetrahedral finite element model is obtained by mesh generation methods [31]. The structured mesh is composed of horizontal and vertical lines that cross orthogonally at intersections called nodes. Generally, the structured mesh generation process is best achieved by discretizing the physical domain defined by the square or rectangular boundary; thus, it is mainly suitable for simple geometric shapes. In unstructured meshes, the nodes used to discretize the two-dimensional physical domain change from being ordered orthogonally in the structured mesh to seemingly randomly placed. These nodes are connected to other nodes via triangular or quadrilaterally shaped subdomains or elements. In the three-dimensional physical domain, the most common types of elements are tetrahedrons and hexahedrons. The generation of unstructured meshes requires more thought and effort than structured meshes. In general, one puts more nodes (or elements) near surfaces and in regions where activity (or steep gradients) is likely to occur. Therefore, the unstructured mesh is a better choice for irregularly shaped asteroids. The unstructured mesh can be generated by the octree algorithm [32], the Delaunay algorithm [31] or the advancing front algorithm [33]. The Octree algorithm uses a cube to cover the entire shape area that needs to be meshed, and then continuously subdivides the cube according to the requirements of the mesh size. The advancing front algorithm generates tetrahedral elements from the boundary of the meshed area and continues to advance to the center until the entire area is covered by tetrahedral elements. The Delaunay algorithm requires the tetrahedron to satisfy the Delaunay criterion [31]; that is, there are no other points in the circumscribed ball of each tetrahedron, except for its own four vertices. A significant advantage of the Delaunay algorithm is that it can make the minimum angle of each tetrahedron element in the mesh system constituted by a given point set as large as possible, allowing us to obtain tetrahedron elements with as high a quality as possible.

The asteroid 216 Kleopatra is taken to generate the free-vertex tetrahedral meshes as an example. As shown in Figure 3, the shape model of Kleopatra is modeled as the polyhedron model and the free-vertex tetrahedral FEMs. The meshes of the tetrahedral FEMs are generated by the Delaunay triangulation algorithm [31] with different maximum element growth rates (the difference between the sizes of two adjacent elements).

Since we generate the free-vertex tetrahedral FEMs, if we obtain the prior information about the density distribution, we can calculate the gravitational field of the asteroid. Although the meshes of these models are different, the gravity field calculated by each FEM is valid if the density distribution of the corresponding mesh can be obtained. This also shows that the meshes of the free-vertex tetrahedral FEMs are not restricted by the internal structures of the asteroid. On the contrary, if the internal structures are known, the divided tetrahedral finite elements can be combined into the corresponding structures, which is also the reason that we propose the free-vertex tetrahedral FEM.

3. Transformation for Gravity Field Expressions

3.1. Transformation from Free-Vertex Tetrahedral FEM to Spherical Harmonics with Simplex Integral

The spherical harmonic expansion of the gravitational potential [34] is expressed by:

$$U_{spherical}=\frac{G{M}_A}{r}\left\{\sum _{n=0}^N\sum _{m=0}^n{\left(\frac{r_e}{r}\right)}^n{P}_{nm}\left(sin\phi \right)\left[C_{nm}cosm\lambda +S_{nm}sinm\lambda \right]\right\},$$

(4)

where r is the distance of the field point to the center of the expansion, and $\lambda$ and $\phi$ are the longitude and latitude of the field point, respectively. $P_{nm}\left(sin\phi \right)$ is the associated Legendre polynomial; $r_e$ is the radius of the reference sphere, which represents the convergence range of the series; $M_A$ is the total mass of the asteroid; n and m are the degree and order of the expansion, respectively; and $C_{nm}$ and $S_{nm}$ are the spherical harmonic coefficients that are specified up to degree N. The coefficient $C_{00}=1$ and the total number of gravitational coefficients is ${\left(N+1\right)}^2$, with $2n+1$ coefficients at each degree n.

The fully normalized harmonic coefficients, ${\overline{C}}_{n,m}$ and ${\overline{S}}_{n,m}$, can be calculated analytically by approximating the extended body with a swarm of points [35]:

$$\left[\begin{array}{c}{\overline{C}}_{n,m}\hfill \\ {\overline{S}}_{n,m}\hfill \end{array}\right]=\sum _{i\in \mathrm{point}\mathrm{masses}}\left[\begin{array}{c}{\overline{c}}_{n,m}\left({x^{\prime }}_i,{y^{\prime }}_i,{z^{\prime }}_i\right)\hfill \\ {\overline{s}}_{n,m}\left({x^{\prime }}_i,{y^{\prime }}_i,{z^{\prime }}_i\right)\hfill \end{array}\right]M_i,$$

(5)

where $M_i$ is the point’s mass, and ${\overline{c}}_{n,m}$ and ${\overline{s}}_{n,m}$ are integrands of the point’s normalized spherical harmonic coefficients. $\left(x^{\prime },y^{\prime },z^{\prime }\right)$ are the coordinates that represent the point.

Then, the point is replaced by a tetrahedron. The density ${\rho }_i$ of each tetrahedron is assumed to be constant, and the coefficients ${\overline{C}}_{n,m}$ and ${\overline{S}}_{n,m}$ can be written as the following integral form:

$$\left[\begin{array}{c}{\overline{C}}_{n,m}\\ {\overline{S}}_{n,m}\end{array}\right]={\int \int \int }_{\mathrm{extended}\mathrm{body}}\left[\begin{array}{c}{\overline{c}}_{n,m}\hfill \\ {\overline{s}}_{n,m}\hfill \end{array}\right]\mathrm{d}m={\int \int \int }_{\mathrm{extended}\mathrm{body}}{\rho }_i\left[\begin{array}{c}{\overline{c}}_{n,m}\hfill \\ {\overline{s}}_{n,m}\hfill \end{array}\right]\mathrm{d}{x}^{\prime }\mathrm{d}{y}^{\prime }\mathrm{d}{z}^{\prime }.$$

(6)

The recurrent relationships for integrands ${\overline{c}}_{n,m}$ and ${\overline{s}}_{n,m}$ proposed by Werner et al. [35] are as follows:

$$\begin{array}{c}\mathbf{sectorial}:\hfill \\ n=0:\left[\begin{array}{c}{\overline{c}}_{0,0}\hfill \\ {\overline{s}}_{0,0}\hfill \end{array}\right]=\frac{1}{M}\left[\begin{array}{c}1\hfill \\ 0\hfill \end{array}\right]\hfill \\ n=1:\left[\begin{array}{c}{\overline{c}}_{1,1}\hfill \\ {\overline{s}}_{1,1}\hfill \end{array}\right]=\frac{1}{\sqrt{3}M}\left[\begin{array}{c}{x}^{\prime }/\phantom{x^{\prime }{r}_e}{r}_e\hfill \\ y^{\prime }/\phantom{y^{\prime }{r}_e}{r}_e\hfill \end{array}\right]\hfill \\ n>1:\left[\begin{array}{c}{\overline{c}}_{n,n}\hfill \\ {\overline{s}}_{n,n}\hfill \end{array}\right]=\frac{2n-1}{\sqrt{2n\left(2n+1\right)}}\left[\begin{array}{cc}{x}^{\prime }/\phantom{x^{\prime }{r}_e}{r}_e& -y^{\prime }/\phantom{-y^{\prime }{r}_e}{r}_e\\ y^{\prime }/\phantom{y^{\prime }{r}_e}{r}_e& x^{\prime }/\phantom{x^{\prime }{r}_e}{r}_e\end{array}\right]\left[\begin{array}{c}{c}_{n-1,n-1}\hfill \\ s_{n-1,n-1}\hfill \end{array}\right]\hfill \\ \mathbf{vertical}:\hfill \\ \left[\begin{array}{c}{\overline{c}}_{n,m}\hfill \\ {\overline{s}}_{n,m}\hfill \end{array}\right]=\left(2n-1\right)\sqrt{\frac{\left(2n-1\right)}{\left(2n+1\right)\left(n+m\right)\left(n-m\right)}}\frac{z^{\prime }}{r_e}\times \left[\begin{array}{c}{\overline{c}}_{n-1,m}\hfill \\ {\overline{s}}_{n-1,m}\hfill \end{array}\right]\hfill \\ -\sqrt{\frac{\left(2n-3\right)\left(n+m-1\right)\left(n-m-1\right)}{\left(2n+1\right)\left(n+m\right)\left(n-m\right)}}{\left(\frac{r^{\prime }}{r_e}\right)}^2\left[\begin{array}{c}{\overline{c}}_{n-2,m}\hfill \\ {\overline{s}}_{n-2,m}\hfill \end{array}\right]\hfill \\ \mathbf{subdiagonal}:\hfill \\ \left[\begin{array}{c}{\overline{c}}_{n,n-1}\hfill \\ {\overline{s}}_{n,n-1}\hfill \end{array}\right]=\frac{2n-1}{\sqrt{2n+1}}\frac{z^{\prime }}{r_e}\left[\begin{array}{c}{\overline{c}}_{n-1,n-1}\hfill \\ {\overline{s}}_{n-1,n-1},\hfill \end{array}\right]\hfill \end{array}$$

(7)

where $r_e$ is the radius of the reference sphere, and $r^{\prime }=\sqrt{{x^{\prime }}^2+{y^{\prime }}^2+{z^{\prime }}^2}$.

Although the recurrence relationships of spherical harmonic coefficients are given, they cannot be applied directly to calculate ${\overline{c}}_{n,m}$ and ${\overline{s}}_{n,m}$ of the free-vertex tetrahedral finite elements proposed in this article as the coordinates $\left(x^{\prime },y^{\prime },z^{\prime }\right)$ cannot be obtained directly for a tetrahedron. Furthermore, as the asteroid model used in Werner et al. [35] and Scheeres et al. [16] is the centroid-vertex tetrahedral finite element, we cannot simply refer to their calculation results of the spherical harmonic coefficients. Therefore, we re-derived the calculation of spherical harmonic coefficients based on the free-vertex tetrahedral FEM in this study. The coordinates $\left(x^{\prime },y^{\prime },z^{\prime }\right)$ of the given free-vertex tetrahedral elements $P_0{P}_1{P}_2{P}_3$ are denoted as $\left(x_0,y_0,z_0\right),\left(x_1,y_1,z_1\right),\left(x_2,y_2,z_2\right)$ and $\left(x_3,y_3,z_3\right)$, which are different from those defined in [35]: $\left(0,0,0\right),\left(x_1,y_1,z_1\right),\left(x_2,y_2,z_2\right)$ and $\left(x_3,y_3,z_3\right)$. In order to solve Equation (6), variable conversion is needed, and the coordinates $\left(x^{\prime },y^{\prime },z^{\prime }\right)$ are expressed with standard simplex coordinates $\left(u_0,u_1,u_2,u_3\right)$ by the relation [36]:

$$\left(\begin{array}{c}1\hfill \\ x^{\prime }\hfill \\ y^{\prime }\hfill \\ z^{\prime }\hfill \end{array}\right)=\left(\begin{array}{cccc}1& 1& 1& 1\\ x_0& x_1& x_2& x_3\\ y_0& y_1& y_2& y_3\\ z_0& z_1& z_2& z_3\end{array}\right)\left(\begin{array}{c}{u}_0\hfill \\ u_1\hfill \\ u_2\hfill \\ u_3\hfill \end{array}\right).$$

(8)

Here, the simplex refers to the free-vertex tetrahedral finite element, and the standard simplex is defined with coordinates $\left(u_0,u_1,u_2,u_3\right):\left(0,0,0\right),\left(1,0,0\right),\left(0,1,0\right),\left(0,0,1\right)$.

The expansion of Equation (8) is:

$$\begin{array}{c}1=u_0+u_1+u_2+u_3\hfill \\ x^{\prime }\left(u_0,u_1,u_2,u_3\right)=x_1{u}_1+x_2{u}_2+x_3{u}_3+x_0{u}_0\hfill \\ y^{\prime }\left(u_0,u_1,u_2,u_3\right)=y_1{u}_1+y_2{u}_2+y_3{u}_3+y_0{u}_0\hfill \\ z^{\prime }\left(u_0,u_1,u_2,u_3\right)=z_1{u}_1+z_2{u}_2+z_3{u}_3+z_0{u}_0.\hfill \end{array}$$

(9)

We use Equation (9) to replace $x^{\prime }$, $y^{\prime }$ and $z^{\prime }$ in Equation (7). The integrands ${\overline{c}}_{n,m}$ and ${\overline{s}}_{n,m}$ can be represented as follows:

$$\left[\begin{array}{c}{\overline{c}}_{n,m}\left(x^{\prime },y^{\prime },z^{\prime }\right)\hfill \\ {\overline{s}}_{n,m}\left(x^{\prime },y^{\prime },z^{\prime }\right)\hfill \end{array}\right]\to \left[\begin{array}{c}{\overline{c}}_{n,m}\left(u_0,u_1,u_2,u_3\right)\hfill \\ {\overline{s}}_{n,m}\left(u_0,u_1,u_2,u_3\right)\hfill \end{array}\right]=\sum _{i_0+i_1+i_2+i_3=n}\left[\begin{array}{c}{\overline{\alpha }}_{i_0,i_1,i_2,i_3}\hfill \\ {\overline{\beta }}_{i_0,i_1,i_2,i_3}\hfill \end{array}\right]{u_0}^{i_0}{u_1}^{i_1}{u_2}^{i_2}{u_3}^{i_3}.$$

(10)

Here, ${\overline{\alpha }}_{h,j,k,l}$ and ${\overline{\beta }}_{h,j,k,l}$ are the tetranomial coefficients of coordinates $u_0,u_1,u_2$, and $u_3$. We provide examples of these in Appendix A.

Then, Equation (6) can be rewritten as follows:

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{\overline{C}}_{n,m}\hfill \\ {\overline{S}}_{n,m}\hfill \end{array}\right]& ={\int \int \int }_{\mathrm{extended}\mathrm{body}}\left[\begin{array}{c}{\overline{c}}_{n,m}\hfill \\ {\overline{s}}_{n,m}\hfill \end{array}\right]\mathrm{d}m\hfill \\ \hfill & =\sum _{\mathrm{simplices}}{\rho }_i\left(detJ\sum _{i_0+i_1+i_2+i_3=n}\left[\begin{array}{c}{\overline{\alpha }}_{i_0,i_1,i_2,i_3}\hfill \\ {\overline{\beta }}_{i_0,i_1,i_2,i_3}\hfill \end{array}\right]\times {\int \int \int }_{\mathrm{standard}\mathrm{simplex}}^{1=u_0+u_1+u_2+u_3}{u_0}^{i_0}{u_1}^{i_1}{u_2}^{i_2}{u_3}^{i_3}d{u}_1\mathrm{d}{u}_2\mathrm{d}{u}_3\right).\hfill \end{array}$$

(11)

Here, the matrix J is the Jacobian matrix.

$$\begin{array}{c}J\equiv \frac{\partial \left(x^{\prime },y^{\prime },z^{\prime }\right)}{\partial \left(u_1,u_2,u_3\right)}=\left[\begin{array}{ccc}\frac{\partial x}{\partial u_1}& \frac{\partial x}{\partial u_2}& \frac{\partial x}{\partial u_3}\\ \frac{\partial y}{\partial u_1}& \frac{\partial y}{\partial u_2}& \frac{\partial y}{\partial u_3}\\ \frac{\partial z}{\partial u_1}& \frac{\partial z}{\partial u_2}& \frac{\partial z}{\partial u_3}\end{array}\right]\hfill \\ detJ=\left|\begin{array}{ccc}{x}_1-x_0& x_2-x_0& x_3-x_0\\ y_1-y_0& y_2-y_0& y_3-y_0\\ z_1-z_0& z_2-z_0& z_3-z_0\end{array}\right|=\left|\begin{array}{cccc}1& x_0& y_0& z_0\\ 1& x_1& y_1& z_1\\ 1& x_2& y_2& z_2\\ 1& x_3& y_3& z_3.\end{array}\right|\hfill \end{array}$$

(12)

Since the tetrahedral simplex integral is:

$$\begin{array}{c}{\int }_{\mathrm{simplex}}^{1=u_0+u_1+u_2+u_3}{u_0}^{i_0}{u_1}^{i_1}{u_2}^{i_2}{u_3}^{i_3}d{u}_1\mathrm{d}{u}_2\mathrm{d}{u}_3\hfill \\ =\frac{i_0!i_1!i_2!i_3!}{\left(i_0+i_1+i_2+i_3+3\right)!},\hfill \end{array}$$

(13)

the expression of the spherical harmonic coefficients for the free-vertex tetrahedral finite element model is:

$$\left[\begin{array}{c}{\overline{C}}_{nm}\hfill \\ {\overline{S}}_{nm}\hfill \end{array}\right]=\sum _{\mathrm{simplics}}{\rho }_i\left(\frac{detJ}{\left(n+3\right)!}\sum _{i_0+i_1+i_2+i_3=n}{i}_0!i_1!i_2!i_3!\left[\begin{array}{c}{\overline{\alpha }}_{i_0,i_1,i_2,i_3}\hfill \\ {\overline{\beta }}_{i_0,i_1,i_2,i_3}\hfill \end{array}\right]\right).$$

(14)

The integration of the simplex in Equation (13) is shown in Appendix B. In the process of analytically calculating simplex integrals, the tetranomial $i_0!i_1!i_2!i_3!$ becomes increasingly complex as the degree of the spherical harmonic coefficients grows, and the dimensions of tetranomial coefficients become enormous due to the recurrent relationships in Equation (7), which requires large amounts of computing resources. In order to save computation time, the analytical solution can be approximated by numerical integration methods in finite element analysis [37], such as the Hammer integral [38,39]. This study used the Hammer integral to simplify calculations for the research. The details of the Hammer integral are shown in Appendix C.

3.2. Inverse Transformation from Spherical Harmonics to Free-Vertex Tetrahedral FEM with the Least-Squares Method

From Equation (14), an inverse transformation can be used to obtain the density of each free-vertex tetrahedron. The least-squares method proposed by Scheeres et al. [16] was used to evaluate the density distribution of the asteroid.

To estimate the density, the serial numbers of tetrahedrons and vertex coordinates were obtained. The usual form of the spherical harmonic gravitational potential is shown in Equation (4), and the potential needs to be rewritten by the tetrahedral finite element model:

$$U_{sphe-tetra}=\sum _i^p\frac{G{V}_i{\rho }_i}{r}\left[\sum _{n=0}^N\sum _{m=0}^n{\left(\frac{r_e}{r}\right)}^n{\overline{P}}_{nm}\left(sin\theta \right)\left({\overline{C}}_{nm}^{i}cosm\lambda +{\overline{S}}_{nm}^{i}sinm\lambda \right)\right],$$

(15)

where i is the index of the tetrahedron, p is the number of tetrahedrons, N is the degree of the highest calculated spherical harmonic coefficient, and ${\overline{C}}_{nm}^i$ and ${\overline{S}}_{nm}^i$ are the normalized spherical harmonic coefficients of the corresponding tetrahedron. Since the volume of the tetrahedron $V_i$ can be calculated from the geometric shape, due to the constant density hypothesis for each tetrahedron, the only unknown in the gravitational potential is the density ${\rho }_i$.

The expression of the gravitational potential is different from the usual spherical harmonic gravitational potential. The usual form of Equation (4) is necessary to calculate the spherical harmonic coefficients of the whole asteroid and to consider the total mass of the asteroid, while in the free-vertex tetrahedral FEM, it is necessary to calculate the spherical harmonic coefficients corresponding to each tetrahedron separately, and then the potential of each tetrahedron is summed. Therefore, M refers to the mass of the tetrahedron in Equation (7).

Ideally, for the same asteroid, the gravitational potentials calculated by different models at the same location should be the same. For the inverse transformation, a least-squares method is used.

The loss function of the general least-squares method is given as:

$$J\left(\hat{\mathit{X}}\right)={\left(\mathbf{Z}-\mathit{H}\hat{\mathit{X}}\right)}^T\left(\mathit{Z}-\mathit{H}\hat{\mathit{X}}\right)=min.$$

(16)

Here, $\mathit{X}$ is the state vector, and the density ${\rho }_i$ is the element of the vector. $\hat{\mathit{X}}$ is the estimation of the state vector. $\mathit{Z}$ is the measurement vector, and $\mathit{H}$ is the measurement matrix. In this study, the measurement is the real gravitational potential of the asteroid $U_{spherical}$ in Equation (4), and the calculated measurement $\mathit{H}\hat{\mathit{X}}$ is the gravitational potential $U_{sphe-tetra}$ in Equation (15).

In order to minimize the function of the above equation, it is necessary to satisfy:

$${\left.\frac{\partial J}{\partial \mathit{X}}\right|}_{\mathit{X}=\hat{\mathit{X}}}=-2{\mathit{H}}^T\left(\mathit{Z}-\mathit{H}\hat{\mathit{X}}\right).$$

(17)

Therefore, under the ideal situation, the estimated state vector $\hat{\mathit{X}}$ can be expressed as:

$$\hat{\mathit{X}}={\left({\mathit{H}}^T\mathit{H}\right)}^{-1}{\mathit{H}}^T\mathit{Z}.$$

(18)

We can find that:

$$\mathit{H}=\left[\begin{array}{cccc}\frac{V_i}{M}{\overline{C}}_{10}^1& \frac{V_i}{M}{\overline{C}}_{10}^2& \cdots & \frac{V_i}{M}{\overline{C}}_{10}^p\\ \frac{V_i}{M}{\overline{C}}_{11}^i& \frac{V_i}{M}{\overline{C}}_{11}^2& \cdots & \frac{V_i}{M}{\overline{C}}_{11}^p\\ \frac{V_i}{M}{\overline{S}}_{11}^i& \frac{V_i}{M}{\overline{S}}_{11}^2& \cdots & \frac{V_i}{M}{\overline{S}}_{11}^p\\ ⋮& ⋮& ⋮& ⋮\end{array}\right]$$

(19)

Then, we can use the spherical harmonic coefficients to express Equation (18) and simplify the estimated solution as:

$$\mathbf{\rho }={\left({\left[\frac{V_i}{M}{\tilde{\mathit{C}}}_{nm}^i\right]}^T\left[\frac{V_i}{M}{\tilde{\mathit{C}}}_{nm}^i\right]\right)}^{-1}{\left[\frac{V_i}{M}{\tilde{\mathit{C}}}_{nm}^i\right]}^T\left[{\tilde{\mathit{C}}}_{nm}\right],$$

(20)

where $\mathbf{\rho }$ is the density vector, $V_i$ is the volume of the tetrahedron, M is the mass of the asteroid, and ${\tilde{\mathit{C}}}_{nm}$ is a vector representation of the spherical harmonic coefficients: $\left[{\tilde{\mathit{C}}}_{nm}\right]=\left[{\overline{C}}_{10,}{\overline{C}}_{11},{\overline{S}}_{11},{\overline{C}}_{20},{\overline{C}}_{21},{\overline{S}}_{21},\cdots ,{\overline{C}}_{NN},{\overline{S}}_{NN}\right]$.

If we measure the spherical harmonic coefficients of the asteroid, we can use weighted least-squares to calculate the state, which can be expressed as:

$$\mathbf{\rho }={\left({\left[\frac{V_i}{M}{\tilde{\mathit{C}}}_{nm}^i\right]}^T{\mathit{W}}_{\tilde{C}}\left[\frac{V_i}{M}{\tilde{\mathit{C}}}_{nm}^i\right]\right)}^{-1}{\left[\frac{V_i}{M}{\tilde{\mathit{C}}}_{nm}^i\right]}^T{\mathit{W}}_{\tilde{\mathit{C}}}\left[{\tilde{\mathbf{C}}}_{nm}\right],$$

(21)

where ${\mathit{W}}_{\tilde{\mathit{C}}}$ is the covariance matrix of the measured true spherical harmonic coefficients.

Then, we can use $\mathit{K}$ and $\mathit{y}$ to denote the matrixes:

$$\begin{array}{c}\mathit{K}={\left[\frac{V_i}{M}{\tilde{\mathit{C}}}_{nm}^i\right]}^T{\mathit{W}}_{\tilde{\mathit{C}}}\left[\frac{V_i}{M}{\tilde{\mathit{C}}}_{nm}^i\right]\hfill \\ \mathit{y}={\left[\frac{V_i}{M}{\tilde{\mathit{C}}}_{nm}^i\right]}^T{\mathit{W}}_{\tilde{\mathit{C}}}\left[{\tilde{\mathit{C}}}_{nm}\right].\hfill \end{array}$$

(22)

Equation (21) can be abbreviated as:

$$\mathbf{\rho }={\mathit{K}}^{-1}\mathit{y}.$$

(23)

The dimension of $\mathbf{\rho }$ is $\left[p\times 1\right]$, and the dimension of ${\tilde{\mathit{C}}}_{nm}$ is $\left({\left(N+1\right)}^2-1\right)\times 1$. If the number of tetrahedrons is equal to ${(N+1)}^2-1$, $\mathit{K}$ is an invertible matrix, and there is a unique solution. If the number of constraints (tetrahedrons) is greater than the number of unknowns (the density of each tetrahedron), which is an overdetermined problem, the normal matrix is invertible, and the density of each tetrahedron can be obtained. When the number of tetrahedrons exceeds ${(N+1)}^2-1$, the number of equations is less than the number of unknowns. If the normal matrix is not invertible and the system is under-determined, then the density cannot be solved by normal methods. In order to solve both cases at the same time, singular value decomposition (SVD) was used to solve Equation [16].

The singular value decomposition was performed by [40]:

$$\mathit{K}=\mathit{U}\mathit{S}{\mathit{V}}^T,$$

(24)

where $\mathit{U}$ is an $l\times m$ orthogonal matrix and $\mathit{V}$ is an $m\times m$ orthonormal matrix. Matrix S is an $m\times m$ diagonal matrix containing all the solved singular values, but if the ratio of the largest singular value to the smallest singular value is more than ${10}^{16}$, the condition number of the entire matrix will be abnormally large, which will cause the original problem to be poorly conditioned. Thus, the equation solution is very sensitive to disturbances. We selected condition number ${10}^{13}$ to obtain an appropriate ${\tilde{\mathit{S}}}^{-1}$.

The final solution of the least-squares method is:

$$\mathbf{\rho }=\mathit{V}{\tilde{\mathit{S}}}^{-1}{\mathit{U}}^T\mathit{y}.$$

(25)

From this derivation process, it can be found that the accuracy of the density estimation depends on the number of tetrahedrons and the degree and order of spherical harmonic coefficients.

4. Density Distribution Estimation

In this section, a simulation in 216 Kleopatra is presented to verify the effectiveness of the free-vertex tetrahedral FEM for density distribution estimation. The mass of the asteroid 216 Kleopatra is $4.64\times {10}^{18}\mathrm{kg}$, and the reference radius of Kleopatra is 112 km. The coarse mesh that consists of 1467 tetrahedrons shown in Figure 3c was used. The volume of Kleopatra calculated by the proposed tetrahedral FEM is $6.8144\times {10}^5{\mathrm{km}}^3$; thus, its bulk density (the mass to volume ratio of Kleopatra) is $6.809\mathrm{g}/{\mathrm{cm}}^3$. In these simulations, tetrahedrons can be used alone or assembled into polyhedrons. Then, we assigned densities to each tetrahedron (or polyhedron) to compute the “true” gravity potential, and the coefficients of the spherical harmonic function can be obtained from Equation (14). Then, the least-squares method was employed on the same shape discretization to verify that we recovered the proper densities for each tetrahedron (or polyhedron). Since the true gravitational field was fully known in our tests, the covariance matrix ${\mathit{W}}_{\tilde{\mathit{C}}}$ became an identity matrix.

4.1. Density Distribution Estimation for Each Tetrahedron

We first estimated the density distribution when each tetrahedron acts as an independent unit of the free-vertex tetrahedral finite-element model. This indicates that the densities for elements of a mass distribution are estimated without prior information on the geometry of the distribution. This is of great significance for asteroid exploration missions. At the early stage of these missions, we do not particularly understand the internal structure of the target asteroid. If we can use shape information and the spherical harmonic gravitational field information in the mission to detect the asteroid, we can find out whether the interior of the asteroid is uniform, or whether there are density mutations or hollows.

Since the number of tetrahedrons is 1467, which is larger than the dimensions of the vector of the spherical harmonic coefficients ($n\le 37$), the problem of estimating the density distribution is under-determined. Therefore, we needed to study the performance of the method under this condition. The effectiveness and limitations of this method were investigated by the sensitivity of estimating a density mutation and the estimation accuracy of different numbers of density mutations by different degrees and orders of gravity field coefficients.

4.1.1. Density Mutation Assumptions

Three cases were chosen, namely, single mutations (No. 1426), double mutations (No. 35/1426), and ten mutations (No. 2/35/87/624/675/790/960/1307/1401/1426), to estimate the density mutation. We assumed an a priori density of $10\mathrm{g}/{\mathrm{cm}}^3$ for the mutations in the three cases; the remaining tetrahedrons without density mutations took the original bulk density. These three assumptions are shown in Figure 4a.

Figure 4b shows the location of the three density mutation cases in each tetrahedral model, and different densities correspond to different colors. To display the result better, the tetrahedrons with density mutations selected in this study are on the surface of the asteroid. The tetrahedrons with single and double mutations are distributed at both ends of the asteroid, and the ten mutations are evenly distributed. For example, tetrahedrons No. 2/35/1401 and No. 1426 are at the ends, and tetrahedrons No. 624/790/960 are near the center of the asteroid.

4.1.2. Validation Test

We tested the sensitivity of our approach with the under-determined case (using a single mutation as an example). Tetrahedron No. 1426 was assumed to have an inconstant density from 0 to $10\mathrm{g}/{\mathrm{cm}}^3$, and the densities of the remaining 1466 tetrahedrons were $6.809\mathrm{g}/{\mathrm{cm}}^3$. The true spherical harmonic gravitational field was assigned to degree and order 8. The density of all 1467 tetrahedrons was initialized and assumed to be equal to $6.809\mathrm{g}/{\mathrm{cm}}^3$. We computed the correction to this initial assumption.

Figure 5 shows the estimated density of tetrahedron No. 1426 and the minimum, maximum, and mean of the other 1466 tetrahedrons as a function of the assumed density of that tetrahedron. The red solid line with the red points is the estimated density of tetrahedron No. 1426, and the blue dotted line with the blue circle is the mean value of the remaining tetrahedrons. The yellow and magenta dotted lines with symbols of the yellow star and the magenta asterisk show the minimum and maximum values of the remaining tetrahedrons, respectively.

Table 1. Validation test result of the single-mutation case.

Assumed Density of Tetrahedron No. 1426 (g/${\mathrm{cm}}^3$)	Estimated Density of Tetrahedron No. 1426 (g/${\mathrm{cm}}^3$)	Absolute Error of the Estimated Density of Tetrahedron No. 1426 (g/${\mathrm{cm}}^3$)
0	5.1111	5.1111
1	5.3604	4.3604
2	5.6098	3.6098
3	5.8592	2.8592
4	6.1085	2.1085
5	6.3579	1.3579
6	6.6073	0.6073
7	6.8566	0.1434
8	7.1059	0.8940
9	7.3554	1.6446
10	7.6047	2.3953

lists the absolute error of the estimated density of tetrahedron No. 1426. The results show that, when the density of a tetrahedron is different from the others in an under-determined system, we can separate it from other tetrahedrons, although the density of the special tetrahedron cannot be restored accurately. In particular, when the assumed value is very close to the initial value, the accuracy of density estimation is the best.

4.1.3. Influence of High-Order Spherical Harmonics

In these tests, the influence of the different degrees and orders of the spherical harmonic coefficients on the results of estimating different number density mutations were inspected. Under the premise of three density mutation assumptions (single mutation, double mutations, and ten mutations), we assumed that the densities of the tetrahedrons with density mutations were $10\mathrm{g}/{\mathrm{cm}}^3$, and the density of the other 1466/1465 and 1657 tetrahedrons was assumed to be $6.809\mathrm{g}/{\mathrm{cm}}^3$, respectively. The true spherical harmonic gravitational field was specified up to degree and order 40. Then, we initially assumed that the densities of all 1467 tetrahedrons were equal to the bulk density $6.809\mathrm{g}/{\mathrm{cm}}^3$ and calculated the correction for these initial values under the three density mutation assumptions cases. The density estimation results of the mutations and the other tetrahedrons as a function of the harmonic degree used for the “true” gravity field are shown in Figure 6, Figure 7 and Figure 8, respectively. The solid lines of the three figures show the estimated densities of the mutational tetrahedrons, and the dotted lines with the blue circle, the yellow star, and the magenta asterisk represent the mean, minimum and maximum values of the remaining tetrahedrons, respectively.

Table 2. The error statistics of the mutations for three density mutation assumptions.

Degree of Harmonics Coefficients	Single Mutation	Double Mutations			10 Mutations
	No. 1426 g/${\mathrm{cm}}^{\mathbf{3}}$	No. 35 g/${\mathrm{cm}}^{\mathbf{3}}$	No. 1426 g/${\mathrm{cm}}^{\mathbf{3}}$	Mean g/${\mathrm{cm}}^{\mathbf{3}}$	Mean g/${\mathrm{cm}}^{\mathbf{3}}$	Max g/${\mathrm{cm}}^{\mathbf{3}}$	Min g/${\mathrm{cm}}^{\mathbf{3}}$
8	2.3952	2.7569	2.3852	2.5711	2.8880	3.1015	2.5712
15	1.1614	1.7914	1.1782	1.4848	2.1138	3.0453	1.1184
20	0.6193	1.0616	0.6184	0.8399	1.7505	2.9491	0.6786
40	0.0155	0.0325	0.0161	0.0243	1.2106	2.894	0.0033

shows the error statistics of the mutations for three density mutation assumptions (we chose degrees 8, 15, 20, and 40 as examples). Figure 9 shows the density estimation results for the three mutation cases by the spherical harmonic gravitational field with degree and order 20. The left plot is the density of each tetrahedron, and the right plot shows the locations of the tetrahedrons, in which the density value is higher as the color becomes closer to yellow.

It should be noted that the densities of the mutant tetrahedrons were estimated more accurately as the degree became higher. In contrast, as the number of mutant tetrahedrons increased and their position distribution in the asteroid became increasingly random, the error of density estimation grew. When the spherical harmonic coefficients increased from degree 1 to 40, the density of the single mutation (tetrahedron No. 1426) was estimated from $6.833\mathrm{g}/{\mathrm{cm}}^3$ to $9.984\mathrm{g}/{\mathrm{cm}}^3$, moving closer to the assumed density. Meanwhile, the density of tetrahedron No. 35 with the double mutations increased from $6.815\mathrm{g}/{\mathrm{cm}}^3$ to $9.967\mathrm{g}/{\mathrm{cm}}^3$, and the estimated error of ten mutations also increased as the degree became larger.

Comparing the results in Figure 9 with the prior density distributions in Figure 4, it can be found that, when the least-squares method is applied, the density of the mutational tetrahedron can be well estimated. However, we can also see that the densities of the surrounding tetrahedrons will be significantly affected by the mutant tetrahedron, while the tetrahedrons far away from this mutant tetrahedron are almost unaffected. The maximum density of the remaining tetrahedrons is bigger than tetrahedron No. 35 when the harmonic degree is less than 10, as shown in Figure 7. The reason for this is that, due to the close distance between the mutation and its adjacent tetrahedron, it is difficult to distinguish the mutation during the density estimation process. Furthermore, as the degree of the harmonic coefficients is lower, the “measured” gravity field contains less information. Thus, the estimated densities of some tetrahedrons adjacent to the mutation are similar to those of the corresponding mutant, which could be larger than the other tetrahedrons. In addition, to balance the values of the spherical harmonic coefficients of the whole asteroid, some tetrahedrons near the mutation have density values that should be smaller than the other tetrahedrons.

In particular, from the results of the ten mutations case shown in Figure 8 and Figure 9, we can find that the estimated densities of the tetrahedrons, which are closer to the ends of the asteroid, are more accurate than the others, such as tetrahedrons No. 2/35/1401 and 1426. The densities of the tetrahedrons, which are close to the center of the asteroid, such as No. 624/675/790, were estimated to be no more than $8\mathrm{g}/{\mathrm{cm}}^3$—that is, only about $1\mathrm{g}/{\mathrm{cm}}^3$ larger than the mean of the remaining tetrahedrons without density mutations.

The above results can be explained by the appearance of the dumbbell-shaped asteroid Kleopatra. A tetrahedron close to the center of an asteroid could be called a near-center tetrahedron; on the contrary, the tetrahedrons distributed at both ends could be called distal tetrahedrons. Since the size of the tetrahedral element is relatively uniform in the model used in this study, the gravitational perturbation caused by the tetrahedron around the near-center tetrahedron is very close to it. Furthermore, it is difficult to distinguish whether the difference in the gravitational potential caused by the density mutation is caused by the mutational tetrahedron or its adjacent tetrahedrons. At the same time, this also explains that, even when the degree of the spherical harmonic coefficients is 40, the density of each tetrahedron cannot be solved accurately. The dimension of the equations is ${(40+1)}^2-1$, which is more than the number of tetrahedral elements, at 1467. Due to the existence of near-centered tetrahedrons, Equation (23) is easily linearly related. Therefore, a positive definite system and a unique equation solution cannot be obtained.

The results show that the spherical harmonic coefficients can be used to estimate the density distribution of asteroids. When mutations of density exist inside an asteroid, they can be detected. If the mutation is in a distal tetrahedron, the density can be estimated accurately, while for a near-center tetrahedron, the density anomaly near it can be estimated.

4.2. Density Distribution Estimation for Known Structures

If a structural hypothesis of asteroids can be obtained through other detection methods, estimating the densities of the known structures will be easier than assuming each tetrahedron is an unknown variable. In this section, it is assumed that we are exploring the asteroid 216 Kleopatra and its internal structure has been ascertained. It is, therefore, necessary to estimate the density of each constituent structure.

We refer to the one- and two-core models proposed in [23] and combine them to form a new structure, as shown in Figure 10. Assuming sphere structures with a radius of 10 km at the center and both ends of the target asteroid, the density near the center of mass is $10\mathrm{g}/{\mathrm{cm}}^3$, the density of the end along the positive x-axis is $8\mathrm{g}/{\mathrm{cm}}^3$, and the density of the end along the negative x-axis is 0. Then, the entire asteroid model is divided into four parts with different densities. It is important to state that, when dividing the four parts, they are not strict spheres. Judging the property of the tetrahedron is based on whether the center of the tetrahedron element is within the spherical range of the given structure. For this process, fine meshing can be used to obtain more precise structural divisions in the future.

Since there are only four unknowns, the system is positive definite, and all the densities can be solved accurately using the appropriate degree and order of spherical harmonic coefficients. The errors of the four parts are shown in Figure 11. When the degree is more than 2, the errors can reach the order of ${10}^{-13}\mathrm{g}/{\mathrm{cm}}^3$. However, when the degree of the coefficient is 1, we cannot obtain a precise result. As the dimension of Equation (21) is three, it is less than the number of unknowns (four), leading to an under-determined system.

The density data are not listed, but the density distribution of the hypothetical structure in the model is given, as shown in the cross-sectional view of Figure 12 (different colors correspond to different density values). Therefore, as long as the true spherical harmonic coefficients are obtained during the exploration, the density distribution of the asteroid with a known structure can be effectively estimated.

5. Conclusions

This study proposed a free-vertex tetrahedral finite element representation, which is intended to study the internal density distribution of asteroids. The explicit formulation of the free-vertex representation was obtained, based upon which we derived its transformation to the spherical harmonics, enabling us to estimate the internal structure of a massive body from its external gravitational field. The performance of the density distribution was evaluated with the asteroid 216 Kleopatra as an example, showing that: (1) under the condition of only the gravity field and shape data of the asteroid being available, the free-vertex tetrahedral FEM can be used to randomly model the interior of the asteroid, and the density mutation inside the asteroid can be detected; and (2) when otherwise known for the internal structure of the asteroid, the free-vertex tetrahedral FEM can be associated with the corresponding structure, and the density can be estimated accurately. Although the current results exhibited the potential of this free-vertex tetrahedral FEM for estimating the internal structure of an asteroid, we also remarked that the accuracy of the density mutation estimate is related to the location. For the benchmarking test, the density estimation becomes more accurate as the mutations approach the two ends of the asteroid and the density of the near-center meshes is difficult to distinguish using our method.