text:
# Statistical properties of visual binaries as tracers of the formation and early evolution of young stellar clusters

## 1. Introduction

The current paradigm for star formation is that most stars form in stellar clusters rather than in Taurus-like star-forming regions. Since binary systems are more frequent than single stars, the study of multiple systems in clusters is expected to provide us with some insight on the stellar formation process itself. For instance, it has been suggested that the gravitational collapse of molecular clouds leads to a binary fraction close to 100 %, which would decrease over time to about 60 %, the value observed on the Main Sequence (Ghez et al. 1993). The physical processes leading to such an evolution and the timescale over which they act still need to be asserted. Alternatively, the impact of local physical conditions in the molecular cloud, such as the magnetic field strength or gas temperature, on the output of the star formation process has still to be ascertain and might lead to different binary frequencies in various environments (Durisen & Sterzik 1994).

Observing stellar clusters with a wide range of physical parameters allows to perform a direct comparison between binary properties and model predictions. Clusters with ages ranging from one to several hundred million years are known within 500 pc from the Sun with some hosting many massive stars while other containing only low-mass stars. The stellar density can vary over more than an order of magnitude from cluster to cluster, as well as the physical conditions of the parent molecular cloud. In this contribution, we review the observed properties of visual binaries in clusters and compare them with cluster properties and with model predictions regarding the formation and evolution of stellar clusters.

## 2. Summary of observational results

### 2.1. Visual binary frequency

The advent of several high-angular imaging techniques makes the binary frequency a relatively straightforward property to derive in stellar populations. For instance, $`HST`$ imaging (Prosser et al. 1994), speckle interferometry (Leinert et al. 1993, Patience et al. 1998) and adaptive optics (Bouvier et al. 1997) have been used, mostly in the near-infrared in order to facilitate the detection of relatively faint and cold low-mass companions. The low-mass populations of extremely young clusters (the Trapezium cluster, IC 348) as well as of more evolved ones ($`\alpha `$ Per, the Pleiades, the Hyades, Præsepe) have been studied, which span an age range from $``$1 to several 100 Myr. The range of semi-major axes probed by these surveys usually extends from a few 10 to $``$2000 AU. The Hyades cluster, which is the closest to the Sun, is the only one for which visual binaries could be resolved down to about 5 AU.

Most surveys miss the faintest companions so that a statistical completeness correction has to be applied in order to estimate the true binary frequency from the observed one. This correction can be large, especially at small separations for adaptive optics surveys due to the separation-dependent detection limit inherent to this technique. The correction factor also depends upon assumptions regarding the distribution of mass ratios and/or orbital periods of the binary population, which obviously increases the statistical uncertainties on the derived binary frequency.

With these limitations, all studied clusters are found to exhibit a visual binary frequency consistent with that of field G-dwarfs (Duquennoy & Mayor 1991), to within 10 % or less (e.g. Prosser et al. 1994, Bouvier et al. 1997, Duchêne et al. 1999). This result does not seem to depend upon the age of the cluster between $``$1 Myr and 600 Myr. One possible exception might be the Hyades cluster which appears to host more tight (5–50 AU) binaries than the field population does, though statistical uncertainties could account for this slight excess (Patience et al. 1998). It hence appears that the fraction of binaries in very young clusters, at an age of about a million years, is already similar to that measured for field dwarfs and that it does not significantly evolve as the clusters dynamically relax on the ZAMS.

### 2.2. Orbital period distribution

While the overall binary frequency integrated over a range of semi-major axes already provides clues to the properties of binary populations in clusters, one would ideally wish to compare the distribution of their orbital periods. One is thus led to convert the observed angular separations into orbital periods assuming an average system mass and applying a statistical correction factor to take into account projection effects.

Over a limited range of orbital periods, the distribution of semi-major axes in cluster binaries again appears to be consistent with that of field G-dwarf binaries (Duquennoy & Mayor, 1991). The comparison between cluster and field binaries can be further investigated by considering not only visual binaries but spectroscopic binaries as well. This has been done for ZAMS clusters, where the low-mass stars have already been searched for short-period companions. Within statistical uncertainties, the results reinforce the similarity between cluster and field binaries, as shown for the Pleiades in Figure 1.

In order to improve the statistics, we combined the results obtained for several ZAMS clusters. This yields the orbital period distribution shown in Figure 2 for visual binaries. While the distribution of period in cluster binaries agrees with that of field binaries over the orbital range $`10^{4.4}`$$`10^{6.4}`$days, a paucity of wide (sep $``$ 300 AU) cluster binaries appears relative to the field at the 2.5 $`\sigma `$ level. Such a deficit cannot be related to a detection bias, since the distant companions of wide binaries are the most easily detected. The uncompleteness correction which was applied here affects only the first orbital bin of the distribution.

This difference between field and cluster binaries is reinforced by the evidence that the orbital period distribution peaks at a smaller value in clusters than in the field (Patience et al. 2000). The mode of the distribution is 5 AU for cluster binaries and 40 AU for field binaries, the former being estimated under the assumption that all clusters, including the Hyades, can be combined into a single sample. Furthermore, Scally et al. (1999) have shown that there are no very wide binaries (sep $``$ 2000 AU) in the Trapezium cluster.

Hence, even though clusters at various evolutionary stages appear to host a roughly similar fraction of low-mass binary systems, which itself is similar to the field binary fraction, the distribution of orbital periods appears to exhibit small but significant differences between clusters and field binary populations.

### 2.3. Mass ratio distribution

The mass of the primaries and secondaries of binary systems on the ZAMS can be reliably estimated from their near-infrared flux using main sequence mass-luminosity relationships. In order to increase the statistical significance of the derived mass ratio distribution, the combined results for three clusters observed with similar set-ups are shown in Figure 3. The observed distribution appears to be essentially flat over the range $`q=M_b/M_a=0.1`$–1. This agrees with the $`q`$-distribution derived by Duquennoy & Mayor (1991) for low-mass field binaries if the limited sensitivity of our observations, which prevents us from detecting tight systems with low mass ratios, is taken into account.

Contrary to main sequence binaries, young low-mass systems frequently exhibit significant continuum excesses at near-infrared wavelengths and this may affect the derivation of mass ratios from near-IR flux ratios. However, it is usually assumed on a statistical basis that both stars in a young binary system exhibit similar excesses, and that the overall flux ratio distribution can be safely transformed into a mass ratio distribution. Baraffe et al. (1998)’s models were used to get the mass-luminosity relationship relevant for such young stars. Both in the Taurus star-forming region (Leinert et al. 1993) and in the Trapezium cluster (Prosser et al. 1994), the observed $`q`$-distribution is consistent with being flat down to $`q=0.1`$, and therefore not different from the ZAMS and field estimates.

Current results thus suggest that the mass ratio distribution is similar for all low-mass binary populations studied so far, independent of their evolutionary status and environmental conditions (field or clusters). No peak is found in the distribution towards $`q1`$ nor towards low $`q`$. Furthermore, the ZAMS clusters mass ratio distribution does not seem to vary with separation for visual binaries (Patience et al. 2000).

### 2.4. Massive binaries

Studies conducted so far have mostly focussed on the low-mass binary population of clusters. More recently, however, the investigation of massive binaries has started, with the additional difficulty of a huge contrast between the bright primary and the low-mass companion which is a major limitation of current massive binary surveys.

The recent HST study of thirteen OB stars in the Trapezium cluster by Preibisch et al. (1999) already revealed eight visual companions with separations ranging from 15 to 450 AU. The binary frequency of high-mass stars in this cluster thus appears significantly higher than that of their low-mass counterparts (Prosser et al. 1994), and the difference is further enhanced when completeness corrections are taken into account. Similarly, an adaptive optics survey of about fifty OB stars in the distant cluster NGC 6611 revealed eight companions in the separation range 200–2000 AU (Duchêne 2000), indicative of a high binary fraction compared to G-type field dwarfs if the orbital period distribution is the same for high-mass and low-mass binaries. A study of young, mostly isolated, Herbig AeBe stars conducted by Bouvier & Corporon (2000) yields an uncorrected binary frequency which is significantly higher than that of lower mass Main Sequence stars for both Be and Ae PMS stars.

The mass-ratio distribution of massive binaries can be estimated from the flux ratios measured by imaging surveys. The $`q`$-distribution derived for massive binaries with B0–B7 primaries is shown in Figure 4. The sample includes members of the Trapezium cluster and NGC 6611 as well as isolated Herbig Be stars, with a semi-major axis ranging from 100 to 2000 AU. The mass ratio distribution is strongly peaked towards small values: only about 25 % of the binaries have $`q>0.5`$. Patience et al. (2000) obtained a similar trend for A and F stars in various ZAMS clusters. This distribution is dramatically different from that of low-mass binaries, indicating that the primary’s mass is a relevant parameter when investigating the formation and evolution of binary systems.

## 3. Interpretation

The observed properties of young binary systems hold clues to their formation mechanism and their early evolution. In the following, we consider both classes of models, binary formation and dynamical evolution of clusters, and compare their predictions with observations.

### 3.1. Dynamical evolution of stellar clusters

The fact that all clusters studied so far exhibit a similar fraction of visual binaries regardless of their age rules out a long-term decrease of the binary frequency past the first million years after gravitational collapse. In addition, since similarly-aged loose associations and clusters display different binary frequencies, it is tempting to conclude that environmental effects play a key role in fixing the binary content of a star forming region.

