The planar two-body problem for spheroids and disks

We outline a new method suggested by Conway (2016) for solving the two-body problem for solid bodies of spheroidal or ellipsoidal shape. The method is based on integrating the gravitational potential of one body over the surface of the other body. When the gravitational potential can be analytically expressed (as for spheroids or ellipsoids), the gravitational force and mutual gravitational potential can be formulated as a surface integral instead of a volume integral, and solved numerically. If the two bodies are infinitely thin disks, the surface integral has an analytical solution. The method is exact as the force and mutual potential appear in closed-form expressions, and does not involve series expansions with subsequent truncation errors. In order to test the method, we solve the equations of motion in an inertial frame, and run simulations with two spheroids and two infinitely thin disks, restricted to torque-free planar motion. The resulting trajectories display precession patterns typical for non-Keplerian potentials. We follow the conservation of energy and orbital angular momentum, and also investigate how the spheroid model approaches the two cases where the surface integral can be solved analytically, i.e. for point masses and infinitely thin disks.


Introduction
In celestial mechanics, a classical problem is to model the dynamics of two rigid, extended bodies under mutual gravitational attraction. In the most general case, the bodies have arbitrary shapes and can have both translational and rotational motion, yielding twelve degrees of freedom. Naturally, to model such a system is computationally expensive, and simplifications and approximations are commonly made.
During the last 20 years, there has been renewed interest in the extended two-body problem in astronomy as binary asteroids have been discovered and studied in detail (e.g. Margot et al. 2002). The two components of a binary asteroid can be physically close and have irregular shapes, and their translational and rotational motion is coupled through energy Department of Engineering Sciences, University of Agder, Jon Lilletuns vei 9, N-4879 Grimstad Tel.: +47-37233476 E-mail: margrethe.wold@uia.no and torque transfer. To describe the dynamics of such a system requires the full two-body problem to be solved.
The mutual gravitational potential U between two extended bodies A and B of arbitrary shape is a six-dimensional integral as it involves volume integration over each body: with r being the distance between mass elements dm A and dm B . In the literature, a number of approaches are described to compute the mutual potential, largely depending on the application. It is customary to assume constant density and convert the integral to an integral of the scalar potential of one body over the volume of the other. The integrand is expanded in terms of Legendre polynomials or spherical harmonics, and higher order terms are neglected depending on the required accuracy (e.g. Borderies 1978;Hartmann 1994).
Another common approach is to use Cartesian coordinates instead of spherical coordinates, and convert the integrals over mass elements to inertia integrals (Maciejewski 1995;Ashenberg 2007;Scheeres 2009;Boué & Laskar 2009;Jacobson & Scheeres 2011;Hou, Scheeres & Xin 2017). This approach is motivated especially for non-spherical shapes and artificial satellites. By using spherical harmonics and/or inertia integrals, the gravitational potential will always be an approximation because higher order terms in the series are neglected. Truncation errors can become significant, especially close to the surface of the bodies, or can accumulate in a manner that will cause problems when describing long-term dynamical behaviour. Furthermore, for spherical harmonics, the radius of convergence for the approximating series is a circumscribing sphere around each body, and within this sphere the series does not converge (mitigated to some degree by using spheroidal or ellipsoidal harmonics (Jekeli 1988;Garmier & Barriot 2001;Fukushima 2014;Reimond & Baur 2016)).
Asteroid shapes can also be represented by polyhedra, and several works exist where binaries are modeled by two polyhedra (e.g. Werner & Scheeres 2005;Fahnestock & Scheeres 2006;Hirabayashi & Scheeres 2013). However, in these cases the mutual potential is also expressed as an infinite series where truncations become necessary for practical purposes. A somewhat different approach is presented by Shi, Wang & Xu (2017) where one asteroid is modeled as a homogeneous polyhedron with a closed-form gravitational potential thus avoiding truncation, and the other body is modeled with an arbitrary mass distribution, and the mutual potential is expressed by inertia integrals truncated at second order.
The methods listed above have in common that they all involve series expansions to describe the mutual potential and the dynamics of two extended rigid bodies. In this paper we utilize the approach used by Conway (2016) who explored the use of integral theorems to express the force, torque and mutual potential between two gravitating bodies as twodimensional surface integrals over one or both bodies. We test the application of two of the surface integral equations in a dynamical context for two extended bodies under mutual gravitational interaction. To our knowledge, this scheme has not been tried for two extended bodies before. We have chosen to concentrate on planar motion of spheroids and thin disks as they have well-known analytical expressions for the gravitational potential (MacMillan 1930), and as long as this is the case, this method also avoids series expansion. The equations of motion are solved using a Runge-Kutta integrator, and we show example orbits and address changes in orbital angular momentum and total energy as a probe of numerical accuracy.
The paper is organized as follows: In Sect. 2 we introduce the integral equations based on Conway's (2016) description and in Sect. 3 we outline a surface integration method that can be applied to spheroids and/or ellipsoids, and use it to compute the force and mutual potential for the case of two spheroids in Sect. 3.1. In Sect. 3.2 we outline how to find the force and mutual potential between two coplanar thin disks using an analytical expression with complete elliptical integrals. We thereafter run simulations with two spheroids and two disks, and present the results in Sect. 4, where we also investigate the conservation of energy and angular momentum. In Sect. 4.4 we compare pairs of spheroids having extreme shapes (nearly spherical and significantly flattened) to analytical solutions for point masses/spheres and infinitely thin disks. We discuss the results and some applications, as well as prospects for future work in Sect. 5.
2 Surface integral equations for force, torque and mutual potential For the general extended two-body problem, three integrals are needed for the force, torque and mutual potential U. The force F on a body of constant density ρ is given by the wellknown surface integral where dS is a vector perpendicular to the surface S of the body and r is a position vector to a point on the surface of the body being integrated over. This equation is immediately obtained by applying Gauss's theorem to the volume integral for the force F: where g(r) = ∇Φ(r) is the gravitational field acting on the the body at any point r. The corresponding volume integral for the torque M is Conway (2016) introduced two alternative vector potentials V 1 (r) and V 2 (r) for the vector and V 2 = 1 2 |r| 2 g(r).
Substituting these vector potentials into the well-known theorem gives two alternative surface integrals for M: A third vector potential for r × g(r) can be defined as V 3 (r) = (V 1 + V 2 )/2, which has the property The volume integral for the mutual potential is and employing ∇ · V 3 (r) in the divergence theorem immediately gives a surface integral for U as Eqs. 2, 8/9 and 12 are surface integral equations for force, torque and mutual potential, hence the integral dimension is reduced by one compared to the more common volume integrals. For the planar motion of spheroids considered here, the torque integrals are not needed, as no distribution of matter can induce a torque about an axis of rotational symmetry of a body. In a future paper we will implement torques (Ho et al. in prep.).
In this paper, we test the application of Eqs. 2 and 12 in a dynamical context for two extended bodies under mutual gravitational interaction. In the next section, we outline our method of surface integration for spheroids/ellipsoids.

