Extended two-body problem for rotating rigid bodies

A new technique that utilizes surface integrals to find the force, torque and potential energy between two non-spherical, rigid bodies is presented. The method is relatively fast, and allows us to solve the full rigid two-body problem for pairs of spheroids and ellipsoids with 12 degrees of freedom. We demonstrate the method with two dimensionless test scenarios, one where tumbling motion develops, and one where the motion of the bodies resemble spinning tops. We also test the method on the asteroid binary (66391) 1999 KW4, where both components are modelled either as spheroids or ellipsoids. The two different shape models have negligible effects on the eccentricity and semi-major axis, but have a larger impact on the angular velocity along the $z$-direction. In all cases, energy and total angular momentum is conserved, and the simulation accuracy is kept at the machine accuracy level.


Introduction
Asteroids are remnants from the formation of the solar system. It is hoped that detailed study of asteroids may improve our understanding on how the solar system was formed. One way to obtain a better knowledge of asteroids is to study the dynamics of multi-body systems. Roughly two decades ago, the existence of asteroid satellites was still uncertain (Weidenschilling et al, 1989). Not long after, the first asteroid binary, (243) Ida and its satellite Dactyl, was discovered by the Galileo spacecraft mission in 1993 (Chapman et al, 1995). Since then, many more binaries have been observed and it is believed that roughly 15% of near-Earth asteroids, larger than 200 meters in diameter, are binaries . A few asteroid binaries have also been observed amongst the main-belt asteroids, but the frequency of such systems is generally lower (Merline et al, 2002).
Studying the dynamical evolution of two rigid bodies, including both their translational and rotational motion, is known as the full two-body problem. In the full two-body problem, the rotational and translational motion of both bodies are fully coupled (Maciejewski, 1995). However, studying the full two-body problem is not trivial for irregular bodies, such as asteroids, as no analytical solution exists to compute the gravitational force between two non-spherical objects. Therefore, several numerical approaches have been developed to solve the full-two body problem.
A straightforward method of modelling the gravitational field of an asteroid is the mascon model (Geissler et al, 1996). In this approach, a set of point masses are placed in a grid, forming the shape of the asteroid. The total gravitational field of the asteroid will then be the sum of all the point mass fields. The advantage of the mascon model is that it can produce an accurate shape representation of an arbitrary body. Despite this advantage, the mascon model also has several drawbacks. For instance, the accuracy of the model depends on the number of point masses included in the body, and including a large number of point masses is computationally expensive. Furthermore Werner and Scheeres (1997) have shown that, despite including a large number of point masses, there exists significant errors in the force computation. This is due to the errors in the resolution of the surface of the asteroid, as the topology of the body is replaced with spherical balls. Wittick and Russell (2017) have revisited and optimized the model to make it less computationally demanding and increased the numerical accuracy.
A different, yet common method to model the gravitational field of non-spherical bodies is through the application of series expansions. The most common case is generated by using expansions of spherical harmonics (Konopliv et al, 2011). However, for spherical harmonics, the potential can only be computed outside a given sphere, known as the Brillouin sphere, as the potential will diverge inside this region (Moritz, 1980). Spherical harmonics are also utilized in the full twobody problem to describe the mutual gravitational potential (Scheeres et al, 1996;Boldrin et al, 2016;Boué, 2017).
Alternative parametrizations are also used to mitigate the divergence problem present in the spherical harmonics approach. For example, ellipsoidal harmonics has been used to compute the gravity field on various small solar system bodies (Garmier and Barriot, 2001;Garmier et al, 2002;Dechambre and Scheeres, 2002;Reimond and Baur, 2016). However, the use of ellipsoidal harmonics can be cumbersome, due to the mathematical and numerical complexity. Alternatively, prolate spheroidal harmonics (Fukushima, 2014;Sebera et al, 2016) provides simpler mathematical expressions compared to ellipsoidal harmonics, while providing a good geometric fit for non-spherical objects. Furthermore, Reimond and Baur (2016) have shown that, although ellipsoidal harmonics are far more accurate than spherical harmonics, there are virtually no difference between prolate spheroidal and ellipsoidal harmonics. Nevertheless, because all these methods use mathematical expansions, there exists a limit to the number of terms included in the model. Neglecting higher order terms, therefore, results in truncation errors.
Another way to expand potentials is to expand the inertia integrals. Scheeres (2009) used this approach to study the stability of two ellipsoids restricted in the plane. Hou et al (2017) present a fast method to compute the mutual potential, also using inertia integrals, between two arbitrary bodies using recurrence relations. Compère and Lemaître (2014) utilize the STF tensor formalism, which allows one to determine the coupling between spherical harmonics in a compact way, to determine the mutual potential and applied the method to study the (66391) 1999 KW4 binary system.
Another approach is to model an asteroid as a polyhedron of constant density. Similar to the mascon model, the polyhedron model allows one to include finer geometric details of an asteroid. Werner and Scheeres (1997) present a method to compute the gravitational potential of a polyhedron and use it to model the gravity field of the asteroid (4769) Castalia. Conway (2015) gives an alternative formulation of the gravitational potential of a polyhedron through the use of vector potentials, a method used to study how the dust is emitted and transported around Comet 67P . For two bodies, a method was presented by Werner and Scheeres (2005) to determine the mutual potential between two polyhedra. This approach has been used to study the dynamical evolution of the asteroid binary (66391) 1999 KW4 Fahnestock and Scheeres, 2006). Shi et al (2017) present a different model for the mutual potential between a polyhedron and a rigid body of an arbitrary mass distribution and use this method to also study the (66391) 1999 KW4 system.
For most asteroids, detailed shape models are not available. An alternative to the polyhedron model is to approximate an asteroid with a well-defined shape, such as an ellipsoid. The gravitational potential of such bodies can be expressed analytically (MacMillan, 1930). However, the potential of an ellipsoid requires computation of elliptic functions, which may be computationally demanding. Nevertheless, using ellipsoid shape approximation to model asteroids has been e.g. used to study an ellipsoid-sphere system (Scheeres, 2004;Bellerose and Scheeres, 2008).
A new approach, based on vector potentials, to compute the gravitational potential between two extended bodies is presented by Conway (2016). Here, the force, torque and mutual potential energy are formulated as surface integrals under the assumption of constant density. This could potentially be a faster and more accurate technique, as the force and torque integrals are converted from volume integrals to surface integrals. Moreover, this approach does not rely on series expansions and will therefore not suffer from truncation errors. This approach was recently tested on coplanar spheroids and thin disks by Wold and Conway (2021). In this paper, we extend the work of Wold & Conway by giving a full three-dimensional treatment of the dynamics between the two bodies, including the coupling between translational and rotational motion.
In sec. 2, we introduce the equations used to determine the translational and rotational motion and describe how we treat the rotational kinematics. In Sec. 3, we present the results of two dimensionless test scenarios, while in Sec. 4 we apply our method on the asteroid binary system (66391) 1999 KW4. In this paper, all dotted variables denote the time derivative of the corresponding term. All vectors with a hat (e.g.x) show the vector expressions in the body-fixed frame (local frame), while those without a hat are in an inertial frame.

