ABSTRACT Title of Dissertation: RESONANT AND SECULAR ORBITAL INTERAC- TIONS Ke Zhang, Doctor of Philosophy, 2007 Dissertation directed by: Professor Douglas P. Hamilton Department of Astronomy In stable solar systems, planets remain in nearly elliptical orbits around their stars. Over longer timescales, however, their orbital shapes and sizes change due to mutual gravitational perturbations. Orbits of satellites around a planet vary for the same reason. Because of their interactions, the orbits of planets and satellites today are di erent from what they were earlier. In order to determine their original orbits, which are critical constraints on formation theories, it is crucial to understand how orbits evolve over the age of the Solar System. Depending on their timescale, we classify orbital interactions as either short-term (orbital resonances) or long-term (secular evolution). My work involves examples of both interaction types. Resonant history of the small Neptunian satellites In satellite systems, tidal migration brings satellite orbits in and out of resonances. During a resonance passage, satellite orbits change dramatically in a very short period of time. We investigate the resonant history of the six small Neptunian moons. In this unique system, the exotic orbit of the large captured Triton (with a circular, retrograde, and highly tilted orbit) in uences the resonances among the small satellites very strongly. We derive an analytical framework which can be applied to Neptune?s satellites and to similar systems. Our numerical simulations explain the current orbital tilts of the small satellites as well as constrain key physical parameters of both Neptune and its moons. Secular orbital interactions during eccentricity damping Long-term periodic changes of orbital shape and orientation occur when two or more planets orbit the same star. The variations of orbital elements are superpositions of the same number of fundamental modes as the number of planets in the system. We investigate how this e ect interacts with other perturbations imposed by external disturbances, such as the tides and relativistic e ects. Through analytical studies of a system consisting of two planets, we nd that an external perturbation exerted on one planet a ects the other indirectly. We formulate a general theory for how both orbits evolve in response to an arbitrary externally-imposed slow change in eccentricity. RESONANT AND SECULAR ORBITAL INTERACTIONS by Ke Zhang Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park in partial ful llment of the requirements for the degree of Doctor of Philosphy 2007 Advisory Committee: Professor Douglas P. Hamilton, Chair Professor Michael A?Hearn Professor Joseph A. Burns Professor M. Coleman Miller Professor James A. Yorke c Copyright by Ke Zhang, 2007 Napoleon: They tell me that you have written this huge book on the system of the universe without once mentioning its Creator. Laplace: I have no need for that hypothesis. Preface I had hoped that I could have arranged each chapter more carefully. I had hoped that I could write a better story-telling introduction about the history of celestial dynamics. I had hoped to do many other things in this dissertation until I was at the point of running out of time. But here it is - the work I have spent seven years on. This dissertation is organized in four parts: Part I is an introduction. The history of the development of celestial dynamics and orbital dynamics is covered in Chapter 1. In Chapter 2, I review the basics of perturbation theory, which is the foundation of my dissertation. Part II is on the resonant interaction and evolution of the small Neptunian satel- lites, in which I focus on the small inclinations of the blue planet?s six small satellites and try to build a resonant history of the system based on their current orbits. Chap- ter 3 provides background information about the moons, as well as theoretical prepa- ration necessary for this project. In Chapter 4, two new orbital elements are de ned for resonant analysis in this system. Individual resonance passages are deciphered in Chapter 5. Finally, several physical parameters of the system are constrained based on dynamical evidence in Chapter 6. Part III is on the secular evolution. A linear secular theory is derived to handle slow eccentricity-damping. In Chapter 7, I discuss the motivation of the project and present the standard secular theory. Eccentricity-damping is then added in Chapter 8, and secular theory is adjusted accordingly. Lastly, the theory is applied to extrasolar planetary systems in Chapter 9. Part IV is the conclusion, which also includes possible future research directions as extensions of the two projects. Part II and Part III are only loosely related, thus the order of reading is not important. The theoretical background in Part I, however, is helpful for the later chapters. Ke Zhang August 1, 2007 ii To my parents iii Acknowledgment First of all, I thank my adviser, Professor Douglas P. Hamilton. Without his direction in my research, I would never have done what I did, and without his advice on life in the United States, I may not have been able to adapt to this totally di erent world. I also thank the whole Astronomy Department here in Maryland for their support. I am deeply grateful to Professor Joseph A. Burns of Cornell University and Professor Richard Greenberg of the University of Arizona for valuable discussions. I highly appreciate the careful reading of this manuscript by my thesis committee members and their helpful suggestions. I especially thank my parents and my sister?s family, who have continuously sup- ported me all the time. I am thankful to all my friends, who have accompanied me during the seven years of study and research, shared my happiness and sorrow, and given me immediate emotional support during my di culties. Finally, I thank the University of Maryland and NASA for their nancial support. iv Table of Contents List of Tables vii List of Figures viii I Introduction 1 1 Celestial Mechanics and Orbital Dynamics 2 2 Perturbation Theory 12 2.1 De nition of Orbits . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2 Disturbing Function . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.3 Secular and Resonant Perturbations . . . . . . . . . . . . . . . . . . . 17 2.4 Rotational Deformation . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.5 Tides . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 II Orbital Resonances in the Neptunian System 31 3 Background 32 3.1 The Neptunian Satellites . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.2 Tidal Evolution and Mean-Motion Resonance Passage . . . . . . . . . 35 3.3 Resonant History of the Neptunian System . . . . . . . . . . . . . . . 38 3.4 Computing Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4 Perturbations from Neptune?s Oblateness and Triton 44 5 Proteus Resonances 52 5.1 The Recent 2:1 Proteus-Larissa Resonance Passage (PL 2:1) . . . . . 52 5.1.1 Eccentricity Evolution during and after the PL 2:1 Passage . . 53 5.1.2 Inclination Resonances in the PL 2:1 Resonant Zone . . . . . 58 5.1.2.1 Three-Body Resonances . . . . . . . . . . . . . . . . 61 5.1.2.2 Important Higher-order Resonances . . . . . . . . . . 64 5.2 The Second-order Resonance PD 3:1 . . . . . . . . . . . . . . . . . . 70 5.3 The PG 2:1 and Diverging Capture . . . . . . . . . . . . . . . . . . . 72 5.3.1 Resonant Trapping Condition . . . . . . . . . . . . . . . . . . 74 5.3.2 Evolution in RePi2T . . . . . . . . . . . . . . . . . . . . . . . . 77 5.3.3 Trapping into Ri2P . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.4 The Chaotic PL 3:2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.4.1 Eccentricity Resonances Overlapping Criterion . . . . . . . . . 82 v 6 Constraints on Physical Parameters 88 6.1 Satellite Densities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 6.1.1 Constraints from a Single PL 2:1 Passage . . . . . . . . . . . . 89 6.1.2 Constraints from All Proteus Resonances . . . . . . . . . . . . 92 6.2 Tidal Evolution Timescale and QN . . . . . . . . . . . . . . . . . . . 103 6.3 QP and QL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 III Secular Interactions with Eccentricity Damping 110 7 Background 111 7.1 Extrasolar Planets . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 7.2 Secular Modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 8 Secular Solutions with Eccentricity Damping 121 8.1 Secular Modes with Eccentricity Damping . . . . . . . . . . . . . . . 121 8.2 Apsidal Circulation and Libration . . . . . . . . . . . . . . . . . . . . 126 8.3 Relativistic Correction . . . . . . . . . . . . . . . . . . . . . . . . . . 135 9 Applications to Extrasolar Planetary Systems 138 9.1 A Test of the Theory: the HIP 14810 System . . . . . . . . . . . . . . 140 9.2 Constraints on Possible Companions of \Hot-Jupiters" . . . . . . . . 143 IV Conclusion 147 10 Conclusion and Future Directions 148 10.1 The Inner Neptunian System . . . . . . . . . . . . . . . . . . . . . . . 148 10.1.1 The 4:7 Inclination of Naiad . . . . . . . . . . . . . . . . . . 149 10.1.2 Con nement of Neptunian Ring Arcs . . . . . . . . . . . . . . 151 10.2 Perturbed Secular Interactions . . . . . . . . . . . . . . . . . . . . . . 152 10.2.1 E ects of Adiabatic Changes in a and m . . . . . . . . . . . . 152 10.2.2 Systems with Three or More Planets . . . . . . . . . . . . . . 153 References 156 vi List of Tables 2.1 Gravitational properties of major planets . . . . . . . . . . . . . . . . 24 2.2 Tidal timescales for natural satellites . . . . . . . . . . . . . . . . . . 29 3.1 Inner Neptunian satellites and Triton . . . . . . . . . . . . . . . . . . 34 5.1 Critical resonances for rst-order eccentricity resonance overlapping . 87 6.1 Inclination kicks through di erent resonance passages . . . . . . . . . 92 6.2 Q of giant planets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 9.1 Properties of planets and stars discussed in this chapter . . . . . . . . 139 vii List of Figures 1.1 The concept of epicycles . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.1 De nition of orbital elements . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Two-satellite system . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3 Mean-motion resonance stability . . . . . . . . . . . . . . . . . . . . . 21 2.4 Tidal torque . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.1 Images of the inner Neptunian satellites . . . . . . . . . . . . . . . . 34 3.2 Distances of the small satellites from Neptune . . . . . . . . . . . . . 36 3.3 Resonant history of the inner Neptunian satellites . . . . . . . . . . . 40 4.1 Laplace plane of the Neptune-Triton system . . . . . . . . . . . . . . 45 4.2 De nition of new orbital elements . . . . . . . . . . . . . . . . . . . . 48 4.3 Comparison of traditional and Laplacian elements . . . . . . . . . . . 50 5.1 PL 2:1 resonance passage without Triton . . . . . . . . . . . . . . . . 54 5.2 PL 2:1 resonance passage with Triton . . . . . . . . . . . . . . . . . . 55 5.3 Resonant angle of Ri2L during the PL 2:1 . . . . . . . . . . . . . . . . 59 5.4 Resonant arguments of three-body resonances . . . . . . . . . . . . . 62 5.5 The R resonances . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.6 Chaotic R zone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.7 The PD 3:1 resonances . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.8 The PG 2:1 resonant zone . . . . . . . . . . . . . . . . . . . . . . . . 73 5.9 The PL 3:2 passage . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 6.1 Inclination kicks versus satellite density . . . . . . . . . . . . . . . . . 90 6.2 PL 2:1 resonant passage with heavier planets . . . . . . . . . . . . . . 91 6.3 Three-body resonant kicks versus satellite?s initial free inclinations . . 93 6.4 Cumulative free inclinations of satellites . . . . . . . . . . . . . . . . 97 6.5 Possible initial con gurations of the system . . . . . . . . . . . . . . . 104 7.1 Eccentricities of \hot-Jupiters" . . . . . . . . . . . . . . . . . . . . . 112 7.2 Secular modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 7.3 Eccentricity ratios in secular modes . . . . . . . . . . . . . . . . . . . 118 7.4 Secular solution in a phase plot . . . . . . . . . . . . . . . . . . . . . 119 7.5 Secular changes of orbital eccentricities . . . . . . . . . . . . . . . . . 120 8.1 Orbital con gurations in damping secular modes . . . . . . . . . . . . 123 8.2 Eccentricity damping in secular modes . . . . . . . . . . . . . . . . . 124 8.3 Secular evolution for systems with m21a1 m22a2 . . . . . . . . . . . . 128 8.5 Apsidal state evolution for m21a1 m22a2 case . . . . . . . . . . . . . . 134 9.1 The current apsidal state of the HIP 14810 system . . . . . . . . . . . 141 9.2 Orbital states of possible companions for HIP 14810b . . . . . . . . . 143 9.3 Eccentricities of possible companions for \hot-Neptune" GJ 436b and GJ 674b . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 10.1 Resonance capture of Naiad and Galatea . . . . . . . . . . . . . . . . 150 10.2 Apsidal alignment states in a three-planet system . . . . . . . . . . . 154 ix Part I Introduction 1 Chapter 1 Celestial Mechanics and Orbital Dynamics Celestial mechanics has its origins in the cu-Moon Earth Epicenter Figure 1.1: The concept of epicy- cles: the Moon orbits in a cir- cle around an epicenter, which it- self moves in a circle around the Earth. For complicated orbits, several epicycles may be required and the epicenters may shift. riosity that human beings display towards the mysterious motion of objects in the sky, as well as a practical need for accurately recording the passage of time and predicting the seasons. It is among the oldest elds of modern physics and astronomy, and has been substantially developed even since Sir Issac Newton (1643 - 1727) pub- lished his famous Principia in 1687. Before that, ancient theories about the Sun, Earth, Moon, and planets trace back to the Greek astronomers Aristarchus (310 - 230 BC), Hipparchus (190 - 120 BC), and Ptolemy (90 - 168 AD) some two millennia ago. Aristarchus was the rst to propose a heliocentric model based on his estimation of a much heavier Sun than Earth. His view of the universe was opposed by most of his contemporaries and astronomers after him be- cause the lack of parallax of the Sun, and the absence of any perceived motion of the Earth. Hipparchus, with data from hundreds of years of Babylonian observations, measured the precession of the vernal equinox at a value of 4600 per year, close to Most material about the early history prior to Newton is based on the 1968 version of Ency- clopaedia Britannica. 2 the modern value of 50:2600 per year. He also made an accurate measurement of the length of a year to within 7 minutes, and the distance of the Moon to within 10% error. He tried to create a model of the Moon?s motion with epicycles (Fig. 1.1), but discrepancies with observations existed until the model was re ned by Ptolemy. In his book, Almagest ( 150 AD), Ptolemy detailed the mathematical theory about the motions of the Sun, Moon, and planets around the Earth. Although Ptolemy attributed much of his book to Hipparchus, including the original idea of epicycles, he was the rst to work out the big picture of a geocentric model of the universe. As more observational data accumulated with better accuracy, more epicycles were required for the geocentric model to match observation. By Copernicus? time some 1400 years later (1473 - 1543), each planet could have as many as 40 60 epicycles in order to match observations, which drove Copernicus to reconsider the heliocentric model of the universe (De Revolutionibus Orbium Coelestium, 1543). The heliocentric model did not gain popularity until 50 years after Copernicus? death, not only because it challenged the authority of the Church, but also because it used almost as many epicycles as Ptolemy?s model did in order to precisely agree with observations. In 1609, Johannes Kepler published his famous laws on planetary orbits, which claimed that planets, including the Earth, orbit in ellipses instead of circles, with the Sun at one focus. Although requiring a single extra parameter for each orbit (eccentricity), Kepler?s theory agreed with observations perfectly without the need for epicycles. These ancient works, although accurate enough to predict the location of the known celestial bodies, utilize mathematics no more complicated than simple alge- 3 bra and geometry. The underlying physics responsible for celestial motions was not understood even after Kepler published his laws, which were simply empirical rules based on extremely detailed observations by Tycho Brahe and Kepler himself. Never- theless, it was these empirical works that enabled Newton to understand the universal law of gravity and to found modern physics. Newton?s contribution was a huge tri- umph for astronomy, physics and mathematics. His universal law of gravity is still used today to guide spacecraft ying to the outer Solar System and to model the motion of the planets to exquisite accuracy. The orbital precession of Mercury, which requires a small correction from general relativity, is the only serious shortcoming of the theory. His invention of calculus (or co-invention with Leibniz) opened the door of mathematical analysis and made modern science possible. After Newton, theory on celestial mechanics was rapidly developed and reached its zenith with the works of two 18th century mathematicians: Joseph Louis Lagrange (1736 - 1813) and Pierre-Simon Laplace (1749 - 1827). Besides his foundational work in classical mechanics and his creation of the variational calculus, Lagrange found the \Lagrange" points (potential maxima or saddles, see Murray and Dermott, 1999) while attempting to solve the three-body problem, worked out a method to determine a comet?s orbit with only three observations, and did additional important work on orbital precession and stability. Laplace, through a series of memoirs to the Academy of Science in Paris, addressed the stability of the Solar System by showing that the changes of the orbital mean motions of Jupiter and Saturn were periodic and due to their near-resonance orbits (sometimes referred as the great inequality). He also 4 spent a signi cant amount of time in the study of lunar motion perturbed by a non- spherical Earth, and of the oceanic tides induced by the Sun and the Moon. His most signi cant contribution, however, was the compilation of the ve volume Celestial Dynamics (1799 - 1825), which \o er a complete solution of the great mechanical problem presented by the Solar System, and bring theory to coincide so closely with observation that empirical equations should no longer nd a place in astronomical tables." In these books, Laplace included most of his work on planetary orbits and perturbations, as well as problems solved by earlier astronomers. Research on celes- tial dynamics achieved a real predictive triumph when the British astronomer Adams (1846) and the French astronomer Le Verrier (1846) independently \discovered" Nep- tune by analytical calculation of its perturbation on the orbit of Uranus. Galle (1846) later found the planet only 1 o Le Verrier?s prediction. The next hundred years of advances in the eld were mathematical in nature, and many new studies on di erent kinds of perturbations to the motion of planets and satellites were conducted. Most of the analytical works during this period were focused on three-body problems (e. g., expansion of the disturbing function by Boquet, 1889), or low-order approximations for systems with a few more objects and additional perturbations (e.g., secular frequencies in the Solar System by Brouwer et al., 1950). Darwin (1879, 1880) also began to pioneer the analysis of the lower-order e ects of tides and tidal friction. The use of computers for numerical integration opened a new window on the subject in the 1960s, and made it possible to handle more complicated systems for a long period of time and to study the formation and evolution of the 5 whole Solar System. One key numerical integration of the outer Solar System for 120,000 years was undertaken by Cohen and Hubbard (1965). The development of numerical techniques in celestial dynamics (within the So- lar System, to be speci c) was recently reviewed by Morbidelli (2002), who divided the numerical study of Solar System dynamics into four major periods. During the classical period, when slow but accurate integration algorithms (Runge-Kutta and Bulirsch-Stoer methods) were used, Cohen and Hubbard (1965) veri ed the analyti- cal theory about secular interactions by Brouwer et al. (1950). Sussman and Wisdom (1988) found the chaotic nature of Pluto?s orbit through a 845-million-year integra- tion on a speci cally constructed parallel computer called the Digital Orrery. A great e ort had also been given to understand resonance structure (Wisdom, 1983; Murray, 1986; Wisdom, 1987) and the distribution and stability of asteroids (Milani et al., 1989; Nesvorn y and Ferraz-Mello, 1997). The symplectic period in numerical Solar System studies started with the e cient Hamiltonian-preserving algorithm proposed by Wisdom and Holman (1991). A symplectic scheme for the solution of the equa- tions of motion for a Hamiltonian system is able to bound the energy error with a large timestep, typically 10-20 samplings per orbital period, thus enabling faster and longer integrations. With this new algorithm, research was carried out to understand the evolution and stability of the whole Solar System (Sussman and Wisdom, 1992; Murray and Holman, 1999), as well as small-body dynamics (Holman and Wisdom, 1993; Duncan et al., 1995). The most important discovery during this period is the phenomenon characterized as chaotic di usion. In this regime, particles are stable 6 for billions of years, but are eventually able to escape due to long-term weak pertur- bations (Morbidelli, 1997; Nesvorn y and Roig, 2001), overturning earlier beliefs that particles in the Solar System are either unstable in a short period of time or stable forever. The Wisdom and Holman (1991) symplectic scheme reaches its limit when close encounters are involved, because the perturbing Hamiltonian dominates over the Ke- plerian potential during a close encounter, which violates one of the algorithm?s fun- damental assumptions. This problem was rst overcome for zero-weight particles by Levison and Duncan (1994), opening the \statistics period" when systematic studies of comets and near-Earth asteroids (NEAs) was made possible. Based on statistics of the lifetime of the Jupiter-family comets and the Neptune-encounter rate of Kuiper belt objects, Levison and Duncan (1994, 1997) concluded that the Kuiper belt con- tains 6:7 109 comet-sized bodies. Duncan and Levison (1997) showed further that the Kuiper belt?s scattered disk was 100 times more populated at the beginning of the Solar System. With similar techniques, Morbidelli and Gladman (1998) and Michel et al. (2000) studied the distribution of NEAs from di erent sources, the results of which were later used by Bottke et al. (2000, 2002) to construct a model for orbital and magnitude distribution of NEAs. Finally, integrations of planetary accretion involving encounters and collisions among massive planetary bodies and planetesimals was enabled with the design of two algorithms: Symba (Duncan et al., 1998) and Mercury (Chambers, 1999). For the rst time, the formation and evolution of the Solar System could be studied over its 7 entire 4.5 billion years history. For the inner Solar System, Chambers and Wetherill (1998, 2001) and Agnor et al. (1999) studied the formation of terrestrial planets from lunar-mass \planetary embryos", and found that Earth-sized planets can be formed between 0:5 2:0 AU in several hundred million years, but that original embryos in the asteroid belt are mostly scattered into unstable orbits and only 1% of the population is left behind in the main belt (Petit et al., 2001). Progress was also made on the formation and evolution of the giant planets. It was known (Fernandez and Ip, 1984) that interactions with planetesimals cause giant planet to migrate: Jupiter goes inwards, while all the other three giant planets move outwards. This procedure produces certain observed structures in the Kuiper belt (Malhotra, 1993, 1995) and is responsible for Pluto?s strange orbit. Recent development of the \Nice model" suggested that the four giant planets may have formed in a much tighter group: between 5:5 and 14 AU. This model was the rst one to simultaneously reproduce the current orbits of the giant planet (Tsiganis et al., 2005), the Trojan population of Jupiter (Morbidelli et al., 2005) and Neptune (Sheppard and Trujillo, 2006), the source for the Late Heavy Bombardment of the terrestrial planets (Gomes et al., 2005), and the structure of the Kuiper belt (Morbidelli et al., 2007). Our analytical understanding of the physics of orbital interactions and evolution also saw many advances during the computer era, with symbolic algebra and semi- numerical techniques in many cases. Mean-motion resonances (Section 2.3) are by far the most popular subject. Early reviews on satellite orbital resonances were given by Greenberg (1977) and Peale (1986), who also analyzed resonant encounters dur- 8 ing tidal migration. Two basic types of behavior of the orbital elements during an encounter were identi ed: kicks and trapping. Henrard (1982) and Borderies and Gol- dreich (1984) computed the capture probabilities during a resonant encounter. More recently, Hamilton (1994) revisited the resonant encounter problem, comparing mean- motion resonances arising from gravitational and electromagnetic perturbations. The analysis of the more subtle secondary resonances was rst done by Malhotra (1990), and the e ect may be responsible for breaking resonant trapping in many cases. Tidal e ects in the Solar System were also investigated in great detail by Goldreich (1963) and Goldreich and Soter (1966). A review by Burns (1986) covers most development on how tides a ect orbits. For non-resonant, or secular, evolution (long-term orbital interactions), once the secular frequencies of the Solar System were found (see, e.g., Brouwer et al., 1950), research has focused on the e ects of secular resonances (Ward et al., 1976) on asteroids (Froeschle and Scholl, 1987; Scholl and Froeschle, 1990), on Kuiper belt objects (Nagasawa and Ida, 2000), and on other minor bodies in the Solar System including Trojan asteroids (Marzari and Scholl, 2000). The study of planetary rings makes another important branch of the Solar System dynamics. It has a long history that dates back to 1610 when Galileo rst pointed his telescope toward Saturn. His discovery was later interpreted by the Dutch as- tronomer Christiaan Huygens as \a thin, at ring, nowhere touching, and inclined to the ecliptic". In many regards, Saturn?s ring is the most splendid phenomenon in the Solar System, and remained the only known ring system until the 1970s and 1980s when a urry of ground-based observations and space missions found new ring sys- 9 tems around Uranus (Elliot et al., 1977), Jupiter (Owen et al., 1979), and Neptune (Smith et al., 1989). Rings around Mars were rst suggested by Soter (1971) and their properties were predicted by Hamilton (1996) and others. Ring dynamics expe- rienced renewed interest with the discovery of a wide variety of rings around the giant planets. Saturn?s F-ring and the Uranian rings are narrow with sharp edges, unlike the broad main rings of Saturn. Jupiter?s entire ring system is composed of micron- sized dust, and there are additional examples of dusty rings belonging to each of the other giant planets. A review by Cuzzi et al. (1984) showed that various features in the broad main ring of Saturn result from gravitational interactions between ring particles and embedded moonlets and among the ring particles themselves. Burns et al. (1984) reviewed the e ects of electromagnetic perturbations with a focus on the vertical structure of the Jovian rings. A follow-up by Hamilton (1994) provided more details on how planetary magnetic elds interact with micron-sized ring parti- cles. Showalter and Burns (1982) showed numerically that sharp ring edges could be con ned and wakes and spokes inside the narrow rings might be excited by nearby shepherd moons. Since 2004, the Cassini spacecraft has sent back the most detailed data on the Saturnian rings, which con rms some earlier theories about ring-moon interactions, and also shows a number of interesting new features. This has led to active and ongoing numerical and analytical exploration of ring dynamics. Several recent books (Planetary Rings by Esposito (2006), Planetary Ring Systems by Miner et al. (2007)) and articles (e.g. Porco and Hamilton, 2007), cover both observational and theoretical developments in planetary ring systems. 10 Research on orbital dynamics is not limited to the Solar System. With the discoveries of dusty disks (Aumann, 1985) and extrasolar planets (Wolszczan and Frail, 1992; Mayor and Queloz, 1995), a lot of e ort has been spent in understand- ing the di erent orbital statistics of planets (the Extrasolar Planet Encyclopaedia, http://exoplanet.eu/) and planet-disk interactions (Ward, 1988). The most signif- icant di erence between the orbits of extrasolar planets and planets in the Solar System is that the former usually have much larger eccentricities. Several mecha- nisms were proposed to produce these eccentric orbits, including Kozai resonances (Holman et al., 1997), planet-planet scattering (Ford et al., 2001), planet-disk inter- action (Goldreich and Sari, 2003), and stellar encounters (Zakamska and Tremaine, 2004). In Part III of this dissertation, we will study how very close giant planets (\hot-Jupiters") may retain their eccentricities against tidal circularization. Studies in planet-disk interactions are usually orientated towards understanding the features seen in the spatially-resolved disks (Holland et al., 1998; Wilner et al., 2002; Greaves et al., 2005), and towards predicting possible planetary masses and orbits based on these disk patterns (Kuchner and Holman, 2003). 11 Chapter 2 Perturbation Theory Gravity from a massive central body usually determines planetary and satellite sys- tems, resulting in elliptical orbits with the dominate mass at one focus, as stated by Kepler?s rst law. The largest planet in the Solar System, Jupiter, is a thousand times less massive than the Sun, and the most extreme satellite-planet mass ratio, excluding Pluto-Charon, is about 1% for the Moon-Earth pair. Thus, any e ects from forces other than central gravity can be treated as small perturbations to the otherwise Keplerian orbits. In this chapter, we summarize the basics of perturbation theory, with emphasis on results relevant to the later chapters. 2.1 De nition of Orbits In a system of two perfect spherical bodies (ideal planets, stars, satellites, etc.), the orbit of each object is an ellipse with their common center of mass at one focus, as is determined by the Keplerian potential. In the case of a hierarchical star-planet or planet-satellite system, the orbit of the secondary body is usually measured in a body-centric frame with the dominant mass at the center, and is also a closed ellipse. We will call the dominant mass a planet (mp) and the secondary mass a satellite (ms) throughout this chapter. 12 mp f i W w a(1?e) ms Reference Direction Pericenter Reference Plane Figure 2.1: De nition of orbital elements. A Keplerian orbit is a closed ellipse with the dominate mass at one focus. The location of the secondary in space can be determined by 6 orbital elements: the semi-major axis a, eccentricity e, inclination i, longitude of ascending node , argument of pericenter !, and true longitude f. The longitude of pericenter $ = +! is a bent angle measured in two di erent planes. The elliptical orbit of a satellite is well-de ned geometrically by ve constant Keplerian orbital elements (Fig. 2.1): the semi-major axis a, the eccentricity e, the inclination i, the longitude of ascending node , and the argument of pericenter !. The last one is often replaced by the longitude of pericenter $ = + !, which is a bent angle measured partly in the reference plane and partly in the orbital plane. The semi-major axis a is sometimes replaced by orbital mean motion n, which is the average angular velocity of the satellite: n = r G(mp +ms) a3 ; where G is the gravitational constant. To determine the location of the satellite in space, a sixth element is required to describe its position along the orbit. This element takes many forms; for example, the true anomaly f is the angle between the 13 directions from the central star to the planet and its orbital pericenter (Fig. 2.1). A more commonly used angle, although without a simple geometric representation, is the mean longitude , which is a longitude measured in two planes from the reference direction, much like $. It has the important property of increasing linearly in time. In contrast to the other ve elements, the true anomaly or mean longitude changes periodically over time to represent a dynamical system instead of a static one. With two or more satellites in the system, or if the shapes of the objects deviate from perfect spheres, satellite orbits are not closed eclipses any longer. In the most common cases, the additional forces are much weaker than the Keplerian potential and the orbits di er only slightly from perfect eclipses. For each instant in time, an imaginary orbit for each satellite can still be de ned by the so-called osculating elements transformed from its location and velocity vectors, assuming a Keplerian orbit (Murray and Dermott, 1999, x2.9). As a result of the extra perturbations, the rst ve osculating elements (a, e, I, , and $) vary only slowly with time, while the true anomaly f (and similarly the mean longitude ) still change quickly. In some aspects of our work, we will average the orbits over f to study the long-term evolution of the orbital size, shape, and orientation. 2.2 Disturbing Function If a system contains more than one satellite (Fig. 2.2), the orbits can behave in a rather complicated manner because of mutual perturbations between the two objects, as we pointed out in last section. It is not di cult to write down the equations of 14 motion of the satellites based only on Newton?s laws: r1 = G(mp +m1)r1r3 1 Gm2 r 1 r2 jr1 r2j3 + r2 r32 =rr1K1 +rr1R1; (2.1) r2 = G(mp +m2)r2r3 2 Gm1 r 2 r1 jr1 r2j3 + r1 r31 =rr2K2 +rr2R2: (2.2) Here r1 and r2 are the position vectors of the inner satellite m1 and the outer satel- lite m2, respectively, and K1 and K2 are the Keplerian potentials due to the central planet, which alone would cause each satellite to orbit in an ellipse. The two addi- tional potentials,R1 and R2, are the disturbing functions arising from the satellites? perturbation on each other; they can be written explicitly as R1 = Gm2 1 jr1 r2j r1 r2 r32 ; (2.3) R2 = Gm1 1 jr1 r2j r1 r2 r31 : (2.4) Despite the concise look of Eqs. (2.1-2.4), it r2 r1 r2 r1? m1 m2 mp Figure 2.2: A two-satellite sys- tem. Position vectors are mea- sured in the planet-centric frame. By convention, the inner satellite has subscript \1", and the outer one has subscript \2". is conceptually simpler to consider the evolution of orbits in terms of orbital elements rather than the more rapidly-varying position and velocity vectors. A signi cant e ort has been devoted to- ward this end, starting with Peirce (1849) who derived a series expansion of Eqs. (2.3) and (2.4) to sixth order in the eccentricities and mutual in- clination. Le Verrier (1855) published the most commonly-used expansion to seventh- order, which was later expanded to eighth order by Boquet (1889). All these expan- sions are carried out in terms of the mutual inclination and ascending node, which 15 simpli es the calculation and is useful for many situations. When strong perturba- tions from rotational deformation of the planet (Section 2.4) or from a massive foreign object exist, however, inclinations measured from the planet?s equatorial plane or the Laplace plane (Chapter 4) are more useful. Thus an expansion in individual orbital inclinations is necessary. Such an expansion to fourth order can be found in the appendix of Murray and Dermott (1999), and Murray and Harper (1993) have per- formed an error-free eighth-order expansion relying extensively on two independent computer algebra codes. Once we have written a disturbing function in terms of orbital elements, we can determine the e ect of the corresponding perturbation on the orbit by solving La- grange?s planetary equations: d dt = 2 na @R @a + p1 e2(1 p1 e2) na2e @R @e + tan(i=2) na2p1 e2 @R @i ; (2.5) da dt = 2 na @R @ ; (2.6) de dt = p1 e2 na2e (1 p1 e2)@R@ + @R@$ ; (2.7) di dt = tan(i=2) na2p1 e2 @R @ + @R @$ 1na2p1 e2 sini@R@ : (2.8) d$ dt = p1 e2 na2e @R @e + tan(i=2) na2p1 e2 @R @i ; (2.9) d dt = 1 na2p1 e2 sini @R @i : (2.10) The angle = nt is de ned as the mean longitude at epoch. Brouwer and Clemence (1961) and Danby (1988) both give detailed derivations of these equations. Disturbing functions are not limited to two-satellite systems. For a system with many satellites, each satellite raises perturbation potentials, in the form of Eqs. (2.3) 16 or (2.4), on all other satellites, and the disturbing function of a satellite can be calculated by summing perturbation potentials from all other satellites. In addition, any other potentials leading to additional disturbing forces can be treated in this way, as we shall see in Section 2.4. 2.3 Secular and Resonant Perturbations After expansion in terms of small quantities e and i, a disturbing function can be written as a sum of a series of cosine harmonics: R= X (jk) R(jk) cos(j1 1 +j2 2 +j3$1 +j4$2 +j5 1 +j6 2) (2.11) with the strength of each harmonic given to lowest order by R(jk) = (a1;a2)ejj3j1 ejj4j2 (sini1)jj5j(sini2)jj6j: (2.12) Here j1; ;j6 are integers subject to the d?Alembert relations: j1 +j2 +j3 +j4 +j5 +j6 = 0; (2.13) j5 +j6 = even number: (2.14) Hamilton (1994) provides an intuitive derivation of these rules based on spatial sym- metry. Finally, (a1;a2) is the part of strength that is independent of orbital ec- centricities and inclinations. It, however, depends on the six integers jk in addition to the two semi-major axes. Ellis and Murray (2000) devised a method to nd the strength factor for a certain argument = j1 1 +j2 2 +j3$1 +j4$2 +j5 1 +j6 2: (2.15) 17 For orbital evolution studies, the fast changes of 1 and 2, which specify the locations of the satellites, are less important than the change in orbital size (a), shape (e), and orientation (i, , and $). Hence, we average the disturbing function over the two orbital periods. Consequently, most terms in Eq. (2.11) average to zero with the following two exceptions: i) Secular terms: j1 = j2 = 0; or ii) Resonant terms: j1 6= 0 or j2 6= 0, but j1 1 +j2 2 = constant: (2.16) Terms that meet condition (i) cause secular perturbations, which exist for any two interacting orbits. Since the strength of a perturbation term is proportional to powers of e and sini (Eq. 2.12), for orbits with small eccentricities and inclinations, terms with large jj3j, jj4j, jj5j, or jj6j only have weak e ects. For the lowest-order terms with (j3 = 1; j4 = 1) or (j5 = 1; j6 = 1), the Lagrange equations can be linearized and an analytical solution can be found, leading to the rst-order secular theory. Secular perturbations between planets cause periodic changes of orbital eccen- tricities and inclinations, as well as precessions of the ascending nodes and pericenters (Chapter 7.2). Stockwell (1873) was the rst to calculate the secular variations for all eight major planets, which was later improved by Brouwer et al. (1950). Their method is similar to what we will use in Part III of this dissertation, where we study how external eccentricity damping a ects secular interactions. Terms satisfying condition (ii) lead to mean-motion resonances; taking the time 18 derivative of Eq. (2.16), we nd j1n1 +j2n2 = 0: (2.17) This condition can only be met when the two mean motions have a ratio of integers (n2=n1 = j1=j2), hence the name \mean-motion resonance". In general, a resonance happens when a dynamical system is driven at one of its natural frequencies. Here when a \commensurability" exists among the orbits (Eq. 2.17), the same orbital con guration repeats and small perturbations between the satellites can accumulate in phase. Following Greenberg (1977), we examine the resonance mechanism and stability for a simple case in which two planar orbits are in 2:1 mean-motion resonance (j1 = 1 and j2 = 2). Similar analyses can also be found in Peale (1976, 1986) and Murray and Dermott (1999). The system is shown in Fig. 2.3: for simplicity we assume that i)the inner orbit is circular, ii) the outer one is eccentric, and iii) m1 m2 so that the inner orbit is not a ected. The perturbation from m1 results in a net tangential force f on m2. During the period when the two move from conjunction to opposition, f is in the direction of motion and accelerates m2. During the period from opposition to next conjunction, f is in the opposite direction of motion and decelerates m2. If the conjunction occurs exactly at pericenter or apocenter of the outer orbit, the tangential force averages to zero, and there is no energy or angular momentum exchange. As a result, the orbital con guration does not change and the two satellites always line up at pericenter or apocenter. Conjunctions at any other location destroy this symmetry and result in a non-zero 19 averaged force that adjusts the conjunction location. Fig. 2.3 shows the case in which the two satellites line up just past m2?s pericenter. When the two satellites move from conjunction to opposition, the average distance between them is larger than when they move from opposition to next conjunction since m2 passes its apocenter in the former case. As a result, the tangential force on m2 during the accelerating phase is weaker than during the decelerating phase. Furthermore, the system spends less time moving from conjunction to opposition (slightly less than one period of m1) than from opposition to next conjunction (slightly longer than one period of m1). Hence, during a full period between two successive conjunctions, the weaker accelerating force is applied for a shorter time than the stronger decelerating force, resulting in a net loss of energy and angular momentum for the outer orbit. The outer satellite falls inward and its mean motion, n2, increases. The inner satellite then has more di culty catching up with the outer one, so subsequent conjunctions must move towards the apocenter. Similarly, conjunctions occurring during the second half of the outer orbit, when m2 moves from its apocenter to pericenter, also shift towards m2?s apocenter. Thus, the outer orbit?s apocenter is a stable conjunction point, while its pericenter is an unstable one. It can also be shown that conjunctions at the inner orbit?s pericenter are stable, while those at its apocenter are unstable. Thus, during a resonance, extreme close encounters (e.g. conjunction at the pericenter of the outer orbit and the apocenter of the inner orbit) are avoided and the orbits are stabilized, as the case of the 3:2 resonances between the Plutinos and Neptune. The repetition of the same orbital con guration, however, allow the orbital elements to be 20 systematically perturbed, leading to changes of orbital elements on short timescales. Since mean-motion resonances m p m 1 m 2 Pericenter Conjunction Apocenter f Opposition Figure 2.3: Stability of mean motion resonances. During a resonance, the conjunction point always moves towards the apocenter of the outer orbit and towards the pericenter of the inner orbit. are associated with an expansion term of the disturbing function, the angular parameter for that term, commonly called the reso- nant angle or resonant argument, can be used to characterize a res- onance. We rewrite Eq. (2.15) in the form usually used for mean- motion resonances: = (p+q) 2 p 1 +j3 $1 +j4 $2 +j5 1 +j6 2: (2.18) The new coe cients p = j1, q = j1 + j2, and all integers are still subject to the d?Alembert conditions (Eqs. 2.13 and 2.14). The sum jj3j+jj4j+jj5j+j6j, which is usually equal to jqj, is the order of the resonance since it represents how many small quantities (e1, e2, sini1, sini2) appear in the resonant strength Eq. (2.12). There are two possible behaviors for the resonant angle : circulation through a full 360 when the two orbits are distant from the corresponding resonance, or libration through a restricted range of values when the resonance is nearby. The libration amplitude of decreases to zero as the resonance is approached, and at 21 exact resonance, the resonant argument satis es _ = 0: (2.19) If the orbits do not precess, i.e., _ 1 = _ 2 = _$1 = _$2 = 0, this equation leads to the resonant condition Eq. (2.17): resonances occur when the two orbital mean-motions are an exact ratio of integers. In reality, however, both the oblateness of the planet (Section 2.4) and secular and resonant perturbations from other satellites cause orbits to precess, leading to resonance splitting qualitatively similar to the Zeeman e ect in which the energy levels of an atom separate when a magnetic eld is applied. Since the nodal and apsidal precession rates in Eq. (2.18) are much slower than orbital mean motions, these resonances are still packed into a small region around the location determined by the ratio of the satellite mean motions (Eq. 2.17). In Part II of this dissertation, we will investigate the resonant history of the inner Neptunian system and study how the orbits of these moons have been a ected over time. 2.4 Rotational Deformation The spin of a planet induces a bulge to appear along the planet?s equator { for example, Earth?s equatorial radius is about 22 km bigger than its polar radius { and this deviation from a spherical shape causes perturbations to the Keplerian orbit of a satellite. Roy (2005) used the disturbing function and Lagrange?s planetary equations to derive the low-order e ects for small eccentricities and inclinations, which is usually su cient for satellite dynamics. We follow his approach here. In the ideal situation, 22 the rotational deformation is axisymmetric, and the perturbation potential, or the disturbing function, is given by Rrot = Gmpr 1X k=2 Jk R p r k Pk(sin ); (2.20) where Rp is the planetary radius, r is the distance from the center of the planet, is the latitude measured from the equatorial plane, Pk(sin ) is the Legendre polynomial of degree k, and the dimensionless coe cients Jk are determined by the planet?s shape. An arbitrary axisymmetric potential can be represented by a judicious choice of the dimensionless coe cients Jk, which represent the magnitudes of the di erent harmonics. If a planet is symmetric about its equator plane, we have J3 = J5 = J7 = = 0. To an excellent approximation, this is the case for all gaseous planets. More general expansions are needed for terrestrial planets (Hamilton, 1994) whose mass distributions are usually non-axisymmetric. Assuming that the orbit di ers only slightly from a Keplerian ellipse, r and sin can be converted to orbital elements by r = a(1 e 2) 1 +ecosf; sin = sinisin(f +!): Substituting these expressions into Eq. (2.20), dropping terms of J6 and higher, and averaging over an orbital period, we get, to second-order in e and sini, hRroti= 12n2a2 " 3 2J2 R p r 2 98J22 R p r 4 154 J4 R p r 4# e2 12n2a2 " 3 2J2 R p r 2 278 J22 R p r 4 154 J4 R p r 4# sin2 i: (2.21) 23 Table 2.1: Gravitational properties of major planets Planet R (km) J2(10 5) J4(10 5) k2 Q Mercury 2,440 6 - - 190 Venus 6,052 0.4 0.2 0.25 17 Earth 6,378 108.3 -0.2 0.299 12 Mars 3,394 196 -1.9 0.14 86 Jupiter 71,398 1473.6 -58.7 0.58 (0:6 20) 105 Saturn 60,330 1629.8 -91.5 0.40 16;000 Uranus 26,200 334.3 -2.9 0.36 11,000 - 39,000 Neptune 25,225 341.1 -3.5 0.41 9,000 - 36,000 References: The Love numbers k2 for the giant planets are from Bur sa (1992); their tidal Q?s are from Table 6.2 of this dissertation; Q?s of Mercury and Venus are from Goldreich and Soter (1966); all other data are from Murray and Dermott (1999). Since neither nor $ appears inhRroti, planetary oblateness does not a ect e and i (see Eqs. 2.7 and 2.8). Nevertheless, it causes precessions of both the ascending node and the pericenter. The precession rates come from combining Lagrange?s equations (Eqs. 2.9 and 2.10) with Eq. (2.21): _$ = +n " 3 2J2 R p r 2 98J22 R p r 4 154 J4 R p r 4# ; (2.22) _ = n " 3 2J2 R p r 2 278 J22 R p r 4 154 J4 R p r 4# : (2.23) The signs of the two expressions indicate that, for a prograde satellite, the pericen- ter precesses, while the ascending node regresses. For a retrograde satellite, however, the opposite is true, as we shall see in Section 4 for Neptune?s large satellite Triton. Table 2.1 lists the J2 and J4 parameters for the major planets. Giant planets have much larger J2 and J4 values due to their faster rotation rates. 24 2.5 Tides Planets and satellites are not rigid f 2 f 1 f Satellite Planet d Figure 2.4: Tidal torque. Because the near tidal bulge pulls the satellite more strongly than the far one does, the net force f does not point exactly toward the center of the planet, resulting in a tangential component of the total force along the direction of the satellite?s velocity. bodies. Because the gravitational pull from a satellite di ers at di erent loca- tions around the planet, the shape of the planet is changed: it elongates along the planet-satellite line and forms two tidal bulges (Fig. 2.4). In exactly the same way, a planet also raises tides on its satellites. Tidal deformation perturbs orbits of satellites just as planetary oblateness does; the e ect was rst discussed by Darwin (1879, 1880), and later expanded and systematized by Kaula (1964) and MacDonald (1964). Following Goldreich and Soter (1966) and Burns (1977), we show how tides a ect orbits physically. An excellent summary of tides and tidal interactions is also given by Murray and Dermott (1999). If a satellite moves along a circle at the same angular rate as the spin of its planet, the tidal bulge on the planet always aligns with the planet-satellite line. As a result, the total gravitational torque on the satellite is zero and the satellite?s orbit is not a ected. When this happens, the satellite is synchronized. Given a planet?s spin rate p, its synchronous radius for satellites can be easily calculated: rsyn = Gm p 2p 1=3 : 25 When a prograde satellite is away from the synchronous orbit, however, the plan- etary tidal bulge can apply a net torque on the satellite as illustrated in Fig. 2.4. This gure shows the case when a satellite is outside of the synchronous orbit, where it orbits more slowly than the planet spins. Because of internal friction, tides take some time to develop and the planet?s fast spin carries along the tidal bulge ahead of the satellite. Subsequently, because of the di erence between the gravitational pulls from di erent side of the planet (Fig. 2.4), a tangential component of the gravitational force is applied in the direction of the satellite?s motion and accelerates it. As energy is pumped into the satellite?s orbit, the orbit expands away from the planet and its mean motion decreases. On the other hand, since energy is drained from the planet, lost both to the satellite and to internal frictions, the planet?s spin slows down. The opposite is true if a satellite is inside the synchronous orbit: the tidal bulge lags behind the satellite and energy is transferred from the satellite?s orbit to the planet?s spin and to tidal dissipation. The satellite migrates towards the planet, its orbital velocity speeds up, and the planetary rotation rate increases. Finally, if a satellite is in a retrograde orbit in which it revolves in the opposite direction of the planet?s spin, like Neptune?s Triton, the tidal bulge always lags behind the satellite, and the satellite?s orbit always decays inward. The tidal migration rate of the satellite and the despin rate of the planet have been calculated by many authors. Here we adapt the equations given by Murray and 26 Dermott (1999): _n = sign( p n) 92k2pQ p R p a 5 sns; (2.24) _ p = sign( p n) 15 4 k2p Qp R p a 3 2sn2: (2.25) Here s = ms=mp 1 is the satellite-to-planet mass ratio, and k2p is the Love number of the planet (Table 2.1), which is a dimensionless parameter characterizing a planet?s tides-raising ability and determined only by the planet?s internal structure. It is usually less dependent on planetary composition than Qp is and can be estimated based on planetary models (see Gavrilov and Zharkov, 1977; Bur sa, 1992). The tidal dissipation factor Qp quanti es the ability of the planet to dissipate energy; a smaller Qp means stronger tidal friction and higher energy loss rate. Qp generally depends on the amplitude and frequency of tides (Goldreich and Soter, 1966), but this dependence is thought to be very weak for the low-frequency tides with small amplitudes, expected on most planets and satellites. There is a simple relation between Qp and the tidal lag angle (Fig. 2.4): sin 2 = 1=Qp: Note that Eqs. (2.24) and (2.25) change sign at the synchronous orbit as described earlier. The dissipation factor Qp is usually estimated through dynamical constraints (see Goldreich and Soter, 1966; Yoder and Peale, 1981; Peale et al., 1980; Tittemore and Wisdom, 1989; Ban eld and Murray, 1992). In Chapter 6 we use new dynamical arguments to constrain Neptune?s Q. 27 E ects from satellite tides are similar: _n = sign( s n) 92k2sQ s R s a 5 1 s n 2 (2.26) _ s = sign( s n) 15 4 k2s Qs R s a 3 1 2s n 2; (2.27) where Rs is the satellite radius, and k2s and Qs are the Love number and tidal dissi- pation factor of the satellite, respectively. Among the four rates Eqs. (2.24 - 2.27), _ s is by far the largest (Table 2.2). Thus, satellites usually despin to a synchronous state so that one face is locked toward the planet as is the case for our Moon. All regular satellites in the Solar System are currently spin-synchronized. As a result, satellite tides are not important for circular orbits, and satellite migrations are mostly deter- mined by planetary tides. Table 2.2 shows the tidal migration timescale of satellites. Satellites of all giant planets have probably migrated by less than a few planet radii over the age of the Solar System. Because of small Qp values for terrestrial planets, both the Moon and Phobos have migrated substantially. During tidal migration, the planet continues to transfer angular momentum and energy to the satellite?s orbit un- til the planetary spin is also synchronized fully to the satellite mean motion and the system reaches double synchronization, as is the case for the Pluto-Charon system. Satellite tides, however, can circularize an eccentric orbit very e ectively even when the satellite has reached near-synchronization. This is because tides on the satellite are stronger when it is near pericenter and weaker when it is near apocenter, so tides rise and fall with the orbital period. Because the tidal force is nearly radial, the orbital angular momentum is conserved, but the orbital energy decreases. Thus 28 Table 2.2: Tidal timescales for natural satellites Satellite s=_ s (yr) e=_e (yr) Rp=_a (yr) i=_i (yr) Moon 2 107 2 1010 2 108 2 1011 Phobos 3 105 5 109 1 108 1 109 Io 2 103 6 106 1 1010 2 1011 Europa 4 104 3 108 2 1011 1 1013 Mimas 1 103 3 108 5 1010 5 1010 Titan 3 104 2 109 5 1011 3 1013 Miranda 8 103 3 108 2 1010 4 1011 Ariel 1 104 6 107 1 1010 2 1011 Triton 4 104 9 107 1 1010 7 1011 Proteus 1 103 9 107 1 1010 2 1011 Larissa 2 102 9 107 2 1010 1 1011 These timescales are rough estimates: we use k2p and Qp from Table 2.1; k2s and ~ s are estimated by Eq. 2.29 and ~ s (104 km=Rs)2 given in Murray and Dermott (1999); Qs are all assumed to be 100 except for the Moon, for which Yoder (1995) estimates QM = 27. a satellite?s eccentricity is forced to decrease at a rate given by Goldreich (1963): _e e = 63 4 1 ~ sQs s R s a 5 n; (2.28) where ~ s is the ratio between elastic and gravitational forces, another measure of the internal strength. For a homogeneous solid body, k2s and ~ s are related by k2s = 3=21 + ~ s : (2.29) Tidal circularization timescales for satellites are listed in Table 2.2. The eccentric- ities of most satellites have damped away over the age of the Solar System, except, again, for those of the Moon and Phobos. The former is too far away for Earth to raise strong tides while the latter is simply too small for tides to be e ective. Je reys (1961) found that, in addition to causing outward satellite migration, planetary tides also act to increase the orbital eccentricity. The e ect can be under- 29 stood by considering the tidal force to be applied as an impulse at pericenter { if a rises, e must rise as well. Goldreich (1963) showed that, however, this e ect is usually much weaker than eccentricity damping by satellite tides. Orbital inclination can also be a ected by tides because the rotation of the planet shifts the tidal bulge o the satellite?s orbital plane. Kaula (1964) found that, for small orbital tilts, inclinations change at nearly the migration rate: 1 i di dt = 1 4a da dt: (2.30) Since satellite inclinations are usually very small, this e ect can be ignored in most cases, as the long timescales in Table 2.2 indicate. For most satellites, the changes of inclination are less than a tenth of their current tilts over the history of the Solar System. In Part II, we use this evidence to argue that all past inclination excitations of the small Neptunian satellites are retained until today. 30 Part II Orbital Resonances in the Neptunian System 31 Chapter 3 Background The investigation of the resonant history among the small inner Neptunian satellites is motivated by the study of the resonant excitation of the inclinations of Amalthea and Thebe by Hamilton and Proctor (unpublished). In this chapter, we provide some background information about the Neptunian satellites and their tidal migration history, as well as about the numerical tools and techniques we will use to carry out our simulations. In Chapter 4, we de ne two proper orbital elements that are necessary to analyze resonances in the system. Then, in Chapter 5, we detail various Proteus resonance passages and discuss their immediate implications. Finally, in Chapter 6, we give constraints on several physical parameters for the system. Part of this research is published as Zhang and Hamilton (2007); the remainder has also been submitted to Icarus for publication. 3.1 The Neptunian Satellites Prior to the Voyager 2 encounter, large icy Triton and distant irregular Nereid were Neptune?s only known satellites. Triton is located where one usually nds reg- ular satellites (close moons in circular equatorial orbits, which are believed to have formed together with their parent planets). The moon follows a circular path, but 32 its orbit is retrograde and signi cantly tilted, which is common only among irreg- ular satellites (small distant moons following highly inclined and elongated paths, thought to be captured objects). Triton?s unique properties imply a capture origin followed by orbital evolution featuring tidal damping and circularization. Although di erent capture mechanisms have been proposed (McKinnon, 1984; Goldreich et al., 1989; Agnor and Hamilton, 2006), in all scenarios Triton?s post-captured orbit is ex- pected to be remote and extremely eccentric (e> 0:9). During its subsequent orbital circularization, Triton forced Neptune?s original regular satellites to collide with one another, resulting in a circum-Neptunian debris disk. Most of the debris was probably swept up by Triton ( Cuk and Gladman, 2005), while some material close to Neptune survived to form a new generation of satellites with an accretion timescale of tens of years (Ban eld and Murray, 1992). Among the survivors of this cataclysm are the six small moonlets discovered by Voyager 2 in 1989 (Smith et al., 1989, see Fig. 3.1). Voyager 2 also found several narrow rings interspersed amongst the satellites within a few Neptune radii, and found the ring arcs hinted at by stellar occultation years earlier. Karkoschka (2003) reexamined the Voyager images later, and derived more accurate sizes and shapes of the new satellites. Proteus, the largest one, is only about 400 kilometers in diameter, tinier than even the smallest classical satellite of Uranus, Miranda. Owen et al. (1991) used Voyager data to calculate the orbital ele- ments of these small satellites, which were later re ned by Jacobson and Owen (2004) with the inclusion of recent data from the Hubble Space Telescope and ground-based observations. Both analyses show that all the small moons are in direct near-circular 33 Table 3.1: Inner Neptunian satellites and Triton Name R (km) a (km) e ( 10 3) iLap ( ) ifr ( ) Naiad 33 3 48,227 0:4 0:3 0.5118 4:75 0:03 Thalassa 41 3 50,075 0:2 0:2 0.5130 0:21 0:02 Despina 75 3 52,526 0:2 0:2 0.5149 0:06 0:01 Galatea 88 4 61,953 0:04 0:090 0.5262 0:06 0:01 Larissa 97 3 73,548 1:39 0:08 0.5545 0:205 0:009 Proteus 210 7 117,647 0:53 0:09 1.0546 0:026 0:007 Triton 1353 354,759 0.0157 - 156.83 Average radii of the small satellites ( R) are from Karkoschka (2003); their orbital elements (semi-major axis a, eccentricity e, inclination of local Laplace plane iLap, and free inclination ifr) are from Jacobson and Owen (2004). Neptune?s equator plane is tilted by = 0:5064 from the invariable plane (see Fig. 4.1); these small satellites lie nearly in the equator plane. Orbital data of Triton are from Jacobson et al. (1991). orbits with small, but non-zero, inclinations. Their parameters are listed in Table 3.1. Smith et al. (1989) estimated Figure 3.1: Voyager 2 images of the six small in- ner Neptunian satellites. Their sizes are roughly to scale. Images courtesy of NASA. the cometary bombardment rate near Neptune and pointed out that, of the six small satellites, only Proteus was likely to sur- vive disruptive collisions over the age of the Solar System. The in- nermost and smallest satellite, Naiad, might not last much longer than 2 to 2.5 billion years, while the intermediate objects might have been destroyed during an early period of heavy bombardment. In any case, all six small satellites probably formed only after Triton?s orbital migration and circularization was nearly 34 complete and the large moon was close to its current circular tilted retrograde orbit (Hamilton et al., 2005). Triton?s large orbital tilt induces a strong forced component in the inclination of each satellite?s orbit (see Chapter 4). These forced inclinations de ne the location of the warped Laplace plane, about whose normal the satellite?s orbital plane precesses. Measured from their local Laplace planes, the current free orbital inclinations (ifr) of the small satellites are only a few hundredths to tenths of a degree, with one exception (4:75 for tiny Naiad). It is the contention of this work that these free tilts, despite being very small, arose from dynamical excitations during orbital evolution. Physically, the debris disk from which the satellites formed should have damped rapidly into a very thin layer lying in the warped Laplace plane. Satel- lites formed from this slim disk should initially have free inclinations ifr 0:001 , perhaps similar to the thickness of Saturn?s ring. A reasonable explanation for the current non-zero tilts of the satellites thus requires an examination of their orbital evolution history. 3.2 Tidal Evolution and Mean-Motion Resonance Passage Tidal friction between a satellite and its parent planet determines the satellite?s orbital evolution over a long time span (Darwin, 1880; Burns, 1977). The physical e ects of tides have been discussed in Section 2.5. Satellite orbits migrate due to planetary tides according to _a a = 3k2N QN s R N a 5 n; (3.1) 35 where we have rewritten Eq. (2.24) in terms of _a, and the sign in Eq. (3.1) depends on the satellite?s location relative to the synchronous orbit. For the Neptunian system, the synchronous orbit lies between the orbits of Proteus and Larissa, which means that Proteus? orbit has expanded over time while the orbits of the other satellites have shrunk. The large gap between the orbits of Proteus and Larissa may be evidence for this divergence (Fig. 3.2). The tidal migration timescale is rather Figure 3.2: Semi-major axes of the small inner Neptunian satellites in planetary radii. The two arrows indicate the direc- tion of tidal migration. di cult to estimate. Although k2N = 0:41 is theoretically determined by Bur sa (1992), the uncertainty in QN prevents a precise calculation (Goldreich and Soter, 1966). Ban eld and Murray (1992) es- timated 12;000 < QN < 330;000, leading to timescales uncertain by more than an order of magnitude. Here we note that the distances between the two satellites and the synchronous orbit are 1:3RN for Proteus and 0:4RN for Larissa, implying that they have migrated by no more than RN over the age of the Solar System. Triton, due to its retrograde orbit, spirals slowly inward with a typical timescale 1010 years (Table 2.2); this inward drift can be safely ignored in resonant studies. Due to tidal migration, orbital periods of the small satellites change over time, causing satellites to pass through mean-motion resonances (Section 2.3). The be- havior of the orbital elements during resonance passage depends on whether the two orbits are converging or diverging. For rst- and second-order resonances (jj3j+jj4j+ 36 jj5j+jj6j= 1 or 2 in Eq. 2.18), if the two orbits diverge from each other (as Proteus and Larissa) and pass through a resonance, the orbital eccentricities and inclinations are subject to sharp changes or kicks (Hamilton and Burns, 1993), which can be either positive or negative. The signs and magnitudes of these kicks depend not only on the resonant strength, but also on the exact phase (value of ) when the resonance is encountered (Peale, 1986). However, kick amplitudes are predictable if the two satellites diverge so slowly that the variation of orbital elements is in the adiabatic limit both before and after a resonance encounter. In this case, the kicks to eccentric- ities and inclinations are always positive, and the kick magnitudes can be obtained analytically by a Hamiltonian analysis (Peale, 1976; Murray and Dermott, 1999). In contrast, when two converging orbits pass through a resonance, they can be captured into a resonant state and may remain locked therein unless perturbations from other objects or nearby resonances force them out (Greenberg et al., 1972; Malhotra and Dermott, 1990). If tides continue to act on objects trapped in a resonance, the af- fected eccentricities and/or inclinations keep growing on the tidal migration timescale (Hamilton, 1994). In addition to tidal migration of a, tides a ect other orbital elements as well (Eqs. 2.28 and 2.30). Based on reasonable assumptions for QT and ~ T, Goldreich and Soter (1966) estimated that the circularization timescale for Triton is of order 108 years. Triton, therefore, has followed a nearly circular path for most of Solar System history. The eccentricity-damping timescales for the small satellites are longer because of their small sizes, but are still signi cantly shorter than a billion years (Table 2.2), 37 which is consistent with the zero eccentricities of Naiad, Thalassa, Despina, and Galatea (Table 3.1). The eccentricities of Proteus and Larissa, however, stand out as signi cantly larger than zero. We will explain in Section 5.1.1 that excitations from the recent 2:1 Proteus-Larissa mean-motion resonance passage are responsible. For a satellite with a small tilt, the tidal e ect on the inclinations is very weak (Eq. 2.30). The inclinations of the inner Neptunian satellites should change by less than a tenth of their current values over the age of the Solar System (Eq. 2.30). As a result, any inclinations acquired by other means, e.g., resonant excitations, are preserved throughout tidal evolution. This leads directly to the main idea of this work: the small satellites? current orbital tilts contain clues to dynamical events of the past. We begin our investigation by an overview of possible resonances that may have once been active among the Neptunian satellites. 3.3 Resonant History of the Neptunian System As we mention earlier, Proteus migrates away from Neptune, while all other satellites spiral inwards. A satellite?s migration rate is proportional to its mass, and has a steep inverse dependence on its orbital semi-major axis (Eq. 3.1). Since Larissa, Galatea, and Despina have comparable masses (within a factor of 2), the innermost of these migrates most rapidly, and, hence, their orbits all diverge from one another. Diverging orbits usually lead to resonant kicks during which satellite orbital elements change sharply (Hamilton and Burns, 1993). Due to the much smaller masses of the innermost satellties, Naiad and Thalassa, their orbits evolve more slowly and are 38 approached by those of the next three satellites. Converging orbits typically lead to resonant trapping, as has been suggested as the explanation for the large inclination of Naiad by Ban eld and Murray (1992). Over the history of the Solar System, the inner Neptunian satellites migrated slowly enough that most of the rst- and second-order mean-motion resonances were traversed in the adiabatic limit. Magnitudes of resonant kicks on satellite inclinations are then predictable both analytically (Peale, 1976) and numerically. Higher-order resonances are often not transversed adiabatically, but their kicks are smaller by about an order of magnitude, and thus add negligibly to the total inclination growth. Because satellite mean motions are much larger than orbital precession rates, resonances cluster in discrete narrow zones near where the ratio of satellite orbital periods is a rational number ((p + q)=p in Eq. 2.18). We integrate Eq. (3.1) for each satellite and locate all the rst- and second-order resonant zones (q = 1 and q = 2), as shown in Fig. 3.3. We stop the integration when Larissa is fairly close to the synchronous orbit and the unphysical discontinuity in _a=a becomes problematic. Because of this over-simpli cation for satellites near synchronous orbit, the left-hand- side of the plot is less accurate than the right, especially for Larissa. As the time axis indicates, the integration is much longer than the age of the Solar System, but the evolutionary timescale depends on Neptune?s Q and the satellite masses, all of which are poorly constrained. Here we have usedQN = 20;000 and a common satellite mean density = 0:6g=cm3 based on estimates for the giant planets and icy satellites. Under these assumptions, Proteus was 0:28RN closer to Neptune 4 billion years ago 39 Figure 3.3: Possible past rst- and second-order mean-motion resonances among Pro- teus (P), Larissa (L), Galatea (G), and Despina (D). We integrate Eq. (3.1) backward until Larissa is fairly close to the synchronous orbit (SO), assuming QN = 20;000 and a uniform satellite mean density of = 0:6g=cm3. Black solid curves show the migration tracks of the four satellites, and the dashed horizontal line denotes the synchronous orbit. The vertical lines represent strong resonant zones for di erent pairs of satellites. The time scale along the bottom axis depends on QN and . For di erent values of these parameters, multiply all times by a factor of QN=20;000 =(0:6g=cm3). ( age of the small satellites) than it is today and Larissa, Galatea, and Despina have migrated towards the planet by 0:24RN, 0:39RN, and 0:49RN, respectively. A larger QN or a lower satellite density would result in a slower evolution, as described in the caption of Fig. 3.3. Thus the origin of the system could occur anywhere along the horizontal axis of Fig. 3.3, depending on the actual values of QN and . With our assumptions of QN = 20;000 and = 0:6g=cm3, the satellites go through approximately 16 resonant zones involving rst- and second-order mean- motion resonances (Fig. 3.3) since the small satellites formed 4 billion years ago. 40 For di erent QN and , a di erent number of past resonances may have occurred. Our strategy here is to work backwards in time from the present to nd out which of these resonances actually occurred and determine what their e ects were. We focus on orbital inclinations which are particularly una ected by tides and hence accumulate over time, leaving a \fossil record" visible in today?s orbits. 3.4 Computing Techniques Our simulations were carried out with the HNDrag module in the HNBody package (Rauch and Hamilton, 2002). HNBody is a general purpose, hierarchical N-body in- tegrator, which implements both the symplectic mapping algorithms (Wisdom and Holman, 1991) and the classical Bulirsch-Stoer and Runge-Kutta algorithms. HNDrag expands the functionality of the original HNBody code by allowing additional drag forces to act on the satellites, which can simulate a wide range of gravitational and non-gravitational perturbations. Since our interest lies in long-term orbital evolution, we use the symplectic integrator for better performance. The integration stepsize is chosen so that there are at least 20 steps during each orbital period. We have per- formed convergence tests for several of our simulations with the number of sampling points per orbit ranging from 1 to 100. The results are consistent for all tests with greater than ve steps per orbit. In the results presented here, we use a cautious 20 steps per orbit to guarantee convergence. We also con rmed the stability of the code by performing a series of simulations with slightly di erent initial conditions. We conducted the test with a system consisting of a planet and two satellites, with 41 an arti cial drag force pulling the satellites through several mean-motion resonances. These tests are similar to the physical problem that we are exploring, and they lead to consistent results after 108 years. The output of HNDrag can be set to either osculating orbital elements or Carte- sian positions and velocities. The osculating elements are a set of projected Keplerian orbital elements for each instant, calculated with the assumption of no extra pertur- bations. However, perturbations from both Neptune?s oblateness and Triton cause the osculating elements to vary arti cially over a single orbital period. We minimize this e ect by using geometric elements, which de ne the actual shape of the orbit. Following Greenberg (1981), we take the position and velocity output from HNDrag and convert it to geometric orbital elements, correcting for rst-order J2 perturba- tions with our conversion program cj2. This procedure greatly reduces unphysical oscillations in the orbital elements. To determine the evolutionary history of the two Neptunian satellites, it would be best to follow their orbits for 4 billion years. As this is not practical with current computing technology, we take advantage of the fact that mean-motion resonance passages take place only at discrete locations. During most of the evolution when the moons are not in resonance, we apply the tidal evolution equations (Eqs. 3.1 and 2.28) to damp eccentricities and move satellites away from the synchronous orbit. Typical resonance passage times, with the slowest migration rate that we use, are on the order of 10 million years; we only simulate these 10-million-year segments, which greatly reduces the computational burden. 42 The simulated system consists of Neptune (with an equatorial bulge), Proteus, and Larissa, with Triton included (for our inclination study) or excluded (for our eccentricity study). We ignore the Sun in our simulations because its perturbation on the small satellites is much smaller than Triton?s. For simplicity, we x the semi- major axis of Larissa, and apply a simple drag force on Proteus to move it slowly outward across the resonant zone. In reality, both satellites are moving at time- dependent rates. But since most of the strong resonances are traversed slowly (in the adiabatic limit), the kicks to the orbital eccentricities and inclinations are independent of whether one satellite or both are migrating, the rate of migration, and even the nature of the drag force. Most of our simulations are performed on the Borg Beowulf cluster of 85 pro- cessors in the Astronomy Department at the University of Maryland. HNDrag is a single-thread program, and di erent simulations are dispatched to di erent nodes of the cluster through the Condor job control system. With these resources, typical simulation times range from 2 days to 2 weeks. 43 Chapter 4 Perturbations from Neptune?s Oblateness and Triton The orbital con guration of the Neptune system is di erent from other satellite sys- tems in that a massive moon, Triton, orbits far from the equator plane, and its perturbation on the inner orbits is not negligible. Hence, a mathematical preparation for this special system is necessary before we proceed to analyze the resonances. If Neptune were perfectly spherical, the rotational angular momentum of Neptune (LN) and the orbital angular momentum of Triton (LT) would both be constant with xed directions in space. In reality, however, the oblateness of Neptune resulting from spin deformation causes Triton?s orbital plane to precess slowly (Section 2.4). For a circularly-orbiting Triton, the nodal precession rate is (Danby, 1988, x11.15) _ T = 3 2J2nT R N aT 2 cosiT goblT : (4.1) Compared to Eq. (2.22), we only include the second-order J2 term, and the additional cosiT is necessary because Triton?s large inclination and retrograde orbit. Neptune?s J2 = 0:003411. Triton?s ascending node precess rather than the more common re- gression given in Eq. (2.22) because iT > 90 . The precession period 2 _ T is about 600 years, signi cantly longer than Triton?s 5.88-day orbital period. Although LT is no longer a constant vector due to the precession of Triton?s orbital plane, the system still conserves its total angular momentum Ltot = LN + LT and, as 44 Figure 4.1: Laplace plane of the Neptune-Triton system. The plot shows a side view of the invariable plane, Neptune?s equatorial plane, Triton?s orbital plane, and the local Laplace plane of a small satellite. Here, iT and are the inclinations of Triton?s orbit and Neptune?s equator, respectively. Note that they are measured from di erent sides of the invariable plane due to Triton?s retrograde orbit. The inclination of the small satellite?s local Laplace plane is given by iLap. The thin curve de nes the shape of the warped Laplace plane for satellites at di erent distances, or for a debris disk inside Triton?s orbit. The whole Laplace plane precesses together with Triton?s orbit and Neptune?s equator. a result, the plane perpendicular to Ltot is xed in space, which makes it a natural reference plane for the measurement of orbital elements. This plane is usually referred to as the invariable plane. In the Neptune-Triton system, it is tilted by = 0:5064 from Neptune?s equatorial plane (Jacobson and Owen, 2004). Neptune?s equatorial plane is always locked with Triton?s orbital plane and the two precess together. We ignore the spin angular momentum of Triton and the orbital angular momenta of the other satellites since they are much smaller than jLNj and jLTj. Small inner Neptunian satellites (ms mT mN) experience secular perturba- tions both from Neptune?s oblateness and from Triton. The overall e ects of these two perturbing components force the orbit of a small moon to precess about the moon?s 45 local Laplace plane, which is distinct from both the invariable plane and Triton?s or- bital plane. Fig. 4.1 shows the warped Laplace plane in the Neptunian system. Near Neptune, the Laplace plane is close to the planet?s equatorial plane, near Triton it is close to the large moon?s orbital plane, and in between the plane is tilted at di erent angles. The nodes of Laplace planes at di erent distances, however, all lie along a line and move slowly with Triton?s secular precession rate. Thus the whole warped disk precesses as a rigid body along with Triton?s orbit and Neptune?s equator. The location of the local Laplace plane at di erent distances from the central planet can be determined by analyzing the two competing perturbations. We undertake this analysis here, as it will lead to the identi cation of the new resonances that we will encounter in Chapter 5. Neptune?s oblateness causes the orbit of a small satellite to precess with a rate gobl given by an expression similar to Eq. (4.1). Triton, as an external perturber, also causes both the satellite?s pericenter and node to precess through secular interactions. We will study the eccentricity secular perturbations in great detail in Part III. In the Neptune-Triton system, however, the eccentricity e ects are trivial due to Triton?s nearly-circular orbit, but the inclination e ects are important. Following similar calculations laid out in Section 7.2, we can obtain the corresponding nodal precession rate: gsec = 14 Tn 2b(1)3=2( ): Here T is the Triton-Neptune mass ratio, = a=aT is the semi-major axis ratio of the satellite and Triton, and b(1)3=2( ) is one of the Laplace coe cients, which depend 46 only on (Murray and Dermott, 1999, x6.4). Combining both the perturbations from Neptune?s oblateness and the secular e ects of Triton, we obtain the disturbing function for a small satellite: R= na2 1 2(g sec +gobl)i2 gsec( iT)icos( T ) ; where the extra symbols are due to Triton?s retrograde orbit. The solution to the corresponding Lagrange?s planetary equations for inclination (Eq. 2.8) and ascending node (Eq. 2.10) with the above disturbing function is isin = ifr sin fr +iLap sin Lap; (4.2) icos = ifr cos fr +iLap cos Lap; (4.3) where fr = (gsec +gobl)t+ fr0 : The free inclination ifr and the free node at the epoch fr0 are constants deter- mined by the initial state. The angles iLap and Lap de ne the local Laplace plane of the satellite, as illustrated in Fig. 4.2a. The inclination of the local Laplace plane, also called the forced inclination, is iLap = g sec gsec +gobl ( iT); (4.4) and the node of the local Laplace plane, or the forced node, is Lap = T + ; (4.5) both of which are independent of the initial inclination and node of the satellite. Once the satellite?s semi-major axis is given for a nearly-circular orbit, the satellite?s local 47 a b Figure 4.2: De nition of key orbital elements. iLap, Lap: inclination and longitude of ascending node of the local Laplace plane; ifr, fr: free inclination and node of the satellite?s orbit measured relative to its local Laplace plane; i, : inclination and node of the satellite?s orbit measured relative the invariable plane. The longitude of ascending node of the orbit is de ned as the bent angle ~ = Lap + fr measured in two separate planes. a: the physical representation of the planes and orbital elements. b: the phase diagram showing the solutions Eqs. (4.2-4.3). Laplace plane is determined. This plane precesses together with Triton?s orbit and Neptune?s equator. Our solution for the Laplace plane, Eqs. (4.4-4.5), is consistent with that derived by Dobrovolskis (1993) in the case of a solar perturbation on satellite orbits. However, his solution is simpli ed based on the fact that the external perturber is much further away from the planet than the perturbed satellite, which is not the case in the Neptune-Triton system. Fig. 4.2b illustrates the solutions Eqs. (4.2 - 4.3) in a phase diagram of isin versus icos . Perturbations on Triton by Neptune?s rotational bulge cause !OO0 to precess (rotate counterclockwise) about the origin at rate jgoblT j, and perturbations on the small satellite by both Neptune?s oblateness and Triton cause !O0A to regress (rotate clockwise) around O0 at the rate jgsec + goblj. The vector sum of !OO0 and 48 !O0A represents the inclination i and the longitude of the ascending node of the small satellite relative to the invariable plane and an arbitrary reference direction. Measuring the direction of !O0A from the reference direction, we nd the angle ~ = Lap + fr; (4.6) which we rede ne as the longitude of the satellite?s ascending node. Fig. 4.2a shows its physical meaning: a bent angle partially in the invariable plane and partially in the Laplace plane, much like the longitude of pericenter $. The free inclination is the tilt of the satellite?s orbit with respect to its local Laplace plane, and the free node is measured from the node of the Laplace plane on the invariable plane. Fig. 4.3 illustrates how the histories of (i, ) and (ifr, fr) di er. We simulate the orbital evolution of a satellite at 8RN (the satellite illustrated in Fig. 4.1) in the Neptune-Triton system. Measured relative to the Laplace plane, ifr 3:5 is a constant over time and fr regresses at a constant rate. However, measured relative to the invariable plane, i oscillates around iLap 8:5 , and is forced to precess at nearly the same rate as the Laplace plane. If a small satellite is initially in its local Laplace plane (ifr = 0 ), it always stays in the plane and its inclination remains constant relative to the invariable plane. However, if it starts out of its local Laplace plane, it precesses about this plane and its inclination measured from the invariable plane oscillates. Hamilton (1996) noticed similar behavior when studying the orbit of a dust grain around Mars subject to strong solar perturbations (his Fig. 7). The concept of the bent angle ~ can be more intuitively understood through a direct comparison to $ = + !. For an inclined orbit, ! is measured in the orbital 49 Figure 4.3: Inclination and node of a satellite at a = 8RN measured relative to the invariable plane (i; ) and the satellite?s local Laplace plane (ifr; fr). Measured from the Laplace plane, the free inclination, ifr, is nearly constant, and the free node, fr, regresses uniformly with a 50 year period. Orbital elements measured from the invariable plane display a more complicated evolution: here both i and oscillate due to the satellite?s orbital regression, while is also dragged along with Triton?s 600-year orbital precession. plane, while is measured in a reference plane (here the invariable plane). With the addition of Triton, however, there are two dynamically important planes in addition to the orbital plane { the invariable plane about which Triton?s orbit precesses, and the local Laplace plane about which the small satellite?s orbit regresses (Fig. 4.2). Because the local Laplace plane determines the dynamics, fr is measured in that plane, and we require an additional angle Lap to specify the location of the Laplace plane. As with $, we are led to a bent angle de nition (Eq. 4.6). Although not necessary for 50 this work, the de nition of $ must also be updated in the Neptune-Triton system to: ~$ = ~ +!, a perverse bent angle measured in three planes (represented with an equally perverse symbol). Here ! is measured from the ascending node of the orbital plane on the Laplace plane rather than on any other reference plane. For the orbits of Proteus and Larissa, the di erences between ~$ and $ are tiny because their free inclinations are so small. It is safe to replace ~$ with $ in most cases. With the new de nition of the longitude of the ascending node ~ (Eq. 4.6) replac- ing , as well as the new longitude of pericenter ~$ replacing $, resonant arguments de ned in Eq. (2.18) hold the same form for resonances among the small satellites in the Neptune-Triton system as we shall see in Chapter 5. The resonance strengths (Eq. 2.12), however, depend on ifrj rather than ij. 51 Chapter 5 Proteus Resonances The magnitudes of resonant kicks on orbital inclinations depended strongly on the mass of the perturber. Since Proteus is by far the largest satellite, resonances between it and the other satellites are much stronger than those among the three smaller satel- lites. As a rst approximation, therefore, we examine only the Proteus resonances, neglect the weaker ones, and see if we can form a consistent story from this subset of Fig. 3.3. We will return to consider resonances between Larissa, Galatea, and Despina in Section 6.1. In addition to the very recent 2:1 resonant zone between Pro- teus and Larissa, Proteus might have gone through seven other resonant zones that we list backwards from the present: PD 3:1, PG 2:1, PL 5:3, PL 3:2, PD 2:1, PG 5:3, and PL 7:5. In this chapter, we detail these resonance passages, determine the inclination excitation provided by each resonance, and analyze new features found in our simulations. 5.1 The Recent 2:1 Proteus-Larissa Resonance Passage The 2:1 mean-motion resonance between Proteus and Larissa (PL 2:1) is located only about 900 km inside Proteus? current orbit or 600 km outside Larissa?s, implying that the satellites passed through the resonance in the recent past (a few hundred million 52 years ago). The proximity of this resonance suggests a resonant origin for the larger- than-average eccentricities of these two satellites (Table 3.1). In Fig. 5.1, we simulate the passage of Proteus and Larissa through this resonance at roughly the correct tidal migration rate. As our rst step in the investigation, Triton has been excluded from the system so that a typical two-body resonance passage can be observed. We plot the orbital semi-major axes, eccentricities, and inclinations of Proteus and Larissa as the satellites diverge slowly through the resonant zone. The orbital elements of the two moons jump at several locations where di erent individual resonances occur. We name the resonances after the orbital elements they a ect with a capital R to signify the appropriate term in the disturbing function (Murray and Dermott, 1999, x6.9), and mark all of the 1st- and 2nd-order ones in Fig. 5.1. Depending on which orbital elements are most strongly a ected, the resonances can be classi ed as eccentricity- type, inclination-type, or mixed-type. 5.1.1 Eccentricity Evolution during and after the PL 2:1 Passage The eccentricities of the two satellites are shown in the middle panels of Fig. 5.1. The two rst-order eccentricity-type resonances, ReL and ReP, dominate the satel- lites? eccentricity growth. Second-order resonances Re2L and Re2P occur at exactly the same locations, respectively, while ReLeP falls between the two. Larissa?s semi-major axis drops while that of Proteus grows with each eccentricity kick to conserve the energy and angular momentum of the system. If aL is not signi cantly altered by the resonances, then ReLeP would be midway between Re2L and Re2P; we derive a similar 53 Figure 5.1: Proteus and Larissa diverge through PL 2:1. Triton is excluded from this system. Plots show the semi-major axes, eccentricities, and inclinations of the two small satellites. Larissa has a xed semi-major axis at aL = 2:93RN, and Proteus migrates outward with a rate of 3:6 10 10RN=yr. As only the relative divergence rate is important in most cases, it is an excellent approximation to move Proteus alone. We assume satellite densities = 0:8g=cm3. Both satellites begin on circular orbits with inclinations of 0:5590 and 1:0667 measured relative to the invariable plane, respectively. These inclinations are the same as would be forced by Triton were it included in the system. Orbital element kicks due to rst- and second-order resonances are marked in the plots. The unmarked small kicks are due to higher-order resonances. 54 Figure 5.2: Proteus and Larissa diverge through PL 2:1. The system consists of Neptune, Triton, and the two small satellites. Plots show the semi-major axes, ec- centricities, and free inclinations (measured relative to the Laplace plane) of Proteus and Larissa. Larissa?s semi-major axis aL is xed at 2:93RN, while Proteus migrates outward with _aP = 1:8 10 10 RN=yr. The density of the satellites is = 0:8g=cm3. The rst- and second-order resonances are identi ed, including a few strong three- body resonances (ReLiPiT, R , RiPiT, and RiLiT); smaller kicks are due to higher-order resonances. 55 result in Section 5.1.2 for the inclination-type resonances. The amplitudes of resonant kicks depend on the strengths of the resonant per- turbations, which are functions of satellite masses and the instantaneous values of the orbital elements. Since the mass of Proteus is about 10 times Larissa?s, a given resonance (e.g., RePeL) gives a stronger kick to Larissa than to Proteus. The strength of the second-order resonance ReLeP depends on two small eccentricities, so it is much weaker than the rst-order resonances and contributes only about 1/6 of the growth of eP. The tiny kicks to eP before RePeL in Fig. 5.1 are due to higher-order resonances. Our additional simulations with di erent tidal migration rates suggest that the tidal migration rate is slow enough that the rst- and second-order resonances are tra- versed in the adiabatic limit. Higher-order resonances are not traversed adiabatically, so their eccentricity and inclination kicks depend on the drag rate and are di cult to predict. For the 2:1 passage, though, higher-order resonances are weaker by about an order of magnitude and their contributions are minimal (Fig. 5.1). We do not include Triton in our eccentricity studies since its orbit is nearly circular and its perturbation to the small satellites? eccentricities is minimal. We verify this assertion with a direct comparison between simulations with and without the large moon (Figs. 5.1 and 5.2). The masses of Proteus and Larissa are not well-constrained observationally. The higher the masses, the stronger the resonances, and in turn, the larger the eccentricity excitation. Since the small satellites formed from the same circum-Neptunian debris disk, we might expect that they should have similar compositions and densities. We make the simple assumption that both satellites have the same density, and calculate 56 their masses based on their observed sizes. In the simulation shown in Fig. 5.1, we use a mean density of = 0:8g=cm3. The satellites might have a higher or lower density, depending on their composition and porosity. The current eccentricities of the two satellites, 0.00053 for Proteus and 0.00139 for Larissa, place a lower limit on the resonant excitation, which then limits the minimum density of the two satellites. We simulate the resonance passage with a number of di erent assumed mean densities for Proteus and Larissa. These simulations show > 0:05g=cm3 in order for Proteus to acquire an eccentricity eP > 0:00053. With this density, Larissa?s eccentricity is excited to a value signi cantly higher than its current 0.00139. A density of > 0:05g=cm3 for satellites is not particularly a good constraint, but we will derive much better limits after later inclination studies. After the resonance, the satellite orbits must migrate outward while simultane- ously circularizing. If we know the eccentricities of the satellites immediately after the resonance, this provides a constraint on the satellite Q?s, which we now explore. Since tidal migration is determined by planetary tides (Eq. 3.1) and eccentricity damping is mostly accounted for by satellite tides (Eq. 2.28), the ratio between a satellite?s Qs and Neptune?s QN can be estimated based on the satellite?s migration distance and the change of its eccentricity subsequent to the resonant passage: Qs QN = 21 4 1 k2N ~ s N s 2 R N Rs ln(af=ai) ln(ef=ei) ; (5.1) where N and s are the densities of Neptune and the small satellite; the subscripts \i" and \f" indicate initial and nal values of the semi-major axis and eccentricity, respectively. We defer the estimation of Qs to Section 6.3 after we have obtained 57 better constraints on QN. 5.1.2 Inclination Resonances in the PL 2:1 Resonant Zone In addition to the eccentricities of the Proteus and Larissa, Fig. 5.1 also shows the change of the satellite inclinations in the bottom two panels. First-order inclination- type resonances do not exist due to the constraints on resonant arguments (Eqs. 2.13 and 2.14, Hamilton, 1994). The three second-order inclination resonances, Ri2P, RiPiL, and Ri2L, are equally-spaced in time, which can be explained by considering the cor- responding resonant arguments: 0i2 P = 4 P1 2 L 2 P1; (5.2) 0iPiL = 4 P2 2 L L P2; (5.3) 0i2 L = 4 P3 2 L 2 L; (5.4) where the subscripts, 1, 2, and 3, denote the three di erent locations of Proteus. We use 0 instead of here to distinguish these arguments from the new de nitions to be introduced later in this section. Since the three resonant locations are very close, we can safely neglect the di erence between _ P1 and _ P2. Applying Eq. (2.19) and subtracting pairs of equations yield nP1 nP2 nP2 nP3; which, for closely-spaced resonances, is equivalent to aP2 aP1 aP3 aP2: 58 Figure 5.3: Ri2L transit during the PL 2:1 resonance passage and the corresponding resonant angles. Top: the free inclination of Larissa; middle: the traditionally-de ned resonant argument 0i2 L = 4 P 2 L 2 L; bottom: resonant argument with new de nition i2L = 4 P 2 L 2~ L. In the Neptune-Triton system, 0i2 L fails to librate in the vicinity of resonance; instead, i2L is the true resonant argument. Furthermore, since the migration rate of Proteus is nearly constant during the reso- nance passage, these locations are equally spaced in time as well (Fig. 5.1). We continue our investigation by running a simulation that includes Triton (Fig. 5.2). Compared to Fig. 5.1, the eccentricity histories in this simulation show similar fea- tures, with only a few very weak additional kicks arising from high-order, mixed-type resonances. This justi es our neglect of Triton in the previous section. In addition, the tidal migration rate used in Fig. 5.2 is half of that of Fig. 5.1, which demonstrates that the strong resonances of this resonant passage are traversed in the adiabatic limit. The inclinations shown in Fig. 5.2 are free inclinations with superscript \fr", which are de ned in Chapter 4 and directly comparable to those listed in Table 3.1. 59 The pattern of inclination kicks is quite di erent from what is shown in Fig. 5.1. We identify the three traditional second-order inclination-type resonances (Ri2P, RiPiL, and Ri2L) by their positions and spacing (compare with Fig. 5.1). In addition, there are several new and stronger resonances that appear near the standard ones. Ev- idently, Triton has a signi cant impact on the tilts of the small satellites? orbits. Its secular perturbation slightly augments the moonlets? orbital precession rates, but more importantly, it alters the inclination resonant pattern itself. When two satellites pass through a mean-motion resonance, the corresponding resonant argument has a stationary value at the exact resonant location (Eq. 2.19). In our simulations with Triton, however, we notice that the resonant angles of the three second-order inclination-type resonances, as de ned by the standard Eqs. (5.2- 5.4), are not stationary even when the resonant kicks occur. For example, Fig. 5.3 shows the inclination of Larissa during the Ri2L traverse. The traditional resonant angle 0i2 L is plotted in the middle panel, which shows no sign of libration. This problem motivated our theoretical consideration of new orbital elements in Chapter 4. With the mathematical results therein, we now generalize the resonant argument (Eq. 2.18) to = (p+q) 2 p 1 +j3 ~$1 +j4 ~$2 +j5 ~ 1 +j6 ~ 2 +jT T (5.5) for the Neptune-Triton system. The new resonant angles have stationary values at the exact resonant location (bottom panel in Fig. 5.3), supporting our arguments. Note that in Eq. 5.5, we also include Triton?s node to cover three-body resonances in this system, as discussed in next section. Accordingly, the second d?Alembert condition 60 Eq. (2.14) should be changed to j5 +j6 +jT = even number: 5.1.2.1 Three-Body Resonances In a two-satellite system, the inclination evolution during the 2:1 resonance passage is dominated by three equally-spaced, second-order resonances: Ri2L, RiLiP and Ri2P (Fig. 5.1). In the actual Neptunian system with Triton included, however, several stronger kicks appear near the traditional second-order kicks (Fig. 5.2) What are these new resonances? A careful examination of their resonant locations shows that the strongest kicks (labeled RiLiT and RiPiT in Fig. 5.2) are shifted the same distance to the left of RiPiL and Ri2P, respectively, which implies that the resonant arguments of the two new resonances, iLiT and iPiT, can be derived by adding a common term to the corresponding second-order resonant arguments. Because RiLiT only a ects Larissa and RiPiT only a ects Proteus, ~ L cannot appear in iPiT, and ~ P not in iLiT. The locations of the new kicks thus suggest the following resonant arguments: iLiT = 4 P 2 L ~ L T; (5.6) iPiT = 4 P 2 L ~ P T; (5.7) which we verify by noticing their forced librations (Fig. 5.4) immediately prior to the resonant kicks. The node of the Laplace plane appears in both arguments through T, which means that the resonances can be considered to be amongst Proteus, 61 Figure 5.4: Resonant arguments ( iPiT and iLiT) of the three-body resonances RiPiT and RiLiT. The top panels show the free inclinations of Proteus (ifrP ) and Larissa (ifrL ) as they traverse the two resonances; the bottom panels show the corresponding resonant arguments from Eqs. (5.6-5.7). These simulations use similar parameters as in Fig. 5.2, except that Proteus migrates at a slower rate (3:6 10 11 RN=yr). Larissa and the warped rotating plane. When the system is close to RiLiT and RiPiT, the associated angles iLiT and iPiT begin to oscillate around equilibrium points at 180 and 0 , respectively. The libration amplitude decreases and the a ected inclination rises as each resonance is approached. When the resonance is crossed, the free inclination of the a ected satellite is kicked up sharply and the corresponding semi-major-axis jump brings the two out of resonance. The resonant angle ceases to librate and begins to circulate again. Since the Laplace plane is only a mathematical description of Triton?s secular e ects, these new resonances can also be interpreted as three-body resonances among Proteus, Larissa, and Triton, which is why we use T rather than Lap in Eqs. (5.6- 5.7). Three-body resonances are usually weaker than two-body ones because the involvement of Triton as a resonant perturber introduces an extra factor of mT=mN 62 in the expression for resonant strengths. However, this e ect is counter-balanced by Triton?s large orbital tilt. Speci cally, the strengths of the three-body resonant kicks on Proteus and Larissa scale as: RiPiT / mTm N mL mN sini fr P siniT; (5.8) RiLiT / mTm N mP mN sini fr L siniT; (5.9) while those of the respective two-body resonances obey Ri2P / mLm N sin2 ifrP ; RiLiP / mPm N sinifrL sinifrP : The rst pair di er from the second only by a factor of mT mN siniT sinifrP 0:2; implying comparable resonant kicks. Note that this is only a rough estimate, the exact ratio between the strengths of the resonances also depends on the numeric coe cient in each resonant term. This type of resonance is di erent from previously-studied three-body resonances (e.g., the Laplace resonance among the three Jovian satellites: Io, Europa, and Ganymede) in that the third body?s orbital longitude does not appear in the resonant arguments. Nevertheless, Triton?s node is involved in both arguments, implying that its inclination should also be kicked during resonance crossing. This e ect is, however, extremely weak due to Triton?s huge mass. Since resonant locations are mostly deter- mined by the coe cients of the orbital longitudes appearing in the resonant angles, the new resonances are located close to the standard two-body resonances. 63 In general, the resonant argument of a three-body resonance has the form = p1 1 +p2 2 +p3 3 +j1 1 +j2 2 +j3 3 +j4$1 +j5$2 +j6$3; where the integers pi and ji still need to satisfy the two d?Alembert constraints men- tioned before. With di erent integer coe cients, three-body resonances should be thickly packed throughout the region of the inner satellites. In our simulations, how- ever, we fail to locate any that involve the longitude of Triton (i.e., p3 6= 0), from which we conclude that these resonances are very weak. It is unclear why they are so weak since their strengths should scale similarly with satellite masses and inclina- tions as RiLiT and RiPiT. A de nitive explanation would require a Taylor expansion of a 3-body disturbing function similar to what has been done for two interacting satellites (Murray and Dermott, 1999, x6.4), and an examination of the relevant res- onant terms. This, however, is a monumental undertaking beyond the scope of this dissertation. 5.1.2.2 Important Higher-order Resonances By de nition, RiLiT and RiPiT are second-order resonances since their strengths de- pend on inclinations of both Triton and a small satellite. Generally, however, \order" should refer to an expansion over small quantities. Since siniT is not small here, these three-body resonances (Eq. 5.8-5.9) should really be considered as rst-order in incli- nations. But they are much weaker than the rst-order eccentricity resonances due to the extra dependence on mT=mN, and it is better to consider these resonances to be 64 second-order in the small quantities ifrP , ifrL , and mT=mN. We adopt this de nition here. In addition to the two second-order three-body resonant kicks, a few fairly strong higher-order kicks also contribute signi cantly to satellite inclination histories, two of which are identi ed in Fig. 5.2. The strong resonance ReLiPiT occurs right after ReL, and has a resonant argument eLiPiT = 2 P L ~$L ~ P + T: It is a third-order resonance that a ects the eccentricity of Larissa, the free inclination of Proteus, and the inclination of Triton. We expect the strength of this resonance to be of order eL 0:01 times the strength of the second-order RiPiT, but simula- tions show that the two are comparable. Thus, ReLiPiT must have a large numerical coe cient in its strength term that could be derived through Taylor expansion of the three-body disturbing function. Another interesting resonance is marked as R in Fig. 5.2. It occurs almost ex- actly at the location where 2nP nL = 0. Since this resonance a ects the inclinations of both satellites, both nodes, ~ P and ~ L, should appear in the resonant argument. Nodal precession normally should displace the resonant location from the precise 2:1 commensurability. R , however, is not displaced, suggesting that the satellites? peri- centers ($P and $L) must also be involved in the resonant argument. The pericenters are required to explain the lack of o set, since, to rst-order in small eccentricities and inclinations, _ P = _$P and _ L = _$L. A single resonance with all of these 65 properties would need to be at least fth-order, e.g., = 4 P 2 L + ~$P + ~ P ~$L ~ L 2 T: This is surprising, as a fth-order resonance should not be as strong as the second- order resonance RiPiL (Fig. 5.2). A careful examination of resonances in the vicinity of R reveals two pairs of third-order resonances: eLiLi3T = 2 P L + ~$L + ~ L 3 T; ePiPi3T = 2 P L + ~$P + ~ P 3 T; and eLiLiT = 2 P L ~$L ~ L + T; ePiPiT = 2 P L ~$P ~ P + T: Although each individual resonance a ects the orbit of only one small satellite, the two resonances in either pair occur almost on top of each other, and the two pairs themselves are so close that we cannot resolve them in Fig. 5.2. The rst pair is weaker than the latter pair by a factor of sin2 iT 0:15, although the exact factor again depends on the numerical coe cients in their strength expressions. A magni ed look at R with a slower migration rate shows the slightly di erent locations of these four resonances (Fig. 5.5). The tiny o sets between the locations of the resonances in each pair are due to higher-order eccentricity and inclination e ects on the nodal and pericenter precession rates. At the beginning of the simulation, eLiLiT shows large amplitude libration because ReLiLiT is the strongest resonance 66 Figure 5.5: Detail of the R resonance in Fig. 5.2. Larissa?s orbit is xed at 2:931RN and Proteus migrates outwards at 3:6 10 11 RN=yr. Satellite densities are taken to be 0:8g=cm3. The simulation covers a very small vicinity around the location where 2nP = nL, which occurs here at t 1:04 105 year. The plots show free inclinations of the two satellites, together with resonant arguments of four third-order resonances which are marked in the inclination plots and detailed in the text. 67 Figure 5.6: Similar simulation as shown in Fig. 5.5, but the satellites have a slightly higher density = 1:0g=cm3. With a larger density, the resonances are stronger and the inclination excitations behave rather stochastically. The resonance angle plots show that multiple resonances are active simultaneously. The nal inclinations are impossible to predict. 68 in the vicinity. However, the weaker resonances ReLiLi3T and RePiPi3T are traversed rst. As the orbits approach these two resonances, eLiLiT becomes out of phase for libration, and the arguments of two earlier resonances circulate even more slowly. At the resonant locations, these angles reverse their direction of circulation. The argu- ments of the latter pair of resonances behave similarly. Due to their weak resonance strengths, none of the four arguments strongly librates as shown for RiPiT and RiLiT in Fig. 5.4. The two Larissa resonances, which should be stronger due to higher values of e and i as well as mP >mL, display long-range e ects visible as concentration in the resonant arguments around t = 0. The overlap of these resonance e ects is a recipe for chaos, especially if the reso- nances are a little bit stronger, e.g., with larger satellite masses, or if orbits linger in the region due to a slower Proteus migration rate. When two orbits diverge through an isolated resonance, the semi-major axis of the inner satellite decreases, while that of the outer satellite increases. The resulting jump causes the two orbits to diverge from each other more quickly than during tidal migration. If another resonance is in the immediate vicinity, however, the system can be a ected by it before completely leaving the rst resonance, resulting in stochastic behavior. In other words, all res- onances have e ective widths { near resonance e ects emerge before, and continue after, the exact resonant location. Stronger resonances have broader widths. If two resonances are located very close to each other, and if they are strong enough that their widths overlap, temporary capture can occur and the kicks to orbital elements behave somewhat like a random walk. Fig. 5.6 shows the same resonances as Fig. 5.5 69 does, but with satellites just 25% more massive. As the system steps into the rst pair of resonances, the two resonant angles start to librate. But the system does not exit the resonant region quickly and cleanly as in Fig. 5.5. Instead, the two res- onant angles alternate between libration and circulation in a complicated way, and the inclinations are kicked up or down randomly until the system escapes these reso- nances. The second pair of resonances interacts chaotically in a similar manner. The random behavior of the inclinations throughout this region makes it impossible to predict their total excitation. However, given certain migration rates and low-enough satellite densities, these temporary captures only continue for a limited time. In the simulation shown in Fig. 5.6, the maximum inclination gains of Proteus and Larissa are of the same order as the RiPiT and RiLiT kicks. Similar chaotic interactions have also been noticed in simulations of the orbital resonances among the Uranian satel- lites by Tittemore and Wisdom (1988). The existence of these chaotic zones puts an intrinsic limit on how well the orbital histories of Neptune?s small satellites can be reconstructed. 5.2 The Second-order Resonance PD 3:1 The previous resonant zone Proteus traversed consists of the 3:1 resonances with Despina (cf. Fig. 3.3). The 3:1 resonant zone is simpler than the 2:1 one as it contains only the second-order, fourth-order, and other even-order resonances. Compared to the PL 2:1 zone (Fig. 5.2), which has rst-order eccentricity resonances, there are not any strong eccentricity kicks and the lack of third-order resonances makes the 70 Figure 5.7: Proteus and Despina diverge through PD 3:1. The satellites have a density = 0:8g=cm3, and QN = 33;000. Both satellites are initially on circular orbits in their local Laplace planes. Two of the three usual second-order resonances are identi ed (RiPiD and Ri2D), as well as the two three-body resonances, RiPiT and RiDiT. Note that although the e ects of the third second-order resonance Ri2P is too weak to be visible, we indicate its approximate location. inclination evolution much simpler and cleaner as well (Fig. 5.7). The two second- order three-body resonant kicks, RiPiT and RiLiT, dominate the inclination growth. Two of the three traditional second-order two-body resonances can be easily identi ed, and fourth- and higher-order kicks are mostly too weak to have noticeable e ects. The overall PD 3:1 inclination kick on Despina is a little smaller than the overall PL 2:1 kick on Larissa, and the Proteus kick through the PD 3:1 is weaker by a factor of 3. This is due both to the smaller mass of Despina (mD 0:5mL) and to the lack of contributions from odd-order resonances in the PD 3:1. The two sets of strong third-order three-body resonances (ReLiPiT and R ) that contributed signi cantly to the growth of ifrP do not exist in the PL 3:1 and other second-order resonant zones, 71 including the PL 5:3, PG 5:3 and PL 7:5 from Fig. 3.3. 5.3 The PG 2:1 and Diverging Capture The 2:1 resonant zone between Proteus and Galatea is located between PD 3:1 and PL 5:3. This resonance is similar to the rst-order PL 2:1 discussed above; the only di erences come from the di erent mass ratios of the two satellite pairs and their di erent distances from Neptune. The two rst-order eccentricity resonances (ReG and ReP) kick the orbital eccentricities of the two satellites strongly, as can be seen in Fig. 5.8. The inclination kicks due to three-body resonances (RiPiT, RiGiT, ReGiPiT, and R ) are similar in magnitude to those due to the PL 2:1 resonances discussed in Section 5.1. In a small fraction of our simulations such as the one shown in Fig. 5.8, however, we see a new e ect that is surprising at rst glance: resonance capture. In addition to resonant kicks when satellite orbits are subject to sudden and sharp changes, two satellites may also be captured into a mean-motion resonance during tidal migration (Greenberg, 1977). When this occurs, the mean motions and orbital precession rates of the two satellites vary in such a way that the associated resonant angle librates around a stationary point rather than cycling through full 360 rotations. The a ected eccentricity or inclination grows until nearby resonances or other perturbations break the system out of resonance. What is surprising is that previous studies have shown that resonant captures during tidal migration occur when the orbits of two satellites approach each other (e.g. Hamilton, 1994), while our satellites are diverging. Earlier papers, however, all focus on strong rst- and second- 72 Figure 5.8: The PG 2:1 passage with an isolated eccentricity-type resonance capture. After being captured into the RePi2T resonance, Proteus is later captured into the second-order resonance (Ri2P) which a ects its inclination. The second capture is enabled by extremely slow changes to _~$P and _~ P induced by the rst resonance. In this simulation = 0:8g=cm3 and QN = 22;000. 73 order resonances. For some higher-order resonances, temporary capture is possible for diverging orbits. We have found several such resonant trappings in our simulations when the two orbits diverge slowly enough. Nevertheless, trappings are still very rare even at slow migration rates, probably because of the inherent weakness of higher- order resonances. We have run 10 simulations through the PG 2:1 zone with QN =(g=cm3) ranging between 25,000 and 35,000 for di erent satellite densities ranging between 0:4g=cm3 and 1:5g=cm3, but there is only one capture event at = 0:8g=cm3, shown in Fig. 5.8. In this simulation, after the system has gone through the two rst-order eccentricity resonances (ReG and ReP), the two third-order three-body resonances (ReGiPiT and R ), and the rst second-order three-body resonance (RiPiT), the two orbits are captured into a three-body resonance RePi2T, second-order in the small quantities eP and T, with a critical argument ePi2T = 2 P G + ~$P 2 T: (5.10) 5.3.1 Resonant Trapping Condition We now derive the condition in which resonant trapping into RePi2T is possible. We assume that away from resonance, tides force Proteus and Galatea to migrate at rates _adP and _adG, respectively, and that Triton?s orbit is xed in space. Following Hamilton (1994), we derive the rates of change of aP, aG, eP, and iT due to resonant and tidal 74 perturbations to the lowest order in e and i: _aP = _adP 4 G T eP i 2 T nP aP sin ePi2T; (5.11) _aG = _adG + 2 P T eP i 2 T nGaG sin ePi2T; (5.12) _eP = G T i 2 T nP a2P sin ePi2T; (5.13) _iT = 2 G P eP iT nT a2T sin ePi2T: (5.14) Here is the resonant strength of RePi2T. It has the same units as n2a2 and de- pends only on satellite semi-major axes. When the two orbits are in any of the 2:1 resonances, their semi-major axes must follow aP aG = _aP _aG 2 2 3: (5.15) With Eqs. (5.11, 5.12), and Eq. (5.15), we can solve for sin ePi2T at resonant equilib- rium: sin ePi2T = nP a 2 P 2 G T ePi2T _adP=aP _adG=aG 2 + 3p2 P= G : (5.16) For _adP = _adG = 0, the two equilibrium points are at ePi2T = 0 and ePi2T = 180 . They are shifted slightly when dissipation is present. Substituting Eq. (5.16) into Eq. (5.13), we obtain 2eP _eP = _a d P=aP _a d G=aG 2 + 3p2 P= G ; which has the solution eP = e2P0 + _a d P=aP _a d G=aG 2 + 3p2 P= G t 12 ; (5.17) 75 where t is the time elapsed since entering the resonance, and eP0 is the initial eccen- tricity of Proteus. If resonant trapping occurs for a low eccentricity, eP must increase with time and Eq. (5.17) thus requires _adP aP _adG aG > 0; or that the two orbits diverge from each other. Note that this is opposite the usual requirement that the orbits converge for trapping to be possible, yet consistent with the behavior of Fig. 5.8. This interesting result { trapping for diverging orbits { is a direct consequence of the sign of the ~$ term in the resonant angle (Eq. 5.10), as we shall see below. The resonance RePi2T a ects not only the eccentricity of Proteus, but also the inclination of Triton. We can solve for the evolution of iT with Eqs. (5.14) and (5.16): iT = i2T0 2 P T ra P aT _ad P=aP _a d G=aG 2 + 3p2 P= G 12 ; (5.18) where the second term on the right hand side is negative for trapping. Thus, if resonant trapping occurs, eP increases with time while iT decreases. In typical cases when the satellites have comparable masses, the trapping will cease when iT drops to zero. Because of Triton?s huge mass, however, the change of iT is negligible compared to that of eP: iT eP 2 P T ra P aT 0:002: Hence, Triton?s inclination decreases 500 times more slowly than Proteus? eccen- tricity increases and the resonance trapping is stable for a long time. 76 The resonant trapping condition can be easily generalized for an arbitrary reso- nance with an argument de ned by Eq. (5.5). Following similar procedures, we nd that, in order for any relevant eccentricity or inclination to increase, the integer in front of the corresponding node or pericenter angle in Eq. (5.5), ji, must obey ji _ad 2 a2 _ad1 a1 > 0; (5.19) For standard rst-order (q = 1) and second-order (q = 2) two-body resonances, ji must be negative, and thus we recover the standard result: resonant trapping is possible only if _ad2=a2 _ad1=a1 < 0, i.e., for converging orbits. But for resonances with three or more node and pericenter angles involved, like RePi2T shown in Eq. (5.10), some of the ji can be positive, which makes capture into these resonances possible only for diverging orbits, as we have seen in Fig. 5.8. Thus for any isolated mean- motion resonances, capture for converging orbits requires a negative node or pericenter coe cient. For capture from diverging orbits, a positive coe cient is needed. 5.3.2 Evolution in RePi2T We now look at the resonant angle during the resonant trapping (the fth panel of Fig. 5.8). Assuming moderate eccentricities so that the pericenter term in Eq. (5.10) is negligible, ePi2T satis es e Pi2T = 2 _nP _nG; which can be easily transformed into a harmonic equation using Eqs. (5.11) and (5.12): e Pi2T = ! 2 0 sin ePi2T ( _n d G 2 _n d P); (5.20) 77 where _ndP = (3=2)nP _adP=aP and _ndG = (3=2)nG _adG=aG are the tidal drag rates of the two orbits expressed in terms of mean motions. Eq. (5.20) is a standard harmonic oscillator with a libration frequency !20 = 6 T ePi 2 T a2P 2 G + 3p2 P : As shown in Fig. 5.8, ePi2T begins to librate around a stable equilibrium point at 0 as the orbits are rst trapped into the RePi2T resonance. In fact, the divergence forced by tidal dissipation causes the equilibrium point to shift to a small positive value determined by Eq. (5.16). Averaging over a libration, Eqs. (5.11) and (5.12) simplify to _aP = _adP P; _aG = _adG G; where P and G are small positive quantities. These terms resist the tidal divergence and allow the approximate resonance condition (Eq. 5.15) to be maintained. The libration of ePi2T continues to decrease as exact resonance is approached. For slow changes of !0, the system has an adiabatic invariant (Landau and Lifshitz, 1976, x49) !0 cos = constant; where is the libration width, or the amplitude of the oscillating resonant angle ePi2T. As eP increases slowly, !0 grows and the libration width decreases as shown in Fig. 5.8 from 1:2 107 years until 1:6 107 years when a second resonance is activated. 78 5.3.3 Trapping into Ri2P Shortly after capture into RePi2T, the two satellites are trapped into the traditional two-body resonance Ri2P (Fig. 5.8), with a resonant argument i2P = 4 P 2 G 2~ P: The interplay between the two active resonances is quite interesting. With both resonances active, the orbits still appear to diverge, as can be seen from Eq. (5.17) and the fact that eP continues to rise in Fig. 5.8. But, coincident with the second capture, there is an abrupt change in the behavior of ePi2T ( fth panel of Fig. 5.8) whose libration changes from decreasing with time (as expected for an isolated resonance) to increasing with time. This state of a airs continues until t = 2:6 107 years at which point the second resonance ceases to be active (note the attening of the iP curve) and the libration amplitude of ePi2T begins decreasing again. The sharp-eyed reader might have noticed a very subtle change in the density of points for the i2P history (sixth panel of Fig. 5.8) - an oval-shaped feature between t = 1:6 107 and 2:6 107 years that indicates that the Ri2P resonance is active. The oval feature has dark edges (turning points - note the similar dark edges in the fth panel) and a lighter center because the system spends more time near i2P = 0 than near i2P = 180 . The Ri2P resonance is prevented from cleanly librating about i2P = 0 by the more powerful ReGi2T resonance, but the asymmetry in the density of points that it produces is enough to cause the systematic rise in Proteus? free inclination (Fig. 5.8, fourth panel). The libration of the resonant angle i2P initially decreases 79 (from 1:6 107 yr < t < 2:2 107 yr) and then increases again (from 2:2 107 yr i for all i0. Thus, although the amplitude of the inclination kick decreases for increasing initial inclination i0, the energy input actually increases. Since kicks from these three-body resonances dominate the inclination growth and all obey the same equation, the overall change in inclination attained from passage 94 through a full resonant zone roughly follows the same rule. Eq. (6.1) shows that the contribution to inclination growth by a resonant zone passage is more signi cant if those resonances occur earlier in time when the satellites? free inclinations are still small. For example, the PL 2:1 passage can kick Larissa?s inclination up to 0:085 ( = 0:8g=cm3) if it is the only resonant zone the satellite has gone through, but if the satellite already has a 0:1 free tilt excited by earlier resonances, then the kick is only 0:046 according to Eq. (6.1). Now we know how to combine resonant kicks from di erent zones, and are set to compute satellite masses by equating the total resonant kicks to the satellites? current free inclinations. But an immediate di culty is that we do not know how many resonant zones the satellites have passed through since the tidal evolution timescale itself is not well-constrained. Fortunately, however, the current free inclinations of Proteus, Larissa, Galatea, and Despina provide strong constraints on the number of past resonant encounters. The observed free tilts of Galatea and Despina are both about 0:06 , Larissa?s free inclination, at 0:2 , is three times as large, and Proteus has a smaller inclination of 0:026 (Table 3.1). While the rst three satellites have similar masses (within a factor of 2), Proteus is 10 times more massive, making it simultaneously the strongest perturber in the system and the hardest to perturb. For this reason, its free inclination is signi cant although it is the smallest among the four. We now combine the e ects of all resonances discussed in last section, and attempt to construct a migration scenario that can account for the free tilts of all four satellites. 95 We have shown in last section that the single PL 2:1 resonance passage acting alone provides too little current free tilts for Larissa and probably for Proteus as well. Furthermore, the inclinations of Galatea and Despina likely require that these two satellites once had resonant interactions with Proteus, so the earlier resonances, PD 3:1 and PG 2:1, must also have occurred. These three resonances have similar strengths, and the inclination kicks on the smaller satellite are about the same, as discussed previously. If these were the only resonant zones that the system has traversed, the three smaller satellites would have comparable free tilts today. Galatea and Despina do actually have similar inclinations, and Table 6.1 shows that, for an average density 0:4g=cm3 < < 0:8g=cm3, the free inclinations of these two satellite can be excited to close to their current values. Larissa, however, has an inclination three times larger, which requires additional resonance passage(s). In our tidal evolution model, the next strong Proteus resonance is the PL 5:3, which involves Larissa. Thus, these four Proteus resonances might provide a consistent solution. We now test the above hypothesis with our numerical data. We combine the results from Table 6.1, scaled by Eq. (6.1) as appropriate, and plot the nal free inclinations acquired by the satellites for several possible past histories in Fig. 6.4. The curves for each satellite are represented by step functions with steps occurring at resonant zones that involve that satellite. It is important to realize that these curves are not evolutionary tracks, but are rather nal status plots: inclinations at any given point in the past represent the predicted orbital tilts of the satellites today assuming that the system formed at that past time. We only calculate kicks for the 96 Figure 6.4: The cumulative free inclinations of Proteus (ifrP ), Larissa (ifrL ), Galatea (ifrG ), and Despina (ifrD ). The curves show the nal free inclinations that the satel- lites acquire versus their formation time represented in terms of resonant zone pas- sages (c.f. Fig. 3.3). Resonances to the left of the formation time have actually occurred, while those to the right have not. We assume all satellites have the same densities; solid curves represents di erent bulk densities: from top to bottom: = 1:5; 1:2; 1:0; 0:8; 0:6; 0:4g=cm3. The single exception to this is the Despina 0.4 curve which is actually above the 0.6 curve due to stochastic variations. The dashed horizontal lines represent the current free inclinations of the satellites. 97 most recent four Proteus resonances since the chaotic behavior of PL 3:2 prevents any useful estimates of prior inclination growth. The curves for Galatea and Despina are very simple (with only one step) since they each can only be involved in at most one major resonant zone passage. Larissa may go through two major zones, and there are four possible ones for Proteus. If PL 2:1 is the only resonant zone the system has gone through, Fig. 6.4 shows that, with a mean density of 1:5g=cm3, ifrP can be kicked up to near its current value, but Larissa can only obtain about half of its current orbital tilt, as we have discussed in the last section. If the system starts from an earlier con guration to include both the PD 3:1 and the PG 2:1 resonant zones, the inclinations of both Galatea and Despina can be pushed high enough with a density satisfying 0:4g=cm3 < < 0:8g=cm3. For ifrP to reach 0:026 in this case, however, a density of 0:6g=cm3 < < 1:0g=cm3 is required. This three-resonant-zone scenario leaves Larissa with less than half its current tilt. With an earlier formation time so that all four Proteus resonances are traversed, the density required by Proteus?s free inclination drops to 0:4g=cm3 < < 0:8g=cm3 due to the extra kick from Larissa. Although ifrL is more excited, however, Larissa still attains just over a half of its current inclination for this density range. In fact, ifrL can only be excited to 0:14 even if the mean density is as high as 1:5g=cm3. Thus the three or four Proteus resonance passages provide a consistent solution for the free tilts of Proteus, Galatea, and Despina with 0:4g=cm3 < < 0:8g=cm3, a density range that is physically plausible. But the current inclination of Larissa cannot be matched with the assumptions of equal satellite densities in any of the 98 scenarios considered so far. What happened to Larissa? The rst possible solution to the Larissa problem is to allow satellites to have di erent densities. But this leads immediately to the same problem that we faced in last section: in order for ifrL to be kicked more, we need a more massive Proteus, or a less massive Larissa. The mass of Proteus cannot be increased much since this would result in larger inclinations for Galatea and Despina than currently observed. Reducing Larissa?s mass helps a little, but the dependence of the PL kick strengths on Larissa?s mass is very weak since Proteus is the dominate mass. Our simulations show that the cumulative inclination kick on Larissa only increase from 0:1 to 0:13 if we drop L from 0:6g=cm3 to 0:05g=cm3 while keeping P at 0:6g=cm3. Even this unrealistically-low satellite density does not solve the problem. A second possible solution is to allow Larissa to pass through more resonances. The next Proteus-Larissa resonant zone is the chaotic PL 3:2. For this zone, our simulations indicate that chaotic behavior becomes signi cant only for density > 0:8g=cm3, in which case Proteus usually gets an overall kick > 0:025 through the random walk process ). Adding kicks from later resonances, this results in too large a tilt. However, for density 0:6g=cm3, the chaotic behavior is weak and the orbital inclination growths are reasonable. In our example shown in Fig. 5.9 ( = 0:6g=cm3), the kicks on Proteus and Larissa are 0:009 and 0:06 , respectively, which bring the overall inclination growths of the two satellites through the 5 resonance passages to ifrP = 0:029 and ifrL = 0:128 , where we have used Eq. (6.1) to model the later kicks. If = 0:4g=cm3, the PL 3:2 kicks on the two satellites are 0:008 and 0:05 , 99 which can only promote the inclination growth of Larissa through the ve resonant zones to a total of ifrL = 0:107 . Although the actual inclination kicks may vary, our simulations show that for a speci c density, they all fall in roughly this range. Hence, Larissa?s orbit cannot be e ectively tilted even if the troublesome PL 3:2 resonant zone is included. Perhaps there is a rare outcome of the chaotic interactions in which Larissa?s inclination is highly excited, but we have seen no evidence that this is the case. Including earlier resonances is problematic: the preceding one, PD 2:1, excites the inclination of Despina, forcing Proteus? density to < 0:4g=cm2 to keep that satellite?s inclination low and making it less likely to produce similar tilts for Despina and Galatea. The earlier PG 5:3 would then be required. The trend here is clear { additional PD and PG resonances force Proteus? density down, making additional PL resonances less e ective. Tweaking the tidal model might change the order of some resonances (e.g. the PL 3:2 and PD 2:1, cf. Fig. 3.3), but our basic conclusion is unaltered. There is no set of Proteus resonances that can simultaneously make the inclination of Larissa large while keeping those of Galatea and Despina low. A third possibility includes invoking the weaker, but more numerous, resonances among the three smaller satellites (Fig. 3.3), which we have neglected until now. Could the inclusion of these resonances solve the problem of Larissa?s inclination excess? Since the masses of Larissa, Galatea, and Despina are 10 times smaller than Proteus?, these resonances are very weak even though the satellites are closer to one another than to Proteus. Our simulations show that typically inclinations of both satellites show a 0:01 overall growth through a single zone passage, assuming zero 100 initial free inclinations. The cumulative e ect of these kicks could potentially be large given that there are so many of them. However, due to the strong dependence of the kick magnitude on the initial free inclinations, these weak kicks, albeit numerous, do not add much to the satellites? free tilts, especially for the ones occurring after any of the more powerful Proteus resonances. Furthermore, Fig. 3.3 shows that Larissa, Galatea, and Despina have all traversed a similar number of these weak resonant zones { in fact, Galatea receives more kicks than the other two because of its central location { it is nearly impossible to increase Larissa?s inclination signi cantly while keeping those of the other two small. Inclusion of these weak resonances, however, does systematically drive our solution for towards slightly lower values. A fourth possible cause of Larissa?s inclination excess is that the satellite might actually have been captured into a resonance. We have seen an unusual inclination capture following a three-body eccentricity capture with the result that both Proteus? eccentricity and inclination are forced to increase (Fig. 5.8). Although the chance is low because of Larissa?s smaller mass, it is possible that the satellite was once captured into a similar resonance (Ri2L) with either Galatea or Despina. Fig. 3.3 shows that there are several possibilities, with earlier resonances more likely to capture Larissa than later ones because its inclination was smaller in the past. Expanding our search to include the tiny inner moons, Naiad and Thalassa, we note that captures are more likely since Larissa?s orbit is approaching the inner two orbits and resonant trapping into strong low-order inclination resonances is possible without a prior capture. Na- iad?s extremely high free inclination also points to a previous capture (Ban eld and 101 Murray, 1992), and Larissa would be a natural candidate. Assuming that Larissa?s semi-major axis migrated from the synchronous orbit to its current orbit, Larissa?s in- terior 2:1 resonance would move from 2:08RN to 1:83RN, which brackets the current orbit of Naiad at 1:91RN. Although Naiad was somewhat further away from Neptune in the past, there is a good chance that the 2:1 Larissa-Naiad resonance did actu- ally occur. Similarly, Larissa?s 2:1 or 5:3 internal resonances may have swept across Thalassa?s orbit, implying possible strong interactions with that satellite. The odds of capture into a resonance involving Larissa?s inclination is a strongly decreasing function of its initial tilt, and it is hardly possible if Larissa?s orbit is tilted more than some critical value (Borderies and Goldreich, 1984). Hence, the capture probably occurred earlier on. Two main possibilities exist: i) capture occurred prior to the PL 5:3, with two subsequent Larissa kicks from the PL 5:3 and PL 2:1, and ii) capture occurred prior to the PG 2:1 with only one subsequent Larissa kick, the recent PL 2:1. For the density range 0:4g=cm3 < < 0:8g=cm3, scenario (i) requires Larissa to have a free inclination of 0:14 after escaping from the hypothetical resonant capture, while in scenario (ii), the satellite needs to have a 0:16 free tilt after escape. Al- though requiring an additional resonant capture, we think that this is the most likely scenario and will pursue the details in a future publication. Densities in the range 0:4g=cm3 < < 0:8g=cm3 implies large porosities, but similar values are measured for both Saturnian and Jovian satellites. Nicholson et al. (1992) measured the densities of Saturn?s co-orbital satellites, Janus and Epimetheus, at 0:6g=cm3. More recently, Renner et al. (2005) measured the densities of 102 Prometheus (0:4g=cm3) and Pandora (0:49g=cm3), and Porco et al. (2005) estimated the densities of Atlas and Pan to be 0:5g=cm3. All those Saturnian satellites are made of nearly-pure porous water ice. The density of the Jovian satellite Amalthea ( = 0:9g=cm3 Anderson et al., 2005) is a little bit higher, but the satellite is made of a mixture of water ice and rock, so high porosity is also expected. For such low densities, the Roche limit of the Neptunian system is located between 3:1RN and 3:9RN. Accretion for small satellites, however, is possible within the Roche limit through non-gravitational amalgamation. Or perhaps the satellites are fragments from larger bodies that were forced inside the Roche radius and were then disrupted. 6.2 Tidal Evolution Timescale and QN The four satellites have most likely passed though three or four strong Proteus res- onant zones: de nitely the PL 2:1, the PD 3:1, and the PG 2:1, and possibly the PL 5:3 before those. If the system formed with a con guration between PL 3:2 and PG 2:1 in Fig. 3.3, this provides a natural explanation for the inclinations of at least three of the four satellites. We enlarge this region of the plot in Fig. 6.5. Due to the observational error in satellite size measurements (Table 3.1), which leads to mass uncertainties, the evolution tracks of Proteus and Larissa fall somewhere within the lightly-shaded area in Fig. 6.5, resulting in uncertainties of the location of PL 3:2 and PG 2:1 at 4:99Gyr 12;000 by requiring Proteus to form outside of the synchronous orbit and an upper limit QN < 330;000 from the Naiad capture event. As our determined satellite density of 0:4 0:8g=cm3 is half that assumed by Ban eld and Murray (1992), we scale their QN estimates to 4;000 9;000. Our biggest improvement, a factor of 6 times 105 Table 6.2: Q of giant planets Planet Author Q k2 Q=k2 GS66 100;000 1.5 66;000Jupiter YP81 (0:6 20) 105 0.379 (GZ77) 1:6 105-5 106 GS66 60;000 1.5 40;000Saturn P80 16;000 0.341 (GZ77) 45;000 GS66 72;000 1.5 48;000Uranus TW89 11,000 - 39,000 0.104 (GZ77) 105,000 - 375,000 BM92 4,000 - 220,000 0.39 (D88) 10,000 - 560,000 Neptune This work 9,000 - 36,000 0.41 (B92) 22,000 - 90,000 GS66: Goldreich and Soter (1966); GZ77: Gavrilov and Zharkov (1977); P80: Peale et al. (1980); YP81: Yoder and Peale (1981); D88: Dermott et al. (1988); TW89: Tittemore and Wisdom (1989); BM92: Ban eld and Murray (1992); B92: Bur sa (1992); reduction in the upper limit of QN, arises from the constraint that Proteus traversed at least 3 resonant zones. Several authors have previously estimated Q for other giant planets based on similar dynamical constraints; we collect these results in Table 6.2 and compare them to our own. It is best to compare the values of Q=k2 in the nal column because this is the quantity directly constrained by all dynamical studies. Goldreich and Soter (1966) were the rst to systematically investigate the planetary Q?s in the solar system. Based on the known satellites at that time, they estimated lower limits of Q for Jupiter (QJ 100;000), Saturn (QS 60;000), and Uranus (QU 72;000). They ignored the internal structure of the planets and assumed a uniform k2 = 1:5. They also assumed that these satellites could initially migrate from the surface of the planet. The existence of synchronous orbits, however, reduces the amount of tidal migration and lifts their lower limits by a factor of 2 3. Yoder and Peale (1981) 106 addressed this and obtained a lower bound for Jupiter?s Q: QJ 64;000, using k2J = 0:379 computed by Gavrilov and Zharkov (1977). Yoder and Peale (1981) also estimated the upper limit of QJ based on tidal heating required by Io, and obtained QJ 2 106. Peale et al. (1980) corrected Saturn?s Q with more reliable k2S from Gavrilov and Zharkov (1977), and estimated the lower bound to be QS 16;000. The lower bounds of both QJ=k2J and QS=k2S are larger than what we nd for Neptune. Recall that Q is an empirical measurement of the energy dissipation properties of planets that is not well understood physically. For Uranus, the closest sibling of Neptune, Tittemore and Wisdom (1989) placed QU between 11,000 and 39,000 based on the resonant history of the Uranian satellites. However, their QU=k2U value is about 4 times larger than our QN=k2N. Although one might expect these two similar-sized planets to have similar Q and k2, di erent internal structures of the planets { Neptune is denser and has a much stronger internal heat source { may lead to signi cantly di erent Q=k2 values. 6.3 QP and QL The current eccentricities of Proteus and Larissa have signi cant non-zero values (Table 3.1), which we have interpreted as a signature of the recent PL 2:1 resonance passage. In Chapter 5, we derived a constraint on satellite Q based on this interpre- tation (Eq. 5.1). With density between 0.4 and 0:8g=cm3, our numerical simulations show that Proteus and Larissa can attain eccentricities of 0:0011 0:0014 and 0:008 0:01 dur- 107 ing the PL 2:1 passage, respectively. After the resonance, their eccentricities damp to current values of 0:00053 and 0:00139 according to Eq. (2.28). Immediately after PL 2:1, the satellites? semi-major axes must satisfy a3P=a3L ? 4; they then evolve following Eq. (3.1), and the two satellites migrate to their current orbits simultane- ously. Based on these constraints, we calculate the semi-major axes displacements of the two satellites after the PL 2:1 encounter: Larissa has migrated 0:014 0:016RN inward, while Proteus? semi-major axis has increased by 0:010 0:013RN. The un- certainties are primarily due to the observational error of satellite sizes. Substituting these measurements into Eq. (5.1), we nd 0:004 < QPQ N < 0:02; 0:002 < QLQ N < 0:006: We have adopted Neptune?s Love number k2N = 0:41, as computed by Bur sa (1992). The internal strength, ~ s, is unknown for most satellites, but it is not expected to be as sensitive to satellite composition and shape as Qs is. We estimate ~ P and ~ L based on the formula ~ s (104 km=Rs)2 given by Murray and Dermott (1999). For 9;000 0 while s < 0. Thus, Eq. (7.15) states that the pericenters of the two orbits are always aligned in the \+" mode (cos( $) = 1), and always anti-aligned in the \{" mode (cos( $) = 1). In fact, in an eigenmode, the system behaves as a rigid body with the shapes and relative orientation of the elliptical orbits remaining xed while the entire structure rotates in space (Fig. 7.2). The 2-planet secular system is in many ways reminiscent of a double pendulum system. The equations for a simple double pendulum system have the same form as Eq. (7.9), only with a symmetric coe cient matrix A. In both 117 systems, the anti-aligned mode has the faster frequency (g > g+) due to a larger restoring force provided by the spring in the pendulum problem, and the more closely spaced orbits in the planetary problem. We now explore how the eccentricity Figure 7.3: Eccentricity ratios versusp =q in the secular modes as given by Eq. (7.12). Di erent curves represent dif- ferent values. The dashed line represent- ing q =p is a turning point; these orbits have e1 = e2 in both modes. ratios in secular modes depend on the physical parameters q and . We plot 2 as a function of p =q for di erent values in Fig. 7.3. Although rare in real systems, q =p makes an interest- ing special case. When this condition is met, s = 1 (Eq. 7.12), and therefore, the inner and outer orbits have the same eccentricity (Eq. 7.14). When q >p , either the outer planet is more massive (larger q), or the two planets are further apart (smaller ). Then the outer planet has a higher eccentricity in the aligned mode and a lower eccentricity in the anti-aligned mode. The opposite is true for q

p ); b) A 0.3 Jupiter-mass companion is located at 0.2 AU (q

p , or m21a1 < m22a2. The top panel shows results from secular equations, and the bottom panel shows the corresponding N-body simulations. Each panel plots three curves: i) the single planet case, in which the eccentricity of the inner planet damps at rate 1, ii) a two-planet aligned mode with damping rate +, and iii) a two-planet anti- aligned mode (damping rate ). For q>p , eccentricities damp much faster in the anti-aligned mode than in the aligned mode, as predicted by Eq. (8.4). A compari- son between the top and bottom panel shows close agreement (within 2%) between full-scale N-body simulation and integration of the approximate secular equations. Damping rates predicted by Eq. (8.4) also agree well with these observed values. Fig. 8.2b shows a system with m21a1 > m22a2, for which the aligned mode damps faster than the anti-aligned mode. The faster damping rate of the aligned mode in the N-body simulation is still within 2% of the prediction, but that of the slow anti- aligned mode, however, is 10% o . This is probably due to the fact that the drag force that we used in the N-body simulations damps the inner body?s semi-major axis slightly, weakening the coupling between the planets a little. The di erent damping rates for the two modes are particularly interesting, espe- cially for well-separated nearly-decoupled orbits for which and are small. In this case, the 4qp 2 term under the square root of Eq. (8.4) is much smaller than the other term. As a result, if the eccentricity damping on one orbit is much faster than on the other, as the case for tides ( 1 2), one mode damps much faster than the other and the system quickly evolves to a single mode. The more rapid damping rate 125 is near the single-planet tidal rate 1, while the other mode decays substantially more slowly. We now carry out secular integrations and N-body simulations for the same two- planet system (q >p ) shown in Fig. 8.2a, but with initial conditions that lead to a superposition of the two eigen-modes at the beginning. The results are depicted in Fig. 8.3. The two numerical methods agree with each other very well, except for some blurring of the orbital elements in the N-body simulations arising from fast variations at the orbital periods. These oscillations are partially due to impulses at conjunction, and partially due to the subtle di erence between the osculating and geometrical orbital elements as discussed by Greenberg (1981). Before about 4 106 years, the system is in a combination of the two modes, so both eccentricities, as well as their ratio, oscillate (cf. Figs. 7.4 and 7.5). As the short-lived anti-aligned mode damps away, the orbits begin to librate around $ = 0 , the two eccentricities oscillate less and less, and in the end, the eccentricity ratio settles at the aligned mode ratioj +jpredicted by Eq. (8.5). Fig. 8.3 shows the corresponding plots for q

