Signs for Heavy Bombardment in Debris Disks

DEBRIS DISKS AND THE FORMATION OF PLANETS: A SYMPOSIUM IN MEMORY OF FRED GILLETT ASP Conference Series, Vol. 324, 2004 L. Caroff, L.J. Moon, D. Backman and E. Praton, eds.

Signs for Heavy Bombardment in Debris Disks C. Dominik Sterrenkundig Instituut ‘Anton Pannekoek’, Kruislaan 403, NL-1098 SJ Amsterdam, The Netherlands J. Bouwman Service d’Astrophysique, CEA/DSM/DAPNIA, C.E. Saclay, F-91191 Gif-sur-Yvette, France Abstract. We explore the decay of the IR luminosity of debris disk with time as it was observed by observations with the ISO satellite. The observations show a steep decrease of the fractional luminosity, and after about 400 Myrs, the excess becomes undetectable for most stars. We describe a collisional model which predicts a slope of –1 in a collisionally dominated disk, and a slope of –2 in a radiation force dominated disk. Since debris disks are generally collisionally dominated, the observed slope of –1.7 is too large. Additional planetesimal removal processes are required. Finally we discuss the spectrum of the Herbig Ae star HD 100546. This spectrum contains crystalline silicate emissions which are different from those seen in other similar stars. We speculate that this material may not have been thermally annealed as small grains, but may be the result of the destruction of differentiated bodies.

1.

Introduction

Immediately after the discovery of the dust shell around Vega (Aumann et al. 1984) it was realized that this may be the first detection of material directly related to planetary systems around nearby stars. This discovery happened more than 10 years before the first detection of an extrasolar planet (Mayor and Queloz 1995). Additional importance of this early discovery is that it detected material related to small bodies in the system like comets and asteroids, bodies so small that they render attempts of direct detection hopeless. The IRAS discovery, in which Fred Gillett played the leading role, has opened a window into foreign planetary systems. Of course, comparison with our own solar system is one of the most intriguing possibilities. The task put to me by the organizers of this symposium was to describe possible links of the IRAS and later ISO observations of circumstellar dust to the heavy bombardment phase in the solar system. Such a link has been suggested by Habing et al. (1999), because we found that disks around nearby main-sequence stars escape detectability after about 400Myrs, a time scale comparable to the duration of the Heavy Bombardment documented by the cratering record on the Moon. Often the term Late Heavy Bombardment is used to describe a late peak in the Bombardment of the inner Planets (see Levison, this volume) - there is no clear indication for such an event in our ISO data. Therefore, when we use the term Heavy Bombardment, we refer to the 121

122

Dominik and Bouwman

Figure 1. Cululative number of disk detections as a function of stellar age. Age is represented by age rank within the sample. Slopes in this diagram indicate directly detection probabilities for a given age. From Habing et al. (1999).

extended bombardment which started in the accumulation phase of the planets and ended about 600Myrs later. In this contribution we will briefly summarize the ISO results which were obtained by the Habing et al. consortium and discuss these results in terms of a collisional model. Finally, we will discuss the ISO spectrum of the Herbig star HD 100546 in which we find evidence which may be interpreted in terms of the breakup of planetesimals in a highly developed circumstellar disk.

2.

Results from ISO Observations

The IRAS discovery of main-sequence stars with circumstellar dust did leave many questions unanswered. In particular, any dependence on stellar age remained hidden by the rather random sampling caused by the IRAS detection limits. With the ISO satellite, several groups set out to create better samples and to derive the age dependence of the amount of dust seen around these stars. There are two main approaches. The proposal from Becklin et al. tried to probe different ages by observing stars in different clusters (Spangler et al. 2001). Beckwith et al. did similar work with very young clusters (Meyer & Beckwith 2000). The main advantage is that the ages of clusters can be determined quite accurately. The disadvantage of this approach is that most clusters are far away, so that only disks with large excesses can be securely identified while small ex-

