Stability of binaries. Part II: Rubble-pile binaries

Icarus 277 (2016) 125–140

Contents lists available at ScienceDirect

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

Stability of binaries. Part II: Rubble-pile binaries Ishan Sharma a,b,∗ a b

213B Northern Laboratories I, Department of Mechanical Engineering, IIT Kanpur, Kanpur 208016, UP, India Mechanics and Applied Mathematics Group, IIT Kanpur, Kanpur 208016, India

a r t i c l e

i n f o

Article history: Received 16 May 2015 Revised 10 April 2016 Accepted 12 April 2016 Available online 30 April 2016 Keywords: Asteroids Dynamics Satellites of asteroids Rotational dynamics Asteroids Composition Near-Earth objects

a b s t r a c t We consider the stability of the binary asteroids whose members are granular aggregates held together by self-gravity alone. A binary is said to be stable whenever both its members are orbitally and structurally stable to both orbital and structural perturbations. To this end, we extend the stability analysis of Sharma (Sharma [2015] Icarus, 258, 438–453), that is applicable to binaries with rigid members, to the case of binary systems with rubble members. We employ volume averaging (Sharma et al. [2009] Icarus, 200, 304–322), which was inspired by past work on elastic/fluid, rotating and gravitating ellipsoids. This technique has shown promise when applied to rubble-pile ellipsoids, but requires further work to settle some of its underlying assumptions. The stability test is finally applied to some suspected binary systems, viz., 216 Kleopatra, 624 Hektor and 90 Antiope. We also see that equilibrated binaries that are close to mobilizing their maximum friction can sustain only a narrow range of shapes and, generally, congruent shapes are preferred.

1. Introduction In this paper we analyze the stability of rubble binary systems. This work constitutes the second part of a two-part paper. In part I (Sharma, 2015, henceforth S2015) the stability analysis was restricted to binaries with rigid members. The motivation for the current work lies both in the discovery of an ever-increasing number of binary systems amongst the minor planet population, and the fact that a fair number of them are suspected have members that are granular aggregates; see e.g., Richardson and Walsh (2006) and Merline et al. (2002). This latter belief rests upon the low estimated densities of these objects – for example (Abe, et al., 2006) find Itokawa’s density to be 1.9 g cm−3 – and the hypothesis that binaries are often the result of accretionary processes; see Walsh and Richardson (2006), Walsh et al. (2008) and Michel et al. (2001, 2004). S2015 expands further upon the above motivations. The equilibrium of rubble binary systems was investigated by Sharma (2010), henceforth S2010. This constituted an appropriate extension of the second Darwin problem for fluid binaries to granular aggregates. Stability analyses of binaries has typically focussed on orbital stability, ignoring the structural response of the binary’s members. In Solar System studies, Scheeres (2007, 2009) were two of the first works to include structural considerations in the dy-

∗ Correspondence to: 213B Northern Laboratories I, Department of Mechanical Engineering, IIT Kanpur, Kanpur 208016, UP, India. Tel.: +91 5122596152. E-mail address: [email protected]

http://dx.doi.org/10.1016/j.icarus.2016.04.017 0019-1035/© 2016 Elsevier Inc. All rights reserved.

© 2016 Elsevier Inc. All rights reserved.

namics of binaries with rigid members. These papers investigated, respectively, rotational fission of contact binaries, and the stability of the full1 planar two-body problem. In astrophysical applications, structural stability of ellipsoidal fluid binaries has, however, been investigated; see Chandrasekhar (1969, Ch. 8). This classical body of work was revisited by Lai et al. (1993, 1994) who employed energy minimization to probe the stability of compressible inviscid fluid Roche and Roche-Riemann ellipsoids. Recently, Descamps (2015) has considered the equilibrium figures of “dumb-bell” like fluid objects for application to contact binary asteroids. In the context of granular bodies, Goldreich and Sari (2009) studied tidal evolution of rubble piles through scaling arguments and order-of-magnitude approximations. They considered structural failure very briefly – without introducing a yield criterion, as is usually done for granular materials – to estimate the critical diameter for rubble bodies in space. Only orbital stability of such systems was probed, and structural stability was not investigated. Finally, the effect of the body’s shape on stability was ignored. Scheeres (2012) explored the equilibrium and stability of rubble bodies by minimizing the total energy of a collection of self-gravitating, smooth, rigid spheres. However, rate-independent Coulomb friction, which plays an essential role in the constitutive response of granular aggregates, was not explicitly included in the energy minimization. Thus, it is possible that several equilibria were missed during energy minimzation. Indeed, Holsapple

1

This indicates that each body’s mass moment of inertia is accounted for.

126

I. Sharma / Icarus 277 (2016) 125–140

Fig. 1. The general configuration of binary system with ellipsoidal members. The centers of mass of the binary, its primary, and its secondary are identified by C, P, and S, respectively.

(20 01), Harris et al. (20 09) and Sharma et al. (2009) obtained many more equilibrium shapes for frictional rubble asteroids than are possible for fluid asteroids. Finally, as Hill (1957), Chakrabarty (1969), Storåkers (1977) and Sharma (2012) show, stability of rateindependent materials is fundamentally different from that of energy preserving or rate-dependent materials. Here, we investigate the stability of granular binaries. We will, for the first time, include the effect of dry friction, shape, structural deformation, and three-dimensional and finite perturbations on the stability of a binary. The constitutive response, developed through plasticity theory in Section 3, explicitly incorporates rateindependent frictional behavior. A binary system will be considered stable only if the orbit and the structure of each member is stable to both orbital and structural perturbations. The structure designates the collection of material points that constitute a binary member, and structural stability refers to this collective staying close to its equilibrium configuration. We discuss stability more precisely in Section 5. We limit ourselves to affine velocity perturbations of the equilibrium state.2 This consists of homogeneous structural deformations of a member about its possibly moving mass center. The present framework may be utilized for the stability analysis of any binary system with ellipsoidal members moving about the binary’s mass center on circular orbits. However, for simplicity, final calculations will be restricted to binaries with prolate members. We first derive relevant dynamical equations for a binary system. We generally follow the notation of S2015. A short primer on tensors is available in Sharma (2009, App. A), and more information in texts such as Knowles (1998). 2. Binary dynamics A binary system will consist of two ellipsoids – a primary of mass m and a secondary of mass m  – orbiting about the common mass center; see Fig. 1. Three coordinate frames will be employed, 2

This approximation is motivated by work in elasticity and fluid mechanics. Love (1946, Ch. XI) showed that shape of an elastic, self-gravitating, rotating ellipsoid changes homogeneously. At the same time, there is a large body of research on the shapes, dynamics and stability of rotating, self-gravitating, fluid ellipsoids that successfully assumed homogeneous deformation; see Chandrasekhar (1969) for details and older references. Newer work addresses compressible, polytropic fluids (Lai et al., 1993), tidal flybs (Sridhar and Tremaine, 1992), and close binaries (Lai et al., 1994). The extent to which rubble-pile binaries will deform homogeneously remains to be fully tested against simulations, although past work on tidal flybys of rubble asteroids (Sharma et al., 2006) obtained good match with discrete-element (DE) simulations of Richardson et al. (1998), while assuming that asteroids deformed in an affine manner. Similarly, Holsapple (2010) assumed homogenous deformations to study YORP-induced shape changes in rubble asteroids and achieved results comparable to the DE simulations of Walsh et al. (2008). See also the last paragraph of Section 3.



viz., the principal axes system (P, eˆ i ) and (S, eˆ i ) of, respectively, the primary and the secondary, and the frame (O, eˆ i ), in which we will investigate stability.3 The frames P, S and O are taken to be rotating at, respectively, ωP (t ), ωS (t ) and ω (t ). Fig. 1 shows various other geometric parameters that will be utilized. In particular, we note that the mass centers of the primary and the secondary are separated by R. Further, rS and RS = Reˆ S , locate the secondary with respect to, respectively, the primary and the binary’s center of mass. The vectors rP and RP have similar interpretations. We recall that every angular velocity ω (t ) has a corresponding anti-symmetric angular-velocity tensor Ω (t ), so that, for any vector x,

ω × x = Ω · x;

(1)

ω is the axial vector of Ω. Unless otherwise stated, all time derivatives will henceforth be with respect to an observer in the rotating frame O. We now investigate the dynamics of a binary whose members execute affine motions. 2.1. Affine motion During affine motion each member deforms homogeneously with respect to its mass center. The center of mass itself may move rectilinearly. Relative to a rotating frame O, the material velocity within each member may then be written as vP = r˙ P + x˙ P = r˙ P + LP · xP

and

vS = r˙ S + x˙ S = r˙ S + LS · xS ,

(2)

where L is a member’s velocity gradient measured in O, and x is the position of a material point within a member with respect to its mass center see Fig. 1. The tensor L estimates the local rate of change of relative displacement between two neighboring material points. It is helpful to split L as

L = D + W,

(3)

where the symmetric strain/stretching-rate tensor D and the antisymmetric spin-rate tensor W capture local deformation and rotation rates, respectively. Thus, if D were zero, we would recover rigid-body motion of the kind considered in S2015 with W playing the role of an angular-velocity tensor associated with a member’s angular velocity. It is possible to consider structural and orbital motions more or less separately, as we do next.

3 The frame O’s origin is at the barycenter C. However, as shown in Fig. 1, O may not be aligned with the typical barycentric frame, one of whose axis is typically normal to the orbital plane, while the other points from the primary to the secondary; cf. O’s definition in Section 5.

I. Sharma / Icarus 277 (2016) 125–140



2.2. Structural motion During affine motion the primary and the secondary deform homogeneously with respect to their respective mass centers. During such motions, initially ellipsoidal objects remain ellipsoidal, though they may rotate, stretch, shrink and shear. S2010 derives dynamical equations governing such structural motions. These may be written in the frame O as, respectively,