Force and mutual gravitational potential: surface integration method
Assume we have two homogeneous, rigid bodies A and B with constant mass densities ρ A and ρ B and that their respective body-fixed coordinate systems are (x, y, z) and (x , y , z ), see Fig.1. Body B has gravitational potential Φ B (r) where r = [x, y, z] is a position vector in A's body-fixed frame (for now, we choose to express B's gravitational potential in A's coordinate system). The force on body A from body B is given by the surface integral: where dS is an outwardly pointing vector normal to the bounding surface S A of A. Assume that the potential Φ B (r) is known and can be evaluated at all positions on A's surface. The next step is to do the surface integration over body A. We wish to do the surface integration over A in A's body-fixed frame, and use azimuthal angle φ and z-coordinate as integration variables, see Fig.1. We assume that A has the shape of an ellipsoid with semiaxes a > b > c. The surface of A is therefore described by: The cross section of a plane parallel to the xy-plane and the ellipsoid surface is an ellipse with equation where α and β are: α = a c c 2 − z 2 (16) Fig. 1 Illustration showing our use of variables and coordinate systems. A plane parallel to the xy-plane of A intersects the ellipsoid and forms an ellipse traced out by azimuthal angle φ and radius vector r. We integrate over the surface of A, and the area element dS is spanned by the two vectors dr 1 and dr 2 (obtained by substituting the equation for a plane parallel to the xy-plane into the ellipsoid equation). Hence, lines of constant z on the ellipsoid surface can be expressed by Eq. 15, see also Fig. 1. We parametrise these ellipses in terms of azimuthal angle φ in the usual manner as: For finding the area element vector dS, we start with small displacements dr 1 and dr 2 on the ellipsoid surface along constant z and constant φ, respectively. We refer to Fig. 1 in the following. The small displacement dr 1 along a constant z is given by where we have differentiated Eq. 18 and 19. A corresponding displacement dr 2 along a line of constant φ is given by In general the line elements dr 1 and dr 2 are parallel to the sides of a parallelogram, and for a triaxial ellipsoid they are orthogonal only along certain symmetric curves. For the surface of a spheroid and a sphere, the line elements are orthogonal everywhere. The area element dS is evaluated as the cross product between dr 1 and dr 2 : With dS from Eq. 22 we can evaluate the surface integral for the force in Eq. 13 if the gravitational potential of body B is known. The force on body A can thus be found, and by Newton's third law also the force on body B: If the gravitational potential of B is a function that can be evaluated at each position r on A's surface, the force can be computed with this double integral. For the potential energy U given by Eq. 12 we need to evaluate the two terms rΦ B (r)·dS and 1 2 |r| 2 g B (r) · dS. The first term is straightforward, and to compute the second term, the gravitational field is found from the gradient of the potential, g B (r) = ∇Φ B (r). In cases where the potential of B can be expressed analytically, exact values of Φ B and g B can be computed at each point r on A's surface during the integration.