Heavy Bombardment in Debris Disks

123

cesses will be difficult to distinguish from photospheric fluxes. For a discussion of the results from these programs, see Silverstone (this volume). The Habing et al. sample was both volume and flux limited to ensure both random sampling of the solar vicinity and secure detection of even photospheric fluxes (Habing et al 2001). The disadvantage of this approach is of course that determining ages of single field stars is very difficult and uncertain. We observed 84 stars at 25 and 60 µm and used only the 60 µm measurements to determine the fractional dust luminosity fd =

Ldust L⋆

(1)

which is a standard measure of the amount of dust in these systems. We accepted only measurements where both the detection itself and the excess are measured with S/N> ∼ 3. Therefore all detections reported are secure IR excesses, even though chance alignment with background sources cannot be excluded. We determined stellar ages based on isochrones (for A and F stars) and based on a variety of methods (metallicity, rotation, kinematics) for late type stars (Lachaume et al. 1999). The ages for the A stars carry errors of about a factor of two, the ages for late-type stars are significantly more uncertain (Stauffer, this volume). A simple test for age dependency is to determine if there is a dependence of disk detection frequency on age. In order to reduce the noise produced by the uncertainty in age, we plotted the cumulative number of detected disks as a function of age rank (i.e. the number of stars in an age-sorted table). The results are shown in Fig. 1. The curve has a slope of 1 for the youngest five stars, and an average slope of 0.6 for all stars younger than 400 Myrs and a slope of 0.09 for stars older than this limit. These numbers indicate that all very young stars have an IR excess (i.e. probably a disk), while 60% of stars below 400 Myrs and only 9% of the older stars show similar excesses. These results become also clear when the fractional luminosity is directly plotted against the stellar age (Fig. 2). This figure shows an initial steep decrease of fd with age, with a slope of about 1.7, in accordance with the results found by Spangler et al. (2001) for their cluster data. This decreasing line intercepts our detection limit of log fd ≈ −5 at an age of about 400Myr, producing the slope change in Fig. 1. However, we would like to point out, that there are several caveats to this simple interpretation. 1. The determination of the slope of 1.7 in the single-star diagram depends largely on very few very young stars. In our sample it is only β Pic, in the Spangler et al. sample also HR 4796A. 2. Similarly, the slope determination in the cluster data depends largely on the data for the Chameleon cluster which contains many stars with optically thick disks, where fd should not be representative of the total dust mass, but rather of the solid angle covered by the (flaring) disk surface. 3. Our sample contains a significant number of stars with rather old ages, but excesses corresponding to fd values of 10−4 , similar to Vega. If there is indeed an intrinsic decay of the disk, a major problem remains to explain these old disks.

124

log fd

Dominik and Bouwman

Figure 2. Fractional luminosity of the IR excess as a function of age. Dots indicate 3-σ detections, crosses indicate 1-σ upper limits. From Habing et al. (2001), but with new age for β Pic (20 Myr).

One possible explanation for the “old” disks may be that these are either misidentified excesses, or due to incorrect ages. For example, the IR excess of the star 55 Cnc (Dominik et al. 1998) may well be due to a background object. The submm excess originally assigned to this object (Jayawardhana et al. 2000) has now been found to be due to several background sources (Jayawardhana 2002). Similarly, the star HD 207129 with a significant excess (Walker & Wolstencroft 1988, Jourdain de Muizon et al. 1999) and a rotational age of about 5 Gyrs has recently been proposed as a member of a young moving group (Montes et al. 2001). However, these are only two cases, more Vega-like stars exist with suggested high ages in the Gyr regime. More work on studying the ages of these stars more closely is very important in order to understand the mechanisms which create and maintain dust disks around main sequence stars. 3.

Collisional Model