$`N`$-body simulations of the dynamical evolution of stellar clusters were ran by Kroupa et al. (1999) and compared with observations of the Trapezium cluster. They found that close encounters between cluster members could efficiently disrupt binaries in less than a million years while the number of binaries would remain constant afterwards. This results from the fact that the number of encounters falls off dramatically after one relaxation time. If the initial low-mass binary fraction is of the order of 100 % at the end of the protostellar collapse, such a scenario would naturally explain the high binary frequency observed in T Tauri associations: in these regions, the stellar density is several orders of magnitudes too low to generate numerous enough close-by interactions so as to impact on the initial binary content (Clarke & Pringle 1991). In contrast, frequent gravitational encounters in dense clusters would rapidly erode the initially high binary fraction and lead to a binary frequency similar to that of the field after a million years of evolution or so.

The observed properties of binary populations described in Section 2 further support the hypothesis that dynamical evolution of clusters may largely drive the binary properties. In dynamical models, wide binaries are preferentially disrupted, because their binding energy is small compared to the average kinetic energy of the cluster (see the review by P. Kroupa in this volume). This prediction agrees with the observed distribution of orbital periods for binaries in ZAMS clusters which exhibits a paucity of wide systems, as shown in Figure 5. The higher binary fraction found for massive stars is also in qualitative agreement with this model since the potential well of massive primaries is deeper than that of solar-mass stars, making it more difficult to disrupt more massive binaries. The “softening boundary” which delineates hard and soft systems, the latter being prone to disruption, is thus moved towards wider separations.

A number of issues remain however to be addressed before accepting this paradigm. For instance, dynamical models predict that low-$`q`$ systems are preferentially disrupted, introducing a non-negligible modification in the shape of this distribution. This constrasts with the derivation of similar mass ratio distributions for low-mass binaries in clusters and associations. Quantitatively, Kroupa (1995) estimated that 70 % of the initial binaries must have a mass ratio $`q<0.4`$ in order to account for the observed distribution of field binaries. This is unlike the $`q`$-distribution observed for T Tauri associations. Another puzzling result is the fact that all clusters seem to have very similar binary frequency and orbital period distributions. Hence, binary properties appear to be weakly dependent on small differences in initial conditions. Dyamical models have focussed on reproducing the binary population of the Trapezium cluster and it remains to be seen whether the same results would be obtained for slightly less dense clusters such as IC 348.

According to these models, wide binaries found in the solar neighbourhood cannot have formed in clusters, since they would have been rapidly disrupted. One thus has to assume that they were born in loose associations that have subsequently dissipated in the field. However, the field binary population does not appear to be a simple combination of cluster and association members: a weighted mixture of these 2 populations that reproduces the observed fraction of tight and wide binaries also yields a large excess of systems in the 5–50 AU separation range compared to field binaries. The issue of the efficiency of star formation in clusters as opposed to associations is still open.

### 3.2. Binary formation models

If the dynamical evolution of clusters is responsible for the observed binary properties, then the formation process itself must yield a large initial binary fraction. Furthermore, the masses of the components must be taken at random from the initial mass function, as suggested by Kroupa et al. (1995). We note that the observations of high-mass binaries are qualitatively consistent with random pairing, with a limited number of $`q>0.5`$ pairs, though this result is still preliminary due to the small sample size. We will now discuss whether binary formation mechanisms do predict such properties for initial binaries, prior to cluster’s evolution.

The formation of a stellar cluster probably results from the fragmentation of the parent molecular cloud into a large number of individual cores. In this environment, wide binary systems can form either by capture or through the fragmentation of a small collapsing core. Models considering capture in small-$`N`$ ($`N<10`$) subclusters predict that the two most massive members are usually associated in a binary system, while the other stars are ejected (McDonald & Clarke 1993). This leads to a mass ratio distribution which either peaks towards unity or is roughly flat, depending on the value of $`N`$. However, taking into account the dissipative effect of protostellar disks during close encounters leads to a secondary mass distribution that mimics the IMF (McDonald & Clarke 1995). Furthermore, the model predicts a higher binary frequency for more massive stars. Similar simulations made by Sterzik & Durisen (1998) additionally predict that the mass ratio distribution is more heavily weighted towards low-$`q`$ values in higher mass systems. All these predictions are consistent with the observations, indicating that stars in clusters may form in compact subclusters. However, such a scenario does not end up with a 100 % initial binary frequency, since some single stars are ejected from each subcluster.

Capture during the relaxation of small-$`N`$ subgroups appears to be the model that best fits current observations, and especially the random pairing of secondary masses. The alternative model of scale-free fragmentation (Clarke 1996) predicts that binary properties do not depend on the primary mass, which conflicts with the observations. Binary properties could also result, at least partly, from the accretion of residual material onto two gravitationally bound seeds (Bate & Bonnell 1997, Bate 2000). However, these models predict a larger number of $`q1`$ binaries for tighter systems as well as for more massive binaries, which is not verified in the observational trends. Finally, binaries formed through disk fragmentation should preferentially have small separations and low mass ratios, a prediction which contrasts with the apparent lack of a dependence of the mass ratio distribution upon separation.

## 4. Summary

The advent of high-angular resolution imaging techniques have offered the opportunity to systematically survey binaries in young clusters, thus complementing spectroscopic studies. We reviewed some statistical properties of visual binaries such as their frequency, orbital period distribution and mass ratio distribution; both low-mass and high-mass stars were considered. Observational results can be summarized as follows: the low-mass binaries in cluster have similar properties as Main Sequence (field) binaries, except for a significant deficit of systems with a separation larger than 300 AU. High-mass stars host many more companions than do solar-type stars, and most massive binaries have a small mass ratio. These observational constraints are consistent with the results of numerical simulations describing the early dynamical evolution of stellar clusters provided that the initial binary frequency is of the order of 100 % for low-mass stars. Disk-assisted capture in small subclusters is the formation model that best approaches the initial conditions required to account for the observations. In both classes of models, however, several issues still need to be investigated and additional simulations over a wider range of physical conditions have to be run before these scenarios can be considered as offering an acceptable new paradigm for the star and binary formation process. In parallel, future observations of young embedded clusters will help to derive the properties of protobinaries at an even earlier stage of evolution than discussed in this review, which will undoubtly lead to new clues to their formation.

## References

Baraffe, I., Chabrier, G., Allard, F., Hauschildt, P. H., 1998, A&A, 337, 403

Bate, M. R. & Bonnell, I. A. 1997, MNRAS, 285, 33

Bate, M. R. 2000, MNRAS, 314, 33

Bouvier, J., Rigaut, F., & Nadeau, D. 1997, A&A, 323, 139

Bouvier, J. & Corporon, P. 2000, IAU Symposium 200: The Formation of Binary Stars, eds. Matthieu & Zinnecker

Clarke, C. J., & Pringle, J. E., 1991, MNRAS, 249, 584

Clarke, C. J., 1996, MNRAS, 283, 353

Duchêne, G., 1999, A&A, 341, 547

Duchêne, G., Bouvier, J., & Simon, T. 1999, A&A, 343, 831

Duchêne, G. 2000, Ph.D. thesis, University of Grenoble

Duquennoy, A., & Mayor, M. 1991, A&A, 248, 485

Durisen, R. H., & Sterzik, M. F. 1994, A&A, 286, 84

Ghez, A. M., Neugebauer, G., & Matthews, K. 1993, AJ, 106, 2005

Kroupa, P. 1995, MNRAS, 277, 1491

Kroupa, P., Tout, C. A., & Gilmore, G. 1993, MNRAS, 262, 545

Kroupa, P., Petr, M. G., & McCaughrean, M. J. 1999, New Astronomy, 4, 495

Leinert, C., Zinnecker, H., Weitzel, N., Christou, J., Ridgway, S. T., Jameson, R., Haas, M., & Lenzen, R. 1993, A&A, 278, 129

McDonald, J. M. & Clarke, C. J. 1993, MNRAS, 262, 800

McDonald, J. M. & Clarke, C. J. 1995, MNRAS, 275, 671

Mermilliod, J.-C., Rosvick, J. M., Duquennoy, A., & Mayor, M. 1992, A&A, 265, 513

Patience, J., Ghez, A. M., Reid, I. N., Weinberger, A. J., & Matthews, K. 1998, AJ, 115, 1972

Patience, J., Ghez, A. M., & Reid, I. N. 2000, in preparation

Prosser, C. F., Stauffer, J. R., Hartmann, L., Soderblom, D. R., Jones, B. F., Werner, M. W., & McCaughrean, M. J. 1994, ApJ, 421, 517

Preibisch, T., Balega, Y., Hofmann, K., Weigelt, G., & Zinnecker, H., 1999, New Astronomy, 4, 351

Scaly, A., Clarke, C. J., & McCaughrean, M. J. 1999, MNRAS, 306, 253

Sterzik, M. F., & Durisen, R. H. 1998, A&A, 339, 95


format: text

text:
# 1 Introduction

## 1 Introduction

