The quadrupole model for rigid-body gravity simulations

Icarus 225 (2013) 623–635

Contents lists available at SciVerse ScienceDirect

Icarus journal homepage: www.elsevier.com/locate/icarus

The quadrupole model for rigid-body gravity simulations Anthony R. Dobrovolskis a,⇑, D.G. Korycansky b a b

UCO/Lick Observatory, University of California, Santa Cruz, 245-3 NASA Ames Research Center, Moffett Field, CA 94035-1000, United States CODEP, Dept. of Earth Sciences, University of California, Santa Cruz, CA 95064-1077, United States

a r t i c l e

i n f o

Article history: Received 15 March 2012 Revised 19 April 2013 Accepted 22 April 2013 Available online 3 May 2013 Keywords: Asteroids, Dynamics Celestial mechanics Comets, Dynamics

a b s t r a c t We introduce two new models for gravitational simulations of systems of non-spherical bodies, such as comets and asteroids. In both models, one body (the ‘‘primary’’) may be represented by any convenient means, to arbitrary accuracy. In our first model, all of the other bodies are represented by small gravitational ‘‘molecules’’ consisting of a few point masses, rigidly linked together. In our second model, all of the other bodies are treated as point quadrupoles, with gravitational potentials including spherical harmonic terms up to the third degree (rather than only the first degree, as for ideal spheres or point masses). This quadrupole formulation may be regarded as a generalization of MacCullagh’s approximation. Both models permit the efficient calculation of the interaction energy, the force, and the torque acting on a small body in an arbitrary external gravitational potential. We test both models for the cases of a triaxial ellipsoid, a rectangular parallelepiped, and ‘‘duplex’’ combinations of two spheres, all in a point-mass potential. These examples were chosen in order to compare the accuracy of our technique with known analytical results, but the ellipsoid and duplex are also useful models for comets and asteroids. We find that both approaches show significant promise for more efficient gravitational simulations of binary asteroids, for example. An appendix also describes the duplex model in detail. Ó 2013 Elsevier Inc. All rights reserved.

1. Introduction Simulations of the dynamics of N > 1 gravitating objects are one of the basic tools of astrophysics and planetary science, and have been employed for several decades as computers have increased in speed and power. In planetary science, some typical recent applications include orbital mechanics, planetary rings, planet formation, and planetesimal interactions, including asteroid collisions (Leinhardt et al., 2000; Leinhardt and Richardson, 2002; Roig et al., 2003; Michel et al., 2004; Richardson et al., 2005; Walsh and Richardson, 2006; Korycansky and Asphaug, 2006, 2009). The simplest and most common approach to modeling the gravitational interaction of two bodies is to treat them both as gravitational monopoles: either as point masses, or as ideal spheres (i.e., objects with perfect spherical symmetry). Then by Newton’s Law of Universal Gravitation, their mutual gravitational potential energy W is simply GMm/r, where G is Newton’s constant of gravitation, M and m are the masses of the two bodies, and r is the distance between them. For many applications, though, this two-monopole representation is inadequate. A more general approach is to allow one of the two bodies to be non-spherical; examples include a small satellite orbiting an oblate planet or an irregular asteroid, as Dactyl or⇑ Corresponding author. Fax: +1 6506046779. E-mail address: [email protected] (A.R. Dobrovolskis). 0019-1035/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.icarus.2013.04.029

bits Ida. Without loss of generality, say that m is spherical but M is not. Then W = mU, where U is the ‘‘gravitational potential function’’ of M, which depends on its internal mass distribution. (Throughout this paper, we employ the physical convention where U = W/m rather than the geodetic convention U = W/m, as we find the former more intuitive.) Generally U is a more complicated function of the positions (locations and orientations) of both bodies than the simple GM/ r of Newton’s Law, and more accurate methods must be used. Sometimes the potential can be expressed as an analytic formula in closed form, as for ellipsoids or polyhedra, for which U can be calculated analytically or semi-analytically (e.g. Werner, 1994; Werner and Scheeres, 1996; Korycansky, 2004). More often, though, U is expanded in terms of spherical harmonic functions, with terms of order r1, r2, r3, r4, etc. When M is approximated as a gravitational quadrupole, with terms of order r1 and r3 only, U is given by MacCullagh’s formula (e.g. Danby, 1962). However, no closed-form analytic formula for the mutual gravitational energy W seems to be known when both bodies are nonspherical, as in a binary asteroid, for example. In such cases W can be expressed as a double expansion in spherical harmonics (e.g., Borderies, 1978; Schutz, 1981; Paul, 1988; Werner and Scheeres, 2005; Ashenberg, 2007; Tricarico, 2008; Boué and Laskar, 2009), although such expansions tend to be awkward analytically. For some methods, such as Smooth-Particle Hydrodynamics (SPH;

624

A.R. Dobrovolskis, D.G. Korycansky / Icarus 225 (2013) 623–635

e.g., Benz and Asphaug, 1999), individual objects are themselves made of assemblages of spherical elements, so that a sum of their potentials approximates the required non-spherical field. Unfortunately, such high-accuracy models are computationally expensive to implement for numerical applications, where many thousands of bodies might interact for many thousands of timesteps. In many cases such accuracy is not feasible; but it may not actually be required in situations where many small, non-spherical objects (such as comets, asteroids, small satellites, ring particles, ejecta blocks, rubble, spacecraft, etc.) all interact with each other, as well as with a single, much larger non-spherical body (such as an oblate planet, its ring system, an asteroid, etc.). In this paper we introduce two compromise approaches suited to the above cases. In both methods, one object (the ‘‘primary’’) may be described by essentially any means. In our first method, all of the other bodies are represented by small gravitational ‘‘molecules’’ consisting of a few point masses, rigidly linked together. In our second model, all of the other bodies are treated as point quadrupoles, with gravitational potentials including spherical harmonic terms up to the third degree (rather than only the first degree, as for ideal spheres or point masses). Our quadrupole method may be regarded as a generalization of McCullagh’s formula for the interaction between a monopole and a quadupole (Danby, 1962). It also is the gravitational analog of the quadrupole expansion for the energy of a charge distribution in an external electric field (see Jackson, 1975). The novelty of our method is the application to celestial mechanics of such a technique, apparently well known in classical electrodynamics. The next section provides a description of the gravitational molecule model. Then Section 3 gives a simple physical derivation of the gravitational quadrupole model. In Section 4, we examine the accuracy and efficiency of our methods in a number of tests, comparing the quadrupole approximation with accurate analytic and numerical calculations of the gravitational forces and torques between two bodies. One test body is a point mass, while the other is an ideal sphere, a triaxial ellipsoid, a cuboid (rectangular parallelepiped), or a ‘‘duplex’’ consisting of two spheres locked rigidly together. Finally, the Appendix describes the duplex model itself in detail. 2. Molecule model To derive our model, first we describe a simple symmetrical ‘‘molecule’’ consisting of a small number of point masses (or ideal spheres), but which reproduces the total mass and moments of inertia of a given object, as well as its gravitational field up through the quadrupole moments. Then we characterize the interaction of this molecule with an arbitrary external gravitational field. 2.1. MacCullagh’s formula Consider a rigid body of arbitrary shape and density distribution, with mass m and principal moments of inertia A 6 B 6 C. Let its center of mass be the origin of Cartesian coordinates (x, y, z) aligned with its principal axes A, B, C, respectively. MacCullagh’s approximation for its gravitational potential / (energy per unit mass) at some location r = (x, y, z) is

/¼ 

Gm G Gm G  3 ½A þ B þ C  3I ¼   5 ½ðB þ C  2AÞx2 r 2r r 2r þ ðA þ C  2BÞy2 þ ðA þ B  2CÞz2 ; ð1Þ

where

I ¼ ðAx2 þ By2 þ Cz2 Þ=r2

ð2Þ

is the moment of inertia of body m about the axis r (e.g., Ramsey, 1940, Section 4.6; Danby, 1962; Murray and Dermott, 1999).

Note that / above does not depend on the exact shape of m, but only on its mass and moment-of-inertia tensor, along with r. Thus the gravitational field of m can be approximated to quadrupole order O(1/r3) by the field of any other object with the same mass and inertia tensor. For example, the principal moments A, B, C of any solid body are equal to those of a homogeneous ‘‘equivalent ellipsoid’’ of the same mass m and principal radii a P b P c, where

a2 ¼

5 ½B þ C  A; 2m

2

b ¼

5 ½A þ C  B; 2m

c2 ¼

5 ½A þ B  C; 2m ð3Þ

and 2

A ¼ m½b þ c2 =5;

B ¼ m½a2 þ c2 =5;

2

C ¼ m½a2 þ b =5

ð4Þ

(Dobrovolskis, 1996). Then MacCullagh’s formula (1) may be rewritten as

/¼

Gm Gm 2 2 2  ½ð2a2 b c2 Þx2 þð2b a2 c2 Þy2 þð2c2 a2 b Þz2  r 10r 5 ð5Þ

In both Eqs. (1) and (5) above, note that the coefficient of x2 is always non-negative, while that of z2 is always non-positive (the coefficient of y2 may have any sign). Note also that both formulae (1) and (5) reduce to the monopole potential / = Gm/r for an equidimensional object where A = B = C and a = b = c. For example, a planet in the shape of an ideal cube would have no quadrupole moment, like an ideal sphere; but such a bizarre world would be quite artificial. For a non-equidimensional body, MacCullagh’s potential (1) or (5) is a good approximation in the ‘‘far field’’, distant from the body, but it becomes poorer in the ‘‘near field’’ close to the body, and actually becomes positive at some points near the origin 2