We will report here some results of a study of the collisional evolution of debris disks (Dominik & Decin 2003). We consider the simplest case of a distribution of planetesimals which develops collisionally. For simplicity we further assume that the velocity distribution of the planetesimals remains constant. So the basic picture is that of a disk where after the formation of planets the remaining planetesimal population is quickly stirred and then left alone without further stirring or damping. We assume that the dust produced in the system is the tail of a collisional equilibrium size distribution which is ultimately due to the destruction of the planetesimals which

125


we will call comets for convenience. Let Np be the number of such bodies with radius ap , collisional cross section σp and mass Mp . Let these bodies occupy a volume V where they move around and collide with relative velocities vcoll = νvK with vK being the Kepler velocity in the regions where the comets exist. ν will be zero for circular orbits in a single plane and approach one for highly excited orbits. Note that significant stirring will increase the total volume visited by the comets, so V will not be independent of ν. If all collisions are destructive, the time derivative of the number of comets in the system is given by N˙ p = −

Np2 ts

(2)

where ts is the sweeping time, the time for the cross section of a comet moving at velocity vcoll to sweep the entire volume V . The solution of eq. 1 is Np (t) =

N0 1 + Nts0 t

(3)

where N0 is the initial number of comets. For times t ≫ becomes independent of N0 . Np (t) →

ts t

for t ≫

ts N0

.

ts N0 ,

the solution

(4)

Thus, the number of comets is given by a powerlaw, but only after an initial period in which the number is in fact almost constant. Since the breakpoint, ts /N0 depends upon N0 , the time needed to reach the powerlaw behavior is dependent on the initial mass in the planetesimal disk. The higher the mass, the faster the powerlaw takes over. We may now proceed to compute the amount of dust grains produced by the collisional cascade of the comets. It has been show that collisional equilibrium produces a size distribution f (m) ∝ m−11/6 (Dohnanyi 1969). Therefore, the production rate of small particles in this cascade is equal to the destruction rate of the large bodies. Thus, if the gains we are interested in (because they dominate the visibility of the dust) have a radius avis , then the production rate of these grains is given by Rprod

= N˙ p

µ

ap avis

¶3

.

(5)

At the small end of the size distribution, particles are removed. The main removal processes are radiation pressure for grains which are small enough to receive an acceleration close to or higher than the gravitational acceleration of the star. Somewhat bigger grains can still be removed from the disk by Poynting Robertson drag. Which process dominates for the visible is dependent upon the dust density in the disk. If the density is high, collisional destruction dominates all the way down to the smallest particles. If the disk contains less mass, particles of a certain size have Poynting Robertson timescales which are smaller than the collisional time and are being pulled out of the cometary region towards the star. It turns out that in a system with as little mass as the solar Kuiper

126

Dominik and Bouwman

Figure 3. Fractional luminosity decrease as a function of age in the model calculation. The different curves are for different initial masses: from bottom to top 0.01, 0.1, 1, 10, 100 M⊕ . Note the slope change at fd ≈ 5 × 10−8 .

Belt (estimates indicate about 0.1 M⊕ ), Poynting Robertson drag dominates for 1 µm sized grains. In the denser Vega-like systems collisions always dominate (Artymowicz & Clampin 1997). In a collisionally dominated disk, the main removal process of small grains is collisions among the grains themselves. The removal rate is therefore proportional to the number of these grains squared coll Rloss ∝ n2vis

.

(6)

In the radiation dominated regime, the removal rate is only linearly proportional to the number of grains, thus coll Rloss ∝ nvis

.

(7)

Inserting eq. (4) into eqs. (6) and (7) and solving for steady state (Rloss = Rprod ) shows n ∝ t−1 n ∝ t−2

in the collisionally dominated regime in the radiatively dominated regime

(8) (9)

By including more details about the collisions and about radiation forces on grains, one can use these results to compute the expected fd values. We show the results of such a calculation in figure 3 where we have assumed a Kuiper-Belt like structure with different initial mass. The size of the comets was assumed to be ap = 1 km, and ν = 1, i.e. Keplerian collision velocities corresponding to highly excited orbits. The plot shows that the time to reach collisional equilibrium depends strongly on the initial masses. With only 0.01M⊕ , it takes several Gyrs before the comet population develops like a powerlaw. For 10M⊕ , however,


