text:
# The populations of comet-like bodies in the Solar system

## 1 Introduction

There is a need to improve the taxonomy for comet-like bodies in the Solar system. First, as a result of recent discoveries, the populations of comet-like bodies in the Solar system are known to be much more extensive than previously thought. The last decade has seen the discovery of more than 600 trans-Neptunian objects, beginning with (15760) 1992 QB1 (Jewitt & Luu 1993), as well as the identification of $``$100 Centaurs (e.g. Scotti 1992) to supplement the early serendipitous finding of (2060) Chiron (Kowal 1979). Moreover, many unusual objects have now been found circulating essentially within the known planetary system, for example (5335) Damocles (e.g. Asher et al. 1994), and the ‘asteroids’ (4015) Wilson-Harrington = 107P/Wilson-Harrington, (7968) Elst-Pizarro = 133P/Elst-Pizarro and (2060) Chiron = 95P/Chiron. These objects blur the traditional clear distinction between comets and asteroids, either on physical or dynamical grounds.

Secondly, the spatial distribution of objects in orbits largely beyond Neptune, but with aphelion distances $`Q`$ much less than the conventional inner edge of the dynamically active Oort cloud ($`Q\mathrm{20\hspace{0.17em}000}`$au, say), is believed to comprise a flattened disc or belt-like structure. These objects have orbits that in some cases extend hundreds of astronomical units (au) from the Sun, and although resembling asteroids in appearance are widely thought to be cometary in composition. This trans-Neptunian zone includes the region referred to as the Kuiper or Edgeworth-Kuiper belt (EKB), and probably merges in the form of an extended trans-Neptunian disc into the inner core of the Oort cometary cloud. Trans-Neptunian objects are a possibly significant source of low-inclination ‘Jupiter-family’ comets (Quinn, Tremaine & Duncan 1990), again muddying the traditional distinction between comets and asteroids.

In fact, the Centaur population appears to represent an important link between trans-Neptunian objects and Jupiter-family comets (Dones, Levison & Duncan 1996, Stern & Campins 1996, Levison & Duncan 1997). Centaur orbits are typically planet-crossing and have relatively short dynamical lifetimes ($``$$`10^6`$ yr). Chiron, which is one of a number of exceptionally large minor bodies with perihelia close to or within the orbit of Saturn, exhibits cometary activity (e.g. Luu & Jewitt 1990) and even has a periodic comet designation: 95P/Chiron. At least one recently discovered comet, namely C/2001 T4 (NEAT), has a very similar Chiron-like Centaur orbit. Taken together, the evidence suggests a picture in which comets and distant minor planets are dynamically reprocessed from one population to another, for example from the Edgeworth-Kuiper belt, through Centaurs to Jupiter-family comets, and so on. Similarly, it is possible for objects to be scattered out of the main asteroid belt and evolve into similar areas. So, it is possible that some Centaurs may be rocky or asteroidal, others may be icy or cometary. This is analagous to the examples given by Fernández, Gallardo & Brunini (2002) of mixing between the Jupiter-family and near-Earth asteroid populations, and again highlights both the difficulty of separating physically distinct populations of objects using purely dynamical characteristics, and the need for a unified dynamical classification scheme capable of describing the full range of observed orbital types.

In parallel with these observational advances, the march of computer power now ensures that the orbits of most Solar system bodies can be routinely integrated for millions of years. Given the wealth of simulation data and the diversity of new discoveries, the taxonomy of Solar system objects assumes great importance. For a classification scheme to be useful, it should allow us to place objects with similar physical or dynamical characteristics in sets, and to examine the bulk statistics of the sets in detail. However, in making sense of the complicated evolutionary pathways followed by objects such as short-period comets and Centaurs, the historical classification scheme seems obsolete. A new taxonomy should clarify the principal dynamical paths followed by different classes of object through interplanetary space, and facilitate the definition of dynamical lifetimes and the flux and transfer probabilities of objects from one dynamical class to another.

Traditionally, comets have been classified as short-period or long-period according to whether their periods of revolution $`P`$ are less than or greater than 200 years. This unphysical division (Levison 1996) has frequently been justified in terms of the time over which comets have been observed with scientific precision (e.g. Weissman 2001). But, the canonical upper limit for ‘periodic’ comets has now been spectacularly broken by the return of 153P/Ikeya-Zhang with an orbital period in excess of 350 years, and also by recent discoveries such as asteroid 2002 RP120 with a period of about 420 yr.

Centaurs are conventionally defined as asteroids or cometary nuclei circulating largely between the orbits of Jupiter and Neptune and usually crossing the orbit of at least one giant planet (e.g. Jedicke & Herron 1997, Larsen et al. 2001). Jedicke & Herron’s definition restricts Centaur semi-major axes to less than that of Neptune (i.e. $`a<30`$au), and therefore allows objects defined as Centaurs to pass significantly beyond Neptune’s orbit into the Edgeworth-Kuiper belt. A straightforward extension of this definition would allow Centaurs to be regarded simply as a continuation of the ‘scattered disc objects’, perhaps even encompassing those resonant objects that are not protected in the long term from close approaches to Neptune. The scattered disc was the name introduced by Duncan & Levison (1997) following the discovery of 1996 TL66 (Luu et al. 1997) to describe a population of relatively high-inclination, high-eccentricity trans-Neptunian objects. The idea that Centaurs and the scattered disc may be drawn from the same underlying population has been advocated recently both by Marsden (1999a) and by Emel’yanenko, Asher & Bailey (2003).

This paper introduces a new classification scheme for comet-like bodies of the Solar system. There is much evidence from previous work (e.g. Kazimirchak-Polonskaya 1972; Kresák 1972, 1980, 1983; Everhart 1977; Rickman & Froeschlé 1988; Manara & Valsecchi 1991; Dones et al. 1996; Levison 1996) that the evolution of comets and Centaurs often takes place under the control of one or another major planet, usually with near-constant aphelion or near-constant perihelion, or via aphelion-perihelion interchanges. This motivates us to consider the planets controlling the perihelion $`q`$ and aphelion $`Q`$ as fundamental to the classification of any comet-like bodies. There is also evidence that low-inclination objects tend to evolve inwards, whereas the reverse is true for high-inclination objects. This suggests that we use the inclination – or better still, the Tisserand parameter $`T_\mathrm{P}`$ with respect to the controlling planet – as a third criterion. The Tisserand parameter is an approximate constant of the motion, at least when the orbital perturbations are dominated by a single planet. Given the set of values ($`q,Q,T_\mathrm{P}`$), our new taxonomy provides an instantaneous classification for any comet-like body in the Solar system.

## 2 An Outline of the Classification Scheme

### 2.1 The zone of control

The Hill radius is defined as

$$R_\mathrm{H}=a\left[\frac{M_\mathrm{P}}{3M_{}}\right]^{1/3},$$
(1)
where $`M_\mathrm{P}`$ is the mass of the planet, $`M_{}`$ is the mass of the Sun and $`a`$ is the semi-major axis of the orbit of the planet about the Sun (e.g. Murray & Dermott 1999). It corresponds to the position of the Lagrange points in the restricted three-body problem (e.g. Arnold et al. 1987) and so marks the largest distance at which the planet may possess a moon. A Solar system object can be classified according to the planets which have the greatest effect at perihelion and at aphelion. For typical eccentricities, a perihelion lying 4 or 5 Hill radii beyond an inner planet is often sufficiently distant for the dynamics to be controlled by the outer planet. So, we conclude that for a planet to control an object, then the perihelion or aphelion must lie closer than a distance of $``$3 Hill radii to the planet’s orbit (see e.g., Charnoz, Thébault & Brahic 2001).

In our classification scheme, each planet’s zone of control extends to $`3k`$ Hill radii, with the parameter $`k1`$. In this way, the Solar system is banded into concentric zones of control of first Jupiter, then Saturn, Uranus and Neptune. Over the next few million years the mean semi-major axes and maximum aphelia of the four major planets Jupiter, Saturn, Uranus and Neptune are (5.2, 5.5), (9.6, 10.4), (19.2, 20.8) and (30.1, 31.0) au respectively (Applegate et al. 1986). Three times their Hill radii are 1.1, 1.3, 1.4 and 2.3 au respectively. These distances provide a rough basis for estimating the respective boundaries between the zones of control. For example, for Jupiter to control the object at perihelion, the perihelion must lie between 4 and 6.6 au approximately. If the aphelion also lies in the range, then we classify the minor body as a J object. For Saturn to control the object at aphelion, then $`6.6<Q<12.0`$ when the object is denoted by JS. We will give detailed examples shortly, but the general idea should now be obvious to the reader.

### 2.2 The Tisserand parameter

The Tisserand parameter is an approximation to the Jacobi constant, which is an exact integral of motion in the circular restricted three-body problem. It is defined as (e.g. Murray & Dermott 1999)

$$T_\mathrm{P}=\frac{a_\mathrm{P}}{a}+2\mathrm{cos}i\sqrt{\frac{a(1e^2)}{a_\mathrm{P}}}$$
(2)
where $`a_\mathrm{P}`$ is the planet’s semi-major axis.

Ignoring the positive orbital eccentricity of the planet, objects with $`T_\mathrm{P}`$ greater than 3 are in principle confined to regions wholly interior or wholly exterior to the planet’s orbit. There is a long history of usage of the Tisserand parameter with respect to Jupiter $`T_\mathrm{J}`$ to divide the cometary families (e.g., Kresák 1972, 1980, 1982, 1983, 1985; Vaghi 1973 a,b; Carusi & Valsecchi 1987; Levison & Duncan 1994; Levison 1996). Fig. 1 shows all short-period comets in the plane of semi-major axis $`a`$ and eccentricity $`e`$. The comets have been colour-coded according to their Tisserand parameter with respect to Jupiter. Motivated by these plots, we propose a four-fold division according to Tisserand parameter as indicated in Table 1.

In our classification scheme, we use the Tisserand parameters with respect to the major planets ($`T_\mathrm{J},T_\mathrm{S},T_\mathrm{U}`$ and $`T_\mathrm{N}`$ with obvious notation) to classify objects whose perihelion lies within the zone of the planet’s control. So, for example, an SUIV object is a minor body whose perihelion is controlled by Saturn and aphelion by Uranus and whose Tisserand parameter with respect to Saturn is larger than approximately 2.8. Here and henceforth, we always denote the Tisserand parameter class as a subscript.

## 3 Comets

### 3.1 Historical introduction

Objects classified as comets display the widest variety of dynamical behaviour in the Solar system. Unfortunately, there is little or no agreement as to the precise definition of different cometary types. Historically, the major cometary division has been between those classified as ‘long-period’ ($`P>200`$ yr) and ‘periodic’ or ‘short-period’ ($`P<200`$yr) comets (e.g. Duncan, Quinn & Tremaine 1988, Nakamura & Kurahashi 1998, Levison 1996, Crovisier 2001). The majority of the short-period group belong to the so-called Jupiter-family, most of which have very short periods indeed (a mean value close to 8 years). For this reason, the periodic comets are usually divided simply on the basis of orbital period into two classes: Jupiter-family comets with $`P<20`$yr (e.g. Fernández 1985, 1994), and Halley-type comets with $`20<P<200`$yr (e.g. Carusi et al. 1987b).

