Effects of thermal radiation on the dynamics of binary NEAs

Icarus 176 (2005) 418–431 www.elsevier.com/locate/icarus

Effects of thermal radiation on the dynamics of binary NEAs ´ a,∗,1 , Joseph A. Burns a,b Matija Cuk a Department of Astronomy, Cornell University, Ithaca, NY 14853, USA b Department of Theoretical and Applied Mechanics, Cornell University, Ithaca, NY 14853, USA

Received 15 November 2004; revised 3 February 2005 Available online 7 April 2005

Abstract The Yarkovsky force, produced when thermal radiation is re-emitted asymmetrically, causes significant orbital evolution of asteroids in the 10 m–10 km size range. When acting on a non-spherical body, the momentum carried by this radiation generally produces a torque, called the YORP effect, which may be important in re-orienting asteroidal spins. Here we explore a related effect, the “binary YORP” (BYORP), that can modify the orbit of a synchronously rotating secondary in a binary system. It arises because a locked secondary is effectively an asymmetric appendage of the primary. It should be particularly important for Near-Earth Asteroids (NEAs) owing to their small sizes, proximity to the Sun, and benign collisional environment. To estimate BYORP’s strength, we subjected 100 random Gaussian spheroids to the thermal radiation model of Rubincam (2000, Radiative spin-up and spin-down of small asteroids, Icarus, 148, 2–11). For most shapes, a significant torque arose on the secondary’s orbit, typically modifying the orbit’s size, eccentricity and inclination in less than 105 years, for components of 1 and 0.3 km radii, separated by 2 km, at 1 AU, each of density 1750 kg m−3 . Together YORP and BYORP are capable of synchronizing secondaries and circularizing orbits, making tidal dissipation unnecessary to explain the evolved state of observed NEA pairs. However, BYORP’s rapid timescale poses a problem for the abundance of observed NEA binaries, since their formation rate is thought to be much slower. We consider and reject the following resolutions of this quandary: (i) the approximation using Gaussian spheroids inadequately models YORP; (ii) most secondaries are not synchronous, but inhabit other spin-orbit resonances (very unlikely); (iii) tidal dissipation is much more efficient than previously estimated, and thus capable of stabilizing observed systems; and (iv) moderately close encounters with planets can re-orient secondaries, turning BYORP into a slower, random-walk process. Finally, we speculate that most observed binary NEAs are in a stable state in which the obliquity-changing torques of YORP (acting on the primary) and BYORP cancel out, and that those systems must be close to 55◦ inclination, where the momentum-changing torques of both YORP and BYORP tend to be very small. Some retrograde systems might develop such that the nodes precess at a Sun-synchronous rate, while some prograde ones might move into the “evection” resonance. All three of these hypotheses can be tested directly by comparison with the i, Ω and  observed for NEA binaries.  2005 Elsevier Inc. All rights reserved. Keywords: Asteroids; Asteroids, dynamics; Asteroids, rotation; Orbits; Satellites, general

1. Introduction In recent years the field of asteroid dynamics has been profoundly changed by wide acceptance of the importance of the Yarkovsky effect. This effect is based on the anisotropic re-emission of infrared radiation by an asteroid, * Corresponding author. Fax: +1 604 822 5324.

´ E-mail address: [email protected] (M. Cuk). 1 Now at: Department of Physics and Astronomy, University of British

Columbia, Vancouver, BC V6T 1Z1, Canada. 0019-1035/$ – see front matter  2005 Elsevier Inc. All rights reserved. doi:10.1016/j.icarus.2005.02.001

due to its non-zero thermal inertia (for a review, see Bottke et al., 2002). The Yarkovsky effect has so far been used to explain the dynamical diffusion of km-sized asteroids into resonances and the delivery of meteoroids to near-Earth space (Farinella et al., 1998; Bottke et al., 2000), and the dynamical spreading of asteroid families (Bottke et al., 2001; Carruba et al., 2003). It is also crucial for making predictions of NEA hazards even just a few centuries in the future (Giorgini et al., 2002). The Yarkovsky effect requires objects to have a non-zero thermal inertia, but it does not make any demands on their

Radiation effects on binary NEAs

shape. While the shape of the body has to be taken into account for more precise applications (Spitale, 2001), the Yarkovsky effect does not require that the asteroid has to be irregularly shaped, and generally acts even on spherical bodies. More recently, a related effect, usually referred ˇ to as YORP (Rubincam, 2000; Vokrouhlický and Capek, 2002), has been discovered to affect an asteroid’s rotation. While also relying on the “Yarkovsky force” (recoil from thermal radiation), YORP is effective only on highly irregularly shaped bodies, and operates even in the absence of thermal inertia. YORP has been credited with producing the clustering of pole positions for the Koronis family members (Slivan, 2002; Slivan et al., 2003; Vokrouhlický et al., 2003), as well as creating an overabundance of fast and slow rotators among small NEAs (Pravec and Harris, 2000). Another development that has dramatically changed the field of asteroid studies within the last decade is the confirmation of the existence of binary asteroids. While the possibility that some asteroids might have satellites had been considered for many years (Weidenschilling et al., 1989), only in 1993, when Galileo spacecraft discovered the satellite of 243 Ida, did the asteroidal satellites become reality. Since then, many binary asteroids have been discovered by ground-based observers using radar, adaptive optics, lightcurve analysis, as well as by the Hubble Space Telescope (Merline et al., 2002). Binaries are known to exist among almost all dynamical classes of minor planets: NEAs, main-belt asteroids, Trojans and Trans-Neptunian Objects (TNOs). However, moving from one class of object to the

419

other, the physical size of components in known binary systems varies widely: for example, the binary NEAs’ primaries are typically only a kilometer or so across, while those of binary TNOs have radii in excess of 100 km (Burns, 2002; Margot, 2002). Since this work is concerned with the effects of radiation on the mutual orbit of the binary members, we will naturally concentrate on binary NEAs, since they are the only known population of binaries that falls into size range for which Yarkovsky and YORP effects are important. Table 1 (P. Pravec, private communication, 2004) summarizes some physical and orbital characteristics of all known binary NEAs. The typical NEA binary has a fast spin period (P = 2–3 h), and a kilometer-sized primary that is often nearly spheroidal in shape. The secondary is usually a few times smaller in radius, orbits at 2–4 primary diameters on an almost-circular path, and almost always rotates synchronously (or is, in case of radar observations, consistent with such a state). The fast rotation and regular shape of most primaries have usually been ascribed to their being “rubble piles” spinning close to the rotational disruption limit (Richardson et al., 2002). Such bodies are conducive to break-up, caused either by YORP-induced spin-up or (more likely) tidal forces of the terrestrial planets during close approaches (Bottke and Melosh, 1996a, 1996b; Richardson et al., 1998). The low orbital eccentricity and synchronous rotation of the secondary are usually explained by tidal dissipation, by analogy with planetary satellites (Murray and Dermott, 1999), although the efficiency of tidal processes decreases sharply with decreasing size of the components