(0, 0, 0); for example, where z2 < 2CAB ¼ ða2 þ b  2c2 Þ=10 along 2m the axis of greatest inertia x = 0 = y. Of course, a net positive gravitational potential is quite non-physical, but this occurs only at dispffiffiffi tances r < a= 5  0:45 a, less than the longest dimension of the body. 2.2. Gravitational ‘‘molecules’’ There have been several previous attempts to represent nonspherical bodies as gravitational molecules consisting of a few point masses, rigidly linked together. Then the net energies, forces, and torques on the bodies are simply given by the sums of those on the individual masses, rather than obtained from MacCullagh’s formula. In his monograph on attitude stability of satellites, Beletskii (1965, pp. 92 and 98) introduced an example which he called a ‘‘three-dimensional dumbbell’’, consisting of six point particles of masses m1, m1, m2, m2, m3, and m3, symmetrically located at (0, 0, +l1), (0, 0, l1), (0, +l2, 0), (0, l2, 0), (+l3, 0, 0), and (l3, 0, 0), respectively. In related studies, Wang et al. (1991, 1992) adopted an asymmetric ‘‘molecule’’ of six point masses m1 through m6 at (x1, 0, 0), (x2, 0, 0), (0, x3, 0), (0, x4, 0), (0, 0, x5), and (0, 0, x6). Chauvineau et al. (1993) also used a similar model for a triaxial ellipsoid of mass m and principal radii a, b, c. Their ‘‘sixpoints’’ mass distribution consisted of six point particles, each of mass m/6, located at (±a/2, 0, 0), (0, ±b/2, 0), and (0, 0, ±c/2) in the body-fixed frame. Their model reproduces the ellipsoid’s total mass m, but the resulting moments of inertia (A00 = m[b2 + c2]/12, B00 = m[a2 + c2]/12, C00 = m[a2 + b2]/12) and the differences among them underestimate those of the ellipsoid (see above)

A.R. Dobrovolskis, D.G. Korycansky / Icarus 225 (2013) 623–635

by a factor of 12/5 = 2.4, which casts some doubt on their results. A simple model which correctly reproduces the quadrupole moments of an arbitrary body consists of a molecule built from seven point masses. Let such a molecule consist of a central point mass m7 at r = 0, rigidly attached to six other point masses m6 at (±a, 0, 0), (0, ±b, 0), and (0, 0, ±c) in the body-fixed frame ( is a dimensionless positive scale factor, initially unity, which will enable us to shrink the molecule later). This molecule has mass 6m6 + m7 and moments of inertia 2m62[b2 + c2], 2m62[a2 + c2], 2m62[a2 + b2]. These reproduce the mass m and moments A = m[b2 + c2]/5, B = m[a2 + c2]/5, C = m[a2 + b2]/5 of the original body and of its equivalent ellipsoid, provided that m6 = 2m/10 and m7 = m  6m6. From Eq. (1), it is clear that adding or subtracting any constant from A, B, and C leaves / unchanged. Consequently, the gravitational model need not even reproduce the original body’s moments of inertia, but only their differences C  B and B  A. This extra degree of freedom makes it possible to collapse the two point masses at (0, 0, ±c) onto the center, and to simplify the above ‘‘heptatomic’’ molecule into a planar ‘‘pentatomic’’ molecule consisting of only five point masses: a central point particle of mass m5 = m  4m4 at r = 0, attached rigidly to four outrigger point particles, each of mass m4 = 2m/10 = m6, located at (±a0 , 0, 0) and 0 0 (0, ±b0 , 0), where a 2 = a2  c2 = 5(C  A)/m and b 2 = b2  c2 = 5(C  B)/m. These reproduce the mass and moment differences of the original body, but the new moments of inertia are reduced to 0 0 0 0 A0 = mb 2/5 = C  B, B0 = ma 2/5 = C  A, and C0 = m[a 2 + b 2]/ 0 0 5 = A + B = 2C  A  B. To illustrate, Fig. 1 shows a cross-section through the a  b plane of a homogeneous pffiffiffi triaxial ellipsoid of mass m with principal radii a ¼ 2; b ¼ 2  1:414, c = 1, in the ratio typical of many asteroids and satellites (e.g., see Dobrovolskis, 1995). Its

625

principal moments of inertia are A = 3mc2/5, B = mc2, and C = 6mc2/5, while its reduced moments of inertia are A0 = mc2/ 0 5, B0 = 3mc2/5,pand = 4mc2/5, and its reduced equivalent radii ffiffiffi C 0 0 become a ¼ c 3; b ¼ c, and c0 = 0. Fig. 1 shows the equivalent pentatomic molecule (with  = 1) embedded inside the ellipsoid, like seeds in a fruit sliced through its core. The dotted circles represent equipotentials of an external point mass located at the plus sign. For  = 1, the six ‘‘outrigger’’ masses m6 of the heptatomic molecule lie on the surface of its equivalent ellipsoid. However, we may change  to any positive value without affecting the molecule’s total mass or inertia tensor, as long as we maintain m6 = 2m/10 and m7 = m  6m6 (m4 = 2m/10 and m5 = m  4m4 for the pentatomic molecule). If we shrink  below unity, the outrigger masses m6 (or m4) sink into the interior of the equivalent ellipsoid, while growing in mass to keep the same moments of inertia. At the same time, the central mass m7 (or m5) must decrease to maintain the p same ffiffiffiffiffiffiffiffi total mass m. This implies that m7 vanishes when  ¼ 3=5  0:77460, so that m6 becomes m/6, and the heptatomic molecule reduces to a hexatomic molecule. Simpffiffiffiffiffiffiffiffi ilarly, m5 vanishes when  ¼ 2=5  0:63246, so that m4 becomes m/4, and the pentatomic molecule reduces to a tetratomic molecule. Note that both the hexatomic and tetratomic molecule models as defined above are unique for a given object. It is possible to reduce the number of point masses in the gravitational molecule even further, to four or even three, in a triangular pattern; but the Cartesian symmetries of our rectangular pentatomic or heptatomic arrangements are convenient for taking the limit of a point quadrupole.

3. Quadrupole model

3

2

1 b’ a’

0

a’

b’

-1

-2 -2

-1

0

1

2

3

Fig. 1. A ‘‘pentatomic’’ molecule of five point masses. The solid ellipse is the x–y section pffiffiffi of a homogeneous triaxial ellipsoid of mass m, and principal radii a = 2, b ¼ 2  1:414, c = 1 aligned with x, y, z-axes originating at its center. The masses m4 = m/10 and m5 = 6m/10 in the molecule reproduce the same total mass and 0 moment pffiffiffi of inertia differences as the ellipsoid, provided b = c and a0 ¼ c 3  1:733c. From innermost to outermost, the dotted circles represent equipotentials of values 4, 3, 2, 1, 1/2, 1/3, 1/4, 1/5, 1/6, and 1/7 times GM, surrounding a point particle of mass M located at the plus sign.

For small molecules far from the primary, the potential is nearly the same at each elementary point mass. This can cause numerical inaccuracies from subtracting nearly equal quantities, particularly in evaluating torques. This is the same problem as that of numerical differentiation. In order to avoid this issue, we pass to the limit of infinitesimal molecules, as described below. In order to reduce the scale factor  even further, however, m7 (or m5) must go negative, while its gravity must become repulsive to compensate for the attraction of the six positive masses! This need not trouble us, however, because the central mass does not contribute to the moments of inertia, while the total mass of the molecule remains fixed. As  approaches zero from above, our gravitational molecule shrinks down to a point object, but with the same mass as the finite object, as well as its same gravitational quadrupole moments. To verify this, it is straightforward to expand the individual potentials of the five or seven point masses in power series of . Then their net sum recovers MacCullagh’s formula (1), to order 2. This is equivalent to finding the gravitational interaction of a quadrupole with a gravitational monopole (sphere or point mass). However, our formulation is much more general, and enables us to express the interaction of a quadrupole with an arbitrary gravitational field in terms of the external potential and its derivatives at a single point. Consider the gravitational interaction of a molecule of mass m, as described above, with some extended object of mass M, which we may call the ‘‘primary’’. The total potential energy W of their interaction is just the integral of m’s gravitational potential / over the mass distribution of M. However, we get the same answer whether we integrate / over M, or integrate the gravitational

626

A.R. Dobrovolskis, D.G. Korycansky / Icarus 225 (2013) 623–635

potential U of M over m. For our gravitational molecule, the latter integral is just the sum of the interaction energies of M with a small number of point masses. For the heptatomic molecule described above, this sum is

W7 ¼ m7 Uð0; 0; 0Þ þ m6 Uðþa; 0; 0Þ þ m6 Uða; 0; 0Þ þ m6 Uð0; þb; 0Þ þ m6 Uð0; b; 0Þ þ m6 Uð0; 0; þcÞ þ m6 Uð0; 0; cÞ;

ð6Þ

acceleration @ U/@x, @ U/@y, and @ U/@z. Explicitly (in vector shorthand), " # 2 2 m 2 @2 2 @ 2 @ F 7 ¼  mr0 U  a þ b þ c rU @x2 @y2 @z2 10 0 " # 1 @2 @2 @2 ¼  mr0 U  ðB þ C  AÞ 2 þ ðA þ C  BÞ 2 þðA þ B  CÞ 2 rU 4 @x @y @z 0 # " 2 2 1 @ @ @2 ¼  mr0 U  ðB þ C  2AÞ 2 þ ðA þ C  2BÞ 2 þðA þ B  2CÞ 2 rU 6 @x @y @z 0

ð11Þ

2

where again m7 = m  6m6 and m6 =  version, this reduces to

m/10. For the pentatomic for the heptatomic molecule, and

W5 ¼ m5 Uð0; 0; 0Þ þ m4 Uðþa0 ; 0; 0Þ þ m4 Uða0 ; 0; 0Þ 0

0

þ m4 Uð0; þb ; 0Þ þ m4 Uð0; b ; 0Þ;

ð7Þ

where m5, m4, a0 , and b0 are as defined above. The latter interaction is illustrated by the dotted equipotentials in Fig. 1. Provided that U(x, y, z) is an analytic function, we can expand W7 and W5 as MacLaurin series in . To lowest order, we thus obtain