Since 1999, however, one-apparition comets with periods in the range $`30<P<200`$yr have been given a standard cometary ‘C/’ designation, rather than the usual periodic ‘P/’ identifier (Marsden 1999b, Green 2000). The official range for short-period comets therefore now extends to only 30 yr instead of the former 200-year limit, and comets with $`30<P<200`$yr are now described as ‘intermediate period’ (Marsden & Williams 1999). The principal disadvantages of this change are that one-apparition intermediate-period comets cannot be easily identified as periodic simply from their designation, while the phrase ‘intermediate period’ in fact has a rather wide range of previous connotations, having been used to describe objects with periods both in the Halley-type $`20<P<200`$yr range (e.g. Fernández 1980, Fernández & Gallardo 1994) and with much longer periods up to $`10^3`$yr (e.g. Everhart 1974).

In fact, the situation is more confused still. Whipple (1978) appears to have chosen $`P=25`$yr to separate his ‘Class V’ short-period comets (essentially the classical Jupiter-family) from his Class IV intermediate-period group ($`25<P<1000`$yr); while Everhart (e.g. 1972, 1973) adopted 13 years instead of 20 years as a convenient boundary between short-period and intermediate-period comets, the latter being defined to have periods in the range $`13<P<1000`$yr. Others (e.g. Rickman & Froeschlé 1988, Stagg & Bailey 1989, Wetherill 1991) have instead adopted various aphelion distances $`Q`$ in the range 8–10 au to distinguish Jupiter-family comets from Halley types, while a few authors (e.g. Nakamura & Kurahashi 1998) have simply considered all comets with $`P<1000`$yr as ‘periodic’.

Contrary to the view expressed by Fernández (1994), therefore, a classification based on orbital period alone is not very useful. The range of definitions used by different investigators to describe essentially the same cometary sub-groups appears greatly, and unnecessarily, to complicate any comparison of observational or theoretical results obtained by different authors. We hope that our historical summary has convinced the reader that there is a clear need to improve the taxonomy for comet-like bodies in the Solar system!

### 3.2 Taxonomy of comets

Jupiter’s zone of control extends between approximately 4 and 6.6 au. In our classification scheme, a ‘comet’ has a perihelion $`q<4`$ au, and so lies formally outside Jupiter’s zone of control. Nonetheless, Jupiter dominates the lives of all the short-period comets. For example, 1P/Halley has perihelion within the orbit of Venus and aphelion just beyond the orbit of Neptune, yet Jupiter plays the most important role in its evolution (Carusi et al. 1988). We propose a four-fold division of comets with perihelion distances less than approximately 4 au, to include Encke (E), short-period (SP), intermediate-period (I) and long-period (L) types. The division is based on aphelion distance, with – as usual – further sub-division into the Tisserand parameter classes (Table 1) denoted by Roman-numeral subscripts. Further embellishments on this basic scheme could be introduced to accommodate the resonant and non-resonant main-belt asteroids (e.g. the Hildas, in the 3:2 mean-motion resonance with Jupiter), and objects on Earth-crossing or Earth-approaching orbits with periods much less than that of Jupiter, such as near-Earth objects. The approximate boundaries between the E, SP, I and L comet types are listed in Table 2.

Working outwards from the Sun, we first come to the E-type comets. These have aphelion closer to the Sun than approximately 4 au, i.e. $`Q`$ within the main belt. The prototype is 2P/Encke ($`Q4.1`$ au) which lies on an orbit that does not approach closer than a few Hill radii to Jupiter, making it stable within the inner Solar system for periods of time far greater than that of typical comets. (The orbit is not, however, stable for indefinitely long periods, as it is subject to long-term secular perturbations which drive the perihelion into the Sun on a time-scale of $``$100 000 yr; see Levison & Duncan 1994, Farinella et al. 1994, Valsecchi et al. 1995 for details.) Comet 107P/Wilson-Harrington, with a similar orbit, is a largely inactive body also known as asteroid (4015), and noted by Marsden (1992) to be identical to the low-activity comet C/1949 W1 (=1949g) following Bowell’s identification of precovery observations of the asteroid on Palomar Schmidt images. Another object, (2201) Oljato, has displayed outgassing in the past, and although currently classified as an Apollo-type asteroid could well be an evolved Encke-type comet (e.g. Weissman et al. 1989, McFadden et al. 1993). Many similar bodies in this class appear to be bona fide asteroids, rocky bodies perhaps originating via collisions in the main asteroid belt.

Next we come to the SP comets, which include most objects conventionally thought of as Halley-type and Jupiter-family comets. Many authors have used the Tisserand parameter with respect to Jupiter to separate these classes. Levison (1996) suggested that comets with $`T_J<2.0`$ should be defined as Halley-types and those with $`T_J>2.0`$ should be Jupiter-family. However, the precise value of $`T_J`$ separating the two types of behaviour is somewhat arbitrary and may depend on period. For example, Rickman & Froeschlé (1988) and Bailey (1992) identified a value of the order of 2.5 as a better marker of a dynamically significant separation. Rather than worry about the precise boundary between Jupiter-family and Halley-types, we use the four-fold division according to Tisserand parameter introduced in Section 2.2, namely, Class I with $`T_\mathrm{J}<2.0`$, Class II with $`2.0<T_\mathrm{J}<2.5`$, Class III having $`2.5<T_\mathrm{J}<2.8`$ and Class IV having $`T_\mathrm{J}>2.8`$. Very roughly, these can be thought of as corresponding to Halley-types, transition-types, loosely bound Jupiter-family comets and tightly bound Jupiter-family comets respectively. We call all objects with $`q<4`$ au and $`4<Q<35`$ au short-period (SP) objects, and differentiate between objects within the SP class using the Tisserand parameter with the sub-class denoted by a Roman numeral subscript. So, for example, an SPI object corresponds to a Halley-type comet under the conventional classification.

The third class consists of I-type comets, which have perihelion distances less than 4 au and aphelia between approximately 35 and 1000 au. Comets with aphelion distances less than about 35 au frequently librate around high-order resonances with Jupiter (Carusi et al. 1987a,b, 1988; Chambers 1997) and show evidence for systematic secular changes of their orbital elements (Bailey & Emel’yanenko 1996). Recent examples of I-type comets are C/1995 O1 (Hale-Bopp) and 153P/Ikeya-Zhang. They have Tisserand parameters with respect to Jupiter of 0.879 and 0.040 respectively, and so both fall into the the II category.

The fourth class is the L-type comets, which have perihelion distances less than 4 au and aphelia greater than 1000 au. Comets classified under the conventional system as long-period fall into either the L or the I class in our system. L-types – such as C/1910 A1 (the Great January comet) – have made no more than their first few passages through the inner Solar system. By contrast, I-types are expected to have made rather more perihelion passages on average. It is sometimes useful to sub-divide L-types into new and young (e.g. Fernández 1981; Wiegert & Tremaine 1999; Horner & Evans 2002). New comets, with $`a>\mathrm{10\hspace{0.17em}000}`$ au, are first time entrants to the Solar system from the Oort Cloud.

We are of course aware that many comets display activity beyond 4 au from the Sun, and similarly that there is an increasing number of known objects with perihelion distances beyond 4 au that have been classified as active comets (e.g. C/2002 V2 \[LINEAR\] and C/2003 A2 \[Gleason\]) or which have distinctly cometary orbital characteristics (e.g. 2002 VQ94). We seek to develop a coherent classification scheme that is useful for all comet-like bodies in the Solar system.

From now on, as we move outwards in perihelion, we are able to exploit the idea that comet-like bodies lie under the control of the giant planets. After the ‘comets’, the next class of objects we encounter are those whose perihelion lies within the zone of control of Jupiter. In the conventional scheme, these bodies are still classified either as Jupiter-family comets, or as asteroids if they show no outgassing. In our picture, the designations are listed in Table 2. Here, and henceforth, our classification scheme proceeds with the first letter designating the planet controlling the perihelion and the second letter designating the planet controlling the aphelion or the region in which the aphelion lies. The Tisserand parameter is defined with respect to the planet which controls the perihelion. So, for instance, a JNIV object has a perihelion under Jupiter’s control, an aphelion between 22.5 and 33.5 au and so under Neptune’s control, and a Tisserand parameter with respect to Jupiter $`>2.8`$; while a JEIII type has an aphelion between 33.5 and 60 au and so lying in the Edgeworth-Kuiper belt and a Tisserand parameter with respect to Jupiter in the range 2.5 to 2.8.

An example of a JIV object is 29P/Schwassmann-Wachmann 1, which has a near-circular orbit just beyond Jupiter. Its Tisserand parameter is $`T_\mathrm{J}=3.0`$ and so it is not quite Jupiter-crossing. The object itself is unusual in a couple of respects. The low eccentricity of its orbit coupled with its unusually great size ($``$20 km radius) allows it to be followed with ease through its entire orbit. It undergoes regular outbursts, with its visual magnitude rising from its usual 17 (at perihelion) by a number of magnitudes, reaching as high as 10 at its brightest.

An example of a JSIV type object is 39P/Oterma. It was originally discovered in 1942 orbiting in a Hilda-type 3:2 mean-motion resonance with Jupiter, with a low-eccentricity orbit and a period of a mere 8 years. At that time, its orbit lay between those of Mars and Jupiter and it would have been classified as an SPIV comet in our picture (although clearly as a resonant ‘Hilda’ in a more complete scheme). However, it was later noted that it had entered the 3:2 resonance in 1937 from the 2:3 resonance, and returned to the 2:3 in 1963 (e.g., Kazimirchak-Polanskaya 1972). Continuous periods in any one resonance are short owing to the orbit’s unstable nature, but it is noteworthy that $`T_\mathrm{J}3.0`$ was approximately conserved during this evolution. The present orbit lies outside that of Jupiter, with a current perihelion distance of 5.46 au and an aphelion of 9.02 au. This object seems to be a good example of a recent ‘handing down’ of a Centaur from Saturn’s control to Jupiter, with further periods as an SP comet seeming quite probable. Finally, recent discoveries like C/2002 A1 (LINEAR) and C/2002 A2 (LINEAR), with semi-major axes close to that of Uranus, furnish examples of JNIII objects.

### 3.3 Resonances

A key feature of the orbital evolution of Centaurs and many other comet-like objects in the outer Solar system is the tendency for the semi-major axis to lie close to one or another mean-motion resonance with a planet. This behaviour is exemplified by the so-called Plutinos (objects in the 2:3 mean-motion resonance with Neptune, i.e. around $`a39.5`$ au), which appear clearly in Fig. 2, and also by the evolution of 2000 FZ53 shown in the top left plot of Fig. 3. The occurrence or otherwise of a long-lived resonance will clearly affect the transfer probability from one part of the $`(q,Q,T_\mathrm{P})`$ phase space to another and is an important consideration in the dynamical evolution of outer Solar system objects. There are recently discovered objects with semi-major axes close to essentially all the significant Neptunian resonances, from 1:1 \[e.g. 2001 QR322 and 2002 PN34\] to 1:8 \[e.g. (54520) 2000 PJ30\]. However, detailed discussion of the classification and role of such resonances is beyond the scope of this paper.

