Long-Term Evolution of Objects in the Kuiper Belt Zone—Effects of Insolation and Radiogenic Heating

Icarus 160, 300–312 (2002) doi:10.1006/icar.2002.6976

Long-Term Evolution of Objects in the Kuiper Belt Zone—Effects of Insolation and Radiogenic Heating Young-Jun Choi, Merav Cohen, Rainer Merk, and Dina Prialnik Department of Geophysics and Planetary Sciences, Tel Aviv University, Ramat Aviv 69978, Israel E-mail: [email protected] Received November 11, 2001; revised July 21, 2002

The Kuiper Belt zone is unique insofar as the major heat sources of objects a few tens of kilometers in size—solar radiation on the one hand and radioactive decay on the other—have comparable power. This leads to unique evolutionary patterns, with heat waves propagating inward from the irradiated surface and outward from the radioactively heated interior. A major radioactive source that is considered in this study is 26 Al. The long-term evolution of several models with characteristics typical of Kuiper Belt objects is followed by means of a 1-D numerical code that solves the heat and mass balance equations on a spherically symmetric grid. The free parameters considered are radius (10–500 km), heliocentric distance (30–120 AU), and initial 26 Al content (0–5 × 10−8 by mass). The initial composition assumed is a porous mixture of ices (H2 O, CO, and CO2 ) and dust. Gases released in the interior are allowed to escape to the surface. It is shown that, depending on parameters, the interior may reach quite high temperatures (up to 180 K). The models suggest that Kuiper Belt objects are likely to lose the ices of very volatile species during early evolution; ices of less volatile species are retained in a surface layer, about 1 km thick. The models indicate that the amorphous ice crystallizes in the interior, and hence some objects may also lose part of the volatiles trapped in amorphous ice. Generally, the outer layers are far less affected than the inner part, resulting in a stratified composition and altered porosity distribution. These changes in structure and composition should have significant consequences for the short-period comets, which are believed to be descendants of Kuiper Belt objects. c 2002 Elsevier Science (USA) Key Words: comets; comets, composition; ices; Kuiper Belt objects.

1. INTRODUCTION

The region of the Solar System beyond Neptune’s orbit, known as the Kuiper (or Edgeworth–Kuiper) Belt, is populated by small bodies, which represent the surviving remnants of the once more massive protoplanetary disk of planetesimals. Having been formed farthest out from the Sun, they are believed to constitute the most primitive, least thermally processed matter. Since 1992, when the first object in this class was discovered, over 560 more have been detected, with semi-major axes between 35 and over 100 AU. Recently, evidence of H2 O was discovered in the

spectrum of a Kuiper Belt object (KBO), indicating a kinship to comets (Brown et al. 1999). In fact, it has long been suggested by Duncan et al. (1988) and others that the Kuiper Belt is the source of short-period comets and of Centaurs (small bodies characterized by dynamically unstable orbits between Jupiter and Neptune). The latter are probably transition objects from the Kuiper Belt reservoir to the inner Solar System (Gladman and Duncan 1990). Their nature—whether asteroidal or cometary— is uncertain; they show evidence of complex organic molecules, and water ice was detected on the surfaces of several such objects (Brown and Koresko 1998, Brown et al. 1998, Luu et al. 2000). Observed Kuiper Belt objects are larger than typical comets, their radii ranging between a few tens and a few hundred of kilometers (Jewitt et al. 1998, Gladman et al. 1998, Chiang and Brown 1999, Sheppard et al. 2000, Trujillo et al. 2001). The largest objects have radii of about 500 km (e.g., Varuna, 450 km, according to Jewitt et al. 2001). Whether KBOs are pristine, or else, have been altered either thermally or in terms of structure, is still an open question. Their sizes may have been significantly altered by collisions (Davis and Farinella 1997), but the low collision velocities in the Kuiper Belt could not significantly affect their thermal structure (Orosei et al. 2001). Thus the main heat sources that could affect them are insolation and radioactive decay, which determine the evolution of comets. In the region of the inner Solar System, the Sun’s radiation is the dominant heat source, and the outer layers of a comet nucleus are those most affected by it. Far from the Sun, where comets spend most of their lifetimes, internal heating by radionuclides may play an important role. Besides 40 K, 235 U, 238 U, and 232 Th, which constitute the major heat sources in terrestrial planets, the radioactive isotope 26 Al has long been recognized as a potential heat source capable of melting bodies with radii between 100 and 1000 km (Urey 1955). Its lifetime (half-life of 7.4 × 105 years) is just long enough for it to outlive the formation of comets and small icy satellites, and short enough to generate sufficient heat power for overcoming cooling by conduction. Although radioactive heating is expected to have affected comets to a much lesser extent than satellites and planets, it was found to cause an internal temperature rise of tens of kelvins (e.g., Whipple and Stefanic 1966, Wallis 1980,

300 0019-1035/02 $35.00 c 2002 Elsevier Science (USA)  All rights reserved.

301

THERMAL EVOLUTION IN THE KUIPER BELT

TABLE I Radioactive Isotope Characteristics Isotope 40 K 232 Th 238 U 235 U 26 Al

centric distance dH , radius R, and the abundance of radioactive species X rad

τ (yr)

X0

H (erg/g)

Q 0 (erg/g/yr)

1.82(9) 2.0(10) 6.50(9) 1.03(9) 1.06(6)

1.1(−6) 5.5(−8) 2.2(−8) 6.3(−9) 6.7(−7)

1.72(16) 1.65(17) 1.92(17) 1.86(17) 1.48(17)

10.7 0.454 0.645 1.13 9.4(4)

4π 3 X rad Hrad 1 R2 , (1 − A)L  2 = R ρ 4 3 τrad dH

(1)

where A is the albedo, L  is the solar luminosity, τrad is the characteristic radioactive decay time, and Hrad is the energy released per unit mass. If we consider only 26 Al, which is the most powerful source, adopt A = 0.04, ρ = 0.7 g cm−3 , and normalize the parameters,

Consolmagno and Lewis 1978, Schubert et al. 1981, Elsworth and Schubert 1983, Prialnik et al. 1987, Prialnik and Bar-Nun 1990, Haruyama et al. 1993, Yabushita 1993, Prialnik and Podolak 1995, and Podolak and Prialnik 1997). Moreover, when very low thermal conductivities are assumed, a runaway increase in the internal temperature occurs and most of the water ice crystallizes (if it is initially amorphous), although for a sufficiently large conductivity, the amorphous ice may be almost completely preserved. Typical abundances—mass fractions of the dust component—of radioactive isotopes in chondritic meteorites (MacDonald 1959), and other relevant characteristics, are given in Table I. It is instructive to compare the power supplied by these two major energy sources, in terms of the leading parameters: helio-

