The spherical harmonic gravity field is widely adopted with coefficients associated with Legendre functions of degree n and order m (Kaula, 1966). U r , ϕ , λ = G M r 1 + ∑ n = 1 ∞ r 0 r n ∑ m = 0 n P n m sin ⁡ ϕ × C n m ⁡ cos m λ + S n m ⁡ sin m λ (1) where r is a spacecraft’s radial position, ϕ is its latitude, λ is the longitude, G is the gravitational constant, M is an asteroid’s total mass, r 0 is the reference radius, P n m is the normalized associated Legendre function of degree n and order m , and C n m and S n m are the normalized spherical harmonic coefficients. C n m = 2 δ l − m ! 2 l + 1 l + m ! 1 M r 0 n ∫ ∫ ∫ Vol r ′ n P n m sin ϕ ′ cos m λ ′ d M (2) S n m = 2 δ l − m ! 2 l + 1 l + m ! 1 M r 0 n ∫ ∫ ∫ Vol r ′ n P n m sin ϕ ′ sin m λ ′ d M (3) where r ′ , ϕ ′ , λ ′ are the radius, latitude, and longitude of a distributed mass inside an asteroid. δ = 0 if m = 0 ; δ = 1 if m ≠ 0 . r 0 is the same reference radius in Eq. 1. As we discuss the gravity on the surface and inside or outside a cube, it is important to note that the reference radius can be chosen arbitrarily, as usual, but here it is chosen as the mean radius of the asteroid—that is, the radius of a corresponding sphere with same mass as the asteroid. This is closely related to the convergence and divergence regions of the harmonic expansion method.

Spherical harmonic expansion can only be used to calculate gravity potential directly outside a celestial body, but it is not suited to evaluating gravity on the surface and inside, or even near but outside, a celestial body. Examples are given below. See Figure 1 as an example, which shows the gravity potential on a sphere outside a cube. Note that the potential unit in this paper can be normalized as any value. Here we pay more attention to the relative value—the patterns of the gravity.

For a polyhedral gravity field, the potential of a spacecraft at position r is (Werner and Scheeres, 1997): U r = 1 2 G σ ∑ e ∈ e d g e s r e ⋅ E e ⋅ r e ⋅ L e − ∑ e ∈ f a c e s r f ⋅ F f ⋅ r f ⋅ ω f (4) where density σ is assumed to be constant. Suffixes e and f denote the edge and face of the polyhedron. The face’s outward-pointing normal vector is n ̂ f , and its dyad F ̂ f = n ̂ f n ̂ f and n ̂ e f is outward edge normal perpendicular to both n ̂ f and the edge. For the edge connecting vertices i and j shared by face A and B, the edge dyad is E ̂ i j = n ̂ A n ̂ i j A + n ̂ B n ̂ j i B , with other similar definitions of E e s. If r i is the vector form spacecraft to polyhedron vertex P i , r i is its magnitude. e i j is the length between vertices P i and P j ; an edge factor L e is defined L e = ∫ e 1 r = ∫ P i P j 1 r d s = ln r i + r j + e i j r i + r j − e i j . (5) ω f is a dimensionless face factor bounded by the triangular face f : ω f = ∫ ∫ triangle Δ z r 3 d S = 2 ⁡ arctan r i ⋅ r j × r k r i r j r k + r i r j ⋅ r k + r j r k ⋅ r i + r k r i ⋅ r j . (6)

as the spacecraft’s gravitational acceleration is generally no longer on the radius directions for a non-spherical central body. In these calculations, some transformations are needed from spherical coordinates to Cartesian coordinates: tan ⁡ θ = y / x , sin ⁡ ϕ = z / r , r = x 2 + y 2 + z 2 .

The polyhedral method can calculate gravity outside a celestial body, and it can also be used to calculate gravity near the surface on a celestial body; however, here we try to apply it to calculate the gravity on and inside the surface of a celestial body without much modification of its algorithm. This is one of the major contributions of this paper. Many examples and analysis will be given below.