Tests for two spheroids
In this paper, we wish to test the use of Eqs. 12 and 23 in a relatively simple case where both bodies are spheroids with a = b. This allows us to use the integration scheme outlined above, and to utilize the known analytic expressions for the gravitational potential of oblate (a > c) spheroids (prolate spheroids with a < c were also tested). In order to clarify our use of symbols, we add subscripts to the spheroid shape parameters so that a A (= b A ) becomes the equatorial radius of spheroid A, and c A the distance from the centroid of A to its pole along the symmetry axis (z-axis). Similarly for spheroid B. For an oblate spheroid B, the gravitational potential at an exterior field point is given by the following analytic expression (MacMillan 1930): where κ is defined by the equation: (see MacMillan (1930) for the corresponding equation for a prolate spheroid). As mentioned above, Φ B (r) appearing in the double integral is the gravitational potential of B expressed in A's coordinate system. However, MacMillan's expression above uses r = [x , y , z ] which is the position vector of a field point in the body-fixed frame of B with B's centroid at the origin. During the integration, when calling the functional form given in Eq. 24, we therefore replace r with r + r AB , where r AB is the vector from B's to A's centroid (see Fig. 1).
For evaluating the integral for U in Eq. 12, the gravitational field has to be evaluated as g B (r) = ∇Φ(r). Differentiating Eq. 24 yields: where and Again, when calling the function g B in the surface integral, the position vector changes from r to r. By using the analytical expression for Φ B (r ) as well as the surface element dS defined in Eq. 22 in the two integral equations, we obtain expressions for the force and mutual potential that are exact in the sense that they do not involve approximations that normally follow when gravitational potentials are written as truncated series. With this method, the expressions for the force and potential energy are integral equations and therefore not analytical. The integrals have to be solved numerically, hence with "exact" we mean exact to within the limits of numerical integration. Here, we solve the surface integrals numerically by using the Gaussian quadrature integration scheme for double integrals implemented in Python (integrate.dlbquad) with the relative and absolute tolerance set to the default value of 1.49 × 10 −8 .
As far as the authors are aware, the surface integrals for gravitational force and mutual potential energy between two non-spherical bodies can be analytically expressed in just one case, which is that of two coplanar, thin, non-coaxial disks. We investigate this case in the next section.
3.2 Tests for two coplanar disks Conway (2016) has also derived the gravitational force between two rigid, coplanar, noncoaxial, infinitely thin disks with masses m A and m B and radii R A and R B . As opposed to the case with two spheroids, the surface integrals for force and mutual gravitational potential have analytic forms. The expression is in closed form in terms of complete elliptic integrals, and we also wish to test this solution. The force on one disk from the other is given by: where r AB is the distance between the axes of the disks, and f (R A , R B , r AB ) is a dimensionless shape factor representing the deviation from an inverse square law: Here K and E are the complete elliptic integrals of first and second kind, respectively. The two factors k + and k − are given by: For the two disks, the mutual gravitational potential has the following analytic expression (Eq. 151 and 152 by Conway (2016)): where u is another dimensionless shape factor: and K, E and k +/− have the same meaning as in Eq. 31 and 32. It is straightforward to balance the gravitational and centrifugal forces for rotation of two such disks in circular coplanar orbits at constant distances from the centre of mass of the disks. This seems to be the only known truly analytical solution for the motion of two extended bodies when neither body has spherical symmetry.