Mathematical model
The main objective of the method is to calculate the force F, torque M and potential U of a solid body of uniform density ρ in a gravitational field g (r) in terms of surface integrals over the surface of the body.
The force of an extended body with constant density ρ, due to an external gravitational field g(r), can be expressed through a volume integral where r is a position vector to a point on the surface of the body being integrated over. If the gravitational field is given by the gradient of a scalar potential Φ(r), so that g(r) = ∇Φ(r), we can apply the divergence theorem such that the force can be computed through a surface integral It can also be shown that the torque about the mass center of an extended body is given byM and the mutual potential as (Conway, 2016). In these formulas n is a unit normal to the body surface. The above surface integrals are thus alternative expressions for the force, torque and potential between two bodies. In the following when we consider two interacting bodies, we replace the scalar potential Φ(r) with known analytical formulae for spheroids and ellipsoids as the first body, and integrate over either a spheroid or an ellipsoid as the second body. The surface integration method is outlined in detail by Wold and Conway (2021). Analytical formulas for g(r) and Φ(r) in terms of elementary functions are given for spheroids and triaxial ellipsoids in (MacMillan, 1930). The Macmillan formulas for g(r) were validated using Eq. (2) by taking Φ(r) to be the scalar potential of a unit point mass and integrating this over a spheroid. Exact agreement was obtained between Eq. (2) and the MacMillan formulas. The analytical expressions for the moments of inertia of spheroids and ellipsoids are already available. Therefore, it is straightforward to set up surface integration schemes in a local body-fixed coordinate system for spheroids and ellipsoids.
Here, we analyze the motion of the bodies in a global inertial coordinate system using the forces and torques acting on each body and applying Newton's laws. As we consider a two-body problem, if F A and F B are the global forces acting on the two bodies A and B, then from Newton's third law It should be noted thatM A , which is the torque ofF A about the origin of the body-fixed frame of body A, may not necessarily be equal toM B , which is the torque ofF B about the origin of the body-fixed frame of body B. However, the torque of all forces about any arbitrary predefined point in the system must be zero in order to maintain the conservation of the angular momentum. When the forces and torques are computed, the equations of motion can be solved as a standard initial value problem, where the velocities v and positions r can be integrated as and the angular velocityω integrated from where J is the angular momentum of the body. We use the embedded Runge-Kutta method of order 9(8) by Verner (Verner, 2010) to solve the equations of motion. The coefficients of the Butcher tableau are taken from Verner's web page 1 , using the 'most efficient' coefficients provided in the web page.