## 4 Centaurs and trans-Neptunian objects

Table 3 gives the classification scheme for objects under the control of giant planets other than Jupiter. As before, the notation is that the first letter designates the planet controlling the perihelion, and the second letter the planet controlling the aphelion or the region in which the aphelion lies. So, a UN object has a perihelion controlled by Uranus and an aphelion controlled by Neptune, and so on. The letters EK denote objects with perihelia and aphelia lying close to or within the classical Kuiper or Edgeworth-Kuiper belt. Trans-EK belt or T objects, on the other hand, may have much larger aphelion distances, and (perhaps depending on their perihelion distances) may be stable ‘outer’ objects as described by Emel’yanenko et al. (2003), or dynamically active objects similar to those conventionally known as scattered disc objects. Both EK and T classes have perihelia beyond Neptune’s control (33.5 au) but less than the conventional outer edge of the Edgeworth-Kuiper belt ($``$ 60 au). They differ in that the aphelion $`Q`$ is less than approximately 60 au for the former and greater than 60 au for the latter. As for the comets, the information on the zone of control must be supplemented with the Tisserand parameter with respect to the controlling planet to give a complete classification. The Tisserand parameter class is given as a subscript. It is interesting to note the recent discovery of several high-inclination Centaurs, for example 2002 VQ94 and 2002 XU93, which we classify as STI and UTI respectively.

The first Centaur discovered, (2060) Chiron, has a perihelion distance of 8.4 au and an aphelion of 18.8 au. It is classified as an SUIV object, as its Tisserand parameter $`T_\mathrm{S}=2.9`$. This is one of the largest Centaurs and is active at perihelion. Hahn & Bailey (1990) showed that Chiron’s orbit is rapidly evolving and that it may have been a short-period comet in the past, and will probably become one again in the future. The second Centaur discovered, (5145) Pholus, has a perihelion of 8.7 au and an aphelion of 32.1 au, together with a Tisserand parameter $`T_\mathrm{S}`$ of 2.6. It is classified as an SNIII object. Pholus’s orbit is also unstable with a characteristic lifetime of $`10^6`$ yr. As another example, (10199) Chariklo has a perihelion of 13.1 au which is just within Uranus’s control and an aphelion of 18.6 au. Its Tisserand parameter with respect to Uranus is 2.9, making its classification UIV. Chariklo’s orbit is more stable than either Chiron’s or Pholus’s, as the dynamical lifetime is roughly a factor of 10 times longer. Not all orbits in this part of the Solar system are evolving. For example, Holman (1997) found that a few low-eccentricity objects between 24 and 27 au ($``$0.3% of an initial low-eccentricity primordial sample) may survive for 4.5 Gyr without significant dynamical evolution. These objects lie in the NIV class. Others (e.g. 2001 QR322, also classified as NIV) are recognized as long-lived Neptune-Trojans (Marsden 2003), and it is probably only a matter of time before similar examples of Trojans are confirmed around Uranus as well.

Almost all the known Centaurs lie within the Tisserand parameter classes III and IV. There are only a few exceptions – for example, 2000 FZ53 has a Tisserand parameter with respect to its controlling planet (Uranus) of 2.4. However, the nearly isotropic influx from the Oort Cloud must inevitably produce at least a few high-inclination Centaurs. Hence, the other classes are needed, and their apparent emptiness at present is largely a selection effect.

Fig. 2 shows the plane of semi-major axis and eccentricity. The red circles show the current location of the 52 objects listed as Centaurs and Scattered Disk Objects in 2002 June by the Minor Planet Center111http://cfa-www.harvard.edu/iau/lists/Centaurs.html. A further 16 objects lie at semi-major axes too great to be seen on the plot. These objects were chosen from the full list of Centaurs by requiring the observable arc to be greater than 30 days. Objects with shorter arcs are excluded as their orbits often change substantially from the original estimates when the arc is extended. The known population depends both on the true numbers in the various regions and on discovery selection effects. Orbits under the control of Uranus and Neptune are expected to be more long-lived than those under the control of Saturn and Jupiter. However, this is counterbalanced by the fact that objects nearer the Sun are more easily discovered. The question can be studied further by modelling the efficiency of observational surveys (Jedicke & Herron 1997). For the present Centaur population, there is a suggestion in Fig. 2 of a substantial population in the UE category. These may be being scattered down from planet to planet having originated in the Edgeworth-Kuiper belt.

The green circles in Fig. 2 show the current classification of those 373 objects with arcs greater than 30 days belonging to the Edgeworth-Kuiper belt. Again, our data were downloaded in June 2002 from the Minor Planet Center222http://cfa-www.harvard.edu/iau/lists/TNOs.html. The two most populous categories are NE and EK. The Plutinos, a population of objects trapped in the 2:3 mean-motion resonance with Neptune (e.g. Jewitt & Luu 1996), fall into either the NE or EK category depending on eccentricity. The resonance acts to prevent close approaches to Neptune, although it is unlikely that such a mechanism would generally protect all the objects from eventual encounters with Neptune over time-scales comparable to the age of the Solar system. The EK category contains most of what are conventionally called the ‘classical’ Kuiper or Edgeworth-Kuiper belt population, while the T category contains most of what are conventionally called scattered disc objects. These last two classes can be further subdivided (Emel’yanenko et al. 2003).

## 5 Numerical Examples

We now illustrate the usefulness of the classification scheme with some examples drawn from a suite of numerical simulations. Clones of the current Centaur population are integrated in the forward and backward direction for up to 6 Myr using the hybrid integrator in the Mercury program (Chambers 1999). The orbits are integrated under the gravitational effects of the Sun and the four major planets only. The time-step is 120 days. Once objects passed beyond 1000 au, their orbits are no longer followed.

### 5.1 From UEII to NEII object

Fig. 3 shows the evolution of a clone of 2000 FZ53. For most of the first 1 Myr its perihelion is under the control of Uranus and its aphelion lies in the Edgeworth-Kuiper belt. The Tisserand parameter with respect to Uranus is $`T_\mathrm{U}=2.4`$. The clone starts off as an SEII object and moves rapidly into the UEII zone.

During the course of the simulation, the clone is driven gradually towards a stable orbit with perihelion close to Neptune (an NE object). This illustrates one route by which an object can approach the Edgeworth-Kuiper belt from the Centaur region. There are two sharp transitions visible in Fig. 3. The first change, at 2.0 Myr from UN to UE, is caused by a close encounter at aphelion with Neptune; the second, at 2.9 Myr within the UE category, is caused by an encounter at perihelion with Uranus. At the end of the 6 Myr simulation, the object has a perihelion just beyond the orbit of Neptune, a high inclination, and a low eccentricity. The value of its Tisserand parameter with respect to Neptune is $`T_\mathrm{N}=2.4`$. The orbit is similar to those of a number of objects currently residing in the inner regions of the Edgeworth-Kuiper belt, such as 1999 CP133.

The stability of the semi-major axis towards the end of the 2000 FZ53 simulation suggests that it lies in a mean-motion resonance with Neptune. This is borne out by calculation of the orbital period of the clone, which is $``$202 yr, suggesting that the object probably lies in the 4:5 mean-motion resonance with Neptune (period $``$164 yr). It is possible that the orbit will continue to drift outward and stabilize after the end of the simulation, illustrating a rather rare example of outward rather than inward evolution. The ultimate destiny of such an object may therefore be to join the Edgeworth-Kuiper belt or the region beyond it.

### 5.2 From SUIV object to SPIV comet

Fig. 4 shows the evolution of a clone of Chiron (cf. Hahn & Bailey 1990). This object illustrates one of the supply routes for the short-period comets. The clone starts with its perihelion under the control of Saturn, but its aphelion wanders between the control of Saturn, Uranus and Neptune respectively for the first $`2\times 10^4`$ yr. After this, the clone falls under Jupiter’s control through a succession of aphelion-perihelion interchanges. The Tisserand parameter $`T_\mathrm{J}`$ falls below 3 and the clone spends the next part of its lifetime ($`7\times 10^4`$ yr) as a short-period comet before leaving the SP region. This is longer than the typical fading time ($`1.2\times 10^4`$ yr) for Jupiter-family comets found by Levison & Duncan (1997). It is noteworthy that this 200 km diameter object has a significant probability of evolving to an orbit with perihelion distance less than that of the Earth (Hahn & Bailey 1990). This emphasizes the potential significance for the inner solar system of so-called ‘giant’ comets and their disintegration products when time-scales of the order of $`10^5`$ yr or more are considered (e.g., Napier 2001). From the characteristic shape of the trajectory in the rightmost panel of Fig. 4, we see that the evolution proceeds along lines of roughly constant perihelion or aphelion both before and after the object comes under Jupiter’s control. The part of the clone’s life when it is controlled by Jupiter is characterized by it being highly chaotic and showing rapid dynamical evolution.

Note that the clone starts with a Tisserand parameter with respect to Saturn $`T_\mathrm{S}=2.9`$. During its lifetime as a short-period comet, its Tisserand parameter with respect to Jupiter $`T_\mathrm{J}=2.8`$. Just as in our previous example, the Tisserand parameter with respect to the controlling planet is nearly preserved during the whole evolution.

The fact that the clone stays in an area in which it will be active as a comet for such a long time suggests that cometary bodies can be captured into the inner Solar system for long enough to have a reasonable chance of decoupling from Jupiter, whether by non-gravitational forces such as outgassing (cf. Kresák 1982, Harris & Bailey 1998, Asher, Bailey & Steel 2001) or collisional mechanisms. Integrations occasionally provide examples of Encke-type orbits evolving to become Jupiter-family short-period comets, but the reverse evolution by gravitational means is so rare that it has never been seen (Carusi & Valsecchi 1987). It is interesting to note that the stability of an orbit approaching Saturn is lower than one approaching merely Uranus or Neptune. Similarly, an orbit in the vicinity of Jupiter is yet more unstable. This adds weight to our argument that classifying all objects in the trans-Jovian region merely as ‘Centaurs’ is insufficient to explain the full richness of dynamical behaviour in the region.

### 5.3 From UEIV object to SPIV comet