m22a2. Similar to Fig. 8.3, but the planetary system in Fig. 8.2b is shown. The aligned mode damps more rapidly and the system evolves to the anti-aligned mode. 128 apsidal state of the orbits, we plot the aligned and anti-aligned components of the eccentricities for Fig. 8.3 in the top two panels of Fig. 8.5. Recall that the total eccentricities of the orbits at any time can be obtained from the components as illustrated in Fig. 7.4. The lower panels in the gure show the phase diagrams of both orbits (ecos( $) versus esin( $)) at ve di erent points in time indicated in the top two panels. The two orbits travel along the phase curves, which themselves change slowly over time. Two critical instants, labeled s and s+, are circulation- libration separatrices, which represent the transitions of the apsidal state from anti- aligned libration to circulation (s ), and from circulation to aligned libration (s+). These two points divide the evolution curves into three regions. In region I (t< 106 yrs), the anti-aligned components are stronger than the aligned ones for both orbits (e >e+ andj s je > s+e+). In the phase plots, both the e1 and the e2 curves are closed and stay on the left side of the esin( $) axis, indicating the libration of $ about 180 . With the decrease in the amplitudes of all components, especially the faster damping of the anti-aligned ones, the curves move closer toward the origin, resulting in an increased libration width of $. The anti-aligned separatrix s crossing occurs at t = 106 yrs, when the two components for the outer orbit are equal (j s je = s+e+) and e2 may drop to zero, resulting in a phase curve for the orbit that is tangential to the esin( $) axis at the origin (cf. Fig. 8.5). The phase curve for e1 at s is a half-oval whose straight edge includes the origin. Accordingly, when e2 drops to zero, $ jumps from 90 to 90 for the largest possible full libration amplitude of 180 . 129 I a II b III c a b c a b c Figure 8.5: Evolution of the apsidal states during eccentricity damping. The top two panels show the time evolution of the eccentricity components for the system in Fig. 8.3 for which m21a1 < m22a2. The inner orbit has an aligned component e+ and an anti-aligned one e , while s+e+ and j s je are the components for the outer orbit. Two circulation-libration separatrices (s and s+) divide the evolution curves into three parts: anti-aligned libration (region I), circulation (region II), and aligned libration (region III). The bottom panels show phase diagrams of the inner and outer orbits on the complex eexp(i $) plane at the corresponding points indicated in the top panels. The shape of the phase curves depends on the relative strength of the two components for each orbit. Banana shapes result from the large di erence between the two components: e e+ at (a), while j s je s+e+ at (c). 130 As the system moves past s , the phase curves for both orbits enclose the origin and circulation results. The circulation region (II) is located between the two separa- trices (i.e., when 106 yrs e+), while it is weaker for the outer orbit (j s je < s+e+). With the continuous fast damping of the anti-aligned mode, the system crosses the aligned separatrix s+ at t = 3:95 106 yrs, when the two components for the inner orbit are equal (e = e+). The two separatrices occur at those times when each phase curve in Fig. 8.5 touches the origin. After s+, both phase curves are to the right of the esin( $) axis, indicating libration about the aligned mode (region III). The two anti-aligned components are both signi cantly damped and the system now has both e >e+ andj s je > s+e+. The geometry of the orbits can also be illustrated with a component diagram similar to Fig. 7.4, but in a frame rotating at the same rate as the anti-aligned mode (Fig. 8.6). Note that in this rotating frame, the aligned component vectors rotate clockwise because their precessions are slower than those of the anti-aligned ones. Evolution in these coordinates can be visualized as circles whose radii and distances to the origin shrink at the di erent rates + and . Since the anti-aligned components initially dominate the aligned ones (region I), the e1 vector stays on the right side of the esin$ axis, and the e2 vector on the left side. When the aligned components are parallel to the vertical axis, the angle between e1 and e2,j $j> 90 , is at minimum; thus the orbits librate about $ = 180 (anti-aligned libration). When the e2 circle moves to enclose the origin and the e1 circle is still con ned in the rst and fourth 131 vcose e1e 2 sin e v vcose sin e v e1 e2 vcose sin e v e2 e1 Dv Dv (II) (I) (III) Dv Figure 8.6: Eccentricity component diagrams for di erent regions in Fig. 8.5. These diagrams are similar to Fig. 7.4, but now shown in a frame rotating at the same rate as the anti-aligned mode. Here the anti-aligned mode damps faster than the aligned mode ( > +) so that the circles move toward the origin faster than their radii shrink. 132 quadrants, the system reaches circulation region II. This geometry enables $ to cycle through a full 360 (Fig. 8.6). Finally, when the anti-aligned components are su ciently damped and both circles contain the origin, the system goes to aligned mode libration (region III). Now the maximum value of $ < 90 occurs when the two aligned components are parallel to the esin$ axis, and the orbits librate about $ = 0 . Fig. 8.7 shows the case of Fig. 8.4, where the two aligned components are initially stronger and the system starts in the aligned libration region III. The two orbits evolve to cross the aligned separatrix s+ into the circulation region II, and then pass the anti-aligned separatrix s to reach the nal anti-aligned libration region I. The equivalent of Fig. 8.6 for this system would show the radii of the circles shrinking faster than the distances of their center from the origin. In conclusion, the apsidal state of a two-planet secular system depends on the relative strengths of the mode components in the following way: Libration: (e+ e )( s+e+ + s e ) > 0; (8.6) Circulation:(e+ e )( s+e+ + s e ) < 0: (8.7) Libration occurs when the same mode components are stronger for both orbits, and circulation occurs when one mode is stronger for the inner orbit, but weaker for the outer one. Barnes and Greenberg (2006b) studied the non-dissipative secular apsidal state and found the distinction between circulation and libration can be expressed in 133 III a II b I c a b c a b c Figure 8.7: Evolution of the apsidal state during eccentricity damping. Similar to Fig. 8.5, but using data from Fig. 8.4. The system moves from aligned libration to circulation, and nally to anti-aligned libration. 134 terms of the quantity S = ( s + + s )e+e s+e2+ + s e2 : They found that the pericenter di erence $ circulates if jSj> 1, and $ librates if jSj< 1, which is equivalent to our conclusion. Eccentricity damping is e ective in changing the apsidal state of the orbits because the two modes damp at di erent rates. In fact, not only eccentricity damping, but ec- centricity excitation is also able to move the two orbits across the libration-circulation separatrices. This can be easily visualized by running the plots in Fig. 8.5 and 8.7 backwards in time. For systems with m21a1 m22a2 (Fig. 8.7). All mechanisms that change eccentricities slowly cause planetary systems to move toward apsidal libration. 8.3 Relativistic Correction Extrasolar \hot-Jupiters" are close enough to their host stars that general relativis- tic e ects are important. In particular, it is well-known that the post-Newtonian potential produces an apsidal precession of the orbit which, to the lowest order in eccentricities, is (Danby, 1988) _$GR = 3a 2n3 c2 ; 135 where c is the speed of light. We de ne a dimensionless quantity to measure the relativistic e ect: = _$ GR1 ; which is the ratio of the relativistic precession to the characteristic secular precession. General relativity adds an extra term to Eq. (7.3), which now reads _$j = + 1n ja2jej @Rj @ej + _$ GRj : Since _$GR is independent of e for small eccentricities, the extra terms do not change the form of Eq. (7.9), and so all discussion of the general secular modes still holds. In particular, the system still has aligned and anti-aligned modes and the two modes damp separately. The mode frequencies, damping rates, and eccentricity ratios, how- ever, need to be revised. Now the diagonal terms of the coe cient matrix Ajk should be adjusted to A = q + +i 1 q p p (1 + 2 ) +i 2 ; which gives the new mode frequencies and eccentricity ratios: gs = 12 (q + ) +p (1 + 2 ) q [q + p (1 + 2 )]2 + 4qp 2 ;(8.8) s = q + p (1 + 2 ) p[q + p (1 + 2 )]2 + 4qp 2 2q ; (8.9) g = gs q p 2 f[q + p (1 + 2 )]2 + 4qp 2g3=2 ( 1 2) 2; (8.10) = 12 " 1 + 2 p (1 + 2 ) (q + ) p[q + p (1 + 2 )]2 + 4qp 2 ( 1 2) # ; (8.11) = s ( 1 i 1 2p[q + p (1 + 2 )]2 + 4qp 2 ) : (8.12) 136 Note that these equations can be obtained from Eqs. (7.11, 7.12, 8.3, 8.4, and 8.5) with the substitution (q p )![q + p (1 + 2 )]. The relativistic e ect increases both secular rates gs (Eq. 8.8) since it causes the orbits to precess in the same direction as the secular interaction does. Both mode eccentricity ratios are also increased (Eq. 8.9), because the eccentricity of a faster precessing inner orbit is less forced by the outer planet. As for the mode damping rates (Eq. 8.11), relativistic precession decreases the aligned mode damping rate, but increases that of the anti-aligned mode. It also decreases the deviation of the mode apsidal lines from perfect alignment (Eq. 8.12). 137 Chapter 9 Applications to Extrasolar Planetary Systems The purpose of the Wu and Goldreich (2002) study of secular theory and tidal damp- ing was to explain the eccentricity of the planet HD 83443b. At the time there was thought to be a second planet in the system, a claim that has since been retracted. In addition to tidal eccentricity damping, they also considered the small change to the semi-major axis that arises from energy dissipation. They found that, even in a sin- gle apsidal aligned or anti-aligned mode, the eccentricities of the two orbits damp at di erent rates { e1 damps faster because of the semi-major axis migration. For small eccentricities, however, the inward migration, is slower than the eccentricity damping by a factor of e21. Thus we have neglected the semi-major axis migration, which is consistent with the low e assumption made within linear secular theory, and we have derived a simpler theory to study the apsidal state of orbits and to estimate eccen- tricities. Our expression for the mode damping rates, Eq. (8.9), are consistent with that of Wu and Goldreich (2002). In this chapter, we rst use the theory developed in Chapter 8 to analyze the current status of a two-planet system, and then discuss how our equations can be used to predict possible companions of \hot-Jupiters". 138 Table 9.1: Properties of planets and stars discussed in this chapter Planet Orbital Period (days) msini (M J ) a (AU) e ! ( ) m (M ) Age (Gyrs) HIP 14810b 6:674 0:002 3:9 0:6 0:069 0:004 0:147 0:006 159 2 0:99 4 HIP 14810c 95:285 0:002 0:8 0:1 0:41 0:02 0:409 0:006 354 2 0:99 4 GJ 436b 2:64385 0:00009 0:0713 0:006 0:029 0:002 0:16 0:02 351 1 0:44 0:04 > 3 GJ 674b 4:6938 0:007 0.037 0.039 0:20 0:02 143 6 0.35 0.1-1 Data for HIP 14810 system are from Wright et al. (2007); GJ 436 from Maness et al. (2007); GJ 674 from Bon ls et al. (2007). 139 9.1 A Test of the Theory: the HIP 14810 System Several multi-planet systems that host one or more \hot-Jupiters" have been discov- ered, but most are not good candidates upon which to test our theory. The Gliese 876, 55 Cnc, and And systems all have three or more planets, which makes apsidal analysis more complicated (Barnes and Greenberg, 2006a). The HD 217107 system has two planets, but the outer planet is 60 times further away from the star than the inner one, resulting in weak secular interactions that are insu cient to force any in- teresting eccentricity on the inner orbit. The recent discovery of a two-planet system HIP 14810 (Wright et al., 2007), however, provides a good test for our theory. Star HIP 14810a is a Solar-type star with a mass of 0.99 Solar mass. The orbital parame- ters of the two planets are listed in Table 9.1: the 3.9 Jupiter-mass inner planet (HIP 14810b) has an orbital period of 6.674 days, and its orbital eccentricity is as large as 0.147; the outer HIP 14810c is about 20% as massive, and about 6 times further away. It has a larger eccentricity of 0.4. At 0:147 and 0:4, the eccentricities of the two orbits are larger than that assumed by our linear secular theory. But this brings merely a 20% error as shown in Fig. 9.1. Here we depict the orbital histories of the two planets on a phase diagram, similar to the ones shown in Figs. 8.5 and 8.7. We plot curves computed both from secular theory (red) and from N-body simulation (black) and nd reasonable agreement. It appears as if the system librates about $ = 180 , but is near the anti-aligned separatrix so that the libration amplitude is large. This is an example of a general observation that extrasolar multi-planet systems tend to be near a secular separatrix 140 Figure 9.1: The current apsidal state of the HIP 14810 system. We show the orbits of both planets on the eexp(i $) plane, similar to the lower plots in Fig. 8.5 and 8.7. The red curves represent the solution of secular equations, while black curves are obtained through an N-body simulation. The orbits librate about $ = 180 , with the libration width 100 predicted by the secular theory, and 140 measured from N-body simulation. (Barnes and Greenberg, 2006b). The near-separatrix state of the system is surprising; perhaps it is due to insu cient tidal circularization. HIP 14810b is about 0.07 AU from its star, which is relatively far compared to other \hot-Jupiters" so that tidal dissipation may simply too weak to signi cantly alter the apsidal state. How does the current state of the system compare with the expectation from tidal-damping theory as developed in Chapter 8? In Section 8.2, we found that if tides are strong enough, every system will rapidly damp to either apsidally aligned or anti-aligned. Which mode the system nally ends up in depends on the damping rates given by Eqs. (8.4) and (8.11). In a secular eigen-mode, the outer and inner 141 orbits have a predictable eccentricity ratio (Eqs. 8.5 and 8.12). In Fig. 9.2, we show a contour plot of the eccentricity ratio between the two orbits for the more slowly damped eigen-mode in the parameter space (q, ), assuming a xed inner planet b and varying the parameters q and of planet c. The dashed line divides the space into a slow aligned mode region (top left) and a slow anti-aligned region (bottom right). In the blue area, the lifetime of the aligned mode is longer than the age of the systems; in the red region, the anti-aligned mode can last the age of the system; and in the white middle area, both modes should have already damped away. Hence, we expect to nd that the two orbits are apsidally aligned or in aligned libration in the blue area, and they are anti-aligned or in anti-aligned libration in the red region. The boundaries of these areas depend only on the age of the system and the tidal damping timescale. We have assumed an age of 4 billion years and a tidal circularization timescale at 500 million years. An older system age and/or faster damping rate will expand the white area and further constrain the outer planet. The actual location of planet c in the HIP 14810 system is marked in both plots of Fig. 9.2. It resides in the more slowly damped anti-aligned mode region, as expected from Fig. 9.1. A comparison between the two plots in Fig. 9.2 shows the e ect of relativistic pericenter precession. The anti-aligned area (red) diminishes while the aligned area (blue) expands into the region of low-mass distant companions (small q, small , the bottom left corner in the plots), systems in which the relativistic precession rate competes with the planet-induced ones ( _$GR). The shrinking of the red area agrees with Eq. (8.11) which shows that relativistic precession tends to increase the 142 X X a) without general relativity b) with general relativity Figure 9.2: Orbital states of possible companions for HIP 14810b without (a) and with (b) general relativistic precession. We have assumed the tidal-damping timescale to be 500 million years. Each point in these q - plots represents an outer companion with corresponding mass and semi-major axis. The dashed curve divides the plane into regions in which the aligned mode damps faster (top-left) and the anti-aligned mode damps faster (bottom-right). Furthermore, in the colored regions, either the aligned mode (blue) or the anti-aligned mode (red) can survive tidal dissipation and last longer than the age of the system (4 Gyrs). The solid contour lines represent the eccentricity ratio s of the long-lived (slower) mode. The location of HIP 14810c is marked by the black cross. damping rate of the anti-aligned mode and to decrease that of the aligned mode. 9.2 Constraints on Possible Companions of \Hot-Jupiters" If tidal dissipation is strong enough, an extrasolar two-planet system should have evolved into either an aligned or an anti-aligned apsidal-lock state. Equations (8.9) and (8.11) thus provide a constraint on these three parameters: the mass ratio q of the two planets, the semi-major axis ratio , and eccentricity ratio of the two orbits. For a lonely \hot-Jupiter" with con rmed parameters and non-zero eccentricity, we can predict possible companion that may force its orbital eccentricity. We illustrate 143 our method with a few examples. Figure 9.3a is similar to Fig. 9.2, but now for the planetary system around the M dwarf star GJ 436; here we plot contours of the companion?s eccentricity rather than the eccentricity ratios. With only 1.3 Neptune masses, GJ 436b is the rst Neptune- sized extrasolar planet detected through a radial velocity survey (Butler et al., 2004). It has recently been found to transit its host star (Gillon et al., 2007), which removes the uncertainty in the determination of its mass. The planet is 0:023AU away from the star (orbital period 2:6 days), and has an orbital eccentricity of about 0.16. If this eccentricity is due to another planet and the system has damped to an eigen- mode, then Fig. 9.3a excludes almost the entire anti-aligned region and all companions that are less that 1=10 as massive as GJ 436b, or more than 5 times further away { these planets are simply too weak to force the inner planet?s eccentricity. In fact, if we consider the age of the system, which is estimated to be older than 3 billion years (Gillon et al., 2007) and probably much older, most possible outer companions are much more massive than GJ 436b. Larger planets in this region are very likely within the radial-velocity detection window, but none have been found. Most planets outside this region, however, cannot be observationally excluded. Recall that the boundaries between the long-lived mode zones and the short-lives one (dotted lines in Fig. 9.3a) depend on both the tidal timescale and how long the eccentricities have been damped. We have assumed that tides have acted to circularize the orbits throughout the evolution history of the system. It is possible that a later event, such as a resonance crossing, excites the eccentricities and tides have worked for a much 144 Anti?aligned Libration Zone Aligned Libration Zone GJ 436 Anti?aligned Libration Zone Aligned Libration Zone GJ 674 (a) (b) Figure 9.3: Eccentricities of possible companions for the \hot-Neptunes" a) GJ 436b, and b) GJ 674b. These plots are similar to Fig. 9.2b, except that contours show the eccentricity of the outer orbit, and dotted lines replace boundaries between the colored areas; between the dotted lines, both modes should damp away within the 3 Gyr age of the system. shorter time. This pushes the two dotted lines toward the boundary between the faster aligned zone and faster anti-aligned zone (dashed line in Fig. 9.3a), and makes more parameters available to the outer planet. Furthermore, the secular theory is accurate only for small eccentricities. For very eccentric orbits, e.g., e> 0:4, further analysis is necessary to include higher-order terms in the disturbing functions. Another recently discovered \hot-Neptune", GJ 674b in shown in Fig. 9.3b. This system is very similar to the GJ 436 system, except that the star is much younger. With an age of only about half a billion years, most of the parameter space in the q plane is available for a possible secular companion. The slightly larger eccentricity of GJ 674b, however, pushes the eccentricity contours towards the upper-right in Fig. 9.3b, which means that all possible companions with small masses or large semi- major axes are excluded. 145 The assumption of an apsidal lock that we made at the beginning of this section, however, is probably too strong. What we have seen for the HIP 14810 system may be true for most systems: although one mode has been signi cantly damped, its amplitude is not yet near zero. In this case, for a \hot-Jupiter" with known mass, semi-major axis, and eccentricity, the eccentricity of the companion is still free to vary somewhat even for a given mass and semi-major axis. We have shown that the existence of a companion certainly slows the tidal circularization of \hot-Jupiters" and can maintain non-zero orbital eccentricities. 146 Part IV Conclusion 147 Chapter 10 Conclusion and Future Directions Modern computing technology has made it possible to study the orbital motion of the planets over their 4.5 billion year histories, and we are rapidly approaching the point where this will also be possible for satellite systems. Numerical simulation also proves to be an e ective tool to test analytical theories. In this dissertation, we have combined numerical methods with mathematical techniques to address several problems in both resonant and secular orbital dynamics. In the last chapter of this dissertation, we summarize our results, tie up some loose ends, and indicate some promising future directions. 10.1 The Inner Neptunian System In Chapters 3 - 6, we have shown that the dynamics of the small inner Neptunian satellites are strongly a ected by Neptune?s large moon, Triton. We derived mathe- matical tools (Chapter 4) to analyze resonances in this system, and we numerically studied the details of resonance passages between Proteus and other satellites. Tri- ton?s huge mass and large orbital tilt lead to strong three-body resonant kicks in the traditional two-body resonant zones between a pair of small satellites (Chapter 5). We showed that these resonant kicks can explain the existing non-zero orbital tilts of 148 three of the four largest inner moons. If the satellites have a common mean density between 0.4 and 0.8 g=cm3, the inclinations of Proteus, Galatea, and Despina can be excited by three or four Proteus resonance passages. The tilt of Larissa, however, is twice as large as can be explained with these resonant events alone. In Chapter 6, we suggested that a resonant capture event might have tilted Larissa?s orbit, and that this event might also be responsible for the large inclination of Naiad (Table 3.1). 10.1.1 The 4:7 Inclination of Naiad The origin of Naiad?s large tilt was rst investigated by Ban eld and Murray (1992), who undertook a resonant-history analysis for the inner Neptunian satellites similar to our study in Section 3.2. Since the importance of three-body resonances was not known at the time, the authors suggested that Naiad might have been trapped into one of the standard two-body resonances with another satellite. Because two-body resonances are so weak (cf. Fig. 5.2), the probability of resonant capture is low, typically only a few percent. The number of possible resonance passages, however, is large and the total capture probability could reach 76% with multiple resonance passages invoked (Ban eld and Murray, 1992). With our discovery of the strong three-body resonances with Triton, resonance trapping is much more probable. Our numerical observations show that the capture probability for a zero-degree-inclined Naiad approaches one hundred percent. In other words, Naiad would likely be trapped into the rst three-body resonance it passes across, as shown in Fig. 10.1. This gure shows that Naiad and Galatea are 149 Figure 10.1: A numerical simulation in which Naiad and Galatea are trapped into their 5:3 mean-motion resonance with Triton. The corresponding resonant angle iNiT = 5 G 3 N ~ N T is shown in the bottom panel. trapped into a 5:3 inclination resonance with Triton when Galatea migrates inward and the two orbits converge. During the resonant trapping, the free inclination of Naiad grows in time, and its semi-major axis decreases with Galatea?s to keep the satellites in resonance. In order to explain Naiad?s current tilt, the satellites must exit the resonance when ifrN 4:7 . Multiple excitations are unlikely because the capture probabil- ity decreases dramatically once the satellite?s inclination is above a critical value (Borderies and Goldreich, 1984). Several mechanisms can break the resonance lock, including additional mean-motion resonances, secondary resonances (resonances be- tween the libration frequencies of nearby mean-motion resonances, Malhotra, 1990), 150 and possible non-resonant perturbations from other satellites. A systematic study is necessary to understand which resonance, with which satellite, might have been responsible for the capture, and to determine how resonance evolves after capture. Besides providing an explanation for the large inclinations of Larissa and Naiad, a successful theory would also further constrain tidal models. 10.1.2 Con nement of Neptunian Ring Arcs The puzzles in the inner Neptunian system arose even before the discovery of the six small satellites. After the rst evidence for the existence of Larissa (Reitsema et al., 1982), several observations of stellar occultations (e.g. Hubbard et al., 1986) implied possible narrow ring arcs around Neptune, which were later con rmed during the Voyager 2 y-by (Smith et al., 1989). Goldreich et al. (1986) suggested that the sharp edge of the rings could be con ned by a Lindblad resonance ( = (p+1) 2 p 1 $1) with a hypothetical satellite, and further proposed that the ring arcs were azimuthally con ned by an inclination co-rotation resonance ( = (p + 2) 2 p 1 2 1). This shepherd satellite turns out to be Galatea, whose 43:42 resonance multiplet sits near the ring arcs (Porco, 1991). More recent observations (Sicardy et al., 1999; Dumas et al., 1999), however, showed that the rings are not centered on the inclination co- orbital resonance. This led Namouni and Porco (2002) to seek a new explanation in which the non-zero mass of the ring arcs moves the location of the resonances. These authors invoked an eccentricity co-rotation resonance to con ne the ring arcs and estimated that the total mass of the ring arcs is 0:002 0:2 times the mass of 151 Galatea, which is large for such a narrow ring system. The requirement for a massive ring might be relaxed if Triton?s e ect were to be considered. Although none of the strong Triton resonances are located especially close to the ring arcs, and the standard two-body resonances are not shifted signi cantly, the system is complicated and warrants a more careful investigation. 10.2 Perturbed Secular Interactions Additional perturbations that a ect planetary orbits, such as mass accretion and ra- dial migration, may interfere with secular interactions. In Chapter 7 - 9, we examined the case of slow eccentricity decay due to tidal drag. Because of secular interactions, eccentricity damping imposed on one orbit can be transmitted to the other, and the timescale for total damping can be substantially slowed. The current eccentricities of \hot-Jupiters" might be secularly forced by unknown companions. Here we expand our attention to consider the e ects of additional perturbations on the system. 10.2.1 E ects of Adiabatic Changes in a and m Although tidal circularization is important for \hot-Jupiters", additional perturba- tions on orbits may be important in other situations. Here we consider both planet migration and mass loss or accretion. For the simplest case in which the mass of the outer planet is much less than that of the inner planet, or q p , Eqs. (7.11) and (7.12) can be simpli ed to: gs+ = 0;gs = p ; s+ = ; s = 1; 152 and the general solution of the two orbits are (Eq. 7.16) h1 = e+; h2 = e+ +e exp[i( p t+? )]: Thus, the inner orbit is not a ected, and the eccentricity of the outer orbit consists of a \forced" component eforced = e+ and a \free" component efr = e , similar to the discussion in Chapter 4 where Triton forces the inclinations of small Neptunian satellites. Note that for a small free eccentricity, these are aligned orbits. Similar arguments for a system with an inner planet of negligible mass (q p ) also lead to a near aligned con guration: h1 = e+ +e exp[i( qt+? )]; h2 = 1 e+: For slow adiabatic changes of any of a1, a2, m1, and m2, we nd that the free eccentricity is an adiabatic invariant of the system, so it remains nearly constant over time. For more general cases where the two masses are comparable, we expect that similar adiabatic invariants exist, but to date we have been unable to prove this assertion. If these more general adiabatic invariants can be found, we would be able to determine the response of a system to a broad class of slow perturbations. 10.2.2 Systems with Three or More Planets Approximately 10 planetary systems with 3 or more planets have been discovered to date, which motivates us to extend our approach to N planets. Since the system 153 in Eq. (7.9) is linear, a solution exists for any N. The solution to the characteristic equation for the eigenmodes requires solving a polynomial of order N, which is best done numerically. The identi cation of the eigenmodes, however, can be explored qualitatively. Recall that for a two-planet system, there are two natural apsidal co- precession states (aligned and anti-aligned as in Fig. 7.2) and they comprise the two eigenmodes. With three planets, it can be shown that all eigenmodes consist of aligned and/or anti-aligned pericenters. There are four possible states (Figure. 10.2) while the system can have only three eigenmodes. We have found that for di erent systems, sets of three eigenmodes are selected from the four natural states: the fully aligned and fully anti-aligned states are always eigen-modes, while the remaining one depends on mass ratios and orbital spacings. For N (> 3) planets, it is even more complicated to determine N eigen-modes from 2N 1 candidates, but is interesting to pursue. The application of this study would Figure 10.2: The natural apsidal co- precession states for a three-planet system. Arrows show the direction from the star to the orbital pericenters for each planet, with the inner-most planet on the top. The left- most (fully aligned) and rightmost (fully anti-aligned) sets are always eigenmodes. not be limited to extrasolar planetary systems, but could also be used to fur- ther understand the dynamics of our own Solar System. One long-standing puzzle is the origin of the inclinations and ec- centricities of the giant planets. Much attention has been paid to the large eccentric- ities seen in other planetary systems, leading to speculation that large eccentricities might be the norm and that our Solar System is an exception. On the other hand, 154 the simplest disk formation scenarios suggest that planetary orbital planes should be precisely parallel to one another (as the case of inner Neptunian satellites discussed in Part II), and planets should all orbit along perfect circles. This may indeed be the primitive state, with all non-zero eccentricities and inclinations generated by subse- quent planetary interactions (Tsiganis et al., 2005) and/or interactions with the gas and planetesimal disks. 155 References Adams, F. C., Laughlin, G., 2006. E ects of secular interactions in extrasolar plane- tary systems. ApJ 649, 992{1003. Adams, J. C., 1846. Explanation of the observed irregularities in the motion of Uranus, on the hypothesis of disturbance by a more distant planet. MNRAS 7, 149{152. Agnor, C. B., Canup, R. M., Levison, H. F., 1999. On the character and consequences of large impacts in the late stage of terrestrial planet formation. Icarus 142, 219{ 237. Agnor, C. B., Hamilton, D. P., 2006. Neptune?s capture of its moon Triton in a binary-planet gravitational encounter. Nature 441, 192{194. Anderson, J. D., Johnson, T. V., Schubert, G., Asmar, S., Jacobson, R. A., Johnston, D., Lau, E. L., Lewis, G., Moore, W. B., Taylor, A., Thomas, P. C., Weinwurm, G., 2005. Amalthea?s density is less than that of water. Science 308, 1291{1293. Aumann, H. H., 1985. IRAS observations of matter around nearby stars. PASP 97, 885{891. Ban eld, D., Murray, N., 1992. A dynamical history of the inner Neptunian satellites. Icarus 99, 390{401. Barnes, R., Greenberg, R., 2006a. Behavior of apsidal orientations in planetary sys- tems. ApJ 652, L53{L56. Barnes, R., Greenberg, R., 2006b. Extrasolar planetary systems near a secular sepa- ratrix. ApJ 638, 478{487. Beaug e, C., Ferraz-Mello, S., Michtchenko, T. A., 2003. Extrasolar planets in mean- motion resonance: Apses alignment and asymmetric stationary solutions. ApJ 593, 1124{1133. Bon ls, X., Mayor, M., Delfosse, X., Forveille, T., Gillon, M., Perrier, C., Udry, S., Bouchy, F., Lovis, C., Pepe, F., Queloz, D., Santos, N. C., Bertaux, J. ., 2007. The HARPS search for southern extra-solar planets. X. a msinI = 11ME planet around the nearby spotted M dwarf GJ 674. ArXiv e-prints, astro{ph/0704.0270. Boquet, F., 1889. Developpement de la fonction perturbatrice : calcul des termes DU huitie ME ordre. Annales de l?Observatoire de Paris 19, B1{B75. Borderies, N., Goldreich, P., 1984. A simple derivation of capture probabilities for the J+1:J and J+2:J orbit-orbit resonance problems. Celestial Mechanics 32, 127{136. Bottke, W. F., Jedicke, R., Morbidelli, A., Petit, J.-M., Gladman, B., 2000. Under- standing the distribution of near-Earth asteroids. Science 288, 2190{2194. 156 Bottke, W. F., Morbidelli, A., Jedicke, R., Petit, J.-M., Levison, H. F., Michel, P., Metcalfe, T. S., 2002. Debiased orbital and absolute magnitude distribution of the near-Earth objects. Icarus 156, 399{433. Brouwer, D., Clemence, G. M., 1961. Methods of celestial mechanics. New York: Academic Press. Brouwer, D., van Woerkom, A., Jasper, J., 1950. The secular variations of the orbital elements of the principal planets. Washington, U. S. Govt. Print. O . Burns, J. A., 1977. Orbital evolution. In: Burns, J. A. (Ed.), Planetary Satellites. Univ. of Arizona Press, Tuscon, AZ, USA, pp. 113{156. Burns, J. A., 1986. The evolution of satellite orbits. In: Burns, J. A., Matthews, M. S. (Eds.), Satellites. Univ. of Arizona Press, Tuscon, AZ, USA, pp. 117{158. Burns, J. A., Showalter, M. R., Mor ll, G. E., 1984. The ethereal rings of Jupiter and Saturn. In: Greenberg, R., Brahic, A. (Eds.), Planetary Rings. Univ. of Arizona Press, Tuscon, AZ, USA, pp. 200{272. Bur sa, M., 1992. Secular Love numbers of the major planets. Earth, Moon, and Planets 59, 239{244. Butler, R. P., Marcy, G. W., Fischer, D. A., Brown, T. M., Contos, A. R., Korzennik, S. G., Nisenson, P., Noyes, R. W., 1999. Evidence for multiple companions to Andromedae. ApJ 526, 916{927. Butler, R. P., Vogt, S. S., Marcy, G. W., Fischer, D. A., Wright, J. T., Henry, G. W., Laughlin, G., Lissauer, J. J., 2004. A Neptune-mass planet orbiting the nearby M dwarf GJ 436. ApJ 617, 580{588. Chambers, J. E., 1999. A hybrid symplectic integrator that permits close encounters between massive bodies. MNRAS 304, 793{799. Chambers, J. E., Wetherill, G. W., 1998. Making the terrestrial planets: N-body integrations of planetary embryos in three dimensions. Icarus 136, 304{327. Chambers, J. E., Wetherill, G. W., 2001. Planets in the Asteroid Belt. Meteoritics and Planetary Science 36, 381{399. Chiang, E. I., 2003. Excitation of orbital eccentricities by repeated resonance cross- ings: Requirements. ApJ 584, 465{471. Cohen, C. J., Hubbard, E. C., 1965. Libration of the close approaches of Pluto to Neptune. AJ 70, 10{13. Cuk, M., Gladman, B. J., 2005. Constraints on the orbital evolution of Triton. ApJ 626, L113{L116. 157 Cuzzi, J. N., Lissauer, J. J., Esposito, L. W., Holberg, J. B., Marouf, E. A., Tyler, G. L., Boishchot, A., 1984. Saturn?s rings - properties and processes. In: Greenberg, R., Brahic, A. (Eds.), Planetary Rings. Univ. of Arizona Press, Tuscon, AZ, USA, pp. 73{199. Danby, J. M. A., 1988. Fundamentals of Celestial Mechanics (2nd ed.). Willmann-Bell Inc., P.O.Box 35025, Richmond, VA, USA. Darwin, G. H., 1879. A tidal theory of the evolution of satellites. The Observatory 3, 79{84. Darwin, G. H., 1880. On the secular e ects of tidal friction. Astronomische Nachrichten 96, 217{222. Dermott, S. F., Malhotra, R., Murray, C. D., 1988. Dynamics of the Uranian and Saturnian satellite systems - a chaotic route to melting Miranda? Icarus 76, 295{ 334. Dermott, S. F., Nicholson, P. D., 1986. Masses of the satellites of Uranus. Nature 319, 115{120. Dobrovolskis, A. R., 1993. The Laplace planes of Uranus and Pluto. Icarus 105, 400{407. Dumas, C., Terrile, R. J., Smith, B. A., Schneider, G., Becklin, E. E., 1999. Stability of Neptune?s ring arcs in question. Nature 400, 733{735. Duncan, M. J., Levison, H. F., 1997. A scattered comet disk and the origin of Jupiter family comets. Science 276, 1670{1672. Duncan, M. J., Levison, H. F., Budd, S. M., 1995. The dynamical structure of the Kuiper Belt. AJ 110, 3073{3081. Duncan, M. J., Levison, H. F., Lee, M. H., 1998. A multiple time step symplectic algorithm for integrating close encounters. AJ 116, 2067{2077. Elliot, J. L., Dunham, E., Mink, D., 1977. The rings of Uranus. Nature 267, 328{330. Ellis, K. M., Murray, C. D., 2000. The disturbing function in Solar System dynamics. Icarus 147, 129{144. Esposito, L., 2006. Planetary Rings. Cambridge University Press, The Edinburgh Building, Cambridge CB2 2RU, UK. Fernandez, J. A., Ip, W.-H., 1984. Some dynamical aspects of the accretion of Uranus and Neptune - the exchange of orbital angular momentum with planetesimals. Icarus 58, 109{120. 158 Ford, E. B., Havlickova, M., Rasio, F. A., 2001. Dynamical instabilities in extrasolar planetary systems containing two giant planets. Icarus 150, 303{313. Froeschle, C., Scholl, H., 1987. Orbital evolution of asteroids near the secular reso- nance 6. A&A 179, 294{303. Galle, J. G., 1846. Account of the discovery of Le Verrier?s planet Neptune, at Berlin. MNRAS 7, 153{163. Gavrilov, S. V., Zharkov, V. N., 1977. Love numbers of the giant planets. Icarus 32, 443{449. Gillon, M., Pont, F., Demory, B. ., Mallmann, F., Mayor, M., Mazeh, T., Queloz, D., Shporer, A., Udry, S., Vuissoz, C., 2007. Detection of transits of the nearby hot Neptune GJ 436 b. ArXiv e-prints, astro{ph/0705.2219. Goldreich, P., Murray, N., Longaretti, P. Y., Ban eld, D., 1989. Neptune?s story. Science 245, 500{504. Goldreich, P., Sari, R., 2003. Eccentricity evolution for planets in gaseous disks. ApJ 585, 1024{1037. Goldreich, P., Soter, S., 1966. Q in the Solar System. Icarus 5, 375{389. Goldreich, P., Tremaine, S., Borderies, N., 1986. Towards a theory for Neptune?s arc rings. AJ 92, 490{494. Goldreich, R., 1963. On the eccentricity of satellite orbits in the solar system. MN- RAS 126, 257{268. Gomes, R., Levison, H. F., Tsiganis, K., Morbidelli, A., 2005. Origin of the cata- clysmic Late Heavy Bombardment period of the terrestrial planets. Nature 435, 466{469. Greaves, J. S., Holland, W. S., Wyatt, M. C., Dent, W. R. F., Robson, E. I., Coulson, I. M., Jenness, T., Moriarty-Schieven, G. H., Davis, G. R., Butner, H. M., Gear, W. K., Dominik, C., Walker, H. J., 2005. Structure in the Eridani debris disk. ApJ 619, L187{L190. Greenberg, R., 1977. Orbit-orbit resonances among natural satellites. In: Burns, J. A. (Ed.), Planetary Satellites. Univ. of Arizona Press, Tuscon, AZ, USA, pp. 157{168. Greenberg, R., 1981. Apsidal precession of orbits about an oblate planet. AJ 86, 912{914. Greenberg, R. J., Counselman, C. C., Shapiro, I. I., 1972. Orbit-orbit resonance capture in the solar system. Science 178, 747{749. 159 Hamilton, D. P., 1994. A comparison of Lorentz, planetary gravitational, and satellite gravitational resonances. Icarus 109, 221{240. Hamilton, D. P., 1996. The asymmetric time-variable rings of Mars. Icarus 119, 153{172. Hamilton, D. P., Burns, J. A., 1993. Lorentz and gravitational resonances on circum- planetary particles. Adv. Space Res. 13 (10), 241{248. Hamilton, D. P., Zhang, K., Agnor, C., 2005. Constraints on Triton?s orbital evolu- tion. AAS/Division of Dynamical Astronomy Meeting 36, 11.04. Henrard, J., 1982. Capture into resonance - an extension of the use of adiabatic invariants. Celestial Mechanics 27, 3{22. Holland, W. S., Greaves, J. S., Zuckerman, B., Webb, R. A., McCarthy, C., Coulson, I. M., Walther, D. M., Dent, W. R. F., Gear, W. K., Robson, I., 1998. Submillimetre images of dusty debris around nearby stars. Nature 392, 788{790. Holman, M., Touma, J., Tremaine, S., 1997. Chaotic variations in the eccentricity of the planet orbiting 16 Cygni B. Nature 386, 254{256. Holman, M. J., Wisdom, J., 1993. Dynamical stability in the outer Solar System and the delivery of short period comets. AJ 105, 1987{1999. Hubbard, W. B., Brahic, A., Sicardy, B., Elicer, L.-R., Roques, F., Vilas, F., 1986. Occultation detection of a Neptunian ring-like arc. Nature 319, 636{640. Jacobson, R. A., Owen, W. M., 2004. The orbits of the inner Neptunian satellites from Voyager, Earth-based, and Hubble Space Telescope observations. AJ 128, 1412{1417. Jacobson, R. A., Riedel, J. E., Taylor, A. H., 1991. The orbits of Triton and Nereid from spacecraft and earthbased observations. A&A 247, 565{575. Je reys, H., 1961. The e ect of tidal friction on eccentricity and inclination. MN- RAS 122, 339{343. Juric, M., Tremaine, S., 2007. The eccentricity distribution of extrasolar planets. ArXiv Astrophysics e-prints, astro{ph/0703160. Karkoschka, E., 2003. Sizes, shapes, and albedos of the inner satellites of Neptune. Icarus 162, 400{407. Kaula, W. M., 1964. Tidal dissipation by solid friction and the resulting orbital evolution. Rev. Geophys. 2, 661{685. Kuchner, M. J., Holman, M. J., 2003. The geometry of resonant signatures in debris disks with planets. ApJ 588, 1110{1120. 160 Landau, L. D., Lifshitz, E. M., 1976. Course of Theoretical Physics, Volume 1: Mechanics (3rd ed.). Pergamon Press, Elmsford, NY, USA. Le Verrier, M., 1846. Sur la recti cation des orbites des com etes au moyen de l?ensemble des observations faites pendant leur apparition. Astronomische Nachrichten 23, 183{188. Le Verrier, U.-J., 1855. Recherches astronomiques: Chapitre IV. - d eveloppement de la fonction qui sert de base au calcul des perturbations des mouvements des plan etes. Annales de l?Observatoire de Paris 1, 258{342. Lee, M. H., 2004. Diversity and origin of 2:1 orbital resonances in extrasolar planetary systems. ApJ 611, 517{527. Levison, H. F., Duncan, M. J., 1994. The long-term dynamical behavior of short- period comets. Icarus 108, 18{36. Levison, H. F., Duncan, M. J., 1997. From the Kuiper Belt to Jupiter-family comets: The spatial distribution of ecliptic comets. Icarus 127, 13{32. MacDonald, G. J. F., 1964. Tidal friction. Rev. Geophys. 2, 467{541. Malhotra, R., 1990. Capture probabilities for secondary resonances. Icarus 87, 249{ 264. Malhotra, R., 1993. The origin of Pluto?s peculiar orbit. Nature 365, 819{821. Malhotra, R., 1995. The origin of Pluto?s orbit: Implications for the Solar System beyond Neptune. AJ 110, 420{429. Malhotra, R., Dermott, S. F., 1990. The role of secondary resonances in the orbital history of Miranda. Icarus 85, 444{480. Maness, H. L., Marcy, G. W., Ford, E. B., Hauschildt, P. H., Shreve, A. T., Basri, G. B., Butler, R. P., Vogt, S. S., 2007. The M dwarf GJ 436 and its Neptune-mass planet. PASP 119, 90{101. Marzari, F., Scholl, H., 2000. The role of secular resonances in the history of Trojans. Icarus 146, 232{239. Mayor, M., Queloz, D., 1995. A Jupiter-mass companion to a Solar-type star. Na- ture 378, 355{359. McKinnon, W. B., 1984. On the origin of Triton and Pluto. Nature 311, 355{358. Michel, P., Migliorini, F., Morbidelli, A., Zappal a, V., 2000. The population of Mars- crossers: Classi cation and dynamical evolution. Icarus 145, 332{347. 161 Milani, A., Carpino, M., Hahn, G., Nobili, A. M., 1989. Dynamics of planet-crossing asteroids - classes of orbital behavior. Icarus 78, 212{269. Miner, E. D., Wessen, R. R., Cuzzi, J. N., 2007. Springer and Praxis Publishing, Chichester, UK. Morbidelli, A., 1997. Chaotic di usion and the origin of comets from the 2/3 resonance in the Kuiper Belt. Icarus 127, 1{12. Morbidelli, A., 2002. Modern integrations of Solar System dynamics. Annual Review of Earth and Planetary Sciences 30, 89{112. Morbidelli, A., Gladman, B., 1998. Orbital and temporal distributions of meteorites originating in the Asteroid Belt. Meteoritics and Planetary Science 33, 999{1016. Morbidelli, A., Levison, H. F., Gomes, R., 2007. The dynamical structure of the Kuiper Belt and its primordial origin. ArXiv Astrophysics e-prints, astro{ ph/0703558. Morbidelli, A., Levison, H. F., Tsiganis, K., Gomes, R., 2005. Chaotic capture of Jupiter?s Trojan asteroids in the early Solar System. Nature 435, 462{465. Murray, C. D., 1986. Structure of the 2:1 and 3:2 Jovian resonances. Icarus 65, 70{82. Murray, C. D., Dermott, S. F., 1999. Solar System Dynamics. Cambridge University Press, The Edinburgh Building, Cambridge CB2 2RU, UK. Murray, C. D., Harper, D., 1993. Expansion of the planetary disturbing function to eighth order in the individual orbital elements. QMW Math Notes No. 15, Queen Mary and West eld College, Mile End Road, London E1 4NS, U.K. Murray, N., Holman, M., 1999. The origin of chaos in the outer Solar System. Science 283, 1877{1881. Nagasawa, M., Ida, S., 2000. Sweeping secular resonances in the Kuiper Belt caused by depletion of the solar nebula. AJ 120, 3311{3322. Namouni, F., 2007. Origin theories for the eccentricities of extrasolar planets. ArXiv Astrophysics e-prints, astro{ph/0702203. Namouni, F., Porco, C., 2002. The con nement of Neptune?s ring arcs by the moon Galatea. Nature 417, 45{47. Nesvorn y, D., Ferraz-Mello, S., 1997. On the asteroidal population of the rst-order Jovian resonances. Icarus 130, 247{258. Nesvorn y, D., Roig, F., 2001. Mean motion resonances in the trans-Neptunian region part II: The 1:2, 3:4, and weaker resonances. Icarus 150, 104{123. 162 Nicholson, P. D., Hamilton, D. P., Matthews, K., Yoder, C. F., 1992. New observa- tions of Saturn?s coorbital satellites. Icarus 100, 464{484. Owen, T., Danielson, G. E., Cook, A. F., Hansen, C., Hall, V. L., Duxbury, T. C., 1979. Jupiter?s rings. Nature 281, 442{446. Owen, W. M., Vaughan, R. M., Synnott, S. P., 1991. Orbits of the six new satellites of Neptune. AJ 101, 1511{1515. Peale, S. J., 1976. Orbital resonances in the solar system. ARA&A 14, 215{246. Peale, S. J., 1986. Orbital resonances, unusual con gurations and exotic rotation states among planetary satellites. In: Burns, J. A., Matthews, M. S. (Eds.), Satel- lites. Univ. of Arizona Press, Tuscon, AZ, USA, pp. 159{223. Peale, S. J., Cassen, P., Reynolds, R. T., 1980. Tidal dissipation, orbital evolution, and the nature of Saturn?s inner satellites. Icarus 43, 65{72. Peirce, B., 1849. Development of the perturbative function of planetary motion. AJ 1, 1{8. Petit, J.-M., Morbidelli, A., Chambers, J., 2001. The primordial excitation and clearing of the Asteroid Belt. Icarus 153, 338{347. Porco, C. C., 1991. An explanation for Neptune?s ring arcs. Science 253, 995{1001. Porco, C. C., Baker, E., Barbara, J., Beurle, K., Brahic, A., Burns, J. A., Charnoz, S., Cooper, N., Dawson, D. D., Del Genio, A. D., Denk, T., Dones, L., Dyudina, U., Evans, M. W., Giese, B., Grazier, K., Helfenstein, P., Ingersoll, A. P., Jacobson, R. A., Johnson, T. V., McEwen, A., Murray, C. D., Neukum, G., Owen, W. M., Perry, J., Roatsch, T., Spitale, J., Squyres, S., Thomas, P., Tiscareno, M., Turtle, E., Vasavada, A. R., Veverka, J., Wagner, R., West, R., 2005. Cassini imaging science: Initial results on Saturn?s rings and small satellites. Science 307, 1226{ 1236. Porco, C. C., Hamilton, D. P., 2007. Planetary rings. In: McFadden, L. A., Weissman, P. R., Johnson, T. V. (Eds.), Encyclopedia of the Solar System (2nd ed.). Academic Press, San Diego, US, pp. 503{518. Rauch, K. P., Hamilton, D. P., 2002. The HNBody package for symplectic integration of nearly-Keplerian systems. Bulletin of the American Astronomical Society 34, 938. Reitsema, H. J., Hubbard, W. B., Lebofsky, L. A., Tholen, D. J., 1982. Occultation by a possible third satellite of Neptune. Science 215, 289{291. Renner, S., Sicardy, B., French, R. G., 2005. Prometheus and Pandora: masses and orbital positions during the Cassini tour. Icarus 174, 230{240. 163 Roy, A. E., 2005. Orbital Motion (4th ed.). Institute of Physics Publishing, Dirac House, Temple Back, Bristol BS1 6BE, UK. Scholl, H., Froeschle, C., 1990. Orbital evolution of known asteroids in the 5 secular resonance region. A&A 227, 255{263. Sheppard, S. S., Trujillo, C. A., 2006. A thick cloud of Neptune Trojans and their colors. Science 313, 511{514. Showalter, M. R., Burns, J. A., 1982. A numerical study of Saturn?s F-ring. Icarus 52, 526{544. Sicardy, B., Roddier, F., Roddier, C., Perozzi, E., Graves, J. E., Guyon, O., North- cott, M. J., 1999. Images of Neptune?s ring arcs obtained by a ground-based tele- scope. Nature 400, 731{733. Smith, B. A., Soderblom, L. A., Ban eld, D., Barnet, C., Basilevsky, A. T., Beebe, R. F., Bollinger, K., Boyce, J. M., Brahic, A., Briggs, G. A., Brown, R. H., Chyba, C., Collins, S. A., Colvin, T., Cook, A. F., Crisp, D., Croft, S. K., Cruikshank, D., Cuzzi, J. N., Danielson, G. E., Davies, M. E., De Jong, E., Dones, L., Godfrey, D., Goguen, J., Grenier, I., Haemmerle, V. R., Hammel, H., Hansen, C. J., Helfenstein, C. P., Howell, C., Hunt, G. E., Ingersoll, A. P., Johnson, T. V., Kargel, J., Kirk, R., Kuehn, D. I., Limaye, S., Masursky, H., McEwen, A., Morrison, D., Owen, T., Owen, W., Pollack, J. B., Porco, C. C., Rages, K., Rogers, P., Rudy, D., Sagan, C., Schwartz, J., Shoemaker, E. M., Showalter, M., Sicardy, B., Simonelli, D., Spencer, J., Sromovsky, L. A., Stoker, C., Strom, R. G., Suomi, V. E., Synott, S. P., Terrile, R. J., Thomas, P., Thompson, W. R., Verbiscer, A., Veverka, J., 1989. Voyager 2 at Neptune - Imaging science results. Science 246, 1422{1449. Soter, S., 1971. The dust belts of Mars. Cornell CRSR Report 472. Stockwell, J. N., 1873. Memoir on the secular variations of the elements of the oribts of the eight principal planets. Smithsonian Contributions to Knowledge 18, 1{199. Sussman, G. J., Wisdom, J., 1988. Numerical evidence that the motion of Pluto is chaotic. Science 241, 433{437. Sussman, G. J., Wisdom, J., 1992. Chaotic evolution of the Solar System. Sci- ence 257, 56{62. Tittemore, W. C., Wisdom, J., 1988. Tidal evolution of the Uranian satellites. I - passage of Ariel and Umbriel through the 5:3 mean-motion commensurability. Icarus 74, 172{230. Tittemore, W. C., Wisdom, J., 1989. Tidal evolution of the Uranian satellites. II - an explanation of the anomalously high orbital inclination of Miranda. Icarus 78, 63{89. 164 Tsiganis, K., Gomes, R., Morbidelli, A., Levison, H. F., 2005. Origin of the orbital architecture of the giant planets of the Solar System. Nature 435, 459{461. Ward, W. R., 1988. On disk-planet interactions and orbital eccentricities. Icarus 73, 330{348. Ward, W. R., Colombo, G., Franklin, F. A., 1976. Secular resonance, solar spin down, and the orbit of Mercury. Icarus 28, 441{452. Wilner, D. J., Holman, M. J., Kuchner, M. J., Ho, P. T. P., 2002. Structure in the dusty debris around Vega. ApJ 569, L115{L119. Wisdom, J., 1980. The resonance overlap criterion and the onset of stochastic behavior in the restricted three-body problem. AJ 85, 1122{1133. Wisdom, J., 1983. Chaotic behavior and the origin of the 3/1 Kirkwood gap. Icarus 56, 51{74. Wisdom, J., 1987. Urey Prize Lecture - Chaotic dynamics in the Solar System. Icarus 72, 241{275. Wisdom, J., Holman, M., 1991. Symplectic maps for the N-body problem. AJ 102, 1528{1538. Wolszczan, A., Frail, D. A., 1992. A planetary system around the millisecond pulsar PSR1257 + 12. Nature 355, 145{147. Wright, J. T., Marcy, G. W., Fischer, D. A., Butler, R. P., Vogt, S. S., Tinney, C. G., Jones, H. R. A., Carter, B. D., Johnson, J. A., McCarthy, C., Apps, K., 2007. Four new exoplanets and hints of additional substellar companions to exoplanet host stars. ApJ 657, 533{545. Wu, Y., Goldreich, P., 2002. Tidal evolution of the planetary system around HD 83443. ApJ 564, 1024{1027. Yoder, C. F., 1973. On the Establishment and Evolution of Orbit-Orbit Resonances. PhD Thesis. Yoder, C. F., 1995. Astrometric and geodetic properties of Earth and the Solar System. In: Ahrens, T. (Ed.), Global Earth Physics: A Handbook of Physical Constants. American Geophysical Union, Washington, pp. 1{31. Yoder, C. F., Peale, S. J., 1981. The tides of Io. Icarus 47, 1{35. Zakamska, N. L., Tremaine, S., 2004. Excitation and propagation of eccentricity disturbances in planetary systems. AJ 128, 869{877. Zhang, K., Hamilton, D. P., 2007. Orbital resonances in the inner Neptunian system. Icarus 188, 386{399. 165