Rotation angles and rotational motion
The force and torque integrals given in Eqs. (2) and (3) are computed in the body-fixed frames. However, because we are interested in the motion in the inertial frame, the equations of motion must be projected back to the inertial frame. To move between different reference frames, we use a rotation matrix. In this paper, we adopt the Tait-Bryan convention where the rotation matrix R takes the form The angles (φ, θ, ψ) are the Tait-Bryan angles, and correspond to rotations around the body-fixed x, y, z axes respectively 2 . The rotation matrix given in Eq. (12) is then used to project the force in the body-fixed frame, from Eq. (2), back to the inertial frame using the equation The computed torque in Eq.
(3) is used to determine how the angular velocity changes over time, using the equations of motion for each body in its body-fixed frame, given as (Curtis, 2013), where (ω x ,ω y ,ω z ) are the angular velocity components in the body-fixed frame, I is the moment of inertia tensor and (M x ,M y ,M z ) are the components of the torque computed in the body-fixed frame. The angular velocity (ω x ,ω y ,ω z ) is related to the angular velocity in the inertial frame as The rotation angles (φ, θ, ψ) change over time, and the following kinematic equations relate the time rate of change of these angles to the angular velocity of the bodyφ (Fossen, 2011). These kinematic equations become singular when θ = nπ/2 for any odd integer n. This singularity is a mathematical problem and is normally resolved by using Euler parameters (see e.g. Kane et al, 1983). Nevertheless, this singularity is generally not a problem for our cases in question.

Surface integration
The expressions for the surface elements in the surface integrals given by Eqs.
(2)-(4), as well as the components of the moment of inertia tensor, vary depending on the shape of the bodies. This can be generalized to any arbitrary ellipsoidal shapes. Consider a general ellipsoid with semiaxes (a, b, c). The surface elements ndS and n × rdS, used in Eqs. (2)-(4) respectively, for a general ellipsoid, are given as (Wold and Conway, 2021), where α is the angle of the cylindrical coordinates (written as φ in Wold & Conway) and z is the generalized latitude line on the Fig. 1 An illustration on how the gravitational field is computed after the bodies have rotated. The gravitational potential Φ B (blue dotted lines) is rotated according to the rotation of spheroid B. The axes x, y, z denotes the inertial frame of the system. ellipsoid. For a general ellipsoid with density ρ, the only non-zero components of the moment of inertia tensor I ij can be expressed analytically as