Fig. 5 shows the evolution of a clone of 1998 TF35. It starts off as a UE object with Tisserand parameter $`T_\mathrm{U}=2.9`$. This object is of interest since it passes through a number of different categories (UE, UN, NE, N, U, SU and S) before finally becoming a short-period comet and then suffering ejection from the Solar system. Generally, it moves through the categories sequentially rather than frequently jumping back and forth between classes. The clone begins in an orbit influenced both by Neptune and Uranus, and distant encounters provide gradual perturbations for most of its lifetime. Eventually, encounters with Neptune move the object into a near-circular orbit controlled only by that planet. This period of evolution ends with a close encounter which moves the clone’s perihelion back to the control of Uranus. The perihelion drops further until the object encounters Saturn, at which point it is rapidly handed down to Jupiter’s control, where it becomes a short-period comet for a brief time ($`10^4`$ yr) before being ejected from the Solar system. Whilst a short-period comet, $`T_\mathrm{J}=2.9`$, so there is again excellent conservation of the Tisserand parameter with respect to the controlling planet.

### 5.4 Comets Lagerkvist-Carsenty and Jäger

As a final example, let us consider two objects conventionally classified as short-period Jupiter-family comets, namely P/1997 T3 (Lagerkvist-Carsenty) and P/1998 U3 (Jäger). Numerical integrations by Lagerkvist et al. (2000) showed that both objects suffered recent close encounters with Saturn that drastically changed their orbits. Fig. 6 shows the evolution of the orbits of both comets over the last century. Comet Lagerkvist-Carsenty started out as an SUIV object in January 1900. A close encounter with Saturn in October 1954 transferred its osculating elements through almost instantaneous SN and SE phases into its present JSIV orbit. This proceeded through a perihelion-aphelion interchange, with the old perihelion becoming the new aphelion. By contrast, Comet Jäger started as an SIII object in January 1900 and underwent a direct transition to its present SPII orbit in July 1991. The evolution proceeded at nearly constant aphelion, as can be deduced both from the plot of the orbit and from its evolutionary track in the plane of semi-major axis and eccentricity. In this case, the close approach caused the Tisserand parameter class to change, the controlling planet switching from Saturn to Jupiter.

These orbital integrations provide illustrations of why the classification scheme is useful. The scheme is based around the aphelion and perihelion distances. Close encounters often preserve aphelion or perihelion to lowest order, or a perihelion-aphelion interchange may take place. Very often, the Tisserand parameter with respect to the controlling planet is conserved to a good approximation as well.

## 6 Conclusions

The main aim of the paper is to introduce a new classification system for comet-like bodies. Minor bodies between Saturn and Neptune are often described simply as ‘Centaurs’ and those beyond Neptune simply as ‘Kuiper Belt Objects’. This is not very enlightening as the histories and fates of such bodies may be very different. For example, simulations show that some Centaurs are in long-lived and stable orbits, whereas others are dynamically highly active and will evolve rapidly.

The main focus of this paper has been on the Centaurs and short-period comets. Our proposition is to classify the comet-like objects beyond Jupiter according to the planets that control the evolution of their perihelion and aphelion. For example, an SN object has a perihelion under the control of Saturn and an aphelion under the control of Neptune, while a UE object has a perihelion under the control of Uranus and an aphelion in the Edgeworth-Kuiper belt. This provides 20 dynamically distinct categories of outer Solar system objects in the Jovian and trans-Jovian regions. The evolutionary tracks of bodies like Centaurs often show periods in which the aphelion or perihelion distances are individually rather well conserved, or encounters in which the old perihelion becomes the new aphelion or vice versa (‘perihelion-aphelion interchanges’).

The objects with smaller perihelion distances, which we designate as ‘comets’, are dominated by Jupiter and not by the planets nearest at perihelion or aphelion. In our scheme, comets have perihelion distance less than 4 au. They are subdivided into Encke-type, short-period, intermediate-period and long-period according to aphelion distance. Following a succession of authors beginning with Kresák (1972), we favour further sub-divisions based on the Tisserand parameter. Specifically, we subdivide the categories of comets into: Class I which has $`T_\mathrm{J}<2.0`$, Class II which has $`2.0<T_\mathrm{J}<2.5`$, Class III having $`2.5<T_\mathrm{J}<2.8`$ and Class IV having $`T_\mathrm{J}>2.8`$. Very roughly, these can be thought of as corresponding to Halley-types, transition-types, loosely bound Jupiter-family comets and tightly bound Jupiter-family comets respectively. This idea is then extended to all comet-like bodies by identifying the Tisserand parameter of the planet controlling the perihelion with a subscript.

Given the aphelion and perihelion distance, together with the Tisserand parameter of the planet controlling the perihelion, it is straightforward to find the instantaneous classification of any Solar system object. Our new classification scheme extends the existing taxonomy for comets to cover all comet-like bodies in the Solar system.

## Acknowledgments

This research was supported by the Particle Physics and Astronomy Research Council (JH), the Royal Society (NWE) and the Northern Ireland Department of Culture, Arts and Leisure (MEB, DJA). DJA acknowledges the hospitality of the Japan Spaceguard Association during work on this paper. We thank the referee for very constructive comments.


format: text

text:
# Inverse temperature dependence of the dust submillimeter spectral index

## 1 Introduction

Characterizing dust emissivity properties is one of the major challenges of nowadays submillimeter astronomy. In this spectral range, big grains at thermal equilibrium (see for instance Désert et al. (1990)) dominate the dust emission. Their emission is characterized by a temperature and a spectral dependence of the emissivity which is usually simply modelled by a spectral index.

The temperature, density and opacity of a molecular cloud are key parameters which control the structure and evolution of the clumps, and therefore, star formation. The spectral index ($`\beta `$) of a given dust grain population is directly linked to the internal physical mechanisms and the chemical nature of the grains.

It is generally admitted from Kramers-König relations that 1 is a lower limit for the spectral index (Emerson (1988)). A classical value is 2, following the investigations of Gezari et al. (1973) and the calculations of, for instance, Draine & Lee (1984). This value is particularly invoked for isotropic crystalline grains, and it is thought to be an upper limit (Wooten (1972), Emerson (1988)). For amorphous silicate or graphitic grains, a value of $`\beta =2`$ is also favoured for different physical reasons (Rowan-Robinson (1986), Tielens & Allamandola (1987), Mennella et al. (1995)). However, it is not the case for amorphous carbon, which is thought to have a spectral index equal to 1 (Koike et al. (1980)), as well as aggregates of silicates and graphite in a porous structure (Mathis & Whiffen (1989)). Silicate grains embedded in ice mantles are thought to have spectral indices between 1.5 and 2 (Miyake & Nakagawa (1993), Preibisch et al. (1993), Krügel & Siebenmorgen (1994), Pollack et al. (1994)). In general, amorphous rather than crystalline structure and increase of the grain size (Mannings & Emerson (1994)) are the main reasons to predict submillimeter spectral indices lower than 2. Spectral indices above 2 may exist, according to several laboratory measurements on grain analogs (Koike et al. (1995), Mennella et al. (1995), Agladze et al. (1996), Mennella et al. (1998)). Also, Agladze et al. (1996) and Mennella et al. (1998) showed an anticorrelation existing for some types of grains between the temperature and the emissivity spectral index.

Observations of the diffuse interstellar medium at large scales favour $`\beta `$ around 2 (see for example Boulanger et al. (1996), Dunne & Eales (2001)). In the case of molecular clouds, spectral indices are usually found to be between 1.5 (Walker et al. (1990)) and 2 (e.g. Wright et al. (1992) in Orion). However, low values (0.2-1.4) of the spectral index have been observed in circumstellar environments (Weintraub et al. (1989), Beckwith et al. (1990), Knapp et al. (1993)). This can be attributed to grain growth in dense stellar envelopes. Low indices have also been observed in molecular cloud cores (e.g. Blake et al. (1996)). Spectral indices larger than 2 have also been observed in the millimeter range (Schwartz (1982), Wilson et al. (1995), Kuan et al. (1996)), in particular on the whole sky at large scales by the WMAP satellite (Bennett et al. (2003)).

In this context, submillimeter multi-band imaging of the interstellar medium provides a lot of useful information about dust properties and interstellar medium structure, especially if the dust emission parameters, namely both the temperature and the spectral index, can be properly derived on top of submillimeter intensities.

## 2 Observations

PRONAOS (PROgramme NAtional d’Observations Submillimétriques) is a French balloon-borne submillimeter experiment (Serra et al. (2001), Ristorcelli et al. (1998), Lamarre et al. (1994)). Four bolometers cooled at 0.3 K measure the submillimeter flux with sensitivity to low brightness gradients of about 4 MJy/sr in band 1 (200 $`\mu `$m) and 0.8 MJy/sr in band 4 (580 $`\mu `$m). The effective wavelengths are 200, 260, 360 and 580 $`\mu `$m, and the angular resolutions are 2 in bands 1 and 2, 2.5 in band 3 and 3.5 in band 4. The data that we analyze here were obtained during the second flight of PRONAOS in September 1996, at Fort Sumner, New Mexico. The data processing method, including deconvolution from chopped data, is described in Dupac et al. (2001), and the calibration procedure is detailed in Pajot et al. (2003). This experiment has observed various phases of the interstellar medium, from diffuse clouds in Polaris (Bernard et al. (1999)) and Taurus (Stepnik et al. (2003)) to massive star-forming regions in Orion (Ristorcelli et al. (1998) and Dupac et al. (2001)), Messier 17 (Dupac et al. (2002)), Cygnus B, and the dusty envelope surrounding the young massive star GH2O 092.67+03.07 in NCS (Bernard et al. (2003)). The $`\rho `$ Ophiuchi low-mass star-forming region has also been observed, as well as the edge-on spiral galaxy NGC 891 (Dupac et al. (2003)). In this letter, we analyze constraints that PRONAOS data establish on the dust spectral index and its relation to the temperature. Ristorcelli et al. (1998) and Bernard et al. (1999) gave tentative evidences of a spectral index larger than 2, respectively in Orion and the Polaris flare. Dupac et al. (2001) and Dupac et al. (2002) showed large variations of the spectral index, including values above 2, and the presence of an anticorrelation between the temperature and the spectral index in massive star-forming regions. In this letter, we extend the analysis to other regions observed with PRONAOS.

## 3 Analysis

### 3.1 Bases of the analysis

For Orion (Dupac et al. (2001)), M17 (Dupac et al. (2002)) and Cygnus (Boudet et al. (2003)), we analyze the intensity maps using the following procedure. With the PRONAOS and IRAS intensity maps, we produce temperature ($`T`$) and spectral index ($`\beta `$) maps following the processing described in Dupac et al. 2002 (see also Dupac et al. (2001)). We fit a modified black body law to the spectra: $`I_\nu =ϵ_0B_\nu (\lambda ,T)(\lambda /\lambda _0)^\beta `$, where $`I_\nu `$ is the spectral intensity (MJy/sr), $`ϵ_0`$ is the emissivity at $`\lambda _0`$ of the observed dust column density, $`B_\nu `$ is the Planck function, $`T`$ is the temperature and $`\beta `$ is the spectral index. In most of the areas, we use either only PRONAOS data or PRONAOS + IRAS 100 $`\mu `$m data. The reasons why we do not use IRAS data everywhere are that first, some very intense regions are saturated in the 100 $`\mu `$m IRAS band, second, many faint regions which are clearly detected with PRONAOS are very noisy with IRAS, or a large-scale warm component dominates in the IRAS bands. In the non-saturated warmest areas ($`>`$ 70 K), we use IRAS 60 $`\mu `$m data too (see Dupac et al. (2002) for more details).