127

equilibrium is reached after 108 years. Still higher masses are already fully on the powerlaw slope at 107 years. The powerlaw slope is indeed only -1 (collisionally dominated grains). But at fd = 5 × 10−8 , the radiation starts removing grains faster that collisions would, so the slope changes to -2. If we compare this result with the observed slope of -1.7 (Spangler et al. 2001), it is clear that this slope is steeper than a pure collisional evolution does suggest. Therefore, additional comet removal processes like ejection due to planetary encounters must be operating at the same time in order to speed up comet removal. In this sense, the processes in the debris disks seem to be indeed similar to the development in the solar system during the Heavy Bombardment phase. 4.

Crystalline Silicates in HD 100546

Herbig AeBe (HAEBE) stars are the massive counterparts of T Tauri stars. They are pre-main-sequence stars of spectral type A and B, and they possess an IR excess which is thought to originate from a disk. The disk must still contain gas since the SED of these stars shows that about 50% of the stellar luminosity is absorbed and reprocessed by the CS material Only a gas-rich disk can reach and sustain the necessary scale heights to intercept 50% of the stellar radiation. HAEBE stars are usually brighter than T Tauri stars and therefore ideal objects for studying the IR emission spectrum in great detail with ISO. Because of its greater sensitivity, SIRTF (Spitzer) will allow similar studies also for T Tauri stars. HAEBE stars, in particular the Be stars, are often still surrounded by molecular cloud material, so disentangling the IR emission due to that material from the emission from the disk can be difficult. However, there exists a class of HAe stars which is more isolated. These stars share all important properties with the other HAEBE stars, but the IR emission is compact and very likely associated with a disk (Waters and Waelkens 1998). A sample of 14 such stars has been observed with ISO (Meeus et al. 2001). The possibly most suspicious of the sources is the B9V star HD 100546. This is the oldest star of the sample (10 Myr, van den Ancker et al. 1997). The spectrum (Fig. 4) of the star stands out because of several reasons (Bouwman et al. 2003): 1. HD 100546 shows the strongest crystalline silicate (Forsterite) features of all observed Herbig stars. At least 10% of the small grains in the disk surface are crystalline. 2. Modelling of the spectrum shows that the crystalline silicates are cold, with a typical temperature of only 200K, much too cold for thermal annealing to occur. 3. The crystalline silicates in HD 100546 are chemically different from those seen in other sources (Bouwman et al. 2001). In most HAEBE stars, any 11.3 feature of Forsterite is accompanied by an 8.6 µm shoulder of silica (the strength of both features is proportional). The 8.6 µm shoulder is completely absent in HD 100546, suggesting that the crystalline dust grains in this source have a different chemical and/or thermal history.

128

Dominik and Bouwman 300 10

200

Fν [Jy]

8 6 4

F ν [Jy]

2

100

0 2

Ol 3

4

5 6 λ [µm]

7

8

0

Fo −100 10

100 λ [µm]

Figure 4. The spectrum of HD 100546, decomposed into the contributing dust materials. Note the extremely strong rise at 20 µm due to the amorphous Olivine component (Ol). Note also the very strong features of crystalline Forsterite (Fo). From Bouwman et al. 2003.