Although variability plays a key role in constraining the models of the central engine of QSOs, little is known about the physical origin of luminosity variations. The most diverse variability mechanisms have been proposed in the past, including supernovae explosions (Aretxaga et al. 1997), instabilities in the accretion disk (Kawaguchi et al. 1998), and gravitational lensing due to intervening matter (Hawkins 1993). So far, most of the information about the characteristics of variability derives from the statistical analysis of single band light curves of magnitude limited QSO samples. The correlation of variability with either intrinsic luminosity ($`v`$-$`L`$) or redshift ($`v`$-$`z`$) is affected by the strong correlation between luminosity and redshift ($`L`$-$`z`$), present in these samples. The results of these analyses also depend on the specific variability index adopted, as shown by Giallongo et al. (1991, GTV), who found a positive $`v`$-$`z`$ correlation through a variability index based on the rest-frame structure function, later confirmed by Cristiani et al. (1996). The increase of variability with the rest-frame observing frequency found by Di Clemente et al. (1996), supports the suggestion of GTV that QSOs at high redshift appear more variable, on average, since they are observed in a higher rest frame frequency. The dependence of variability on frequency is associated with the hardening of the spectral energy distribution (SED) during the bright phases (Cutri et al. 1985, Kinney et al. 1991, Edelson et al. 1990). Trèvese et al. (1999) have shown that a hardening of the SED in the bright phase occurs, on average, in the statistical, magnitude limited, sample of AGNs of the SA 57 (Trèvese et al. 1989, 1994, Bershady et al. 1998). They also found that the slope $`\alpha `$ and its variations associated with the luminosity changes, are consistent on average with temperature changes of a black body.

## 2 Analysis and results

In a recent paper, Giveon et al. (1999) provide the light curves of a subsample of PG quasars consisting of 42 nearby and bright QSOs ($`z<0.4`$, $`B<16`$ mag) which were monitored for 7 years in the Johnson Cousins B and R bands, with the 1 m telescope of the Wise Observatory, with a typical sampling interval of 39 days. We define the structure function as (see Di Clemente et al. 1996) $`S=\sqrt{\pi /2|m(t+\tau )m(t)|^2\sigma _n^2}`$, where $`m(t)`$ is either the B or the R magnitude, t is the rest-frame time, $`\tau `$ is the time lag between the the observations, $`\sigma _n`$ is the relevant r.m.s. noise and the angular brackets indicate the average taken over all the pairs of observations separated by a time interval $`\tau \pm \mathrm{\Delta }\tau `$. The value of $`\mathrm{\Delta }\tau `$ is the result of a trade-off between time resolution and statistical uncertainty. In the following we define, for each QSO and for each band, two different variability indexes $`v_1`$ and $`v_2`$ defined as $`S`$, computed for $`\tau =0.3\pm 0.09`$ yr and for $`\tau =2\pm 0.6`$ yr respectively. The average values of $`v_{1B}`$, $`v_{1R}`$, $`v_{2B}`$, $`v_{2R}`$ over the ensemble of 42 QSOs, versus the relevant average rest-frame frequency, are reported in Figure 1a, which is adapted from Di Clemente et al. (1996).

The new points are consistent with the general trend and confirm the value of the slope $`S_1/\mathrm{log}\nu _{rest}=S_1/\mathrm{log}(1+z)0.250.3`$ which accounts for the $`v`$-$`z`$ correlation found by GTV. Although this correlation is not present for the whole sample (Giveon et al. 1999), due to the small redshift range ($`z<0.4`$), the $`v`$-$`z`$ correlation appears (see Fig. 1b), if we restrict to an absolute magnitude bin $`23.5<M_B<22.5`$, to reduce the effect of the (positive) $`L`$-$`z`$ and (negative) $`v`$-$`L`$ correlations. The correlation coefficient $`r_{v,z}=0.39`$ is only marginally significant ($`P(>r)=0.09`$), since the number of QSOs in the selected magnitude bin is only 19. Therefore, the $`v`$-$`z`$ correlation and its interpretation in terms of spectral variability are confirmed, despite the poor statistics and the narrow redshift range. For each QSO we computed the instantaneous spectral slope $`\alpha (t)\mathrm{log}(f_{\nu _B}/f_{\nu _R})/\mathrm{log}(\nu _B/\nu _R)`$, regarding as synchronous the observations within a time lag of 9 hours. In Figure 2 $`\alpha (t)`$ is shown for each QSO at each observing time as a function of the monochromatic luminosity $`L_\nu `$ at a fixed rest-frame frequency. The regression lines $`\alpha `$ vs. $`\mathrm{log}L_\nu `$ are also reported for each QSO. With a few exceptions, a general trend of increasing $`\alpha `$ for increasing $`\mathrm{log}L_\nu `$ appears. Moreover, the average $`\alpha `$ of each QSO is also positively correlated with the average luminosity. This fact suggests that both the increase of $`\alpha `$ during luminosity variations and the average increase of $`\alpha `$ from faint to bright QSO have similar physical origins, like e.g. the increase of temperature of the emitting gas. For comparison in Figure 2 are reported the $`\alpha (t)`$ versus $`\mathrm{log}L_\nu (t)`$ lines computed for three emitting black bodies of different area.

This admittedly oversimplified model accounts qualitatively for both the intra-QSO and the inter-QSO $`\alpha (t)\mathrm{log}L_\nu (t)`$ correlation. Future comparison with different emission models will indicate their consistency or inconsistency with SED variability data.


text:
# The formation of Encke meteoroids and dust trail 1footnote 11footnote 1Based on observations with ISO, an ESA project with instruments funded by ESA Member States (especially the PI countries: France, Germany, the Netherlands and the United Kingdom) with the participation of ISAS and NASA.

## 1 Introduction

The discovery of cometary dust trails by the Infrared Astronomical Satellite (Davies et al. 1984, Sykes et al. 1986) revealed that comets emit far more dust than expected based on visual wavelength observations (Sykes and Walker 1992). The trails consist of large, dark particles (in the millimeter to centimeter size range), and they were found to represent the principal mass loss mechanism for the comets with which they were associated. From the IRAS observations, it was inferred that dust trails were a phenomenon common to short-period comets in general, and that the ‘dirty snowball’ of Whipple (1950) was more accurately a ‘frozen mudball’ (Sykes 1993). Even the classical ‘gassy’ comet P/Encke was observed to have a trail, and it was shown that the dust to gas mass ratio of the material ejected from the comet was 3.5, as opposed to the canonical 0.1 to 1 (Sykes and Walker 1992).

We will use the terms ‘meteoroid’ and ‘dust particle’ synonymously in this paper, but it should be evident from our results that the ‘dust particles’ we are talking about are larger than those often referred to by observers of other comets. In fact, the ‘dust particles’ are the same size as the particles that give rise to meteors when they enter the Earth’s atmosphere (Ceplecha 1998). Thus when we describe the ‘Encke dust trail,’ we are describing the origin of meteoroids from comet Encke. Comet Encke is the parent body of the Taurid meteor stream complex, so we infer that the ‘Encke dust trail’ is the set of recently-produced particles that would become a meteor stream if their orbit crossed the Earth’s orbit (as the Taurid meteoroids’ orbits now do). Meteoroids and dust particles are ejected from comets due to pressure from sublimating ices, and after being accelerated to some terminal velocity ($`v_{ej}`$), the principal forces on optically large particles are solar gravity and radiation pressure (parameterized by the force ratio $`\beta =F_{rad}/F_{grav}`$). The width of the dust trails perpendicular to the orbital plane, observed by IRAS, constrained $`v_{ej}`$ and $`\beta `$. For P/Encke, the trail spanned a $`90^{}`$ of mean anomaly, or about 2 AU of physical length. The length of a dust trail ahead of and behind the comet allows an estimate of the ages of the trails, which are typically many orbital periods of the comets (Sykes et al. 1990, Sykes and Walker 1992). However, the relatively low angular resolution of IRAS limited the information that could be easily extracted from the data.

Frequent opportunities exist to study comets from ground-based telescopes, but low albedos and optical depths of $`\tau 5\times 10^9`$ make trails extremely difficult to observe at visible wavelengths. In the thermal infrared they are only available from space-based platforms. The apparition of P/Encke in 1997 provided a special opportunity to study the dynamics and distributions of various size ranges of particles ejected from the nucleus, because we could see the comet from a favorable viewing geometry above its orbital plane. The favorable geometry, increased sensitivity, higher angular resolution of the Infrared Space Observatory camera relative to IRAS, and proximity of Encke to the Earth, made the 1997 apparition of Encke a prime opportunity to advance our understanding of the nature and origin of cometary dust.

## 2 ISO Observations

We observed comet Encke using the Infrared Space Observatory (ISO; Kessler et al. 1996) on July 14, 1997, after Encke had passed perihelion on May 23, 1997. At the time of observation Encke was at a distance $`\mathrm{\Delta }=0.25`$ AU from the Earth and $`R=1.15`$ AU from the Sun. It was observed from a vantage point well above the orbital plane, with the angle between the line of sight and the orbital plane being $`35^{}`$. The phase angle (Sun-Encke-Earth) was $`53^{}`$ for this observation. Figure 2 shows the locations of Encke and Earth in their orbits at the time of observation. For reference, Encke’s orbit has semimajor axis 2.21 AU, eccentricity 0.85, and inclination $`11.9^{}`$.

For our observations, we used the mid-infrared camera ISOCAM (Cesarsky et al. 1996) with the $`6^{\prime \prime }`$ pixel-field-of-view lens and the widest filter (LW10, nominal wavelength 11.5 $`\mu `$m, range 8–15 $`\mu `$m). ISOCAM consists of a $`32\times 32`$ detector array, so in this configuration the instantaneous field of view is $`3.2^{}\times 3.2^{}`$. The ideal diffraction limit is $`5^{\prime \prime }`$, approximately the same as the pixel size, so a typical point source will have a full width at half maximum of less than 2 pixels.