In this work, we restrain the analysis to all fully independent (3.5 side) pixels for which both relative errors on the temperature and the spectral index are less than 20%. This procedure is motivated notably by the fact that there exists a degeneracy between the temperature and the spectral index. This effect is mainly due to the relatively large error bars on the spectral index in cold regions and on the temperature in warm regions (Rayleigh-Jeans limit). Dupac et al. (2001) and Dupac et al. (2002) have shown by fitting simulated data that this artificial anticorrelation effect was small compared to the effect observed in the data. For this work, we have also done other simulations, in a slightly different way from the previous ones: we consider temperature values in the same temperature range as the observations, but with a constant spectral index (1.5), we produce simulated data with a realistic noise level, and we fit them with the usual procedure. With PRONAOS simulated data alone, we find that the relative residual error on the spectral index is 22% rms. With PRONAOS + IRAS 100 $`\mu `$m data, the same error is only 7% rms on the spectral index and 9% on the temperature. The induced inverse correlation effect is very limited. Thus, this noise effect is clearly not capable to explain the $`T`$-$`\beta `$ distribution found in the data. Here we still reduce it by restraining the analysis to pixels for which errors on both the temperature and the spectral index are lower than 20%. Also, we disregard pixels for which at least one of the bands has a spectral intensity below 10 MJy/sr. Moreover, all spectra were visually checked in order to exclude pixels which could be contaminated by noise features. This provides us with a set of 31 pixels for M17, 19 for Orion and 60 for Cygnus. For the cirrus clouds in Polaris and Taurus, we have, respectively, two and one data points available (Bernard et al. (1999) and Stepnik et al. (2003)). For the NCS data point, we fit the PRONAOS fluxes as given in Bernard et al. (2003). We also present five data points from $`\rho `$ Ophiuchi (Ristorcelli et al. (2003)) and three from the NGC 891 edge-on spiral galaxy (Dupac et al. (2003)). We present typical examples of the obtained spectra in Fig. 1.

### 3.2 Dependence on the temperature of the spectral index

We present in Fig. 2 and 3 temperature versus spectral index plots of the above-defined pixel set. For clarity, we present all points in Fig. 3 without the error bars, and a plot with the error bars in Fig. 2 but without the numerous points from Cygnus and M17. These data points have similar error bars to the Orion data points. The temperature in this data set ranges from 11 to 80 K, and the spectral index also exhibits large variations from 0.8 to 2.4. One can observe an anticorrelation on these plots between the temperature and the spectral index, in the sense that the cold regions have high spectral indices around 2, and warmer regions have spectral indices below 1.5. In particular, no data points with $`T>`$ 35 K and $`\beta >`$ 1.6 can be found, nor points with $`T<`$ 20 K and $`\beta <`$ 1.5. All cold regions below 18 K have spectral indices above 1.8, except one point with a large error bar, which is well compatible with a value above 1.8.

This anticorrelation effect is present for all objects in which we observe a large range of temperatures, namely Orion, M17, Cygnus and $`\rho `$ Ophiuchi. It is also remarkable that the few points from other regions are well compatible with this general anticorrelation trend. Following the type of the observed object, the temperature distribution of the points is different, but the points seem to follow this same general trend of an inverse dependence of the spectral index on the temperature. Hereafter, we therefore consider that this effect can be simply investigated with the different objects altogether.

Fitting the data under the assumption of a constant spectral index leads to $`\chi ^2`$ / degrees of freedom (d.o.f. hence on) values of 1835/121 for $`\beta `$=2 and 377/121 for $`\beta `$=1.5. These high values of the $`\chi ^2`$/d.o.f. exclude the possibility of matching the data with an invariant spectral index. The linear regression fit of $`\beta `$($`T`$), considering the temperature as an independent variable, gives a slope of $``$0.014 K-1, and $`\beta `$ = 2.01 at $`T`$ = 0. Then the linear correlation coefficient is $``$0.79, and the $`\chi ^2`$/d.o.f. is 137/120 (=1.14).

Since the points seem quite well distributed along a hyperbolic-shape line, we perform a hyperbolic fit of the data using: $`\beta =\frac{1}{\delta +\omega T}`$, where $`\delta `$ and $`\omega `$ are free parameters. Considering that the temperature is the independent variable, the best fit gives $`\delta `$=0.40 $`\pm `$ 0.02 and $`\omega `$=0.0079 $`\pm `$ 0.0005 K-1, with $`\chi ^2`$/d.o.f. = 120/120. The temperature dependence of the emissivity spectral index is thus very well fitted with an hyperbolic approximating function. Note that since the number of free parameters is the same as in the linear regression, the lower $`\chi ^2`$ is significant, so that we can consider the hyperbolic fit to be better.

If now we consider that the spectral index is the independent variable, fitting the $`T`$($`\beta `$) function by the linear regression gives a large $`\chi ^2`$/d.o.f. value of 1099/120. In this case, the hyperbolic fit gives $`\chi ^2`$/d.o.f. = 1123/120. These bad fits are mainly due to the large dispersion of the spectral indices that we observe at low temperatures.

## 4 Discussion

In order to interpret this inverse temperature dependence of the spectral index, several hypotheses can be investigated. Mixtures of dust components with different temperatures, but having the same spectral index, can be put forward to explain apparent low indices. However, this effect was shown to be small compared to the amplitude of the observed spectral index variations (Dupac et al. (2002)).

As stated in the introduction, low indices have been observed in active environments such as circumstellar disks and warm molecular cloud cores. PRONAOS observations of Orion, M17, Cygnus B and $`\rho `$ Ophiuchi show such low spectral indices around 1. This could be due to the increase of the grain size in dense environments, as predicted by Miyake & Nakagawa (1993), whose models with different grain size distributions show that the maximum grain radius adopted is crucial for the long-wavelength emissivity slope. This is also the reason invoked by Goldsmith et al. (1997) to explain their results, which showed an inverse dependence of the spectral index on the optical depth. However, we do not find such correlation in our data. Moreover, the range of variations that they find for the spectral index (from 1.6 to 3.2) is significantly different from ours (from 1 to 2.4), though with a similar variation amplitude. Also, these grain growth effects are invoked for very peculiar dense environments around stars rather than for large scale observations such as the PRONAOS ones. Therefore, it is difficult to consider that size distribution effects are dominant in the $`T`$-$`\beta `$ anticorrelation that we observe, though they may play a role.

Another slightly different explanation could be that one finds different chemical composition or physical state of the grains following the physical conditions of the medium, in particular following the temperature. Cold environments could harbor fluffy silicate grains including ice compounds, having spectral indices of 2 (e.g. Krügel & Siebenmorgen (1994)), or simply silicate or graphitic particles (e.g. Draine & Lee (1984)), or even olivine core or fused quartz grains with ice mantles (Aannestad (1975), $`\beta `$ 3) to explain high indices above 2. Warm regions could harbor aggregates of silicates, porous graphite or amorphous carbon and therefore have a spectral index around 1 (Mathis & Whiffen (1989)). Nevertheless, more observations are needed to constrain the chemical structure of the large grains in various media and confirm this quite speculative hypothesis.

Finally, an explanation may be searched for in a temperature dependence of the intrinsic optical properties of the grains (see Henning & Mutschke (1997) for a review of the different physical processes which could be involved). This can be linked to laboratory results (Agladze et al. (1996) and Mennella et al. (1998)) which showed a T-$`\beta `$ anticorrelation on interstellar grain analogs. Mennella et al. (1998) measured the absorption coefficient of crystalline and amorphous grain analogs between 20 $`\mu `$m and 2 mm wavelength, in the temperature range 24-295 K. They derived an anticorrelation between the temperature and the spectral index, and attributed it to two-phonon difference processes (Sparks et al. (1982)). Agladze et al. (1996) measured absorption spectra of silicates between 0.7 and 2.9 mm wavelength and deduced an anticorrelation between the spectral index and the temperature in the temperature range 10-25 K, for two different silicate precursors. They attributed the underlying absorption process to a resonant tunneling effect between ground states of two-level systems (Phillips (1972)). The temperature dependence of the spectral index which they show is in a range of temperatures and spectral indices close to our measurements: they measure spectral indices between 1.8 and 2.6 for temperatures ranging from 10 to 25 K, while we measure spectral indices roughly between 1.4 and 2.4 in the same temperature domain. However, the wavelength range is not the same, and there is no evidence that this effect is simply transposable from the millimeter range to the PRONAOS submillimeter range, since experimental results suggest a cut-off around 500 $`\mu `$m for this phenomenon (Fitzgerald et al. (2000)). Furthermore, at submillimeter wavelengths and higher temperatures, thermally-activated relaxation processes and temperature-dependent absorption associated to transitions to excited levels of two-level systems (Fitzgerald et al. (2000), Hutt et al. (1989)) might dominate. These issues will be treated in a forthcoming paper (Mény et al., in prep.) This possible discrepancy between the submillimeter and the millimeter spectral indices could explain the small difference observed between the WMAP (millimeter) spectral index of 2.2 and the PRONAOS one, which is around 2 at DIRBE (rather low) large-scale temperatures.

## 5 Conclusion

We have presented a compilation of PRONAOS-based results concerning the temperatures and spectral indices measured in a variety of regions of the interstellar medium. The temperatures observed range from 11 to 80 K, which provides a large field of investigation for measuring variations of the dust emissivity. We give strong evidence of large variations of the submillimeter spectral index between 0.8 and 2.4, as well as of the existence of an inverse dependence of the emissivity spectral index on the temperature. This dependence is well fitted by the function $`\beta =\frac{1}{0.4+0.008T}`$. Several interpretations are possible for this effect: one is that the grain sizes change in dense environments, another is that the chemical composition of the grains is not the same in different environments and that this correlates to the temperature, a third one is that there is an intrinsic dependence of the spectral index on the temperature, due to quantum processes such as two-level tunneling effects. Additional modeling, as well as additional laboratory measurements and astrophysical observations, are required in order to discriminate between these different interpretations. However, the most likely explanation could be that there is an intrinsic dependence of the emissivity on the temperature. A new dust semi-empirical model is needed to take into account these variations of the spectral index.

## 6 Acknowledgements

We thank very much H. A. Sabat for fruitful discussions on statistics and B. Parise for her careful reading of the manuscript. We are indebted to the French space agency Centre National d’Études Spatiales (CNES), which supported the PRONAOS project. We are very grateful to the PRONAOS technical teams at CNRS and CNES, and to the NASA-NSBF balloon-launching facilities group of Fort Sumner (New Mexico).


text:
# Dynamical ionization ignition of clusters in intense and short laser pulses

## I Introduction