4. Modelling also shows that the degree of crystallinity is increasing outwards in the disk. Therefore, annealing close to the star with subsequent radial mixing cannot be the explanation for the cold grains. Radial mixing can produce a constant fraction of crystalline material at best. 5. The distribution of energy in the SED between near IR (< 7 µm) and mid IR (> 7 µm) is very different from other HAEBE stars like AB Aur. While AB Aur has similar amounts of emission in both regimes, HD 100546 emits 70% of it reprocessed radiation in the mid-IR and only 30% in the near IR. This points to an unusually low mass content of the inner disk region (Dominik et al. 2002). We therefore conclude that the production of crystalline silicates in this source must have been a local event. Candidates for this event are processes similar to the ones which have produced Chondrules in the Solar Nebula, for example shock heating (Boss and Graham, 1993). Another possibility may be the destruction of large planetesimals. These planetesimals would have to be large enough to melt and become differentiated, in order to account for the low Mg content in the grains (< ∼ 5%, Malfait et al. 1999). The grains produced in a collisional cascade would be collected in the disk surface where they become readily observable. A total of about one Earth mass of dust would be needed. The high collisional velocities required for destruction of large bodies could be due to a planet stirring the local population of planetesimals. The gap created by such a planet would also provide a natural explanation for the low mass content of the inner disk

129


10 AU

0.3 AU

380 AU

Figure 5. Cartoon of the disk of HD 100546. A planet at 10AU causes destruction of planetesimals, and starvation of the inner disk. The dust created is distributed over the outer disk and seen as crystalline silicates in the spectrum. From Bouwman et al. (2002)

which is reflected in the large ratio of mid-to-near IR emission. Figure 5 shows a cartoon of the structure of the disk for this scenario. Acknowledgments. CD would like to thank the USRA and the SOFIA program for travel support. References Artymowicz P., Clampin M. 1997, ApJ 490, 863 Aumann H.H., Gillett F.C., Beichmann C.A., et al. 1984, ApJ 278, 23 Boss A.P., Graham J.A. 1993, Icarus 106, 168 Bouwman J., Meeus G., de Koter A., et al. 2001, A&A 375, 950 Bouwman J., de Koter A., Dominik C., Waters L.B.F.M. 2003 A&A 401, 577 Dohnanyi J.W. 1969, J. Geophys. Res. pp. 2531–2554 Dominik C., Decin G. 2003, ApJ, 598, 626 Dominik C., Dullemond C., Waters L., Walch S. 2002, A&A, 398, 607 Dominik C., Laureijs R.L., Jourdain de Muizon M., Habing H.J. 1998, A&A 329, L53 Habing H.J., Dominik C., Jourdain De Muizon M., et al. 1999, Nature 401, 456 Habing H.J., Dominik C., Jourdain de Muizon M., et al. 2001, A&A 365, 545 Jayawardhana R., Holland W.S., Greaves J.S., et al. 2000, ApJ 536, 425

130

Dominik and Bouwman

Jayawardhana R., Holland W.S., Kalas P., et al. 2002, ApJ 570, L93 Jourdain de Muizon M., Laureijs R.J., Dominik C., et al. 1999, A&A 350, 875 Lachaume R., Dominik C., Lanz T., Habing H.J. 1999 A&A 348, 897L Malfait K., Waelkens C., Bouwman J., de Koter A., Waters L.B.F.M.: 1999, A&A 345, 181 Mayor M., Queloz D. 1995, Nature 378, 355 Meeus G., Waters L.B.F.M., Bouwman J., et al. 2001, A&A 365, 476 Meyer M.R., Beckwith S.V.W. 2000, In: ISO Survey of a Dusty Universe, Proceedings of a Ringberg Workshop Held at Ringberg Castle, Tegernsee, Germany, 8-12 November 1999, Edited by D. Lemke, M. Stickel, and K. Wilke, Lecture Notes in Physics, vol. 548, p.341 Montes D., L´ opez-Santiago J., G´ alvez M.C., et al. 2001, MNRAS 328, 45 Spangler C., Sargent A.I., Silverstone M.D., Becklin E.E., Zuckerman B. 2001, ApJ 555, 932 van den Ancker M.E., The P.S., Tjin A Djie H.R.E., et al. 1997, A&A 324, L33 Walker H.J., Wolstencroft R.D. 1988, PASP 100, 1509 Waters L.B.F.M., Waelkens C. 1998, ARA&A 36, 233