Owing to its rapid motion ($`7^{}`$/hr) and large angular size ($`>20^{}`$), comet Encke was a challenging target for ISO observations. We devised a strategy that allowed us to cover a large area with high sensitivity, without allowing the comet to smear or the stars to trail. To do this, we made a series of pointed observations, taking 11 exposures of 2.1 sec duration at each position. Encke’s motion during the 23 sec spent at each position is less than half a pixel, so each set of 11 frames can be coadded without smearing the comet. To build a mosaic large enough to show extended features, we pointed the telescope in a $`21\times 11`$ raster, with step sizes of $`90^{\prime \prime }\times 76^{\prime \prime }`$, oriented perpendicular to the comet’s motion. The raster leg length was chosen to be long enough to cover the desired area, yet fast enough that the comet could not move more than half of the ISOCAM field of view during a raster leg. This ensures that we saw the comet on two raster legs. The spacing between raster legs was chosen to be equal to the motion of the comet during an individual raster leg. This is the optimal strategy for covering the widest possible area in comet-centered coordinates, while maintaining enough overlap between frames to allow accurate removal of the pixel-to-pixel gain variations.

Because the ISO operation system did not allow raster observations of moving targets, we negotiated with the ISO team for an observing date, close to Encke’s passage near the Earth, that fit within the ISO schedule and viewing constraints. For each of the possible observing dates, we calculated the position and orientation of motion of Encke for that observing date, and we specified a window for scheduling. The observations spanned the period between 00:07 and 02:03 UT. If the observation had been scheduled earlier or later, part or all of the desired field would have been missed. (Unfortunately, the other comet in this observing program, 81P/Wild 2, was missed entirely.) Because the Encke observation took place exactly at the time for which we calculated the ephemeris, Encke is centered in the final map. In celestial coordinates, the map covers $`15.7^{}\times 33^{}`$, and in cometocentric coordinates, the map covers $`30^{}\times 33^{}`$.

The set of 2541 exposures were reduced as follows. First, cosmic rays were identified using a running median on the time history of the brightness seen by each pixel. For each cosmic ray, the affected pixel was flagged for 4 exposures beginning with the cosmic ray event; this helps remove pixels affected by transient gain variations induced by the cosmic ray. Then the dark current was subtracted from each exposure using a library image of the dark current from the routine calibration observations. Then the time history of each pixel was corrected for the transient gain response, by deconvolving the the time series from the response curve, which consists of a fast response to 60% of a brightness change followed by a slower response to the remaining 40% with a time constant $`\tau =320/I`$ sec for a brightness $`I`$ MJy sr-1; the transient corrections make only small changes to the image because the comet and zodiacal light were very bright for this observation (cf. Coulais and Abergel 1999). The 10 exposures at each position for which the telescope was not slewing were coadded to give 242 images. The pixel-to-pixel gain response and the optical vignetting were corrected by dividing each frame by a flat-field image, constructed from the median of all images excluding images taken within $`10^{}`$ of the comet nucleus.

Figure 2 shows a mosaic of the ISOCAM images in the rest frame of comet Encke. Each of the 242 images was shifted to match the comet’s motion relative to its position halfway through the observation. Therefore, this image simulates a ‘snapshot’ of the comet on July 14, 1997 at 01:00 UT. In Fig. 2, stars appear as a pair of point sources spread along the direction of the comet’s motion (roughly right to left). The spacing between the star images depends on the interval between crossings of the star during the raster (which was performed with legs of alternating direction).

Several structures can be distinguished in our ISOCAM observation of comet Encke. Figure 2 shows an unsharp-masked version of the ISOCAM image, with several features labeled. The structures include a coma, consisting of a lobe extending to the south and a spike extending to the northwest, and a thin, linear trail that stretches across the image, running directly through the nucleus from the upper right (following the comet) to the lower left (leading the comet). The trail is brighter following the comet than leading, and there is a steeper gradient of brightness with distance leading the comet than following the comet. The vertical profile of the trail is also somewhat different for the portions of the trail leading and following the comet. Figure 2 shows two perpendicular brightness profiles of the dust trail. The two cuts are at angular separations of $`\pm 12^{}`$ from the nucleus, corresponding to a separation from the nucleus by an orbital mean anomaly of $`\pm 0.014^{}`$ and a physical distance of approximately $`\pm 1.5\times 10^5`$ km. The profile leading the comet is consistent with a single Gaussian component, while the profile following the comet is very different. In addition to a nearly Gaussian core (with width similar to that of the trail ahead of the comet), there is a wider, underlying component, as well as significant emission from the tail, which curves over the trail behind the comet. For the purpose of this paper, we fit the trail as the sum of two Gaussian components, labeled ‘core’ and ‘skirt’ in Fig. 2. Table 1 shows the parameters of Gaussian fits to the brightness profiles. The optical depth was calculated for the peak brightness through the trail: $`\tau =I_\nu /B_\nu (T)`$, where $`B_\nu `$ is the Planck function. We assume a temperature $`T=270`$ K, based on the time of the ISOPHOT spectrophotometry of the coma (Lisse et al. 2000). For reference, $`B_\nu (270\mathrm{K})=2.6\times 10^8`$ MJy sr-1 at 11.7 $`\mu `$m (the center of the response curve for the LW10 filter). As discussed below, the particles in the dust trail are likely to be larger than those in the coma and tail. But even the coma and tail particles appear to be large enough that they emit as greybodies: their temperature is only 4% warmer than expected for rapidly-rotating grains larger than 10 $`\mu `$m radius, and they have a grey emissivity even out to 100 $`\mu `$m wavelength (Lisse et al. 2000). Very large grains may preserve a day-night temperature difference across their surface, which could elevate the temperature of the trail particles somewhat. Sykes and Walker (1992) determined that Encke trail particles had temperatures $`<6`$% in excess of a greybody. The temperature excess is likely to be less than 15 K, meaning our optical depths are unlikely to be wrong by more than 25%.

Much of the ISOCAM image is filled with emission from Encke dust. To check whether the background emission in our image is partially due to a very extended distribution of dust spilling over the edges of the image, we use the COBE/DIRBE zodiacal light model (Kelsall et al. 1998) and the ISOCAM filter transmission (Pérault et al. 1994) to estimate the zodiacal light in our filter and on our observing date to be $`13\pm 2`$ MJy sr-1. Because our observation was at a relatively low galactic latitude, we also expect some emission from the diffuse interstellar medium. Using the DIRBE-corrected IRAS dust template of Schlegel et al. (1998) to get the interstellar surface brightness at 100 $`\mu `$m wavelength, and a nominal interstellar dust spectrum (Reach and Boulanger 1998), we estimate the brightness of interstellar dust in our image to be $`1.2`$ MJy sr-1. Thus the total predicted background for our image is $`14\pm 2`$ MJy sr-1. The faint regions in our image, near the bottom, have a brightness of $`13.7\pm 2`$ MJy sr-1, which is entirely consistent with the predicted background. Therefore, it appears that our image does extend far enough to catch the southern edge of the emission from Encke dust. The dust trail clearly extends off the image toward the west, and the coma clearly extends off the image toward the northwest.

## 3 Nucleus and inner coma

Figure 3 shows a close-up of the inner coma made using data from a single raster leg, with logarithmic contours. The individual pixel that pointed at the nucleus was saturated, meaning that the brightness of the $`1100\times 1100`$ km area centered on the nucleus cannot be determined. An approximate fit to the coma surface brightness profile between $`6^{\prime \prime }`$ (1100 km) and $`6^{}`$ (65,000 km) is

$$I_\nu 70(\theta /30^{\prime \prime })^1\mathrm{MJy}\mathrm{sr}^1,$$
(1)
where $`\theta `$ is angular distance from the nucleus. Note that the coma is highly asymmetric at $`\theta >30^{\prime \prime }`$, and this equation only gives the azimuthally-averaged falloff of brightness with distance from the nucleus. (Also, the coma profile flattens close to the nucleus; a better fit to the coma from $`1^{}`$ (11,000 km) outward has a radial power-law exponent of -1.15.) A background of 13.7 MJy sr-1 (the brightness in the edges of the image) due to zodiacal light and diffuse interstellar medium, was subtracted from the image (see previous section). Using the dust temperature of 270 K, measured by Lisse et al. (2000) with broad-band 5–100 $`\mu `$m ISOPHOT observations taken only 4 days later than our observations, the projected optical depth of the coma is

$$\tau 2.6\times 10^7(\theta /30^{\prime \prime })^1.$$
(2)
Because the profile is flatter than $`\theta ^2`$, the integrated flux from the comet is not strongly dominated by the inner coma. The dust coma of Encke was very bright, with an integrated flux within $`6^{}`$ (65,000 km) diameter of 40 Jy.