Clusters bridge the gap between bulk material and single atoms. When placed into a laser field all atoms inside (not too big) clusters experience the same laser field, contrary to the bulk where a rapidly created surface plasma and the skin effect prevents the laser from penetrating deeper into the target. In rarefied gases, on the other hand, the laser can propagate but the density is too low to yield high absorption of laser energy and high abundances of fast particles. Hence, the absorption of laser energy is expected to be optimal in clusters. In fact, highly energetic electrons shao , ions ditmirenature ; ditmirePRA ; springate , photons mcpherson ; ditmireJPB ; teravet , and neutrons originating from nuclear fusion zweibackDD  were observed in laser cluster interaction experiments. A prerequisite for clusters as an efficient source of energetic particles and photons is the generation of high charge states inside the cluster. Several mechanisms for the increased ionization (as compared to single atoms in the same laser field) have been proposed. Boyer et al. boyer  suggested that collective electron motion (CEM) inside the cluster is responsible for inner shell vacancies the corresponding radiation of which was observed in experiments schroeder . Rose-Petruck et al. rosepetru  introduced the so-called “ionization ignition” (IONIG) where the combined field of the laser and the ions inside the cluster leads to increased inner ionization, i.e., electrons are more easily removed from their parent ion as if there was the laser field alone. The removal of electrons from the cluster as a whole is called outer ionization. An “outer-ionized” cluster will expand because of the Coulomb repulsion of the ions while a heated, quasi-neutral cluster expands owing to the electron pressure. According to experiments the latter viewpoint of a nanoplasma (NP) ditmirePRA ; milchberg  seems to be appropriate for big clusters of about $`10^4`$ atoms or more kim  while numerical simulations indicate that the Coulomb explosion dynamics prevail for smaller clusters ishika ; rusek . For recent reviews on the interaction of strong laser light with rare gas clusters see posthumus ; krainsmirn .

In this paper we are aiming at clarifying the early ionization dynamics in laser cluster interaction. Because the plasma has to be created first the NP model is not applicable. Although we shall find CEM to be indeed an important ingredient for the early ionization dynamics of clusters, the origin of the CEM should be explained rather than assumed beforehand. The IONIG model predicts that “as the ionization proceeds beyond the first ionization stage, strong electric fields build up within the cluster that further enhance ionization” rosepetru . Here we are interested in how these strong inner fields are generated. Once the cluster is charged after the first electrons are removed, the remaining bound electrons experience the attraction of neighboring ions. This attraction will be largest for bound electrons of ions sitting near the cluster boundary. The other ions will pull these electrons into the cluster interior. This force might be supported by the electric field of the laser, thus leading to further ionization. The IONIG model was successfully called on for interpreting the experimental results in snyder .

Although this IONIG mechanism is appealing, the details of the ionization dynamics remain unclear. If the force exerted by the other ions is strong enough to ionize further it should be also strong enough to keep the already freed electrons inside the cluster. Hence, it remains to be explained how the electrons are removed from the cluster as a whole (outer ionization). The electrons, after inner ionization, may as well shield the space charge of the ions so that the IONIG mechanism would come to an end, and the NP model would take over even in small clusters. So, why does a strong electric field build up inside the cluster whose interplay with the electric field of the laser is constructive for inner as well as outer ionization?

Classical particle methods are frequently applied to investigate the electron and ion dynamics of clusters in laser fields rosepetru ; lastI ; lastII ; ishika ; saalmann ; siedschlagI ; siedschlagII  since a full three-dimensional quantum treatment is out of reach with nowadays computers. In these classical accounts inner ionization was either accounted for by sampling the quantum mechanical electron distribution of the bound state by a microcanonical, classical ensemble of electrons, or by using ionization rates so that the electron dynamics was simulated only after inner ionization. Semi-classical Thomas-Fermi theory was employed in rusek  to study the explosion dynamics of clusters consisting of up to 55 atoms. In our work we use time-dependent density functional theory rungegross ; burkegross ; grossdobson  since we are mainly interested in the early ionization dynamics of clusters were quantum mechanics is important.

The article is organized as follows. In Section II the numerical model is introduced. In Section III results are presented concerning the groundstate properties of the model cluster (III.1), the electron dynamics (III.2), the formation of collective electron motion inside the cluster and outer ionization (III.3), the effect of the collective electron motion on inner ionization (III.4), the frequency dependence of outer ionization (III.5), and the influence of the ionic motion (III.6). Finally, we conclude in Section IV.

Atomic units (a.u.) are used throughout the paper unless noted otherwise.

## II Numerical model

Time-dependent density functional theory (TDDFT) is employed to study the ionization dynamics of small and medium size rare gas clusters in intense and short laser pulses. To that end the spin degenerate time-dependent Kohn-Sham equation (TDKSE) is solved. However, solving the TDKSE in three spatial dimensions (3D) for rare gas clusters in laser fields is too demanding for todays computers. One reason for this is that, contrary to metal clusters calvay  or fullerenes bauerc60 , the electrons in the ground state of rare gas clusters are not delocalized so that jellium models where the ionic background is smeared out and assumed spherical are not applicable. In order to make a numerical TDKS treatment feasible two simplifications were made. First, as in previous studies of clusters in laser fields brewI ; brewII ; veniard ; grigor , the dimensionality of the problem was restricted to 1D, namely to the direction of the linearly polarized laser field described by the vector potential $`𝑨(t)=A(t)𝒆_x`$ in dipole approximation dipoleremark . To that end the Coulomb interactions were replaced by soft-core Coulomb interactions, i.e., $`|𝒓𝒓^{}|^1[(xx^{})^2+a_{ee}^2]^{1/2}`$ and $`|𝒓𝑹_k|^1[(xX_k)^2+a_{ei}^2]^{1/2}`$ for the electron-electron interaction and the electron-ion interaction, respectively. The smoothing parameters $`a_{ee}`$ and $`a_{ei}`$ may be chosen to yield ionization energies similar to real 3D systems. The second simplification was the use of the exchange-only local density approximation (XLDA) so that $`V_{\text{xc}}[n(x,t)]=V_{\text{XLDA}}[n(x,t)]=\alpha (3n(x,t)/\pi )^{1/3}`$ where $`n(x,t)`$ is the electron density. The pre-factor $`\alpha `$ would be unity in full 3D XLDA calculations. In the 1D model it may be chosen to yield satisfactory ground state properties (we used $`\alpha =3/4`$). The ionic potential was $`V_{\text{ion}}(x)=_{k=1}^{N_{\text{ion}}}Z[(xX_k)^2+a_{ei}^2]^{1/2}`$ with constant nearest-neighbor distances $`X_{k+1}X_k=d2r_{\text{WS}}`$ where $`r_{\text{WS}}`$ is the Wigner-Seitz radius. One may look at the 1D ion chain as representing those $`N_{\text{ion}}`$ ions of a 3D spherical cluster which are situated along the diameter parallel to the linearly polarized laser field. The cluster radius then is $`R(N_{\text{ion}}1)d/2`$, and the number of ions in the real 3D cluster $`N_{\text{3D}}=R^3/r_{\text{WS}}^3`$ would be $`(N_{\text{ion}}1)^3`$.

The TDKSE for the $`N_{\text{KS}}`$ Kohn-Sham (KS) orbitals $`\mathrm{\Psi }_j(x,t)`$, $`j=1,\mathrm{},N_{\text{KS}}`$ reads

$$\text{i}\frac{}{t}\mathrm{\Psi }_j(x,t)=H_{\text{KS}}\mathrm{\Psi }_j(x,t)$$
(1)
where

$$H_{\text{KS}}=\frac{1}{2}\mathbf{}^2+U+V_{\text{XLDA}}+V_{\text{ion}}+V_{\text{laser}}.$$
(2)
From the doubly spin-degenerate KS orbitals the total electron density $`n(x,t)=2_{j=1}^{N_{\text{KS}}}|\mathrm{\Psi }_j(x,t)|^2`$ is calculated. $`U(x,t)=\text{d}x^{}n(x^{},t)[(xx^{})^2+a_{ee}^2]^{1/2}`$ is the 1D Hartree potential (accounting for the electron-electron repulsion), and $`V_{\text{laser}}=\text{i}A(t)_x`$ governs the interaction with the laser (taken in velocity gauge with the purely time-dependent term $`A^2`$ transformed away). Eq. (1) was solved using the Crank-Nicolson method with the accuracy of the spatial derivatives boosted to fourth order. Since the Hamiltonian (2) itself depends (through the density) on the KS orbitals a predictor-corrector method should be used in combination with the usual Crank-Nicolson time-propagation. In practice, however, the evaluation of the Hamiltonian (2) using the density from the previous time step is usually accurate enough.

For the ion-ion interaction a soft-core potential $`Z^2[(X_jX_k)^2+a_{ii}^2]^{1/2}`$ was assumed as well. The ions $`j=1,\mathrm{},N_{\text{ion}}`$ were moved according to their classical, non-relativistic equations of motion,

$`M\ddot{X_j}`$ $`=`$ $`Z^2{\displaystyle \underset{k=1}{\overset{N_{\text{ion}}}{}}}{\displaystyle \frac{X_jX_k}{[(X_jX_k)^2+a_{ii}^2]^{3/2}}}`$ (4)
$`Z{\displaystyle \text{d}x\frac{n(x,t)(X_jx)}{[(X_jx)^2+a_{ei}^2]^{3/2}}}`$
$`Z{\displaystyle \frac{}{t}}A(t),`$
using the Verlet algorithm.

The results which will be presented in the Sections III.1III.5 were obtained for an ion mass $`M=1311836`$ (Xe atom) and $`a_{ii}=1`$. Due to the short laser pulse durations ($`<40`$ fs) and the modest charge states $`<4`$ created, the ionic motion did not appreciably affect the ionization dynamics of our 1D model. However, since in a real 3D cluster consisting of $`(N_{\text{ion}}1)^3`$ ions the Coulomb explosion would be more violent than in our 1D model with $`N_{\text{ion}}`$ ions (of the same ion mass and charge state) the issue of ionic motion will be discussed separately in Section III.6.

Finally, before presenting our results, we briefly relate our approach to the work in Refs. brewI ; veniard ; grigor . Instead of the simple XLDA, Véniard et al. veniard  used in their 1D TDDFT approach the more advanced exchange potential proposed by Krieger, Li, and Iafrate (KLI) kli  but restricted their studies to ion chains up to 7 fixed atoms with one active shell only. Besides the extensive study of the harmonic radiation emitted by the 1D cluster, Véniard et al. focussed on the dependence of ionization on the internuclear distance and on the number of ions in the chain. This is different from our goal to highlight the dynamics of IONIG. The methodological approach instead is very similar. The only differences are the mobile ions and the choice of simple XLDA instead of the KLI exchange potential. The latter simplification gave us the freedom to simulate bigger clusters and more active shells per atom.

The 1D models in grigor  and brewI  were called on to study the explosion dynamics of Xe clusters. Since for this purpose a full TDDFT treatment with many KS orbitals is too demanding even in 1D, further simplifications had to be adopted. In grigor  a single orbital, representing all electrons, was introduced while in brewI  time-dependent Thomas-Fermi theory was used. Our main emphasis in the present paper is complementary to this work since we are not dealing so much with the Coulomb explosion of the cluster but are predominantly interested in the early electron dynamics.