Rotated gravitational potential
In this section, we will describe how to incorporate the change in the gravitational potential Φ(r) and gravitational field g(r) when the bodies are rotated. Consider two bodies A and B whose mass centers are respectively positioned by r A and r B from the origin of the inertial reference frame and with their respective rotation matrices R A and R B . The vector linking the two centroids therefore becomes Let us consider the case where we compute the force on body A. In the inertial frame, the position vector of an arbitrary point located on the surface of A is expressed as wherer p,A is a vector describing the surface points of A in its body-fixed frame. For a general ellipsoid, this becomeŝ The term r used as input argument in the Φ(r) and g(r) functions is a vector that points from the centroid of B to any point on the surface of A. This can be written as which can be expressed in the body frame of reference aŝ Equation (30) is then used as the input argument of both the gravitational potential Φ(r) and the gravitational field g(r). As the bodies rotate, so must the relative gravitational field g(r). This can be achieved using the following equation

Total energy
When all the differential equations are solved, we use the computed velocities, positions, angular velocities and rotation angles to determine the energy of the system. The translational kinetic energy of each body is computed as where m and v are the mass and the center of mass velocity of the body, respectively. For an ellipsoid with constant density, the kinetic energy of rotational motion is The mutual potential energy U is computed by Eq. (4). The total energy of the system is the sum of the kinetic, rotational and potential energies.

Angular momentum
Because we do not include external forces and moments to model the system, the angular momentum and total energy of the system must be conserved. The total angular momentum, in the inertial frame, is where the angular momentum of each body is given by As previously stated, for a general ellipsoid, the only non-zero components of the moment of inertia are I 11 , I 22 and I 33 . Therefore, the total angular momentum becomes

Gravitational potential
To compute the force and torque given in Eqs. (2) and (3), we require an expression for the gravitational potential Φ(r). This potential can take different forms depending on the shape of the body. In this paper, we consider spheroids and triaxial ellipsoids, for which analytical expressions for Φ(r) are available. For a general ellipsoid with semiaxes a > b > c and constant density ρ, the gravitational potential is given by , 1930), where F (ω κ , k) and E(ω κ , k) are the elliptic integrals of the first and second kind respectively, κ is the largest root of the equation and Table 1 The initial parameters used for the two dimensionless test cases. The second and third columns show initial positions r 0 and velocities v 0 in the inertial frame, respectively. The fourth column indicates the initial rotation angles in radians. The fifth column shows the initial angular velocities in the body-fixed frames, in units of radians per unit time. The sixth and seventh columns show the semiaxes (a, b, c) of the spheroids and their masses, respectively. The top and bottom rows, for each case, correspond to body A and B respectively.
For an oblate spheroid with semiaxes a = b > c, the gravitational potential can be expressed as (MacMillan, 1930), where κ still satisfies Eq. (38). It should be noted that the gravitational potentials given in Eqs. (37) and (41) are exterior potentials. The expression for the gravitational field g(r), which is required to compute the mutual potential energy in Eq. (4), is derived in Appendix A.

Dimensionless test scenarios
We first test our method for two dimensionless cases, where the gravitational constant is set to unity, i.e. G = 1. In these scenarios, we let both bodies take spheroidal shapes. The initial conditions and spheroid parameters for the simulations are summarized in Tab. 1. The time span of the simulations is t ∈ [0, 200]. The actual simulation times on an ordinary laptop computer for these test scenarios varied from 20 seconds to two minutes.