The radius of the nucleus of P/Encke was recently measured by Fernandez et al. (2000) using ESO and HST observations; they found $`R=2.4\pm 0.3`$ km. Using our observations combined with the ISOPHOT observations, we can make an independent estimate of the size of the nucleus. The single ISOCAM pixel containing the nucleus was saturated, which means that the flux in that pixel exceeded 2 Jy. If we extrapolate the coma surface brightness inward using eq. 1, the coma contributes about 0.9 Jy to the central pixel. The remaining $`>1.1`$ Jy/pixel could be due to the nucleus or to enhanced inner coma brightness (relative to our extrapolation). For a point source, about half of the flux falls in the central pixel (based on integrating the ISOCAM point spread function; Okumura 1998), so the lower limit to Encke’s flux is 2 Jy. An upper bound on the nuclear flux is set by an unsaturated measurement of the flux in a $`3^{}`$ aperture using ISOPHOT (Lisse et al. 2000). Interpolating between the ISOPHOT 10 $`\mu `$m and 12.8 $`\mu `$m fluxes, the flux at our wavelength would be 19 Jy. Integrating our ISOCAM image over a $`3^{}`$ diameter aperture centered on the nucleus (but excluding the central, saturated pixel), the flux is 21 Jy. From controlled observations of calibration stars, the relative calibration between ISOCAM and ISOPHOT is accurate to $`10`$% (Moneti 1998); also, the observations themselves have $`10`$% uncertainties. Therefore the nuclear flux from the ISOPHOT-ISOCAM difference is $`2\pm 3`$ Jy, which we can interpret as an upper limit of 4 Jy. The approximate confidence interval of nuclear fluxes, 2–4 Jy, constrains the radius of P/Encke as follows. If the nucleus were a absorbing, rapidly-rotating sphere with a subsolar latitude of zero, its equatorial temperature would be 260 K, and the radius must be in the range 1.5–3 km. If the nucleus were a slow-rotating sphere with a single temperature on the sunlit hemisphere and zero on the other, the sunlit hemisphere temperature would be 307 K, and the radius would be in the range 1–2 km. Thus the radius of P/Encke must be in the range 1–3 km. This compares well with the direct radar measurement of Encke’s radius, which yielded a confidence interval 0.5–3.8 km (Kamoun et al. 1982). Our limit is similar to previous limits determined from infrared and optical observations (Campins 1988; Gehrz et al. 1989; Luu and Jewitt 1990; Fernandez 2000).

## 4 Dynamics of the Encke Dust Trail

### 4.1 Dynamical principles

Space-based infrared observations of short-period comets reveal a more extensive dust environment than what is easily observable from ground-based telescopes. To understand the origin of this extended dust environment, we need to compare it to predictions for particles of a given ejection velocity ($`v_{ej}`$) and size (parameterized by $`\beta `$). A dynamical analysis of cometary particles yields only $`\beta `$ and $`v_{ej}`$, from which may be inferred the size (and mass). Assuming that the particles are spherical with radius $`a`$ ($`\mu `$m) and mass density of $`\rho `$ (g cm-3), and they are large compared to the wavelength of sunlight (i.e. $`a1`$), then

$$\beta =\frac{K}{\rho a},$$
(3)
where $`K=5.7\times 10^5`$ g cm-2 (Burns et al. 1979). Note that the particles size for a given $`\beta `$ is rather uncertain, because it depends on the particles’ mass density and geometry. Collected interplanetary dust particles and chondritic meteorites appear rather solid, with $`\rho 13`$ (Corrigan et al. 1997), but cometary dust may be rather porous and irregular (cf. Lien 1992).

If the ejection velocity is assumed to be zero, then lines of constant $`\beta `$ (syndynes) can be used to determine the range of $`\beta `$s ejected by a comet. Figure 4.1 shows the syndynes overlaid on a greyscale version of the ISOCAM image. This figure can be used to define the canonical three main parts of the comet: coma, tail, and trail. Once the dust is ejected it will end up in either the tail or the trail, and the ‘dust coma’ is only a transition region. We consider the dust coma to be a roughly spherical region around the nucleus whose radius is equal to the distance at which sunward emitted dust with intermediate values of $`\beta `$ ‘turns around’ in the reference frame of the comet due to radiation pressure. For our observations of comet Encke, the ‘dust coma’ is the peanut-shaped region of intense emission with a diameter of about $`9^{}`$ ($`10^5`$ km) centered on the nucleus. Dust with larger values of $`\beta `$ would form a comet tail; particles ejected at perihelion will have unbound orbits if $`\beta >(1e)/2`$ (Burns et al. 1979), which means $`\beta >0.075`$ for Encke. Based on Figure 4.1 and the more detailed dynamical models presented below, there is no evidence for a dust tail due to small particles (less than 200 $`\mu `$m in radius) in comet Encke. The syndynes for particles with very small $`\beta <10^4`$ are not shown in Fig. 4.1 because they crowd into the dust trail, as is found in the simulations below.

### 4.2 Numerical simulations

To better understand the observed shapes of the coma and trail, we created a series of simulations using a Monte Carlo dynamical model (cf. Campins et al. 1990). Our simulations assume a differential size distribution

$$\frac{dn}{d\beta }\beta ^k,$$
(4)
ranging from $`10^2`$ to $`10^6`$, with $`k=1`$ as the initial guess. The ejection velocity depends on size and heliocentric distance as

$$v_{ej}=1.35\beta ^{0.5}R^{0.5}\mathrm{km}\mathrm{s}^1$$
(5)
(cf. Whipple 1951). Lisse et al. (1998) found that this velocity distribution reproduced the COBE observations of C/Austin, C/Levy, and C/Okazaki-Levy-Rudenko. Over the range of $`\beta `$s used in the simulations, our ejection velocities (eq. 5) reproduce Sekanina’s ejection velocities (1988a,b) to within 10%.

See f7.jpg

Figure 4.2 shows the results of three simulations. Each panel shows the number of particles per pixel weighted by the time a particle spends in that pixel. The images were smoothed to $`15^{\prime \prime }`$ resolution to improve their appearance. The first simulation, shown in panel Figure 4.2a, assumes isotropic emission over the sunlit hemisphere. The dust production rate scales with heliocentric distance as $`\dot{M_d}R^2`$, with dust ejection beginning at aphelion ($`R=4.1`$ AU; October 1995). The isotropic emission model in Figure 4.2a fails to reproduce any of the coma structure observed in our ISOCAM image. This confirms the previous work that demonstrated that Encke’s mass loss is anisotropic, based on the non-gravitational forces and orientation of the visible fan-like tail (Sekanina 1986, 1988a,b), that showed that the visual appearance of Encke’s fan-like dust tail could be understood in terms of a vents on the northern and southern hemispheres of a rotating and precessing nucleus. The isotropic emission model with smooth $`\dot{M}`$ also fails to reproduce the observed time-variation of the OH production rate (A’Hearn et al. 1985). An active area on the northern hemisphere was suggested by Fanale and Salvail (1984) to explain observed variations in OH production, and by Newburn and Spinrad (1985) to explain an anomalous decrease in dust loading of the coma near perihelion.

The second simulation, shown in Figure 4.2b, assumes emission from a single jet at a latitude of $`15^{}`$ with the north pole pointing at at RA=102.8, Dec=32.6. This polar orientation corresponds to an obliquity of nearly $`90^{}`$, so that Encke’s pole is in its orbital plane. The rate of particle ejection from the jet was proportional to the cosine of the zenith angle of the Sun. A rotational period of 24 hours is used in this analysis. This period differs from the 15 hours found by Luu and Jewitt (1990). However, our results are completely insensitive to rotational periods of less than a few days, because any rotational structure is smoothed out by variations in $`\beta `$ and the low ejection velocities. The single-jet, rotating model in Figure 4.2b includes dust emission occurring since the most recent aphelion (October 1995). This model reproduces the gross features of the ISO image: the southern lobe of material in the coma, the ‘spike’ to the NW, and the orientation and width of the trail (illustrated in Fig. 2). However, the extent of the trail is not reproduced. Figure 4.2c shows the same model as Figure 4.2b, but with dust ejection beginning at the aphelion of the previous apparition (June 1992). Clearly, reproducing all of the observed structures requires not only that the dust is ejected anisotropically, but also that the dust is emitted over multiple orbits.

To better understand the types of particles that give rise to the dust coma and trail, we repeated the rotating jet model used in Figure 4.2b for particles in four separate decades of $`\beta `$: $`10^2`$$`10^3`$; $`10^3`$$`10^4`$; $`10^4`$$`10^5`$; and $`10^5`$$`10^6`$. Figure 4.2 shows the simulated images for particles in these four size ranges. The ISOCAM image must be produced by particles with $`\beta <10^3`$, because the higher-$`\beta `$ particles create a diffuse coma that does not look at all like the one we observed. The bright lobe of emission extending $`2^{}`$ to the southwest of the nucleus, and the ‘spike’ to the northwest, are due to the projection of the expanding material from the jet. The ‘gap’ in the images for intermediate-$`\beta `$ particles is created when the jet is pointed away from the Sun, near perihelion, due to the tilt of Encke’s rotational pole—that is, when it becomes ‘winter’ for the jet (late April 1997). The ‘turn-around’ seen coming from the ‘spike’ is due to the first material ejected just after the comet turned on again during this most recent apparition (mid June 1997). Figure 2 shows the locations of Encke and Earth for these events. All of the dust observed in these simulations eventually ends up in the dust trail of comet Encke.

See f8.jpg