## III Results and discussion

### III.1 Groundstate properties

Let us consider a chain of $`N_{\text{ion}}=9`$ ions with nearest-neighbor distance $`d=8`$ and nuclear charge per ion $`Z=4`$. Hence, the ground state of the neutral cluster consists of $`N=36`$ electrons. The smoothing parameters for the soft-core Coulomb interaction were simply chosen to be $`a_{ee}=a_{ei}=a_{ii}=1`$. In Fig. 1 the ground state density and the various contributions to the total potential are plotted: the ion potential $`V_{\text{ion}}`$, the Hartree potential $`U`$, and the exchange potential $`V_{\text{XLDA}}`$. Although the ionic potential alone has its absolute minimum at the central ion the total potential, including the classical Hartree-repulsion and XLDA, consists of nine almost identical and equidistantly separated potential wells, each locally similar to that of the corresponding individual atom. Consequently, the ground state density displays nine almost equal, well localized peaks, and the energy levels of the cluster are approximately at the same positions as those of the atom, namely around $`1.17`$ for the $`2N_{\text{ion}}=18`$ inner electrons, and around $`0.23`$ for the $`18`$ outer electrons (note, that in the 1D model there are only 2 electrons per shell). Since in XLDA Koopman’s theorem is usually not well fulfilled the energy of the highest occupied orbital does not equal the ionization energy for removing the outermost electron. Calculating the ionization energy from the difference of the total energy of the neutral and the singly ionized cluster $`\mathrm{\Delta }_{\text{cluster}}=0.26`$ was obtained. Doing the same for the single atom led to $`\mathrm{\Delta }_{\text{atom}}=0.49`$ which is a reasonable value for rare gas atoms.

### III.2 Electron dynamics

Let us start by comparing the electron motion in the cluster with that of a single atom. In Fig. 2 contour plots of the logarithmic density are shown vs. space and time. The results for the full cluster (a), the single atom (b), and an artificial cluster made of noninteracting atoms (c) in a laser pulse of frequency $`\omega =0.057`$ ($`\lambda =800`$ nm) and rather modest field amplitude $`\widehat{E}=0.033`$, corresponding to $`3.9\times 10^{13}`$Wcm<sup>-2</sup>, are shown. The laser field was ramped up linearly over $`3`$ cycles, held $`8`$ cycles constant, and ramped down again over $`3`$ cycles (hereafter called a $`(3,8,3)`$-pulse) so that the pulse duration was $`37`$fs. The density in the contour plot (c) was calculated by assuming that all of the $`N_{\text{ion}}=9`$ atoms behave as the single atom in plot (b).

It is seen that the electron dynamics of the cluster (a) and the noninteracting atoms (c) differ significantly from each other already at $`t150`$ because more electron density leaves to the right at that time instant in (a) than it does in (c). This behavior of stronger ionization of the cluster than of the individual atoms continues during the subsequent half laser cycles.

A qualitatively different electron dynamics inside the cluster emerges for $`t400`$. While in the cluster (a) an accumulation of electron density bounces from one boundary of the cluster to the other such an electron dynamics, of course, cannot build up in the ensemble of independent atoms (c). The formation of an electron wave packet which travels through the entire cluster in step with the laser field is remarkable in view of the fact that the excursion of a free electron in the laser field $`\widehat{E}=0.033`$ amounts to $`\widehat{x}=10.2`$ only while the diameter of the cluster is $`2R64`$. By varying the cluster size (results for $`N_{\text{ion}}=17`$ will be presented in the following subsection) it was found that the formation of the bouncing wave packet at laser intensities where $`\widehat{E}/\omega ^2<R/2`$ is a robust phenomenon, not sensitive to the cluster parameters. Test runs with more advanced exchange potentials such as the Slater potential with self-interaction corrected XLDA and the KLI potential kli  were performed to ensure that the coherent electron motion is not an artifact of plain XLDA. However, the wave packet formation is sensitive to the laser frequency, as will be shown in subsection III.5.

### III.3 Formation of collective electron motion inside the cluster and outer ionization

It is useful to study the phase relation between the oscillating electron density inside the cluster and the laser field in order to understand the formation of the electron wave packet with unexpected large excursion amplitude. Given a laser field $`𝑬(t)`$ of optical (or lower) frequency, the polarization of an atom is $`𝑬(t)`$ because the bound electrons are able to follow adiabatically the force exerted by the field. Hence, the phase lag of the polarization with respect to the laser field is $`\mathrm{\Delta }\varphi =\pi `$. Free electrons, on the other hand, oscillate perfectly in phase $`𝑬(t)`$ so that $`\mathrm{\Delta }\varphi =0`$. Energy absorption from the laser field is low in both cases because $`\text{d}t\dot{𝒙}(t)𝑬(t)0`$ (with $`𝒙(t)`$ the expectation value for the position of an electron). During ionization there is necessarily a transition from $`\mathrm{\Delta }\varphi =\pi `$ to $`0`$ where energy absorption can take place. During the ionization of atoms this transition occurs rapidly while in clusters, after inner ionization, the electrons may be still bound with respect to the cluster as a whole. Hence, the free motion of electrons $`𝑬(t)`$ comes to an end at latest when they arrive at the cluster boundary. There, they either escape from the cluster, contributing to outer ionization, or they are reflected so that their phase relation with the driving laser is affected, leading on average to an enhanced absorption of laser energy. Although collisions of the electrons with ions are included in our TDDFT treatment the effect of the boundary on the electron dynamics clearly dominates. The unimportance of electron ion collisions in medium size and small clusters was pointed out in ishika  while the relevance of boundary effects was recently affirmed in megi  within the NP model remarkII .

In Fig. 2a it is seen that some electrons enter about ten atomic units into the vacuum before they are pulled back by the cluster charge. This is reminiscent of what in laser plasma physics is called Brunel effect brunel , “vacuum heating” vacheat , or, more expressively, “interface phase mixing” mulser . Thanks to the fact that the fast electrons leave the cluster (and slow electrons are accelerated) a filter effect comes into play so that a wave packet can form that oscillates with the laser frequency and an excursion amplitude of about the cluster radius $`R`$.

In order to underpin this scenario the dipole of the electron density inside the cluster $`x_{\text{inner}}(t)=_{(R+d/2)}^{R+d/2}\text{d}xxn(x,t)`$ was calculated for several laser and cluster parameters. The results are shown in Fig. 3. In panel (a) the laser and cluster parameters were the same as in Fig. 2. One sees that during the first few laser periods the electrons indeed move $`E(t)`$ (green, dotted curve) as indicated by the first, green bar at the top of the $`x_{\text{inner}}`$-plot. Then, the electrons inside the cluster get out of phase with the laser for about nine cycles (red bar) so that $`\mathrm{\Delta }\varphi \pi /2`$. During this period the dipole amplitude $`\widehat{x}_{\text{inner}}`$ is particularly high remark , and the wave packet bouncing inside the cluster is clearly visible in Fig. 2. Finally, towards the end of the laser pulse the phase relation of the few electrons which were removed from their parent ions but did not make it to leave the cluster becomes that of free electrons (blue bar), i.e., $`x_{\text{inner}}(t)E(t)`$. In the lower plot the number of electrons inside the cluster $`N_{\text{inner}}`$ is plotted and compared with $`N_{\text{ion}}`$ times the result for the single atom. It is seen that during the first phase (green bar) ionization of the cluster proceeds similar to the single atom case. However, when the phase lag is shifted to $`\mathrm{\Delta }\varphi =\pi /2`$ the cluster continues to ionize while the single atom ionization comes to an end.

In Fig. 3b the same is shown for a higher laser intensity. Essentially, the phase $`\mathrm{\Delta }\varphi `$ behaves in the same way but this time, due to the stronger laser field, ionization of the outer shell is almost completed during the first few laser cycles. Hence, the final average charge state is almost the same for the cluster and the single atom. The period where the phase lag is about $`\mathrm{\Delta }\varphi \pi /2`$ lasts only a few laser cycles (red bar) and so does the wave packet motion inside the cluster. In Fig. 3c the result for a bigger cluster ($`N_{\text{ion}}=17`$, $`Z=2`$) is presented, revealing a qualitatively similar scenario as in (a) with the bouncing wave packet surviving for about 9 cycles. Finally, in panel (d) a higher laser frequency was used ($`\omega =0.18`$) while keeping the cluster parameters as in (a) and (b). The $`N_{\text{inner}}`$-plot reveals that the single atom ionizes more efficiently than the cluster for these laser parameters. We will come back to the frequency dependence of outer ionization in subsection III.5.

### III.4 Effect of the collective electron motion on inner ionization: dynamical ionization ignition

The formation of collective electron dynamics as exposed in the previous subsection explains how the absorption of laser energy is increased due to a phase shift into resonance with the driving field, and how the electrons, after inner ionization, are efficiently transported out of the cluster (outer ionization). The increased inner ionization still remains to be analyzed. IONIG states that the presence of the other ions is responsible for the more efficient removal of bound electrons. This is because two neighboring ions form a potential barrier (cf. the $`V_{\text{ion}}`$ curve in Fig. 1) through which an electron may tunnel when the whole cluster is submitted to an electric field so that the entire potential is tilted. However, we found from our numerical studies that this energetic advantage of a bound electron inside the cluster as compared with an electron in the corresponding single ion is not very pronounced for medium size and small clusters at moderate charge states. Instead, we propose a dynamical version of IONIG where the previously introduced collective electron motion plays an important role. Coherent electron motion was suggested to be responsible for inner shell vacancies in Xe clusters, leading to x-ray emission mcpherson ; schroeder . In our numerical model there are only two shells and we did not find particularly high line emission from the cluster as compared to the single atom. However, it is possible that, owing to the lack of dynamical correlation in XLDA, the interaction of the electrons in the wave packet with the still bound electrons is underestimated in our model. Thus, dynamical IONIG might be the mechanism behind the experimental results reported in mcpherson .

Despite the fact that in our mean field approach there are no “hard” collisions of the electrons in the wave packet with the still bound electrons, the electric field associated with the oscillating electron packet already enhances inner ionization. In Fig. 4 two snapshots of the total effective potential and the electron density are presented for the same laser parameters as in Fig. 2. In panel (a), electron density and total potential are shown for a time where the electron wave packet is close to the left cluster boundary (the red bar at the bottom of the density plot indicates $`x_{\text{inner}}(t)`$). The contour plot (b) shows for each ion at position $`X_i`$ the “difference density” $`_{X_id/2}^{X_i+d/2}\text{d}xn(x,t)N_{\text{inner}}/N_{\text{ion}}`$ which indicates whether there is a lack of electron density at that position inside the cluster (black and dark colors) or whether there is excess density (yellow and light colors) compared to the average density $`N_{\text{inner}}/N_{\text{ion}}`$.