Due to the symmetry in a rectangle’s shape, the gravitational coefficients of its spherical harmonic expansion taken about its center of mass have a relatively simple form. First of all, S l m coefficients are identically equal to zero. Secondly, all coefficients C l m are such that either l or m are odd and are also equal to zero. Thus, the only non-zero gravity coefficients are those of the form C 2 l , 2 m , l , m = 0,1,2 , … . Some explicit formula for the rectangle’s gravitational field coefficients are derived here, and the first few with the dimension 2a × 2b × 2c terms without normalization are specified as: C 20 = 1 6 r 0 2 2 c 2 − a 2 − b 2 (7) C 22 = 1 r 0 2 a 2 − b 2 (8) C 40 = 1 15 r 0 4 9 a 4 + 9 b 4 + 24 c 4 + 10 a 2 b 2 − 40 a 2 c 2 − 40 b 2 c 2 (9) C 42 = 1 2 r 0 4 a 2 − b 2 10 c 2 − 3 a 2 − 3 b 2 (10) C 44 = 7 r 0 4 3 a 2 − b 2 a 2 − 3 b 2 (11) We can find that there is a similarity of the gravitational coefficients between a rectangle and an ellipsoid (Balmino, 1994; Scheeres, 2012), especially about C 20 = ( 2 γ 2 − ( α 2 + β 2 ) ) / ( 10 r 0 2 ) , C 22 = ( α 2 − β 2 ) / ( 20 r 0 2 ) in the form of its semi-axis ( α ≥ β ≥ γ ) . But other coefficients are different, there being no simple relation like C 40 = 15 / 7 ( C 20 2 + 2 C 22 2 ) , C 42 = 5 / 7 C 20 C 22 , C 44 = 5 / 28 C 22 2 for an ellipsoid. It is obvious that for a cube ( a = b = c ) , C 20 = C 22 = C 42 = 0 .

One simple way to calculate the spherical harmonic coefficients is to evaluate the gravity potential by Eq. 1 on a sphere outside the irregular central body using the polyhedral method. Then, by a usual least square algorithm, the coefficients can be obtained. See the compact form of the linear Eq. 12 about the coefficients. Since the r in Eq. 1 is here a constant on a sphere, the calculation is simplified. This is one of the important reasons that we use spheres as calculating examples here–see the potential example on an exterior sphere in Figure 1. The larger potential (yellow) parts are directly above the center of the planes, and smaller potential (blue) parts are above the corners of the cube. The patterns show that the physical meaning is quite clear, which can justify the algorithm and the code to some degree. This is one reason why a cube has been chosen first for calculating examples in this paper. One drawback of this method is that outliers may appear in the results. For example, S 32 = − 0.1777 is not trivial when calculating 6 × 6 coefficients, but S 54 = 2.0173 appears when calculating 8 × 8 coefficients, and other coefficients are the same as in Table 1. It is unexpected that only one outlier exists for each of these two cases; however, because of this, it is easy to remove the outlier from the data. Luckily, the drawback of the results by this algorithm gives us a special strategy to remove the outlier. … U r , ϕ i , λ i G M − 1 r … = … P 10 sin ϕ i … … P 11 sin ϕ i cos λ i … … P 11 sin ϕ i sin λ i … … P 20 sin ϕ i sin ϕ i … … P 21 sin ϕ i cos λ i … … P 21 sin ϕ i sin λ i … … P 22 sin ϕ i cos 2 λ i … … P 22 sin ϕ i sin 2 λ i … … … … … P n m sin ϕ i cos m λ i … … P n m sin ϕ i sin m λ i … T r 0 r 2 C 10 r 0 r 2 C 11 r 0 r 2 S 11 r 0 2 r 3 C 20 r 0 2 r 3 C 21 r 0 2 r 3 S 21 r 0 2 r 3 C 22 r 0 2 r 3 S 22 … r 0 n r n + 1 C n m r 0 n r n + 1 S n m (12)

