UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
DYNAMICAL EVOLUTION OF ASTEROIDAL AND COMETARY PARTICLES AND THEIR CONTRIBUTION TO THE ZODIACAL CLOUD BY JERCHYI LIOU DISSERTATION PRESENTED TO THE GRADUATE SCHOOL THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1993 4jwfl, twJAiw o To my parents, brother, sister in law, and my beloved wife for their love and support. m^^wm.sf~~q~~ ACKNOWLEDGMENTS I would like to thank my adviser, Dr. Stan Dermott, and the other members in my committee, Dr. Campins, Dr. F. E. Dunnam, and Dr. H. Smith for their support and guidance through out the work of this dissertation. Special thanks go to Dr. A. Smith for reading the thesis and providing me with useful suggestions. Special thanks to Dr. Yulin Xu for helping me understand SIMUL and providing me with the observational data from IRAS. I would like to thank Dr. B. Gustafson for many useful discussions. I would also like to thank Dan Durda for his comment on the collisional evolution of the asteroids. Finally, I would like to thank my fellow graduate students, Jaydeep Mukherjee for his help over the years, Sandra Clements for showing me how to be a good observer at RMH Observatory, David Kaufmann and Damo Nair for making me the best backgammon player in the department. TABLE OF CONTENTS ACKNOWLEDGMENTS LIST OF TABLES . LIST OF FIGURES ABSTRACT iii v S . . . . . . . v n . .. .. .. .. ... S. S. S XV CHAPTERS INTRODUCTION Li~l * * * * * 1 The Zodiacal Clo Objective of this Dissertation 2 SECULAR PERTURBATION THEE C S . S S S S S S S 3 S S S S S S S S S S S S C S S S S 5 6J Loworder Secular Perturbation Theory SS S S S S . S .. 6 The Solar System. Small Bodies in the Solar System Numerical Approach Comparison between the Analytical and Numerical Calculations SunJupiterSaturn System S S S S S S S S S S 10 . 14 S C 5 5 5 5 5 5 S* S l. T14 Asteroidal Particles in the Solar System. . . S S S S 15 ORBITAL EVOLUTION OF ASTEROIDAL PARTICLES Radiation Pressure , PoyntingRobertson Drag, and the Corpu scular Drag. Static Theory and Dynamical Theory S S S S S S S S S 28 Numerical Comparisons SS . C 32 S S S S . S S S S S S 5 6 5 5 5 6 S. S . 12 . C S S 5 25 25 Passage Through Resonance . 4 4 C 4 4 C 9 4 C 34 Orbital Evolution of Dust Particles from the Major Hirayama Asteroidal Families S36 9 pm Diameter Particles. . ...... . 38 4 pm and 14 pm Diameter Particles. . . . 41 NeEarth Dust Grains. 4 ORBITAL EVOLUTION OF COMETARY PARTICLES Orbital Evolution of Encketype Dust Particles .. Gauss' * a a a a .77 .* a a 77 Equations Can Comets be the Source of the Solar System Dust Bands? Distribution of Encketype Dust Particles in the Solar System . 8 STRUCTURE OF THE ZODIACAL CLOUD ......... 106 Contribution from Asteroidal Dust Particles . 108 Contribution from Encketype Dust Particles . . . 110 Discussion BIBLIOGRAPHY . 4 a a a C . S S S S 111 129 BIOGRAPHICAL SKETCH . . . . . . . f f .. . ..80 132 LIST OF TABLES Comparison between numerical and analytical results for the real Jupiter and Saturn . . S S S S S S 4 S S S 17 Comparison between numerical and analytical results for small m, e, I, Jupiter and Saturn 18 LIST OF FIGURES Variation in h, k, p, and q of Jupiter subscriptt J) with time obtained from numerical integration of the SunJupiterSaturn system 19 Variation in h, k, p, and q of Saturn subscriptt S) with time obtained from numerical integration of the SunJupiterSaturn system 20 Distribution in inclination space of Koronislike particles at t=0 and t=280,000 years . .. .*..ft ... .... f ft2 1 Distribution in eccentricity space of Koronislike particles at t=0 and t=280,000 years, where w is the longitude of pericenter Variation in the forced p and q over a period of 300,000 years. dots are from direct numerical integration while the solid line is from secular perturbation theory . . 23 Variation in the forced h and k over a period of 300,000 years. dots are from direct numerical integration while the solid line is from secular perturbation theory .. 24 Variation in the total subscriptt t), forced subscriptt f), and proper subscriptt p) eccentricities and longitudes of pericenter from t=0 to Variation in the osculating inclination vectors of 9 um diameter Eos type dust particles as they spiral towards the Sun from 1.8 AU Each panel shows the vectors at various mean distances from the Sun Variation in the osculating eccentricity vectors of 9 pm diameter Eos type dust particles as they spiral towards the Sun from 1.8 AU. Each panel shows the vectors at various mean distances from the Sun 22 Proper inclination as a function of semimajor axis of 9 pm diameter Eostype dust particles as they spiral towards the Sun from 1.8 AU Error bars are the dispersions in the radii from the mean radius obtained from the least squares fit . S S S S S 50 Comparison of the forced inclination obtained from the dynamical theory (bold line) and from numerical integration (points) of 9 pm diameter Eostype dust particles as they spiral towards the Sun from 1.8 AU. Error bars are the dispersions in the forced inclinations from the mean forced inclination obtained from the least squares fit 51 Comparison of the forced node obtained from the dynamical theory (bold line) and from numerical integration (points) of 9 pm diameter Eostype dust particles as they spiral towards the Sun from 1.8 AU. Error bars are the dispersions in the forced nodes from the mean forced node obtained from the least squares fit . 52 Comparison of the forced eccentricity obtained from the dynamical theory (bold line) and from numerical integration (points) of 9 pm diameter Eostype dust particles as they spiral towards the Sun from 1.8 AU. Error bars are the dispersions in the forced eccentricities from the mean forced eccentricity obtained from the least squares fit . . . S S S S S S S 53 Comparison of the forced longitude of pericenter (Wforced) obtained from the dynamical theory (bold line) and from numerical integration (points) of 9 pm diameter Eostype dust particles as they spiral towards the Sun from 1.8 AU. Error bars are the dispersions in the forced pericenters from the mean forced pericenter obtained from the least squares fit S . . 54 Forced and proper eccentricities as functions of the semimajor axis. The two straight reference lines have slopes of 1.25 3.10: Variation in the osculating inclination vectors of 9 pm diameter Eos particles as they spiral towards the Sun. Each panel shows the vectors at various mean distances from the Sun. *. C 55 56 Variation in the osculating eccentricity vectors of 9 jpm diameter Eos particles as they spiral towards the Sun, where w is the longitude of pericenter. Each panel shows the vectors at various mean distances from the Sun. 3.12: S S S S S S C * C S 57 Proper inclination as a function of semimajor axis of 9 jm diameter Eos dust particles as they spiral towards the Sun. Error bars are the dispersions in the radii from the mean radius obtained from the least squares fit 3.13: S. S S S S 58 Comparison of the forced inclination obtained from the dynamical theory (bold line) and from numerical integration (points) of 9 jm diameter Eos dust particles as they spiral towards the Sun. Error bars are the dispersions in the forced inclinations from the mean forced inclination obtained from the least squares fit . S S 59 3.14: Comparison of the forced node obtained from the dynamical theory (bold line) and from numerical integration (points) of 9 pm diameter Eos dust particles as they spiral towards the Sun. Error bars are the dispersions in the forced nodes from the mean forced node obtained from the least squares fit S S S C S 4 5 60 3.15: Comparison of the forced eccentricity obtained from the dynamical theory (bold line) and from numerical integration (points) of 9 /m diameter Eos dust particles as they spiral towards the Sun. Error bars are the dispersions in the forced eccentricities from the mean forced eccentricity obtained from the least squares fit 3.16: S S S C f 61 Comparison of the forced longitude of pericenter (Wforced) obtained from the dynamical theory (bold line) and from numerical integration (points) of 9 pm diameter Eos dust particles as they spiral towards the Sun. Error bars are the dispersions in the forced pericenters from the mean forced pericenter obtained from the least squares fit 3.11: . 62 Comparison of the forced eccentricity obtained from the dynamical theory (bold line) and from numerical integration (points) of 9 pm diameter Eos dust particles as they spiral towards the Sun. Error bars are the dispersions in the forced eccentricities from the mean forced eccentricity obtained from the least squares fit. The masses of Jupiter and Saturn have been reduced to 0.1 of their real values 3.18: Variation in the osculating inclination vectors of 9 pm diameter Eos particles, in the real solar system, as they spiral towards the Sun. Each panel shows the vectors at various mean distances from the Sun .................................... 64 3.19: Variation in the osculating eccentricity vectors of 9 pm diameter Eos particles, in the real solar system, as they spiral towards the Sun, where w is the longitude of pericenter. Each panel show vectors at various mean distances from the Sun. St 65 3.20: Forced inclination as a function of semimajor axis for 9 pm diameter Eos particles in 1983. These are the results of integrating waves of particles from the past such that they end up at various distances from the Sun in 1983 0 5 S S S 9 S S 5 66 3.21: Forced node as a function of semimajor axis for 9 pm diameter Eos particles particles in 1983. These are the results of integrating waves of from the past such that they end up at various distances from the Sun in 1983 3.22: . 67 A vertical crosssection of a model dust band produced from Eos particles alone 3.23: 68 Forced inclination as a function of semimajor axis for 4 different systems. 4 G.P. stands for four giant planets, Jupiter, Saturn, Uranus, and Neptune. , E, M, stand for Venus, Earth, and Mars, respectively 3.24: Forced inclinations as a function of semimajor axis for 4, 9, and 14 rpm diameter Eos particles in 1983. These are the results of integrating waves of particles from the past such that they end up at 'mr4 nr i n A Clntnin 4n rnrv. *4l9 Cnn 4 n. 1 AQ '2 rn . . . . . . . . 69 3.25: Forced nodes as a function of semimajor axis for 4, 9, and 14 pm diameter Eos particles in 1983. These are the results of integrating waves of particles from the past such that they end up at various distances from the Sun in 1983 3.26: . a a a a a a a a a a 7 1 Distribution in inclination space of particles with 4 to 9 pm diameters from Eos, Themis, and Koronis families in the Earthcrossing region. The outer circle contains Eos particles while the inner group contains Themis and Koronis particles a a a a a 72 Seasonal variation of the relative numbers of asteroidal particles that will encounter the Earth. These are particles from the Eos, Themis, and Koronis families with diameters in the range 4 to 9 pm 3.28: .. 73 Relative numbers of particles that will encounter the Earth in November. These are particles with diameters in the range 4 to 9 pm originating from Themis (bold line), Eos (bold line), and all other asteroids (thin line). 3.29: a a a a a a a 4 a a a a a 7 4 Variation in eccentricity and inclination of 9 pm diameter Eos particles as they approach and pass the Earth. upper right hand corners of the panels The numbers at the indicate the time example, 19833100 means 3100 years before 1983 3.30: a a a a a /75 Variation in the semimajor axes of fifty 9 pm diameter Eos particles as they approach and pass the Earth. A few of them are trapped in mean motion resonances with the Earth. While they are trapped, their semimajor axes remain constant and produce horizontal lines in the diagram. Sa a a a a a a a 76 Variation in the osculating eccentricity vectors of 9 pm diameter Encketype dust particles as they spiral towards the Sun, where cw is the longitude of pericenter. panels is 600 years . The time interval between consecutive 90 Variation in the osculating inclination vectors of 9 pm diameter Encketype dust particles as they spiral towards the Sun. interval between consecutive panels is 600 years The time * S 4 a a 9 1 The variation in (a,e) for every 1000 years of one wave of 9 pm diameter Encketype dust particles as they spiral towards the Sun. A few of them have their semimajor axes increased up to 4 AU due to the process of release from the comet Relative rate of change of semimajor axi 92 s (in arbitrary units) as a function of eccentricity due to the action of drag force 93 Distribution in (a,e) space for 9 pm diameter Encketype dust particles in 1983 from the second model. These are the results of integrating particles released continuously (once every 300 years) starting from the year (19835,400). S94 Distribution in (a,I) space for 9 pm diameter Encketype dust particles in 1983 from the second model. These are the results of integrating particles released continuously (once every 300 years) starting from the year (19835,400). . . S S S S 95 Longitude of the ascending node (Q) vs. longitude of pericenter (?w) for 9 pm diameter Encketype dust particles in 1983 from the second model. These are the results of integrating particles released continuously (once every 300 years) starting from the year (19835,400) . 96 Distributions of the osculating eccentricity vectors in different semimajor axis ranges for 9 pm diameter Encketype dust particles in 1983, from the second model. represented by The longitude of pericenter is These are the results of integrating particles released continuously (once every 300 years) starting from the year (19835,400) S* 97 Distributions of the osculating inclination vectors in different semimajor axis ranges for 9 pm diameter Encketype dust particles 1983, from the second model. These are the results of integrating particles released continuously (once every 300 years) starting from the year (19835,400). 98 4.10: Distribution of the osculating inclination vectors in different semimajor axis ranges for 9 pm diameter Encketype dust particles in 1983, from the first model. These are the results of integrating particles released continuously (once every 300 years) starting from the year (19835,400) . . . . . . 99 Distribution of the inclinations in 1983 of 9 pm diameter 4.11: Encketype dust particles from the first model. 4.12: a 100 Distribution of the inclinations in 1983 of 9 pm diameter Encketype dust particles from the second model 4.13: . 101 Relative number of orbits per unit semimajor axis for 9 pm diameter Encketype dust particles 4.14: . * 9 4 5 102 Relative number of orbits per unit semimajor axis for 9 pm diameter Eos dust particles. 4.15: . S 4 103 Forced eccentricity as a function of semimajor axis for 9 pm diameter Encketype dust particles in 1983. Error bars are the dispersions in the forced eccentricities from the mean forced eccentricity obtained from the least squares fit 4.16: . . 104 Proper eccentricity as a function of semimajor axis for 9 pm diameter Encketype dust particles in 1983. Error bars are the dispersions in the radii from the mean radius obtained from the least squares fit Observed shape of the zodiacal cloud obtained from IRAS data in the 25/pm waveband with an 900 elongation angle (SOP406) 105 118 Observed latitudes of the peak fluxes of the zodiacal cloud around the sky obtained from IRAS data in the 25pm waveband with an 900 elongation angle (SOP406) Comparison between IRAS' .. . . . 5 9 5 119 observations (bold line) and the allasteroid model (thin line) at SOP406, 900 elongation angle, leading direction . . . . a ... 120 Comparison between IRAS' observations (bold line) and the allasteroid model (thin line) at SOP392, 970 leading direction . elongation angle, S. S S S * . 121 Latitudes of the peak fluxes around the sky from the allasteroid model (thin curves) and the best fit observational data (bold curves) S S S S S S S S S S S 5 S S S 122 The zodiacal cloud at SOP392; the bold lines are the IRAS observations. The thin lines are from the allasteroid model (top), the first all comet model (central), and a combination of asteroid and comet model (bottom) . . . S S S 123 The zodiacal cloud at SOP406; the bold lines are the IRAS observations. The thin lines are from the allasteroid model (top), the first all comet model (central), and a combination of asteroid and comet model (bottom) Latitudes of the peak fluxes arou nd the sky for the 27% asteroi124 nd the sky for the 27% asteroid 125 plus 73% comet model Comparison between the second allcomet model (thin lines) and the observations (bold lines) at two different SOP positions . 126 5.10: Inclination distribution of Jupiter Family short period comets . 127 5.11: Distribution in (e, I) space of Jupiter Family short period comets which have their apocenters inside the orbit of Jupiter 128 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy. DYNAMICAL EVOLUTION OF ASTEROIDAL AND COMETARY PARTICLES AND THEIR CONTRIBUTION TO THE ZODIACAL CLOUD By JerChyi Liou December, 1993 Chairman: Major Department: anley F. Dermott Astronomy Asteroids particles comets that populate are the zodiacal major cloud. sources how the interplanetary particles evolve and in exactly what proportion asteroids and comets contribute to this cloud are questions that have never been convincingly answered. In this dissertation we try to answer these questions studying the dynamical evolution asteroidal and cometary particles and constructing a physical model to compare with observations from the Infrared Astronomical Satellite (IRAS). We first study the orbital evolution of asteroidal and cometary particles under the influence of planetary perturbations, radiation pressure, PoyntingRobertson drag, t 1 A 1 11 ** ,m .... secular perturbation theory for asteroidal particles. We use a numerical integrator, RADAU, on an IBM ES/9000 supercomputer to study orbital evolution thousands of particles, from both asteroids and comets, systematically, from their origins to a heliocentric distance of around 0.1 AU. We then analyze the results and apply the distributions of the orbital elements from the dynamical study to a three dimensional model, SIMUL, which generates the observed flux from a given number of orbits with a known distribution of orbital elements, and compare the results with the IRAS observations. We conclude that it is possible to construct a zodiacal cloud model to fit the shape of the observed cloud by a combination of approximately 27% asteroidal particles and 73% cometary particles. CHAPTER INTRODUCTION The Zodiacal Cloud The zodiacal cloud is a feature of the night sky that can be seen by the naked eye just before sunrise and just after sunset. brightness increases towards the ecliptic plane and is roughly symmetric with respect to that plane. Since the all sky survey in 1983 by the Infrared Astronomical Satellite (IRAS), the structure of the zodiacal cloud in four infrared wavebands (at 12, 25, 60, and 100 pm) has been known in great detail, including the overall shape of the cloud and the latitudes of the peak fluxes around the sky (e.g. Hauser et al., 1985). The infrared flux is mainly due to thermal emission by micronsized or larger dust particles. What are the origins of these particles? How do they form and how do they evolve? What are the global and finescale structures of the cloud? Is there any time variation in the cloud structure? The observational data from IRAS has provided us with a detailed look at the structure of the cloud, not only increasing our understanding of the cloud but also setting critical constraints on the models which will be used to explain the formation and structure of the cloud. Previous models of the structure of the zodiacal cloud are based on some assumed n nll r t r A n 'ii 4tit Cl.nno *di\tPlt A a.. I A : ,^1.. na n.A an1 ..2 ^1_ la.. Giese et al., 1985, 1986, Kneissel and Mann, 1991). These models do not consider the formation, source, and evolution of the particles; they are restricted to a description of the observations rather than an explanation of the physical structure and evolution of the cloud. When compared with IRAS observations, none of these models are capable of explaining several key features in the data. The dynamical evolution of individual micronsized dust particles in the solar system been studied previously Gustafson et al., 1987 Jackson Zook, 1992). A more recent study Dermott et al. (1992) has shown that it is possible to study the orbital evolution of dust particles coming from different sources systematically by incorporating the concept of the socalled proper forced elements from the secular perturbation theory (to be explained in Chapter However, the work of Dermott et al. needs some improvement; in particular, it did not, for example, include the effects of solar wind corpuscular drag and the changes in the orbital elements of dust particles due to the process of release. we will show, these are effects that significantly affect the orbital evolution of dust particles. In summary, our aim is to study the orbital evolution of dust particles coming from different sources systematically, and including important dynamical effects. the structure of a model zodiacal cloud based on kind of dynamical study and compare our results with the IRAS observations. 3 Objective of this Dissertation dissertation, we consider dynamical evolution particles, systematically, from source to sink, as described by Dermott et al. (1992) and use the results from the dynamical study to construct a zodiacal cloud model. know that asteroids and comets are the two major sources of the interplanetary dust particles that populate the zodiacal cloud. But in what way these particles evolve and in exactly what percentage they contribute to this cloud are questions that have not been answered before. These are the questions that we would like to answer. In Chapter 2 we examine the classical secular perturbation theory which describes the orbital evolution of planets due to their mutual gravitational perturbations and the orbital evolution of asteroidal particles due to the gravitational perturbations of the planets. We first study the evolution of the solar system numerically and find the eigenfrequencies, amplitudes, and phases associated with the planets in the system that describe the variations of their eccentricities and inclinations and compare them with predictions from secular perturbation theory. We also examine numerically the orbital evolution of asteroidal particles in a (SunJupiterSaturn) system and, again, compare the results with those from secular perturbation theory. In Chapter 3 we revise the classical secular perturbation theory to include the effect of drag. The idea that we use was first developed by Dermott et al. in 1992. complete the theory and apply this dynamicall" theory to the evolution of asteroidal dust particles starting from several different heliocentric distances until they reach thp nnPr rt nfrt nrf thp crliar cnotprm x1z7 thsn r*mnrpa thpcp' Qnlhrtrlnl lrailto u;th th numerical results. We next integrate numerically the orbits of thousands of 4, 9, and 14 pm diameter asteroidal dust particles systematically from their origin until they reach about 0.1 AU from the Sun in the real solar system. All the planetary perturbations (except those of Mercury and Pluto), radiation pressure, Poynting Robertson light drag, and solar wind corpuscular drag effects are included in our calculations. We analyze the data and obtain the distributions of their orbital elements as functions of semimajor axis, and also examine the orbits of these dust particles when they approach the Earth to gain information useful for future Earthorbiting dust collection facilities. In Chapter 4 we study the evolution of Encketype dust particles using procedures similar to those in Chapter 3. This is the first time that this has been attempted for particles in cometary type orbits. We construct two models in which the dust particles are given different initial inclination distribu numerically as they spiral towards the Sun. to those in Chapter We follow their orbital evolution The data are analyzed in a similar way 3 in order to obtain the distributions of their orbital elements as functions of semimajor axis. The difference in the initial inclinations turns out to be a key element in constructing the zodiacal cloud model. We also show that, because of their high eccentricities, comets cannot be the source of any of the solar system dust bands discovered by IRAS. We use the results from Chapters 3 and 4 to construct several physical zodiacal cloud models (Chapter 5). A threedimensional model, SIMUL, is used to calculate 5 cloud and latitudes of the peak fluxes from different zodiacal cloud models with the observations from IRAS. Finally we show that by combining cometary dust particles and asteroidal dust particles, shape of the zodiacal cloud. we can construct a model that describes the observed We also discuss alternative models and some unsolved problems. CHAPTER 2 SECULAR PERTURBATION THEORY In this chapter we will first summarize the theory of secular perturbation, which is the classical method to study the longterm orbital evolution of planets as well as minor planets in the solar system. We apply this theory to a (SunJupiterSaturn) system to find the eigenfrequencies, amplitudes, phases associated planets in the system and compare the results with those from the direct numerical integration. We then use the results from the numerical integration to calculate the initial forced elements for asteroidal type particles at 2.6 AU and follow their orbital evolution in this (SunJupiterSaturn) system both numerically and analytically and compare the results. Loworder Secular Perturbation Theory The Solar System The classical way to study the long term rates of change of orbital elements, eccentricities (e), longitudes of pericenter (w), inclinations (I), and longitudes of ascending node (St) of small bodies in the solar system without any mean motion resonance involved is by using secular perturbation theory (Brouwer and Clemence, 1961, Dermott et al., 1986). It is based on the assumption that the eccentricities and 7 and Q of a small body are largely determined by the secular terms in its disturbing function (to be defined later). These are terms independent of the mean longitudes of any objects involved. Considering the solar system with the Sun at the center, the Lagrange equations which describe the rates of change of the orbital elements of any given planet are n na A OR  1 e 2 \/~ e2R na2e tanI ,OR na2V12 OA na2e OR S) na2 vT sin I 01 (2.1) na2 ~/ 2 sinI l V'STe2R tan I 2 na2e na2 .527 dt nada I /J C VY) OR na2e tan 11 na21/F7e2 where n23 n a = G(M + m) (2.2) E+ and G, M, m, n, R are the universal gravitational constant, mass of the sun, mass of the given planet, the mean motion of the given planet, and the disturbing function, respectively. The disturbing function is the potential acting on this given planet due to the existence of other planets in the system. The semimajor axis, eccentricity, inclination, longitude the ascending node, longitude pericenter, mean 1 an n.. 4...,) a nH tAn ..a tan n a 4 1. n n F fl .. ..,,,..,l.. ..1 9e + ,,^, \ 8 In order to avoid singularities that arise in Equation 2.1 for small eccentricities and inclinations, we define four quantities (h,k,p,q) to replace e, w, I, and Q: = e sin w = e COS W (2.3) =tan I = tan I sinf cosfl . (2.4) If the eccentricities and inclinations of the planets are small, we can ignore the second and higher order terms in e and I in Equation 2.1 and the resultant Lagrange equations involving e, w, I, and t become 1 ORj nya^ Okj k  ' 1 nja3 (2.5) ORj 9h, . 1 OR, p} = nja 9qj 2q q. = 1 nya 2 nja1 (2.6) ORj dpj By using similar assumptions, the disturbing function to the second order in e and I experienced by the jth planet has the following form (Dermott and Nicholson, 1986): Rj =nja 2us mi f ol .. .b(l) mrc ,8 3/2( m 1 2)" 1 (2) 4 ji ,' 3/2 e3e cos (wj Wi 2. m ri 1 (1 7n 33/2 3 aj I I 2 =nj a3 1 k) + < A,,(h + kc) + 2 I3 ! +njaG 1Bjj1(p4 + q2) + Aji(hjhii + kjki) ( P ) where j is not equal to i and Si/aj 1aj aji/ag A3~ = n .7 .74 (internal (external (internal (external perturber) perturber) perturber) perturber) mi i _ (1 )3 M&JbOJb3/ 2 (2.9) (2.10) (2.11) 1M^ b(2) 4M J1 3/'2 1  n 4 j mi ..(1) M O3 3 3/2 (2.12) S1 mi ( ii =~ nl jioj 3/2 and b1') 3/2' b arethesocalledLaplacecoefficients(e.g.2) 3/2 are the socalled Laplace coefficients (e.g. Brouwer and Clemence, 1961). Mercury and Pluto are not included in our calculation due to their small masses. The solutions of hj, kj, qj are hIj klj  ej sin(git+ fli) (2.13) cos(git + ypi) Ij, sin(fit + yi) (2.14) pj = qi j Ijicos(fjt+ 7i) (2.8) Otj  aji = A ^in  B = r\ ** 10 where gi and fi are the eigenfrequencies of the matrices A and B whose elements are defined by Ajj, Aji, jj, and Bji above. The amplitudes, eji, Iji, and phases, %i, are determined from the initial conditions of the planets. Equations 2.13 and 2.14 are the results from the classical loworder secular per turbation theory. These two equations describe analytically the respective variations in the eccentricity, inclination, longitude of pericenter, and longitude of ascending node of the jh planet. They are based on the assumption that the eccentricities and inclinations of all the objects involved are small. When there is no mean motion res onance phenomenon involved, this theory describes fairly accurately the longterm variations in e, zt, I, and Q of the planets. Small Bodies in the Solar System For small bodies moving under the gravitational perturbations of planets, can follow similar procedures to find the longterm rates of change in their (h,k), and (p,q). If we use A to represent the nodal precession rate of the small body subscriptt s), (2.15) where n, mi, and M are the mean motion of the small body, the mass of the it planet, and the mass of the Sun, respectively, then the solutions for (h,k) and (p,q) are h = epsin(At + p) + ho k = ep cos(At + f3) + ko Sp = Ip sin(At + 7) + po (2.16) /) 17\ A= aj .(1) A 4 2 M ^^3/ 1=1 where ha = 7 i= =1 i= i=1 Po = qo = i sin(git + /i) (2.18) Si cos (it + /3i) A + sin(fit + yi) A+ f; (2.19) +i fcos(fit + 7i) A+ fi  mi (2) M>_c0^1s l3/2 ti (2.20) Pi = _ lb Iii M 3/2 and gifi, fe, eii, and pi, 7i are the eigenfrequencies, amplitudes, and phases associated with the planets of the solar system from the previous section. In the previous equations, ep, Ipt, and 7 can be determined from the initial conditions of the small body. The forced eccentricity and inclination are defined as ef = If= (2.21) The facts that (h,k), and (p,q) are decoupled from each other and, for each of them, we can use the proper and forced elements (their characteristics under the influence of drag force will be described in Chapter 3) instead of the osculating elements to describe their evolution make it possible to study dust particles coming from the same source systematically (Dermott et al., 1985). The precession term, A, depends on the masses of the planets if the semimajor aC a C C C C + qo2 * fi r 1 A due to the existence of planets. If all the planets have circular orbits and the same inclinations, then there will be no forced eccentricity and longitude of pericenter (due to zero eccentricities of planets) and the orbit of the minor planet will process along the common plane of the planets (due to the same inclinations of planets). However, this is not the case in the real solar system. The forced components in Equations 2.16 and 2.17 will be present due to the eccentricities and different inclinations of the planets' orbits. The forced components introduce a plane of symmetry (different from any of the planets' orbital planes), with respect to which the orbit of the minor planet processes, along with the displacement of the pattern of the orbit of the minor planet away from the Sun (Dermott et al., 1985). The equations of (h,k) and (p,q) describe the time variation in orbital elements of small bodies in the solar system. This is a result from the loworder theory. When the eccentricity or inclination of the small body becomes larger, the loworder theory can no longer be applied to study the evolution of this minor planet' orbit. In this case, direct numerical integration is the only way to study its orbital evolution. In the next section we will test this loworder theory on the real solar system and compare the results with those from direct numerical integration. Numerical Approach After the invention of high speed computers, astronomers found another powerful to study celestial mechanics: numerical integration. Many effects that can mtegration. For example, if we include the gravitational forces due to the Sun and planets and the radiation pressure force as well as the drag terms in the equation of motion of a dust particle and numerically integrate the equation, passage through mean motion resonance, planetary scattering, and trapping in mean motion resonance associated with planets can be seen from the numerical output of the orbit of this particle. However, the numerical error built up in the process needs to be considered carefully. Besides, not all problems can be solved numerically. In some cases, for example, when the number of particles involved becomes too large or the integration time becomes too long, even the fastest machine in the world is not sufficient to obtain the desired results. The numerical integrator we chose is RADAU. It was developed by Everhart (1985). The major computer we use to run the numerical integration is an IBM supercomputer ES/9000 in Northeast Regional Data Center (NERDC) at the University of Florida. IBM ES/9000 is equipped with a vector facility. For each processor, it can take up to 256 objects and perform the same operation at the same time. In order to improve the performance, we have vectorized RADAU. In addition to IBM ES/9000, we have also used IBM RS/6000 350, IBM RS/6000 340, IBM RS/6000 220, and several Sun workstations. The numerical integration we describe in this chapter include two parts. first part is the integration of the solar system with only two planets, Jupiter and Saturn, with their real masses and orbital elements. This is a simple threebody gvgtem Wp hnve. sl1 and inclinations of Jupiter and Saturn. Both numerical results were analyzed and compared with analytical predictions from secular perturbation theory. second part is the integration of the orbital evolution of 100 asteroidal type particles in the (SunJupiterSaturn) system. Again, the comparison between analytical and numerical results was made. Comparison between the Analytical and Numerical Calculations SunJupiterSaturn System We first test the secular perturbation theory on the (SunJupiterSaturn) system. This threebody system was integrated from 1983 backward for five hundred thousand years. The variations in h, k, p, and q of Jupiter and Saturn are shown in Figure 2.1 and Figure 2.2. We analyzed the data by using the Fast Fourier transform method to find the eigenfrequencies, amplitudes, and phases of the system from the power spectra and compared the results with the analytical calculation. shown in The comparison is Table 2.1. As we can see from the table, the two results agree reasonably well. there are certain discrepancies for the eigenfrenquencies, gi and g2. are 15.2% and 16.8%, respectively. motion resonance between Jupiter an However, The differences A possible reason for this might be the mean d Saturn. These two planets are very close to 5:2 mean motion commensurability (a third order resonance). The analytical third order resonance theory has not been worked out, and its development is beyond the 15 We can also check this mean motion resonance effect by reducing the masses, eccentricities, and inclinations of Jupiter and Saturn and compare the results again. The comparison is in Table The masses of Jupiter and Saturn are 10% of their original masses. Their eccentricities and inclinations are ej = 0.0004948, Ij = 0.003273 degree, es = 0.0005085, and Is = 0.008872 degree. The largest difference is only 3.3 All the other differences are less than 1%. The secular perturbation theory predicts very accurately the results from the numerical integration in this case. From these comparisons we can conclude that the loworder secular perturbation theory does describe the behavior of planets accurately except when a mean motion resonance between planets is involved. We will further examine how good this theory is when we apply it to the asteroidaltype particles in the next section. Asteroidal Particles in the Solar System In the calculations described in this section we put 100 asteroidaltype particles in the SunJupiterSaturn system and integrate their evolution for 280,000 years. Again, we compare the results with the analytical calculation from secular perturbation theory. The semimajor axes for the particles are chosen to be 2.6 Astronomical Units (AU). This is to avoid several mean motion resonance locations Jupiter at around 3 AU. We use the "particles in a circle" method (Dermott et al., 1985, 1992) to set up the initial orbital elements of the particles. Their proper eccentricities and inclinations are 0.049 and 2.12 degrees, respectively, which are the proper elements a' rr / , r\ I 1 randomly chosen between 00 and 360 . The initial forced elements at 2.6 AU are calculated by using the eigenfrequencies, amplitudes, and phases from the numerical integration of the (SunJupiterSaturn) system from the previous section. The positions of particles in the (p,q) and (h,k) phase spaces at t = 0 and t = 280,000 years are shown in Figures 2.3 and 2.4. Due to the difficulty in generating w in the diagram, we use w to represent the longitude of pericenter in all diagrams. fact that after 280,000 years all the particles still remain in pretty welldefined circles in both diagrams shows that the method of using the forced and proper elements to describe the evolution of asteroidaltype particles is a good one. The variations in the forced elements are shown in Figures 2.5 and 2.6. dots are results from the direct numerical integration while the solid lines are results from the analytical calculation from the secular perturbation theory. are about 0.20 The differences in the forced inclination and less than 0.01 in the forced eccentricity. These results are actually not too bad. This shows that when the particles are not in mean motion resonance with any planet, the secular perturbation theory does describe the evolution of asteroidaltype particles to a certain accuracy. However, when small particles are spiralling towards the Sun due to PoyntingRobertson and corpuscular drag, this classical secular perturbation theory does not apply any more. chapter we use the theory of Dermott et al. In the next (1992) to include the effect of drag. 1 0 26  0 0 18 z 0 CQtf cd& *~a iO n no m O in Co t 'n 'n ^o 0^ ON^^ 00 * S0o n 00 00 O ^ o o 00 a U oQo o o0000 5 4, _O I 00 o r0 N 00 0600V 'n enO\ 00Oo ^ ^ O O o o d '< m T o o o o o < S CtogOOief0 Coi en o 14 O0  2 C 00 *a E I ,z rt14 o In Iin In 000 (NI ON 00 0 00 00 cr) en en 0 0o en mi en \o o o in tf 0g 0 0 0 0 O0 N ' 0 0o 0 o 0 0000 ^ O c g so O i; O O O; i ooo i< 0 U C.) F o 0 '.. 2.# l 22 oo oo ??ccc^^c ^ ,. '^ ,( ,^ ~ ^^^ ^ ^^.c S^r 19 5 X *S I 1 m I' i i 11) 1 I' '' 1 11"i' l I' o :0 0 S*** D m O N m a I I F'' C C *" n0 '0 I Cr, 0 a N 0 N 0 0 C Cf g "\ g \ g C \ 8 B Sc oE / /  *M *e 5.. SE 1.7.. I ( / ^ 3 20 S 0 U 8 O O oo _ 13~~ c'F o ^ 0^ 0 0)g  Co aggg I< I< ^^ ~ I s^ "* 1f '* a .^. P2 I C P2m aCa 0o^^" 0 * 0 N ? " N N^_ V1 <^ ." N ^CT '~ .S 1 I ^ U a a C S oN ~ 0 o c (. 0/ 0.) N 0 0 3cru aE *U cu 'f >E S. (M~i) vS C aC O I (U)uis r 0 o ea *. *..* * (m)uis aa C) 9l " o 00 *r 0 Qg 'eg bC) a .g U 0  c.2 _ paoJJo d * 0 wr(o i S O p~Jaojq CHAPTER ORBITAL EVOLUTION OF ASTEROIDAL PARTICLES Radiation Pressure, PoyntingRobertson Drag, and the Corpuscular Drag The orbital evolution of micronsized dust particles in the solar system is affected by several forces in addition to the gravitational forces due to the Sun and planets. These include radiation pressure, PoyntingRobertson (PR) drag, and solar wind corpuscular drag. The radiation pressure force is the force exerted on a dust particle due to the interception of sunlight. It varies as the inverse square of the distance from the Sun and depends on the properties of the particle (size, chemical composition, etc.). When a micronsized dust particle is moving around the Sun, it absorbs the sunlight and reradiates the energy away isotropically in its own reference frame. However, when viewed from the rest frame this particle radiates more energy along the same direction as its velocity vector. This will cause a net loss in its orbital energy and angular momentum. It acts like a retarding force on the particle. This deceleration drag causes the particle to spiral towards the Sun. It is called PoyntingRobertson drag (Poynting, 1903 , Robertson, 1937). This effect is present even if the particle scatters the sunlight instead of absorbing and reradiating the .. . 26 The solar wind corpuscular drag is due to the interaction between the solar wind protons and the dust particles. drag. The physical interpretation is similar to that of the PR The classical paper dealing with radiation pressure, PR drag, and corpuscular drag is the paper by Burns et al. (1979). The equation of motion of a dust particle with geometrical cross section A and mass m moving under the influence of the gravitational forces of the Sun, planets, and the radiation force from the solar radiation can be written as mv = rGMm GMnm (3.1) A/c)Qpr (1 r\  C c/ where v and r are the velocity and position vectors of the particle in the heliocentric coordinate system, while Mn, rn are the mass of the nh planet and the position vector of the particle with respect to that planet, respectively. S, c, and Qpr are the solar energy flux density, the speed of light, and the radiation pressure coefficient, respectively. The particle's radial velocity and the unit vector of the incident radiation are represented by t and S, respectively. The nonvelocity dependent part of the third term at the right hand side of Equation 3.1 is the radiation pressure force term, while velocity dependent parts are the drag terms. Qpr depends on the properties (density, shape, size, etc.) of the particle. It can be calculated from the Mie theory (Burns et al., 1979, Gustafson, 1994). Traditionally, a dimensionless quantity, P/, is introduced to specify the effect of radiation pressure and PR drag. It is defined as radiation "1 1I force r (3.2) 27 The major topic in this chapter is to study the dynamical evolution of asteroidal dust particles which contribute most to the 25 pm waveband IRAS observations. we adopt the solar wind drag as being 35% of the PR drag, the major contributors to pm waveband are particles with y3 equal to 0.05037, corresponding to 9 pm in diameter. These results are obtained by assuming 2.7 gm/cm3 perfectly spherical astronomical silicate and using theory to integrate over the solar spectrum (Gustafson, 1994). The corpuscular drag is 35% that of the PR drag according to the latest average solar wind data from Helios 1 and 2 (Schwenn, 1990). When a dust particle is released from a large parent body, it immediately "sees" a less massive Sun, due to the radiation pressure. Because the dust particle has identical position and velocity vectors as its parent body (assuming it leaves with zero relative velocity), this reduction in the centripetal force will change some of its orbital elements right away. Since radiation force is a radial force, it has no effect on I and fl of the dust particle; only the semimajor axis, eccentricity, and longitude of pericenter need to be recalculated. The new semimajor axis, an, and eccentricity, en, are given by an=a  2a'3/r (3.3) (1 2a/3/r)(1 p3)2 (3.4) X where a and e are the semimajor axis and eccentricity of the parent body and r is the radial distance from the sun at which the dust particle is released. From the first equation we can see that the new semimamir axis of the dust article is always larger  e2 en than that of its parent body. However, this is not true for the new eccentricity. new semimajor axis and eccentricity of the dust particle also depend on where it is released (i.e., r). It is not too difficult to show that both an and en have their maxima at perihelion release and have their minima at aphelion release. When the quantity X in Equation 3.3 becomes negative, the new eccentricity is the square root of the absolute value of X, and there is an 1800 shift in the longitude of pericenter of the dust particle compared with that of its parent body. Static Theory and Dynamical Theory The classical secular perturbation theory mentioned in Chapter 2 deals particles having constant semimajor axes. This is the socalled static theory. Micron sized dust particles in the solar system will spiral in towards the Sun due to the drag force in about 50,000 years (Wyatt and Whipple, 1950, Burns et al., 1979). Thus, it is necessary to develop a new secular perturbation theory which includes the effect of drag. We use the theory developed in Dermott et al. 1992. We test the theory on asteroidal dust particles in this section. In the loworder static theory, the rates of change of p, q, h, and k are N E vj cos(gjt + yj) j=l (3.5) Ah + Aq + I = An  vy sin(gyt + yj) PJ cos(fyt + yj) (3.6) "' u." fsin( f;t 4+ 'y h = Ak where A, v gj, fj, Pj, and 7j, are defined in Chapter When there is no drag, the precession rate, A, is spirals in towards the Sun. a constant. When there is a drag force, a dust particle Its semimajor axis becomes a function of time. So, A, vY, and pj (through their dependence on semimajor axis) all become functions of time. The eigenfrequencies and phases, gj, fj, and 7j, are constants regardless of the presence of drag force. In order to solve the problem, we define a complex equation to combine the equations involving p and q (or h and k), namely z=q + zp . (3.7) By using z, Equation 3.6 becomes N z = iAz+i ei(ft+ . j=l In the static theory, the solution of the above equation is (3.8) pei(At+Y) + Sj ei(tft+1+) = A+ f j~l (3.9) where j and 7 are the proper inclination and node, respectively, of the particle at time equals to zero. When A and /j become functions of time, Equation 3.8 becomes a standard first order inhomogeneous differential equation. Its solution is z(t) = e e'rdt (3.10) where iAdt (3.11) 13jei(fut+rY) (3.12) The constant C in Equation 3.10 depends on the initial conditions of the dust particle. Before we find the solution of Equation 3.10, we need to know the rates of change of semimajor axis and eccentricity due to PR drag. They are (Wyatt and Whipple, 1950, Burns et al., 1979) (yr/a)Qpr( + 3e2) (3.13) 1 e2)3/2 5'Qpre (3.14) 2a2(1 e_ 2)1/2 where 17 = 2.53 x l01 /ps for spherical particles in heliocentric orbits with density p, radius s in (c.g.s.) units, and Qpr is the radiation pressure coefficient. Now let us take a look at how the drag force affects the forced and proper eccentricities. We first consider a group of dust particles from a given family having a distribution in the (h,k) phase space as shown in Figure 3.1. circle is The equation of the given 2 = epO (3.15) where he and kc define the forced elements (i.e., the center of the circle) at time equal to 0 and epo is the initial proper eccentricity. The total eccentricity of a given particle r = i da dt) de dtJ + (h hc)2 center of the circle moves to ni forced eccentricity changes from efo to en. while the radius of the circle changes to ep1. The change in the proper eccentricity the radius of the circle) is due to drag, while the change in the forced eccentricity and longitude of pericenter is due to both gravitational perturbation and drag. From Equation 3.14 we particle see that the drag force changes the total eccentricity of any given by the amount 6e= Cleto (3.16) where CI is a positive constant. in the circle, etl The new total eccentricity of the previous particle then becomes = (1 C1)eto The equation of the new circle after S [(1 C1)k (1 1 C1) kcj (3.17) becomes 2 + [(1 C1)h (1 C)h]2 (3.18) = [(1 Ci)epo] This means that the drag force moves the center of the circle (which is defined by the forced eccentricity and longitude of pericenter) towards the origin by the factor of (1C1) and, at the same time, makes the circle (defined by its radius, the proper eccentricity) shrink by the same factor. This factor, can of course be calculated analytically from Equation 3.14. The overall picture of the dust evolution thus becomes very clear. The gravita tional perturbation changes the forced eccentricities of the dust particles while having Sl .* m i I i . + h l 32 towards the Sun and changes their proper and forced eccentricities the same way it changes their total eccentricities. By combining Equations 3.13 and 3.14 we can obtain the following equation: 5e(1 e2) (3.19) 2a(2 + 3e2) For asteroidal dust particles we can drop the e2 terms. log(ep) plot during their evolution is then 1.25. This v The slope in a log(a) vs. vill be verified numerically the next section. Numerical Comparisons In this section we first test the dynamical theory by placing dust particles at 1.8 AU initially in a (SunJupiterSaturn) system. This is to avoid the passage through mean motion resonances associated with Jupiter. Then we place dust particles initially at 3 AU and observe how they are affected when they do go through the Kirkwood Gaps (to be defined later). Drag Effect We place 254 Eostype dust particles initially at 1.8 AU to start the numerical integration. Again, we use the "particles in a circle" set up their initial orbital elements. and inclinations as those of Eos. method (Dermott et al., 1992) to The particles have the same proper eccentricities Their initial forced elements are calculated using the eigenfrequencies, amplitudes, and phases from the numerical integration of the (a, e, and w) due to the releasing process. The reason for this is to simplify the problem and study the effect of drag on the orbital evolution alone. The dust particles are integrated for 20,000 years. We plot their (p,q) and (h,k) once every 2,200 year Figure 3.3. They end up around 0.5 AU. s as shown in Figure 3.2 and The fact that they remain in nice circles in both (p,q) and (h,k) phase spaces during their evolution proves that the method of using forced and proper elements to describe their evolution in the presence of drag force is valid. Figure 3.4 shows that the proper inclinations of the dust particles are unchanged during their evolution from 1.8 AU to 0.5 In Figures 3.5 through 3.8 we compare the results from the numerical integration (shown as dots) and the prediction from the dynamical theory (solid lines). The forced inclination and node will be shown in Chapter 5 to be very important parameters when we construct the zodiacal cloud model. difference After 20,000 years of evolution, the the forced inclination between the dynamical theory and the actual numerical integration is about 0.02 degree. less than 1 degree. The difference in the forced node is This shows that the dynamical theory can describe properly the evolution of dust particles under the influence of both gravitational perturbations and PR drag as long as the effects of mean motion resonance and point mass scattering are not important. In Figure 3.9 we plot the forced and proper eccentricity in a log(a) vs. frame. log(e) Two reference lines with slopes equal to 1.25 are also shown in the diagram. A uxe^ mentirnndP intbn th lcl etrtln rQn frrst ;0 ts etnwr foutr tlhat hhanao c 34 the proper eccentricity. It changes the proper eccentricity according to Equation 3.19, which predicts an 1.25 slope in the log(a) vs. log(e) plot. The forced eccentricity does not approach the reference line until the particles are very close to the Sun. This is because both the gravitational perturbation and drag act on the forced eccentricity. When particles are far away from the Sun, the gravitational perturbation dominates the change in the forced eccentricity. It is only when the particles are very far away from the planets (Jupiter and Saturn), that the drag becomes the dominant factor and the forced eccentricity in the diagram approaches the straight line. Passage Through Resonance In the real solar system, asteroid dust particles produced from the main asteroid belt will go through several Kirkwood Gaps and secular resonance regions while they are spiralling towards the Sun. Kirkwood Gaps are the regions at which a particle is in mean motion resonance with Jupiter (e.g. Dermott and Murray, 1983). For example, resonance is at a=2.46 resonance at AU for particles having /=0.05037 . The locations of the Gaps are slightly shifted due to the radiation force. Secular resonance regions are the places at which a particle's orbital precession rate matches one or a simple combination of several of the eigenfrequencies of the solar system. The major secular resonance region is located near Williams, J. 1969, Scholl et al., 1989). Passage through different Kirkwood Gaps has some complicated effects on the orbits of the 35 is beyond the scope of this dissertation. However, one can observe their effects from direct numerical integration, as we have done here. Again, this is a SunJupiterSaturn system. The initial forced elements calculated similar to the methods mentioned in the previous section. and fifty four Eos dust particles are placed at 3.015 AU at t=0. The Two hundred ir other orbital elements are generated by the "particles in a circle" method (Dermott et al., 1985, 1992). The changes in orbital elements due to the releasing process are not included. These dust particles are integrated for 50,000 years. They end up at about 0.7 The distributions in the (p,q) and (h,k) phase spaces of these particles while they are moving towards the Sun are shown in Figure 3.10 and Figure 3.11. time interval between two consecutive panels is 6,000 years. The first impression from the figures is that they remain in a welldefined circle in the inclination space while going through large scattering in the eccentricity space. In Figure 3.12 we plot the proper inclination vs. semimajor axis of the particles. The proper inclination jumps up and down somewhat (particularly at the 3:1 resonance) during the particles' passage through Kirkwood Gaps even though the magnitude of change is not very significant (from the observational point of view). Figure and Figure show the variation in the forced inclination and longitude of ascending node from the numerical integration (dots) and the results from the dynamical theory (solid line). The difference is quite small. This implies the dynamical theory works well in inclination even in the case of mean motion resonance passage. The variations z C .... r .. ..! L ji 1 2 C A  r " _ F,  2 1^ ___ 3.16. The sudden drop at the 3:1 resonance region is very obvious. The difference between the numerical results and the dynamical theory is around 30% near the end of the integration. In all the diagrams there is no dramatic change that can be seen when particles pass through the secular resonance region. We have also performed another test run here. We changed the masses of Jupiter and Saturn to 10% of their real masses and repeated the same numerical integration. This was to reduce the strength of each mean motion resonance (Dermott and Murray, 1983). The difference in the forced eccentricity was dramatically reduced this time as we can see in Figure . This also supports what we have just suggested, namely that the discrepancies in Figures 3.15 and 3.16 are actually due to passage through mean motion resonance zones. From the above results we can conclude that passage through Kirkwood Gaps does produce a significant effect on particles' eccentricities and longitudes of peri center. Orbital Evolution of Dust Particles from the Major Hirayama Asteroidal Families The orbital evolution of dust particles in the real solar system is more complicated than the cases presented previously. First of all, even if the parent bodies which produce dust particles constitute a welldefined group (i.e. they form nice circles in both (p,q) and (h,k) spaces), the releasing process changes some of the orbital 37 the dust particles from a welldefined family do not form a welldefined circle in the eccentricity space from the very beginning. This makes the attempt to study analytically the variation in forced and proper eccentricities of dust particles coming from the same family very difficult. Secondly, inside the orbit of Jupiter, there are three important planets, Mars, Earth, and Venus. When dust particles pass right by these planets, there are singularities (due to zero distance to the planet) which cannot be handled analytically. The third reason is that some of the dust particles will be trapped in the nearEarth mean motion resonance (and to a lesser extent, nearMars and nearVenus resonances). Due to these difficulties, the practical way to study the orbital evolution of micronsized dust particles is direct numerical integration. In this section we perform the numerical integration for dust particles coming from three major Hirayama asteroid families, Eos, Themis, and Koronis. All the planetary perturbations (Mercury and Pluto are not included), radiation pressure, PoyntingRobertson drag, and corpuscular drag are included in the calculation. changes in the orbital elements due to the process of release are also included. Particles are assumed to leave their parent bodies with zero relative velocity. In order to utilize the vector facility on the IBM ES/9000, 7 planets for each integration. we set up 249 dust particles and The solar system is wound back in time numerically to an epoch such that dust particles being released at that epoch will reach a chosen heliocentric distance 1983. initial forced elements of particles calculated from the orbital elements of the 7 nrnnsrar nrnAsC anrl nrinti'h1A f narlrn^a tr planets at the beginning epoch. Anrra n r n n n, A ,,1.i, /^iic'i nnr between 00 and 3600 Dermott et al. The proper eccentricities and inclinations of families are from (1985). 9 um Diameter Particles We first integrate one wave of particles from each family. Each wave ends up inside 1 AU in 1983. For Eos particles, the variations in (p,q) and (h,k) spaces during their evolution are shown in Figures 3.18 and 3.19. The first panel in each figure is the distribution of the parent bodies which produce the dust particles. The time interval between two consecutive panels is 6,500 years. The label above each panel is the average semimajor axis of the dust particles at that time. In the inclination space, Eos dust particles remain in a circle even after they pass the Earth while Themis and Koronis dust particles have large scatterings when they approach the Earth. This is because Themis and Koronis dust particles have much smaller inclinations compared to that of Eos dust particles. The probability of low inclination objects being scattered by terrestrial planets is higher than that of high inclination objects. In the eccentricity space, PoyntingRobertson drag reduces the eccentricities of dust particles while they are moving towards the Sun. However, when they approach the Earth some of them are trapped in the nearEarth mean motion resonances. These particles gain energy and lose angular momentum during the process. Eventually they will escape from the resonance zone and continue their journey towards the Sun. This phenomenon will be discussed in more detail in the 39 Obtaining the particle distribution at different heliocentric distances in 1983 (for the purpose of constructing the zodiacal cloud model) is different from studying the dynamical evolution of the particles. Different waves (each "wave" is the numerical integration of the orbits of 249 dust particles) of particles must be sent out at different epochs such that they end up at different heliocentric distances in 1983. This task requires considerable CPU time. data for three families took about 3 months on an IBM ES/9000 supercomputer. For Eos dust particles, we sent in 23 different waves dating back to the year 57,500 years before 1983. between each wave is 2,500 years. integration and obtained the orbital distances. The time interval Then we analyzed the data at the end of each elements distribution at different heliocentric Thus, we have the forced elements as a function of heliocentric distance from around 0.5 AU to 3.2 AU The forced inclinations and nodes for Eos particles are shown in Figures 3.20 and 3.21. Themis and Koronis particles have similar forced inclinations and nodes until they reach around 1 AU where they no longer remain in a circle due the gravitational scattering by the Earth. In the (h,k) space, particles are scattered even from the very beginning as we can see in Figure 3.19. Hence, we did not attempt to find the corresponding forced elements by fitting a circle to the points in eccentricity space. Instead, we used a tabular form from the data to generate the eccentricities and longitudes of pericenter for the dust particles when we constructed the zodiacal cloud model. Once the orbital element distributions of dust particles from a known family have l ^ ^ ^ ^ J ^  ._ _ _ _ l.. i* 1. 1 i i cloud from the results. We did just that for Eos dust particles. Figure 3.22 is the vertical crosssection picture of the Eos cloud. This is the cloud which may responsible for the dust band. Actually the socalled "solar system dust band" is not a band. Because of the drag effects, dust particles go all the way from their origin to the Sun (the reason there is a hole in the figure inside 0.5 AU is that we terminate the integration when dust particles reach there). The actual "band" is the enhancement of the observed flux at certain ecliptic latitudes due to the concentrations of dust particles at the edges. This concentration can be seen in Figure 3.22. In the real solar system, this cloud is embedded in the zodiacal cloud which includes dust particles from many different sources. However, this concentration is strong enough such that when the whole zodiacal cloud was scanned by IRAS, this "band" (more precisely, a tiny peak on top of the background) at around 100 was observed. The forced inclination and longitude of the ascending node of the dust particles determine the plane of symmetry of the zodiacal cloud. This plane is a warped plane as the forced inclination and longitude of the ascending node vary with heliocentric distances. When we observe the zodiacal cloud, the latitudes of the peak fluxes around the are determined warped plane. These positions are very welldefined in the IRAS observations. One of the most important achievements of the present work is being able to use the forced inclination and longitude of the ascending node from the dynamical study and the SIMUL model to predict the latitudes of the peak fluxes as observed by IRAS. Comparison of those predictions 41 In Figure 3.20 the data show that the forced inclination increases towards small heliocentric distances. This increase had been reported from observations prior to the launch of IRAS (e.g. Misconi and Weinberg, 1978, Leinert et al., 1980). Is this increase due to terrestrial planets, especially Venus (which has a very high inclination, 3.40) as some early papers suggested (Misconi and Weinberg, 1978, Gustafson and Misconi, 1986)? We examine this by taking out the terrestrial planets one at a time and integrating the orbits of Eos dust particles in each case. inclinations from these cases are in Figure 3.23. can be seen clearly. The resulting forced The effect of the terrestrial planets The increase in the forced inclination has nothing to do with the existence of the terrestrial planets. It is the passage through the secular resonance zone that increases the forced inclination. 4 pm and 14 pm Diameter Particles We assume that 9 pm diameter dust particles are the major contributor to the 25 Cpm waveband observation of IRAS. So far we have studied the dynamical evolution of 9 pm diameter dust particles and their distribution in the solar system. In Chapter 5 we will show a zodiacal cloud model based on particles of this one size. However, particles of different sizes do contribute to the 25 waveband observations as well. Eventually the model zodiacal cloud must include results from particles with a range of sizes. In this section we present the dynamical study for particles with diameters equal to 4 pm and 14 pm. This covers a range of 3.5 in sizes (a factor nf 43 in m ~ p. if the.v have the same den itv. The integration procedures here are similar to those in the previous section. corresponding CPU time for 14 pm diameter dust particles of course is much longer than that of 9 pm diameter dust particles. The forced inclinations and nodes as functions of semimajor axes for the three different size particles are shown in Figures 3.24 and 3.25. Among them, the 14 pm diameter dust particles have the smallest drag rate. They pick up the most increase in the forced inclination while they pass through secular resonance regions and end up with the largest forced inclination in the inner solar system. The difference in inclination at 1 AU between these three different size particles is about 1.5 degree. Obviously this difference must be considered when a better zodiacal cloud model (comprising particles of 4, 9, and 14 pm in diameters) is to be constructed. How to combine these results is very important in constructing a successful zodiacal cloud model in the future. NearEarth Dust Grains In this section we will describe some properties of asteroidal dust particles when they are in the vicinity of the Earth. These include seasonal variations in the number of particles that intersect the Earth, how to identify the origins of nearEarth dust particles, and the phenomenon of trapping in nearEarth mean motion resonances. The forced inclination and longitude of ascending node determine the plane of symmetry of the zodiacal cloud. They also produce another interesting phenomenon collected in nearEarth orbit. The distributions in inclination space for 4 pm and 9 pm dust particles are shown in Figure 3.26. These are particles in the Earth crossing region in 1983 found from our numerical integration. Notice the offcenter distribution for Themis and Koronis particles. There are more Themis and Koronis particles in the first and fourth quadrants than in the second and third quadrants of the diagram. When the Earth moves around its orbit over one year, it will encounter more Themis and Koronis dust particles at certain places. These places correspond to the areas in the inclination space where there are more particles (first and fourth quadrants). Figure 3.27 This kind of seasonal variation for 4 and 9 pm particles is shown in Here we have a histogram showing the variation in the relative number of particles as a function of the time of the year. The Earth encounters more Themis particles at their ascending nodes in November than in April. The difference is more than a factor of 8. There is no significant seasonal variation for particles with large proper inclina tions, like Eos. In the vicinity of the Earth the proper inclination of Eos particles is much larger than the forced inclination. Consequently, even though the circle is shifted from the origin in the inclination space (see Figure 3.26), the difference in terms of number of particles in each quadrant is not very large. Hence the seasonal variation for Eos particles in Figure 3.27 is not as dramatic as that of Themis particles. Themis, Koronis, and Eos are the most prominent families in the asteroid belt. The ratio of dust particles produced from these three families to nonfamily asteroids 'ten~~ 4.l a n 4..l .~ r A Ct r. e 4 t fl a * n 1t I .. Sc *n~ nA *Lr Themis : Eos : Koronis = 0.61 : 0.12 (Durda, 1993). This means Themis and Eos produce 17% and 10% of all the asteroidal dust particles, respectively. This, together with the facts that (1) particles from certain families do not arrive uniformly in time (e.g. most Themis particles with ascending nodes encounter the Earth around November) and (2) particles from one family have a welldefined range of their inclinations, makes it possible to identify the origins of some family dust particles. In Figure 3.28 we show the relative numbers of particles as a function of inclination. These are 4 and 9 pm dust particles at their ascending nodes which will encounter the Earth in November. Themis and Eos particles are represented by bold lines in the histogram while the whole nonfamily asteroidal particles are shown as thin lines. Much useful information can be obtained from this figure. First of all, the low inclination dust particles in November are dominated by Half of the particles around 120 inclination are of Eos origin. Themis members. If we utilize these results to interpret the samples from an Earthorbiting dustcollection probe, we can obtain detailed information (e.g. chemical composition) about some asteroid families without having to go to the asteroid belt. For example, from this dynamical study we know that half of the dust particles being collected in November with 12 inclination are Eos particles. Then, we can collect, say, 20 particles with inclinations around and analyze their composition. Half of them will have similar compositions (if we assume particles from the same family have the same compositions). must be Eos particles. family. These Their chemical compositions are those of the Eos asteroid Also, most of the low inclination (around 30) dust particles being collected in November must have similar compositions. These are particles from the Themis family. An important objective for any future dust collection mission is to establish a relationship between the dust particles being collected and their parent bodies. results here have shown that it may be possible to achieve such a goal. With modem techniques, intact capture of particles with velocities less than 9 km/sec is possible (these are most likely to be asteroidal dust particles). From our study, we can provide not only the best time but also the orientation of the orbits of the dust particles that can be collected and related to several known asteroidal families. When dust particles approach the Earth from outside, some of them will be trapped in the nearEarth mean motion resonances. on many factors. The trapping probability depends These include the drag rate, the orbital elements of the particles, the strength of the resonance, and the geometry of how a particle approaches the resonance. This trapping in mean motion resonance will produce some interesting observational consequences (e.g. Jackson and Zook, 1989 , Reach, 1991, Dermott et 1993). Here we will describe qualitatively the behavior of particles when they approach the Earth based on our numerical results. Similar phenomena happen near Venus and Mars, but to a lesser extent. The orbital evolution of particles being trapped in a mean motion resonance with a planet is, actually, very easy to understand. When a particle is spiralling towards the Sun, the drag force takes away its energy and angular momentum. When a particle 10 fronnaA in a oa nnai n no anno ono m Crrnn *tTfl annnnnnn fal nn,.ndaa *n i *n 4 loss due to drag force. Whether a particle will be trapped or not depends on how much energy it gains when it passes through a given resonance zone. If it does not obtain enough energy, it simply passes right through (with some minor changes in its orbital elements). Otherwise it will get trapped. However, this particle will not be trapped forever. Eventually it will have a series of close encounters or a very deep close encounter with the planet. particle from Then, the close encounter removes this the resonance. In our numerical integration, diameters equal to 4, 9, and 14 pm. than that for the big particles. Hard we have studied the evolution of particles The drag rate for small dust particles is larger ly any of the 4 pm dust particles are trapped in any resonances. A few of the 9 pm dust particles are trapped. show a wave of 9 pm Eos particles moving towards the Earth. show the semimajor axes vs. In Figure 3.29 we In Figure 3.30 we time for 50 of them when they approach the Earth. Some of them do get trapped, but the overall trapping in resonance is not very significant. 14 pm dust particles, the trapping percentage is higher because of their lower drag rate. ^4 ^.,.."a.g ..*  \ B 44 tI3 f t I Q.I 0 0 = o * C Ua ** e *wU #rar U* *... S S S S U .. * (1  \., * *.. ** egggggggst4 ^*. **... c *a pa (u)UTs *  /'" N I 1 iN (c'uis 4 1~~ 0 C02 jadosidI E O co C 1n Q paosoJI I I I I I paoJoJLj 'a  M =< C a0  rn al cI. C C0 00 CO a " V 'Sc *M " E .CS) S2g I ^ CO w T3C3 *5 &D C *oE a  Cl? 0 paojoj a r; cn a = 9 Co 0O 0 0 C 0 CD g .9 V1k ""^ /^ M O * O O L  OC pa~oj.oj (IoiTUa u oog) (O)Uis  * \ i t^I 44 p ^*w ' 4. i 4, .a 'i.. * ' 'I i I x~~r /1 _ f ' t/ ..* t Goc S0 c1  (m)uIs II I I I I*I 'ml I*I 'ml !*t I*! I~* !*! I  ' I  I  I   l e  I I I I I f II I I I Il nadosdI T I Il I 10 C1 paiJOJI I I I I I I II I I I IIIIII I 1 I 1 I ,, I ,, I ,,, I paaojO paojoj9 paDaojo0 o I. Cu3 6 c e ow paoJoj9 C ) . ET (u)mus ** rN \* (T 5)UIS 1I 0 rI II gJ. II I I I I I I I I I I I I I I I I I I I I I I I I I paojoJI I I ll I I I I Il 1 1 =^ c> "o *i^ 0E i.. paojojz Cs) * I ~r (nv) !aH o .d.p t paaojI ^s do, ^S ca Sa '0 14 paojojI xe p9aJOJJ V1 ^ ^i CM CL 0 0 0 4 I O 01 S o 0 . 2o\ O h4 *.!M . S .  C *  J** (U)uis O C\2 j98qujnf nCJ " ^...^ ^^ ^ GAnT.PIG "r 103 Cr3 "r CO T l 1 O1 s^ 0' ' 08 o o 2 CO 0 0 aoq unN ahiii3Ta'ti 5 E  8...' .4,*** * *  * * .*  *  * .. * I I ' ... U . a ' . * . *C W I I 8 *,  I 8 *  * II we *^. * .b * '3 L  V1 (nv) 'x x * ~ ug c;EE &bw. ^ fl *A ^ o Cumimu aS CHAPTER ORBITAL EVOLUTION OF COMETARY PARTICLES major difficulty studying the dynamical evolution cometary particles sources. and their contribution to the zodiacal cloud is that of identifying their Unlike asteroidal dust particles, which are mostly from the main asteroid belt, cometary dust particles have no welldefined parent bodies. We do not know whether they are from a highly active comet, a group of comets, or all short period comets. This makes the attempt to find the contribution of cometary particles to the zodiacal cloud seem impossible. Nevertheless, we feel that, even if we don't know the exact source of cometary dust particles, by studying the characteristics of the model zodiacal cloud produced from a typical cometary source, we will gain some insight into the overall problem. The typical cometary source we choose is comet Encke. In this chapter we study the orbital evolution of Encketype dust particles numerically and find their distribution in the solar system. We also show that comets cannot be the source of the solar system dust bands observed by IRAS. Orbital Evolution of Encketype Dust Particles Encke has long been proposed to be the major source of the zodiacal cloud *^ I 4 4 S S S C C S S C C C C xI 4 ^\ /^jrv 78 against the gravitational perturbation of Jupiter (Opik, 1963) and because of its strong correlation with several meter showers and the socalled Taurid Complex (e.g. Steel, et al, 1991a,b, Hartung, 1993). It is thus natural to study the orbital evolution of Encketype dust particles when there are no other better choices. Our definition for Encketype dust particles is those with initial semimajor axes, eccentricities and inclinations similar to those of Encke. The mean semimajor axis, eccentricity, and inclination for Encke are 2.2 AU 0.85, and 12, respectively. The initial longitudes of ascending nodes, longitudes of pericenters, and mean longitudes of dust particles are randomly chosen between 00 and 360 . This means these micronsized dust particles are not released by means of direct sublimation from the parent body when it is near its pericenter. The PR drag lifetime for Encketype micronsized dust particles is about 5,000 years. If they are produced from a comet at its pericenter, they will not have time to process and spread out their orbits to produce the axial symmetry structure of the zodiacal cloud. We assume that it is the disruption of large Encketype cometary meteoroids (large enough such that their longitudes of ascending nodes and longitudes of pericenters have been dispersed between 00 3600) that produces these dust particles. Griin et al. (1985a,b) concluded that this is the most efficient way to deposit cometary dust particles into the zodiacal cloud. We first set up the initial condition for five hundred 9 pm diameter dust particles by the "particles in a circle" method (Dermott et al., 1985, 1992). We start the integration at five thousand years before 1983. Dust particles are then integrated for 79 variations in the (h,k) and (p,q) spaces when the particles are spiralling towards the Sun are shown in Figures 4.1 and 4.2. panels is six hundred years. The first The time difference between two consecutive : panel is the distribution of the parent bodies at time equal to zero. The major difference between Figure 4.2 and Figure 3.18 is that cometary dust particles do not remain in a welldefined circle in inclination space right after they are released. If we examine the inclination of each individual particle, we will find that it goes through huge variations as the particle moves towards the Sun. As a whole, cometary particles form a doughnutlike distribution with welldefined maxima and minima in the inclination space. of cometary type dust particles. These features are the most important characteristics They may play crucial roles in combining with the asteroidaltype dust particles to form the zodiacal cloud. detail in Chapter 5. This will be discussed in The physical interpretation for these inclination variations is discussed in the next section. The evolution of Encketype dust particles in (a,e) space is shown in Figure These five groups correspond to data from one thousand (large filled circles) to five thousand (crosses) years after they were released. Very few particles are trapped in mean motion resonances with Jupiter or the Earth. This is due to their high eccentricities which give them much higher drag rates than asteroidal particles with the same pf and thus very much lower probabilities of capture. We demonstrate, particle. The drag rate of the dust particle is given by (see Equation 3.13) const.(2 + 3e2) const. a (4.1) where a and e are the semimajor axis and eccentricity of the dust particle, respec tively. With the same starting semimajor axis, the drag rate of the dust particle increases dramatically towards the high eccentricity range. For example, da/dt of an eccentricity 0.85 particle is nearly particle. 14 times larger than that of an eccentricity 0.1 Comparing Encketype dust particles with typical asteroidal dust particles, we notice that this high drag rate makes them move much faster through the mean motion resonance zone in such a way that they do not have time to gain enough energy from the resonance to counterbalance the loss due to PR drag orbital decay. Gauss Equations The reason why Encketype dust particles go through huge variations in their inclinations can be understood from Gauss' equations. These equations provide a better intuitive understanding of the gravitation perturbations from planets (see, for example, Danby, 1988) than the Lagrange equations. The rate of change of the inclination of the particle is given by  1Wcos(w + f) na2vTC where r is the heliocentric distance of the dust particle, (4.2) w and f are the argument pericenter and anomaly the dust particle, respectively. W is the e2 of the particle. When we consider the longterm average perturbation, it is obvious that as the eccentricity of a dust particle' W increase. orbit increases, both the factor 7 and This will lead to an increase not only in the rate of change of inclination but also in the maximum and minimum inclinations that particle can reach. Particles starting with the same orbital elements will reach the same maxima and minima in their inclinations. Can Comets be the Source of the Solar System Dust Bands? IRAS discovered several dust bands during its allsky survey in 1983. central band (later it was found that there are actually two central bands) and the 100 band are very clear in all four wavelengths. Several theories have been proposed to explain their origins. Dermott et al. (1984, 1985) first pointed out the possible correlation between these bands and three major Hirayama asteroidal families and later confirmed it from their dynamical study and modeling results (Dermott et al., 1990, 1992, 1993). The central bands are probably due to the Themis and Koronis families while the Eos family may be responsible for the 100 band. However, the possibility still exists that some comets may produce these bands (Dermott et al., 1984, Clube and Asher, 1990). Among all the possible candidates, Encke has been proposed to be the source of the 100 band. This is primarily due to the fact that Encke's mean inclination is very close to 10 degrees. Here we will discuss, from Before we show why there is no cometary origin of dust bands, we must first understand how to form a dust band. In the (p,q) phase space, the osculating inclination a dust particle is the vectorial sum two components, forced and proper components. proper inclinations of asteroidal dust particles remain unchanged during their orbital evolutions Figure 3.18, a key feature to form a dust band) while their proper nodes process about a point in inclination space defined by the forced inclination and node. From the geometric point of view, this means that the orbits of asteroidal dust particles process about the axis of a common plane, the plane of symmetry, which is defined by the forced inclination and node. Dust particles produced from the same family will form a doughnutlike structure if the family members are old enough such that their nodes have been spread all over the sky. This doughnutlike structure is tilted to the ecliptic with the inclination and node defined by the forced inclination and node. When we observe this structure from the Earth, we will see two maxima in flux at the two edges of the doughnut at any given ecliptic longitude. This is the socalled "dust band" . These maxima still exist even when we consider the fact that dust particles are spiralling towards the Sun and fill up the space between their origin and the Sun (see Figure 3.22). With this band forming mechanism in mind, we can now examine the orbital evolution of cometary dust particles and see if they are capable of forming a dust band. As we mentioned previously, while asteroidal dust particles remain in a well variations in their inclinations. From the geometric point of view, this means there will be no welldefined doughnut like structure, hence no welldefined edges. Thus, it is impossible for cometary particles to form a dust band. Another observational consequence from this huge variation in inclination is that cometary dust particles will contribute more (relatively speaking, when compared with asteroidal dust particles) to the observed flux at high ecliptic latitudes. Chapter 5, we will show that this may be the key feature needed to account for the observed shape of the zodiacal cloud. Distribution of Encketype Dust Particles in the Solar System In order to find the contribution of Encketype dust particles to the zodiacal cloud we must have an idea of the distribution of these particles in the solar system at the time of the IRAS observations. We need the distributions of all the orbital elements. We have constructed two models, to be described in this section. The only difference between these two models is in their initial inclinations. Tl used to obtain the dust distributions are similar to those in Chapter ie procedures We send out waves of particles once every 300 years starting from year 5,400 years before 1983. Most of the particles being released before the year (19835,400) will have semimajor axes less thai contributions to the IRAS flux data orbits are integrated numerically. AU by the year 1983 and thus make no Each wave contains 500 dust particles. All the planetary perturbations from Their 7 planets, the numerical integration. The 3 for the dust particles is assumed to be 0.05037 The corpuscular drag is assumed to be 35 % that of PR drag. The initial semimajor axes and proper eccentricities of particles are 2.2 AU and 0.85, respectively. In the first model, the initial proper inclinations of all the dust particles are set equal to Encke's mean inclination, 120. In the second model, the initial proper inclinations of the particles are generated according to Encke's past. variation in inclination in the The longitudes of ascending node, the longitudes of pericenter, and the mean longitudes of the dust particles are generated randomly between 0 and 3600 . The forced elements at 2.2 AU are calculated according to the orbital elements of the planets at the time of particle release. All the calculations are performed on an IBM ES/9000. We exclude particles when their semimajor axes become less than 0.1 AU or when their orbits become unbounded. For each integration, all six orbital elements (a, e, I, fl, w, A) of the dust particles are recorded every 100 years until the end of the integration. Finally we combine the results from the end of each integration to get the orbital elements of particles at different heliocentric distances in 1983. The final distributions of dust particles from the second model (nearly 8,000 of them) in (a,e), (a,I), and (f2,w) phase spaces are shown in Figures 4.5 to 4.7 Few particles are found to be trapped in mean motion resonances with Jupiter or the Earth, as we expected. Figure 4.6 shows the distribution in inclination as a function of semimajor axis in 1983. Particles are distributed quite uniformly between 20 and 350 Particles with 1*;U ~ ..11 ;*,^1:./ ..,, ,;1 LL, t!,t _l 1* ^ _~( ^ l ^ * ^ i  