Two coplanar spheroids
Various tests were run with two spheroids of varying shapes and initial parameters. With the forces determined, the equations of motion for the two bodies become a system of six first-order differential equations which can be solved as a standard initial value problem. Table 1 Parameters for the two models discussed in the text, one with two identical spheroids ('2S') and one with two identical thin disks ('2D'). The mass of the bodies is m, spheroid axes are (a, c) and disk radii are R (denoted R A and R B in Eqs.31 and 33). The last column is the time step used in the simulation. We solve the equations of motion in an inertial frame using a fourth-order Runge-Kutta integrator (RK45) with a fixed time step. There are two issues with this integrator: First, Runge-Kutta integrators are known to cause a drift in energy, but as we have restricted ourselves to only shorter simulations covering a few orbits, not much energy drift is expected to occur. Secondly, by using a fixed time step integration errors will increase at periapsis, in particular for systems where the bodies are close. However, since the simulations span only a few orbits, we have chosen a small time step which will contribute toward minimizing integration errors. We include an error analysis in Sect. 4.3 by following the conservation of energy and angular momentum and find that the errors are acceptable for the demonstration purposes of this paper. In the following, we make the simplifying assumption that the two spheroids have a common equatorial plane, i.e. they share a common xy-plane. Hence, there will be no exchange of angular momentum between the two bodies and all torques vanish. The two spheroids can be given arbitrary initial spins about their symmetry axes, but these spins will remain constant throughout the motion.
We probe the conservation of total energy and total orbital angular momentum as a function of time. The kinetic energy is readily computed from the velocities, and the gravitational potential energy from Eq. 12. The magnitude of the total orbital angular momentum J is where R and V are positions and velocities of A and B in the inertial frame (we denote all quantities relative to the inertial frame with capital symbols). We tested several cases of spheroid pairs with different shapes and initial parameters, but show only one case with two identical, oblate spheroids here. The parameters for the two spheroids are given in Table 1, and their tracjectories are shown in the top row of Fig. 2. The initial centroid-to-centroid distance (for the chosen initial velocity) for the two bodies is rather small, corresponding to 7 length units. Applying this, the bodies reach a minimum distance of r AB ≈ 2.55 length units in the simulation, i.e. almost touching since they both have radii of 1.0. All quantities are kept unit-less.
The orbits of A and B around the system centre of mass (at the origin) are confined between an inner radius, the pericentre distance, and an outer radius, the apocentre distance. The precession pattern is characterized by a radial and an azimuthal period, where the radial period, T r , is the time it takes for one of the bodies to travel from apocentre to pericentre and back to apocentre, and the azimuthal period T ψ is the period it takes for one of the bodies to travel 2π radians (Binney & Tremaine 2008). The non-closed orbits shown in Fig.2 have T r > T ψ , hence the bodies, upon passing the apocentre the second time, have completed more than 2π radians around their common centre of mass. The non-closed orbits are approximate ellipses where the major axis precesses with an angle ψ p equal to the angle overshooting 2π radians for each radial period. For the case displayed in Fig. 2, the Table 2 The characteristics of the ∆E and ∆J distributions for the simulated cases, with ∆E and ∆J being the difference in total energy and angular momentum between successive time steps. The columns named 'slope' contain the slope of a straight line fit to ∆E and ∆J versus time −3 × 10 −11 3 × 10 −08 5 × 10 −10 −4 × 10 −10 2D 7 × 10 −09 6 × 10 −06 3 × 10 −08 −8 × 10 −07 1 × 10 −10 7 × 10 −08 1 × 10 −09 1 × 10 −09 2D −1 × 10 −18 2 × 10 −16 9 × 10 −19 −9 × 10 −17 precession angle is ψ p = 18 • , and the precession rate defined as Ω p = ψ p /T r is 23 per time unit. It is positive, indicating that the major axis of the ellipse rotates in the same direction as the bodies rotate around the centre of mass. For two identical prolate spheroids with (a, c) = (0.25, 1.0) and with the same masses and initial parameters as the two oblate ones, the precession pattern rotates in the opposite direction with ψ p = −11 • (T r < T ψ ) and Ω p = −14 per time unit.