X rad,0 = 6.7 × 10−7

x = X rad / X rad,0 y = dH /dH,0

dH,0 = 30 AU

z = R/R0

R0 = 10 km,

we obtain the following relation in three variables x y 2 z = 0.55,

(2)

which is illustrated by the shaded surface in Fig. 1. Above this surface radioactive heating dominates, while solar heating dominates below it. The interesting conclusion of this exercise is that the Kuiper Belt zone is quite unique within the Solar System in

x y2z=0.55

50 45 40 35

z

30 25 20 15 10 1

5 0 0

2 −0.5

3 −1

4 −1.5

−2

5

y

log (x) 10

FIG. 1. space.

Solution of Eq. (2), showing the plane where solar heating and radiogenic energy release have equal power, in the distance—radius—26 Al-content

302

CHOI ET AL.

the sense that the major heat sources—solar radiation on the one hand and radioactive decay on the other—have comparable power for objects of the observed KBO size. This should lead to unique evolutionary patterns. When radioactive heating dominates (e.g., in the Oort cloud), evolution (or change in structure and composition) proceeds from the center outward. By contrast, in the inner Solar System, where insolation dominates, structural changes propagate from the surface inward. Kuiper Belt objects may experience both, and the resulting structure may be quite unexpected. In this context, we should keep in mind that although the two energy sources are comparable in power, they differ in character: radioactivity is a body source, homogeneously distributed, and time-dependent; solar radiation is a surface source and, for KBOs which have close to circular orbits, it is almost constant in time, although its effect should decrease as the surface temperature reaches equilibrium. The thermal evolution of distant comets (dH = 100 AU) under the effect of radiogenic heating, with special emphasis on 26 Al, was studied by Prialnik and Podolak (1999) for a range of radii and porosities, but only for relatively short time spans—until the rising central temperature attained maximum. Only water ice was considered in that study, and volatiles trapped in the amorphous ice. Changes in structure and composition thus resulted only when the temperatures were sufficiently high for crystallization and gas release to occur. In those cases a layered structure emerged, with ices of more volatile species overlying ices of less volatile species. Recently, the thermal evolution of KBOs was considered by De Sanctis et al. (2001), who evolved models of two such objects for up to 10 Myr, including the long-lived radioactive species and, in one case, 26 Al as well, and assuming a composition of mixed ices and dust. They found that the surface layers became depleted of the most volatile species. However, the evolution was not followed long enough for the effect of radiogenic heating to become fully developed and mainly the effect of insolation was, in fact, examined. In the present study we follow, numerically, the long-term evolution of KBO models over the [R, dH , X 0 (26 Al)] parameter space. The 26 Al issue is further discussed in the next section. The model, assumptions, and parameters are addressed in Section 3; the results of numerical computations are described and discussed in Section 4, focusing on volatile depletion and loss of homogeneity in structure as well as in composition. Finally, our main conclusions are summarized and compared to other studies in Section 5. 2.

26

AL AND THE FORMATION OF SMALL BODIES IN THE SOLAR SYSTEM

The radionuclide 26 Al is already widely used in thermal modeling of asteroidal bodies (Miyamoto et al. 1981, Akridge et al. 1998, Ghosh and McSween 1998). For silicate small bodies such as asteroids, accretion models predict growth times in the range between 104 years and 1 Myr, depending on whether or not conditions allow for runaway accretion (Weidenschilling 1988,

Wetherill and Stewart 1989). Thus, the comparably short lifetime of 26 Al and the growth time of asteroidal bodies are not mutually contradictory (Ghosh and McSween 1998). Furthermore, the decay product 26 Mg can be found experimentally in Ca–Al inclusions of meteorites, the best-known among them being Allende (Lee et al. 1976). Srinivasan et al. (1999) detected for the first time 26 Mg in a differentiated meteorite and thus could confirm the role of 26 Al in the differentiation of meteoritic parent bodies. Renewed interest in this radionuclide followed the detection of interstellar 1.809 Mev γ -rays from the decay of 26 Al (Mahoney et al. 1984, Share et al. 1985, Clayton and Leising 1987). Accretion models for KBOs (Kenyon and Luu 1998, Farinella et al. 2000) or for cometary icy bodies (Weidenschilling 1997) show that accretion times are in the range between 105 years (Weidenschilling 1997) and several million years (Kenyon and Luu 1998, Farinella et al. 2000). An upper limit for accretion time scales in the Kuiper Belt region seems to be the formation time of Neptune, since it is assumed that the formation of Neptune efficiently terminated growth in the Kuiper Belt region (Farinella et al. 2000). Kenyon and Luu (1998) assumed a value of 50–100 Myr (Pollack et al. 1996) as an upper limit of this time scale. Models based on Wetherill’s particle-in-abox method (Wetherill and Stewart 1989), such as that used by Kenyon and Luu (1998), lead to longer accretion times than models that include the transition between drag-dominated and gravity-dominated regimes (Weidenschilling 1997). The former models start with a seed body and exclude sticking effects, whereas the latter models describe the growth of particles from micrometer to kilometer size. Weidenschilling (1997) showed that it is possible to grow large icy bodies of a radius of about 40 km within 2.5 × 105 years in the region at 30 AU. These investigations show that the short-lived isotope 26 Al can be considered an important contributor to internal heat generation within KBOs. All the observational evidence points toward an interstellar isotopic ratio of 26 Al/ 27Al ≈ 5 × 10−5 , implying an initial mass fraction X 0 (26 Al) ≈ 7 × 10−7 in dust (rock) and presumably an order of magnitude less in objects whose time of aggregation did not exceed a few million years. Therefore, an initial 26 Al mass fraction of a few 10−8 should be a reasonable working assumption. 3. MODEL, ASSUMPTIONS, AND PARAMETERS