˙ + Ω 2 + 2 Ω · LP ) · I P L˙ P + L2P · I P = −σ P V  + M TP − (Ω I˙ P = LP · I P + I P · LTP ,

and and





˙ + Ω 2 + 2 Ω · LS ) · I S L˙ S + L2S · I S = −σ SV  + M TS − (Ω

127

(4a)

(4b)

(5a)



M QP = −2π ρ  GI P · BS(0) + eˆ P · BS(1)







(9)

 (1 )

M QS = −2π ρ  GI S · BP(0) + eˆ S · BP

and

,

(10)

where B(0 ) and B (1 ) are, respectively, second- and third-order tensors. These tensors are obtained by expanding the tidal shape tensor B associated with an ellipsoid about its companion’s mass center, and S2015 provides further details. The total moment on a member is obtained as a sum of the moments due to internal and external gravity. Thus, the total moment acting on the primary and the secondary are, respectively,

˜S M P = −2π ρ  GI P · AP − 2π ρ  GI P · B

(11a)

˜ P, M S = −2π ρ  GI S · AS − 2π ρ  GI S · B

and

(11b)

where

I˙ S = LS · I S + I S · LTS ,

and

(5b)

where the dot (·) indicates time derivative with respect to an ob server in O, σ = V σ dV/V is the volume-averaged stress tensor within a binary member,



I=

V

ρ x  xdV 

and

M=

V

(6a)

˜ S = B(0) + eˆ P · B (1) B S S

(6b)

are, respectively, a member’s inertia tensor and the external moment tensor acting on it, and ρ and V are the density and volume of a member, respectively. We note that in (4a) and (5a), the last three terms on the right-hand side are, respectively, angular, centripetal and Coriolis’ accelerations, and act as external moment tensors in the rotating frame. The first of (4) and (5) follow, respectively, the evolution of LP and LS in O by balancing internal stresses, external moments and inertial effects. The second of (4) and (5) describe, respectively, the changing inertia tensors I P and I S . Together these two sets of equations follow the structural deformation relative to frame O of the binary system’s ellipsoidal primary and secondary members. Recalling (3), it is clear that these equations would match the corresponding ones of S2015 were the strain rates DS and DP zero. 2.3. Moments

The secondary and the primary are located at rS and rP , respectively. These are related to RP and RS by

M G = −2π ρ GI · A,

(7)

where A is the gravitational shape tensor that captures the effect of the ellipsoidal shape. The tensor A is diagonalized in the principal coordinate system of a member, and Sharma et al. (2009) provides its components. To compute the tidal moment tensor due to the action of one ellipsoid on another, we first note that the force per unit mass bQ exerted by a homogeneous ellipsoid of density ρ at an external point at X is bQ = −2π ρ GB · X,

m m + m

(8)

where B is the tidal shape tensor that contributes the influence of the ellipsoid’s shape. The tensor B is diagonalized in the principal coordinate system of the ellipsoid, and a detailed description is provided in S2015. Employing the above formula with definition (6), S2015 computes the tidal moment tensors due to the primary’s action on the secondary, and vice versa to be, respectively,

R S = μ R S

and

rP =

m m + m

RP = μ RP .

(13)

S2015 shows that orbital motions in the frame O of the primary and the secondary are governed by, respectively,

˙ · RP + 2Ω · R˙ P = ¨ P + Ω · RP + Ω R 2

1 FP m

(14a)

2 ˙ · RS + 2Ω · R˙ S = 1 FS , and R¨ S + Ω · RS + Ω m

(14b)

where

m = μ m = μ m =

m m , m + m

(15)

while FP = −2π ρ  Gm bS · Reˆ P

and

FS = −2π ρ  Gm bP · Reˆ S ,

(16)

are, respectively, the total force exerted by the secondary on the primary, and vice versa, in terms of

bS = BS(0) +

The moment tensor M includes the effect of a member’s own gravity and that of its companion. The self-gravity moment tensor for an ellipsoidal body is

(12)

2.4. Orbital motion

rS =

ρ x  bdV

˜ P = B(0) + eˆ S · B (1) . and B P P

 1  (1) BS + eˆ P · BS(2) : I P  eˆ P  2 mR

and bP = BP(0) +

(17a)

 1  (1) BP + eˆ S · BP(2) : I S  eˆ S , m R2

(17b)

where B(2 ) is a fourth-order tensor. This tensor is obtained after expanding B, as was done in (9) to find B(0 ) and B (1 ) , and more details are available in S2015. At equilibrium, RP = REP and RS = RES are fixed, and (14) reduces to FEP − μ m ΩE · REP = 0 2

and

FES − μ m ΩE · RES = 0, 2

(18)

where FEP and FES are obtained after evaluating (16) at equilibrium. It is important to realize here that there are several rotations involved in the current analysis: the orbital angular velocities of the mass centers of the primary and the secondary as they orbit the binary’s center of mass; the rotation rates of the principal axes of inertia of the two members; the material spins W of each binary member when they deform, as introduced in (3); and the rotation rate Ω of the coordinate system O in which we will pursue our stability analysis. All these rotation rates will differ in a general motion wherein the binary members are allowed to deform.

128

I. Sharma / Icarus 277 (2016) 125–140

3. Rheology We utilize the rheology introduced by Sharma (2013) that we summarize here for completeness. Rubble-piles are modeled as rigid-plastic cohesionless materials obeying the pressuredependent Drucker–Prager yield criterion and an associative flow rule. We enforce associativity, as non-associative flow rules are trivially secularly4 unstable. In associative materials the yield surface is also the plastic potential which generates the flow rule. This is not the case in non-associative materials; see, e.g., Chen and Han (1988). In our model, a granular aggregate remains rigid until it violates the Drucker–Prager yield condition

|s|2  k2 p2 ,

(19)

prescribed in terms of the pressure

1 p = − tr σ , 3

(20)

the magnitude

|s|2 = si j si j of the deviatoric stress tensor

s = σ + p1, √ 2 6 sin φF , 3 − sin φF

(22)

which is known in terms of the aggregate’s internal friction angle φ F . The internal friction angle depends on interfacial friction due to particle interaction, and a geometric friction due to interlocking and rearrangement of finite-sized grains. The latter is typically dominant in dense aggregates, but reduces with density. Similarly, φ F depends on confining pressure, and grain crushing lowers φ F when the confining pressure is 1 MPa or beyond. Dense soils at these pressures typically have 30o ≤ φ F ≤ 40o ; see Sharma (2013) for more detail. Post-yield, an associative flow rule relates stress and strain increments:

σ



D = −p(1 +  )1 + p k2 + 3 2 , |D|

(23)

where D is defined by (3) and  = k2 /9. The stress σ and strain rate D are called compatible whenever they satisfy (23), σ lies on the yield surface, and D is along the outward normal to the yield surface at σ ; see Fig. 2. During a stability analysis we will also investigate perturbations that impose strain rates incompatible to the stress state. This may mean that (a) the material’s stress state σ is at yield, but the perturbing D is not normal to the yield surface at σ or (b) the material’s stress state is below yield, but is subjected to a non-zero strain rate perturbation leading to plastic flow. In both cases, the material’s current stress state changes suddenly to one that is compatible with the imposed strain rate where (23) holds. This is explained in Fig. 2. The pressure p is related to the dilation rate

3|D| trD = √ , k2 + 3 2

(24)

by

p˙ = −κ trD,

(25)

where κ is the granular aggregate’s plastic bulk modulus, i.e., bulk modulus post–yield. Because the dilation rate is positive, the material expands after yielding causing a drop in pressure. Sharma 4

Fig. 2. The Drucker–Prager yield surface in principal– stress/strain space is a circle when viewed on a plane perpendicular to the pressure axis σ1 = σ2 = σ3 . The origin is at ‘O’. The point ‘E’ corresponds to an equilibrium stress state. The strain rate DE compatible to the stress state σ E is normal to the yield surface at ‘E’. If the incompatible strain rate D0 is superimposed, the stress state shifts immediately to a location σ 0 where the imposed strain rate will be compatible, introducing a stress jump σ . The equilibrium stress σ A of another configuration is not at yield but, when perturbed, may move to a state σ B with corresponding compatible strain rate DB ; cf. Section 5.2.

(21)

and the parameter

k=

DE

Systems found stable by the energy criterion are secularly stable, cf., Section 5.2.

(2013) discusses possible values of κ for aggregates such as sand in detail. We note that the bulk modulus for sands confined at hundreds of kilo-Pascals is of the order of 1 MPa. We assume isotropic homogeneous bodies and homogenous deformation. Therefore, strain rates and internal stress fields will be uniform. Consequently, post-yield, we will replace σ and p in the constitutive relations (23) and (25) by, respectively, their volume averages σ and p. The average density ρ and volume V of a homogeneously deforming body change as

V˙ ρ˙ = − = trD. V ρ

(26)

From the above, it is clear that we look for yielding on average, rather than locally. Clearly, the latter will precede the former, though it may not be as important. Local yielding should be contrasted from global collapse in which the body cannot support imposed loads and loses integrity. Most bodies yield locally, but continue to retain their shape and support loads, e.g., an internally pressurized cylinder will locally yield from its inner surface, but will collapse only when the yielded zone extends to its outer surface. When investigating shapes and stability of entire bodies (like asteroids), it appears more important to estimate when global collapse will occur; see, e.g., Holsapple (2001). Such estimates for rigid-plastic bodies may be arrived at through the application of theorems in limit analysis; see Chen and Han (1988, Ch. 8). Following such an analysis, Holsapple (2001) showed that a rubble, cohesionless asteroid (modeled as a rigid-plastic continuum) collapsed globally only when it yielded everywhere and, moreover, the stress field that led to global collapse was unique and coincided with the stress field that caused yielding on average. This result could be extended easily to rubble satellites (Holsapple and Michel, 20 06, 20 08; Sharma, 20 09) and rubble binaries (Sharma, 2010), again when modeled as rigid-plastic bodies. Thus, there appears to be close link between the ability of a rigid-plastic body to retain its integrity and shape (global collapse), and yielding on average. In fact, Holsapple (2008) proved that in any rigid-plastic body the stress field that produces yielding on average will also cause global collapse in the body. However, the accuracy to which yielding on average will capture the actual failure of a granular aggregate will depend on our success in modeling rubble-piles as

I. Sharma / Icarus 277 (2016) 125–140

129

Fig. 3. A binary with prolate members at equilibrium. The primary and the secondary are aligned along their long axes with their shortest principal axes eˆ 3 and eˆ 3 , respectively, normal to the orbit. The binary’s mass center lies at C. The coordinate system O in which stability is evaluated is shown, and is parallel to the primary’s principal axes P. The principal axes S of the secondary is rotated from P by 180° about eˆ 3 . The unit vector eˆ P (eˆ S ) locates the center of mass of the primary (secondary) with respect to that of the secondary (primary).

rigid plastic materials. This may only be ascertained by extensive finite-element or discrete-element (DE) simulations. 4. Equilibrium shapes The equilibrium configuration of the binary system is shown in Fig. 3 along with the coordinate systems P, S and O. The primary and the secondary orbit the binary’s center of mass on circular paths, while being tidally locked with their long axes parallel. Each member rotates at ωE about its shortest axis that is parallel to the normal eˆ 3 of the orbital plane. This is prompted by the prediction of Burns and Safronov (1973), Efroimsky and Lazarian (20 0 0), Sharma et al. (20 05) and Breiter et al. (2012) that energy dissipating rotors ultimately align into a state of pure rotation about their shortest axis. Orbital evolution of energy dissipating rubble binaries has been investigated recently by Goldreich and Sari (2009) and Jacobson and Scheeres (2011), although they did not consider the evolution of the orientation of the binary’s members. The rotation rate at equilibrium of the primary and the secondary equals their common orbital angular velocity ωE = ωE eˆ 3 . The rotation rate of the frame O at equilibrium is also taken to be ωE eˆ 3 , that may now be obtained from (16) and first of (18):

 m + m 1 ωE2 = 2π ρ  G eˆ P · BS(0) + 3 B˜ S  RP · eˆ P .  R

m

(27)

Analogously, from (16) and the second of (18) yields

 m + m 1 ωE2 = 2π ρ  G eˆ S · BP(0) + 3 B˜ P  RS · eˆ S .  R

m

(28)

The two formulae for ωE differ as tidal forces are approximated by a series expansion as detailed in Section 2.5 of S2010. In general, ωE is better approximated by the formula that utilizes the B of the more spherical of the two members. We emphasize that the difference in the above formulae arises because we approximate the tidal shape tensor B in (8) through a series consisting of B(0 ) , B (1 ) , B(2 ) , etc. This approximation is done to avoid computing B directly, which is computationally intensive. If we had retained BS and BP in their entirety when computing FS and FP , then, indeed FS and FP match5 as they should.

5 From Maciejewski (1995), the total force exerted by the primary on the sec  ondary, FS = −ρ  ρ  G V  V  XS −XP 3 dV  dV  , where X locates a material point in an

|XS −XP |

ellipsoid with respect to the mass center C in Fig. 1. The force FP has the same formula, but with minus sign, so that FS = −FP . From Chandrasekhar (1969, p.  48, Th. 7), −ρ  G V  XS −XP 3 dV  = −2πρ  GBP · XS , as it is the force per unit mass

At equilibrium the members of the binary appear stationary in the frame O, although this may not be the case following a perturbation; cf., Section 5. At equilibrium, the tensors LP and LS and their derivatives vanish, as does Ω’s derivative, so that (4a) and (5a) provide the average stress at equilibrium in the primary and the secondary as, respectively,

−σ P V  + M TP − ΩE · I P = 0 2

and

− σ SV  + M TS − ΩE · I S = 0, 2

(29) where, we recall, ΩE is the angular velocity tensor related to ωE through (1). From (29), and one of (27) or (28), we may estimate the volume-averaged stresses at equilibrium in the binary’s members as functions of their shape and their separation. These estimates are then combined with the yield criterion (19) to obtain equilibrium regions in the six-dimensional space spanned by the scaled separation q∗ = R/(a1 + a1 ), the ratios6 χ = m /m and η = ρ  /ρ  , and the shapes of the primary and the secondary characterized by the axes ratios α  = a2 /a1 , β  = a3 /a1 , α  = a2 /a1 and β  = a3 /a1 , wherein a rubble-pile binary may exist in equilibrium; see S2010. These equilibrium zones are parameterized by the friction angle φ F . This six-dimensional region is explored through appropriate planar sections. An example is displayed in Fig. 4 for the case of a binary with equally massive (χ = 1 ) prolate members (α  = β  , α  = β  ) separated by q∗ = 1.5. The equilibrium zones corresponding to various choices of φ F are delineated in the twodimensional β  − β  space. At the boundaries of these zones, either the primary or the secondary are at yield. The primary yields at the vertically aligned boundaries, while the secondary does so at boundaries that are horizontally oriented. The two members simultaneously yield at the meeting point of these boundaries, and these intersections lie along the plot’s diagonal passing through (0.1, 0.1). The binary may remain at equilibrium as long as the axes ratios of its prolate members are such that they lie within the region corresponding to the internal friction angle φ F of the binary. For example, the equal-mass rubble binary with β  = 0.7 and β  = 0.54, shown in Fig. 4 by a filled circle, will necessarily require an internal friction angle greater or equal to 10° to survive as a granular aggregate. Similarly, a binary with equal-mass members will necessarily have to lie within the shaded region in Fig. 4 were its φF = 5◦ , else it would not exist as a rubble-pile. Analogous equilibrium plots may be generated for other choices of q∗ , χ , α  , η and α   ; see S2010. We now probe the stability of equilibrated rubble-pile binaries.

|XS −XP |

exerted by the primary at the location XS with respect to the primary. Simi  larly, −ρ  G V  XP −XS 3 dV  = −2πρ  GBS · XP . Thus, FS = ρ  V  −2πρ  GBP · XS dV  , |XP −XS |    and FP = ρ  V  −2πρ  GBS · XP dV  , which are exactly the formulae employed by Sharma (2010, 2015). Clearly, FS = −FP .

6

In S2010 the ratio χ is denoted by κ . Here κ is the plastic bulk modulus.

130

I. Sharma / Icarus 277 (2016) 125–140

where

J S = (trI S )1 − I S

and J P = (trI P )1 − I P

(32)

inertia7

are the familiar Euler moment of of, respectively, the primary and the secondary about their own mass centers, while

j S = m r 2 1 − m rS  rS

and

j P = m r 2 1 − m rP  rP .

(33a)

The above are the moments of inertia about the binary’s mass center of, respectively, the secondary and the primary, considering them to be point masses. At equilibrium, the frame O rotates at the equilibrium rate ωE . We restrict ourselves to perturbations that do not affect Hrel , so that post-perturbation ω0 = ωE . The rate of change of the rotation rate of the frame O immediately after the perturbation is obtained by evaluating

J B · ω˙ = −J˙ B · ω − ω × (J B · ω )

(34)

just after the perturbation, i.e., at t = Fig. 4. Equilibrium landscape in β  –β   space for a binary with equal-sized members (χ = 1 ) separated by q∗ = 1.5. The equilibrium zones are parameterized by the binary’s internal friction angle φ F , which are indicated next to their associated critical curves. This plot follows Fig. 5 of S2010.

5. Stability S2015 discusses definition of stability due to Lyapunov, as well as, the spectral and energy stability tests. As per Lyapunov, ‘a system is locally stable if small perturbations of the system lead to small departures of its coordinates from its equilibrium values’; see LaSalle and Lefschetz (1961, p. 28). The spectral method decides stability by investigating the eigenvalues of the system’s linearized equations of motion. Thus, this method requires that the governing equations of the system be amenable to linearization. Consequently, it cannot be employed to investigate the stability of non-smooth systems, such as the one considered here. The energy method, on the other hand, predicts stability if the system’s equilibrium state lies at the minima of an appropriately defined potential energy. Depending on the method employed, a system is said to be, respectively, spectrally/dynamically or energetically/secularly stable or unstable. We consider a binary to be stable only if the structures and the orbits of both its members are stable to orbital and structural perturbations. Because the material model of Section 3 is non-smooth, we follow Sharma (2012, 2013) and employ an incremental formulation of the energy criterion. For this we consider the total kinetic energy, relative to an appropriate coordinate system, that is stored in the orbital and structural motions of the binary. S2015 suggests that an appropriate coordinate system O in which to consider the stability of a binary is one in which, relative to O, its total angular momentum

 Hrel =

V 

ρ  xS × x˙ S dV  +

 V

ρ  xP × x˙ P dV  + m rS × r˙ S + m rP × r˙ P (30)

vanishes. The origin of the frame O is attached to the center of mass of the binary. This definition generalizes the Tisserand’s mean axes of the body. We now list some of the formulae from Sharma (2015, Section 4.1) that will be relevant to the subsequent development. The total Euler moment of inertia of the binary about its center of mass is

J B = J S + J P + jS + jP ,

(31)

0+ :

  −1 ˙

˙ ω˙ |t=0+ = ω˙ 0 = − J −1 B · J B · ω t=0+ = −J B · J 0B · ω E ;

(35)

here and everywhere the subscript ‘0’ denotes evaluation at t = 0+ . Note that (34) follows from the conservation of the total angular momentum of the isolated binary; see Section 4.1 of S2015. We ¨ 0 , and this is obtained by differentiating (34), and also require ω evaluating the resultant at t = 0+ . We now formulate the energy criterion for a binary system. We begin by defining admissible perturbations in a manner similar to Section 4.2 of S2015. 5.1. Admissible perturbations Each binary member is restricted to affine motions, and so, too, will be their perturbations. Thus, structural perturbations of the secondary and the primary will be limited to homogeneous perturbations: x˙ 0S = LoS · xS

and x˙ 0P = LoP · xP ,

(36)

respectively, where LoS and LoP are arbitrary tensors. The center of mass of each member may be rectilinearly perturbed. Therefore, orbital perturbations of the secondary and the primary, with respect to the binary’s mass center C, will be r˙ oS = r˙ iS eˆ i and r˙ oP = r˙ iP eˆ i , respectively. We prefer to phrase orbital motions in terms of RS and RP . To this end, we invoke (13), to write R˙ oS = R˙ 1S eˆ 1 + R˙ 2S eˆ 2 + R˙ 3S eˆ 3

and

R˙ oP = R˙ 1P eˆ 1 + R˙ 2P eˆ 2 + R˙ 3P eˆ 3 . (37)

Affine velocity perturbations of the primary and the secondary are, respectively, voS = r˙ oS + x˙ oS = μ R˙ oS + W oS · xS + DoS · xS

(38a)

and voP = r˙ oP + x˙ oP = μ R˙ oP + W oP · xP + DoP · xP ,

(38b)

where, again, (13) is employed, and we have split LoS and LoP as in (3). The tensors Do and W o represent initial structural deformation, and rigid-body spin or misorientations, respectively. In addition to viewing the above perturbations as a combination of orbital and structural parts, it is helpful to view them as being composed of rigid-body modes, of the form r˙ + W · x, and strain-rate driven deformation modes, of the kind D · x. 7 The Euler moment of inertia J enters into the equations of motion governing the rotation of a rigid body; see, e.g., Greenwood (1988, p. 301). The inertia tensor I , defined by (6) as the moment of the mass distribution, arises naturally during volume averaging.

I. Sharma / Icarus 277 (2016) 125–140

Stability will be tested for perturbations (38) generated by all admissible choices of R˙ oS , R˙ oP , DoS , DoP , W oS and W oP . However, these perturbations cannot all be independently specified, and we now discuss the conditions that these perturbations must satisfy. First, because O’s origin is attached to C,

m r˙ oS + m r˙ oP = 0,

so that from (13), R˙ oS + R˙ oP = 0.

 



+ m rES × r˙ oS + m rEP × r˙ oP = 0,

where ‘sk’ represents the anti-symmetric part of a tensor, while ‘ax’ is the axial vector associated with an anti-symmetric tensor; cf. (1). The first term in the equation above represents the angular momentum stored in the structural motions of the secondary and the primary calculated about their respective mass centers. Employing (13) and its derivative, we rewrite the above equation as

 

−2 ax sk I S · LToS + I P · LToP



+m



   ˙ 0 + Ω 2 · I S : LT δ (1) Es = − σ SV  + M TS − Ω E oS and

 μ RES × R˙ oS + μ REP × R˙ oP = 0.

Utilizing RES + REP = 0 and (39), we may simplify the above considerably to

 



+ m RES × R˙ oS = 0.

Among all choices of perturbations R˙ oS , R˙ oP , DoS , DoP , W oS and W oP , only those that satisfy conditions (39) and (40) are admissible. The class of perturbations (38) so generated constitute admissible perturbations. 5.2. Energy criterion We follow S2015 to compute the total kinetic energy of a binary relative to the rotating frame O as

Ek = Es + Es + Eo + Eo , where Es and Es represent, respectively, the kinetic energy measured in O of the structural motions of the primary and the secondary relative to their corresponding mass centers, while Eo and Eo are the orbital kinetic energies, again as seen in O, carried by the centers of mass of the primary and the secondary, respectively. A binary is secularly stable as per the energy criterion if, over small time δ t, the change δ Ek in Ek is negative for all admissible perturbations (38), i.e.,

δ Ek =

 δt 0

E˙ k dt < 0.

(41)

We estimate δ Ek by expanding it in a Taylor series about t = 0+ :

δ Ek = δ ( 1 ) Ek δ t + δ ( 2 ) Ek

δt 2 2

+O



 δt 3 ,

(42)

where, for i = 1, 2, ...,

δ (i) Ek = δ (i) Es + δ (i) Es + δ (i) Eo + δ (i) Eo .

(43)

The terms on the right hand side above are evaluated as in Sharma (2013, 2015). For homogeneous structural motions the secondary’s structural energy is

δ Es =

 δt  0

− σ SV  + M TS −



 ˙ + Ω2 · I S : LT dt. Ω S

This is expanded in a Taylor series to yield8 8 We recall that the subscript ‘0’ refers to evaluation of a tensor at t = 0+ , while ‘E’ indicates the equilibrium value. With this, ΩE is the rotation rate tensor at equilibrium, which equals Ω0 as perturbations don’t change the binary’s angular mo˙ 0 is the angular acceleration tensor evaluated at t = 0+ . mentum. The term Ω

(44a)

   d  ˙ + Ω 2 · I S : LT . − σ SV  + M TS − Ω S 0 dt (44b)

Similarly for the primary, we find

   ˙ 0 + Ω 2 · I P : LT δ (1) Es = − σ PV  + M TP − Ω E oP and

δ (2) Es =

(45a)

   d  ˙ + Ω2 · I P : LT . (45b) − σ P V  + M TP − Ω P 0 dt

For the orbital kinetic energy terms, we obtain from S2015 that

δ (1) Eo = μ R˙ 0S · and

δ (1) Eo = μ R˙ 0P · and



FES − μ m



 Ω2E + Ω˙ 0 · RES

     δ (2) Eo = μ ddt R˙ S · FS − μ m Ω2 + Ω˙ · RS 0 ,

and

(40)

−2 ax sk I S · LToS + I P · LToP

δ (2) Es =

(39)

Thus, orbital perturbations of the secondary and the primary are not independent. Second, the structural and orbital perturbations must be such that they leave Hrel defined by (30) unchanged at zero. Employing (36), this leads to Hrel = −2 ax sk I S · LToS + I P · LToP

131



FEP − μ m



 Ω2E + Ω˙ 0 · REP

     δ (2) Eo = μ ddt R˙ P · FP − μ m Ω2 + Ω˙ · RP 0 .

(46a) (46b)

(47a) (47b)

For small δ t, the sign of δ Ek follows that of the first non-zero term in its expansion above. A binary is stable if this term is negative for all admissible structural and orbital perturbations Do , Do , W o , W o , R˙ 0 and R˙ 0 . S

P

S

P

S

P

Now, the leading order estimate δ (1) Ek is identically zero. We show this in a manner similar to S2015. Adding the first-order contributions, assuming that the stress post-perturbation (σ 0 ) equals that at equilibrium (σ E ), and invoking (18) and (29), we obtain

 ˙ 0 : I P · LT + I S · LT + μ2 m RE  R˙ o δ (1) Ek = −Ω oP oS P P   2  ˙ + μ m RES  RoS .

Utilizing the tensor identity in footnote 4 of S2015, (13), its derivative, (1), and (40), the above becomes

    δ (1) Ek = −ω˙ 0 · − 2 ax sk I S · LToS + I P · LToP    + m rP × r˙ oP + m rS × r˙ oS

˙ 0 · Hrel |t=0+ , = −ω which vanishes, because Hrel remains zero in frame O. The foregoing calculation showing that δ (1) Ek vanishes relied on the equilibrium stress remaining unchanged due to the perturbation (38). As discussed in detail by Sharma (2012, 2013) in the context of an isolated body, this happens only if the body is critically equilibrated9 and is perturbed by compatible structural perturbations. Compatible strain rates were defined in Section 3, and ensure that the stress post-perturbation σ 0 equals the equilibrium stress σ E , so that there is no stress jump. Sharma (2013) proved that for incompatible structural perturbations, the maximum dissipation postulate ensures that the first-order change in structural energy is negative. Briefly, this postulate reflects the impossibility of obtaining useful work from plastic materials by cyclically applying and removing external forces, i.e., plastic materials dissipate, or at best conserve energy during closed cycles; see also Lubliner (1990, p. 117). Later, Sharma (2014) employed this fact to prove stability of satellites to all incompatible structural perturbations. In order make arguments, similar to the ones outlined in the preceding paragraph, for a granular binary system, we need to define critical equilibrium and compatible structural perturbations 9 A body is critically equilibrated if its volume-averaged equilibrium stress state is on the yield surface.

132

I. Sharma / Icarus 277 (2016) 125–140

for a rubble-pile binary. A binary is said be critically equilibrated if it lies on the boundary of its equilibrium zone corresponding to its internal friction angle; see Fig. 4. In contrast to the case of asteroids and satellites, a binary may be critically equilibrated even while the stress state of one of its member is below yield. In fact, as the discussion accompanying Fig. 4 shows, at the boundary of a binary’s equilibrium zone, typically only one of the two members will be critically equilibrated; both members yield simultaneously only at a point. The stability of a binary is tested by imposing affine perturbations of the kind (38) that affect both the orbit and the structure of the primary and the secondary. The structural perturbations are a superposition of deformation and rigid-body modes. A structural perturbation of the binary will be said to be compatible if, either (a) the binary is critically equilibrated and the critically equilibrated member is perturbed by a velocity gradient with a compatible strain rate – in the sense of discussion following (23) and the previous paragraph – while the member whose stress state is below yield is subjected to a rigidbody motion, or (b) both members are below yield and their structures are imparted only rigid-body perturbations; no strain rate perturbations are imposed. An incompatible structural perturbation of the binary indicates that at least one of its members is perturbed structurally in an incompatible manner. This could mean that we suddenly expose a member below yield to strain rate perturbations, or that a critically equilibrated member is imparted an incompatible strain rate. As discussed above, this will cause a negative increment in the structural kinetic energy of that member at first order, i.e., δ (1 ) Es or δ (1 ) Es will be negative. This, in turn, will translate into a negative change at first order in the total relative kinetic energy of the binary, i.e., δ (1) Ek < 0, thereby proving the stability of the binary to incompatible perturbations. The above proves that a granular binary is stable at first order to those velocity perturbations (38) that are incompatible. It now remains now to investigate the stability of the binary to admissible compatible perturbations. This is decided at the second order, i.e., by δ (2) Ek ’s sign. We now apply the stability test to rubble-pile binaries with prolate members. 6. Example: rubble-pile binaries We model binaries as self-gravitating rubble-pile prolate ellipsoids orbiting each other on circular paths, with Fig. 3 depicting the situation at equilibrium. The rubble-pile follows the rheology of Section 3. As mentioned, a binary is stable only if it is structurally and orbitally stable to structural and orbital perturbations of both the primary and the secondary. This is now analyzed through the energy criterion. For this, as we saw above, it remains to consider the sign of the second-order estimate δ (2) Ek of the change in the relative kinetic energy Ek . We now obtain necessary formulae. We employ the subscript ‘0’ – indicating evaluation post-perturbation at t = 0+ – for only those quantities that change, or can potentially change, abruptly due to a perturbation, e.g., stresses, O’s rotation rate, etc. It is assumed that equilibrium values are employed for quantities without a subscript or with the subscript ‘E’. These clearly include quantities that depend only on a material point’s location or density, and cannot change suddenly due to velocity perturbation, e.g., mass, volume or inertia. We first summarize appropriate non-dimensionalization, which follow S2015.

stresses at equilibrium in the primary and the secondary are found from (29) as, respectively,

 −2/3 1 2 ˜ S · Q P α  β  Ω + A + B P E η ηγ12   −2/3 1 2 ˜ P · Q S α  β  and σ ES = − ΩE + AS + B , η σ EP = −



Time is scaled by 1/ 2π ρ  G, stress by (3/20π ) (2π ρ   Gm  )(4π /3V  )1/3 , energy changes such as δ (2) Ek by (2π ρ  G )2 a1 2 m /5, and inertias I P and I S by, respectively, m a12 /5 and m a1 2 /5. Employing these, the non-dimensional average

(48a) (48b)

where η = ρ  /ρ  , γ1 = a1 /a1 , and Q P and Q S are non-dimensional tensors obtained from I P and I S whose forms are given by Eq. (51) in S2015. We will also employ the scalings RS

qS =

a1

RP

, qP =  , q = |qP |, a1

and

q = |qS |,

(49)

and introduce the scaled separation

q∗ =

1

1 +  q q

−1

R =   1; a1 + a1

(50)

the inequality prevents interpenetration. Finally, from (27) and (28), the non-dimensional equilibrium rotation rate ωE is either

ωE2 =

1

μ

eˆ P · bS · eˆ P

or

1

ωE2 =

μ

eˆ S · bP · eˆ S ,

(51)

where employing (17) and the scalings introduced above, we compute

 1  (1) BP + eˆ S · BP(2) : Q S  eˆ S  2 5q  1  (0 ) and bS = BS + 2 BS(1) + eˆ P · BS(2) : Q P  eˆ P . 5q bP = BP(0) +

For brevity, we will retain the same symbols for non-dimensional quantities. 6.2. Orbital kinetic energy The second-order estimate of the orbital kinetic energy follows S2015 with minor changes. For example, δ (2 ) Eo is found by expanding (46b), and replacing for FS and its derivative from (16) and R¨ 0S from (14b), to find

  2 δ (2) Eo = −m 2π ρ  GbP + μ Ω2E : R˙ oS  R˙ oS − μ mΩ˙ 0 : RES  RES    − m 2π ρ  G b˙ P − bP trDoP   ˙ 0 : RE  R˙ o , ¨ 0+Ω ˙ 0 · Ω E + 3Ω E · Ω + μ Ω S S

where definition (15) and equilibrium condition (18) have been utilized. The difference above from the corresponding expression in Section 5.2 of S2015 is that the density ρ  may also change and this contributes the term −2π ρ  GbP trDoP . A similar formula is obtained for δ (2 ) Eo . Following Section 6.1, we non-dimensionalize δ (2 ) Eo and δ (2 ) Eo to obtain:

δ (2) Eo = −

5 μ

 1

γ 1

η

2 1

+

bP + μ ΩE

2

b˙ P − bP trDoP



2

˙ : qE  qE : q˙ oS  q˙ oS + μ Ω 0 S S



η    ˙ 0 ¨ 0+Ω ˙ 0 · Ω E + 3Ω E · Ω + μ Ω : qES  q˙ oS ,

and

δ (2) Eo = −

6.1. Non-dimensionalization



1

5μ

+



χ 

bS + μ ΩE

b˙ S − bS trDoS

+ μ



2





(52a)

2

˙ : qE  qE : q˙ oP  q˙ oP + μ Ω 0 P P

¨ 0+Ω ˙ 0 · Ω E + 3Ω E · Ω ˙ 0 Ω





: qEP  q˙ oP .

The tensor b˙ required above is computed in App. B of S2015.

(52b)

I. Sharma / Icarus 277 (2016) 125–140

6.3. Structural kinetic energy We now compute the second-order corrections δ (2 ) Es and δ (2) Es to the structural energy terms given by, respectively, (44) and (45). We limit ourselves to a rubble-pile binary with prolate members. Consider δ (2 ) Es given by (44). We expand (44) and replace for ˙LS from (5a), M P and M ˙ P from (11) and its derivative, derivative of the post-yield stress σ˙ S after differentiating (23) and invoking (25), and, finally, ρ˙  , ρ˙  and V˙  by adapting (26). Simplifying further, we obtain

δ

(2 ) 

 

κ ( k2 + 3 2 |DoS | − (1 +  )trDoS )1 − σ ES trDoS : DoS   T 2T + V  σ ES · I −1 S : Lo S · I S · Lo S + I S · Lo S    ˙  ˜ P trDoP ˜P − B − 2π ρ  G A˙ 0S − AS trDoS + 2π ρ  G B 2 ¨ 0+Ω ˙ 0 · Ω E + 3Ω E · Ω ˙ 0 · I S : LT − Ω ˙ : IS , +Ω oS 0

Es = V 

where σ ES is given by the second of (29). The above expression is similar to the one obtained in Section 5.3 of S2015, except for the presence of terms related to internal gravity and contributions due to strain rates in the primary and the secondary. In arriving at the form above, we have restricted ourselves to compatible velocity perturbations that obey (23). A similar formula may be found for the structural kinetic energy δ (2 ) Es of the primary. The non-dimensionalized δ (2 ) Es and δ (2 ) Es are:

δ (2) Es =

 η   2 κ k + 3 2 |DoP | − (1 +  )trDoP 1 − σ EP trDoP : χ  3 / 2 DoP α  β      3/2 η T 2T + σ EP · Q −1 α β P : Lo P · Q P · Lo P + Q P · L o P χ 1  ˙ 1 ˙ 2 1 ˜ S trDoS ˜ 0S − B − Ω0 : Q P − A˙ 0P − AP trDoP + B 2 2 χγ1 χ γ1 η  ˙ 0+Ω ¨ 0 + 3Ω E · Ω ˙ 0 · Ω E · Q : LT +Ω (53a) P oP

and

 

 δ (2) Es = κ k2 + 3 2 |DoS | − (1 +  )trDoS 1 − σ ES trDoS :  3 / 2 DoS α  β      3/2 T 2T + σ ES · Q −1 α β S : Lo S · Q S · Lo S + Q S · Lo S    2 ˙ : Q − A˙ 0 − AS trDo + 1 B ˜ P trDoP ˜˙ 0P − B −Ω S 0 S S η  ˙ 0+Ω ¨ 0 + 3Ω E · Ω ˙ 0 · Ω E · Q : LT , +Ω S oS

133

tem P of the primary. Section 5.4 of S2015 lists the forms of the various tensors in O. Additionally, we have



[σ E P ] =

σ 1P

0

σ 2P

0 0

0

0 0



σ 3P

and [σ ES ] =

 σ 1S 0 0

0

σ 2S 0

0 0

σ 3S



,

(54)

where σiP = σi and σiS = σi , because P, S and O are aligned; P

S

see Fig. 3. The components of the stress tensors will be found in Section 6.5. The next section will show that the compatible strain rates will be diagonal, so that



[D o P ] =

D1P 0 0

0 D2P 0

0 0 D3P





and [DoS ] =

D1S 0 0

0 D2S 0



0 0 , (55) D3S

where, again, because P, S and O are aligned, DiP = Di and DiS = P

Di . The spin rates W oP and W oS will correspond to full threeS

dimensional rotations and their forms are given by eq. (59) of S2015. 6.5. Equilibrium The previous section notes that the components of the equilibrium stress in a binary member in its principal axes coordinate system, and in O, coincide. The same holds for tensors associated with self-gravity and tidal moments. Thus,

Bi jS = BiS δi j (no sum), and Bi jP = Bi jP δi j (no sum), where we have recalled that BS and BP are diagonal in, respectively, S and P. Similarly, the components in O of BS(0 ) and BS(1 ) ,

1) i.e., Bi(j0 ) and Bi(jk , respectively, and also of BP(0 ) and BP(1 ) in O, will S

S

be the same as their corresponding components in S. These components for a prolate primary, i.e., a1  a2 = a3 , are tabulated in App. A of S2015. Adopting the scaling of Section 6.1, the non-dimensional components in O of the equilibrium stresses in the secondary and the primary may be found from (29) and definition (12). They are, respectively,

  −4/3 1  (0 ) (1 ) σ 1S = ωE2 − A1S − B11 + B β 111 P P η  1 (0) −4/3 β σ 2S = β 2 ωE2 − A3S − B11 η P  1 (0 ) β −4/3 , and σ 3S = −β 2 A3S + B33 η P and (53b)

  (0 )  −4/3 1 (1 ) 2 ω − A − B + B β 1 11S 111S ηγ12 E η P β 2  2 1 (0 ) ω β −4/3 = − A − B 3 E P 11S η ηγ12 β 2  1 (0 ) β −4/3 , =− A3P + B33 2 S ηγ1 η

σ 1P =

(56)

1

where κ is now the dimensionless post-yield bulk modulus. The tensor A is available in Sharma et al. (2009), and its derivative in ˜ and its derivatives are provided in Apps. A and B Sharma (2012), B of S2015, respectively. In estimating the second-order corrections above, we require the equilibrium stress state of the binary’s members in the frame O and the structural perturbations compatible to these stress states. These we find in Section 6.5. First we evaluate various tensors in the frame O.

and

6.4. Components

6.6. Perturbations

Equilibrium and stability computations will be performed in the rotating coordinate system O defined in Section 5. We will require to evaluate foregoing formulae in O. At equilibrium, and just after perturbation, O is aligned with the principal axes coordinate sys-

Structural perturbation of a member are compatible with the corresponding stress state only if their associated strain rates satisfy the flow rule (23). Employing components from Section 6.4, we find that the strain rate perturbations in (38) are compatible only

σ 2P σ 3P

(57)

where ωE is given in (51). The shear stresses are zero in both members. The equilibrium pressures pE and pE in, respectively, the secondary and the primary may now be obtained from (20).

134

I. Sharma / Icarus 277 (2016) 125–140

when the strain-rate perturbations DoS and DoP are of the form (55), and their diagonal components satisfy

D2S = D1S and

σ 2S + pE (1 +  ) D3S σ 3S + pE (1 +  ) and = ,  D1S σ 1S + pE (1 +  ) σ 1S + pE (1 +  ) D2P = D1P

(58a)

σ 2P + pE (1 +  ) D3P σ 3P + pE (1 +  ) and = .  D1P σ 1P + pE (1 +  ) σ 1P + pE (1 +  ) (58b)

Thus, the strain rate tensors DoP and DoS are determined up to an arbitrary constant multiplier by their associated equilibrium stresses σ EP and σ ES . As discussed, the non-dimensional perturbation rates q˙ oS , q˙ oP , DoS , DoP , W oS and W oP cannot be chosen independently. For compatible perturbations, DoS and DoP are diagonal in O. Employing this fact, solving (39) and (40), and non-dimensionalizing, yields

χ χ β 2 q˙ , q˙ 3P = − q˙ 3S , W1P = −γ12 2 W1S γ1 2S γ1 β     1 q˙ 2S = − 1 + β 2 W3P + χ γ12 1 + β 2 W3S ,   5μ q    1  and q˙ 3S = 1 + β 2 W2P + χ γ12 1 + β 2 W2S ,   5μ q q˙ 2P = −

(59a) (59b) (59c)

which are the same as the relations found for rigid binaries in Section 5.5 of S2015. Furthermore, DoS and DoP are known

up to an arbitrary positive scaling of the magnitudes DoS and DoP through, respectively, σ ES and σ EP . Finally, as in S2015, a prolate binary’s stability will be unaffected by W1S and W1P . Thus, when probing stability, we may independently select only q˙ 1S , W2S , W3S , W2P , W3P and the magnitudes of the strain rate tensors DoS and DoP . It is convenient to visualize a general perturbation as a vector in a seven dimensional by the ‘unit vectors’

space spanned

q˙ 1S , W2S , W3S , W2P , W3P , DoP and DoS . Thus, we may express any admissible perturbation of the binary as

c1 q˙ 1S + c2W2S + c3W3S + c4W2P + c5W3P + c6 |DoP | + c7 |DoS |,

(60)

with the added restriction that c6 and c7 not be negative. It is clear that δ (2) Ek ’s sign depends only on the relative ratio of the various unit vectors q˙ 1S ,..., etc. Thus, we will assume that the ‘energy’ of  the perturbation is fixed, i.e., we impose 7i=1 ci2 = 1. All admissible perturbations may now be explored by varying ci appropriately. We finally turn to the stability analysis of rubble-pile binaries. 6.7. Local stability We investigate the local stability of a binary through the energy criterion. We emphasize that our stability analysis accounts for structural and orbital stability of both binary members, considering them as frictional rubble piles of comparable mass. Past stability analyses of granular satellites (Sharma, 2014) or asteroids (Sharma, 2013) may be recovered as special cases when the mass ratio χ → 0. Indeed, at very small χ , the primary’s orbital motion is negligible, and its structural dynamics is unaffected by the secondary, i.e., δ Eo  δ Es , and Es depends only on the primary’s structural deformation, analogous to the case of an isolated asteroid. Similarly, for χ  1, the secondary behaves as a rubble satellite of the (relatively) massive primary, and its stability analysis will parallel the one developed in Sharma (2014). For simplicity, and also because little is known about densities of known binaries, we henceforth take the density ratio η = 1. S2015 investigated the stability of prolate binaries with rigid members. It was seen that a rigid binary was stable to W3S and W2S ; these were labeled there as, respectively, in-plane and out-of-plane

Fig. 5. The stability of a granular binary with equal-sized (χ = 1) prolate members separated by q∗ = 1.5 is explored in β  − β  space. Stability regions are shaded and identified by numbers corresponding to non-dimensional plastic bulk moduli κ . Boundaries of equilibrium zones for various internal friction angles φ F are shown by dashed lines that follow Fig. 5 of S2015 I. The binary ‘B’ is an example discussed in the text.

misorientations. Furthermore, a rigid binary was unstable to radial perturbations q˙ 1S whenever the binary’s members were closer than a critical distance q∗c (β  , β  ), which depended on the shapes of the primary and the secondary. It is clear from our earlier discussion that rigid-body modes – as defined by q˙ 1S , W2S and W3S – form a subset of affine perturbations (38). Investigating the stability of a rigid binary is, therefore, equivalent to considering the stability of a granular binary to rigid-body modes alone. Consequently, the current stability analysis recovers the results of S2015 as a special case in the absence of strain-rate perturbations. Thus, we need only investigate those rubble binaries that are stable to rigid body perturbations. These will necessary be binaries whose members are separated by more than q∗c . Consider, first, binaries whose members are not critically equilibrated, i.e., their average stress states are below yield, and whose separation is greater than q∗c . These are, therefore, stable to rigidbody modes that are, also, the only compatible perturbations. We have already seen, in the discussion following (46) and (47), that such binaries will be locally stable at first order to any perturbation that deforms their structure. Thus, such binaries are always stable. We turn then to critically equilibrated granular binaries, again more than q∗c apart. We need only to test stability to compatible affine perturbations. We will explore this on appropriate twodimensional sections of the binary’s equilibrium landscape. For each critical equilibrium configuration, which entails a choice of β  , β   and q∗ , and each choice of plastic bulk modulus κ , we will explore all possible combinations of orbital and structural perturbations of the primary and the secondary, as defined by (60). Figs. 5–7 plot the stability zones for binaries with prolate member for mass ratios χ = 0.5, 0.7 and 1, and separations q∗ = 1.5 and 1.4. At these values S2015 shows that the binary system is stable to rigid-body modes. In each case, several plastic moduli κ are considered. For each κ , there is associated a shaded region within which all critically equilibrated binaries are secularly stable following affine perturbations, and, hence, locally stable. Equivalently, binaries that lie on their corresponding equilibrium zone’s boundary outside the shaded stability region associated with their plastic bulk modulus will be secularly, and so

I. Sharma / Icarus 277 (2016) 125–140

135

Fig. 6. Continued from Fig. 5. Stability regions of prolate granular binaries displayed in β  − β  space for two different mass ratios χ at q∗ = 1.5. See also Fig. 5’s caption.

locally, unstable. In all cases it is observed that the stability zone is most significantly determined by structural perturbations, while orbital – specifically radial (q˙ 1S ) – perturbations lead to relatively minor reshaping of the stability zones. We also see that the stability of a granular binary is unaffected by the presence/absence of in-plane (W2S ) out-of-plane (W3S ) misorientations. Finally, following the discussion on the structure of the equilibrium landscape in Fig. 4, we may similarly identify the boundaries of the stability region. In general lower/horizontally aligned stability boundaries are obtained by testing the correspondingly oriented critical equilibrium curves. These in turn correspond to secondary, so that we may say that the lower/horizontally aligned stability boundaries are dictated by the secondary. Analogously, the upper/vertically directed stability boundaries relate to a critically equilibrated primary. As an example, consider a equal-sized (χ = 1) rubble binary B with prolate members having β  = 0.7 and β  = 0.4, and separated by q∗ = 1.5. Fig. 9a of S2015 confirms that q∗c (0.75, 0.4 ) ≈ 1.19, so that B is stable to rigid-body modes. The binary may be

Fig. 7. Stability regions of prolate granular binaries separated by q∗ = 1.4 displayed in β  − β  space for three different mass ratios χ . See also Fig. 5’s caption.

136

I. Sharma / Icarus 277 (2016) 125–140

located in Fig. 5 as shown. Let us assume that B’s internal friction angle φF = 20o. Its location in Fig. 5 indicates that the binary’s members are not critically equilibrated. Thus, the discussion above shows that B will be locally stable to affine velocity perturbations. However, in case the secondary’s shape was more slender, with β   ≈ 0.3, then B would lie on the boundary of the φF = 20◦ equilibrium region as shown by the open circle in Fig. 5. The binary would then be critically equilibrated. However, this location lies outside the stability regions corresponding to any choice of plastic bulk modulus κ . Our analysis will then find the binary B to be locally unstable, and would raise doubts about B being granular. Alternatively, it could be that B’s internal friction is higher than 20°, so that even when β  = 0.3 it is not critically equilibrated, and, therefore, locally stable at first order itself. It is interesting to note that for χ = 1, except for binaries both of whose members are very slender, the equilibrium curves corresponding to friction angles φ F ≥ 25° typically lies outside any of the shaded local stability regions. This indicates that equal-mass rubble-pile binaries separated by q∗ = 1.5 and with a φ F ≥ 30° cannot typically be critically equilibrated, no matter what is their plastic bulk moduli. Fig. 6 repeats Fig. 5, but this time for members with unequal masses. The essential structure of the stability regions in Fig. 5 is preserved, although we see that stability zones reduce with χ , i.e., the bigger the primary is relative to the secondary, the smaller the stability regions. We also see that stability is lost more due to the structural instability of the secondary than the primary. This may be seen by the skewed nature of stability regions about the diagonal through (0.1, 0.1). In Fig. 7, we report the stability of a prolate granular binary with three different mass ratios, whose members are separated by q∗ = 1.4. When compared to the stability plots of Figs. 5 and 6, we observe that the stability regions reduce dramatically, even though the separation between the members has only decreased by about a tenth of the primary’s size. In fact, our analysis reveals that critically equilibrated binaries that are separated by q∗  1.4 are not secularly stable. A similar result holds for binaries with mass ratios χ  0.5 and separated by q∗  1.5. At the same time, recall that any binary that is close enough – q∗c ≈ 1.1 for χ = 1 – is unstable to rigid body perturbations, whether or not it be critically equilibrated. Thus, granular binaries between q∗ = q∗c and q∗ about 1.4 cannot be critically equilibrated. The above highlights the crucial role played by the separation q∗ and mass ratio χ on the stability of granular binaries. For each case considered above, we observe that as the plastic bulk modulus is lowered, the stability zones expand, with stability regions associated with lower κ containing the one corresponding to higher κ as a subset. This was also observed in the case of asteroids and satellites by Sharma (2013, 2014), and the explanations provided there in may be adapted easily for the present case. 6.8. Stability to finite structural perturbations The local analysis pursued above was limited to infinitesimal perturbations of critically equilibrated binaries. However, given the improbability of finding critically equilibrated binaries, it is of interest to investigate binaries both of whose members have equilibrium stress states lying within the yield surface. The previous section’s local analysis found such a binary to be locally stable whenever they are separated by more than q∗c ; at lower separations the binary is unstable to rigid-body modes. Thus, for q∗ > q∗c , the maximum dissipation postulate ensured stability at first order in the presence of structural perturbation of either member. This was because the leading order estimate of the change in the binary’s structural kinetic energy Eks , i.e., δ (1 ) Eks = δ (1 ) Es + δ (1 ) Es , and so also δ (1) Ek , was negative. However, if structural perturbations are finite, then it may happen that Eks increases following an initial decline, and this may lead to instability. We make a first at-

tempt to analyze this following Sharma (2013, 2014). The primary and the secondary are taken more than q∗c away. Consider a binary the equilibrium stresses of both of whose members are within yield, so that the binary’s location in β  − β  space is within the boundary of its associated equilibrium region. For such binaries, δ (1 ) Eks was shown to be negative following any perturbation that deformed the structure of either the members. However, if structural perturbations are large enough to ensure that motion persists for times longer than δt ∗ ∼

the subsequent

2 δ (1 ) Eks /δ (2 ) Eks , where

δ (2) Eks = δ (2) Es + δ (2) Es ,

(61)

then the contribution δ (2 ) Eks δt 2 /2 of δ (2 ) Eks would overcome that of

δ (1) Eks . Consequently, if δ (2) Eks > 0, it may lead to δ Eks , and consequently δ Ek , becoming positive. In that case, a locally stable binary may become unstable at longer times. Here, for definiteness, we will assume that binaries lying within their associated equilibrium zones will be unstable if δ (2) Ek > 0. This assumption will be better the greater δ (2 ) Eks is in comparison to δ (1 ) Eks , as will generally hold for binaries that are close to being critically equilibrated. The expression for δ (2 ) Eks for binaries that are not critically equilibrated will involve finding δ (2 ) Es and δ (2 ) Es . These may be found in a manner similar to the derivation of (53). We find

 

δ (2) Es = κ k2 + 3 2 |DoS | − (1 +  )trDoS 1 − σ ES trDoS   3 / 2  + 2σ ES · DoS : DoS α  β  − A˙ 0S − AS trDoS  1˙ ˙ 0+Ω ˙ 0 · Ω E · Q : Do ˜ P trDoP + 3ΩE · Ω ˜ −B + B S S η 0P   2 ˙ : Q , +σ S : D2 + Q −1 · σ S (α  β  )2/3 −Ω S 0 oS S × (α  β  )2/3 , and

δ (2) Es =



1

χ γ12

(62a)

 

 ηγ12 κ k2 + 3 2 |DoP | − (1 +  )trDoP 1

− σ EP trDoP + 2σ EP · DoP : DoP −

1 η





α  β 

2 / 3

˙ 0 ˜ S trDoS + 3ΩE · Ω ˜˙ 0S − B A˙ 0P − AP trDoP + B



2

˙ 0 · ΩE · Q : Do − Ω ˙ :Q +Ω P P 0 P







  2/3 (α  β  )2/3 , + ηγ12 σ P : D2oP +ηγ12 Q −1 P · σ P ( α β )

(62b) where the jumps in the stress states of the secondary and the primary from their equilibrium to their post-perturbation values, are, respectively,

σ S = σ 0S − σ ES = sES +



k pE − |sES | sES , |sES |

and

σ P = σ 0P − σ EP = sEP +



k pE − |sEP | sEP . |sE |

The above formulae are obtained by assuming that perturbations do not change the equilibrium pressure. The equilibrium deviatoric stress sE simply scales to a value s0 that lies on the yield surface, and where the strain rate Do is compatible with σ 0 . This choice guarantees the smallest possible values for δ (1 ) Eks and δ (1 ) Eks . Note that sES and sEP may be obtained by combining (21) with, respectively, (57) and (56), whenever the binary’s shape and rotation rate ωE are known. In the above, the states ‘E’ and ‘0’ correspond, now, to locations A and B in Fig. 2, respectively.

I. Sharma / Icarus 277 (2016) 125–140

137

Similar analyses may be pursued to investigate the effect of even higher-order corrections to δ Ek in the presence of finite structural perturbations. 7. Application: near-Earth binaries

Fig. 8. Regions in in β  − β  space where δ (2) Ek < 0 for prolate binaries with q∗ = 1.5, χ = 1 and φF = 30o . The shaded regions correspond to different nondimensional plastic bulk moduli κ . See also Fig. 5’s caption.

S2010 explored the equilibria of the binary systems 216 Kleopatra, 25143 Itokawa, 624 Hektor and 90 Antiope. The stability of these systems considering their members to be rigid objects was analyzed by S2015. We will now utilize the framework developed above to investigate the stability of these binaries as granular aggregates. We begin by first discussing in the next section appropriate values for the material parameters employed to model rubble piles. We note that we will apply our results to binaries whose primary and secondary are of comparable mass, i.e., the mass ratio χ is O(1). Systems with smaller χ , such as 1999 KW4 that has χ ≈ 0.04, are much more easily analyzed as follows: (a) the primary is investigated as an isolated rubble pile, as per the equilibrium and stability analyses of Sharma et al. (2009) and Sharma (2013), respectively, and (b) the secondary is modeled as a rubble satellite and probed through an appropriate extension of the methods of Sharma (2009, 2014). 7.1. Material parameters

Fig. 9. Location of 216 Kleopatra in β  –β   space assuming it to be an equal-mass (χ = 1) prolate binary separated by q∗ = 1.5. The two shape models discussed in the text are shown as I and II. Stability zones for various κ are shaded. Critical equilibrium curves for various friction angles φ F are displayed as dashed lines; cf. Figs. 4 and 5.

Fig. 8 displays regions in β  − β  space where δ (2) Ek , with

δ (2) Eks given by (61) and (62), is negative for an ellipsoidal equal-

mass rubble binary separated by q∗ = 1.5. We take the internal friction angle φ F of the binary to be 30°. The regions are parameterized by the non-dimensional plastic bulk modulus κ . Because we have taken the members of the binary to have φF = 30◦ , these regions do not extend beyond the φF = 30◦ equilibrium zone. We observe that there are no regions for κ  1 where δ (2) Ek < 0. For κ = 0 and 1, only a small corner of the φF = 30◦ region has δ (2) Ek < 0. This is a generic feature for other choices of χ , q∗ and φ F , with the regions getting smaller at lower q∗ and greater φ F . This indicates that while most binaries with internal friction angle φF = 30◦ that are close to being critically equilibrated may be stable at first order, as per the maximum dissipation postulate, they will be unstable to finite perturbations.

In Section 3 we modeled granular aggregates like soils through a two-parameter model in terms of its internal friction angle φ F and its post-yield plastic bulk modulus κ . To apply our stability results to granular binaries we require estimates of both φ F and κ . Both these parameters depend significantly on the pressure confining these binaries. This pressure in typical binaries is on the order of a few tens of Pascals, making rubble binaries extremely lightly-confined granular aggregates. Sharma (2013) discusses in detail tests by which these parameters are measured for terrestrial soils, and highlights the paucity of data for lightly-confined granular aggregates. For granular binaries, we employ the parameter choices made by Sharma (2013) for asteroids. Thus, we will assume that a rubble binary may have an internal friction angle φ F over the broad range of 20◦ − 50◦ , and allow κ to be of the order of a kPa. This would correspond to a non-dimensional of κ = 10 for a body like Itokawa. We refer the interested reader to Sharma (2013) for further details. Finally, we note that Scheeres et al. (2010) suggests that in granular bodies with meter-sized or smaller constituent grains, cohesion due to van der Waals forces may have an effect comparable to internal and external gravitational interactions. In contrast to our rheological model, that is suitable for cohesionless grains, cohesive aggregates are able to accommodate tensile stresses to some extent. It is possible to extend the rigid-plastic material model utilized here to include cohesive effects; see, e.g., Chen and Han (1988, Section 2.3.4, p. 94). However, we do not pursue this here. 7.2. Local stability Scheeres (2009) and S2015 have shown that contact10 rigid binaries are unstable. Thus, Itokawa will require a force-carrying ‘neck’ in order to survive. The present analysis provides the additional information that Itokawa is also locally unstable to deformation modes. Indeed, the equilibrium analysis of (Sharma, 2010, Fig. 18) found Itokawa to be at the φF = 50◦ boundary, which placed it tantalizingly close to being critically equilibrated; see 10 A contact binary was taken to be a special case of an orbital binary with zero separation. The primary and the secondary were assumed to be resting upon each other, and the contact force, if any, was ignored

138

I. Sharma / Icarus 277 (2016) 125–140

Section 7.1. The analysis of Section 6.7 now shows that critically equilibrated binaries having members separated by q∗  1.4 are secularly unstable to strain-rate perturbations alone, irrespective of the plastic bulk modulus κ . Hektor’s stability was probed by S2015, considering it to be comprised of prolate, rigid members. It was seen that the most probable scenario is for Hektor to be a contact binary connected by a force-carrying ‘bridge’. Other combinations of mass ratio, separation and shape that would allow a separated Hektor to exist as a stable binary could not be found. Kleopatra has recently been studied by Hirabayashi and Scheeres (2013) and Descamps (2015). Hirabayashi and Scheeres (2013) constrained Kleopatra’s shape by considering surface shedding and structural failure through volume-averaging. They performed linear elastic analysis and concluded that structural failure was more important than surface shedding. As per their linear elastic analysis, the possibility that Kleopatra’s ‘bridge’ could fail plastically (global collapse) arises at the largest admissible size scale; see Hirabayashi and Scheeres (2013, Figs. 9–11). Now, when modeling granular aggregates as linear elastic materials it is important to note the following: (a) In granular aggregates, local plastic yielding precedes global collapse, and these cause local slipping and readjustments. The consequent shape changes are often large enough to not be captured by a linear analysis. (b) The continuum response of granular materials is generally not linear elastic. Firmer conclusions will, thus, require either a large deformation, elastoplastic analysis, or, preferably, a DE simulation as in, e.g., Sanchez and Scheeres (2012). Given these caveats, we believe that there is merit in exploring the Kleopatra’s ‘bridge’ through an alternative approach. Our stability analysis below shows that, contrary to the conclusions of Hirabayashi and Scheeres (2013), Kleopatra’s ‘bridge’ may not be near failure. Descamps (2015) investigated Kleopatra by modeling it as a fluid contact binary and, so, obtained estimates of its physical characteristics. Both Hirabayashi and Scheeres (2013) and Descamps (2015) analyses do not, however, consider the stability of Kleopatra as a frictional rubble pile. Earlier, S2010 and S2015 explored two self-similar shape models of Kleopatra that had either β  = β  = 0.5 or β  = β  = 0.4 due to, respectively, Tanga et al. (2001) and Hestroffer et al. (2002). S2015 found that a minimum scaled separation q∗ of about 1.2 is required for the survival of Kelopatra as a rigid binary, unless its primary and secondary too were linked by a force-carrying ‘bridge’. Recently, Descamps et al. (2011) reported that Kleopatra does appear to have a bridge around 90 km long connecting its two lobes. As a first approximation,11 we explore the stability of Kleopatra assuming that the ‘bridge’ does not contribute to the structural and orbital dynamics of Kleopatra, i.e., we ignore the inertia of the ‘bridge’ and any possible stress that it may carry. In doing so, we arrive at a model for Kleopatra that consists of two members separated by q∗ ≥ 1.5. A benefit of such an analysis will be to reveal the dynamical significance of Kleopatra’s ‘bridge’. For example, if our stability analysis predicts that a separated Kleopatra is unstable, then the ‘bridge’ is crucial to Kleopatra’s survival. We now explore Kleopatra’s stability assuming, as in S2010, that its lobes are equal-mass prolate ellipsoids, each with dimensions 72 × 36 km. The separation is taken to be about 70 km, so that q∗ = 1.5. Fig. 9 locates both models of Kleopatra in Fig. 7a. We see that a thinner Kleopatra with β  = β  = 0.4, and an internal friction close to 20°, could be critically equilibrated. In that case, Kleopatra would require a low non-dimensional plastic modulus of about 1 – equivalent to 0.5 kPa – to be stable. Larger plastic

11 The degree to which this approximation is justified will only be known after the availability of detailed simulations that estimate the stress carried by such ‘bridges’.

Fig. 10. Stability region in β  –β   space for an equal-mass (χ = 1) binary with members separated by q∗ = 2. Stability zones for various κ are shaded. Critical equilibrium curves for various friction angles φ F are displayed as dashed lines; cf. Figs. 4 and 5. This plot follows Fig. 20 of S2010.

bulk moduli may be admitted for more widely separated Kleopatra, as stability regions enlarge with separation; cf. Section 6.7. In case Kleopatra’s internal friction was higher, say 30°, then, because it will be located within the corresponding equilibrium regions, it will be stable to infinitesimal perturbations. The conclusion then is that Kleopatra’s ‘bridge’ is not crucial to its local stability, so that the possibility of the ‘bridge’ being under high stresses – and hence close to failure – is low. This prediction is further reinforced by the analysis of the next section. Finally, 90 Antiope was studied through the “Roche formalism” by Descamps (2008). Antiope was modeled in S2010 as a samesized binary (χ = 1) with separation q∗ = 2. At this separation Antiope is stable as a rigid binary. Fig. 10 investigates the stability regions in the context of a rubble Antiope. It is seen that a critically equilibrated Antiope with internal friction angles greater than about 20° will require a very low plastic modulus κ to be stable. On the other hand, in case the non-dimensional κ of a granular Antiope was about, say, 500, it would almost certainly not be critically equilibrated, provided that we accept that φ F  20° for granular aggregates. Of course, much more information about Antiope’s shape is necessary to draw firmer conclusions. 7.3. Stability to finite structural perturbations We now explore what the analysis of Section 6.8 may tell us about the stability to finite perturbation of some of the binaries considered above. Fig. 11 plots both proposed shape models of Kleopatra assuming two different values for Kleopatra’s internal friction φ F . From Fig. 11a we find that a Kleopatra with φF = 20◦ will be stable to finite perturbations only if its scaled post-yield bulk modulus is less than 1. The previous section’s local stability analysis too reached a similar conclusion. Thus, a Kleopatra with φ = 20◦ and κ ≤ 1 is stable irrespective of the stresses carried by its ‘bridge’. On the other hand, from Fig. 11b we observe that, while both shape models are locally stable if Kleopatra had φF = 30◦ , they are potentially unstable to finite perturbations for all plastic bulk moduli. This would indicate that Kleopatra’s ‘bridge’ will now play an important role. However, this conclusion must be balanced by the fact that both models of Kleopatra in Fig. 11b lie some distance within the φF = 30o critical curve. It is thus possible that, in this case, the perturbations required to reshape Kleopatra are large

I. Sharma / Icarus 277 (2016) 125–140

Fig. 11. Location of the two shape models of Kleopatra in regions in β  − β  where δ (2) Ek < 0; cf. Fig. 9. The separation q∗ = 1.5 and χ = 1. Two internal friction angles are explored. See also Figs. 8’s caption.

enough to require a stability analysis of even higher order than the one developed in Section 6.8. Finally, Fig. 12 shows that in case the internal friction angle of an equal-mass Antiope was about 20°, then to ensure stability to finite structural perturbations, the shape of its members would need to satisfy 0.4  β  , β    0.8. On the other hand, in case the friction angle of Antiope was closer to 30°, its members would need to be much more slender in order to be stable to finite structural perturbations. In each case, the conclusions assume that the effect of finite perturbations are well estimated by the higher-order analysis of Section 6.8. This may not be the case for binaries located well within their corresponding equilibrium regions, which will demand a still higher-order stability analysis. 8. Conclusion This paper completes the stability analysis of rubble-pile binary systems begun in S2015. A binary was defined to be stable provided the orbit and the structure of its members were stable to

139

Fig. 12. Regions in β  − β  where δ (2) Ek < 0 for binaries with q∗ = 2 and χ = 1. Two different internal friction angles φ F are explored. See also Fig. 8’s caption.

both orbital and structural perturbations. As the primary and the secondary were of comparable size, and close together, it was not possible to decouple orbit and structural stability investigations, as in the case of satellites in Sharma (2014). We, thus, had to test stability by considering change in the total kinetic energy of the binary, relative to some appropriate frame O. The total kinetic energy comprised of orbital and structural energies of both members. Also necessary was an appropriate generalization, to the case of binaries, of the coordinate system O in which to test stability through the energy criterion. The stability criterion was employed to test the stability of rubble binaries to both infinitesimal and finite perturbations. We saw that stability was sensitively dependent on the separation of its members. Critically equilibrated binaries could not be locally stable for separations q∗ ≤ 1.4. Similarly, it was observed that binaries that were close to being critically equilibrated had a very narrow region of stability, in which congruent shapes were preferred. We applied our results to four binary systems: 216 Kleopatra, 25143 Itokawa, 624 Hektor and 90 Antiope. Contact binaries like Itokawa and Hektor were known to be unstable to rigid-body modes from Scheeres (2009) and S2015, and we showed how those

140

I. Sharma / Icarus 277 (2016) 125–140

results are recovered as a special case of the more general stability analysis of rubble binaries. We also modeled Kleopatra as a separated binary, in order to explore the dynamical significance of the ‘bridge’ reported by Descamps et al. (2011). We saw that at friction angles of about 20° and low plastic bulk moduli, a separated Kleopatra is stable to infinitesimal and finite perturbations, in which case the ‘bridge’ is not crucial to Kleopatra’s survival. This, in turn, indicates a low possibility of the ‘bridge’ carrying high stresses and, so, being susceptible to failure. This is in contrast to the predictions of Hirabayashi and Scheeres (2013). On the other hand, were Kleopatra to have φ F ≥ 30°, then the ‘bridge’ could be important; however at these φ F we may require a stability analysis of even higher order than the one presented here. Finally, only very broad constraints could be placed on Antiope, and more shape information is needed to obtain firmer conclusions. Looking ahead, we see that there is a need to investigate the force-carrying ‘neck’ or ‘bridge’ that links contact binaries. In this context, various relative alignments of the two members should also be explored. Triaxial shapes for the primary and the secondary may be included, as and when more shape information becomes available. Finally, it will be interesting to explore the dynamical passage into equilibrium of such rubble-pile binaries. As in the case of rubble piles in pure spin (Sharma et al., 2009), this exercise will provide insight into the accretionary process that formed/shaped them, and will also be helpful for comparisons with n-body simulations of these granular binaries. Acknowledgments I would like to thank two anonymous referees and Dr. Michael Efroimsky of US Naval Observatory for insightful reviews that have made this work better. I am grateful to the members of the Mechanics & Applied Mathematics Group at IIT Kanpur (www.iitk.ac.in/mam) for helpful discussions. Finally, I thank Khyati Makkar, Ashish Bhateja, and Arvind of Yukti Lab at IIT Kanpur, for help with the plots. This research was supported by grant (PRL/ADMAC/PB/2013-14) under the PLANEX program of the Physical Research Laboratory, India. References Abe, S., et al., 2006. Mass and local topography measurements of Itokawa by Hayabusa. Science 312, 1344–1347. Breiter, S., Rozek, A., Vokrouhlicky, D., 2012. Stress fields and nutation damping for inelastic three-axis ellipsoids. Mon. Not. R. Astron. Soc. 165, 403–411. Burns, J.A., Safronov, V.S., 1973. Asteroid nutation angles. Mon. Not. R. Astron. Soc. 165, 403–411. Chakrabarty, J., 1969. On uniqueness and stability in rigid/plastic solids. Int. J. Mech. Sci. 11 (9), 723–731. Chandrasekhar, S., 1969. Ellipsoidal Figures of Equilibrium. Yale University Press, New Haven, CT. Chen, W.F., Han, D.J., 1988. Plasticity for Structural Engineers. Springer-Verlag, New York. Descamps, P., 2008. Roche figures of doubly synchronous asteroids. Planet. Space Sci. 56, 1839–1846. Descamps, P., 2015. Dumb-bell-shaped equilibrium figures for fiducial contact-binary asteroids and EKBOs. Icarus 245, 64–79. Descamps, P., et al., 2011. Triplicity and physical characteristics of asteroid (216) Kleopatra. Icarus 211, 1022–1033. Efroimsky, M., Lazarian, A., 20 0 0. Inelastic dissipation in wobbling asteroids and comets. MNRAS 311, 269–278. Goldreich, P., Sari, R., 2009. Tidal evolution of rubble piles. Astrophys. J. 691, 54–60. Greenwood, D.T., 1988. Principles of Dynamics. Prentice-Hall, Englewood Cliffs, NJ. Harris, A.W., Fahnestock, E.G., Pravec, P., 2009. On the shapes and spins of “rubble pile” asteroids. Icarus 199, 310–318. Hestroffer, D., Berthier, J., Descamps, P., et al., 2002. 216 Kleopatra: Tests of the radar-derived shape model. Astron. Astrophys. 392, 729–733.

Hill, R., 1957. Stability of rigid-plastic solids. J. Mech. Phys. Solids 6 (1), 1–8. Hirabayashi, M., Scheeres, D.J., 2013. Analysis of asteroid (216) Kleopatra using dynamical and structural constraints. Astrophys. J. 780, 160–170. Holsapple, K.A., 2001. Equilibrium configurations of solid cohesionless bodies. Icarus 154, 432–448. Holsapple, K.A., 2008. Spinning rods, elliptical disks and solid ellipsoidal bodies: Elastic and plastic stresses and limit spins. Int. J. Nonlin. Mech. 43, 733–742. Holsapple, K.A., 2010. On YORP-induced spin deformations of asteroids. Icarus 205, 430–442. Holsapple, K.A., Michel, P., 2006. Tidal disruptions: A continuum theory for solid bodies. Icarus 183, 331–348. Holsapple, K.A., Michel, P., 2008. Tidal disruptions II. A continuum theory for solid bodies with strength with applications to the Solar System. Icarus 193, 283–301. Jacobson, S.A., Scheeres, D.J., 2011. Long-term stable equilibria for synchronous binary asteroids. Astrophys. J. Lett. 736, L19:1–5. Knowles, J.K., 1998. Linear Vector Spaces and Cartesian Tensors. Oxford University Press, Oxford. Lai, D., Rasio, F., Shapiro, S.L., 1993. Ellipsoidal figures of equilibrium: compressible models. Astrophys. J. Suppl. S. 88, 205–252. Lai, D., Rasio, F., Shapiro, S.L., 1994. Equilibrium, stability, and orbital evolution of close binary systems. Astrophys. J. 423, 344–370. LaSalle, J., Lefschetz, S., 1961. Stability by Lyapunov’s Direct Method. Academic Press, New York. Love, A.E.H., 1946. A Treatise on the Mathematical Theory of Elasticity, fourth ed. Dover, New York. Lubliner, J., 1990. Plasticity Theory. Macmillan, New York. Maciejewski, A.J., 1995. Reduction, relative equilibria and potential in the two rigid bodies problem. Celest. Mech. Dyn. Astron. 63, 1–28. Merline, W.J., Weidenschilling, S.J., Durda, D.D., et al., 2002. Asteroids do have satellites. In: Bottke Jr., W.F., Cellino, A., Paolicchi, P., Binzel, R.P. (Eds.), Asteroids III. University of Arizona Press, Tucson, AZ, pp. 289–312. Michel, P., Benz, W., Richardson, D.C., 2001. Collisions and gravitational reaccumulation: Forming asteroid families and satellites. Science 294, 1696–1700. 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. Richardson, D.C., Bottke Jr., W.F., Love, S.G., 1998. Tidal distortion and disruption of earth-crossing asteroids. Icarus 134, 47–76. Richardson, D.C., Walsh, K.J., 2006. Binary minor planets. Annu. Rev. Earth Planetary Sci. 34, 47–81. Sanchez, P., Scheeres, D.J., 2012. DEM simulation of rotation-induced reshaping and disruption of rubble-pile asteroids. Icarus 218, 876–894. Scheeres, D.J., 2007. Rotational fission of contact binary asteroids. Icarus 189, 370–385. Scheeres, D.J., 2009. Stability of the planar full 2-body problem. Celest. Mech. Dyn. Astron. 104, 103–128. Scheeres, D.J., 2012. Minimum energy configurations in the N-body problem and the celestial mechanics of granular systems. Celest. Mech. Dyn. Astron. 113, 291–320. Scheeres, D.J., Hartzell, C.M., Sanchez, P., et al., 2010. Scaling forces to asteroid surfaces: the role of cohesion. Icarus 210, 968–984. Sharma, I., 2009. The equilibrium of rubble-pile satellites: the Darwin and Roche ellipsoids for gravitationally held granular aggregates. Icarus 200, 636–654. Sharma, I., 2010. Equilibrium shapes of rubble-pile binaries: the Darwin ellipsoids for gravitationally held granular aggregates. Icarus 205, 638–657. Sharma, I., 2012. Stability of rotating non-smooth complex fluids. J. Fluid Mech. 708, 71–99. Sharma, I., 2013. Structural stability of rubble-pile asteroids. Icarus 223, 367–382. Sharma, I., 2014. Stability of rubble-pile satellites. Icarus 229, 278–294. Sharma, I., 2015. Stability of rubble-pile binaries. Part I: Rigid binaries. Icarus 258, 438–453. Sharma, I., Burns, J.A., Hui, C.-Y., 2005. Nutational damping times in solids of revolution. Mon. Not. R. Astron. Soc. 359, 79–92. Sharma, I., Jenkins, J.T., Burns, J.A., 2006. Tidal encounters of ellipsoidal granular asteroids with planets. Icarus 183, 312–330. Sharma, I., Jenkins, J.T., Burns, J.A., 2009. Dynamical passage to approximate equilibrium shapes for spinning, gravitating rubble asteroids. Icarus 200, 304–322. Sridhar, S., Tremaine, S., 1992. Tidal disruption of viscous bodies. Icarus 95, 86–99. Storåkers, B., 1977. On uniqueness and stability under configuration-dependent loading of solids with or without a natural time. J. Mech. Phys. Solids 25 (4), 269–287. Tanga, P., Hestroffer, D., Berthier, J., et al., 2001. HST/FGS observations of the asteroid (216) Kleopatra. Icarus 153, 451–454. Walsh, K.J., Richardson, D.C., 2006. Binary near-earth asteroid formation: rubble pile model of tidal disruptions. Icarus 180, 201–216. Walsh, K.J., Richardson, D.C., Michel, P., 2008. Rotational breakup as the origin of small binary asteroids. Nature 454, 188–191.