Two thin, non-coaxial disks
We also include a model with two thin disks to test the equations in Sect. 3.2. The two disks are identical with the same masses as the spheorids, see the details in Table 1. As the force and mutual potential are given by analytic expressions, the two disks compute much faster than two spheroids, hence we chose a shorter time step (set to 0.01) for the RK integrator.
The resulting trajectories are shown in the bottom two panels of Fig. 2. As can be seen, the trajectories for the two disks are very similar to that for the two spheroids, except that the precession angle and precession rate are larger, with ψ p = 30 • and Ω p = 38 per time unit.

Conservation of energy and orbital angular momentum
Because no external forces or torques are applied to the system, the total energy and orbital angular momentum are conserved. In reality however, these two quantities are affected both by the finite numerical precision of the simulations (numerical errors), and by the finite timestepping of the integration scheme (integration errors). In order to analyze the accuracy and energy drift in our simulations, we probe the change in total energy and angular momentum between successive time steps, termed ∆E and ∆J. In the absence of errors, both these quantities should be zero. The results are displayed in Fig. 3, where also a figure of how the distance r AB between the two centroids changes throughout the simulation is included. As expected, the errors increase as r AB approaches a minimum (pericentre) where the two bodies undergo the most rapid changes in kinetic and potential energy. The maximum value of ∆E and ∆J for the simulation with the two spheroids is 1 − 2 × 10 −7 . In other simulations (not showed here) where the initial distance between the two bodies is three times larger, the maximum errors decrease by 2-3 orders of magnitude. For the two disks, where an analytical expression is used, the errors in ∆J are at the machine accuracy level. The distributions of ∆E and ∆J are displayed as histograms in Fig. 4, and the accompanying Table 2 lists characteristics of these distributions such as the mean, the standard deviation We have also fitted a straight line to ∆E and ∆J as a function of time to measure if there is some drift in energy or angular momentum over the simulation, and the slope is also listed in Table 2. We expect both numerical errors and integration errors in the simulation, however the integration errors are likely small as they accumulate over time and our simulations are relatively short (just a few orbits). The numerical errors are random in nature and should have a mean of zero, whereas integration errors accumulating over time should have a nonzero average (Eastman & Pande 2016). If integration errors are negligible (which we argue that they are given the small time steps and the short duration of the simulations), then ∆E and ∆J should be dominated by random errors and their means should be close to zero. The standard deviations in ∆E and ∆J will in that case reflect the magnitude of the random errors. Fig. 4 and the numbers in Table 2 show that ∆E and ∆J have mean values of < 10 −8 , hence integration errors are negligible. The standard deviations σ are two to four orders of magnitude larger than the mean, also suggesting that the simulations are dominated by Error analysis of the simulations with two spheroids (left) and two disks (right). The centroid-tocentroid distance is r AB , and ∆E and ∆J are the change in energy and angular momentum between successive time steps. Simulation time is t random rather than integration errors. There is also negligible drift over the duration of the simulations as the slopes are smaller than the standard deviations. Table 2 also lists the standard error, σ/ √ N, which is a measure of the uncertainty in the mean, and we find values that are typically an order of magnitude larger than the mean, indicating that the mean cannot be determined properly. This happens because the changes in energy and angular momentum are correlated at short time scales (smaller variations on shorter time scales superimposed on variations happening on longer time scales).