" # 2 2 m 2 @2 2 @ 2 @ W7 ¼ mU0 þ a þb þc U 10 @x2 @y2 @z2 0 " # 1 @2 @2 @2 ¼ mU0 þ ðB þ C  AÞ 2 þ ðA þ C  BÞ 2 þ ðA þ B  CÞ 2 U 4 @x @y @z 0

ð8Þ 2

(plus terms of order  ) for the heptatomic molecule. By Laplace’s equation r2U = 0, we may subtract (a2 + b2 + c2)r2U/ 30 = (A + B + C)r2U/12 from Eq. (8) above, to obtain

"

" # 2 m 02 @ 2 02 @ F 5 ¼ mr0 U  a þb rU 10 @x2 @y2 0 " # 1 @2 @2 ¼ mr0 U  ðC  AÞ 2 þ ðC  BÞ 2 rU 2 @x @y

for the pentatomic version. Like Eqs. (8)–(10), Eqs. (11) and (12) above are equivalent by Laplace’s equation r2U = 0. To order (s/r)3, formulae (8)–(12) approximate the interaction energy W and force F which the primary M exerts on any object m whose size s is much smaller than its distance r from the center of the primary. Of course, the force f which m exerts on M is just F by Newton’s third law. The simplest expression for the torque s which m exerts about M’s center of mass is rM  f =rM  F, where rM is the location of M with respect to m. By a similar consideration of the individual forces on the gravitational molecule, taking due care over signs, we find that the torque T which M exerts about the center of m has Cartesian components

   @ 2 U  @ 2 U  @ 2 U  T x ¼ ðC  BÞ  ; T y ¼ þðC  AÞ  ; T z ¼ ðB  AÞ  ; @y@z @x@z @x@y 0

2

W7 ¼ mU0 þ

0

ð9Þ Formula (9) above is arguably more intuitive than Eq. (8) for nearly equidimensional objects (A  B  C); Eq. (9) also is formally equivalent to Eq. (4.24) of Jackson (1975), although the latter is expressed in terms of electrostatic fields. For the pentatomic molecule, Eq. (8) reduces to

" # 2 m 02 @ 2 02 @ W5 ¼ mU0 þ a þb U 10 @x2 @y2 0 " # 1 @2 @2 ¼ mU0 þ ðC  AÞ 2 þ ðC  BÞ 2 U þ Oð2 Þ: 2 @x @y

0

0

ð13Þ

2

m @ @ 2 2 ð2a2 b c2 Þ 2 þð2b a2 c2 Þ 2 30 @x @y # 2 2 @ þð2c2 a2 b Þ 2 U @z 0 " # 1 @2 @2 @2 ¼ mU0 þ ðBþC 2AÞ 2 þðAþC 2BÞ 2 þðAþB2CÞ 2 U: 6 @x @y @z

ð12Þ

0

along the x, y, and z axes, respectively. All of the torques on m are proportional to the differences among its moments of inertia, or its quadrupole moments. Thus the quadrupole model is the simplest one capable of gravitationally affecting the rotational energy or angular momentum of the small bodies. Formulae (11) and (13) above are formally equivalent to the solutions given for Problem 4.5 of Jackson (1975), although the latter are expressed in terms of electrostatic fields. Note that the interaction energy W depends on the value of the gravitational potential U of M, and on two or all three of its homogeneous second derivatives, while the torque T which M exerts about the center of m depends on all three of the mixed second derivatives of U. The force F which M exerts on m depends on all three first derivatives of U, and on six (or nine) of its 10 distinct third derivatives (@ 3U/ @x@y@z does not appear). In each case, U and its derivatives all are evaluated at x = y = z = 0, the center of mass of the quadrupole m. Note also that the position of each quadrupole is described by its attitude (orientation) in addition to its location, unlike the position of a monopole such as a point mass or an ideal sphere.

ð10Þ

0

In each case, the subscript 7 refers to the heptatomic molecule, and the subscript 5 to the pentatomic version, while the subscript zeroes indicate that U and its derivatives all should be evaluated at (0, 0, 0), the center of the quadrupole. In the limit as the scale factor  approaches zero, each molecule shrinks into a point quadrupole. Likewise, we obtain the force F which M exerts on m to lowest order in  by applying the same expansion to the components of

4. Analytical examples To illustrate, the most trivial application of the quadrupole method is to a constant potential field U = U0, such as that found inside a hollow spherical shell of constant thickness and density. U is trivially harmonic (i.e., r2U = 0, as for any valid gravitational potential in free space), since all of its derivatives vanish. Then from Eqs. (10)–(12), the force F and torque T on the quadrupole both vanish too, while the interaction energy W is just mU0 from

627

A.R. Dobrovolskis, D.G. Korycansky / Icarus 225 (2013) 623–635

Eq. (8) or (9). Thus the quadrupole acts just like a gravitational monopole of the same mass m in this field; in fact, any extended body would behave the same. 4.1. Constant acceleration Another very simple example of our method is the interaction of a point quadrupole with the gravitational field of a ‘‘flat Earth’’, or on either side of an infinite planar sheet of uniform mass density r per unit area. Such a field might approximate that of a nearby disk, such as a broad flat planetary ring, for example. The gravitational potential at a distance Z from such a sheet is just U = gZ, where the constant g = Gr is the magnitude of the gravitational acceleration g. In the rotated (x, y, z)-frame originating in the quadrupole, Z = d + ax + by + cz, where d is the distance of the quadrupole from the sheet Z = 0, and a, b, and c are the direction cosines of the Z-axis on the x-, y-, and z-axes, respectively. Then U is of first degree in x, y, and z, and the corresponding gravitational acceleration g = rU = grZ = gai  gbj  gck is a constant vector. Because all of the second (and higher) derivatives of this potential vanish, it is harmonic outside the mass sheet, while the net interaction energy of a point quadrupole with this field is simply W = mU0 = mgd from Eq. (8) or (9). The corresponding force on the quadrupole is just F = mg from Eq. (10) or (11), while the torque T on it vanishes by Eq. (12). Again, the quadrupole acts just like a gravitational monopole of the same mass m. In fact, any extended body of mass m outside the mass sheet would experience the same force mg, torque 0, and net interaction energy mgd in this field, where d is now the distance of its center of mass from the mass sheet. 4.2. Tidal potential Another simple application of the quadrupole method is to a tidal potential. Say that a perturber of mass M lies at a distance d along the Z-axis of an (X, Y, Z)-frame centered at the quadrupole. Then to lowest degree in 1/d, its tidal potential takes the form 3

2

2

3

2

2

2

U ¼ GMd ½X =2 þ Y =2  Z  ¼ GMd ½R  3Z =2; 2

2

2

Homogenous derivatives

First

@u/@x = x  3a[ax + by + cz] @u/@y = y  3b[ax + by + cz] @u/@z = z  3c[ax + by + cz]

Second

@ 2u/@x2 = 1  3a2 @ 2u/@y2 = 1  3b2 @ 2u/@z2 = 1  3c2

3

¼ GMd ½x2 ð1  3a2 Þ þ y2 ð1  3b2 Þ þ z2 ð1  3c2 Þ  6bcyz  6acxz  6abxy=2:

@ 2u/@y@z = 3bc @ 2u/@x@z = 3ac @ 2u/@x@y = 3ab

3

T x ¼ þ3GMðC  BÞbc=d ;

3

T y ¼ 3GMðC  AÞac=d ;

3

T z ¼ þ3GMðB  AÞab=d :

ð17Þ

Note that this situation is opposite from the preceding case of constant acceleration, where the torque vanishes but the force does not. 4.3. Monopole potential An important application of the quadrupole method is to the gravitational potential of a point mass or ideal sphere. In this case qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi U = GM/d, where d  ðx  xM Þ2 þ ðy  yM Þ2 þ ðz  zM Þ2 is the distance between the test point and the center of the sphere at rM = (xM, yM, zM). Because this potential is isotropic, it is not necessary to transform coordinates from some (X, Y, Z) inertial reference frame; but rather, the derivatives can be taken in the (x, y, z)-frame fixed in the body of the quadrupole. Table 2 displays the first, second, and third derivatives of u  1/d in this frame. Substituting the prescribed derivatives into either formula (8) or (9) then gives

W¼ 

GMm GM  5 ½ðB þ C  2AÞðx  xM Þ2 þ ðA þ C  2BÞðy  yM Þ2 d 2d

þ ðA þ B  2CÞðz  zM Þ2 

ð18Þ

(although the pentatomic form is somewhat simpler to manipulate). Evaluating Eq. (17) above at (x, y, z) = (0, 0, 0) then leaves

W¼

GMm GMK  5 ; rM 2r M

ð19Þ

where we define K  ðB þ C  2AÞx2M þ ðA þ C  2BÞy2M þ ðAþ B  2CÞz2M for convenience. From Eq. (10) or (11), the force F which M exerts on m is

GMmxM GMðB þ C  2AÞxM 5GMKxM  þ ; r 5M 2r7M r 3M GMmyM GMðA þ C  2BÞyM 5GMKyM  þ ; Fy ¼ r 5M 2r 7M r3M GMmzM GMðA þ B  2CÞzM 5GMKzM Fz ¼  þ : r 5M 2r 7M r 3M

Fx ¼ ð15Þ

U is still a quadratic function of x, y, and z, so its second derivatives are constants, while its third (and higher) derivatives all vanish. Table 1 displays the first and second derivatives of U, normalized by GMd3. Note that U and its three first derivatives all vanish at the origin (0, 0, 0). Then F also vanishes from Eq. (10) or (11), so this tidal potential exerts no net force on the quadrupole. From Eqs. (8) and (9), the net interaction energy is 2

½a2 ð1  3a2 Þ þ b ð1  3b2 Þ þ c2 ð1  3c2 Þ 3 10d GM ¼ 3 ½ðC  AÞð1  3a2 Þ þ ðC  BÞð1  3b2 Þ; 2d