Case 1: Small initial rotation of the spheroids
We first consider a case where the spheroids have a small initial rotation angle φ 0,A = −φ 0,B = π/32 about the x-axis and zero angular velocity. Their physical parameters and initial conditions are listed in Tab. 1 under "Case 1 (Rotated)". Fig. 2a shows the spheroid orbits projected into the xy-plane, while Fig. 2b shows their z-positions as functions of time. Because forces are acting parallel to the z-axis, as seen in Fig. 3, the motions of the spheroids are no longer restricted to one plane. The motion along the z-axis is 11 orders of magnitude smaller than the motions along the x and y-axes. The motion in the z-direction grows with simulation time, and becomes of the order 10 −1 at t ≈ 2200. The motion along the x and y-directions, unlike the motion in z-direction, did not have any significant changes in the same time duration.
The rotation angles (φ, θ, ψ), plotted as sine functions, for spheroids A and B, are shown in Fig. 2c and Fig. 2d, respectively. Both spheroids rotate multiple times about their x and z-axes, and therefore have tumbling motion. The angular velocities, projected onto the respective body-fixed frames, are shown in Fig. 2e and 2f for spheroids A and B. Due to the rotational symmetry of the bodies about their z-axes, no gravitational forces can change the rotational speeds about the z-axes. Hence, the body-fixed componentω z remains constant for both spheroids.
The top panel of Fig. 4 shows the energies in the system. While both spheroids have no initial rotational motion, the rotational energy starts to make a noticeable contribution near t ≈ 75, which is consistent with Figs. 2e and 2f, indicating that the angular velocities start to increase around t = 75. The middle panel of Fig.  4 shows the relative error in the total energy as a function of time. The relative error is computed as where N s is the number of data points. Because the relative error is smaller than 10 −12 , the total energy can be considered as conserved in our simulation. Finally, the absolute error of the three total angular momentum components are shown in the bottom panel of Fig. 4. The absolute error is computed as Similar to the total energy, the error in the total angular momentum is smaller than 10 −12 , and we conclude that also the total angular momentum is conserved in the simulation.

Case 2: Adding initial angular velocity to the spheroids
We now consider a system where both spheroids are initially spinning around their z-axes, withω 0,z = 2, while the other angular velocity components are initially set to zero. We also let both spheroids start with a slight tilt, in which θ 0,A = −θ 0,B = π/32, while the other angles are initially zero. The other parameters are the same as those in case 1. Fig. 5a shows the spheroid orbits projected onto the xy-plane. The orbits for this case closely resemble the case for two spheroids in co-planar motion (Wold and Conway, 2021). Fig. 5b shows their z-positions as a function of time. As in Case 1, a force along the z-axis can occur, but the motion in the z-direction is 11 orders of magnitude smaller than the motion in the x and y-directions for both spheroids.
Figs. 5c and 5d show the rotation angles φ and θ, plotted as sine functions, of bodies A and B respectively. As the spheroid spins around its z-axis, the angle ψ increases linearly and will appear as a straight line with constant slope in the figures. Meanwhile, the angles (φ,θ) oscillate between −π/32 and π/32. As ψ is increasing linearly and (φ, θ) are both oscillating, the rotational motion of the   Fig. 4 For case 1. Top: The potential energy U (green), translational kinetic energy E k (red), rotational kinetic energy Erot (orange) and the total energy Etot (blue). Middle: The relative error of the total energy. Bottom: The absolute errors in the components of the total angular momentum. The errors in both the total energy Etot as well as the total angular momentum Jtot are smaller than 10 −12 , which demonstrates that these quantities are conserved in our simulations. The angular velocity components (ω x ,ω y ) are shown in Figs. 5e and 5f for bodies A and B, respectively. Because of the rotational symmetry of the spheroid, ω z will remain constant at the valueω z = 2 throughout the simulation, and is therefore excluded from the figures. The values of theω x andω y components peaks when the two spheroids are in close proximity. By increasing the initial angles by a factor of four, i.e. θ 0 = π/8, the amplitude of both theω x andω y components nearly doubles, indicating that the amplitude of these two angular velocity components is sensitive to initial angle.
The different parts of the total energy are shown in the top panel of Fig. 6. Due to the spheroids' rotations, the rotational energy now makes a significant contribution to the total energy, being greater than the kinetic energy from the translational motion. However, becauseω x andω y are roughly 3 orders of magnitude smaller thanω z , andω z is constant, the rotational energy is nearly constant. Furthermore, both the total energy and total angular momentum remain constant in the simulation, as their respective errors are smaller than 10 −12 .