We use a numerical code that solves the equations of mass and energy conservation for a spherically symmetric porous nucleus. The nucleus is assumed to be composed of dust and a mixture of volatiles, which may be found in a solid or gaseous state. The dust is assumed to include radioactive elements in abundances typical of meteorites. The water ice is assumed to be initially amorphous, but crystallization is taken into account according to the temperature-dependent rate λ(T ) = 1.05 × 1013 exp (−5370/T ) s−1 , given by Schmitt et al. (1989). In some test

303

THERMAL EVOLUTION IN THE KUIPER BELT

cases we also consider gases that are trapped in the amorphous ice and released upon crystallization. The equations of mass conservation for a volatile species α in gaseous (subscript g) and solid (subscript s) phases are, respectively, ∂ρα,g + ∇ · Jα = q α , ∂t ∂ρα,s = −qα , ∂t

(3) (4)

where ρ, with the appropriate subscript, denotes the bulk density (mass per unit volume of cometary material) of a species in a given phase, J is the gas flux (for details, see Prialnik et al. 1993), and q is the sublimation rate. The porosity is given by p =1−

 α

ρα,s / α − ρdust / dust ,

(5)

where α is the characteristic solid density. The energy conservation equation is     ∂  ρα u α + ∇ · F + u α Jα = λρa Hac + Q˙ − qα Hα , ∂t α α α (6) subject to the following relations  qα = S(Pα (T ) − Pα ) Q˙ =



µα 2π Rg T

X 0,j Hj τj−1 exp(−t/τj ),

(7) (8)

j

where u denotes specific heat, H is the latent heat, P is saturated vapor pressure as given by the Clausius–Clapeyron equation, P is gas pressure, µ is molecular mass, Rg is the gas constant, S is the surface-to-volume ratio (depending on porosity and pore size), Q˙ is the rate of radioactive energy release, including all the isotopes of Table I, and τj is the decay time of the j  th radioactive isotope; X 0,j is its mass fraction, which for the long-lived radionuclides is obtained by the values listed in Table I, multiplied by the dust mass fraction, while for 26 Al it is a free parameter. Finally, Hj is the energy released per unit mass upon decay. The heat flux is F = −ψ( p)K (T )∇T , where K (T ) is the thermal conductivity of the solid material (as given by Klinger 1980 for the ice component, while for dust K = 106 erg s−1 cm−1 K−1 ), corrected by a factor ψ( p) = 1 − p 2/3 < 1 to take account of the porosity. The boundary conditions are F(0, t) = Jα (0, t) = 0 Pα (R, t) = 0 F(R, t) = σ T 4 (R, t) − (1 − A)L  /16πdH2 (t).

The system of non-linear partial differential equations is solved numerically on a one-dimensional spherical grid. The radial resolution is higher near the surface, where gradients are steeper. The time steps are adjusted by the code to keep the changes in temperature during one time step confined to a few percent. Further details on the modeling of the porous medium, and on the method of computing the sublimation rate, may be found in Mekler et al. (1990); details of the implicit, iterative numerical code are given by Prialnik (1992). The working assumption of the present calculations is that gases released in the interior can easily escape, that is, can channel easily to the surface and out of the nucleus. This is justified by the very low strength of cometary material, on the order of 104 dyn cm−2 (Klinger et al. 1989, Kochan et al. 1989, Greenberg et al. 1995, Sirono and Greenberg 2000). For illustration consider an element of porous solid material containing a fraction X of volatile ice, which absorbs heat continually. Eventually, the ice will start sublimating and if the vapor is not allowed to escape, the resulting gas pressure in the pores may rise up to Rg XρT / pµ (the temperature will rise as well, controlled by the saturated vapor pressure, as in a pressure-cooker). This gas pressure exceeds the estimated material strength of cometary nuclei already for X ≈ 5 × 10−5 . Thus even a moderate buildup of internal pressure is bound to open channels that would release the pressure, allowing the gas to escape. As a result of the easy escape of gas, the gas pressure in the porous interior will be lower than the saturated pressure and therefore most of the heat released will be absorbed in sublimation. The effect of a few such scattered wide channels cannot, however, be accounted for by a 1-D code. We could, of course, adopt a very large average channel width, but this would also imply an unrealistic surface-to-volume ratio, as well as an unrealistically small Knudsen number, and it would amplify heat advection. Another approach is to assume that the effective permeability is sufficiently high to permit a quasi-steady-state approximation, where the first term on the left-hand side of the mass balance equation for the gas phases (Eq. (3)) is neglected, while at the same time adopt a reasonable average pore-size. Thus Eqs. (3)– (4) are replaced by ∇ · Jα = qα ,

∂ρα,s = −qα . ∂t

(9)

In this way we have to strictly solve only one time-dependent equation—although gas densities and production rates change as the temperature distribution changes—supplemented by structure (space-dependent) equations. This constitutes a huge computational advantage, since a long-term calculation including a detailed account of gas flow through the porous medium, coupled with heat transfer, would require a prohibitively large amount of computing time. Combining Eqs. (9) and integrating over volume, we obtain ˙ α,s = Jα (R, t)4π R 2 , −M

(10)

304

CHOI ET AL.

TABLE II Composition Parameters (Given ρ = 0.7 g cm−3 )

Component

Mass fraction

Specific density g cm−3

Partial density g cm−3

Relative volume

H2 O CO CO2 Dust Pores

0.40 0.05 0.05 0.50 —

0.917 1.1 1.56 3.25 —

0.280 0.035 0.035 0.350 —

0.305 0.032 0.022 0.108 0.533