Only the oldest and largest particles from our simulation of the present apparition of Encke are close to joining the dust trail. Thus, we conclude that the core of the dust trail contains mostly particles emitted during previous apparitions. To reproduce the entire dust trail as seen in our image, we ran the same model as Figure 4.2c for the same 4 particle size ranges as Figure 4.2. Dust production begins in June 1992, at the aphelion of the previous apparition. Figure 4.2 shows a number of important features. First, even including particles from the previous apparition, particles with $`\beta >10^3`$ do not contribute significantly to the ISOCAM image. Second, Figure 4.2 shows that the particles that extend off the northwest edge of the ISOCAM image most likely have intermediate size: $`10^4>\beta >10^5`$. Third, the core of the dust trail is produced by very large particles. The width (FWHM) of the simulated trail in Fig. 4.2c is $`4^{}`$, while the width of the trail in Fig. 4.2d is $`2^{}`$. Fig. 4.2d is a very good match to the ‘core’ of the dust trail in Figure 2 and Table 1; therefore, the core of the dust trail is produced by particles with $`\beta <10^5`$. The ‘skirt’ of the dust trail, which appears only trailing the nucleus, is produced by particles with a range of sizes, with a typical value $`\beta 10^4`$.

See f9.jpg

The ISOCAM image shows, for the first time, cometary trail material before it actually enters into the trail. In principle, this allows us to determine the trail mass injection rate unequivocally once we have the correct rotational and jet model. Unfortunately, the model presented in Figure 4.2 is not unique. Figure 4.2 shows the simulation of Figure 4.2c along with a simulation based on Sekanina’s rotational model (two jets with latitudes of $`75^{}`$ and $`+55^{}`$ and a pole position of RA$`=206^{}`$, Dec$`=+3^{}`$ which corresponds to an obliquity of $`70^{}`$). In the Sekanina model, the NW ‘spike’ is created solely by the $`75^{}`$ jet, which turns off just before perihelion, and the southern lobe is created solely by the $`+55^{}`$ jet. In our model both features are created by the same jet. Sekanina’s analysis assumed that the bisection of the observed visible emission fan indicates the instantaneous projection of the spin vector. This assumption is only correct if the dust is relatively ‘young.’ If the fan is created by particles with $`\beta <10^3`$, as suggested by our results, then the youngest dust particle would be over 1 month old at the outer edge of the wedge. Since most of the observations were made around perihelion, there can exist a large change in the projection of the spin axis over the course of 1 month, in which case the instantaneous projection of the spin axis will not coincide with the bisection of the observed dust wedge. Therefore it is very difficult to make an accurate measurement of the orientation of the rotational pole. Neither our model nor Sekanina’s model are perfectly matched to the observations. However, given the relatively large phase space of possible solutions, it is interesting that such simple models can match the ISOCAM observations as well as they do. Taken together, our work and that of Sekanina show that a high obliquity is required to interpret both infrared and optical observations of dust from Encke.

See f10.jpg

### 4.3 Model results

Using the dynamical simulations, we constrained the particle sizes through their force ratio $`\beta `$ (Eq. 3). Each particle was assigned a surface brightness equal to a Planck function at $`T=270`$ K (see above) multiplied by the solid angle of the particle as seen from Earth. Then the production rate and $`\beta `$ distribution were adjusted in order to reproduce the observed surface brightness. The production rate was adjusted by making simulations for particles ejected over short intervals and determining how their contributions add to make the total image. The initial guess for the production rate was determined by choosing a narrow range of $`\beta `$ and calculating the surface brightness of dust ejected over a small range of heliocentric distance. Dividing the observed image by the initial-guess simulation yields an estimate of the heliocentric dust production rate. The $`\beta `$ distribution was also adjusted by making simulations for size ranges and determining how their relative contributions add to make the total image. The initial guess for the $`\beta `$ distribution was obtained by choosing a heliocentric distance where there is a good separation between the effects of $`\beta `$ and production rate, then calculating a simulated image for a small range of $`\beta `$ at that heliocentric distance. These adjustments were made step by step, and our final solution is not unique; for example, the orientation of the pole and the jet are correlated with the production rate. However, the salient features of the model are robust. The simulations were normalized by the total number of particles emitted over the simulation multiplied by the production rate and $`\beta `$ distribution functions. Because the model assumes a single jet which turns on and off with season, the resultant production rate is a complex function of heliocentric distance.

We found the $`\beta `$-distribution index from equation 4 to have a value $`1.3<k<1.6`$ for the range $`10^7<\beta <10^3`$, with no particles with larger $`\beta `$. This is not a unique solution for the size distribution, and a much more exhaustive modeling effort would be able to determine the size distribution in more detail—including the higher-$`\beta `$ particles. For our size distribution, both the mass and the surface area are in relatively large particles. The mass is in the very largest particles and the area is in particles with $`\beta 10^3`$. The ISOCAM observations and our model clearly shows a paucity of the smallest particles (with high $`\beta `$). This is not due to physical optics alone, which predicts that particles with $`2\pi a/\lambda <1`$ will not emit efficiently. For our observation $`\lambda =11.5`$ $`\mu `$m, so only particles smaller than 2 $`\mu `$m radius, or $`\beta 0.3`$, will emit inefficiently due to physical optics. But we find a paucity of particles with $`\beta >10^3`$, which corresponds to $`a<600`$ $`\mu `$m. If physical optics were the reason for the observed lack of high-$`\beta `$ particles, then the relation between size and $`\beta `$ would have to be modified drastically. Alternatively, there could be a real lower limit to the size of particles ejected from Encke. In either case, our inferred $`\beta `$ distribution differs significantly from those found for other comets such as P/Halley, for which the mass distribution is dominated by large particles while the surface area (and the infrared emission) is dominated by much smaller particles (McDonnell et al. 1991). We find that the infrared emission from P/Encke is dominated relatively large particles, so the size distribution for Encke must be more shallow than that for Halley.

Using the numerical simulations that best match the ISOCAM image, we find that the total mass emitted since aphelion is within the range 2–6$`\times 10^{13}`$ g (assuming particles with mass density $`\rho =1`$ g cm-3). The range in values comes from the different $`\beta `$ distributions that produced acceptable matches to the observed surface brightness of the comet. The model we present here is not unique, and it could be improved by optimizing the parameter values to match the observations in detail. Therefore, our derived mass loss rate remains uncertain, although the factor of 3 range is our best estimate.

To estimate the contribution to the trail from the previous apparition, we ran the numerical model for approximately 1.5 orbits, starting from the aphelion before the previous apparition. In the core of the trail, the model brightness is only about 10% of the observed brightness. This suggests the trail is at least 10 orbits old, which is consistent with the trail age estimated from its total extent (Sykes and Walker 1992). The mass of the trail is therefore at least 2–6 $`\times 10^{14}`$ g, comparable to the estimate of $`1.6\times 10^{14}`$ g estimated from the IRAS observations (Sykes and Walker 1992).

As a check, a model-independent estimate for the mass of the portion of the dust trail present in the ISO image can be measured from the observables as follows. If we assume the trail is a cylinder oriented perpendicular to the line of sight, with cross-sectional angular diameter $`W`$ and length $`L`$, the distance from the Earth is $`\mathrm{\Delta }`$, particle density is $`\rho `$, and the optical depth through the trail is $`\tau `$, then

$$M_{trail}=\frac{\pi }{3}\mathrm{\Delta }WL\rho a\tau ,$$
(6)
where $`a`$ is the particle radius in the dust trail. This equation assumes that both the mass and emitting surface area are produced by large particles of effective size $`a`$. Because only the combination $`\rho a`$ determines the trail mass, and the dynamics constrain $`\beta (\rho a)^1`$, the masses that we derive here are independent of the mass density $`\rho `$ of the particles, and the masses are independent of whether the particles are fluffy or compact. We can set a generous upper limit $`\beta <10^4`$ by comparing the dynamical simulations to the morphology of the dust trail and coma. Using the brightness profile fit from Table 1 we find a trail mass within our image $`M_{trail}1\times 10^{12}`$ g, with comparable amounts coming from the core and skirt. This trail mass is a drastic underestimate for at least three reasons. First, the lower limit to $`\beta `$ determined from the dynamical simulations is far form the typical values that create the observed dust trail; the particles form previous apparitions are in the core of the trail and can have $`\beta `$s that are more than an order of magnitude lower. Second, even if $`\beta =10^4`$ particles were to dominate the brightness, larger particles would also exist and could dominate the total dust mass if the size distribution is at least as shallow as was observed for P/Halley (McDonnell et al. 1991). Third, it is clear that the trail extends off the edge of the ISOCAM image. Using the IRAS observations, Sykes and Walker (1992) found the Encke dust trail extends fully $`90^{}`$ of mean anomaly relative to the nucleus. At the time of the IRAS observation, this corresponded to a length of 2 AU, and at the time of our observation, the same mean anomaly would correspond to a length of 5.5 AU along the orbit. Our ISO observation covered only 0.002 AU of the comet’s orbit. The actual mass of the Encke dust trail is at least an order of magnitude larger that our limit just due to this affect. These arguments suggest the trail mass is at least $`10^{14}`$ g, which is consistent with the results of our numerical simulations and provides a sanity check.

## 5 Impact hazard in Encke’s coma

In November 2003, the CONTOUR spacecraft is scheduled to fly within 200 km of Encke (Veverka 1999). At this time the comet will be at heliocentric and geocentric distances similar to those when our ISOCAM observations were made.222Comet distances for the 2003 CONTOUR encounter were calculated using the JPL Horizons program, http://ssd.jpl.nasa.gov/horizons.html. Using a color temperature of 270 K, as measured measured by ISOPHOT (Lisse et al. 2000), the optical depth at an impact parameter 200 km is $`\tau 7\times 10^6`$. For spherical particles of a single size, the mass of particles encountered by a spacecraft (with area $`A`$) is