Application to (66391) 1999 KW4
In this section, we apply our method to the asteroid binary system (66391) 1999 KW4, where we refer to the primary as Alpha, and to the secondary as Beta. We will study the dynamical evolution of the bodies and see how the outcome changes when the body shapes are changed from ellipsoids to spheroids by setting the semiaxes to be a = b > c. The simulations run with a timespan of 1 year.
The physical parameters, as well as initial conditions, used in the simulations are shown in Tab. 2. We let Alpha start at rest from the origin of the inertial frame in the simulations. We also assume that Beta is located at the pericenter at The total energy and total angular momentum of the system is conserved, as demonstrated by the small relative error in these two quantities. t = 0, and thus set the mean anomaly M 0 to be initially zero.  also find that the angle between the rotation pole of Alpha and the binary orbit is between 0 • and 7.5 • with a nominal separation of 3.2 • . As such, we let Alpha be initially rotated with φ 0 = 3.2 • , whereas Beta initially remains non-rotated. It is assumed that the orbit of Beta is synchronous, and the initial angular velocity of Beta is set so that its rotation period is equal to its orbital period.
The top, middle and bottom panels of Fig. 7 shows the eccentricity, semimajor axis and inclination of Beta respectively. The inclination is calculated with respect to the original plane of orbit. For a spheroidal model, the eccentricity can take slightly higher values, where the eccentricity peaks at e = 0.0140 for the spheroidal model, while it peaks at e = 0.0126 for the ellipsoidal model. This is lower to the eccentricity in the excited state of Fahnestock and Scheeres (2008) (see also Compère and Lemaître, 2014;Hou et al, 2017;Shi et al, 2017), which find that the eccentricity can surpass e = 0.03. The range of the semi-major axis are similar in both simulation types, where it takes values between a s ∈ [2.5418, 2.5526] km and averaging atā s = 2.5471 km for the ellipsoidal model, and a s ∈ [2.5433, 2.5521] km and averages atā s = 2.5477 km for the spheroidal model.
Because the orbital period is proportional to the semi-major axis, the fact that the average semi-major axis is larger in the spheroidal simulation, indicates that the orbital period is also longer. We find that the average orbital period of Squannit isT = 17.4116 hrs for the ellipsoidal simulation andT = 17.4177 for the spheroidal simulation. The period of the inclination is longer for the ellipsoidal model, where the period is approximately 3900 hrs between the two maxima, whereas the period is approximately 3400 hrs for the spheroidal model. Fig. 8 shows the angular velocity components of the secondary for the first 200 hours of the simulation. Here, we find that the range ofω z isω z ∈ [9.55, 11.16]·10 −5 rad/s. Compared to the findings of Fahnestock and Scheeres (2008), the range range is smaller than the excited state, but also larger than the relaxed state, than that of Fahnestock and Scheeres. Another major difference is in the components of bothω x andω y . The work of Fahnestock and Scheeres find that both of these components change very insignificantly for both the excited and relaxed configurations. Our findings, however, show that the components respectively oscillates around the valuesω x ∈ [−1.89, 1.86] · 10 −5 rad/s,ω y ∈ [−3.20, 3.23] · 10 −5 rad/s. As previously mentioned in dimensionless test scenario "spinning system" (see Sec. 3.2), the amplitude of theω x andω y components could be affected by the initial angle. The difference in our result, and the one of Fahnestock and Scheeres, could therefore be due to the difference in the initial angles.
Changing the body shapes from ellipsoids to spheroids significantly affectsω z . As seen in Fig. 8, by allowing the bodies to take a spheroidal shapes,ω z for Beta becomes constant as opposed to oscillating when it had an ellipsoidal shape. This is because spheroids are rotationally symmetric, and no torques can act to change the angular velocity in the z-direction. This is also what was seen in the dimensionless test scenarios in Sec. 3. Theω x andω y components, on the other hand, still oscillates between the values seen for an ellipsoidal shape. The angular velocity of Beta in the spheroidal case is similar to the result that was previously shown in Fig. 5e and 5f, in whichω z was constant and bothω x andω y were oscillating over time. Fig. 9 shows the relative error in the total energy and the relative error of the components of the total angular momentum in the top and bottom panels respectively. The blue and red curves correspond to simulations where the bodies take ellipsoidal and spheroidal shapes. The relative error of the total energy is smaller than 10 −14 for both the ellipsoidal and spheroidal simulations. For the total angular momentum, we only show the relative error in the ellipsoidal simulation. Here, the relative errors in each component are smaller than 10 −11 . The errors in the spheroidal simulation are similar to the ellipsoidal simulation. Although, throughout the simulation, there is a drift in the total angular momentum, causing the errors in each component to increase over time.