Mixed derivatives

while from Eq. (12), the torques which M exerts about the center of m are

ð14Þ

U ¼ GMd3 ½x2 þ y2 þ z2  3ðax þ by þ czÞ2 =2

GMm

Order

2

where we define R  X + Y + Z . Note that this potential is harmonic, since r2U = GMd3[1 + 1  2] = 0. In the more usual spherical coordinates, U may be written as GMd3R2[1  3 cos2 f]/2, where f = Arccos(Z/R) is the angular distance from the direction of the perturber. In the rotated (x, y, z)-frame aligned with the quadrupole, R2 2 = r  x2 + y2 + z2, while Z = ax + by + cz. Then the tidal potential becomes

W¼

Table 1 Derivatives of normalized tidal potential u  [x2 + y2 + z2  3(ax + by + cz)2]/2. Note that @ 2u/@x2 + @ 2u/@y2 + @ 2u/@z2 = 0, as required by Laplace’s equation; recall that a2 + b2 + c2 = 1.

ð16Þ

ð20Þ

Results (19) above are equivalent to formula (2) of Ramsey (1940, Section 4.61). Note that this force is non-central; that is, not directed along the line between the centers of the monopole and quadrupole. From Eq. (12), the corresponding torque T which M exerts on m is

T x ¼ þ3GMðC  BÞyM zM =r 5M ; ¼ 3GMðC 

AÞxM zM =r 5M ;

Ty T z ¼ þ3GMðB  AÞxM yM =r 5M :

ð21Þ

628

A.R. Dobrovolskis, D.G. Korycansky / Icarus 225 (2013) 623–635 Table 2 Derivatives of normalized monopole potential u  [(x  xM)2 + (y  yM)2 + (z  zM)2]1/2, up to third order. Note that @ 2u/@x2 + @ 2u/@y2 + @ 2u/@z2 = 0, as required by Laplace’s equation. Homogenous derivatives

Mixed derivatives

@u/@x = [xM  x]u3 @u/@y = [yM  y]u3 @u/@z = [zM  z]u3 @2 u @x2 @2 u @y2

¼ ½þ2ðx  xM Þ2  ðy  yM Þ2  ðz  zM Þ2 u5

@ 2u/@y@z = 3[y  yM][z  zM]u5

¼ ½ðx  xM Þ2 þ 2ðy  yM Þ2  ðz  zM Þ2 u5

@ 2u/@x@z = 3[x  xM][z  zM]u5

2

¼ ½ðx  xM Þ2  ðy  yM Þ2 þ 2ðz  zM Þ2 u5

@ 2u/@x@y = 3[x  xM][y  yM]u5

@ u @z2

@ 3u/@x3 = 9[x  xM]u5  15[x  xM]3u7

@ 3u/@y3 = 9[y  yM]u5  15[y  yM]3u7

@ 3u/@z3 = 9[z  zM]u5  15[z  zM]3u7

@3 u @x@y2

¼ 3½x  xM u5  15½x  xM ½y  yM 2 u7

3

@ u @x@z2

¼ 3½x  xM u5  15½x  xM ½z  zM 2 u7

@3 u @y@x2

¼ 3½y  yM u5  15½y  yM ½x  xM 2 u7

@3 u @y@z2

¼ 3½y  yM u5  15½y  yM ½z  zM 2 u7

@3 u @z@x2 @3 u @z@y2

¼ 3½z  zM u5  15½z  zM ½x  xM 2 u7 ¼ 3½z  zM u5  15½z  zM ½y  yM 2 u7

@ 3u/@x@y@z = 15[xM  x][yM  y][zM  z]u7

Thus the monopole potential produces both a net force and torque on the quadrupole. Note that expressions (20) above are equivalent to Eq. (16) for the torque produced by a tidal potential, and to formula (3) of Ramsey (1940, Section 4.61), as well as Eqs. (5.43)– (5.45) of Murray and Dermott (1999). It is easy to verify that Eq. (18) is equivalent to MacCullagh’s formula (1) when M is a monopole (ideal sphere or point mass), but our quadrupole method is much more general, and can be applied to any mass distribution M for which an analytic (or semianalytic) potential U is known. Besides spheres and infinite sheets or slabs, this includes straight rods, prolate spheroids, oblate spheroids, triaxial ellipsoids, polyhedra, quadrupoles and other potentials expressed in terms of spherical harmonics, etc. We can also handle situations where many small objects interact (each with its own quadrupole moments), such as satellite swarms, rings, or rubble piles. Our quadrupole method requires only that we take all of the necessary derivatives of U(X, Y, Z), and then transform them geometrically into the rotated x, y, z system of the quadrupole. In some cases, generating such derivatives analytically may be inconvenient or impossible. Symbolic math packages or numerical calculation may be a viable strategy to deal with those situations. 5. Numerical tests We have also performed preliminary numerical tests in order to compare our gravitational molecule and quadrupole models with more exact results from full potentials. We compute not only the interaction energies W, but also forces F and torques T for numerical integrations of the objects’ dynamics. We present here the results of several calculations involving different types of test objects (spheres, ellipsoids, cuboids, and duplexes) in the potential of a unit monopole. By analogy with ellipsoids and spheroids, we use the term ‘‘cuboids’’ for rectangular parallelipipeds (bricks). We also introduce a useful toy model we call a ‘‘duplex’’, consisting of two non-concentric spheres rigidly connected together (not to be confused with a dipole; see Appendix). Except for spheres, duplexes may be the physically realizable objects with the simplest gravitational fields, with convenient analytic expressions for their potentials (being just the sum of two monopole potentials), as well as for the various derivatives thereof (see Table 2). Therefore we use them below as versatile test bodies for comparison with the quadrupole and molecule models; but the

Appendix develops the duplex model in detail as a valuable representation of certain primaries, such as contact binary asteroids. 5.1. Interaction energy For consistency with the previous discussion, we denote the ‘‘test object’’ by (lower-case) m, with principal moments of inertia A, B, and C. Each of our test objects was oriented so that its principal axes of inertia lay along the coordinate axes, with x along its A (longest) axis and z along its C (shortest) axis. Then we computed the interaction energy W of each test object m with a unit point mass M at (x, y, z) in four different ways: as Wa, by applying the exact analytic potential / of m to M; as Wq, by applying the quadrupole approximation (8); and as W4 and W6, by applying the tetratomic and hexatomic molecule models, respectively. For the triaxial ellipsoid and triaxial cuboid cases, we evaluated W all four ways for M located at nk = 50 latitudes p/2 < ki < p/2 by nh = 100 distinct longitudes p < hj < p at 10 different distances 2 6 r 6 500. In contrast, the symmetry of the sphere and duplexes about the x-axis means that all the information about their interactions with M can be found by placing M in the equatorial xy plane, so we evaluated W only at zero latitude (nk = 1). Our figure of merit DWq for the relative accuracy of Wq is given by the RMS difference between Wq and Wa, taken at the same distance r, but summed over latitude and longitude, and normalized by the RMS average of Wq and Wa:

DWq 

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi .Xn ;n Xnk ;nh k h 4 ij ½Wq ði; jÞ  Wa ði; jÞ2 ½Wq ði; jÞ þ Wa ði; jÞ2 : ij ð22Þ

Figures of merit DW4 and DW6 for the molecule models were determined analogously. Fig. 2 plots DW6, DW4, and DWq for six different test objects. The middle panel on the left-hand side refers to a homogenous triaxial ellipsoid of unit mass and principal semi-axes a = 2.0, b = 1.4, and c = 1.0 units, with corresponding moments of inertia A = 0.592, B = 1.000, and C = 1.192 from Eq. (4). Similarly, the bottom left panel refers to a homogenous triaxial cuboid of unit mass and semiedges a = 2.0, b = 1.4, and c = 1.0 units. Its corresponding moments of inertia are A = m[b2 + c2]/3  0.987, B = m[a2 + c2]/3  1.667, and C = m[a2 + b2]/3  1.987 (comparepffiffiffiffiffiffiffiffi Eq. (4)), whilepits ffiffiffiffiffiffiffiffiequivalent ellipsoid has principal radii of 5=3a  2:582, 5=3b  1:807,

A.R. Dobrovolskis, D.G. Korycansky / Icarus 225 (2013) 623–635

1

10

100

10

100

0

1000 0

-2

-2

-4

-4

-6

-6

-8

-8

symmetric duplex

sphere -10

-10

-2

-2

-4

-4

-6

-6

-8

-8

asymmetric duplex

ellipsoid -10

-10

-2

-2

-4

-4

-6

-6

-8

-8

unbalanced duplex

cuboid -10

-10

-12 1

10

100

10

100

-12 1000

Distance r Fig. 2. The relative accuracy DW of the gravitational energy of a unit monopole at a distance r from bodies described in the text, as labeled. Hexatomic molecules are denoted by , while tetratomic molecules are denoted by +. Squares denote point quadrupoles, while for comparison, the dotted lines plot power laws, as labeled.

pffiffiffiffiffiffiffiffi and 5=3c  1:291 from Eq. (3). Both of these bodies are centered at the origin, and for each one, DW starts out at a few percent for r = 2 (extrapolating down to order unity at r = 1), but drops off approximately as the inverse fourth power of the distance r, until it approaches a limiting accuracy of 1010. This lower limit presumably is due to the accumulation of roundoff error. The right-hand side of Fig. 2 refers to three different duplexes. The uppermost right panel refers to a symmetric duplex consisting of two spheres of equal unit masses and radii, centered at (1, 0, 0) and (+1, 0, 0). For this simple body, DW behaves almost the same way as it does for the ellipsoid and cuboid. In contrast, the middle right panel refers to an asymmetric duplex consisting of one sphere of mass 1 and radius 1, centered at (1, 0, 0), and a second, denser sphere of mass 2 but radius 1/2, centered at (+1/2, 0, 0). Again, the accuracy starts out at a few percent for r = 2, but now DW declines only as the inverse cube of the distance. This difference is explained as follows. The net potential field / at a distance r from the center of mass of an arbitrary object of mass m is dominated by the monopole term Gm/r. However, / also contains terms varying as r3 proportional to the quadrupole moments of m, other terms varying as r4 proportional to its hexapole moments, others varying as r5 proportional to its octopole moments, and so on. Our quadrupole approximation includes the quadrupole terms, but neglects hexapole and higher terms. Therefore its absolute accuracy should be of order r4 in the potential, while the corresponding relative accuracy should be of order r4/r1 = r3. This is just what we find for the asymmetric duplex. In contrast, the ellipsoid and cuboid have no hexapole moment, because they each possess three mutually perpendicular planes of