In panel (a), at time $`t=662.5`$ the wave packet is close to the left boundary of the cluster while there is a lack of electron density near the right boundary. This charge distribution leads to a force $`F_{\text{int}}`$ on the other electrons (pointing to the right) and therefore increases the ionization probability. The electric field of the laser instead is close to zero so that $`F_{\text{laser}}`$ is small. The situation approximately a quarter of a laser cycle later is shown in panel (c). The wave packet is at the center of the cluster, moving with maximum velocity to the right and repelling bound or slow electrons in front of it. The force $`F_{\text{laser}}`$ is close to its maximum value at that time, pointing into the direction in which the wave packet moves. The ionization probability is, again, greater than with the laser field alone. Thus, during the course of a laser cycle the total force $`F_{\text{laser}}+F_{\text{int}}`$ clearly leads to higher ionization as if there was the laser field only.

The fact that the electron wave packet dynamics indeed increases inner ionization is underpinned by Fig. 5 where the number of electrons in the region $`\pm d/2`$ around the central ion in the cluster was subtracted from the initial value $`Z=4`$ and is compared with the corresponding single atom result. The laser parameters were the same as in Figs. 2 and 4. While for the single atom the ionization is completed for $`t>700`$ the average electron density around the central ion in the cluster is still decreasing. When the electron wave packet sweeps over the central ion the density is temporarily increased, leading to local minima in the curve of Fig. 5. When the wave packet is closest to one of the two cluster boundaries the lack of electrons around the central ion is maximal. The absolute increase of this maximum each half cycle means ongoing inner ionization. In contrast, the single atom, where only the laser field is present but no wave packet can form, does not ionize any further.

### III.5 Dependence of outer ionization on the laser frequency

The interaction of the model cluster with $`N_{\text{ion}}=9`$ was investigated for the two different laser frequencies $`\omega _l=0.057`$ and $`\omega _h=0.18`$ (corresponding to $`800`$ and $`254`$ nm, respectively) and laser intensities between $`4\times 10^{12}`$ and $`10^{16}`$Wcm<sup>-2</sup>. The pulses were of $`(3,8,3)`$-shape for both frequencies, that is, the pulse durations were $`T_l37`$ fs and $`T_h=12`$ fs, respectively. After the laser pulse, the average charge state in the cluster $`Z_{\text{av}}=ZN_{\text{inner}}/N_{\text{ion}}`$ was calculated. Note that $`Z_{\text{av}}`$ only yields information about outer ionization for it does not distinguish between electrons that are still bound to their parent ions and those which move inside the cluster.

In Fig. 6 the average charge state is plotted vs. the laser intensity for the two frequencies $`\omega _h`$ and $`\omega _l`$. The results for the single atom are also shown. As discussed in the previous subsections, it is seen that in the low frequency case the atoms in the cluster are stronger ionized than an individual atom in the same laser field. Both charge states come close only for $`Z_{\text{av}}=2`$, that is when the two electrons of the first shell are removed but the two electrons of the next shell are still strongly bound. The stepwise increase of the charge state due to the electronic shell structure is very pronounced in the low-frequency cluster case as well as for the single atoms at both, low and high frequency. The cluster in the high-frequency field instead shows a very different behavior: between $`10^{14}\text{Wcm}\text{-2}`$ and the threshold to the inner shell at $`5\times 10^{15}\text{Wcm}\text{-2}`$, the charge state of the single atom is higher than the average charge state in the cluster.

In Fig. 7 the dynamics of the cluster electrons is shown for a $`(3,8,3)`$-pulse of frequency $`\omega _h=0.18`$ and peak field amplitude $`\widehat{E}=0.099`$, corresponding to an intensity of $`3.44\times 10^{14}`$Wcm<sup>-2</sup>. It is seen from Fig. 6 that for this intensity the single atom ionizes more efficiently than the cluster as a whole, contrary to what happens at the lower frequency $`\omega _l=0.057`$. Fig. 7 reveals that the electrons, although removed from their parent ions, mostly remain inside the cluster. From plot (c) one infers that if the atoms inside the cluster were independent there would be a strong electron emission for $`100<t<350`$. The emitted electrons have sufficient high kinetic energy to escape from their parent ion (and, thus, from the “independent atom”-cluster). Contour plot (a), instead, shows that in the real cluster a significant fraction of the electrons near the cluster boundaries return due to the space charge created by all the ions. A wave packet dynamics that could enhance outer ionization, as in the low frequency result of Fig. 2, does not form. We attribute this to the fact that the initial inner ionization occurs less adiabatic (multiphoton instead of tunneling ionization). Moreover, the excursion $`\widehat{x}=3.06<d/2`$ is too small to trigger any collective motion. The dynamics inside the cluster is rather “splash-like,” as can be inferred from the strongly fluctuating electron density between the ions in Fig. 7a. Hence, although at high laser frequencies inner ionization is high, outer ionization remains low since there is not the coherent electron dynamics supporting outer ionization. Consequently, for creating quasi-neutral nanoplasmas and suppressing Coulomb explosion the use of high frequency lasers is favorable. Reduced ionization of clusters in laser fields of, however, many times higher frequency (to be generated by x-ray free electron lasers in the near future) was also found in the numerical simulations of saalmann . On the other hand, in the soft-x-ray FEL experiment performed by Wabnitz et al. wabnitz  increased ionization of Xe clusters (as compared to single atoms) was observed at $`98`$ nm wavelength, $`100`$ fs pulse duration, and intensities up to $`7\times 10^{13}`$Wcm<sup>-2</sup>. The reason for this unexpected behavior at short wavelengths is not yet clear.

For the high frequency $`\omega _h`$ and laser intensities $`I>3\times 10^{15}\text{Wcm}\text{-2}`$ the average charge state in the cluster overtakes the charge state of the single atom (see Fig. 6). It was therefore interesting to check whether under those conditions a wave packet also forms at the higher frequency $`\omega _h`$. This is indeed the case. However, due to the rapid ionization the wave packet dynamics lasts only a few laser cycles (or less).

### III.6 Influence of ionic motion on the wave packet formation

So far the ion mass was set to $`M=1311836`$ (Xe atom). This high mass lead to no appreciable ionic motion in the 1D cluster during the pulse. Even in the “worst case” where all active electrons are removed at $`t=0`$ so that the ionic charge $`Z=4`$ is not screened at all, the $`N_{\text{ion}}=9`$-cluster radius increases only from $`32`$ a.u. to $`33.8`$ a.u. within the 37 fs pulse duration. However, the corresponding spherical 3D cluster with $`(N_{\text{ion}}1)^3`$ Xe-ions of charge state $`Z=4`$ doubles its initial radius within the same time doublingformula . Clearly, the Coulomb explosion is underestimated in the 1D model because in 3D an ion sitting at the cluster surface “sees” not only the charges of the ions aligned along the laser field polarization direction but essentially a charged sphere containing all the other $`(N_{\text{ion}}1)^31`$ ions.

Hence, in order to allow the 1D ion chain to explode similarly to the corresponding 3D cluster one has to reduce the ion mass. To obtain the results shown in Fig. 8 $`M=1836`$ was set. In the “worst case” introduced in the previous paragraph the radius of the 1D cluster now expands from $`32`$ to $`130`$ a.u. within $`37`$ fs. However, in the early stage when the wave packet formation inside the cluster takes place the charge states of the ions are still low and the cluster radius remains close to its original value. This is clearly visible in Fig. 8a where between $`t=300`$ and $`t=800`$ the bouncing wave packet can be easily identified. In Fig. 8b,c the dipole $`x_{\text{inner}}`$ and the number of electrons inside the cluster $`N_{\text{inner}}`$ are presented, as in Fig. 3. For comparison, the results for a computer run with immobile ions is included (blue and broken curves). The differences are small. Only the ionization at the end of the laser pulse is slightly reduced in the case of mobile ions. This is expected because the ions gain their kinetic energy at the expense of the electrons.

At higher laser intensities the higher charge states are generated earlier. In these cases the Coulomb explosion is more violent but the wave packet mechanism contributes little to the total ionization of the cluster anyway. In Fig. 9 the ionization and Coulomb explosion scenario for the highest laser intensity $`I=1.14\times 10^{16}\text{Wcm}\text{-2}`$ considered in this paper is presented for the same cluster parameters as in Fig. 8. The outer shell of all atoms is depleted rapidly already during the first laser cycle both for the cluster and the isolated atom. On average three quarters of the electrons in the next shell are depopulated during the remainder of the laser pulse in the case of the cluster with mobile ions. The ionization probability of the isolated atom is slightly less. The ionization of the cluster with heavier ions is, again, higher.

Due to the rapidly increasing charge states of the ions no formation of an electron wave packet moving $`\mathrm{\Delta }\varphi =\pi /2`$ out of phase with respect to the driving laser field for several laser cycles can be inferred in Fig. 9a. Instead, the typical motion $`𝑬(t)`$ of free electron density in the laser field $`𝑬(t)`$ seems to dominate the dynamics. It is interesting to observe that the removal of the electrons in the first shell is even more efficient for the isolated atom than in the cluster (see Fig. 9c around $`t=180`$). However, shortly after “cracking” the second shell the cluster ionization overtakes the single atom result. At that time the $`x_{\text{inner}}`$-plot in Fig. 9b reveals two extrema (marked by arrows) which are approximately $`\mathrm{\Delta }\varphi =\pi /2`$ out of phase with the driving field. Thus, although a dephased wave packet dynamics over several cycles cannot develop in such intense fields the stronger ionization of the cluster still relies on the fact that the electron density inside the cluster moves temporarily with the appropriate phase lag necessary for efficient absorption.

## IV Conclusion

The ionization dynamics of a one-dimensional rare gas cluster model in intense and short laser pulses was investigated by means of time-dependent density functional theory. An electron wave packet dynamics was found to build up inside the cluster when the laser intensity $`I=\widehat{E}^2`$ was sufficiently high for modest inner ionization but not that high that all electrons of an atomic shell are freed within a few laser cycles. The electron wave packet is driven into resonance with the laser field through the collisions with the cluster boundary. The phase lag between the bouncing electron wave packet and the laser field then is $`\pi /2`$ so that the absorption of laser energy is particularly high. The fastest electrons in the wave packet escape from the cluster (outer ionization). The electric field of the bouncing electron wave packet adds up constructively to the laser field, thus enhancing inner ionization. This effect was called dynamical ionization ignition. It is a robust phenomenon with respect to the cluster size and, since it occurs during the early ionization stage, it is not affected by ionic motion. However, with increasing laser frequency (keeping the laser intensity fixed) the mechanism is less efficient.

We expect the wave packet scenario being valid also in real, three-dimensional rare gas clusters. Due to the spherical cluster-vacuum boundary the wave packet should assume a sickle-like shape in that case. However, in order to verify this, studies with higher-dimensional cluster models will be pursued in the future.

###### Acknowledgements.

This work was supported by the Deutsche Forschungsgemeinschaft (D.B.) and through the INFM advanced research project Clusters (A.M.). The permission to run our codes on the Linux cluster at PC<sup>2</sup> in Paderborn, Germany, is gratefully acknowledged.