Summary and discussion
We have validated and explored an alternative method for simulating the dynamics of a fully three dimensional rigid two-body problem suggested by Conway  (2016). The method, which is based on vector potentials instead of scalar potentials, uses surface integrals to determine the force, torque and mutual potential energy between two bodies. Wold and Conway (2021) outlined the surface integration method in detail, and tested it for a pair of coplanar spheroids and thin disks.
In this work, we extend the work and apply the method to pairs of ellipsoids and spheroids that can be randomly oriented with respect to each other, hence torques and angular momentum exchange is included. Tab. 3 shows a summary on the CPU time used for each simulation. The CPU time used for the dimensionless simulations varied between 20 seconds to two minutes. For the (66391) 1999 KW4 system, the CPU time required varied between one hour to four hours.
Two dimensionless cases were studied, where both test scenarios considered spheroidal body shapes. In the first case, both the spheroids are initially rotated around their body-fixed x-axes. Despite the lack of initial motion along the zdirection, the initial rotational tilt allowed a small force component along the x, 2 ellipsoid y, 2 ellipsoid z, 2 ellipsoid x, 2 spheroid y, 2 spheroid z, 2 spheroid z-direction to take place. This results to a small motion along the z-direction over time, although the motion is 11 orders of magnitude smaller than the motion along the x and y direction. Furthermore, the rotational motion of the spheroids develops into tumbling motion. The second dimensionless case considered is similar to the first one, where both bodies are now initially rotated around their bodyfixed y-axes. However, in this scenario, both bodies also start with an angular velocity about their body-fixed z-axes. The motions of the spheroids' in this case closely resembled two tops spinning in orbit around the common center of mass. By spinning the spheroids about their axes of symmetry, the rotational motion stabilises so that it is no longer tumble-like.
Finally, we apply the method on the asteroid binary system (66391) 1994 KW4. In this scenario, we consider two types of simulations: one where both bodies have ellipsoidal shapes and one where both have spheroidal shapes.
We compare the difference in the dynamical evolution of Beta when the bodies had ellipsoidal and spheroidal shapes. The eccentricity, on average, is larger in the spheroidal simulation. Furthermore, the values of the eccentricity is smaller than the findings of of Fahnestock and Scheeres (2008). The values of the semimajor axis are similar for both simulation types, but the average semi-major axis is slightly larger for the spheroidal simulation compared to the ellipsoidal simulation. This also indicates that the orbital period becomes longer when the bodies take ellipsoidal shapes. We also find that the time period it would take for the inclination to reach its maximum are longer when both bodies took ellipsoidal shapes.
The angular velocity components of Beta is also studied. The results are also compared to the findings of Fahnestock and Scheeres (2008), where we find that thê ω z component, for the ellipsoidal simulation, is similar to the findings of Fahnestock and Scheeres. However, the evolution of bothω x andω y are different, in which we find that these components are oscillating with larger amplitudes compared to the findings of Fahnestock and Scheeres, which is due to the difference in the initial conditions. Studying the errors in the total energy and total angular momentum serves as a check of simulation accuracy. We find that the errors for both the total energy and total angular momentum are smaller than 10 −12 for all simulations presented. The errors, which are numerical in origin, are small enough to demonstrate that our model conserves energy and angular momentum.
While our method has only been demonstrated here for a handful of scenarios, the method can also be generalized to an N -body simulation, which can be used to simulate e.g. an asteroid triple system and even include the gravitational pull from the planets in the Solar system.

A Gravitational field of an ellipsoid and spheroid
The gravitational field is required to compue the mutual potential energy in Eq. (4). We will here derive an analytical expression of g = (gx, gy, gz) = ∇Φ based on the expression in Eqs. (37) and (41). It should be noted that, while κ is a function of (x, y, z), when taking the partial derivatives of the gravitational potential Φ, κ can be treated as a constant (MacMillan, 1930).
For a general ellipsoid, the components of the gravitational field thus become For an oblate spheroid, the components of g are The value of κ, for both the ellipsoid and spheroid cases, still satisfies Eq. (38).