629

symmetry. Likewise, the symmetric duplex has no hexapole moment either, because it is symmetrical about the yz-plane as well as the x-axis. Consequently, the accuracy is dominated by these bodies’ octopole moments. Therefore the absolute accuracy of the potential should be of order r5, while its relative accuracy should be of order r5/r1 = r4. This is just what we find for the symmetric duplex, as well as for the ellipsoid and cuboid. In this respect, these shapes are unrealistically simple, while the asymmetric duplex is more realistic, despite its radial symmetry about its long axis. The bottom right panel of Fig. 2 presents the even more extreme example of an ‘‘unbalanced’’ duplex consisting of one sphere of mass 1 and radius 1, centered at (1, 0, 0), and another sphere of radius 1 but mass 2, centered at (+1, 0, 0). In this case, the center of mass does not lie at the origin, but at (1/3, 0, 0). As a consequence, the duplex possesses a dipole moment (of unit strength), whose potential falls off as r2, so its relative contribution varies as r2/r1 = r1, as the figure shows. Thus small errors in locating a body’s center of mass ultimately can dominate the errors in its gravitational field far from the object. Surprisingly, the accuracy of both molecule models is comparable to that of the quadrupole model in all of the above cases. This implies that the issue of numerical differentiation is not very important after all. As a control, the top panel on the left-hand side refers to an ideal sphere of unit mass and unit radius, centered at the origin. For the hexatomic model, DW6 again starts out at 0.01 for r = 2, but then drops off approximately as the inverse fourth power of the distance r. However, DW4 for the tetratomic molecule and DWq for the point quadrupole are vanishingly small, because for an ideal sphere, both of the latter models reduce to a single point mass. 5.2. Forces and torques As a result of their gravitational interaction, each test object m experiences forces and torques from the unit monopole M located at r = (x, y, z) with respect to m. The force exerted on m by M is F, and is equal but opposite to the force f from m on M. We denote the ‘‘orbital torque’’ (the moment of f about m) by r  f, and the ‘‘spin torque’’ about the center of mass of m by T. The sum of both torques is zero due to conservation of angular momentum, and we take of advantage of that fact to evaluate our results and to check their accuracy. We calculated both F and T at the same points where we evaluated the interaction energy W. But unlike W, we calculated F and T by five different methods rather than only four: in addition to the exact analytic formulae, the quadrupole approximation, and the two molecule models, we also evaluated F and T by means of surface quadratures, as done by Korycansky (2004) for interacting polyhedra. Our figure of merit D for the relative accuracy of the components of F or T is now given by formulae analogous to Eq. (21) for the RMS difference between the exact analytic evaluation and the quadrupole approximation, either molecule model, or the numerical quadrature. We describe the results below, but do not display them to save space. For the ellipsoid from the preceding subsection, we used n = 4, 8, 16, or 32 quadrature points in both latitude and longitude, so that there were a total of n2 = 16, 64, 256, or 1024 grid points over the surface of the ellipsoid. The results show that the accuracy D of the quadrupole approximation and molecule models varies as r4 for the forces (just as for DW in Fig. 2), but only as r2 for the torques, until it begins to degrade between r = 200 and r = 500. At such large radii, D is dominated by accumulating numerical errors in our routines to evaluate the elliptic integrals required for an analytic evaluation of the ellipsoidal potential. The quadrupole approximation and molecule models are more accurate than Gaussian quadrature for the 4  4 grid, but roughly comparable

630

A.R. Dobrovolskis, D.G. Korycansky / Icarus 225 (2013) 623–635