Limiting cases for spheroids: point masses and thin disks
As a final check, we investigate how two cases with 1) nearly spherical and 2) very flattened spheroids compare with two interacting point masses and two interacting thin disks. In other words, we are comparing two cases were we use the surface integrals in Sect. 3 with the two cases with known analytical solutions for the mutual gravitational force, i.e. that of point masses and thin disks. We are thus testing the surface integration at the limit when it approaches the analytical solutions.
First, we set c/a = 0.99 for the two spheroids which make them nearly spherical, and thereafter c/a = 0.01 which flattens them considerably so that they become comparable to two disks. The mass ratio is kept at m A /m B = 1 for all cases, and the systems are given identical initial conditions. A pairwise comparison is made in Fig. 5 where we plot the relative difference in the distance between the two bodies, i.e. ∆r AB /r AB . The dotted lines show the relative differences between point masses and flat spheroids, and between thin disks and the c/a = 0.99 spheroids, and are of course the largest as these are cases that are not really comparable. However, the relative difference decreases when we compare c/a = 0.99 spheroids with point masses, and c/a = 0.01 spheroids with thin disks, confirming that the surface integration model produces results that approach the analytical cases. Finally, we note that  , and how nearly flat spheroids compare with thin disks (blue). The dotted lines compare cases that are clearly expected to be poor comparisons, i.e. nearly spherical spheroids with disks, and flattened spheroids with point masses the relative difference peaks each time the bodies are close, and that the amplitude at periapsis increases over time. From this figure, we also note that using two point masses as an approximation to two nearly spherical spheroids seems to be produce a better approximation than using two thin disks as an approximation to two flattened spheroids.

Summary and conclusions
We explore the application of a surface integration method to compute the force and mutual gravitational potential between two extended, rigid bodies. We first briefly describe how the surface integral expressions are obtained by vector calculus following Conway's publication (Conway 2016), and thereafter how the surface integrals can be computed for cases of two interacting ellipsoids/spheroids. By assuming that the gravitational potential of one body can be analytically expressed as a spheroid (MacMillan 1930), and integrating over a spheroid assumed to be the second body, we solve the equations of motion to test the method in a few simple planar cases. The resulting trajectories are non-closed orbits with either positive or negative precession, typical for bodies moving in non-Keplerian potentials. The surface integration scheme outlined in this paper can be applied to a spheroid or an ellipsoid, and if the other body has a gravitational potential that can be expressed analytically (like for ellipsoids and spheroids), the mutual force, torque and gravitational potential can be computed exactly between the two bodies (to within the limits of numerical integration over the surface of one of the bodies), thus avoiding truncation errors associated with series expansions of potentials. In a forthcoming paper, we apply the method to two ellipsoids with a full three dimensional treatment including torques (Ho et al., in prep.).
If the two bodies are coplanar, non-coaxial, infinitely thin disks, the surface integrals can be solved analytically, hence we also test the analytical expressions given by Conway (2016). Moreover, we demonstrate that in two limiting cases of almost spherical and significantly flattened spheroids, the spheroid solutions approach that of the analytical solutions for point masses and thin disks.
We have used the analytic expression by MacMillan (1930) for the gravitational potential of a spheroid inside the surface integral, but is it also possible to replace this with other forms of arbitrary potentials for studying problems involving e.g. an ellipsoid/spheroid and an arbitrary body.
The method we present is promising for studying dynamics between two solid, nonspherical bodies, for instance two ellipsoids or two spheroids, or a combination of the two. Since the method is exact and avoids using series expansions for gravitational potentials, it is applicable, and will produce exact results, in situations where the two bodies are close to each other. It is therefore interesting to apply it for studying dynamics of asteroids binaries, in particular contact binaries and cases of rotational fission where the two bodies are initially in contact with each other.