$$M_1=\frac{4KA\tau }{3\beta }.$$
(7)
The encountered mass is relatively insensitive to assumptions about the particle properties; in particular, it is independent of the mass density (hence, fluffiness) of the particles. The number of particles encountered by the spacecraft is

$$N=\frac{A\tau \rho ^2\beta ^2}{\pi K^2},$$
(8)
which is sensitive to particle density as well as $`\beta `$. Using the upper limit $`\beta <10^3`$ for coma particles, we find that a 10 m2 spacecraft will encounter $`>0.05`$ g of cometary particles during its traverse. The expected number of particles encountered with $`\beta <10^3`$ is $`N10^2`$.

Using the size distribution from the numerical simulations (eq. 4), the encountered mass increases relative to this lower limit to

$$M_{enc}=M_1\frac{k1}{2k}f^{k2}\frac{1x^{k1}}{1x^{2k}},$$
(9)
where $`x=\beta _{min}/10^3`$ and $`\beta _{min}`$ is the $`\beta `$ of the largest particle (which dominates the mass). For example, for $`\beta _{min}=10^6`$ and $`k=1.5`$, the mass of cometary particles encountered is $`M_{enc}1`$ g. Larger particles exist, but their spatial distribution is more confined. Thus, to make an accurate assessment of the impact hazard, one would need to compare the spacecraft trajetory and the particle distribution as a function of both particle size and position in the coma. For comparison, the dust detectors on the Giotto spaceprobe passing through the coma of comet Halley in 1986 detected only 0.01 g, though an estimated total mass of $`0.15`$ g penetrated the dust shields (McDonnell et al. 1986), at a distance of 2,200 km. Also, several very large particles (1–50 mg each) must have impacted the Giotto spacecraft to cause the observed attitude shifts (Curdt and Keller 1990).

In the dust trail core, there are relatively fewer particles than in the coma, but they are larger. At high flyby velocities of tens of kilometers per second, such large particles near the comet could represent a significant hazard to a spacecraft. Using the trail optical depth of $`\tau 1.4\times 10^8`$ and an upper-limit to $`\beta <10^4`$, there is $`<0.1`$% chance of a spacecraft with an area of 10 m2 hitting a trail particle, but each particle is potentially catastrophic.

## 6 Mass loss and meteoroid production of comet Encke

The total mass of particles ejected by Encke in its 1997 apparition is 2–6 $`\times 10^{13}`$ g. Dividing the mass of the nucleus (using $`R=2.4`$ km \[Fernandez et al. 2000\] and a nuclear mass density $`\rho _N`$ in g cm-3) by the mass production we infer for the 1997 apparition, the fraction of Encke’s mass that was lost in the 1997 apparition is

$$\frac{\mathrm{\Delta }M}{M}=(0.11)\times 10^3\rho _N^1.$$
(10)
For comparison, in his classic papers, Whipple (1950, 1951) found that the non-gravitational forces observed on comet Encke’s orbit could be produced by mass loss at a (recent) rate of $`\mathrm{\Delta }M/M=0.002`$, with certain assumptions about the nature of the mass loss. Our inferred mass loss rate is only smaller by a factor of few, which is a remarkable confirmation of the general principle of Whipple’s comet model, in which the mass loss rate depends on the velocity of the emitted particles and their angular distribution (which determines the direction of the net momentum carried away by the emitted material). A modern version of this model takes into account that the mass loss occurs primarily from small regions (‘jets’) whose illumination depends on season (Sekanina 1986, 1988a). Further, a modern model should take into account that most of the mass of comet Encke is lost in solids as opposed to gas, so there is a component of the non-gravitational force due to the solids that could be somewhat different from that of the gas. If the dust is lost through the same jets as the gas, as appears to be the case from an analysis of the inner coma of P/Halley (Reitsema 1989), then the only effect of our new dust-to-gas estimate is to increase the amount of mass lost per unit ice sublimation. The non-gravitational acceleration, which is the observable quantity, depends only on $`\mathrm{\Delta }M/M`$. Thus the main effect of increasing mass loss rate by a factor of 10–30 times the gas mass loss rate is to increase the mass (and density) of the nucleus as inferred from non-gravitational forces.

The dust mass loss rate we infer from the dust trail observations is indeed significantly larger than the gas mass loss rate. A’Hearn et al. (1985) summarized observations the OH production along Encke’s orbit, and Sekanina (1986) provide a convenient fit. The total mass of sublimating H2O is $`\mathrm{\Delta }M_{gas}=2.1\times 10^{12}`$ g/orbit. Thus the dust-to-gas ratio for comet Encke’s mass loss is $`\mathrm{\Delta }M/\mathrm{\Delta }M_{gas}10`$–30. Clearly our estimate is much larger than the dust-to-gas mass ratios normally inferred for comets based on optical observations. The fraction of the mass loss occurring in non-volatile form is 91–97%. For comparison, Whipple (1955) used 20% as the fraction of mass lost in non-volatile, meteoric, material. Based on the P/Halley encounter data from Giotto, the dust-to-gas mass ratio was found to be $`2`$ (McDonnell et al. 1991). Based on the IRAS dust trail observations, Sykes & Walker (1992) found a dust-to-gas ratio of 3.5 for Encke. In fact their dust mass loss loss rate is $`3\times 10^{13}`$ g/orbit, consistent with the result of this paper, and the gas mass loss rate was taken from an average of historical visual wavelength observations (Kresák and Kresáková 1990).

If Encke continues to lose mass at this rate, its remaining lifetime before shedding essentially its entire mass is 3000–10,000 $`\rho _N`$ yr. Since Encke is presently in a dynamically stable orbit, it appears that mass loss will be its cause of death unless the mass loss is somehow quenched in the next few millenia. If the mass is lost as a uniform layer over the surface, then the fractional decrease in surface area of the nucleus per orbit is $`2\mathrm{\Delta }M/3M`$, which would lead to a reduction in nuclear magnitude by 0.0002–0.0006 $`\rho _N^1`$ mag/orbit. This decrease is significantly less than the average 0.01 mag/orbit secular decrease in comet brightnesses inferred by Kresák and Kresáková (1990), and it is far smaller than the aging of 0.09 mag/orbit estimated by Ferrin and Gil (1988). These conflicts should not be surprising, as the coma brightness is likely a nonlinear function of the nuclear surface area, and the derivation of fading from historical observations is exceptionally difficult. As another application of the mass loss rate, we can estimate the decrease in the nuclear radius with time. Again assuming the mass is lost as a sheet from the upper crust of the nucleus, the fractional decrease in nuclear radius per orbit is $`\mathrm{\Delta }M/3M`$, leading to a change of radius of 0.3-0.9 m/orbit. If Encke is losing mass like Halley, via jets rather than erosion of the entire surface layer, then the jets must either be very wide or they must evolve from place to place on the surface, lest they become so deep that sunlight could not illuminate them.

## 7 Conclusions

The ISOCAM image of the near-nucleus region of P/Encke reveals a widespread region with dust having two principal spatial components: coma and trail. There is no evidence of a tail composed of particles with $`\beta >10^3`$ (or radius $`a<500/\rho `$ $`\mu `$m, with the density $`\rho `$ in g cm-3). The dust structures that are observed within the ISOCAM image require a high obliquity rotational axis and one or more jets. The core of the observed dust trail is dominated by material ejected during previous apparitions of comet Encke; particles from the current apparition particles have not had enough time to spread into the trail.

The dust trail near the nucleus (distinguished from the trail well away from the nucleus as observed by IRAS) appears to have a broad component and a narrow component which tracks the comet’s orbit. The components are ‘fed’ by the southern lobe of the coma, which consists mainly of particles with $`10^5<\beta <10^3`$ (500 $`\mu `$m–5 cm) and a NW ‘spike.’ Particles with $`\beta =10^4`$ (5 mm) appear to be the delimiter between broad and narrow trail components, with all particles larger than 5 mm ($`\beta <10^4`$) ending up in the narrow portion of the trail. The broad wings of the trail are composed of particles with $`10^3<\beta <10^4`$. This result is consistent with the maximum beta of $`\beta 10^3`$ found in P/Tempel 2 dust trail (Sykes et al. 1990) suggesting that the size distributions of large particles may be similar from comet to comet.

The dynamical result that P/Encke’s dust mass loss is dominated by the emission of large particles is consistent with the non-detection of a silicate feature from ISOPHOT observations taken 4 days later (Lisse et al. 2000). In addition, little or no 10 $`\mu `$m silicate feature was detected during the 1987 apparition, when the inner coma of the comet was observed at its brightest, close to perihelion (Gehrz et al. 1989). The lack of a silicate feature, however, does not mean that silicate particles are necessarily absent. Rather, Gehrz et al. conclude that the particles are large enough ($`a>10`$ $`\mu `$m) that their emission spectrum does not contain bright 10 $`\mu `$m (or 20 $`\mu `$m) spectral features. The broad-band spectrophotometry both from the 1987 (Gehrz et al. 1989) and the 1997 (Lisse et al. 2000) apparitions indicates that the temperature of the grains in the coma is similar to that expected for rapidly rotating greybodies. Because the coma emission spectrum appears similar to that of a greybody out to wavelengths $`>100`$ $`\mu `$m, Lisse et al. (2000) conclude that the particles must be larger than 15 $`\mu `$m in radius. Thus, dynamical and spectroscopic evidence all point strongly to a coma dominated by at least mid-sized grains. This distinguishes Encke from new comets such as C/Hale-Bopp, which evidences fragmentation of coma dust and generating a significant high-$`\beta `$ particle population. The dominance of large particles is even more extreme for P/Encke than it is for P/Halley (McDonnell et al. 1991; Fulle et al. 1995) and P/Grigg-Skjellerup (McDonnell et al. 1993), for which in situ observations showed the total mass is dominated by large particles, with mass $`>10^4`$ g or $`a>200\rho ^{1/3}`$ $`\mu `$m.