which means that the total mass of gas ejected through the comet’s surface per unit time is equal to the total amount of gas evaporated throughout the nucleus per unit time, for each species. This approximation is valid for non-abundant species, for which the bulk density is low. It breaks down when the net gas sublimation rate is negative, that is, when recondensation surpasses sublimation, but this is not the case for the conditions considered in this paper. The calculations are terminated if vigorous water ice sublimation takes place. It is worth noting that a different approximation with the same computational advantage—reduction of the number of timedependent equations—has been used in other studies (e.g., De Sanctis et al. 2001). It assumes that when both the ice and gas phases are present, the gas density is equal to the saturated vapor density, which is a function of temperature. Strictly, this would imply that no evaporation/condensation could take place. However, as the temperature changes, the saturated density (pressure) changes with it, and this change can be translated into a rate of sublimation/condensation. This approximation is not valid for minor volatile components and fails close to the surface of the nucleus, but it applies to water ice in the interior of the comet. The two simplifying approximations are thus complementary. The free parameters of our study are the cometary radius, the heliocentric distance (≥30 AU), and the initial content of 26 Al. For all models we assume an initial temperature of 10 K, equal mass fractions of ices and dust, and a bulk density of 0.7 g cm−3 (corresponding to a porosity of 0.5). The ice mixture includes 10% (by mass) of each CO—representative of extremely volatile species—and CO2 —representative of moderately volatile species—and the rest is amorphous water ice. Details of composition-related properties are summarized in Table II. We assume an albedo of 0.04 and an emissivity of 0.96, but the exact albedo value has little effect on the results, so long as it is significantly smaller than unity, which is the case for KBOs (Jewitt et al. 2001 and references therein). 4. RESULTS OF EVOLUTIONARY CALCULATIONS

The eight models that were calculated, for different combinations of the three basic parameters, are listed in Table III. The baseline model is labeled a; it has a radius of 100 km—the typical value for radii of observed KBOs. A circular orbit is as-

sumed, at a heliocentric distance of 30 AU, although KBOs do not move in circular orbits; the exact orbit has no importance in a long-term evolution calculation focused on the interior of the nucleus, and the distance chosen is smaller than the average semi-major axis of ∼40 AU, since the most marked effect of irradiation is expected to occur near perihelion (i.e., closer than 40 AU). The initial 26 Al content is 5 × 10−8 , about 10 times lower than the nominal value (see Section 2), allowing for a few Myr formation time. The models in Table III are grouped according to different criteria: those in the first group differ from the baseline model each in only one parameter. This allows an assessment of the influence of each parameter. Models in the second group share one parameter with the baseline model, while the other two are determined by requirement (1), that is, equal power of the two energy sources. The parameter combinations are selfconsistent: thus, model e has the same 26 Al content, but a smaller radius than a, implying a shorter formation time at the same distance, which is compensated by a larger heliocentric distance, where formation rates are lower. Model f has the same radius as model a, but a larger distance and thus a lower 26 Al content; model g has a larger radius at the same distance as model a, compensated by a lower 26 Al content. The two models in the third group are meant to provide a lower limit to the effect of long-term evolution in the Kuiper Belt region. In what follows we shall refrain from describing in detail the evolution of each model; instead, we shall focus on the broad view of their evolutionary pattern, looking, in particular, for properties and processes that appear to be shared by objects in the Kuiper Belt. We shall compare models, in order to single out the potential effects of individual parameters. The evolution of the central temperature Tc , the surface temperature Ts , and the maximal temperature Tmax attained within the nucleus is shown in Fig. 2 for models a–d. The abscissa in these (and the following) diagrams is time on a logarithmic scale, since changes are largely dictated by radioactive decay, which is exponential. Clearly, 26 Al plays a crucial role, raising the internal temperature to ∼180 K. (Calculations were terminated at this point, when significant sublimation of water in the interior TABLE III Initial Parameters of Computed Models Model

R (km)

dH (AU)

X 0 (26 Al)

a b c d

100 100 100 10

30 30 120 30

5(−8) 0 5(−8) 5(−8)

a e f g

100 10 100 500

30 90 120 30

5(−8) 5(−8) 3(−9) 1(−8)

b h

100 10

30 30

0 0

THERMAL EVOLUTION IN THE KUIPER BELT

305

FIG. 2. Evolution of central, surface, and peak temperatures, Tc (dotted line), Ts (dashed line), and Tmax (solid line) for models a, b, c, d, respectively, as marked. Note that the time scale is logarithmic, since changes are more pronounced during the early stages of evolution. The temporary halt in the rising trend of the central temperature, first at ∼20 K and then again at ∼70 K is due to the sublimation of CO and CO2 , respectively, which absorbs the energy released by radioactive decay. The plateau at ∼140 K is maintained by the crystallization of amorphous ice, which is exoergic.

was in conflict with the assumption of easy escape of gases—a very large amount of vapor could not escape through very cold regions without considerable refreezing (see Section 3).) We note the effect of radioactive 40 K causing a late rise in temperature in model b. We also note the effect of a small radius: in model d the maximal temperature is only ∼140 K. In all cases, the temperature rises rapidly at the surface, due to insolation, and more slowly in the interior. When the surface temperature reaches equilibrium with the solar radiation, it stops rising, and the more slowly rising central temperature becomes the temperature maximum. In models a and c (large radii and 26 Al content), the temperature peak shifts again outward after the 26 Al decay. This is due to another energy source that has now been activated—the exoergic crystallization of amorphous ice. Thus

the temperature peak moves with the crystallization front from the center toward the surface. In all models there is a conspicuous temporary halt in the rising trend of the central temperature, first at ∼20 K and then again at ∼70 K. This is due to the sublimation of CO and CO2 , respectively, which absorbs the energy released by radioactive decay. Similar conclusions are drawn from Fig. 3, which shows the evolution of Tc , Ts , and Tmax for models a and e–g. In addition, we note the effect of the initial amount of 26 Al: lowering the mass fraction by an additional factor of ∼10 (to 3 × 10−9 —model f) reduces its effect significantly—the maximal internal temperature attained is only slightly above 50 K. Even for a large body, a reduction in the 26 Al content (model g) results in considerably lower temperatures, somewhat below 130 K. In this case,

306

CHOI ET AL.

FIG. 3.

Same as Fig. 2 for models a, e, f, g.

however, the effect of the next potent radionuclide, 40 K, is apparent, causing the internal temperature to rise above 130 K on a time scale of a few 108 years. We note that the evolutionary course is little affected by heliocentric distance, except that at distances much larger than 30 AU (models c, e, f), the surface temperatures rise above equilibrium for a while, but remain well below 50 K at all times. The effect of heliocentric distance is more pronounced for a smaller body: this emerges from the comparison of model e with model d of the previous group—the peak temperature is lower in the more distant object and it is maintained for a shorter time span. The reason is that while both models attain sufficiently high internal temperatures (130–140 K) for crystallization to occur, the small temperature difference between them results in a difference of a factor of 20 in the exponential crystallization rate. Keeping in mind that crystallization is exoergic, this explains why the process—as well as the relatively high temperature accompanying it—is main-