Table 1 Orbital and physical parameters of known NEA binaries (P. Pravec, private communication, 2004). The columns list the minor planet number (if assigned), name or temporary designation, primary’s diameter, secondary-to-primary diameter ratio, separation (in primary diameters), eccentricity, orbital period and the primary’s spin period, respectively. A more current version of this table can be found on P. Pravec’s Binary NEA Internet site (http://www.asu.cas.cz/~asteroid/binneas.htm); the table shows the values given on September 24th, 2004 Number 35,107 3671 5407 31,345

66,391

5381 66,063 69,230 65,803

Object

d p (km)

ds /dp

a/dp

e

P orb (h)

P rot (h)

1994 AW1 1991 VH Dionysius 1996 FG3 1992 AX 1998 PG 1999 HF1 2000 DP107 2000 UG11 1999 KW4 1998 ST27 2001 SL9 2002 BM26 2002 KK8 Sekhmet 1998 RO1 2003 SS84 Hermes 1990 OS Didymos 1999 DJ4 2003 YT1 2002 CE26

0.9 1.2 0.9 1.4 4.0 0.9 3.5 0.8 0.23 1.2 0.8 1.0 0.6 0.5 1.0 0.8 0.12 0.4 0.3 0.8 0.4 1.0 3

0.53 0.40 >0.28 0.31
0.30
0.30 0.24 0.38
0.6
0.3 ∼0.13 0.31 ∼0.2 ∼0.2 0.3
0.5 0.5 1.0 0.15 0.2 0.5 0.18 0.07

2.3 2.7 2.6 1.7 (1.7) (1.7) 2.0 3.3 1.8 2.1 (∼5) 1.8 (5)

< 0.05 0.07

22.40 32.69 27.72 16.4 (13.52) (14.01) 14.02 42.2 18.4 17.45 (∼100) 16.40 72

2.52 2.62 2.71 3.59 2.55 2.52 2.32 2.78 (4.44) 2.76 3.1 2.40 ∼2.7

(12 or 24) 14.53 (24) (13.89) (18–24) 11.90 17.72 ∼30 ∼16

(∼3) 2.49

1.5 (1.7)

(1.5) (2.0) (2.7) 1.7

0.05 (<0.05)

0.01 0.03

13.89 Hours 2.26 2.51 2.34 3.29

420

´ M. Cuk, J.A. Burns / Icarus 176 (2005) 418–431

(Weidenschilling et al., 1989). While the rotational dynamics of single bodies in the same size range are today thought to be dominated by the YORP effect, no investigation to date has explored the potential importance of YORP in binary systems. In general, radiation effects on the orbit of a planetary satellite are quite different from those on a body orbiting the Sun. The classical Yarkovsky effect disappears, when averaged over the mean motions of both the planet and the satellite, as there is no preferred angle between the direction to the Sun and the satellite’s orbital velocity. Nevertheless, a number of other radiation-induced perturbations are specific to planetary satellites, both natural and artificial. Before discussing the more general problem, we need to mention the “pathological” case of retrograde quasisatellites, which are bodies that inhabit an interesting coorbital state with their primaries, with the center of their librations coinciding with the primary. In the reference frames fixed on their primaries, the quasi-satellites appear to follow an elongated, retrograde orbit around the primary with a period equal to that of the primary’s heliocentric orbit (see Namouni, 1999). Since their orbits are resonant with the Sun, the Yarkovsky effect on the quasi-satellites is likely non-trivial and may lead to interesting secular effects. Similarly, objects librating around linear Lagrangian points of a primary would also feel a non-zero net Yarkovsky effect, as the Sun appears stationary relative to the system (J. Spitale, personal communication, 2004). Both of these cases are interesting, but lie outside of the scope of this paper. To date, the most detailed studies of radiation effects on a satellite orbit have been those dedicated to explaining the results from the LAGEOS mission in the 1980s (Rubincam, 1987, 1988; Milani et al., 1988). The most important radiation-driven perturbations on LAGEOS’s orbit are the Rubincam and Yarkovsky–Schach effects. The Rubincam effect is just the Yarkovsky effect powered by infrared radiation from Earth rather than by visible light from the Sun. The Yarkovsky–Schach effect is a result of solar eclipses by the planet, as seen from the satellite. Since eclipses always happen in approximately the same configuration of all three bodies, the sunlight is no longer azimuthally symmetric, and their effects end up being additive rather than averaged out, producing a net Yarkovsky force. However, neither the Rubincam nor the Yarkovsky– Schach effect is likely to influence secondaries in binary asteroid pairs, which usually exhibit synchronous rotation. The essential aspect of the Yarkovsky and associated effects is that surface elements must absorb and re-emit radiation in different rotational phases, so that the momenta carried by the incoming and outgoing radiation make different angles with the object’s orbital velocity. Obviously, any given surface element on a synchronously rotating body will always have the same orientation relative to the body’s motion, so there cannot be any of the “thermal lag-angle” that is fundamental to the Yarkovsky effect, as well as to the Rubincam and Yarkovsky–Schach effects. The only caveat to this state-

ment is that we assumed that the secondary has a spherical shape, which is probably unrealistic. The purpose of this work is to show how the Yarkovsky force, despite not producing secular effects for synchronous spherical asteroidal satellites, can be very important for bodies with more realistic shapes. ´ Cuk et al. (2002) have described how the most distant satellites of planets with large Hill spheres should also experience a form of the Yarkovsky effect, due to an appreciable gradient of sunlight over their orbits. This case is relevant for the tiny irregular moons of Jupiter as well as any hypothetical satellites of massive close-in extrasolar planets ´ (Burns and Cuk, 2002), but has no useful applications for asteroids, since their Hill spheres always have small absolute sizes (Hamilton and Burns, 1991).

2. The mechanism To explain the principles behind the “binary YORP effect” (BYORP), we start by describing how the usual YORP effect (Rubincam, 2000) operates. To do so, we illustrate how a body’s geometry can affect the radiation forces that are experienced. In this discussion, for the sake of simplicity, the planes of the “asteroid’s” orbit and its equator are assumed to coincide, so that the Sun always remains in the equatorial plane. Fig. 1 shows an idealized, “irregularly shaped” rotating asteroid, which consists of a sphere with two identical equilateral, right triangular prisms fastened 180◦ apart along the sphere’s equator (Fig. 1a). Each prism is attached to the sphere at its base, which is an equilateral right triangle. Each prism has its shorter sides horizontal and vertical. The hypotenuse of one prism points up and right, whereas the other’s hypotenuse faces up but to the left (see Fig. 1b). Both prisms have the same cross-sections whether viewed from the left or the right, so the amount of sunlight absorbed by either prism from two diametrically opposite directions will always be the same. Each photon of sunlight carries momentum. Thus, by the impulse-momentum theorem, a force is imparted when the photon strikes the body and is absorbed (cf. Burns et al., 1979); depending on the part of the surface that a photon strikes, it likely also creates a torque on our model body. Nonetheless, for a body shaped like our idealized asteroid, the net torque from the light absorbed by the two prisms during one synodic rotational period is zero, because the contributions at phases 180◦ apart cancel each other. The “asteroid’s” spherical part, due to symmetry, obviously contributes torque at no time. Hence the entire asteroid experiences no net torque from the absorbed momentum. Now we consider the radiation emitted by the body; this re-emission is required if the body is in thermal equilibrium. We will assume that the radiation is instantaneously re-emitted according to Lambert’s law, so that the total momentum of the outgoing radiation is always perpendicular to the surface element from which it is emitted. Once again,

Radiation effects on binary NEAs

(a)

(b) Fig. 1. Schematic view of how YORP operates on an idealized “asteroid” consisting of a sphere with two triangular prisms (panel (a)) attached to the equator (see text for discussion). Panel (b) shows how radiation leaves the surfaces of the two prisms that lie on either side of the sphere. This figure combines Figs. 1 and 2 from Rubincam (2000).

due to symmetry, the spherical part of the “asteroid” makes no contribution to the average torque on the body. However, now the prisms behave differently when their different sides are lit. The bottom panel of Fig. 1b shows the prism’s narrower side being struck perpendicular to its face by solar radiation (represented by the wiggly arrows); by our assumption, the re-radiation of the absorbed energy and momentum (represented by the block arrow) occurs directly backwards in the same plane, yielding a net force on the prism from the interaction that lies totally in the orbital (equatorial) plane. The resultant torque on the body only has a component normal to the asteroid’s orbit plane. But 180◦ later, when

421

the Sun illuminates the prism’s slanted side (top panel of Fig. 1b), energy and momentum are re-radiated out of the equatorial plane. Accordingly, some fraction of the incident momentum has been transferred from the equatorial plane to the direction perpendicular to the equator. Of course, since the overall momentum flux is conserved, any modification in the photon’s momentum flux has to be compensated by a force on the prism and thence on the body to which it is an appendage. When averaged over all rotational phases, this force produces a torque capable of changing the “asteroid’s” spin. More generally, when the equator and the orbit do not coincide, the YORP effect produces a vector torque; that vector can be decomposed into one component which changes rotational spin, another which alters obliquity, and a third which causes the rotational pole to precess. Any directly reflected light will also produce a similar effect; this will be ignored here as asteroid surfaces tend to have low albedos, absorbing significantly more energy than they reflect. We now address what would happen if the prism were not attached to the sphere, but rather located at some distance from it, while keeping the same orientation as in Fig. 1, with one base facing the sphere. Obviously the radiative forces and resulting torques would affect only the prism, which would experience a net torque around the sphere’s center. This abstract geometry is actually quite similar to that of a binary asteroid having a nearly spherical primary and a synchronously rotating, tidally locked irregularly shaped secondary. Of course, a real asteroid would not be prismshaped, but nevertheless, due to the asymmetry of its facets, it would similarly experience a net torque like the one described for our toy model. It is obvious that a non-synchronously rotating satellite would, at any given time, have a random orientation relative to its orbital motion. A synchronous one, on the other hand, will have permanent leading and trailing hemispheres, making the “promontory” analogy appropriate. Therefore, the effect discussed here (which we will refer to as the “Binary YORP” effect or simply BYORP) can produce a secular acceleration only for synchronously rotating asteroidal satellites, whereas it averages to zero in non-synchronous cases. The torque on the secondary would alter the angular momentum of its orbit around the system’s center of mass and also its own spin. Once again, we can expect a vector torque to arise, and its three perpendicular components would change the spin period and obliquity, as well as induce precession of the line of nodes, respectively. Since the magnitudes of all three torques depend quite sensitively on the secondary’s shape, no reliable estimate of the effect’s order of magnitude can be obtained through these basic considerations. In order to determine the potential significance of the BYORP effect, we must model it numerically, using more realistic shapes, similar to those of real-world asteroids.

´ M. Cuk, J.A. Burns / Icarus 176 (2005) 418–431

422

3. Numerical experiment In order to estimate how BYORP should affect secondaries in asteroids, we model the objects as Random Gaussian Spheroids (RGS), described by Muionen (1996) and previously used for exploring YORP by Vokrouhlický ˇ and Capek (2002). The RGS technique is based on representing an irregular shape as a superposition of many spherical harmonics, or more precisely, by having the natural logarithm of the radius be the harmonics’ weighted sum s(θ, φ) =

∞
l


Plm (cos θ )(alm cos mφ + blm sin mφ), (1)

l=1 m=1

with log-radius s(θ, φ) defined as   s(θ, φ) = ln 1 + σ 2 r(θ, φ)/a   = β 2 /2 + ln r(θ, φ)/R ,

(2)

where R is the mean radius, σ and β are variances of radius and log-radius, respectively, and Plm are associated Legendre functions. The spherical harmonics’ coefficients alm and blm are generated using a pseudo-random number generator (i.e., one that is not truly random but the output of which is statistically indistinguishable from a truly random sample). They were generated as part of a Gaussian distribution with dispersion βlm , that was chosen so that the resulting log-radii of the surface points have a normal distribution with variance β 2 βlm = (2 − δm0 )

(l − m)! 2 cl β , (l + m)!

(3)

with the autocorrelation coefficients cl derived from modified spherical Bessel functions il , following Muionen (1996) cl = (2l + 1) exp(−κ)il (κ),

of “pole-pointing” triangles were defined by pairs of neighboring points which were on the same parallel, while the third corner was found by linear interpolation between the two closest points on the next parallel from the equator. Then the “equator-pointing” triangles were just defined as being complementary to the “pole-pointing” ones. The equatorial points defined two sets of “pole-pointing” triangles, while the poles themselves were the top corners of multiple “pole-pointing” triangles, without any “equator pointing” ones in between. We then find the radii of the triangles’ center points, their areas, and the surface-normal vectors for all 2430 surface elements. Having the data on our artificial shape in this format allows us to find its center of mass, by dividing the volume into 2430 triangular pyramids with surface elements as bases, and assuming uniform density. This is needed as the center of coordinates in which the body was generated is not necessarily its center of mass for any one of the artificial shapes, but just the center of mass of the whole ensemble (since our technique should create no systematic anisotropy in the resulting shapes). After the center of mass is found, we convert the central radii of our surface elements into the new reference frame, centered on the center of mass (of course, the areas and surface-normal vectors stay the same during translation). Then we find the principal moments of inertia of our bodies, once again choosing the triangular pyramids to have uniform density. For each pyramid whose vertex is at the body’s center of mass, and whose height is much larger than any of the base’s sides, the integral giving the moment of inertia integral is  R

1 3 xi xj r 2 dr dΩ = Xi Xj R 3 Ωa = Xi Xj Va , 5 5

(6)

Ωa 0

(4)

where il is a modified spherical Bessel function, while parameter κ measures the angular irregularity of the object’s surface, and is tied to the correlation angle Γ −1  2 Γ κ = 4 sin (5) . 2 While the process of generation of an artificial shape is somewhat complex, it requires only two parameters as input: σ and Γ , which characterize the scales of the radial and the angular variabilities, respectively. We used σ = 0.274 and Γ = 30.9◦ , following Muionen and Lagerros (1998) and ˇ Vokrouhlický and Capek (2002); and we also introduced a cutoff in the spherical harmonics employed (Eq. (1)) at l = 10, once again to make our results more easily comparaˇ ble to those of (Vokrouhlický and Capek, 2002). The objects are generated in form of 1162 points, distributed approximately equally across all solid angles (the “latitude” interval was constant, but the one in “longitude” increased toward the “poles”). The next step was to convert the points into triangular surface elements. The bases

 = (Xi , Xj , Xk ) is the position of the surface elwhere X ement’s center, Ωa is the solid angle it subtends from the center of mass, and Va is the volume of the pyramid it defines. After computing all the elements in the inertia tensor, we find its eigenvalues (which are the principal moments of inertia C > B > A); in doing so, we have diagonalized the matrix (i.e., we have transformed to the principal axes, the system of coordinates where the inertia tensor has only diagonal components, see Goldstein et al., 2002). In these new axes, we define the z-axis to be that of the maximum moment of inertia, the x-axis to provide the intermediate moment and y-axis to be that of the minimum moment, chosen so that zˆ × xˆ = y. ˆ Since we assume uniform density for our artificial asteroids, the z-, x- and y-axes, typically, lie along the shortest, intermediate and longest diameters of our body (this is not always guaranteed as our body is not an ellipsoid). It was essential to determine the principal axes of inertia because we will assume that the secondary has tidally evolved to a synchronous state rotating about its axis of largest moment of inertia (Peale, 1977) with its long axis pointing toward

Radiation effects on binary NEAs

the primary. Of course, the 180◦ ambiguity in the choice of positive x and z directions is inherent to the physics of the problem; later we will discuss the implications of this degeneracy on the evolution of real NEA binaries. The next step in building our model is to introduce the Sun as the source of radiation. In order to describe the interaction between the Sun and the rotating and revolving binary asteroid, we introduce an “inertial” reference frame centered on the primary. While not strictly inertial (as it is not tied to the center of mass of the Solar System), the axes of this frame do not rotate, but stay fixed relative to the background stars. In this frame, not only is the secondary revolving around the primary while rotating synchronously, but the Sun is also considered to orbit the pair, on an inclined orbit. We decided to account for the yearly motions in this way for simplicity, as no important physics in the problem gets affected: most importantly, the line of nodes between the binary’s mutual and heliocentric orbits (the latter being the plane of the Sun’s apparent motion) stays constantly oriented in space. We define the inertial coordinates ξ , η and ζ as follows: ξ -axis points toward the ascending node of the Sun’s apparent orbit, while ζ points along the angular momentum vector of the binary (we assume that the secondary’s axis of rotation is normal to its orbit, which lies in the ξ –η plane). At t = 0, the x, y and z axes of the coordinate system are fixed to the secondary and are, respectively, parallel to the ξ , η and ζ axes of the inertial frame. The secondary’s orbital velocity is always in the positive x direction, while the primary will always be at y = a in the positive y direction. From the primary’s point of view, the secondary will be at η = −a and have velocity along ξ axis at t = 0. At t = t  , the secondary will be at ξ = sin ωt  and η = − cos ωt  . In order to find the radiation’s net effect on the secondary’s motion during a complete orbit of the binary (but with the fixed position of the Sun), we add together 36 evaluations of the radiative force computed at 10◦ intervals. On the other hand, in order to find how this effect changes in the course of the asteroid’s year, we simply move the Sun along its apparent path. To find the net annual effect, we evaluate the radiative force over the whole binary revolution (as described in the above paragraph) for each of the 36 positions of the Sun, which are equally spaced along the “ecliptic.” We ignore the precession of the binary in the course of one year. The last approximation is not always safe; for example, a binary having a separation of six primary radii, an orbital period of one day and J2 = 0.1 for the primary (this is realistic for barely stable rubble piles) would have its line of nodes precessing once in a year or so. Still, averaged effects of precession should be zero as long as it is not commensurable with other frequencies in the system (see Section 5). To ascertain the total force acting on the asteroid, we will find the differences between the momentum fluxes of the incoming and outgoing radiation, for each of the surface el-

423

ements and then sum over the entire body F (ψ, φ) = −


ΦSi   nˆ i (ψ) · sˆ (φ) nˆ i , c

(7)

i

where F (ψ, φ) is the total Yarkovsky force on the body for binary position ψ and solar position φ, Φ is the solar flux, c the speed of light, Si the area of the ith element, nˆ i (ψ) its surface normal and sˆ (φ) is the direction to the Sun. The program has been instructed to set F = 0 for all elements for which nˆ i (ψ) · sˆ (φ) < 0 (i.e., Yarkovsky force is zero during the night since our model has no thermal inertia). After that, we sum the torque caused by F (ψ, φ) over 36 discrete ψ s (one whole “month”) for each of 36 discrete φ’s (one whole “year”): Ttotal =

35
35


a (ψk ) × F (ψk , φj ),

(8)

j =0 k=0

where ψk = 2πk/36, φj = 2πj/36 and a (ψk ) is the secondary’s radius-vector. The resulting torque has, in general, three components: Tζ parallel to the binary’s angular momentum, Tξ along the line of nodes and Tη normal to the first two. Only Tζ changes the binary’s angular momentum, while Tξ and Tη cause precession and alter the binary’s inclination, respectively. Two approximations made in the above calculation are certainly not completely accurate: we assumed that the body has no thermal inertia and that there is no shadowing. Our rationale for ignoring thermal inertia is that the BYORP effect (just like YORP) does not require it, and thus in this first discussion of the phenomenon, we can ignore it. However, if we would want to make useful predictions about any body’s exact behavior under BYORP, we would need to take thermal inertia into account. In the case of YORP, Vokrouhlický ˇ and Capek (2002) have not treated thermal inertial directly, while Vokrouhlický et al. (2003) have, and their results differ significantly. Vokrouhlický et al. (2003) find that the relative abundance of asteroids exhibiting different types of behavˇ ior (designated “I–IV” by Vokrouhlický and Capek (2002)) changes once thermal inertia has been included. Therefore, in this work, we will not assign any percentages to different types of behavior among Gaussian spheroids, as such assignments are likely to be incorrect. We will instead focus on the order-of-magnitude strength of different processes acting on binary NEAs. ˇ Vokrouhlický and Capek (2002) used a “Rubincam approximation” to treat thermal inertia, which basically means including a factor of 2/3 in Eq. (7). We have not done that as we are interested only in the order of magnitude and wished to keep as simple a model as possible. However, it would be trivial to apply this approximation to the final result if needed, as it would simply require one to multiply all results by 2/3. The other approximation we used—that shadowing could be overlooked—is even less justified. It certainly makes us

424

´ M. Cuk, J.A. Burns / Icarus 176 (2005) 418–431

underestimate the overall diversity of BYORP-related behavior among binary NEAs. However, our computational and personnel constraints have made it impractical to expand our model to incorporate shadowing, as it is noˇ toriously computation-intensive (Vokrouhlický and Capek, 2002). Given that any shadowing is unlikely to change the order of magnitude of BYORP, we decided that the present model without shadowing would be sufficient. Once highresolution models of real secondaries are available from observations, the inclusion of realistic thermal inertia and shadowing may become worthwhile.

4. Results

(a)

Before presenting the results of the numerical experiment for BYORP, we first test our software by modifying it to calculate YORP rather than BYORP (the changes in the code are quite minor). This allows us to then compare ˇ our results with those of Vokrouhlický and Capek (2002). Our program computes the torques in dimensionless units, in which R = 1, ρ = 1 and Φ/c = 1 (ρ is the asteroid’s mean density). In order to convert the result into SI units (meters, seconds and kilograms), we need to multiply the dimensionless results (Ti /C  ) by a combination of constants T Φ Ti = i . C C cρR 2

(9)

We avoided making this step automatic, as it gives us the opportunity to apply the same numerical results to bodies with a range of sizes, mean densities and heliocentric distances. Although later we mostly discuss sub-km NEA secondaries, Fig. 2 displays results with a = 1 km, ρ = 2500 kg m−3 and rSun = 2.5 AU, for purposes of comparˇ ison with (Vokrouhlický and Capek, 2002). Fig. 2a shows the strength of the secular acceleration of the spin rate for ten bodies in our sample (Tζ /C), while Fig. 2b plots the intensity Tη /C of the component of angular acceleration that changes the obliquity. The results can be directly compared ˇ to those in Fig. 11 of Vokrouhlický and Capek (2002), as −18 they are plotted in the same units (10 s−2 ). Considering ˇ that the results of Vokrouhlický and Capek (2002) containing a multiplicative factor of 2/3, the sets of curves appear identical (we plotted all four types of behavior together in both Figs. 2a and 2b). Now that we have assured ourselves that our software is generally accurate for YORP, we move on to computing BYORP. Just like in the case of YORP, our raw results are dimensionless and need to be multiplied by a set of constants T T Φ P =  , L L cρRa 2π

(10)

where a is the binary’s semimajor axis, L is its orbital angular momentum and T is its period. We expressed this result as the torque over angular momentum, rather than torque

(b) Fig. 2. (a) Angular acceleration due to ζ -component of the YORP torque. Each of the lines shows the results for one of ten random Gaussian spheroids. This component of the torque acts to change the magnitude of angular momentum and therefore the spin period of the body. The units of y-axis are 10−18 s−2 . (b) Angular acceleration due to η-component of the YORP torque for the same subset of bodies. This component acts to change the body’s obliquity of the body. The units of y-axis are 10−18 s−2 .

over moment of inertia, as an orbiting body’s moment of inertia changes together with the orbit’s angular momentum. We express the results in terms of observables a and T rather than the primary’s mass, just to make comparison with Table 1 easier. In plotting our results, we assume Rp = 0.5 km, R/Rp = 0.3 and a/Rp = 4, or R = 150 m and a = 2 km. For Φ, we use the solar constant, Φ0 = 1368 W m−2 , since we are dealing with NEAs. We employ an orbital period of 20 h, which implies a density of 1750 kg m−3 for the primary, so we also adopted the same density for the secondary. This density indicates significant porosity for our artificial asteroids, since the density of solid rock is closer to 2500– 3000 kg m−3 . This is probably justified for the primaries, which are widely thought to be rubble piles prone to disruption, but may not always hold for the secondaries, which could conceivably be monolithic. The assumption of equal densities is supported by radar observations of the NEA bi-

Radiation effects on binary NEAs

Fig. 3. Relative rate of change of the binary’s angular momentum due to BYORP. This quantity is the inverse of the angular momentum’s doubling timescale. The units about the y-axis are 10−12 s−1 , and 1 on this scale is equivalent to an evolution timescale of about 3 × 104 yr.

nary 2000 DP107 (Margot et al., 2002) which suggest both significant porosity for the primary and a similar mean density for both components. Fig. 3 plots the inverse orbit modification timescale due to the torque parallel to the orbit normal, Tζ /L, in units of 10−12 s−1 for all possible inclinations of 100 artificial binary systems. We see that the effect is symmetric for prograde and retrograde binaries, as our model assumes that the binary’s orbital motion is much faster than its heliocentric one. Note that in this case reversing the direction of the orbital motion is not equivalent to changing the inclination by 180◦ , because such a reversal would also make the leading and trailing hemispheres exchange places. In Fig. 3, a body with a retrograde inclination still has the same leading hemisphere it had when i < 90◦ . However, as mentioned in Section 3, the choice of leading and trailing hemispheres is arbitrary. So the behaviors of our artificial bodies should be viewed as symmetric relative to the abscissa. Overall, the behavior of BYORP with inclination is very similar to that of YORP. We found only eight bodies for which the acceleration was zero for more than one inclination between 0◦ and 90◦ , while all the others had a single zero, typically in the i = 50◦ –60◦ range, just like YORP (see Figs. 4 and 5 in Rubincam, 2000). Fig. 4 shows the corresponding behavior of the inverse for the inclinationchanging timescale, Tη /L, in units of 10−12 s−1 . As expected, the relevant torque is anti-symmetric around 90◦ . It is easy to see that, in the system in which the Sun is treated as a uniform, inclined ring around the binary, the case with i = β is equivalent to that with i = 180◦ − β, just with the positive direction of ξ -axis reversed. The most important result of the numerical experiment concerns BYORP’s strength. An intensity of 1 × 10−12 s−1 implies a doubling time for the angular momentum of about 3 × 104 yr. If the torque is negative, this gives the approximate timescale on which the secondary will impact the primary. The intensity of the effect is proportional to

425

Fig. 4. Rate of change of the inclination of the binary pair due to BYORP. The units of y-axis are 10−12 rad s−1 , and 1 on this scale is equivalent to an evolution timescale of about 3 × 104 yr.

√ T /a ∼ a, so outward migration should be a runaway process while inward drift should slow asymptotically. So a naive view of BYORP would be that it inevitably causes either collision or separation of the components after just a few tens of thousand years. However, tidal disruption and non-synchronous rotation can put bounds on the inward and outward motions, respectively. It is clear that, once the secondary moves inside the Roche limit of approximately 2.5 primary radii, it will tidally disrupt unless it possesses material strength. Monolithic secondaries are, in principle, able to hold themselves together by material forces, but the lack of known pairs with a/dp < 1.25 in Table 1 suggests that most secondaries avoid being too close to the primary. While observational biases might lead to some very close binaries being misclassified as single elongated objects, such pairs do not appear to be common among NEAs. The breaking of synchronous lock is somewhat more uncertain, but is bound to happen sooner or later as the secondary drifts to larger separations. In addition, small-body impacts and planetary fly-bys are both capable of disrupting the synchronous lock, with the second process usually considered more relevant for NEAs. However, for the large majority of bodies, which show either Type I or Type II behavior, to use the designations ˇ of Vokrouhlický and Capek (2002), the breaking of synchronous lock is not required to stop (or even reverse) a secondary’s outward evolution. All bodies of “Type I,” those that show descending curves in Fig. 2 (i.e., whose BYORP acceleration decreases with inclination), also have positive values of Tη /C (Fig. 3). Similarly, bodies with Type II BYORP behavior have ascending curves in Fig. 2 and negative ones in Fig. 3. Together, these two types account for over 90% of our artificial bodies. In both cases, the obliquitychanging torque Tη tries to move the obliquity into the region where Tζ is negative. So even if the initial obliquity puts the binary in the Tζ > 0 regime, eventually the inclination’s evolution should drive the binary into the state where Tζ < 0,

426

´ M. Cuk, J.A. Burns / Icarus 176 (2005) 418–431

following which the slow inward spiral of the secondary appears unavoidable. This result, taken at face value, appears inconsistent with both observations (Margot, 2002) and theoretical predictions of binary formation rates (Bottke and Melosh, 1996a). About 15% of NEAs appear to be binaries, and this is roughly consistent with the formation rate expected from close passages by Earth and Venus. Without BYORP, the lifetimes of binaries are mostly limited by their dynamical lifetimes as NEAs; hence most binaries should live for millions of years (Bottke and Melosh, 1996a). However, if our predictions of BYORP are accurate, almost every NEA has to go through a binary phase at least once in order for us to observe as many as we do. In the next section, we will show that, when both components and all physical processes acting on them are included, the short lifetime problem might be averted. Before discussing the long-term orbital behavior (semimajor axis and inclination), we mention that both BYORP and YORP also affect the secondary’s rotational period and orbital eccentricity. It is well known that tangential torques acting at one apse of the orbit tend not to affect that distance, but rather alter the position of the other apse. For example, the gas drag felt by a satellite at pericenter tends to decrease its apocenter distance, while leaving the pericenter unchanged (Burns, 1976). Even when gas drag influences both apses, the typically greater density of gas at pericenter tends to shrink the apocenter distance faster and therefore decreases the orbital eccentricity. We find a somewhat similar behavior in the case of BYORP: its force is independent of the distance from the primary, but its torque is stronger at larger distances. So, if Tζ is positive, we would expect it to be greater at apocenter and less at pericenter. Accordingly, the pericenter distance will grow faster than that of apocenter, and the eccentricity will decrease. If Tζ < 0, the pericenter will collapse faster than the apocenter and the eccentricity will increase. If the binary is “fresh,” i.e., just created during a planetary encounter, it is expected to have significant eccentricity and relatively small separation. It is clear that only binaries with Tζ > 0 at that epoch can avoid disaster, and evolve toward orbits having larger semimajor axis and lower eccentricity. To estimate the timescale of eccentricity damping, we use the appropriate equation (Burns, 1976) for the evolution of perturbed orbits  1/2   de   = a 1 − e2 µ−1 FR sin f + FT (cos f + cos E) , dt (11) where e is eccentricity, µ = Gmp is the primary’s mass multiplied by the gravitational constant, FR and FT are radial and tangential perturbing forces per unit mass (i.e., they are accelerations), and f and E are the true and elliptic anomalies, respectively. Using cos f = (cos E − e)/(1 − e cos E),

and dE = n dt/(1 − e cos E), where n is the mean motion (Murray and Dermott, 1999), and integrating over a full orbit, we get 

   1/2 3   3 eFT de = − e a 1 − e2 µ−1 , FT  − (12) dt 2 2 v where v is the orbital velocity, and the eccentricity has been assumed to be small in the last expression. Clearly, due to BYORP, e/e ˙ operates on a similar timescale to a/a. ˙ The binary’s eccentricity should, of course, also evolve due to the tidal forces. The timescale of tidal eccentricity damping due to tides in the secondary (Murray and Dermott, 1999) is given by   e 38 ms a 5 µr Q , τe = − = (13) e˙ 63 mp Rs nρ 2 GR 2 where Q is the material’s quality factor and µr is the secondary’s rigidity. We choose µr Q = 1011 N m−2 , based on the measurements on dissipation in Phobos (Yoder, 1982). Using our standard binary parameters (a/R = 13.3, R/Rp = 0.3, R = 150 m, n = 2π/(20 h) and ρ = 1750 kg m−3 ), we find τe  5 × 1010 yr, or six orders of magnitude longer than that obtained for e/e ˙ due to BYORP. While the numerical value of the secondary’s µr Q is basically unknown, and our estimate could be significantly in error, this comparison suggests that, in most cases, BYORP should damp eccentricity more rapidly than tidal forces. Therefore the observed low eccentricity of many binary NEAs does not necessarily mean that the system underwent any tidal evolution, as previously suggested (Merline et al., 2002), because BYORP alone can achieve the same qualitative result. Another process that may have been mis-attributed to tidal forces is the secondary’s spindown. Now many planetary scientists suspect that the spins of asteroids smaller than 1 km are dominated by YORP, especially if they are in a benign collisional environment like the inner Solar System. If we take 10−18 s−2 to be an average Tζ /C due to YORP for a body with R = 1 km, rSun = 2.5 AU and ρ = 2500 kg m−3 , then for our “standard secondary” (see above) the same quantity should be, according to Eq. (9), Tζ /C = 4 × 10−16 s−2 . Assuming an initial spin period of 3 h (Pravec and Harris, 2000), we compute τY  50,000 yr. The relevant expression for tidal spindown of a satellite is (Murray and Dermott, 1999)   9 ρ 2 GR 2 mp R 3 2 n , Ω˙ = − (14) 38 αµr Q ms a where α = C/(ms R 2 )  0.4. Using the same values as above, we get Ω˙ = 3 × 10−21 s−2 . If the initial spin period were about 3 h, then tides would take about 6 × 109 yr to synchronize the secondary. Even if we assume that the rigidity of rubble piles is much, much lower than that for Phobos (whose density of 1900 kg m−3 is not much higher than the

Radiation effects on binary NEAs

one we adopted), it appears unrealistic that such a process could be relevant on timescales shorter than the dynamical lifetime of the NEAs. Once again, our results strongly suggest that the secondaries’ slow spins are mostly due to YORP rather than tidal forces. Once the spin period becomes very long, the satellite’s permanent quadrupole moment can become locked by the primary’s tides. Unlike tidal dissipation, first-order tidal interactions scale with radius, so there is no reason to doubt the possibility of a tidal lock. Also, Vokrouhlický et al. (2003) find that, once thermal inertia has been taken into account, Type II behavior becomes more common, meaning that the obliquities move to 0◦ or 180◦ , while the spin slows down. Since a tidally created secondary is likely to have its spin aligned with that of the primary (as well as their mutual orbital plane; Richardson et al., 1998), its obliquity (relative to the Sun) is likely to move to 0◦ if the primary is prograde or 180◦ if the primary is retrograde. After tidal locking, internal dissipation in the secondary is likely to damp out both longitudinal and latitudinal librations and to exactly align the secondary’s spin with its orbit (Peale, 1977).

5. Discussion Our results for the secondary’s eccentricity and rotation clearly demonstrate that radiation effects alone can explain the observed abundance of synchronously rotating secondaries on low-e orbits. Therefore, the related effects on semimajor axis and inclination likely operate, too. Thus their rapid timescales must somehow be reconciled with the observed high incidence of binary NEAs. Before analyzing other processes that might affect the long-term orbital evolution of binary NEAs, we must consider the possibility that our model of BYORP is completely inaccurate. Perhaps real secondaries cannot be approximated by Gaussian spheroids. While we cannot prove that our approach is valid, we mention that the little that is known about the shapes of secondaries indicates that they are appreciably irregular (J.-L. Margot, private communication, 2004); thus a significant YORP-like force is very probable. It is conceivable that real secondaries might exhibit more complex BYORP behaviors than those that parallel Vokrouhlický and ˇ Capek (2002) Types I and II for YORP; e.g., these might include stable states toward which the binary might asymptotically evolve. However, Vokrouhlický et al. (2003) find that simple Type II behavior is even more common whenever a more sophisticated model is used, and given how closely the two effects resemble each other, this likely applies to BYORP, too. Another factor that could prevent BYORP from operating is any asynchronous rotation of the secondary. Radar observations of binary NEAs can only indicate that the secondary’s spin rate is slow and compatible with synchronous rotation, but they cannot determine the period accurately. Even though some bodies can evolve into non-synchronous

427

stable rotation states associated with spin-orbit resonances (Mercury being the archetype), for such a state to be stable, the secondary’s orbit must have substantial eccentricity (Murray and Dermott, 1999), which is inconsistent with the orbits of most known NEA binaries (Table 1). Currently, among tens of known binary NEAs, there is just a single good candidate for a spin-orbit resonance (J.-L. Margot, personal communication, 2004). Additionally, whenever lightcurve studies are able to distinguish a secondary’s lightcurve variation, it typically has a period equal to the orbital period deduced from eclipses (Pravec et al., 2005, in preparation). The phase of the secondary’s lightcurve also usually shows maxima halfway between the eclipses, as expected if the long axis of the body pointed toward the primary. So, according to current observations, widespread non-synchronous rotation of secondaries in NEA binaries is extremely unlikely. While it appears that there is no way to make BYORP disappear completely, it might be possible to stabilize some binaries by some other process counteracting BYORP, e.g., tidal dissipation inside the primary (as opposed to that inside the secondary discussed earlier). Tides raised on the primary by the secondary tend to drive the secondary out of synchronous orbit. In the case of most binary NEAs, this means that the secondary moves outward, while the spin of the primary slows down. The rate of change of the secondary’s semimajor axis (Murray and Dermott, 1999) is   9 ρp2 GRp2 ms Rp 5 a˙ = n, (15) a 19 αp µp Qp mp a where the quantities with subscript p have the same meaning as in Eq. (14), but now they describe the properties of the primary, which was taken to have µQ = 1011 N m−2 , like the secondary. Hence, τa = a/a˙ = 2 × 1010 yr. Just like in the case of tidal dissipation within the secondary, radiation effects appear to be about six orders of magnitude more important than tidal dissipation. Once again, we caution that the rigidity and the Q-factor of small bodies are basically unknown. However, since BYORP and YORP can accomplish everything that was previously attributed to tidal forces, we believe that tidal dissipation is indeed unimportant over the lifetime of NEA binaries. Close planetary encounters are the last non-radiationrelated process that could, in principle, restrict the secondary’s radial migration due to BYORP. Tidal forces of a planet could perturb the binary system by either “flipping” the synchronous lock or perturbing the mean motion sufficiently so that the secondary escapes the synchronous lock. In the first case, after the fly-by is over, the secondary will have a 50%–50% chance to end up with the nearside and the farside reversed. In this case the leading and trailing hemisphere flip, meaning that Tζ changes sign. If these reversals happen often enough, BYORP becomes a random-walk process, and the binary’s lifetime would be much longer. Quantitatively, in order for the secondary’s long axis to switch briefly from pointing toward the primary in order

428

´ M. Cuk, J.A. Burns / Icarus 176 (2005) 418–431

to face the planet during the fly-by, the planet’s tidal forces have to overwhelm the primary’s mp m (16)  3 , 3 a r where m is the planet’s mass and r  is the distance from the planet to the secondary. Equation (16) makes evident that the binary would have to pass within a few planetary radii for this to hold (more precisely r  /R   1.5a/Rp , with the factor of 1.5 arising because planetary mean densities (Earth and Venus) are much higher than those of most asteroids). If we demand that r  < 10R  for such a flip to happen, we have an event that is only 100 times more likely than a direct collision (i.e., r  < R  ). If a km-sized asteroid impacts Earth about every 105 years, and about 2 × 103 NEAs are larger than 1 km (Morrison et al., 2002), then, on average, any given large NEA would pass closer than ten Earth radii from our planet once in 2 million years. Presumably the probabilities are about the same for Venus (cf. Gladman et al., 1996), leaving a 106 yr timescale for this process, or about two orders of magnitude too infrequent to produce a random-walk behavior in BYORP. Another way that a planetary fly-by can disrupt synchronous lock is by changing the semimajor axis and, therefore, the secondary’s mean motion. The new mean motion will not be synchronized with the spin, and when synchronization is re-established, it is equally likely that the leading and trailing hemispheres will be exchanged. Lissauer (1986) explored the possibility of impacts breaking the synchronous lock, and derived the criterion for the change in the spin rate needed to rotate the body more than 90◦

3(B − A) 1/2 ω = , (17) Ω C where ω = Ω − n, so in our case ω = −n. Our Gaussian spheroids typically have moments of inertia of a few times ρR 5 , with A, B and C typically differing by about a single ρR 5 . Hence 3(B − A)/C is of order unity, so a neardisruption of the system is necessary to break the synchronous lock in this way. This would require the system’s Hill sphere (bounded by planetary tides) to shrink to about the size of binary’s separation, which once again brings us to Eq. (16). Therefore, this method of randomizing BYORP is no more effective than the direct planetary tides affecting the secondary, as both require close planetary approaches which do not happen more often than once in 106 yr. Note that we used the best-case scenario in computing the probabilities of planetary tidal perturbations, and that we do not imply that most binaries will get disrupted in a million years (see Bottke and Melosh (1996b) for a more detailed analysis of effects of planetary fly-bys on binary NEAs). Apparently, none of the above-mentioned processes are capable of significantly interfering with BYORP. However, we still have not explored the full implications of the orbital evolution of a system consisting of two irregularly shaped components. Apart from BYORP and YORP acting on the

secondary, radiation forces should also affect the primary. Also, the spin axes and the mutual orbit necessarily precess around each other, depending on the system’s exact parameters. The last phenomenon can tie different components of the system to each other, and thereby change the way that BYORP operates. While NEA primaries are often close to having generally spherical shapes (Margot et al., 2002), YORP usually depends on much smaller-scale irregularities of the surface, which are always omnipresent on low-gravity bodies. If we assume that the primary behaves like a Gaussian spheroid, the YORP behavior should just scale as 1/R 2 , and a 1-km primary with an initial spin period of 3 h should have a YORP timescale of about 5 × 105 yr. This timescale is somewhat longer than that for BYORP, but these results are valid only for the case in which the two angular momenta are decoupled and can evolve in parallel without affecting each other. Precession, a dynamical topic that has been studied for centuries longer than YORP, is quite complex. The relative magnitudes of the secondary’s orbital and primary’s rotational angular momenta determine which one effectively precesses around the other (strictly speaking, both vectors precess around their sum). In our model system, L  0.16Lp , so the primary’s angular momentum clearly dominates the system (we will assume αp = 0.4 throughout this section). But once the secondary’s momentum starts precessing around that of the primary, the two become closely connected. Perturbations of the satellite’s orbital plane can be transmitted to the primary. For example, if the secondary has a small inclination to the primary’s equator, but the latter is significantly inclined to the binary’s heliocentric orbital plane, the orbit’s inclination relevant for BYORP is not the free inclination to the Laplace plane (i.e., primary’s equator), but rather the primary’s obliquity. Similarly, the inclinationchanging BYORP torque will ultimately affect the primary’s obliquity. In the opposite case, where the primary’s obliquity is low but the secondary’s free inclination is high, BYORP will affect mostly the secondary’s free inclination (we can basically ignore such precession as long as it is much slower than the binary’s heliocentric motion). In intermediate cases, the evolutions of inclination and obliquity are necessarily more complex, with both quantities being substantially modified. Obviously, the more interesting possibility is that of a high-obliquity, low-inclination system (in our terminology, “inclination” is the angle between the secondary’s orbital plane and the primary’s equator, while “obliquity” is the angle between the primary’s equator and the pair’s heliocentric orbit plane; see Fig. 5). Systems of this character should be common, since both tidal and rotational disruption are likely to produce binaries with approximately aligned angular momenta (Richardson et al., 1998). For such systems, we need to consider that BYORP now modifies the system’s total angular momentum (including that of the primary). The primary’s rotational momentum for most of the pairs listed

Radiation effects on binary NEAs

Fig. 5. Schematic illustration of our definitions of inclination and obliquity of the binary system, as used in Section 5. Inclination (i) is the angle between the binary’s mutual orbit plane and the primary’s equator plane, while obliquity (ε) is the angle between the primary’s equator and the binary pair’s heliocentric orbit. We also assume that the secondary’s rotation axis is perpendicular to the binary’s mutual orbit plane.

in Table 1 is several times larger than the secondary’s orbital angular momentum. For our model system, the new timescale for the obliquity-changing component of BYORP is [Tη /(Lp + L)]  2.2 × 105 yr. This is of the same order of magnitude as our estimate of the timescale for the primary’s YORP-induced obliquity evolution. The implication of this result is that, if primary and secondary exhibit different YORP and BYORP behaviors, respectively (i.e., one being Type I and the other Type II), the two obliquitychanging torques can, in principle, cancel one another. Since a changes due to BYORP, the outward-evolving system can always find the distance at which the two torques will balance out, and it appears that such a separation is similar to the one we actually observe for binary NEAs. Also, if the system’s obliquity is around 50◦ –60◦ , a separation-changing torque can also disappear, and such a system would be stable. It would also, at least locally, be stable to perturbations, as both a˙ and the variation in obliquity ε˙ change signs when the other variable passes through its equilibrium value. Roughly speaking, we can expect just less than half of the freshly formed binary systems to consist of a primary with Type II YORP behavior and a secondary with Type I BYORP behavior, or vice versa. The systems where BYORP and the primary’s YORP belong to the same type are not viable, and such systems are likely to re-impact on the BYORP timescale. If any extant systems do inhabit this intermediateobliquity stable state, their inclinations should reveal their nature. If future observations find that many binary NEAs indeed have inclinations (relative to their heliocentric orbits) in the 50◦ –60◦ range, a balance between YORP and BYORP will be a likely cause. A stabilization of the orbit by balancing YORP and BYORP torques could also account for the small, but non-zero, eccentricity of many NEA binaries. If the system is close to that critical inclination for which the ζ -component of BYORP vanishes, eccentricity-damping should become much slower than we calculated in Section 4. Small perturbations from planets can periodically excite the eccentricity, which

429

would damp out too slowly to reach zero until the next perturbation. However, major perturbations of the orbit due to close fly-bys would also toss about the binary’s inclination, and BYORP would be “re-activated,” which would bring the eccentricity back to low values. In those cases for which i > ε, the system eventually evolves back toward a collision and the precession of the orbital plane around the equatorial plane can be ignored, unless that precession is rather fast. To first order, the precession of the nodes of a satellite’s orbit around an oblate body (Murray and Dermott, 1999) is

  3 Cp − (Ap + Bp )/2 Rp 2 n. Ω˙ n = − (18) 2 a mp Rp2 The expression in square brackets is commonly known as the primary’s J2 gravity moment. The line of nodes always precesses in the direction opposite to that of orbital motion. If we assume that our model binary is retrograde and has a heliocentric orbital period of 1 yr, then the 1:1 commensurability between the nodal regression rate with the Sun’s apparent mean motion would require J2 = 0.024, which is not exceptional for an irregularly shaped body. Actually, our model system is likely to precess in less than one year. In any case, it is clear that resonances between the Sun (and therefore the radiation-caused effects) and the binary’s line of nodes are possible. What would exactly happen in case of such a commensurability cannot be answered without additional numerical modeling that is outside this paper’s scope. In any case, such a rapid precession in principle can be detected by observations, especially those of eclipsing systems. For prograde orbits, there is another interesting, potentially stable state, somewhat similar to the one just discussed. It involves the resonance between the Sun’s apparent mean motion and the precession rate of the longitude of pericenter of the secondary’s orbit ( ˙ ). To first-order,  ˙ = −Ω˙ n , so the conditions for this commensurability are the same as for the nodal one (Eq. (18)), even if the physics is likely to be very different. This type of resonance is known as the “evection” resonance and it is thought to have been important during early history of Earth–Moon system (Touma and Wisdom, 1998; Touma, 2000). The Moon likely evolved through an evection resonance very early, and during its resonant phase the lunar semimajor axis was constant while its eccentricity increased. The resonance is broken once the eccentricity becomes too high for continued stability of the lock. In order for the evection resonance to be permanently stable, e˙ < 0 due to BYORP must balance e˙ > 0 due to the resonance. The strength of the evection perturbation on eccentricity (Brouwer and Clemence, 1961) is roughly 15 n2 de =− e sin(2λ − 2 ), (19) dt 4 n where λ and n are the mean longitude and apparent mean motion of the Sun. Using parameters of our “standard binary,” we get τe
102 yr, which is much shorter than the BYORP timescale. The only way that the evection resonance

430

´ M. Cuk, J.A. Burns / Icarus 176 (2005) 418–431

can be stable is if sin(2λ − 2 ) 1. This would require the pericenter to be either always parallel or always normal to the direction to the Sun. This can help us to decide if an evection resonance is likely in a system that has been observed only for a short time, so no measurement of precession was possible. One binary, 2000 DP107, described by Margot et al. (2002), is especially interesting when considering evection. This system appears to have a low inclination of its mutual orbit to its heliocentric orbit (i = 15 ± 7◦ ); in this case, a J2 of about 0.05 would suffice to produce the evection resonance. Observations were made during a close Earth encounter in 2000 between September 30th and October 7th (λ  185◦ –195◦ ), when Ω = 10±17◦ and ω = 7±39◦ . The error bars are large, and certainly allow that  = Ω + ω  10◦ , which would mean that pericenter is anti-aligned with the Sun. Of course, we cannot exclude many other orientations in view of the large error bars, but the evection resonance is, at least in principle, possible for this binary system. To conclude this chapter, we note that our crude analysis of various processes affecting long-term evolution of the binary NEA systems has identified at least three potentially stable states, the YORP-BYORP cancellation point, which is stable within our model, and the Sun-node and Sun-apse resonances, the effects of which need to be explored in more detail in future work.

(5)

(6)

(7)

6. Conclusions Since the effect of radiation on binary asteroids has not been discussed before, this paper contains several new, if tentative, findings: (1) We describe a new variety of YORP, which we call BYORP (Binary YORP), that should affect the orbits of synchronous asteroidal satellites. (2) We present results of a numerical experiment based on ˇ methods described by Vokrouhlický and Capek (2002), using random Gaussian spheroids. The experiment confirms the existence of BYORP, and suggests that the orbits of typical NEA systems should evolve on 104 to 105 -yr timescales due to BYORP. In terms of geometry, the effects of BYORP are very similar to those of YORP, and our test-systems exhibit types of behavior that are similar to those found by Vokrouhlický and ˇ Capek (2002) for YORP. (3) We also find that BYORP should affect the binary’s orbital eccentricity, damping it on the same 104 –105 -yr timescale. In addition to that, YORP acting on the secondary should also evolve its spin on similar 104 –105 -yr timescales, making it a likely mechanism for synchronization of the secondary’s spin. (4) We estimate the timescale for competing tidal processes and find them to be much longer. We therefore suggest that both the low eccentricities and the synchronous ro-

(8)

tation of many NEA binary pairs almost certainly result from radiation forces, rather than tidal dissipation. We discount several explanations for the apparent survival of binaries beyond the brief BYORP timescale, all involving processes other than radiation. BYORP is to be expected, since it is simply another manifestation of YORP which has arguably been proven to exist within the Koronis family (Vokrouhlický et al., 2003). Binary near-Earth asteroids do appear, by and large, to have synchronous secondaries, and close encounters with Earth and Venus are too rare to disrupt that configurations in times timescales shorter than that for BYORP. Finally, orbit expansion due to tidal dissipation appears to be many orders of magnitude slower than that caused by BYORP. As a solution to the binary lifetime problem, we suggest that many extant systems are in one of several possible stable states, in which the total radiative force is zero. One such state requires an equatorial orbit of the secondary, an obliquity of 50◦ –60◦ for the primary, and separations of several primary’s radii, the last being consistent with observations of known binary NEAs. Another putative stable point might involve a resonance between the precession rate of the line of nodes of a retrograde orbit and the Sun’s apparent motion, both of which should have periods of the order of a year for many of the known binary NEAs. A similar resonance involving the longitude of pericenter (“evection” resonance) could be a stable point for prograde binaries. All the possible stable states described above could be directly verified by observing the binaries’ inclinations and the orientation of the lines of nodes and apsides. It is also possible that other, unidentified, stable states exist, so data on the orientation of a number of NEA binaries’ orbits would be invaluable for theorists trying to model their long-term evolution.

Acknowledgments The authors thank Petr Pravec for allowing use of his statistical data on NEA binaries (Table 1), David Rubincam for permitting us to reproduce Fig. 1 from his paper, and Jean-Luc Margot for very helpful discussions. Joe Spitale and David Rubincam provided helpful reviews.

References Bottke, W.F., Melosh, H.J., 1996a. The formation of asteroid satellites and doublet craters by planetary tidal forces. Nature 381, 51–53. Bottke, W.F., Melosh, H.J., 1996b. Binary asteroids and formation of doublet craters. Icarus 124, 372–391. Bottke, W.F., Rubincam, D.P., Burns, J.A., 2000. Dynamical evolution of main belt meteoroids: Numerical simulations incorporating planetary perturbations and Yarkovsky thermal forces. Icarus 145, 301–331.

Radiation effects on binary NEAs

Bottke, W.F., Vokrouhlický, D., Brož, M., Nesvorný, D., Morbidelli, A., 2001. Dynamical spreading of asteroid families by the Yarkovsky effect. Science 294, 1693–1696. Bottke, W.F., Vokrouhlický, D., Rubincam, D.P., Brož, M., 2002. The effect of Yarkovsky thermal forces on the dynamical evolution of asteroids and meteoroids. In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R.P. (Eds.), Asteroids III. Univ. of Arizona Press, Tucson, AZ, pp. 395–408. Brouwer, D., Clemence, G.M., 1961. Methods of Celestial Mechanics. Academic Press, New York. Burns, J.A., 1976. Elementary derivation of the perturbation equations of celestial mechanics. Am. J. Phys. 44, 944–949. Burns, J.A., 2002. Two bodies are better than one. Science 297, 942–943. ´ Burns, J.A., Cuk, M., 2002. Death to satellites of extrasolar planets: Tidal collapse, vaporization and Yarkovsky decay. Bull. Am. Astron. Soc. 34, 916. Burns, J.A., Lamy, P.L., Soter, S., 1979. Radiation forces on small particles in the Solar System. Icarus 40, 1–48. Carruba, V., Burns, J.A., Bottke, W.F., Nesvorný, D., 2003. Orbital evolution of the Gefion and Adeona asteroid families: Close encounters with massive asteroids and the Yarkovsky effect. Icarus 162, 308–327. ´ Cuk, M., Burns, J.A., Carruba, V., Jacobson, R.A., Nicholson, P.D., 2002. New Yarkovsky effect for irregular satellites. Bull. Am. Astron. Soc. 34, 882. Farinella, P., Vokrouhlický, D., Hartmann, W.K., 1998. Meteorite delivery via Yarkovsky orbital drift. Icarus 132, 378–387. Erratum: Icarus 134 (1998) 347. Giorgini, J.D., 13 colleagues, 2002. Asteroid 1950 DA’s encounter with Earth in 2880: Physical limits of collision probability prediction. Science 296, 132–136. Gladman, B.J., Burns, J.A., Duncan, M., Lee, P., Levison, H.F., 1996. The exchange of impact ejecta between terrestrial planets. Science 271, 1387–1392. Goldstein, H., Poole, C., Safko, J., 2002. Classical Mechanics, third ed. Addison–Wesley, San Francisco, CA. Hamilton, D.P., Burns, J.A., 1991. Orbital stability zones about asteroids. Icarus 92, 118–131. Lissauer, J.J., 1986. Can cometary bombardment disrupt synchronous rotation of planetary satellites? J. Geophys. Res. 90, 11,289–11,293. Margot, J.-L., 2002. Worlds of mutual motion. Nature 416, 694–695. Margot, J.-L., Nolan, M.C., Benner, L.A.M., Ostro, S.J., Jurgens, R.F., Giorgini, J.D., Slade, M.A., Campbell, D.B., 2002. Binary asteroids in the near-Earth object population. Science 296, 1445–1448. Merline, W.J., Weidenschilling, S.J., Durda, D.D., Margot, J.-L., Pravec, P., Storrs, A.D., 2002. Asteroids do have satellites. In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R.P. (Eds.), Asteroids III. Univ. of Arizona Press, Tucson, AZ, pp. 289–312. Milani, A., Nobili, A.M., Farinella, P., 1988. Non-Gravitational Perturbations and Satellite Geodesy. Hilger, Bristol. Morrison, D., Harris, A.W., Sommer, G., Chapman, C.R., Carusi, A., 2002. Dealing with the impact hazard. In: Bottke, W.F., Cellino, A., Paolicchi,

431

P., Binzel, R.P. (Eds.), Asteroids III. Univ. of Arizona Press, Tucson, AZ, pp. 739–754. Muionen, K., 1996. Light scattering by Gaussian random particles. Earth Moon Planets 72, 339–342. Muionen, K., Lagerros, J.S.V., 1998. Inversion of shape statistics for small Solar System bodies. Astron. Astrophys. 333, 753–761. Murray, C.D., Dermott, S.F., 1999. Solar System Dynamics. Cambridge Univ. Press, Cambridge, UK. Namouni, F., 1999. Secular interactions of coorbiting objects. Icarus 137, 293–314. Peale, S.J., 1977. Rotation histories of the natural satellites. In: Burns, J.A. (Ed.), Planetary Satellites. Univ. of Arizona Press, Tucson, AZ, pp. 87– 112. Pravec, P., Harris, A.W., 2000. Fast and slow rotation of asteroids. Icarus 148, 12–20. Richardson, D.C., Bottke, W.F., Love, S.G., 1998. Tidal distortion and disruption of Earth-crossing asteroids. Icarus 134, 47–76. Richardson, D.C., Leinhardt, Z.M., Melosh, H.J., Bottke, W.F., Asphaug, E., 2002. Gravitational aggregates: Evidence and evolution. In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R.P. (Eds.), Asteroids III. Univ. of Arizona Press, Tucson, AZ, pp. 501–515. Rubincam, D.P., 1987. LAGEOS orbit decay due to infrared radiation from Earth. J. Geophys. Res. 92, 1287–1294. Rubincam, D.P., 1988. Yarkovsky thermal drag on LAGEOS. J. Geophys. Res. 93, 13,805–13,810. Rubincam, D.P., 2000. Radiative spin-up and spin-down of small asteroids. Icarus 148, 2–11. Slivan, S.M., 2002. Spin vector alignment of Koronis family asteroids. Nature 419, 49–51. Slivan, S.M., Binzel, R.P., Crespo da Silva, L.D., Kaasalainen, M., Lyndaker, M.M., Krˇco, M., 2003. Spin vectors in the Koronis family: Comprehensive results from two independent analyses of 213 rotation lightcurves. Icarus 162, 285–307. Spitale, J.N., 2001. Detailed study of the Yarkovsky effect on asteroids and Solar System Implications. Ph.D. thesis. Univ. of Arizona, Tucson, AZ. Touma, J., 2000. The phase space adventure of the Earth and the Moon. In: Canup, R.M., Righter, K. (Eds.), Origin of the Earth and the Moon. Univ. of Arizona Press, Tucson, AZ, pp. 165–178. Touma, J., Wisdom, J., 1998. Resonances in the early evolution of the Earth–Moon system. Astron. J. 115, 1653–1663. ˇ Vokrouhlický, D., Capek, D., 2002. YORP-induced long-term evolution of the spin state of small asteroids and meteoroids: Rubincam’s approximation. Icarus 159, 449–467. Vokrouhlický, D., Nesvorný, D., Bottke, W.F., 2003. The vector alignments of asteroid spin vectors by thermal torques. Nature 425, 147–151. Weidenschilling, S.J., Paolicchi, P., Zappalà, V., 1989. Do asteroids have satellites? In: Binzel, R.P., Gehrels, T., Matthews, M.S. (Eds.), Asteroids II. Univ. of Arizona Press, Tucson, AZ, pp. 643–658. Yoder, C.F., 1982. Tidal rigidity of Phobos. Icarus 49, 327–346.