Comparing the numerical simulations to the observations, we find the total mass of particles ejected by Encke in its 1997 apparition is 2–6 $`\times 10^{13}`$ g. Most of the particles from the 1997 apparition are in the extended dust coma around the comet, and all of these particles eventually end up in the dust trail. Most of the dust trail in our image consists of particles from previous apparitions. Comparing the observed surface brightness of the dust trail to a numerical simulation for 1.5 orbits, and considering the total length of the dust trail (Sykes and Walker 1992), we infer that the trail is at least 10 orbits old. Therefore the total mass of the trail is at least 2–6 $`\times 10^{14}`$ g. Dividing the mass of the nucleus (using $`R=2.4`$ km \[Fernandez et al. 2000\] and a nuclear mass density $`\rho _N`$ in g cm-3) by the mass production we infer for the 1997 apparition, the remaining lifetime of comet Encke is 3000–10,000 $`\rho _N`$ yr. The decrease in surface area of the nucleus per orbit is approximately (0.7-2)/$`\rho _N`$ % per century, so it will gradually fade. Therefore, we are witnessing the final life stages of Encke, before it eventually disintegrates into a meteoroid stream or becomes inactive.

## 8 Acknowledgments

We thank Russ Walker for encouraging us to observe comet Encke with ISO and for many useful discussions during this project. We thank the ISO operations team led by Timo Prusti for getting the Encke observation scheduled exactly on time. And we thank Mike A’Hearn for his critical reading of the manuscript and his helpful comments. This paper benefited from the detailed comments of two anonymoous referees whom we thank. MVS acknowledges NASA grant NAG5-3359.

## 9 REFERENCES

A’Hearn, M. F., P. V. Birch, P. D. Feldman, and R. L. Millis 1985. Comet Encke: Gas production and lightcurve. Icarus 64, 1–10.

Burns, J. A., P. L. Lamy, and S. Soter 1979. Radiation forces on small particles in the Solar System. Icarus 40, 1–48.

Campins, H. 1988. The anomalous dust production in periodic comet Encke. Icarus 73, 508–515.

Campins, H., R. G. Walker, and D. J. Lien 1990. IRAS images of comet Tempel 2. Icarus 86, 228–235.

Ceplecha, Z. K., J. I. Borovicka, W. G. Elford, D. O. Revelle, R. L. Hawkes, V. Porubcan, and M. Simek 1998. Meteor phenomena and bodies. Space Sci. Rev. 84, 327–471.

Cesarsky, C. J. and 65 colleagues 1996. ISOCAM in flight. Astron. Astrophys 315, L32–L37.

Corrigan, C. M., M. E. Zolensky, J. Dahl, M. Long, K. Weir, C. Sapp, and P. J. Burkett 1997. The porosity and permeability of chondritic meteorites and interplanetary dust particles. Meteorit. Planet. Sci. 32, 509.

Coulais, A., and A. Abergel 2000. Transient correction for the LW-ISOCAM data for low contrasted illumination. Astron. Astrophys. Suppl 141, 533–544.

Curdt, W., and H. U. Keller 1990. Large particles along the Giotto trajectory. Icarus 86, 305–313.

Davies, J. K., S. F. Green, A. J. Meadows, B. C. Stewart, and H. H. Aumann 1984. The IRAS fast-moving object search. Nature 309, 315–319.

Fanale, F. P. and J. R. Salvail 1984. An idealized short-period comet model, surface insolation, H2O flux, dust flux and mantle evolution. Icarus 60, 476–511.

Fulle, M., L. Colangeli, V. Mennella, A. Rotundi, and E. Bussoletti 1995. The sensitivity of the size distribution to the grain dynamics: simulation of the dust flux measured by GIOTTO at P/Halley. Astron. Astrophys. 304, 622–630.

Gehrz, R. D., E. P. Ney, J. Piscetelli, E. Rosenthal, and A. T. Tokunaga 1989. Infrared photometry and spectroscopy of Comet P/Encke 1987. Icarus 80, 280–288.

Kamoun, P. G., D. B. Campbell, S. J. Ostro, G. H. Pettengill, and I. I. Shapiro 1982. Comet Encke—Radar detection of nucleus. Science 216, 293–295.

Kelsall, T., J. L. Weiland, B. A. Franz, W. T. Reach, R. G. Arendt, E. Dwek, H. T. Freudenreich, M. G. Hauser, S. H. Moseley, N. P. Odegard, R. F. Silverberg, and E. L. Wright 1998. The COBE Diffuse Infrared Background Experiment search for the cosmic infrared background. II. Model of the interplanetary dust cloud. ApJ 508, 44–73.

Kessler, M. F., and 10 colleagues 1996. The Infrared Space Observatory (ISO) mission. Astron. Astrophys. 315, L27–L31.

Kresák, L., and M. Kresáková 1990. Secular brightness decrease of periodic comets. Icarus 86, 82-92..

Lien, D. 1992. Numerical simulations of cometary dust. In Asteroids, Comets, Meteors 1991 (A. Harris and E. Bowell, Eds.), pp. 359-362. LPI, Houston.

Liou, J. C., and H. A. Zook 1996. Comets as a source of low eccentricity and low inclination interplanetary dust particles. Icarus 123, 491–502.

Lisse, C. M., M. F. A’Hearn, M. G. Hauser, D. J. Lien, S. H. Moseley, W. T. Reach, and R. F. Silverberg 1998. Infrared observations of comets by COBE. Astrophys. J. 496, 971–991.

Lisse, C. et al. 2000. Infrared observations of Comet Encke, .submitted to Icarus

Luu, J., and D. Jewitt 1990. The nucleus of P/Encke. Icarus 86, 69–81.

McDonnell, J. A. M., P. Lamy, and G. S. Pankiewicz 1991. Physical properties of cometary dust. In Comets in the Post-Halley Era (R. L. Newburn, Jr., M. Neugebauer, and J. Rahe, Eds.), pp. 1043–1073. Kluwer, Dordrecht.

McDonnell, J. A. M., W. M. Alexander, W. M. Burton, E. Bussoletti, D. H. Clark, J. L. Grard, E. Grün, M. S. Hanner, Z. Sekanina, and D. W. Hughes 1986. Dust density and mass distribution near comet Halley from Giotto observations. Nature 321, 338–341.

Moneti, A. 1998. Quantitative cross-calibrations: I. Theory and first results, preprint (European Space Agency).

Newburn, R. L. and H. Spinrad 1985. Spectrophotometry of seventeen comets. II. The continuum. Astron. J. 90, 2591–2608.

Okumura, K. 1998. ISOCAM PSF report, ISO/ESA.Memorandum

Pérault, M., F. X. Désert, A. Abergel, F. Boulanger, Ch. Dupraz, A. Soufflot, C. J. Cesarsky, and L. G. Vigroux 1994. Infrared Space Observatory camera calibration facility and preflight characterization. Opt. Eng. 33, 762–770.

Reach, W. T., and F. Boulanger 1998. Infrared emission from interstellar dust in the local interstellar medium. In The Local Bubble and Beyond (D. Breitschwerdt, M. J. Freyburg and J. Trümper, Eds.), pp. 353–362. Springer, Berlin.

Reitsema, H. J., W. A. Delamere, A. R. Williams, D. C. Boice, W. F. Huebner, and F. L. Whipple 1989. Dust distribution in the inner coma of Comet Halley—Comparison with models. Icarus 81, 31–40.

Schlegel, D. J., D. P. Finkbeiner, and M. Davis 1998. Maps of dust infrared emission for use in estimation of reddening and cosmic microwave background radiation foregrounds. Astrophys. J. 500, 525–553.

Sekanina, Z. 1986. Effects of the law for nongravitational forces on the precession model of Comet Encke. Astron. J. 91, 422–431.

Sekanina, Z. 1988a. Outgassing asymmetry of periodic comet Encke. I - Apparitions 1924-1984. Astron. J. 95, 911–924.

Sekanina, Z. 1988b. Outgassing asymmetry of periodic comet Encke. II - Apparitions 1868-1918 and a study of the nucleus evolution. Astron. J. 96, 1455–1475.

Sykes, M. V. 1993. Great balls of mire. Nature 363, 696–697.

Sykes, M. V., L. A. Lebofsky, D. M. Hunten, and F. J. Low 1986. The discovery of dust trails in the orbits of periodic comets. Science 232, 1115–1117.

Sykes, M. V., D. J. Lien, and R. G. Walker 1990. The Tempel 2 dust trail. Icarus 86, 236–247.

Sykes, M. V., and R. G. Walker 1992. Cometary dust trails. I. Survey. Icarus 95, 180–210.

Veverka, J. 1999. Contour Mission, website.http://www.contour2002.org

Whipple, F. L. 1950. A comet model. I. The acceleration of Comet Encke. Astrophys. J. 111, 375–394.

Whipple, F. L. 1951. A comet model. II. Physical relations for comets and meteors. Astrophys. J. 113, 464–474.

Whipple, F. L. 1955. A comet model. III. The zodiacal light. Astrophys. J. 121, 750–770.