tained for a longer period of time in the slightly hotter object. Finally, when no 26 Al is present, the other radioactive isotopes have a limited effect on large bodies, as already mentioned, but no effect at all on small bodies (∼10 km), as shown in Fig. 4 for model h (cf. Consolmagno and Lewis 1978). The high internal temperatures that we have found to develop in objects in the Kuiper Belt region are bound to lead to alterations in structure and composition. We proceed to illustrate the thermal and structural evolution of the models in color-coded 3-D plots. The ordinate represents the depth into the nucleus, again on a logarithmic scale, since usually gradients become steeper toward the surface. Figure 5, top and bottom, shows the evolution of temperature for models a–d and a ∪ e–g, respectively. We draw attention to model g, the 500-km object, which appears to maintain a high temperature in the interior, possibly up to the present. Such a structure might have implications for unusual activity patterns, such as are sometimes observed

307

THERMAL EVOLUTION IN THE KUIPER BELT

FIG. 4.

Same as Fig. 2 for models b, h.

(Hainaut et al. 2000, Delsanti et al. 2001), if the relatively hot interior could be tapped (e.g., by a collision). Regarding the composition of the comet models, the changes are considerable, which should be expected in view of the high temperatures attained. The initial content of volatiles besides water corresponds to a bulk mass density of 0.07 g cm−3 (half CO and half CO2 , see Table II). The change of this density throughout the nucleus with time is shown in Fig. 6, top and bottom, for models a–d and a ∪ e–g, respectively. We draw attention to the obvious difference in the nature of the two energy sources: insolation, which is a surface source, generates a sublimation front of CO ice, which advances inward. This is illustrated by the slope of the boundary between the CO-rich and CO-depleted zones in all panels of Fig. 6. Radiogenic heating (of homogeneously distributed radioactives) is a volume source and affects most of the comet’s volume simultaneously; this is shown by the sublimation of CO2 , whose gradually decreasing density is illustrated by the vertical, color-changing fronts in Fig. 6. The CO ice is lost, sooner or later, in all models. The CO2 ice is retained in a surface layer in all models, but the thickness of this CO2 -rich layer varies considerably among them: all the CO2 is retained (everywhere) if no (or very little) 26 Al is included; only the outer few hundred meters retain the CO2 ice in the baseline model. Typically, the CO2 -rich layer is about 1 km thick, almost regardless of the model’s radius, distance, and X 0 (26 Al), > 10−8 ). In some cases, some of the CO2 gas provided X 0 (26 Al ∼ released in the interior refreezes at the lower boundary of the CO2 -rich outer layer, raising the CO2 density there by up to a factor of 4. The region that retains the CO2 ice also retains the water ice in amorphous form. Below it there is a ∼100-m thick layer of amorphous ice, devoid of CO2 ice, and below that only crystalline water ice and dust. The loss of volatiles from the interior naturally affects the porosity. The evolution of porosity throughout the nucleus is

shown in Fig. 7, top and bottom, for models a–d and a ∪ e–g, respectively. The ordinate in these plots is the relative volume (V (r )/V (R)) rather than the logarithm of depth, which means that the surface is now at the top. In addition, each panel has its own gray scale, to emphasize individual details. When very little or no 26 Al is included and only the CO ice is lost, the resulting increase in porosity is small (about 7%, from 0.533 to 0.565—see Table II) and the porosity remains almost homogeneous. When the CO2 ice is lost as well, regions of relatively high and relatively low porosities are obtained due to evaporation of water on the one hand and refreezing of CO2 on the other. Thus the porosity ranges from as low a value as 0.2 to almost 0.8. These changes should, however, be taken as indications of the expected trends, since the assumption of very high permeability is bound to affect these results far more than it has affected the previous results. The range of possible porosities and bulk densities that could arise from evolution of the initial structure assumed (for all models) is summarized in Table IV. It is noteworthy that large variations in porosity arise in a relatively thin region (in terms of radius), at depths ranging among models from a few hundred meters to 1 km, which means that a mechanically “weak” region forms at that depth below the surface. This may TABLE IV Porosity and Density Ranges Composition

Porosity

Density g cm−3

Dust + H2 O + CO2 + CO Dust + H2 O + CO2 Dust + H2 O Dust Dust + H2 O-filled pores Dust + CO2 -filled pores

0.533 0.565 0.587 0.892 0 0

0.70 0.665 0.63 0.35 1.17 1.74

308

CHOI ET AL.

FIG. 5. Color map of the evolving temperature profile for models a, b, c, d—top panels, as marked—and models a, e, f, g—bottom panels, as marked. The change in color along a vertical line represents the radial temperature variation at the corresponding time. The change in color along a horizontal line represents the change of temperature with time at the corresponding (fixed) depth. The time and depth scales are logarithmic, since changes are more pronounced during the early stages of evolution, and sharper, generally, toward the surface of the nucleus.

FIG. 6. Color map of the total mass density of CO and CO2 ices, as it changes with time throughout the nucleus, for the same models (as marked) and on the same time and depth scales as those in Fig. 5. The initial CO + CO2 is 0.07 g cm−3 (see Table II); the density drops to 0.035 g cm−3 , when the CO ice has completely evaporated and only the CO2 ice is left, and drops to zero, when both these volatile species have evaporated. Density values higher than 0.07 g cm−3 result from refreezing of CO2 gas, as it flows outward through cold layers.

THERMAL EVOLUTION IN THE KUIPER BELT

309

different assumptions regarding physical characteristics, but retaining the same initial porosity. The conclusions of these tests may be summarized as follows:

FIG. 7. Gray-scale map of the changing porosity for models a, b, c, d—top panels, as marked—and models a, e, f, g—bottom panels, as marked—where the ordinate is relative volume, rather than depth, providing a different perspective of the models’ structures. It is meant to indicate the mass distribution (since porosity is a volume-related parameter). Individual scale bars are used in this case.

be significant for the effect of collisions, which are relatively frequent in the Kuiper Belt, bearing on the sizes of broken-off fragments. So far, we have performed a systematic parameter study, involving only three parameters. In fact, however, the number of uncertain parameters is larger and it would be impossible to account for all of them in a methodical manner. Therefore, in order to test the robustness of our results regarding the change of internal structure, we have computed additional test models, with

(a) When gases trapped in amorphous ice are considered, in addition to those included in frozen form, the only long-term change is that more volatiles can freeze, eventually, in the outer layers of the nucleus, reducing the porosity. However, this effect depends entirely on the initial amounts of trapped and frozen volatiles and it would be impossible to deduce the origins of the ice in the outer layers of the nucleus based on the final composition there. (b) We have run models for which we solved the full set of time dependent equations, adopting a very high permeability (large average pore size)—an option mentioned in Section 3— for comparison. The grid spacing in the outer part of such models is far coarser, since the mass-shell thickness is limited from below by the pore size. Therefore the surface temperature is less accurate. As expected, advection surpasses conduction, sometimes by a large factor, in this case. Nevertheless, we do not find a marked difference in the thermal evolution of these models as compared to the former models. Our calculations show that among the different terms of the heat equation, conduction and advection are smaller than the other terms during the decay of 26 Al, meaning that the energy released is mainly absorbed locally, in evaporation. This strengthens our conclusions regarding the internal structure of KBOs. (c) For the test models we also find that larger amounts of gas flowing outward from the interior freeze in the cold outer regions, reducing the porosity there to very low values. For example, for a 500-km model (similar to model g) a porosity of only about 0.06 may be reached after 1 Myr of evolution, in an outer layer about 1 km thick. This is because the internal gas densities are higher, meaning that a gas molecule dwells in the nucleus interior for a longer period of time, which increases the probability of its refreezing. (d) Increasing the baseline heliocentric distance from 30 to 40 AU has no significant effect. Although the surface temperature is lower, we find that the volatile CO is still completely lost. We should mention here that the heliocentric distance adopted is, in fact, only a measure of the solar radiation. The true parameter of these models is L  /dH2 . Since the solar luminosity was weaker in the early stages of the Solar System, whereas we have assumed here the present value of L  , the true distances corresponding to our nominal distances may, in fact, be smaller by as much as 15%. Hence the fact that a 25% change in distance has no significant effect is encouraging. A detailed comparison with other similar studies is very difficult, if not impossible, to perform for two reasons: first, the approximations and basic assumptions employed differ among models, and second, the free parameter values adopted are different. Nevertheless, we find that, broadly, the results of different calculations agree and complement each other. As mentioned in Section 1, a calculation similar to the present calculation was

310

CHOI ET AL.

recently performed by De Sanctis et al. (2001), for two particular KBO orbits and radii and for a limited evolution time. As a result, the effect of radioactive heating was less apparent in that study. Depletion of volatiles, differentiation of the initially homogeneous composition, and refreezing are the main results of both these calculations. Our results are also in agreement with those of the extensive parameter study performed by Prialnik and Podolak (1995), although that study was devoted mainly to the process of crystallization and release of gas trapped in amorphous ice, whereas here volatile species are included as ice mixtures. 5. CONCLUSIONS

We have carried out long-term evolutionary calculations for models simulating potential KBOs. A composition of porous ice and dust was assumed, where the ice was a mixture of H2 O—the most abundant species—and CO and CO2 in smaller amounts. The main purpose of the study was to determine whether these more volatile species could survive radioactive heating during the early stages of evolution and whether a presumably homogeneous initial structure and composition could be preserved. To account for large vents or channels that should allow gases to escape—and also in order to facilitate the computations—we have adopted a quasi-steady-state approximation for the gaseous phases. Spherical symmetry was also assumed, but this is perfectly justified in long-term evolution calculations, especially when the major heat source is homogeneously distributed over the mass. Based on these assumptions and on the set of parameters adopted, and allowing for the approximations employed, our main conclusions are as follows: 1. The internal temperature profile of KBOs may have been substantially affected by both short- and long-lived radionuclides, with accompanying changes in composition and structure. The evolution of the temperature profile and the structural modifications are a strong function of the accretion times of KBOs, their sizes, the dust-to-ice mass fraction, and, to some extent, the heliocentric distances. 2. KBOs may have lost entirely all volatiles that sublimate below ∼40–50 K, which were initially included as ices. This can occur even in the absence of 26 Al. However, in this case the conclusion is valid only on the assumption that the entire surface area is—on average—equally heated. Depending on the angle of the rotation axis of the object, it may well be that some areas are always concealed from the Sun and stay at very low temperatures at all times. But it seems unreasonable to assume that rotation of a small (and probably uneven) body would remain undisturbed for billions of years. Thus, it is highly unlikely that KBOs still retain extremely volatile species in icy form. This leads to the further conclusion that if comets originating in the Kuiper Belt emit such volatiles (e.g., CO) these must have been trapped and retained in the amorphous ice, meaning that the ice of KBOs must have been initially amorphous.

26

Molecule H2 O CO CO2 NH3 CH4 H2 CO N2

TABLE V Al Required for Evaporation of Volatiles Tsubl (K)

H (erg g−1 )

X 0 (26 Al)

120 20 70 83 30 73 21

2.8(10) 2.3(9) 5.9(9) 1.7(10) 6.2(9) 8.2(9) 2.5(9)

1.0(−7) 8.2(−9) 2.4(−8) 6.2(−8) 2.2(−8) 3.1(−8) 8.9(−9)

3. KBOs may have partially lost less volatile ices as well. The temperatures at which sublimation of different volatiles sets in (Tsubl ), the corresponding latent heat, and the mass fraction of 26 Al that would be required in order to raise the temperature of cometary material from initially 10 K to Tsubl and sublimate various volatiles, assuming they constitute 10% of the ice, are listed in Table V. These are, of course, lower limits to the necessary amounts of 26 Al, since some of the heat is conducted outward and “wasted.” Thus, CO, as well as N2 and possibly methane, are expected to be lost entirely, but the moderately volatile species, such as CO2 , H2 CO, and NH3 , should be partially retained. As for water, only a relatively small fraction could be evaporated. 4. As a result, the internal structure of KBOs is most probably not uniform; rather, density, porosity, H2 O ice phase, and strength, vary with depth. Generally, the porosity should increase with depth (since compaction is unlikely to occur in small bodies), but not necessarily monotonically. In particular, we have found that weak regions—by which we mean regions of sharp density variations—may form at depths ranging among models from several hundred meters to 1 km. Since breaking occurs preferentially at, or along, a weak zone, it is reasonable to deduce that fragments breaking off KBOs upon collision are likely to have sizes on the order of 1 km. Now, collisions in the Kuiper Belt are believed to be responsible for the Jupiter-family comets and thus, indulging in speculation, we may have found a reason for the relatively small sizes of these comets (see also Farinella and Davis 1996, Davis and Farinella 2001). 5. Not only is the structure not uniform, but also, similarly, the internal composition of KBOs is most probably not homogeneous, but stratified (see also De Sanctis et al. 2001), with the outer layers being less altered by evolution. We have found that regions enriched in volatile species, as compared with the initial abundances assumed, arise due to gas migration and refreezing. This effect should be more pronounced, if the assumption of very high permeability (in different approximations) is relaxed. When such regions are exposed in comets originating in the Kuiper Belt, after the overlying material has been eroded, vigorous evaporation of these volatile species may result in outbursts. Our results should be taken to indicate evolutionary trends and our conclusions have been formulated as qualitative. An uneven shape, rotation, orbital migration, and collisions are only a few

THERMAL EVOLUTION IN THE KUIPER BELT

among the factors that may affect the structure of KBOs, but are difficult—if not impossible—to be accounted for in a systematic manner. Nevertheless, we may state with a reasonable degree of confidence that if KBOs did experience radioactive heating, their structure and composition were altered mainly to the extent of considerable loss of volatiles and significant departure from internal homogeneity. ACKNOWLEDGMENTS We thank Paul Weissman for a very careful reading of the original manuscript and for many valuable comments an suggestions. R.M. acknowledges the support of the Minerva Foundation.

REFERENCES Akridge, G., P. H. Benoit, and D. W. G. Sears 1998. Regolith and megaregolith formation of H-chondrites: Thermal constraints on the parent body. Icarus 132, 185–195. Brown, R. H., and C. C. Koresko 1998. Detection of water ice on the Centaur 1997 CU26. Astrophys. J. 505, L65–L67. Brown, R. H., D. P. Cruikshank, Y. Pendleton, and G. J. Veeder 1998. Identification of water ice on the Centaur 1997 CU26. Science 280, 1430– 1432. Brown, R. H., D. P. Cruikshank, and Y. Pendleton 1999. Water ice on Kuiper Belt Object 1996 TO66 . Astrophys. J. 519, L101–L104. Chiang, E. I., and M. E. Brown 1999. Keck pencil-beam survey for faint Kuiper Belt objects. Astron. J. 118, 1411–1422. Clayton, D. D., and M. D. Leising 1987. 26 Al in the interstellar medium. Phys. Rep. 144, 1–50. Consolmagno, G. J., and J. S. Lewis 1978. The evolution of icy satellite interiors and surfaces. Icarus 34, 280–293. Davis, D. R., and P. Farinella 1997. Collisional evolution of Edgeworth–Kuiper Belt Objects. Icarus 125, 50–60. Davis, D. R., and P. Farinella 2001. Collisional effects in the Edgeworth– Kuiper Belt. In Collisional Processes in the Solar System (M. Ya. Marov and H. Rickman, Eds.), pp. 277–285. Kluwer, Dordrecht. Delsanti, A. C., H. Boehnardt, L. Barrera, K. J. Meech, T. Sekiguchi, and O. R. Hainaut 2001. BVRI photometry of 27 Kuiper Belt objects with ESO/very large telescope. Astron. Astrophys. 380, 347–358. De Sanctis, M. C., M. T. Capria, and A. Coradini 2001. Thermal evolution and differentiation of Edgeworth–Kuiper Belt objects. Astron. J. 121, 2792–2799. Duncan, M., T. Quinn, and S. Tremaine 1988. The origin of short-period comets. Astrophys. J. 328, L69–L73. Ellsworth, K., and G. Schubert 1983. Saturn’s icy satellites: Thermal and structural models. Icarus 54, 490–510. Farinella, P., and D. R. Davis 1996. Short-period comets: Primordial bodies or collisional fragments. Science 273, 938–941. Farinella, P., D. R. Davis, and S. A. Stern 2000. Formation and collisional evolution of the Edgeworth–Kuiper Belt. In Protostars and Planets IV (V. Mannings, A. P. Boss, and S. S. Russell, Eds.), pp. 1255–1282. Univ. of Arizona Press, Tucson. Ghosh, A., and H. Y. McSween 1998. A thermal model for the differentiation of Asteroid 4 Vesta, based on radiogenic heating. Icarus 134, 187–206. Gladman, B., and M. Duncan 1990. On the fates of minor bodies in the outer Solar System. Astron. J. 100, 1680–1693. Gladman, B., J. J. Kavelaars, P. D. Nicholson, T. J. Loredo, and J. A. Burns 1998. Pencil-beam surveys for faint trans-neptunian objects. Astron. J. 116, 2042–2054.

311

Greenberg, J. M., H. Mizutani, and T. Yamamoto 1995. A new derivation of the tensile strength of cometary nuclei: Application to Comet Shoemaker–Levy 9. Astron. Astrophys. 295, L35–L38. Hainaut, O. R., C. E. Delahode, H. Boehnardt, E. Dotto, M. A. Barucci, K. J. Meech, J. M. Bauer, R. M. West, and A. Doressoundiram 2000. Physical properties of TNO 1996 TO66 —Lightcurves and possible cometary activity. Astron. Astrophys. 356, 1076–1088. Haruyama, J., T. Yamamoto, H. Mizutani, and J. M. Greenberg 1993. Thermal history of comets during residence in the Oort cloud: Effect of radiogenic heating in combination with the very low thermal conductivity of amorphous ice. J. Geophys. Res. 98, 15,079–15,088. Jewitt, D., J. Luu, and C. Trujillo 1998. Large Kuiper Belt objects: The Mauna Kea 8K CCD survey. Astron. J. 115, 2125–2135. Jewitt, D., H. Aussel, and A. Evans 2001. The size and albedo of the Kuiper Belt Object (2000) Varuna. Nature 411, 446–447. Kenyon, S. J., and J. X. Luu 1998. Accretion in the early Kuiper Belt. I. Coagulation and velocity evolution Astron. J. 115, 2136–2160. Klinger, J. 1980. Influence of a phase transition of ice on the heat and mass balance of comets. Science 209, 271–272. Klinger, J., S. Espinasse, and B. Schmidt 1989. Some considerations on cohesive forces in sun-grazing comets. Europ. Space Agency Spec. Publ. 12, 197–200. Kochan, H., B. Feuerbacker, F. Joo, J. Klinger, W. Sebolt, A. Bischoff, H. Duren, D. Stuffler, T. Spohn, H. Fechtig, E. Gruen, H. Kohl, D. Krankowsky, K. Roessler, K. Thiel, G. Schwehm, and U. Weishaupt 1989. Comet simulation experiments in the DFVLR space simulators. Adv. Space Res. 9, 113–121. Lee, T., D. Papanastassiou, and G. J. Wasserburg 1976. Demonstration of 26 Mg excess in Allende and evidence for 26 Al. Geophys. Res. Lett. 3, 109–112. Luu, J. X., D. C. Jewitt, and C. Trujillo 2000. Water ice in 2060 Chiron and its implications for Centaurs and Kuiper Belt objects. Astrophys. J. 531, L151– L154. MacDonald, G. J. 1959. Calculations on the thermal history of the Earth. J. Geophys. Res. 64, 1967–2000. Mahoney, W. A., J. C. Ling, Wm. A. Wheaton, and A. S. Jacobson 1984. Heao 3 discovery of 26 Al in the interstellar medium. Astrophys. J. 286, 578–585. Mekler, Y., D. Prialnik, and M. Podolak 1990. Evaporation from a porous comet nucleus. Astrophys. J. 356, 682–686. Miyamoto, M., N. Fujii, and H. Takeda 1981. Ordinary chondrite parent body: An internal heating model. Proc. Lunar Planet. Sci. Conf. 12B, 1145–1152. Orosei, R., A. Coradini, M. C. De Sanctis, and C. Federico 2001. Collisioninduced thermal evolution of a comet nucleus in the Edgeworth–Kuiper Belt. Adv. Space Res. 28, 1563–1569. Podolak, M., and D. Prialnik 1997. 26 Al and liquid water environments in comet. In Comets and the Origin of Life (P. Thomas, C. Chyba, and C. McKay, Eds.), pp. 259–272. Springer-Verlag, New York. Pollack, J. B., O. Hubickyj, P. Bodenheimer, J. J. Lissauer, M. Podolak, and Y. Greenzweig 1996. Formation of the giant planets by concurrent accretion of solids and gas. Icarus 124, 62–85. Prialnik, D. 1992. Crystallization, sublimation, and gas release in the interior of a porous comet nucleus. Astrophys. J. 388, 196–202. Prialnik, D., and A. Bar-Nun 1990. Heating and melting of small icy satellites by the decay of 26 Al. Astrophys. J. 355, 281–286. Prialnik, D., and M. Podolak 1995. Radioactive heating of porous cometary nuclei. Icarus 117, 420–430. Prialnik, D., and M. Podolak 1999. Changes in the structure of comet nuclei due to radioactive heating. Space Sci. Rev. 90, 169–178. Prialnik, D., A. Bar-Nun, and M. Podolak 1987. Radiogenic heating of comets by 26 Al and implications for their time of formation. Astrophys. J. 319, 993– 1002. Prialnik, D., U. Egozi, A. Bar-Nun, M. Podolak, and Y. Greenzweig 1993. On pore size and fracture in gas-laden comet nuclei. Icarus 106, 499–507.

312

CHOI ET AL.

Schmitt, B., S. Espinasse, R. J. A. Grim, J. M. Greenberg, and J. Klinger 1989. Laboratory studies of cometary ice analogues. Europ. Space Agency Spec. Publ. 302, 65–69.

Urey, H. C. 1955. The cosmic abundances of potassium, uranium and thorium and the heat balances of the Earth, the Moon and Mars. Proc. Nat. Acad. Sci. 41, 127–144.

Schubert, G., D. J. Stevenson, and K. Ellsworth 1981. Internal structures of the Galilean satellites. Icarus 47, 46–59.

Wallis, M. K. 1980. Radiogenic heating of primordial comet interiors. Nature 284, 431–433.

Share, G. H., R. L. Kinzer, J. D. Kurfess, D. J. Forrest, E. L. Chupp, and E. Rieger 1985. Detection of galactic 26 Al gamma radiation by the SMM spectrometer. Astrophys. J. 292, L61–L65.

Weidenschilling, S. J. 1988. Formation processes and time scales for meteorite parent bodies. In Meteorites and the Early Solar System (J. F. Kerridge and S. W. Metthews, Eds.), pp. 348–371. Univ. of Arizona Press, Tucson.

Sheppard, S. S., D. C. Jewitt, C. Trujillo, M. J. I. Brown, and M. C. B. Ashley 2000. A wide-field CCD survey for centaurs and Kuiper Belt objects. Astron. J. 120, 2687–2694.

Weidenschilling, S. J. 1997. The origin of comets in the solar nebula: A unified model. Icarus 127, 290–306.

Sirono, S., and J. M. Greenberg 2000. Do cometesimal collisions lead to bound rubble piles or to aggregates held together by gravity? Icarus 145, 230–238.

Wetherill, G. W., and G. R. Stewart 1989. Accumulation of a swarm of small planetesimals. Icarus 77, 330–357.

Srinivasan, G., J. N. Goswami, and N. Bhandari 1999. 26 Al in Eucrite Piplia Kalan: Plausible heat source and formation chronology. Science 284, 1348– 1350.

Whipple, F. L., and R. P. Stefanic 1966. On the physics and splitting of cometary nuclei. Mem. R. Soc. Liege (Ser. 5) 12, 33–52.

Trujillo, C., J. X. Luu, A. S. Bosh, and J. L. Elliot 2001. Large bodies in the Kuiper Belt. Astron. J. 122, 2740–2748.

Yabushita, S. 1993. Thermal evolution of cometary nuclei by radioactive heating and possible formation of organic chemicals. Mon. Not. R. Astron. Soc. 260, 819–825.