Fajardo-Acosta: The silicate spectrum of Beta Pic looks rather strange. Is that from ISO? The ground based spectra show more structure, including PAH features. Dominik : Yes, all the silicate spectra in this plot are from ISO. Liseau: In the plot of age vs f, you said that the data for T-Tauri are not relevant because of the high optical depth in the disk. But farther out the disk must be optically thin? Dominik : I meant optically thick to the stellar radiation. Becklin: But, using mm data you find that there still seems to be a correlation between f and total mass in dust. Dominik : Indeed you can measure the mass in the mm, but the IR number gives only a lower limit to f. Becklin: Yes, but they still agree; that is, they are still correlated. Forrest: Can you and Murray reconcile the broken power law that you show (break at 470 Myr) with the –1.75 power law that he showed? Silverstone: In the relevant age range, our data set (i.e the young end of our data) is very sparse. It’s consistent with a single power law, but not convincing. Our data are not inconsistent with the break that Carsten shows. Dominik : I agree. The power law drops below the detection limit at the same age where we find the break in disk frequency. So, the two plots are actually consistent. Werner : We are seeing collisional destruction of dust in the middle of protoplanetary disks. It’s important to understand how the transition from a protoplanetary to a debris disk occurs. The composition should be telling us something. Grady: The uv data also show an underabundance of Fe relative to Si. This shows that the silicates being evaporated are Mg rich, i.e. Fe poor.


131

Wooden: With respect to mineralogy, the detection of non-crystalline features is very temperature-dependent. The predictions are not yet in agreement with the data. The current models are optically thin. More detailed radiative transfer calculations are just beginning to be made. There is lots of work to be done yet. Dominik : I agree that the optically thin models only provide the relative abundances of particles in the upper layers of the disks. I don’t know how to resolve the discrepancy between your results and ours. Li : How can we get PAH features in these systems? We see PAHs in HD100546 which has cometary silicates, but not in comet Hale-Bopp, which should be similar according to some interpretations. Dominik : It’s possible that the PAHs could be detected farther out in the disk than silicates can, because the PAHs can be excited by single photon processes, whereas the silicate features would not be seen out there. Meyer : Can John Stauffer comment on the ages here? Stauffer : For the high mass stars in the sample there is likely no problem with the ages. But it would still be useful to look at several open clusters in that age range, and use the spread in that data (as described earlier) to set error bars on the ages for the stars in this group. For low mass stars the situation is more problematical. A few detected objects have ages younger than can be inferred from rotation, so you really have determined an upper limit to the age. Song: HD207129 is a member of an association. So you can use that to update the age estimate. As to the uncertainty: the ages of early type stars are heavily dependent on their rotation rates. The rotation affects the colors, hence the assignment of stellar classification. The net effect is that the ages can be off by as much as 50%. Dominik : Even so, our conclusion would probably not be affected. Artymowicz : A possible mechanism for producing crystalline features might be: the large grains are moved in by PR drag and the small grains go out due to radiation pressure. The large grains can get annealed, break up through collisions, and then the resulting small grains move out. Dominik : Interesting idea. Weinberger : In models of the Solar System, in-situ shocks can make crystalline materials far away from the star, where the equilibrium temperature is low. Meyer : What do you see if you sort by spectral type? Dominik : The data are not inconsistent, but are really too sparse to provide any conclusions. Our main conclusion is based on the A stars. Silverstone: I agree with this for our data. SIRTF and SOFIA will help a lot in providing a larger range of detected sources. Grady: How many of the old stars are left with disks? Dominik : Not many in our sample. Two or three. I don’t know of others. Silverstone: One star in our data set, HD151044. But the age determination is a bit dicey.

132

Dominik and Bouwman

Dominik : SIRTF will help to find excesses around lower mass stars, which are older on average. These stars can’t blow away the dust, so the dust just piles up.