to the 8  8 grid. The 16  16 and 32  32 grids are more accurate in the near field, but their accuracy tends to level out at 108, while the errors in the torques actually rise for r J 50. For the cuboid from the preceding subsection, quadratures were done separately on each rectangular face, with n = 4, 8, 16, or 32, so that each face had n2 = 16, 64, 256, or 1024 points, for a total of 6n2 = 96, 384, 1536, or 6144 points. The results show that the accuracy D of the quadrupole approximation and molecule models again varies as r4 for the forces, but only as r2 for the torques, until accumulating numerical errors begin to degrade D between r = 200 and r = 500. For the numerical quadratures, D generally levels out at 109 (or 107 for the 4  4 grid), while the torques again degrade for r J 50. All of these quadratures are now more accurate than the quadrupole approximation or molecule models, except for the 4  4 grid in the far field forces (r J 100). Comparison with the previous results also reveals that the numerical quadratures are significantly more accurate for the cuboid than for the ellipsoid. We also performed Gaussian quadratures over the duplexes of the preceding subsection, as well as the sphere of unit mass and radius. The quadratures over the sphere were done using n = 4, 8, 16, or 32 points in both latitude and longitude, just as for the ellipsoid. The same grid was used for all of the duplexes as well, but for each of the two component spheres, so that there were twice as many grid points in all, for a total of 2n2 = 32, 128, 512, or 2048 integration points on each duplex. The results for the sphere indicate that the errors in the numerical quadratures for the forces mostly lie between those for the ellipsoid and for the cuboid, but the error in the quadrupole approximation and molecule models is again vanishingly small, just as for DW in Fig. 2. The torque on the sphere and its error vanish as well, for all of the models, as expected. For the symmetric duplex, the errors in the forces again generally lie between those for the ellipsoid and for the cuboid. The errors in the torque mostly lie between those two cases in the near field (r [ 50); but in the far field (r J 50), the errors in the torque for the duplex level off and are smallest, because duplexes are not as susceptible to roundoff error. For the quadrupole approximation and molecule models, D varies as r4 for the forces, and as r2 for the torque, just as for both the ellipsoid and cuboid. The quadrupole approximation and molecule models are more accurate than Gaussian quadrature for the 4  4 grid, but roughly comparable to the 8  8 grid, as for the ellipsoid. The 16  16 and 32  32 grids are more accurate in the near field, but their accuracy tends to level out at 109, as for the cuboid. For the asymmetric duplex, in contrast, the accuracy of the quadrupole approximation and molecule models varies as r3 for the forces (comparable to the 8  8 quadrature), and only as r1 for the torque. This is explained as follows. For an arbitrary mass distribution centered at the origin, the net force on a distant point mass should vary as r2, the inverse square of the distance, due to the body’s monopole moment. In contrast, the net torque on a distant point mass should scale as r3, due to the body’s quadrupole

moments. Our quadrupole approximation and molecule models reproduce these moments, but neglect the body’s hexapole moments, introducing absolute errors of order r5 in the force, and of order r4 in the torque. The corresponding relative errors then should vary as r5/r2 = r3 in the force (just as for DW in Fig. 2), but as r4/r3 = r1 in the torque. This is just what we find for the asymmetric duplex. In contrast, the symmetric duplex has no hexapole moment, like the ellipsoid and cuboid. Consequently the accuracy is dominated by these bodies’ octopole moments. The resulting absolute errors vary as r6 in the force, and as r5 in the torque. The corresponding relative errors should scale as r6/r2 = r4 in the force (just as for DW in Fig. 2), but as r5/r3 = r2 in the torque. This is just what we find for the symmetric duplex, as well as for the ellipsoid and the cuboid. In the most extreme case of the unbalanced duplex, the residual dipole moment exerts a net force varying as r3, and a net torque varying as r2. The corresponding relative errors in the force vary only as r3/r2 = r1, just as for DW. In comparison, the relative errors in the torque should be of order r2/r3 = r+1, or directly proportional to the distance! In fact we find that the error in the torque rises somewhat with distance up to r  20, but then levels off at values slightly greater than unity. Again, we find that errors in locating a body’s center of mass can dominate its gravitational field far from the object. 5.3. Speed versus accuracy Of course, accuracy is not the only criterion for numerical simulations; speed also is a factor. In order to compare the speed of the quadrupole technique and molecule models to that of Gaussian quadratures, we benchmarked the simulations described above. Table 3 lists the CPU times required for each evaluation of the forces and torques between the test objects and a point mass at r = 5. We chose this distance because the near field is probably the most interesting range for most problems, especially for physical situations involving collisions, binary formation, or Yarkovskyrelated effects. For each type of object, it is clear that the quadrupole and molecule models are about as fast as evaluating the analytic expression, or even faster; and much faster than all of the numerical quadrature schemes. Note that the tabulated CPU times for the quadratures scale roughly as the total number of evaluation points in the numerical grid, for a given body. The duplex and ellipsoid quadratures require much more CPU time than the cuboid quadrature due to the necessity of evaluating more trigonometric expressions in calculating the force and torque. There is an additional penalty of a factor of two for the duplexes because two bodies are involved in the quadrature. Still, what really matters is the speed of computation for comparable accuracy. For the ellipsoid simulations, the accuracy of the quadrupole and molecule methods is comparable to that of the 8  8 quadrature, but the quadrupole and molecule methods are

Table 3 Mean CPU times (in microseconds on a 1.4 GHz workstation) for each evaluation of the forces and torques between a unit monopole and various test objects. Method

Sphere

Ellipsoid

Cuboid

Symmetric duplex

Asymmetric duplex

Unbalanced duplex

Analytic Quadrupole Tetratomic molecule Hexatomic molecule

0.17 0.64 0.13 0.80

6.62 0.56 0.52 0.79

9.07 0.60 0.58 0.54

0.25 0.63 0.41 0.80

0.25 0.63 1.53 0.79

0.25 0.65 0.40 0.80

4  4 quadrature 8  8 quadrature 16  16 quadrature 32  32 quadrature

72.8 284 1129 4494

258 1028 4092 16,340

10.8 42.3 169 652

145 569 2257 8969

145 569 2256 8969

145 569 2255 8976

631

A.R. Dobrovolskis, D.G. Korycansky / Icarus 225 (2013) 623–635

some three orders of magnitude faster. For the cuboid case, the quadrupole and molecule methods are comparable to the 6  4  4 quadrature, or less accurate; but the quadrupole and molecule methods are still an order of magnitude faster. For the duplexes, the accuracy of the quadrupole and molecule methods generally fall between those of the 4  4 and 8  8 quadratures; but the quadrupole and molecule methods are from two to three orders of magnitude faster.

6. Discussion In this paper we have outlined two formulations for approximating the response of a body to an external gravitational field, including the body’s quadrupole moments: the gravitational molecule and the point quadrupole. Both methods are intended to include one primary, which may be described by essentially any method, along with any number of simpler bodies, which are modeled either as small gravitational molecules, or as point quadrupoles. Both techniques show promise for improved studies of binary and multiple asteroids, for example. Note that, like the primary but unlike gravitational monopoles such as ideal spheres or point masses, the position of each molecule or quadrupole is described by its attitude (orientation) in addition to its location. The point quadrupole method also requires knowledge of the external potential up to its third derivatives, or of the force up to its second derivatives. In contrast, the gravitational molecule method requires knowledge only of the force itself, or of the potential up to its first derivatives. Both methods are simple to implement, and represent substantial improvements in the speed of computation over numerical quadrature techniques, for comparable accuracy. To our surprise, both versions of the molecule model tested here performed about as well as the point quadrupole model in terms of both speed and accuracy. For comparison, the expansion given by Schutz (1981) of the gravitational interaction energy W between two extended bodies (of masses M and m and size scales S and s, respectively) includes the monopole–monopole term GMm/r, terms of order S2/r3 and s2/r3 relative to GMm corresponding to the quadrupole–monopole and monopole–quadrupole interactions, and quadrupole–quadrupole terms of order S2s2/r5. However, Schutz (1981) neglects comparable terms of order S4/r5 and s4/r5 corresponding to octopole–monopole and monopole–octopole interactions, respectively, as well as terms of lower order S3/r4 and s3/r4 corresponding to hexapole–monopole and monopole–hexapole interactions. In contrast, our expressions (6)–(9) for the gravitational interaction energy W between a molecule or quadrupole and the external potential (primary) include terms up to only second degree in the spherical harmonic expansion of the quadrupole’s potential, but they include the external potential to all degrees and orders. Furthermore, our techniques can include dipole moments1 of the external potential (primary); in contrast, Schutz (1981) also neglects hexapole–dipole and dipole–hexapole terms of order S3s/r5 and Ss3/r5, quadrupole–dipole and dipole–quadrupole terms of order S2s/r4 and Ss2/r4, dipole–dipole terms of order Ss/r3, and even dipole– monopole and monopole–dipole terms of order S/r2 and s/r2. For comparison, our quadrupole model includes all terms of orders Sn/ rn+1 and Sns2/rn+3, but neglects all terms of orders Sns3/rn+4 and higher in s, while all terms of orders Sns/rn+2 vanish. In applications of our methods, the interaction between two molecules or point quadrupoles should include only the quadrupole–monopole and monopole–quadrupole terms, for consistency;

including the quadrupole–quadrupole terms would be inconsistent with the neglect of octopole–monopole and monopole–octopole interactions, and especially with the neglect of hexapole–monopole and monopole–hexapole terms. Acknowledgments We thank two anonymous reviewers for their helpful suggestions. This work was supported by NASA Planetary and Geophysics Program Grant NNX07AQ04G. Appendix A. The duplex A.1. Construction We define a duplex as a pair of ideal spheres, connected together rigidly, as sketched in Fig. 3. They may be fully detached or disjoint, as in Fig. 3a; tangent to each other from outside, as in Fig. 3b; interpenetrating or conjoint, as in Fig. 3c; or even nested, as in Fig. 3d. Note that the two spheres need not have the same size, mass, or mean density. Neither one needs to be uniform in density either, as long as each is spherically symmetric. This construct is quite simple, but versatile, and useful both as an ideal test object and as a model for real bodies, such as certain comets and asteroids. A slightly simpler model has been used before, consisting of two point masses connected by a massless rigid rod. Numerous authors (Klemperer and Baker, 1957; Schindler, 1960; Rosner, 1960; Moran, 1961; Beletskii, 1965; Robinson, 1973; Chauvineau et al., 1993; Prieto-Llanos and Gómez-Tierno, 1994; Elipe et al., 2008) refer to this as a ‘‘dumb-bell’’ or a ‘‘dumbbell’’, while Hirabayashi et al. (2010) call it a ‘‘two-connected-mass’’ or a ‘‘rod-connected’’ body.

a

b

c

d

1

Dipole terms correspond to an offset of the origin from the body’s center of mass. They can be useful for asymmetrical objects, or for those whose center of mass is unknown, like most asteroids. Our ‘‘molecules’’ and point quadrupoles have no dipole moments, but the primary may.

Fig. 3. Various configurations of the duplex. (Panel a) Detached, (Panel b) tangent, (Panel c) interpenetrating, and (Panel d) nested. The dotted line indicates the axis of symmetry, while the two crosses on it denote the centers of both spheres.

632

A.R. Dobrovolskis, D.G. Korycansky / Icarus 225 (2013) 623–635

Prieto-Llanos and Gómez-Tierno (1994) also refer to their model as a ‘‘mass dipole’’, but this is a misnomer, as explained below. A true dipole consists of two point masses (or electrical charges, or magnetic monopoles, etc.) of opposite signs +m and m, but separated by a distance d. A pure dipole field is produced in the limit as d approaches zero, but +m and m grow without bound, such that their dipole moment md remains finite. This is analogous to our derivation of the quadrupole field in the main text. In contrast to all of the above models, each sphere in our duplex model may have a finite size, density, and moment of inertia. Consider a duplex consisting of two spheres labeled a and b, with individual masses Ma and Mb, radii Ra and Rb, and central moments of inertia Ia and Ib, respectively, whose centers are separated by a fixed distance d. In order to preserve the simplicity of this model, both spheres must rotate and translate together as a single rigid body, as if connected by a massless rigid rod, even in the case where the two spheres are disjoint, as in Fig. 3a. Marchis et al. (2005) refer to a similar disjoint model for Asteroid 121 Hermione as a ‘‘detached peanut’’ configuration. Such a model might also be appropriate for the tidally locked Pluto-Charon system, or for certain close binary stars. In any such disjoint duplex, the two spheres attract each other with a gravitational force GMaMb/d2. In order to maintain their relative positions, the entire duplex must rotate as a whole with an qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 angular speed x ¼ G½M a þ M b =d about an axis perpendicular to the line connecting the centers of both spheres and through their mutual center of mass. Note that this relation is equivalent to Kepler’s third law. If the two spheres just touch, as in Fig. 3b, they are considered to be welded together at the point of tangency. Marchis et al. (2005) refer to such a configuration as a ‘‘non-detached peanut’’. German and Friedlander (1991) used a similar model they called a ‘‘dumbbell’’ consisting of two identical tangent spheres (although their abstract misleadingly states that they used two ellipsoids). Hudson and Ostro (1995) also used two spheres in contact as a preliminary model for Asteroid 4179 Toutatis. As before, the two spheres of a duplex attract each other with a gravitational force GMaMb/d2, where now d = Ra + Rb if the spheres are tangent from outside. If such a duplex is rotating with angular speed x about an axis perpendicular to the line through the centers of both spheres, the weld at the point of tangency also is subM M

ject to a centrifugal tension lx2d, where l  MaaþMbb is known as the ‘‘reduced mass’’ of the duplex, as for any two-body system. The net force on the weld is the difference between these two forces; if they cancel exactly, we recover Kepler’s third law above. If the two spheres overlap, as in Fig. 3c, they are presumed to interpenetrate in a lenticular ‘‘crush zone’’. Note that the density inside this lens must be the sum of the densities of both spheres separately, again to preserve the simplicity of the model. Marchis et al. (2005) refer to a similar model as ‘‘imbricated’’, or as the ‘‘snowman’’ configuration. The two spheres even may be nested, as in Fig. 3d, to represent a spherical body with a mantle and offset core of different densities. If the core is denser than the mantle (as in the Earth, for example), then the inner sphere has a positive mass and density. Conversely, if the core is less dense than the mantle (as in an icy comet surrounded by a silicate crust, for example), then the inner sphere must have a negative mass and density. Certain nested configurations even can contain void spaces. For example, if both spheres are homogeneous, and the density of the inner sphere is exactly the opposite of the outer sphere’s density, the duplex will contain a spherical void. Interestingly, such a void

contains a gravitational field of constant acceleration, like that described in Section 3.1 above (see Section 3.52 of Ramsey (1940)). A.2. Physical properties The equations of motion for any rigid body are simplest when formulated with respect to its center of mass. Without loss of generality, we assume Ra P Rb, and that the net local density q is nowhere negative. Then Ma and the total mass Ma + Mb are both positive (even in the special case when Mb is negative), and the center of mass of the duplex lies along its symmetry axis, at a distance

da ¼

jM b jd Ma þ Mb

ð23Þ

from the center of a, and at a distance

db ¼

Ma d Ma þ Mb

ð24Þ

from the center of b. An interesting case arises when the two spheres overlap like a ‘‘snowman’’ (Ra  Rb < d < Ra + Rb), as in Fig. 3c. Then their surfaces intersect at an exterior angle 2

d  R2a  R2b d ¼ Arccos 2Ra Rb

! ð25Þ

(from the law of cosines), forming a concave circular ‘‘neck’’ or ‘‘waist’’ of diameter

W¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 4 2R2a R2b þ 2R2a d þ 2R2b d  R4a  R4b  d =d

ð26Þ

(from Heron’s formula). Note that W attains its maximum width of qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R2a  R2b and d = Arccos(Rb/Ra). Note also that if

2 Rb when d ¼

rotation is neglected, the gravitational potential at the surface of the duplex attains its global minimum of GMa/Ra  GMb/Rb along the waistline. Meanwhile, the total outer surface area U and volume V of this conjoint duplex work out to

h i h i h i 2 2 U ¼ 2p R2a þ R2b þ p R2a  R2b þ d Ra =d þ p R2b  R2a þ d Rb =d ð27Þ and

h i p h 2 2 i2 p h 2 2i 3 V ¼ 2p R3a þ R3b =3 þ Ra  Rb =d þ R þ Rb d  pd =12: 4 2 a ð28Þ These two formulae may be useful in computing mean crater den  ðM a þ Mb Þ=V of a contact binary sity and mean mass density q asteroid, for example. In the limit of outer tangency (d = Ra + Rb), as in Fig. 3b, both d and W reduce to zero, while Eqs. (27) and (28) above simplify to h i h i U ¼ 4p R2a þ R2b and V ¼ 4p R3a þ R3b =3. The latter two formulae also apply when both spheres are completely detached, as in Fig. 3a. In the limit of inner tangency (d = Ra  Rb), W again reduces to zero, but now d becomes 180°, while Eqs. (27) and (28) simplify to just U ¼ 4pR2a and V ¼ 4pR3a =3. These last two formulae also apply when b lies entirely inside of a, as in Fig. 3d. For any rigid body (with positive total mass), its moment of inertia in a given direction is least about an axis through its center of mass. The moment of inertia of any duplex is

Ik ¼ Ia þ Ib

ð29Þ

A.R. Dobrovolskis, D.G. Korycansky / Icarus 225 (2013) 623–635

633

about the line containing the centers of both a and b. By symmetry, this is also a principal axis of inertia for the duplex. Furthermore, any line through the center of mass and perpendicular to the symmetry axis is also a principal axis, about which the moment of inertia is 2

2

2

I? ¼ Ik þ Ma da þ M b db ¼ Ik þ

Ma Mb d 2 ¼ I k þ ld : Ma þ Mb

ð30Þ

In our notation, Ik and I\ are always positive, but both Ib and l have the same sign as Mb, while l < Ma and l < Mb whether Mb is positive or negative. Internal dissipation of rotational energy eventually drives an imperfectly rigid body to spin about an axis of greatest inertia through its center of mass. Normally, the center of mass of a duplex lies between the centers of the two spheres, and I\ > Ik, so that the axes of greatest inertia are perpendicular to the symmetry axis. Then setting A = Ik, B = C = I\, and putting the symmetry axis along the x-axis, MacCullagh’s approximation (1) for the potential of the duplex reduces to 2

/¼

G½M a þ Mb  Gld  ½2x2  y2  z2 : r 2r 5

ð31Þ

However, in the exceptional nested case when Mb < 0 and Ib < 0 to account for a sparser core, the centers of both spheres lie on the same side of the center of mass, and I\ < Ik, so that the axis of symmetry is also the axis of greatest inertia. Then setting A = B = I\, C = Ik, and putting the symmetry axis along the z-axis, MacCullagh’s approximation (1) reduces to 2

/¼

G½M a þ Mb  Gld  ½2z2  x2  y2 : r 2r 5

ð32Þ

Fig. 4. Normalized gravitational potential of a duplex consisting of two positive point masses. The large X at (1/3, 0) denotes the location of a particle of mass 2, while the small x at (2/3, 0) denotes a particle of mass 1. The center of mass of the duplex lies between them at the origin (0, 0), while the horizontal dotted line indicates the axis of symmetry. The upper half of the figure plots contours of the exact potential from Eq. (35), while the lower half contours the quadrupole approximation from formula (31). Contour values: 20, 10, 8, 6, 4.5, 3.5, 3, 2.5, 2, 1.5, 1, 0, +20.

Note that both Ia and Ib cancel out of Eqs. (31) and (32) above. A.3. Gravitational field The exact gravitational field of a duplex is quite simple: it is just the sum of the fields of its two component spheres. Indeed, this is the greatest virtue of the duplex model. The external gravitational potential of sphere a is just

Ua ¼ GM a =r a ;

ð33Þ

where ra is the distance of the field point from the center of a. For completeness, the internal potential of sphere a is

Ua ¼

i GMa h 2 2 r a =Ra  3 ; 2Ra

ð34Þ

provided that its density is uniform. The gravitational potential Ub of sphere b is given by the analogous expressions, simply replacing the subscripts a by b. Then the total external potential of a duplex is just

U ¼ GMa =r a  GMb =rb :

ð35Þ

To illustrate, Fig. 4 plots equipotential contours of formula (35) above for a duplex with Ma = 2 and Mb = 1, in units where both G and the separation d between a and b are normalized to unity. The total mass of the duplex is 3, while its reduced mass l is only 2/3. Because the actual radii of a and b are irrelevant to the external field, we set them both to zero without loss of generality, and take both spheres as point masses for simplicity. Then Ia = 0 = Ib and A = Ik = 0 from Eq. (29), while B = C = I\ = ld2 = 2/3 from Eq. (30).pThen ffiffiffiffiffiffiffiffiffiffiffifrom Eq. (3), the equivalent radii of the duplex become a ¼ 10=9d  1:054d, b = 0, and c = 0. The top half of Fig. 4 shows the exact duplex potential from Eq. (35), while the bottom half plots the quadrupole approximation from formula (31). Note how the outer contours match relatively

well, producing somewhat egg-shaped equipotential surfaces approaching spheres centered about the origin in the far field (the 1 contour is just visible in the corners of the figure). In the near field, the exact equipotentials split in two and again become egg-shaped, approaching spheres centered around each mass. In contrast to the outer equipotentials, though, the inner contours of the quadrupole approximation match the exact potential poorly, especially on the right-hand side near the lesser mass; this is partly due to the artificial left–right symmetry of the quadrupole equipotentials. Also in the near field, note the teardrop-shaped white zone of positive potential just below the center of mass, aligned with the axes of greatest inertia. This too is a non-physical artifact of the quadrupole approximation, but it only extends out about 1/3 unit from the origin. However, an equally non-physical zone of repulsive gravity also can be seen extending out almost twice as far below the center. Of course, these artifacts are mitigated in the case of a physical duplex, where both spheres have positive radii. In this case, it is easy to find the gravitational potential, acceleration, and slope over the entire surface of the duplex from the exact formula (35) for the external potential of the duplex. Note that Eq. (35) is simpler than its point quadrupole approximation, Eq. (31). Furthermore, consideration of numerical errors indicates that Eq. (35) is also more accurate than Eq. (31), for distances less than some 105d. A.4. Spherical harmonics It is instructive to express the duplex potential (35) in terms of spherical harmonics, as usual for planetary gravity. Recall that the external gravitational field of sphere a is the same as that of a point particle of the same mass Ma, located at the center of a. For simplicity, let this particle lie along the ‘‘northern’’ axis of spherical

634

A.R. Dobrovolskis, D.G. Korycansky / Icarus 225 (2013) 623–635

coordinates R,k,h, at a distance R = da from the origin (here taken as the center of mass of the duplex). Then its latitude k is 90°, while its longitude h is arbitrary. In these coordinates, the external potential (33) of sphere a is

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

Ua ¼ GMa = R2 þ d2a  2Rda sin k

ð36Þ

(Kaula, 1968, p. 66). Note that the longitude h does not even appear in formula (36) above, because of our choice of poles along the symmetry axis of the duplex. The analogous expression Ub for sphere b is obtained from Eq. (36) by replacing Ma with Mb, but replacing da with db (or +db in the special case when Mb is negative). The spherical harmonic expansion of the point mass potential (36) is

Ua ðR; k; hÞ ¼ 

ð37Þ

(Kaula, 1968, p. 66), where Pn() is the Legendre polynomial of degree n. Again, Ub may be expanded similarly. For the duplex as a whole, the total potential is just the sum Ua + Ub, so the total potential contains two terms for each Legendre polynomial. Note that formula (37) above applies only outside of the sphere R = da; for R 6 da, this series expansion fails to converge in general, so da is its radius of convergence. Likewise the corresponding expansion for Ub diverges for R 6 db, so the radius of convergence for the entire duplex is the greater of da and db. Note that if da – db, it is possible to reduce the radius of convergence to a minimum of d/2 = da/2 + db/2, by placing the origin at the midpoint between the individual centers of a and b. Then the radius of convergence also lies entirely inside the duplex, provided that the two spheres over2 lap such that d 6 R2a þ R2b (so d P 90° from Eq. (25)). Other simple choices of origin would be the center of either sphere, their point of tangency (if just touching, as in Fig. 3b), or the ‘‘center of figure’’ of the whole duplex, here defined as the midpoint of its longest axis. For a body of arbitrary shape, its center of figure can be defined as the center of the smallest sphere enclosing it (often called its Brillouin sphere); sometimes this choice of origin may be the best alternative to its center of mass (which may be unknown anyway). Of course, the axis of symmetry runs through all six points mentioned above: the center of mass, the center of figure, the center of a, the center of b, the midpoint between the latter two, and the point of tangency if a and b just touch. Note that all three of the duplexes studied in the main text are of the latter type, and centered at the point of tangency, for convenience; but this does not coincide with either the center of figure or the midpoint between the centers of a and b for the second duplex, nor with the center of mass for the third. It is illuminating to interpret expansion (37) by looking at its first few terms separately. Since P0(z) = 1 for all z, to zeroth degree Eq. (37) simply reduces to the monopole potential U0(R, k, h) = G Ma/R, as if the offset da of a from the origin of reference were zero. Then adding the contributions of both spheres of the duplex gives

U0 ðR; k; hÞ ¼ G½Ma þ Mb =R;

ð38Þ

where the subscript zero now refers to the degree of the terms. Thus the monopole potential of the duplex is the same as if both spheres were centered at the origin. For the first-degree terms of formula (37), P1(sin k) = sin k, so the a ðda =RÞ sin k. Analdipole contribution U1 of sphere a becomes  GM R GM b R

ðdb =RÞ sin k if

GM  Rb

ðþdb =RÞ sin k if Mb < 0. In either case, the total diMb > 0, or pole potential becomes

U1 ðR; k; hÞ ¼ G½Ma da  jMb jdb  sinðkÞ=R2 :

a contribution U2 of sphere a as  GM ðda =RÞ2 P2 ðsin kÞ, where R

2

P2 ðsin kÞ ¼ 32 sin ðkÞ  12 ¼ 14  34 cosð2kÞ. The analogous contribution GM

of sphere b is  R b ðdb =RÞ2 P2 ðsin kÞ. Then the total quadrupole potential of the duplex is

U2 ðR; k; hÞ ¼ GMa d2b P2 ðsin kÞ=R3  GM b d2a P2 ðsin kÞ=R3 

n 1  GMa X da Pn ðsin kÞ R n¼0 R

ogously, the dipole contribution of sphere b is 

Formula (39) above is valid in general; but in the usual case where the origin lies at the center of mass, Mada  jMbjdb = 0 from Eqs. (23) and (24), and we recover the familiar fact that the gravitational dipole moment of any body vanishes when evaluated relative to its center of mass. Any other choice of origin introduces a dipole term in the harmonic expansion. The second-degree terms of formula (37) give the quadrupole

ð39Þ

2 Mb d P2 ðsin kÞ=R3 Ma þ Mb  2 Ma d  GMb P2 ðsin kÞ=R3 Ma þ Mb " # 2 Ma Mb d P2 ðsin kÞ=R3 ¼ G Ma þ Mb ¼ GMa

2

¼ Gld P2 ðsin kÞ=R3 ;

ð40Þ

where again we have substituted Eqs. (23) and (24) for da and db, and l for the reduced mass. M M b d2 2 Note that the factor ld ¼ Maa þM appearing in formula (40) b above is just the moment of inertia difference I\  Ik from Eq. (30), which shows the relation of the quadrupole potential to MacCullagh’s formula (1) from the main text. In terms of the conventional ‘‘dynamical oblateness coefficient’’ J2 and its reference I I

M a M b d2

radius S, we have S2 J 2 ¼ Mka þM? b ¼ ½M

a þM b 

2

. The sign of the above for-

mula indicates that the duplex is actually prolate with respect to its symmetry axis (or oblate in the special case when one mass is negative). If instead we prefer to situate both masses in the equator plane of our spherical coordinates, we may use the addition theorem of spherical harmonics (Kaula, 1968, pp. 67 and 124) to transform the quadrupole coefficient into the zonal harmonic C2,0   J2, along with the sectoral harmonic C2,2 characterizing the longitude dependence

of

þM M d

the

potential:

then

þM M d2

S2 J 2 ¼ 2½M aþMb a

b

2

and

2

S2 C 2;2 ¼ 4½M aþMb 2 . The addition theorem can be applied similarly a

b

to any other harmonic coefficient. References Ashenberg, J., 2007. Mutual gravitational potential and torque of solid bodies via inertial integrals. Celest. Mech. Dynam. Astron. 99, 149–159. Beletskii, V.V., 1965. Motion of an Artificial Satellite about its Center of Mass. NASA Technical Translation F-429. Benz, W., Asphaug, E., 1999. Catastrophic disruptions revisited. Icarus 142, 5–20. Borderies, N., 1978. Mutual gravitational potential of N solid bodies. Celest. Mech. 18, 295–307. Boué, G., Laskar, J., 2009. Spin axis evolution of two interacting bodies. Icarus 201, 750–767. Chauvineau, B., Farinella, P., Mignard, F., 1993. Planar orbits about a triaxial body: Application to asteroidal satellites. Icarus 105, 370–384. Danby, J.M.A., 1962. Fundamentals of Celestial Mechanics. MacMillan, New York. Dobrovolskis, A.R., 1995. Chaotic rotation of Nereid? Icarus 118, 181–198. Dobrovolskis, A.R., 1996. Inertia of any polyhedron. Icarus 124, 698–704. Elipe, A., Palacios, M., Pre ß tka-Ziomek, H., 2008. Equilibria of the three-body problem with rigid dumb-bell satellite. Chaos Solitons Fract. 35, 830–842. German, D., Friedlander, A.L., 1991. A simulation of orbits around asteroids using potential field modeling. Adv. Astronaut. Sci. 75, 1183–1201. Hirabayashi, M., Morimoto, M.Y., Yano, H., Kawaguchi, J., Bellerose, J., 2010. Linear stability of collinear equilibrium points around an asteroid as a two-connectedmass: Application to fast rotating Asteroid 2000 EB14. Icarus 206, 780–782. Hudson, R.S., Ostro, S.J., 1995. Shape and non-principal axis spin state of Asteroid 4179 Toutatis. Science 270, 84–86. Jackson, J.D., 1975. Classical Electrodynamics, second. ed. Wiley, New York.

A.R. Dobrovolskis, D.G. Korycansky / Icarus 225 (2013) 623–635 Kaula, W.M., 1968. An Introduction to Planetary Physics. Wiley, New York. Klemperer, W.B., Baker Jr., R.M., 1957. Satellite librations. Astronaut. Acta 3, 16–27. Korycansky, D.G., 2004. Orbital dynamics for rigid bodies. Astrophys. Space Sci. 291, 57–74. Korycansky, D.G., Asphaug, E., 2006. Polyhedron models of asteroid rubble piles in collision. Icarus 181, 605–617. Korycansky, D.G., Asphaug, E., 2009. Low-speed impacts between rubble piles modeled as collections of polyhedra. Icarus 204, 316–329. Leinhardt, Z.M., Richardson, D.C., 2002. N-body simulations of planetesimal evolution: Effect of varying impactor mass ratio. Icarus 159, 306–313. Leinhardt, Z.M., Richardson, D.C., Quinn, T., 2000. Direct N-body simulations of rubble pile collisions. Icarus 146, 133–151. Marchis, F., Hestroffer, D., Descamps, P., Berthier, J., Laver, C., de Pater, I., 2005. Mass and density of Asteroid 121 Hermione from an analysis of its companion orbit. Icarus 178, 450–464. Michel, P., Benz, W., Richardson, D.C., 2004. Catastrophic disruption of asteroids and family formation: A review of numerical simulations including both fragmentation and gravitational reaccumulations. Planet. Space. Sci. 52, 1109– 1117. Moran, J.P., 1961. Effects of plane librations on the orbital motion of a dumbbell satellite. Am. Rocket Soc. J. 31, 1089–1096. Murray, C.D., Dermott, S.F., 1999. Solar System Dynamics. Cambridge University Press, Cambridge, UK. Paul, M.K., 1988. An expansion in power series of mutual potential for gravitating bodies with finite sizes. Celest. Mech. 44, 49–59. Prieto-Llanos, T., Gómez-Tierno, M.A., 1994. Stationkeeping at libration points of natural elongated bodies. J. Guid. Control Dynam. 17, 787–794. Ramsey, A.S., 1940. An Introduction to the Theory of Newtonian Attraction. Cambridge U. Press.

635

Richardson, D.C., Elankumaran, P., Sanderson, R.E., 2005. Numerical experiments with rubble piles: Equilibrium shapes and spins. Icarus 173, 349–361. Robinson, W.J., 1973. The restricted problem of three bodies with rigid dumb-bell satellite. Celest. Mech. 8, 323–330. Roig, F., Duffard, R., Penteado, P., Lazzaro, D., Kodama, T., 2003. Interacting ellipsoids: A minimal model for the dynamics of rubble-pile bodies. Icarus 165, 355–370. Rosner, H.R., 1960. Motion of a dumbbell satellite. Adv. Astronaut. Sci. 7, 125–137. Schindler, G.M., 1960. Satellite librations in the vicinity of equilibrium solutions. Astronaut. Acta 6, 233–240. Schutz, B.E., 1981. The mutual potential and gravitational torques of two bodies to fourth order. Celest. Mech. 24, 173–181. Tricarico, P., 2008. Figure–figure interaction between bodies having arbitrary shapes and mass distributions: A power-series approach. Celest. Mech. Dynam. Astron. 100, 319–330. Walsh, K.J., Richardson, D.C., 2006. Binary near-Earth asteroid formation: Rubble pile model of tidal disruptions. Icarus 180, 201–216. Wang, Li.-S., Krishnaprasad, P.S., Maddocks, J.H., 1991. Hamiltonian dynamics of a rigid body in a central gravitational field. Celest. Mech. Dynam. Astron. 50, 349– 386. Wang, Li.-S., Maddocks, J.H., Krishnaprasad, P.S., 1992. Steady rigid-body motions in a central gravitational field. J. Astronaut. Sci. 40, 449–478. Werner, R.A., 1994. The gravitational potential of a homogeneous polyhedron or don’t cut corners. Celest. Mech. Dynam. Astron. 59, 253–278. Werner, R.A., Scheeres, D.J., 1996. Exterior gravitation of a polyhedron derived and compared with harmonic and mascon gravitation representations of Asteroid 4769 Castalia. Celest. Mech. Dynam. Astron. 65, 313–344. Werner, R.A., Scheeres, D.J., 2005. Mutual potential of homogeneous polyhedra. Celest. Mech. Dynam. Astron. 91, 337–349.