The normalized spherical harmonic coefficient ( 6 × 6 ) , C n m , shown in Table 1 together for comparison, are S n m from 0 to the order 6 by calculating a cube with a dimension 100 × 100 × 100 length scale, and rectangles with the dimension 100 × 200 × 200 and 100 × 200 × 300 , and also an ellipsoid whose semi-axes are 100 × 200 × 300 . It is interesting to see that only C 00 , C 40 , C 44 , C 60 , C 64 are not negligible for a cube ( 200 × 200 × 200), C 20 appear also for a rectangle ( 200 × 400 × 400), and C 22 , C 42 , C 62 appear also for a rectangle ( 200 × 400 × 600), which can be validated by Eqs 7–11. It is evident that the number of a cube’s nontrivial coefficients is minimum except for a sphere. On the other hand, it also reveals the physical meaning of C 40 , C 44 , C 60 , C 64 , etc., which represent the shape of a cube. Interested readers can also see a comparison of these coefficients to those of an ellipsoid with semi-axis ( 100 × 200 × 300) length scale.

Figures 2 and 3 visualize the total polyhedral gravity magnitudes on the surfaces of a cube and a rectangle. In Figure 3, one unexpected result is that the maximum gravity appears at a center of the mid-planes with the dimension 200 × 600 . This conclusion differs from the well-known result for an ellipsoid, whose maximum surface gravity appears at the end of its short axis while its minimum surface gravity is at the ends of its long axis (Guibout and Scheeres, 2003; Dobrovolskis, 2019). Figure 4 shows gravity on a cube in its surface direction, which could be helpful for understanding dust motion on surface. It is evident that near the edge, gravity is high but that it is zero at the center of the plane. The careful reader can find the largest tangential gravities appearing near the center of the edge. This gives some clue that a natural cube is rare in the celestial realm. All dust has a tendency to move to its surface center and accumulate to become a sphere. Some trajectory analysis is given in Section 4.2.

Attention can be given to the special sphere with the radius R = 141.4 of a cube with the dimension 200 × 200 × 200. In Figure 5, part of the sphere is tangential to the edge, and part of it is inside the corner or outside the plane center of the cube. Figure 6 shows the 2D map of the gravity of the cube on a sphere whose radius R = 141.4 corresponds to the 3D Figure 5. Special attention should be given to the maxim gravity at the boundary line, which is both on the surface of the sphere and the cube. In this figure, the larger value appears at the middle of each edge of the cube, tangential to the sphere. The yellow boundary around the corners of the cube is also on its surface. Inside the boundary near the corners, the gravity is inside the corners of the cube, which is smaller than the corresponding boundary with same radius on the surface. This is physically correct. This also shows that the polyhedron algorithm works well for all these different cases. Hence, we visualize the gravity inside the cube, which could be helpful for understanding internal structure and the evolution of a celestial body. This sphere includes parts both inside the corner and outside the central plane of the cube. It is obvious for the outside parts of the sphere above the planes of the cube that its gravity is also smaller than the gravity on the boundary both on the surface of the cube and of the sphere. This is good example of how the polyhedral method is not only valid for the gravitational calculation outside a solid body but also works well for the place inside or on its surface; relatively very few researchers have done it this way. The harmonic expansion method only works for the spheres far enough from its central body. Again, we note that the acceleration unit in this paper can be normalized as any value. For example, if the length unit is kilometers, density is about 1g/ cm 3 and the acceleration unit is mm/ s 2 . Here, we pay more attention to the relative value–the acceleration pattern and distribution about a cube.

For this, see Figures 7–9. For even larger spheres outside the cube, their overall gravitational patterns are all similar. Interested readers can compare Figure 7 to Figure 1, in that the place (above corner of the cube) on the sphere of larger gravity corresponds to the place with smaller potential, whereas the place (above the center of a plane) of smaller gravity corresponds to the place with larger potential on the sphere. Thus, the place on the sphere with larger potential is at the same place with smaller gravity, and vice versa is compatible with the general definitions of potential and gravity. Particularly in Figure 9, the overall error is acceptable, and largest errors appear near the corner of the cube at about 10%. The error in the large yellow part is very small. The accuracies can be improved even further by adding more spherical harmonic terms in the model, such as C 60 , C 64 , C 80 , C 84 , C 88 .

The polyhedral method is again applied to investigate the gravity inside a celestial body, which is closely related to the internal stress and strain. This kind of research is helpful for analyzing the internal structure and evolution of an asteroid. The gravity on a sphere r = 100 , which is just inside the cube, is calculated. See Figure 10 and compare with Figure 2 ,which is similar in overall pattern; the maximum value appears just at the center of each plane of the cube. For other even smaller spheres, the gravitational configurations on them are all similar. However, the errors calculated by spherical harmonics are all too large, especially near the corners (Figures 11, 12). We see that the errors are totally unacceptable for the spherical harmonic expansion method, especially near the corners of the cube. The relativity is up to 200–300%. Interior spherical harmonic coefficients are needed (Takahashi et al., 2013). With this situation, some special modifications for this method should be given (Takahashi et al., 2013; Scheeres, 2012; Pearl et al., 2018; Arora and Russell, 2012).

With the polyhedral model, spacecraft landing and ascending trajectories can be calculated, but more interesting theoretical results can be obtained if secular motion is studied. The secular motion of periapsis in the orbital tori is similar to the well-known Mercury apsidal secular perturbation motion around the Sun. See Figures 13 and 14 and consider an initial circular orbit e = 0 in the equatorial plane i = 0 ; it can be seen that the semi-major axis a , eccentricity e , and inclination i are constant on average, the longitudinal of the ascending node Ω increases, and the argument of periapsis ω decreases by averaging the true anomaly ν . The Ω ’s rate is just twice the rate of ω , which can be explained by Eqs 13;–;16; the corresponding curves are on Figure 14. The orbit moves 16 periods, whereas the longitude processes only one period in about 350 time units. This can be derived without difficulty for the averaged equation of six elements using classical averaging theory (Kaula, 1966): a ̇ = e ̇ = i ̇ = 0 (13) Ω ̇ = − 35 C 40 n ⁡ cos ⁡ i 8 p 4 6 / 7 + 9 / 7 e 2 − 3 / 2 + 9 / 4 e 2 sin 2 ⁡ i , (14) ω ̇ = 35 C 40 n 8 p 4 12 / 7 + 27 / 14 e 2 − 93 / 14 + 27 / 4 e 2 sin 2 ⁡ i (15) + sin 4 ⁡ i 21 / 4 + 81 / 16 e 2 .

Short period oscillation is mainly due to the term C 44 , C 64 . In Figure 14, from top to bottom, the calculated classical six elements for the orbit in Figure 13 are semi-major axis a , eccentricity e , inclination i , longitude of ascending node Ω , argument of periapsis ω , and true anomaly ν . The unit for angles are degrees. More attention should be paid to the trend of variation.Other researchers are also interested in the orbital motion around a cube (Liu et al., 2011), which shows a few families of periodic orbits.

Figure 15 shows three examples of motion: a particle moving without friction in a straight line (green), a circle (red), and a general tori trajectory (blue) on a surface of a cube via polyhedral gravity calculations. When there is no initial lateral velocity, the dust or a cart’s motion is on a straight line (green line), like a spring’s harmonic motion. With the proper choice of lateral velocity, the motion can be in a periodic trajectory, similar to a circle (red circle). A random velocity normally leads to a tori trajectory (blue trajectory). The only forces acting on the particle are the gravity from the cube and the supporting force out of the surface. No friction is assumed here.

For a cube with dimension 2 × 2 × 2 unit length, the motion is simple with initial condition x 0 = 0.9 , z 0 = 0 in the plane y 0 = 1 . Even for this simple case, a complete proof is also tedious as too many terms should be included. A simplified analysis is here given. In gravity, these kinds of terms appear by simple formulation. r i = x 2 + 2 x + 6 (16) r j = x 2 − 2 x + 6 (17) k = − r i + r j r i r j r i + r j 2 − 4 (18) b = − r i − r j r i r j r i + r j 2 − 4 (19) x ̈ = k x + b (20) Where the length of the cube is 2, r i , r j is the distance from the attracted mass point to the corners of the adjacent corners i , j of the related edge. It is easy to show that | b / k | ≪ 1 (less than 10%) and k is almost a constant when − 1 < x < 1 . It is only near the edge that a little nonlinear property appears. The motion is then similar to mass spring harmonic motion. These accelerations can also be seen in the arrow representations in Figure 3 for both direction and magnitude.

When friction exists, the tori and circle motion become spiral. The linear oscillation’s amplitude will dwindle. Finally, all the motions will be stopped at the center of the plane. This could be an explanation of why in the natural world a pure cube celestial body with sharp corners is rare.