# The GENGA Code: Gravitational Encounters in N-body simulations with GPU Acceleration

DOI: 10.1088/0004-637X/796/1/23 · Source: arXiv

Abstract

We describe a GPU implementation of a hybrid symplectic N-body integrator,
GENGA (Gravitational ENcounters with Gpu Acceleration), designed to integrate
planet and planetesimal dynamics in the late stage of planet formation and
stability analysis of planetary systems. GENGA is based on the integration
scheme of the Mercury code (Chambers 1999), which handles close encounters with
very good energy conservation. It uses mixed variable integration (Wisdom &
Holman 1991) when the motion is a perturbed Kepler orbit and combines this with
a direct N-body Bulirsch-Stoer method during close encounters. The GENGA code
supports three simulation modes: Integration of up to 2048 massive bodies,
integration with up to a million test particles, or parallel integration of a
large number of individual planetary systems. GENGA is written in CUDA C and
runs on all Nvidia GPUs with compute capability of at least 2.0. All operations
are performed in parallel, including the close encounter detection and the
grouping of independent close encounter pairs. Compared to Mercury, GENGA runs
up to 30 times faster. GENGA is available as open source code from
https://bitbucket.org/sigrimm/genga.

- ... Hamada et al. [16] came up with an efficient GPU implementation of Barnes–Hut algorithm on large data sets in astrophysics and turbulence and Jetley et al. [18] provided a scalable implementation on GPU clusters. Grimm and Stadel [15] designed a hybrid symplectic N-body integrator to analysis planet and planetesimal dynamics in the late stage of planet formation with GPU acceleration. Significant speedups have been observed in the above applications. ...
- ... Third, we use Monte Carlo impact simulations to determine the amount of mass accreted by each satellite. All dynamical simulations were performed using the Gravitational ENcounters with Gpu Acceleration (GENGA) (Grimm & Stadel, 2014). The Monte Carlo (MC) simulations are based on the studies by Brasser et al. (2016) and Brasser & Mojzsis (2017 ...Preprint
- Nov 2018

The intensity and effects of early impact bombardment on the major satellites of the giant planets during an episode of giant planet migration is still poorly known. We use a combination of dynamical N-body and Monte Carlo simulations to determine impact probabilities, impact velocities, and expected masses that collide with these satellites to determine the chronology of impacts during the migration. Volatile loss through bombardment is typically 20% for Miranda, a few percents for the larger Uranian satellites and negligible for the Galilean satellites. Due to its small size and the high impact velocity there is a >99% chance that Miranda suffered a catastrophic impact that shattered the satellite. Subsequent re-accretion from a circum-Uranian ring could account for its peculiar surface morphology and low density. The probability to destroy Ariel and Umbriel is 15% and 1% for Titania and Oberon. Approximately 90% of the mass in planetesimals that passes through the Jovian and Uranian satellite systems (about $4 {\rm \ M_{\oplus}}$ and $2 {\rm \ M_{\oplus}}$ respectively) does so in about 15 Myr. This extremely rapid and intense bombardment causes repeated local crustal melting on all satellites. The combination of these effects results in an entirely different impact chronology than that of the inner solar system. We conclude that the simple extrapolation of the lunar chronology to the outer solar system satellites is not correct. The tail end (after 25 Myr) of the chronology function has an e-folding time of 100 Myr at Jupiter, but follows a cumulative Weibull distribution at Uranus, making direct comparisons between the gas and ice giant planets difficult. Based on our results the surfaces of the Uranian satellites, Callisto, and possibly Ganymede, are all about the same age, and are roughly 150 Myr younger than the timing of the dynamical instability. - ... Nbody simulations are accurate and can automatically handle any complicated phenomena, such as resonant interactions and spatially non-uniform distributions of planetesimals. Gravity calculations are accelerated by such as tree-methods (Richardson et al., 2000) or special hardwares (Kokubo and Ida, 1996; Grimm and Stadel, 2014), and a large time step can be used with sophisticated integrators, such as Mixed Variable Symplectic (MVS) or Wisdom-Holman integrators (Duncan et al., 1998; Chambers, 1999). Even with these novel techniques, N-body simulations are computationally intense and the number of particles or the number of time steps in a simulation is severely limited. ...Article
- Oct 2015
- Icarus

We introduce a new particle-based hybrid code for planetary accretion. The code uses an N-body routine for interactions with planetary embryos while it can handle a large number of planetesimals using a super-particle approximation, in which a large number of small planetesimals are represented by a small number of tracers. Tracer-tracer interactions are handled by a statistical routine which uses the phase-averaged stirring and collision rates. We compare hybrid simulations with analytic predictions and pure N-body simulations for various problems in detail and find good agreements for all cases. The computational load on the portion of the statistical routine is comparable to or less than that for the N-body routine. The present code includes an option of hit-and-run bouncing but not fragmentation, which remains for future work. - Article
- Feb 2019
- NUMER ALGORITHMS

The integrator SSS performs accurate N-body simulations of the Solar System when there is a mix of massive bodies and test particles. The orbital motion of all bodies at all times is integrated using a 12-10 explicit Runge-Kutta Nyström (RKN) pair. The test particles are divided into sets and each set integrated on a different processor. The explicit RKN pair uses an order 12 interpolant for the position and velocity when checking for collisions. We report on two significant improvements to SSS. The first improvement reduced the local round-off error in interpolated values by approximately four orders of magnitude, permitting more accurate modelling of collisions. The technique used to reduce the round-off error can be applied to other high-order interpolants. The second improvement is hand optimization of the implementation of SSS. This optimization increased the speed of SSS by approximately 60%, permitting more accurate modelling through the use of more test particles. We also present a summary of the numerical performance of SSS on a simulation of the Sun, the planets Earth to Neptune, and 500,000 test particles over 100 million years. - Article
- Jul 2018
- ASTRON ASTROPHYS

We report the discovery of an additional substellar companion in the CoRoT-20 system based on six years of HARPS and SOPHIE radial velocity follow-up. CoRoT-20 c has a minimum mass of 17 ± 1 MJup and orbits the host star in 4.59 ± 0.05 yr, with an orbital eccentricity of 0.60 ± 0.03. This is the first identified system with an eccentric hot Jupiter and an eccentric massive companion. The discovery of the latter might be an indication of the migration mechanism of the hot Jupiter, via the Lidov-Kozai effect. We explore the parameter space to determine which configurations would trigger this type of interactions. - Article
- Feb 2018
- ASTRON ASTROPHYS

We report the discovery of four super-Earth planets around HD 215152, with orbital periods of 5.76, 7.28, 10.86, and 25.2 d, and minimum masses of 1.8, 1.7, 2.8, and 2.9 M_Earth respectively. This discovery is based on 373 high-quality radial velocity measurements taken by HARPS over 13 years. Given the low masses of the planets, the signal-to-noise ratio is not sufficient to constrain the planet eccentricities. However, a preliminary dynamical analysis suggests that eccentricities should be typically lower than about 0.03 for the system to remain stable. With two pairs of planets with a period ratio lower than 1.5, with short orbital periods, low masses, and low eccentricities, HD 215152 is similar to the very compact multi-planet systems found by Kepler, which is very rare in radial-velocity surveys. This discovery proves that these systems can be reached with the radial-velocity technique, but characterizing them requires a huge amount of observations. - Article
- Dec 2017
- ASTRON ASTROPHYS

Using high-contrast imaging with the SPHERE instrument at the VLT, we report the first images of a cold brown dwarf companion to the exoplanet host star HD4113A. The brown dwarf HD4113C is part of a complex dynamical system consisting of a giant planet, stellar host and a known wide M-dwarf companion. Its separation of $535\pm3$mas and H-band contrast of $13.35\pm0.10$mag correspond to a projected separation of 22AU and an isochronal mass estimate of $36\pm5$M$_J$ based on COND models. The companion shows strong methane absorption, and through atmospheric model fitting we estimate a surface gravity of $\log g$=5 and an effective temperature of ~500-600K. A comparison of its spectrum with observed T dwarfs indicates a late-T spectral type, with a T9 object providing the best match. By combining the observed astrometry from the imaging data with 27 years of radial velocities, we use orbital fitting to constrain its orbital and physical parameters, as well as update those of the planet HD4113Ab, discovered by previous radial velocity measurements. The data suggest a dynamical mass of $66\pm5$M$_J$ and moderate eccentricity of $0.44\pm0.08$ for the brown dwarf. This mass estimate appears to conflict with the isochronal estimate and that of similar objects, which may be caused by the newly detected object being an unresolved binary brown dwarf system or the presence of an additional object in the system. Through dynamical simulations we show that the planet may undergo strong Lidov-Kozai cycles, raising the possibility that it formed on a quasi-circular orbit and gained its currently observed high eccentricity through interactions with the brown dwarf. Follow-up observations combining radial velocities, direct imaging and Gaia astrometry will be crucial to precisely constrain the dynamical mass of the brown dwarf and allow for in-depth comparison with evolutionary and atmospheric models. - Article
- Oct 2017
- PUBL ASTRON SOC JPN

We have newly developed a parallelized particle-particle particle-tree code for planet formation, PENTACLE, which is a parallelized hybrid N-body integrator executed on a CPU-based (super)computer. PENTACLE uses a fourth-order Hermite algorithm to calculate gravitational interactions between particles within a cut-off radius and a Barnes-Hut tree method for gravity from particles beyond. It also implements an open-source library designed for full automatic parallelization of particle simulations, FDPS (Framework for Developing Particle Simulator), to parallelize a Barnes-Hut tree algorithm for a memory-distributed supercomputer. These allow us to handle 1-10 million particles in a high-resolution N-body simulation on CPU clusters for collisional dynamics, including physical collisions in a planetesimal disc. In this paper, we show the performance and the accuracy of PENTACLE in terms of Rcut and a time-step Δt. It turns out that the accuracy of a hybrid N-body simulation is controlled through Δt/Rcut and Δt/Rcut ∼ 0.1 is necessary to simulate accurately the accretion process of a planet for ≥10 6 yr. For all those interested in large-scale particle simulations, PENTACLE, customized for planet formation, will be freely available from https://github.com/PENTACLE-Team/PENTACLE under the MIT licence. © The Author 2017. Published by Oxford University Press on behalf of the Astronomical Society of Japan. All rights reserved. - Article
- Apr 2017
- ASTRON ASTROPHYS

Context. Low-mass stars are currently the best targets when searching for rocky planets in the habitable zone of their host star. Over the last 13 yr, precise radial velocities measured with the HARPS spectrograph have identified over a dozen super-Earths and Earth-mass planets (m sin i ≤ 10 M⊕) around M dwarfs, with a well-understood selection function. This well-defined sample provides information on their frequency of occurrence and on the distribution of their orbital parameters, and therefore already constrains our understanding of planetary formation. The subset of these low-mass planets that were found within the habitable zone of their host star also provide prized targets for future searches of atmospheric biomarkers. Aims. We are working to extend this planetary sample to lower masses and longer periods through dense and long-term monitoring of the radial velocity of a small M dwarf sample. Methods. We obtained large numbers of HARPS spectra for the M dwarfs GJ 3138, GJ 3323, GJ 273, GJ 628, and GJ 3293, from which we derived radial velocities (RVs) and spectroscopic activity indicators. We searched for variabilities, periodicities, Keplerian modulations, and correlations, and attribute the radial-velocity variations to combinations of planetary companions and stellar activity. Results. We detect 12 planets, 9 of which are new with masses ranging from 1.17 to 10.5 M⊕. These planets have relatively short orbital periods (P < 40 d), except for two that have periods of 217.6 and 257.8 days. Among these systems, GJ 273 harbor two planets with masses close to the Earth's. With a distance of only 3.8 parsec, GJ 273 is the second nearest known planetary system - after Proxima Centauri - with a planet orbiting the circumstellar habitable zone. - Article
- Apr 2016

Collisionless N-body simulations over tens of millions of years are an important tool in understanding the early evolution of planetary systems. We first present a CUDA kernel for evaluating the gravitational acceleration of N bodies that is intended primarily for when N is less than several thousand. We then use the kernel with a variable-order, variable-stepsize Adams method to perform long, collisionless simulations of the Solar System near limiting precision. The varying stepsize means no special scheme is required to integrate close encounters, and the motion of bodies on eccentric orbits or close to the Sun is calculated accurately. Our method is significantly more accurate than symplectic methods and sufficiently fast. - ArticleFull-text available
- Jan 2016

The fates of planetary systems provide unassailable insights into their formation and represent rich cross-disciplinary dynamical laboratories. Mounting observations of post-main-sequence planetary systems necessitate a complementary level of theoretical scrutiny. Here, I review the diverse dynamical processes which affect planets, asteroids, comets and pebbles as their parent stars evolve into giant branch, white dwarf and neutron stars. This reference provides a foundation for the interpretation and modelling of currently known systems and upcoming discoveries. - Article
- Oct 2015
- ASTRON ASTROPHYS

We know now from radial velocity surveys and transit space missions that planets only a few times more massive than our Earth are frequent around solar-type stars. Fundamental questions about their formation history, physical properties, internal structure, and atmosphere composition are, however, still to be solved. We present here the detection of a system of four low-mass planets around the bright (V = 5.5) and close-by (6.5 pc) star HD 219134. This is the first result of the Rocky Planet Search programme with HARPS-N on the Telescopio Nazionale Galileo in La Palma. The inner planet orbits the star in 3.0935 ± 0.0003 days, on a quasi-circular orbit with a semi-major axis of 0.0382 ± 0.0003 AU. Spitzer observations allowed us to detect the transit of the planet in front of the star making HD 219134 b the nearest known transiting planet to date. From the amplitude of the radial velocity variation (2.25 ± 0.22 ms-1) and observed depth of the transit (359 ± 38 ppm), the planet mass and radius are estimated to be 4.36 ± 0.44 M⊕ and 1.606 ± 0.086 R⊕, leading to a mean density of 5.76 ± 1.09 g cm-3, suggesting a rocky composition. One additional planet with minimum-mass of 2.78 ± 0.65 M⊕ moves on a close-in, quasi-circular orbit with a period of 6.767 ± 0.004 days. The third planet in the system has a period of 46.66 ± 0.08 days and a minimum-mass of 8.94 ± 1.13 M⊕, at 0.233 ± 0.002 AU from the star. Its eccentricity is 0.46 ± 0.11. The period of this planet is close to the rotational period of the star estimated from variations of activity indicators (42.3 ± 0.1 days). The planetary origin of the signal is, however, thepreferred solution as no indication of variation at the corresponding frequency is observed for activity-sensitive parameters. Finally, a fourth additional longer-period planet of mass of 71 M⊕ orbits the star in 1842 days, on an eccentric orbit (e = 0.34 ± 0.17) at a distance of 2.56 AU. - Article
- Nov 2014
- MON NOT R ASTRON SOC

Although 25–50 per cent of white dwarfs (WDs) display evidence for remnant planetary systems, their orbital architectures and overall sizes remain unknown. Vibrant close-in (≃1 R⊙) circumstellar activity is detected at WDs spanning many Gyr in age, suggestive of planets further away. Here we demonstrate how systems with 4 and 10 closely packed planets that remain stable and ordered on the main sequence can become unpacked when the star evolves into a WD and experience pervasive inward planetary incursions throughout WD cooling. Our full-lifetime simulations run for the age of the Universe and adopt main-sequence stellar masses of 1.5, 2.0 and 2.5 M⊙, which correspond to the mass range occupied by the progenitors of typical present-day WDs. These results provide (i) a natural way to generate an ever-changing dynamical architecture in post-main-sequence planetary systems, (ii) an avenue for planets to achieve temporary close-in orbits that are potentially detectable by transit photometry and (iii) a dynamical explanation for how residual asteroids might pollute particularly old WDs.

- Article
- Sep 1991
- CELEST MECH DYN ASTR

We discuss the use of symplectic integration algorithms in long-term integrations in the field of celestial mechanics. The methods' advantages and disadvantages (with respect to more common integration methods) are discussed. The numerical performance of the algorithms is evaluated using the 2-body and circular restricted 3-body problems. Symplectic integration methods have the advantages of linear phase error growth in the 2-body problem (unlike most other methods), good conservation of the integrals of the motion, good performance for moderately eccentric orbits, and ease of use. Its disadvantages include a relatively large number of force evaluations and an inability to continuously vary the step size. - Article
- Jan 2011

The standard planetesimal model of terrestrial-planet formation is based on astronomical and cosmochemical observations, and the results of laboratory experiments and numerical simulations. In this model, planets grow in a series of stages beginning with the micrometer-sized dust grains observed in protoplanetary disks. Dust grains readily stick together to form millimeter- to centimeter-sized aggregates, some of which are heated to form chondrules. Growth beyond meter size via pairwise sticking is problematic, especially in a turbulent disk. Turbulence also prevents the direct formation of planetesimals in a gravitationally unstable dust layer. Turbulent concentration can lead to the formation of gravitationally bound clumps that become 10-1000-km planetesimals. Dynamical interactions between planetsimals give the largest objects the most favorable orbits for further growth, leading to runaway and oligarchic growth and the formation of Moon- to Mars-sized planetary embryos. Large embryos acquire substantial atmospheres, speeding up planetesimal capture. Embryos also interact tidally with the gas disk, leading to orbit modification and migration. Oligarchic growth ceases when planetesimals become depleted. Embryos develop crossing orbits, and occasionally collide, leaving a handful of terrestrial planets on widely spaced orbits. The Moon probably formed via one such collision. Most stages of planet formation probably took place in the asteroid belt, but dynamical perturbations from the giant planets removed the great majority of embryos and planetesimals from this region. - Article
- May 2013
- MON NOT R ASTRON SOC

We report on the stability of hypothetical super-Earths in the habitable zone of known multiplanetary systems. Most of them have not yet been studied in detail concerning the existence of additional low-mass planets. The new N-body code genga developed at the University of Zürich allows us to perform numerous N-body simulations in parallel on graphics processing units. With this numerical tool, we can study the stability of orbits of hypothetical planets in the semimajor axis and eccentricity parameter space in high resolution. Massless test particle simulations give good predictions on the extension of the stable region and show that HIP 14810 and HD 37124 do not provide stable orbits in the habitable zone. Based on these simulations, we carry out simulations of 10 M⊕ planets in several systems (HD 11964, HD 47186, HD 147018, HD 163607, HD 168443, HD 187123, HD 190360, HD 217107 and HIP 57274). They provide more exact information about orbits at the location of mean motion resonances and at the edges of the stability zones. Beside the stability of orbits, we study the secular evolution of the planets to constrain probable locations of hypothetical planets. Assuming that planetary systems are in general closely packed, we find that apart from HD 168443, all of the systems can harbour 10 M⊕ planets in the habitable zone. - Article
- Mar 1996
- Science

Orbital histories of ejecta from the terrestrial planets were numerically integrated to study their transfer to Earth. The properties of the lunar and martian meteorites are consistent with a recurrent ejection of small meteoroids as a result of impacts on their parent bodies. Long-range gravitational effects, especially secular resonances, strongly influence the orbits of many meteoroids, increasing their collision rates with other planets and the sun. These effects and collisional destruction in the asteroid belt result in shortened time scales and higher fluxes than previously believed, especially for martian meteorites. A small flux of mercurian ejecta appears possible; recovery of meteorites from the Earth and Venus is less likely. - Article
- Dec 1992
- ASTROPHYS J

We continue our investigation of the dynamical evolution and coagulation process of planetesimals With a numerical N-body scheme, we simulate gravitational scattering and physical collisions among a system of planetesimals. The results of these simulations confirm our earlier analytical results that dynamical equilibrium is attained with a velocity dispersion comparable to the surface escape velocity of those planetesimals which contribute most of the system mass. In such an equilibrium, the rate of energy transfer from the systematic shear to dispersive motion, induced by gravitational scattering, is balanced by the rate of energy dissipation resulting from physical collisions. We also confirm that dynamical friction can lead to energy equipartition between an abundant population of low-mass field planetesimals and a few collisionally induced mergers with larger masses. These effects produce mass segregation in phase space and runaway coagulation. Collisions also lead to coagulation and evolution of the mass spectrum. The mergers of two field planetesimals can provide sufficient mass differential with other planetesimals for dynamical friction to induce energy equipartition and mass segregation. For small velocity dispersions, the more massive planetesimals produce relatively large gravitational focusing factors. Consequently, the growth time scale decreases with mass and runaway coagulation is initiated. Our numerical simulations show that, provided there is sufficient supply of low-mass planetesimals, runaway coagulation can lead to the formation of protoplanetary cores with masses comparable to a significant fraction of an Earth mass. We estimate that, at 1 AU, the characteristic time scale for the initial stages of planetesimal growth is ˜104 yr and ˜105 yr for the growth to protoplanetary cores. At Jupiter's present distance, these time scales are an order of magnitude longer. - Article
- Sep 1992
- ASTRON J

Many problems in solar system dynamics are described by Hamiltomans of the form H =HKep + εHpert, ε≪ where HKep is the usual Hamiltonian for the Kepler two-body problem and εHpert represents (for example) much weaker perturbations from the planets. We review symplectic integrators for Hamiltonians of this kind, focusing on methods that exploit the integrabihty of HKep. We show that the long-term errors in these integrators can be reduced by a factor of order ε by suitable starting procedures, for example, by starting with a very small stepsize and gradually increasing the stepsize to its final value. The resulting integrators are easily the best available for a wide range of solar system problems. - Article
- Feb 1984
- ASTRON ASTROPHYS

The aim of this work is the construction of a fast integration method for differential equations (DE), especially the equations of the motion of celestial bodies. Although a number of integration schemes are available none of them seem to be adequate for treating n-body systems with variable masses, which arise in some cosmogonic problems of the early solar system. As a first step we are now able to present a high-speed numerical integration scheme of the classical n-body system. The basic idea of solving differential equations with Lie-series is due to Grobner (1967) but, unfortunately, he did not elaborate on this method and stopped after some numerically unsatisfactory results. We could simplify the calculation of the Lie-terms and derived finally a recurrence formula for the Lie-terms. Whereas Grobner tried to solve the two-body and three (n-body)problem by two different approaches we solved, at first, in an optimal way the 2-body- problem. Then we were able to derive in a quite similar way the solutions of the 3-body and n-body system. Our integration method for planetary motions has two major advantages: First, it is a relatively fast method (about the factor 3-10 faster in comparison with the n-body program by Schubart-Stumpff, which is commonly used by Astronomers). Second, because larger step lengths can be used, roundoff errors are smaller (e.g. step length 135 d for Jupiter). - Article
- May 1991
- ASTRON J

The equations of motion of the nine planets and the spin axis of the earth are integrated for 3.05 million years into the past. The initial conditions are taken from the JPL DE 102 ephemeris; dominant relativistic corrections and corrections for the earth-moon system's quadrupole moment are included; but solar mass loss is not considered. A low-pass filter is employed to remove short-term variations and isolate variations of more than 2000 years. At the end of integration, the fractional error in the earth's position is shown to be less than 0.03 radian, and larger errors of up to several radians are shown for other planets. Because the adjacent trajectories do not diverge exponentially, the planetary orbits are apparently regular. Between 1 and 2 Myr in the past, a burst of oscillations are found in the semimajor axes of Venus and the earth. The results are compared to other models, and disagreement is shown between the results of secular perturbation theory and this model's predictions beyond 1-1.5 Myr into the past. - Article
- Sep 2012

Since the formulation of the problem by Newton, and during three centuries, astronomers and mathematicians have sought to demonstrate the stability of the Solar System. Thanks to the numerical experiments of the last two decades, we know now that the motion of the planets in the Solar System is chaotic, which prohibits any accurate prediction of their trajectories beyond a few tens of millions of years. The recent simulations even show that planetary collisions or ejections are possible on a period of less than 5 billion years, before the end of the life of the Sun. - ArticleFull-text available
- Dec 2007

Symplectic integration algorithms are well suited for long-term integrations of Hamiltonian systems, because they preserve the geometric structure of the Hamiltonian flow. However, this desirable property is generally lost when adaptive time step control is added to a symplectic integrator. We describe an adaptive time step, symplectic integrator that can be used if the Hamiltonian is the sum of kinetic and potential energy components and the required time step depends only on the potential energy (e.g., test-particle integrations in fixed potentials). In particular, we describe an explicit, reversible, symplectic, leapfrog integrator for a test particle in a near-Keplerian potential; this integrator has a time step proportional to distance from the attracting mass and has the remarkable property of integrating orbits in an inverse-square force field with only "along-track" errors; i.e., the phase-space shape of a Keplerian orbit is reproduced exactly, but the orbital period is in error by O(N-2), where N is the number of steps per period. - Article
- Dec 2007

We present a new symplectic algorithm that has the desirable properties of the sophisticated but highly efficient numerical algorithms known as mixed variable symplectic (MVS) methods and that, in addition, can handle close encounters between objects. This technique is based on a variant of the standard MVS methods, but it handles close encounters by employing a multiple time step technique. When the bodies are well separated, the algorithm has the speed of MVS methods, and whenever two bodies suffer a mutual encounter, the time step for the relevant bodies is recursively subdivided to whatever level is required. We demonstrate the power of this method using several tests of the technique. We believe that this algorithm will be a valuable tool for the study of planetesimal dynamics and planet formation. - Article
- Dec 2008
- ASTROPHYS J

We investigate the formation of protoplanet systems from planetesimal disks by global (N = 5000 and 10,000 and 0.5 AU < a < 1.5 AU, where N is the number of bodies and a is the distance from a central star) N-body simulations of planetary accretion. For application to extrasolar planetary systems, we study the wide variety of planetesimal disks of the surface mass density Σsolid = Σ1(a/1 AU)-α g cm-2 with Σ1 = 1, 10, 100 and α = 1/2, 3/2, 5/2. The results are all consistent with the prediction from the "oligarchic growth" model. We derive how the growth timescale, the isolation (final) mass, and the orbital separation of protoplanets depend on the initial disk mass (Σ1) and the initial disk profile (α). The isolation mass increases in proportion to Σ, while the number of protoplanets decreases in proportion to Σ. The isolation mass depends on a as a(3/2)(2-α), which means it increases with a for α < 2 while it decreases with a for α > 2. The growth timescale increases with a but decreases with Σ1. Based on the oligarchic growth model and the conventional Jovian planet formation scenario, we discuss the diversity of planetary systems. Jovian planets can form in the disk range where the contraction timescale of planetary atmosphere and the growth timescale of protoplanets (cores) are shorter than the lifetime of the gas disk. We find that for the disk lifetime ~108 yr, several Jovian planets would form from massive disks with Σ1 30 with Uranian planets outside the Jovian planets. Only terrestrial and Uranian planets would form from light disks with Σ1 3. Solar system-like planetary systems would form from medium disks with Σ1 10. - Recent results have shown that many of the known extrasolar planetary systems contain regions that are stable for massless test particles. We examine the possibility that Saturn mass planets exist in these systems, just below the detection threshold, and predict likely orbital parameters for such unseen planets. We insert a Saturn mass planet into the regions stable for massless test particles and integrate the system for 100 million years. We conduct 200-600 of these experiments to test parameter space in 55 Cancri, HD 37124, HD 38529, and HD 74156. In 55 Cnc we find three maxima of the survival rate of Saturn mass planets, located in semimajor axis a and eccentricity e space at (a,e) = (1.0 AU, 0.02), (2.0 AU, 0.08), and (3.0 AU, 0.17). In HD 37124 the maximum lies at a = 0.90-98 AU, eccentricity e ~ 0.05-0.15. In HD 38529, only 5% of Saturn mass planets are unstable, and the region in which a Saturn mass planet could survive is very broad, centered on 0.5 < a < 0.6, e < 0.15. In HD 74156 we find a broad maximum at a = 0.9-1.4 AU, e ≤ 0.15. These orbital values are initial conditions, and not always the most likely values for detection. Several are located in the habitable zones of their parent stars and are of astrobiological interest. We suggest the possibility that companions may lie in these locations of parameter space, and encourage further observational investigation of these systems.
- Article
- Sep 2012
- ICARUS

The abundances of elements in the Earth and the terrestrial planets provide the initial conditions for life and clues as to the history and formation of the Solar System. We follow the pioneering work of Bond et al. (2010) and combine circumstellar disk models, chemical equilibrium calculations and dynamical simulations of planet formation to study the bulk composition of rocky planets. We use condensation sequence calculations to estimate the initial abundance of solids in the circumstellar disk with properties determined from time dependent theoretical models. We combine this with dynamical simulations of planetesimal growth that trace the solids during the planet formation process. We calculate the elemental abundances in the resulting planets and explore how these vary with the choice of disk model and the initial conditions within the Solar Nebula. Although certain characteristics of the terrestrial planets in the Solar System could be reproduced, none of our models could reproduce the abundance properties of all the planets. We found that the choice of the initial planetesimal disk mass and of the disk model has a significant effect on composition gradients. Disk models that give higher pressure and temperature result in larger variations in the bulk chemical compositions of the resulting planets due to inhomogeneities in the element abundance profiles and due to the different source regions of the planets in the dynamical simulations. We observed a trend that massive planets and planets with relatively small semi-major axes are more sensitive to these variations than smaller planets at larger radial distance. Only these large variations in the simulated chemical abundances can potentially explain the diverse bulk composition of the Solar System planets, whereas Mercury's bulk composition can not be reproduced in our approach. - Article
- Aug 2012
- NEW ASTRON

We present Swarm-NG, a C++ library for the efficient direct integration of many n-body systems using highly-parallel Graphics Processing Unit (GPU), such as NVIDIA's Tesla T10 and M2070 GPUs. While previous studies have demonstrated the benefit of GPUs for n-body simulations with thousands to millions of bodies, Swarm-NG focuses on many few-body systems, e.g., thousands of systems with 3...15 bodies each, as is typical for the study of planetary systems. Swarm-NG parallelizes the simulation, including both the numerical integration of the equations of motion and the evaluation of forces using NVIDIA's "Compute Unified Device Architecture" (CUDA) on the GPU. Swarm-NG includes optimized implementations of 4th order time-symmetrized Hermite integration and mixed variable symplectic integration, as well as several sample codes for other algorithms to illustrate how non-CUDA-savvy users may themselves introduce customized integrators into the Swarm-NG framework. To optimize performance, we analyze the effect of GPU-specific parameters on performance under double precision. Applications of Swarm-NG include studying the late stages of planet formation, testing the stability of planetary systems and evaluating the goodness-of-fit between many planetary system models and observations of extrasolar planet host stars (e.g., radial velocity, astrometry, transit timing). While Swarm-NG focuses on the parallel integration of many planetary systems,the underlying integrators could be applied to a wide variety of problems that require repeatedly integrating a set of ordinary differential equations many times using different initial conditions and/or parameter values. - Article
- Oct 2008
- NEW ASTRON

We present sixth- and eighth-order Hermite integrators for astrophysical N-body simulations, which use the derivatives of accelerations up to second-order (snap) and third-order (crackle). These schemes do not require previous values for the corrector, and require only one previous value to construct the predictor. Thus, they are fairly easy to implement. The additional cost of the calculation of the higher-order derivatives is not very high. Even for the eighth-order scheme, the number of floating-point operations for force calculation is only about two times larger than that for traditional fourth-order Hermite scheme. The sixth-order scheme is better than the traditional fourth-order scheme for most cases. When the required accuracy is very high, the eighth-order one is the best. These high-order schemes have several practical advantages. For example, they allow a larger number of particles to be integrated in parallel than the fourth-order scheme does, resulting in higher execution efficiency in both general-purpose parallel computers and GRAPE systems. - Article
- Jul 2012
- J COMPUT PHYS

We present a new implementation of the numerical integration of the classical, gravitational, N-body problem based on a high order Hermite's integration scheme with block time steps, with a direct evaluation of the particle-particle forces. The main innovation of this code (called HiGPUs) is its full parallelization, exploiting both OpenMP and MPI in the use of the multicore Central Processing Units as well as either Compute Unified Device Architecture (CUDA) or OpenCL for the hosted Graphic Processing Units. We tested both performance and accuracy of the code using up to 256 GPUs in the supercomputer IBM iDataPlex DX360M3 Linux Infiniband Cluster provided by the italian supercomputing consortium CINECA, for values of N up to 8 millions. We were able to follow the evolution of a system of 8 million bodies for few crossing times, task previously unreached by direct summation codes. The code is freely available to the scientific community. - Article
- Aug 2001
- Icarus

The results of 16 new 3D N-body simulations of the final stage of the formation of the terrestrial planets are presented. These N-body integrations begin with 150–160 lunar-to-Mars size planetary embryos, with semi-major axes 0.3 < a < 2.0 AU, and include per-turbations from Jupiter and Saturn. Two initial mass distributions are examined: approximately uniform masses, and a bimodal dis-tribution with a few large and many small bodies. In most of the integrations, systems of three or four terrestrial planets form within about 200 million years. These planets have orbital separations sim-ilar to the terrestrial planets, and the largest body contains 1/3–2/3 of the surviving mass. The final planets typically have larger ec-centricities, e, and inclinations, i than the time-averaged values for Earth and Venus. However, the values of e and i are lower than in earlier N-body integrations which started with fewer embryos. The spin axes of the planets have approximately random orientations, unlike the terrestrial planets, and the high degree of mass concen-tration in the region occupied by Earth and Venus is not reproduced in any of the simulations. The principal effect of using an initially bimodal mass distribution is to increase the final number of planets. Each simulation ends with an object that is an approximate ana-logue of Earth in terms of mass and heliocentric distance. These Earth analogues reach 50% (90%) of their final mass with a me-dian time of 20 (50) million years, and they typically accrete some material from all portions of the disk. - Article
- May 2006
- ASTROPHYS J

The final stage of terrestrial planet formation is known as the giant impact stage, where protoplanets collide with one another to form planets. As this process is stochastic, in order to clarify it, it is necessary to quantify it statistically. We investigate this final assemblage of terrestrial planets from protoplanets using N-body simulations. As initial conditions, we adopt the oligarchic growth model of protoplanets. We systematically change the surface density, surface density profile, and orbital separation of the initial protoplanet system, and the bulk density of protoplanets, while the initial system radial range is fixed at 0.5–1.5 AU. For each initial condition, we perform 20 runs, and from their results we derive the statistical properties of the assembled planets. For the standard disk model, typically two Earth-sized planets form in the terrestrial planet region. We show the dependences of the masses and orbital elements of planets on the initial protoplanet system parameters and give their simple empirical fits. The number of planets slowly decreases as the surface density of the initial protoplanets increases, while the masses of individual planets in-crease almost linearly. For a steeper surface density profile, large planets tend to form closer to the star. For the param-eter ranges that we test, the basic structure of planetary systems depends only slightly on the initial distribution of protoplanets and the bulk density as long as the total mass is fixed. - Article
- Feb 2009

Abstract— We have examined the fate of impact ejecta liberated from the surface of Mercury due to impacts by comets or asteroids, in order to study 1) meteorite transfer to Earth, and 2) reaccumulation of an expelled mantle in giant-impact scenarios seeking to explain Mercury's large core. In the context of meteorite transfer during the last 30 Myr, we note that Mercury's impact ejecta leave the planet's surface much faster (on average) than other planets in the solar system because it is the only planet where impact speeds routinely range from 5 to 20 times the planet's escape speed; this causes impact ejecta to leave its surface moving many times faster than needed to escape its gravitational pull. Thus, a large fraction of Mercurian ejecta may reach heliocentric orbit with speeds sufficiently high for Earth-crossing orbits to exist immediately after impact, resulting in larger fractions of the ejecta reaching Earth as meteorites. We calculate the delivery rate to Earth on a time scale of 30 Myr (typical of stony meteorites from the asteroid belt) and show that several percent of the high-speed ejecta reach Earth (a factor of 2–3 less than typical launches from Mars); this is one to two orders of magnitude more efficient than previous estimates. Similar quantities of material reach Venus.These calculations also yield measurements of the re-accretion time scale of material ejected from Mercury in a putative giant impact (assuming gravity is dominant). For Mercurian ejecta escaping the gravitational reach of the planet with excess speeds equal to Mercury's escape speed, about one third of ejecta reaccretes in as little as 2 Myr. Thus collisional stripping of a silicate proto-Mercurian mantle can only work effectively if the liberated mantle material remains in small enough particles that radiation forces can drag them into the Sun on time scale of a few million years, or Mercury would simply re-accrete the material. - Article
- Aug 1999
- CELEST MECH DYN ASTR

By Hamiltonian manipulation we demonstrate the existence of separable time‐transformed Hamiltonians in the extended phase‐space. Due to separability explicit symplectic methods are available for the solution of the equations of motion. If the simple leapfrog integrator is used, in case of two‐body motion, the method produces an exact Keplerian ellipse in which only the time‐coordinate has an error. Numerical tests show that even the rectilinear N‐body problem is feasible using only the leapfrog integrator. In practical terms the method cannot compete with regularized codes, but may provide new directions for studies of symplectic N‐body integration. - Article
- Jan 1996
- CELEST MECH DYN ASTR

Large scale chaos is present everywhere in the solar system. It plays a major role in the sculpting of the asteroid belt and in the diffusion of comets from the outer region of the solar system. All the inner planets probably experienced large scale chaotic behavior for their obliquities during their history. The Earth obliquity is presently stable only because of the presence of the Moon, and the tilt of Mars undergoes large chaotic variations from 0 to about 60. On billion years time scale, the orbits of the planets themselves present strong chaotic variations which can lead to the escape of Mercury or collision with Venus in less than 3.5 Gyr. The organization of the planets in the solar system thus seems to be strongly related to this chaotic evolution, reaching at all time a state of marginal stability, that is practical stability on a time-scale comparable to its age. - Article
- May 2012
- MON NOT R ASTRON SOC

We describe the use of Graphics Processing Units (GPUs) for speeding up the code NBODY6 which is widely used for direct $N$-body simulations. Over the years, the $N^2$ nature of the direct force calculation has proved a barrier for extending the particle number. Following an early introduction of force polynomials and individual time-steps, the calculation cost was first reduced by the introduction of a neighbour scheme. After a decade of GRAPE computers which speeded up the force calculation further, we are now in the era of GPUs where relatively small hardware systems are highly cost-effective. A significant gain in efficiency is achieved by employing the GPU to obtain the so-called regular force which typically involves some 99 percent of the particles, while the remaining local forces are evaluated on the host. However, the latter operation is performed up to 20 times more frequently and may still account for a significant cost. This effort is reduced by parallel SSE/AVX procedures where each interaction term is calculated using mainly single precision. We also discuss further strategies connected with coordinate and velocity prediction required by the integration scheme. This leaves hard binaries and multiple close encounters which are treated by several regularization methods. The present nbody6-GPU code is well balanced for simulations in the particle range $10^4-2 \times 10^5$ for a dual GPU system attached to a standard PC. - Article
- Aug 2012
- ICARUS

It has been suggested that the ejection of terrestrial crustal material to interplanetary space, accelerated in a large impact, may result in the interchange of biological material between Earth and other Solar System bodies. In this paper, we analyze the fate of debris ejected from Earth by means of direct numerical simulations of the dynamics of a large collection of test particles. This allows us to determine the probability and conditions for the collision of Earth ejecta with other planets of the Solar System. We also estimate the amount of particles falling back to Earth and colliding with the Moon as a function of time after being ejected. The Mercury-6 code is used to compute the dynamics of test particles under the gravitational effect of the planets in the Solar System and the Sun. A series of simulations are conducted with different ejection speeds, considering more than 10(5) particles in each case. We find that in general, the collision rates of Earth ejecta with Venus and the Moon, as well as the fall-back rates, are within an order of magnitude of results reported in the literature. By considering a larger number of particles than in all previous calculations we have also determined, on the basis of direct numerical simulations, the collision probability with Mars and, for the first time, computed collision probabilities with Jupiter and Saturn. We find that the collision probability with Mars is greater than values determined from collision cross section estimations previously reported. - Article
- Apr 2012
- EUR PHYS J-SPEC TOP

In this short review we present the developments over the last 5 decades that have led to the use of Graphics Processing Units (GPUs) for astrophysical simulations. Since the introduction of NVIDIA's Compute Unified Device Architecture (CUDA) in 2007 the GPU has become a valuable tool for N-body simulations and is so popular these days that almost all papers about high precision N-body simulations use methods that are accelerated by GPUs. With the GPU hardware becoming more advanced and being used for more advanced algorithms like gravitational tree-codes we see a bright future for GPU like hardware in computational astrophysics. - Article
- Feb 2008
- NEW ASTRON

We present the results of gravitational direct N-body simulations using the graphics processing unit (GPU) on a commercial NVIDIA GeForce 8800GTX designed for gaming computers. The force evaluation of the N-body problem is implemented in “Compute Unified Device Architecture” (CUDA) using the GPU to speedup the calculations. We tested the implementation on three different N-body codes: two direct N-body integration codes, using the 4th order predictor–corrector Hermite integrator with block time-steps, and one Barnes-Hut treecode, which uses a 2nd order leapfrog integration scheme. The integration of the equations of motions for all codes is performed on the host CPU.We find that for N > 512 particles the GPU outperforms the GRAPE-6Af, if some softening in the force calculation is accepted. Without softening and for very small integration time-steps the GRAPE still outperforms the GPU. We conclude that modern GPUs offer an attractive alternative to GRAPE-6Af special purpose hardware. Using the same time-step criterion, the total energy of the N-body system was conserved better than to one in 106 on the GPU, only about an order of magnitude worse than obtained with GRAPE-6Af. For N ≳ 105 the 8800GTX outperforms the host CPU by a factor of about 100 and runs at about the same speed as the GRAPE-6Af. - Article
- Aug 2006
- ICARUS

The final stage in the formation of terrestrial planets consists of the accumulation of ∼1000-km “planetary embryos” and a swarm of billions of 1–10 km “planetesimals.” During this process, water-rich material is accreted by the terrestrial planets via impacts of water-rich bodies from beyond roughly 2.5 AU. We present results from five high-resolution dynamical simulations. These start from 1000–2000 embryos and planetesimals, roughly 5–10 times more particles than in previous simulations. Each simulation formed 2–4 terrestrial planets with masses between 0.4 and 2.6 Earth masses. The eccentricities of most planets were ∼0.05, lower than in previous simulations, but still higher than for Venus, Earth and Mars. Each planet accreted at least the Earth's current water budget. We demonstrate several new aspects of the accretion process: (1) The feeding zones of terrestrial planets change in time, widening and moving outward. Even in the presence of Jupiter, water-rich material from beyond 2.5 AU is not accreted for several millions of years. (2) Even in the absence of secular resonances, the asteroid belt is cleared of >99% of its original mass by self-scattering of bodies into resonances with Jupiter. (3) If planetary embryos form relatively slowly, then the formation of embryos in the asteroid belt may have been stunted by the presence of Jupiter. (4) Self-interacting planetesimals feel dynamical friction from other small bodies, which has important effects on the eccentricity evolution and outcome of a simulation. - Article
- Oct 2009
- NEW ASTRON

We present Sapporo, a library for performing high precision gravitational N-body simulations on NVIDIA graphical processing units (GPUs). Our library mimics the GRAPE-6 library, and N-body codes currently running on GRAPE-6 can switch to Sapporo by a simple relinking of the library. The precision of our library is comparable to that of GRAPE-6, even though internally the GPU hardware is limited to single precision arithmetics. This limitation is effectively overcome by emulating double precision for calculating the distance between particles. The performance loss of this operation is small (≲20%) compared to the advantage of being able to run at high precision. We tested the library using several GRAPE-6-enabled N-body codes, in particular with Starlab and phiGRAPE. We measured peak performance of 800 Gflop/s for running with 106 particles on a PC with four commercial G92 architecture GPUs (two GeForce 9800GX2). As a production test, we simulated a 32 k Plummer model with equal-mass stars well beyond core collapse. The simulation took 41 days, during which the mean performance was 113 Gflop/s. The GPU did not show any problems from running in a production environment for such an extended period of time. - Article
- Jan 2000
- ICARUS

We describe a new direct numerical method for simulating planetesimal dynamics in which N∼106 or more bodies can be evolved simultaneously in three spatial dimensions over hundreds of dynamical times. This represents several orders of magnitude improvement in resolution over previous studies. The advance is made possible through modification of a stable and tested cosmological code optimized for massively parallel computers. However, owing to the excellent scalability and portability of the code, modest clusters of workstations can treat problems with N∼105 particles in a practical fashion.The code features algorithms for detection and resolution of collisions and takes into account the strong central force field and flattened Keplerian disk geometry of planetesimal systems. We demonstrate the range of problems that can be addressed by presenting simulations that illustrate oligarchic growth of protoplanets, planet formation in the presence of giant planet perturbations, the formation of the jovian moons, and orbital migration via planetesimal scattering. We also describe methods under development for increasing the timescale of the simulations by several orders of magnitude. - Article
- Nov 2007
- NEW ASTRON

We present the results of gravitational direct N-body simulations using the commercial graphics processing units (GPU) NVIDIA Quadro FX1400 and GeForce 8800GTX, and compare the results with GRAPE-6Af special purpose hardware. The force evaluation of the N-body problem was implemented in Cg using the GPU directly to speed-up the calculations. The integration of the equations of motions were, running on the host computer, implemented in C using the 4th order predictor–corrector Hermite integrator with block time steps. We find that for a large number of particles (N ≳ 104) modern graphics processing units offer an attractive low cost alternative to GRAPE special purpose hardware. A modern GPU continues to give a relatively flat scaling with the number of particles, comparable to that of the GRAPE. The GRAPE is designed to reach double precision, whereas the GPU is intrinsically single-precision. For relatively large time steps, the total energy of the N-body system was conserved better than to one in 106 on the GPU, which is impressive given the single-precision nature of the GPU. For the same time steps, the GRAPE gave somewhat more accurate results, by about an order of magnitude. However, smaller time steps allowed more energy accuracy on the grape, around 10−11, whereas for the GPU machine precision saturates around 10−6 For N ≳ 106 the GeForce 8800GTX was about 20 times faster than the host computer. Though still about a factor of a few slower than GRAPE, modern GPUs outperform GRAPE in their low cost, long mean time between failure and the much larger onboard memory; the GRAPE-6Af holds at most 256k particles whereas the GeForce 8800GTX can hold 9 million particles in memory. - Article
- Nov 1999
- ICARUS

We perform three-dimensional N-body integrations of the final stages of terrestrial planet formation. We report the results of 10 simulations beginning with 22–50 initial planetary embryos spanning the range 0.5–1.5 AU, each with an initial mass of 0.04–0.13M⊕. Collisions are treated as inelastic mergers. We follow the evolution of each system for 2×108 years at which time a few terrestrial type planets remain. On average, our simulations produced two planets larger than 0.5M⊕ in the terrestrial region (1 simulation with one m≥0.5M⊕ planet, 8 simulations with two m≥0.5M⊕ planets, and 1 simulation with three m≥0.5M⊕ planets). These Earth-like planets have eccentricities and orbital spacing considerably larger than the terrestrial planets of comparable mass (e.g., Earth and Venus). We also examine the angular momentum contributions of each collision to the final spin angular momentum of a planet, with an emphasis on the type of impact which is believed to have triggered the formation of the Earth's Moon. There was an average of two impacts per simulation that contributed more angular momentum to a planet than is currently present in the Earth/Moon system. We determine the spin angular momentum states of the growing planets by summing the contributions from each collisional encounter. Our results show that the spin angular momentum states of the final planets are generally the result of contributions made by the last few large impacts. Our results suggest that the current angular momentum of the Earth/Moon system may be the result of more than one large impact rather than a single impact. Further, upon suffering their first collision, the planetary embryos in our simulations are spinning rapidly throughout the final accretion of the planets, suggesting the proto-Earth may have been rotating rapidly prior to the Moon-forming impact event. - Article
- Nov 1990
- PHYS LETT A

For Hamiltonian systems of the form H = T(p)+V(q) a method is shown to construct explicit and time reversible symplectic integrators of higher order. For any even order there exists at least one symplectic integrator with exact coefficients. The simplest one is the 4th order integrator which agrees with one found by Forest and by Neri. For 6th and 8th orders, symplectic integrators with fewer steps are obtained, for which the coefficients are given by solving a set of simultaneous algebraic equations numerically. - Article
- Dec 1998
- ICARUS

We simulate the late stages of terrestrial-planet formation using N-body integrations, in three dimensions, of disks of up to 56 initially isolated, nearly coplanar planetary embryos, plus Jupiter and Saturn. Gravitational perturbations between embryos increase their eccentricities,e, until their orbits become crossing, allowing collisions to occur. Further interactions produce large-amplitude oscillations ineand the inclination,i, with periods of ∼105years. These oscillations are caused by secular resonances between embryos and prevent objects from becoming re-isolated during the simulations. The largest objects tend to maintain smallereandithan low-mass bodies, suggesting some equipartition of random orbital energy, but accretion proceeds by orderly growth. The simulations typically produce two large planets interior to 2 AU, whose time-averagedeandiare significantly larger than Earth and Venus. The accretion rate falls off rapidly with heliocentric distance, and embryos in the “Mars zone” (1.2 <a< 2 AU) are usually scattered inward and accreted by “Earth” or “Venus,” or scattered outward and removed by resonances, before they can accrete one another. The asteroid belt (a> 2 AU) is efficiently cleared as objects scatter one another into resonances, where they are lost via encounters with Jupiter or collisions with the Sun, leaving, at most, one surviving object. Accretional evolution is complete after 3 × 108years in all simulations that include Jupiter and Saturn. The number and spacing of the final planets, in our simulations, is determined by the embryos' eccentricities, and the amplitude of secular oscillations ine, prior to the last few collision events. - Article
- Sep 2006
- ICARUS

We have performed 8 numerical simulations of the final stages of accretion of the terrestrial planets, each starting with over 5× more gravitationally interacting bodies than in any previous simulations. We use a bimodal initial population spanning the region from 0.3 to 4 AU with 25 roughly Mars-mass embryos and an equal mass of material in a population of ∼1000 smaller planetesimals, consistent with models of the oligarchic growth of protoplanetary embryos. Given the large number of small planetesimals in our simulations, we are able to more accurately treat the effects of dynamical friction during the accretion process. We find that dynamical friction can significantly lower the timescales for accretion of the terrestrial planets and leads to systems of terrestrial planets that are much less dynamically excited than in previous simulations with fewer initial bodies. In addition, we study the effects of the orbits of Jupiter and Saturn on the final planetary systems by running 4 of our simulations with the present, eccentric orbits of Jupiter and Saturn (the EJS simulations) and the other 4 using a nearly circular and co-planar Jupiter and Saturn as predicted in the Nice Model of the evolution of the outer Solar System [Gomes, R., Levison, H.F., Tsiganis, K., Morbidelli, A., 2005. Nature 435, 466–469; Tsiganis, K., Gomes, R., Morbidelli, A., Levison, H.F., 2005. Nature 435, 459–461; Morbidelli, A., Levison, H.F., Tsiganis, K., Gomes, R., 2005. Nature 435, 462–465] (the CJS simulations). Our EJS simulations provide a better match to our Solar System in terms of the number and average mass of the final planets and the mass-weighted mean semi-major axis of the final planetary systems, although increased dynamical friction can potentially improve the fit of the CJS simulations as well. However, we find that in our EJS simulations, essentially no water-bearing material from the outer asteroid belt ends up in the final terrestrial planets, while a large amount is delivered in the CJS simulations. In addition, the terrestrial planets in the EJS simulations receive a late veneer of material after the last giant impact event that is likely too massive to reconcile with the siderophile abundances in the Earth's mantle, while the late veneer in the CJS simulations is much more consistent with geochemical evidence. - Article
- Jun 2010
- ICARUS

We present results from a suite of N-body simulations that follow the accretion history of the terrestrial planets using a new parallel treecode that we have developed. We initially place 2000 equal size planetesimals between 0.5--4.0 AU and the collisional growth is followed until the completion of planetary accretion (> 100 Myr). All the important effect of gas in laminar disks are taken into account: the aerodynamic gas drag, the disk-planet interaction including Type I migration, and the global disk potential which causes inward migration of secular resonances as the gas dissipates. We vary the initial total mass and spatial distribution of the planetesimals, the time scale of dissipation of nebular gas, and orbits of Jupiter and Saturn. We end up with one to five planets in the terrestrial region. In order to maintain sufficient mass in this region in the presence of Type I migration, the time scale of gas dissipation needs to be 1-2 Myr. The final configurations and collisional histories strongly depend on the orbital eccentricity of Jupiter. If today's eccentricity of Jupiter is used, then most of bodies in the asteroidal region are swept up within the terrestrial region owing to the inward migration of the secular resonance, and giant impacts between protoplanets occur most commonly around 10 Myr. If the orbital eccentricity of Jupiter is close to zero, as suggested in the Nice model, the effect of the secular resonance is negligible and a large amount of mass stays for a long period of time in the asteroidal region. With a circular orbit for Jupiter, giant impacts usually occur around 100 Myr, consistent with the accretion time scale indicated from isotope records. However, we inevitably have an Earth size planet at around 2 AU in this case. It is very difficult to obtain spatially concentrated terrestrial planets together with very late giant impacts. - Article
- Mar 2003
- ICARUS

Mounting attention has focused on interplanetary transfer of microorganisms (panspermia), particularly in reference to exchange between Mars and Earth. In most cases, however, such exchange requires millions of years, over which time the transported microorganisms must remain viable. During a large impact on Earth, however, previous work (J.C. Armstrong et al., 2002, Icarus 160, 183–196) has shown that substantial amounts of material return to the planet of origin over a much shorter period of time (< 5000 years), considerably mitigating the challenges to the survival of a living organism. Conservatively evaluating experiments performed [by others] on Bacillus subtilis and Deinococcus radiodurans to constrain biological survival under impact conditions, we estimate that if the Earth were hit by a sterilizing impactor ∼ 300 km in diameter, with a relative velocity of 30 km s−1 (such as may have occurred during the Late Heavy Bombardment), an initial cell population in the ejecta of order 103–105 cells kg−1 would in most cases be sufficient for a single modern organism to survive and return to an again-clement planet 3000–5000 years later. Although little can be said about the characteristics or distribution of ancient life, our calculations suggest that impact reseeding is a possible means by which life, if present, could have survived the Late Heavy Bombardment. - We have undertaken a thorough dynamical investigation of five extrasolar planetary systems using extensive numerical experiments on the supercomputer of the Max Planck Institute for Gravitational Physics (Albert Einstein Institute). This work was performed as part of the Helmholtz Institute for Supercomputational Physics Summer School on “Chaos and Stability in Planetary Systems” (2003). The systems Gl 777 A, HD 72659, Gl 614, 47 Uma and HD 4208 were examined concerning the question of whether they could host terrestrial like planets in their habitable zones (=HZ). First we investigated the mean motion resonances between fictitious terrestrial planets and the existing gas giants in these five extrasolar systems. Then a fine grid of initial conditions for a potential terrestrial planet within the HZ was chosen for each system, from which the stability of orbits was then assessed by direct integrations over a time interval of 1 million years. For each of the five systems the 2-dimensional grid of initial conditions contained 80 eccentricity points for the Jovian planet and up to 160 semimajor axis points for the fictitious planet. The computations were carried out using a Lie-series integration method with an adaptive step size control. This integration method achieves machine precision accuracy in a highly efficient and robust way, requiring no special adjustments when the orbits have large eccentricities. The stability of orbits was examined with a determination of the R ́enyi entropy, estimated from recurrence plots, and with a more straight forward method based on the maximum eccentricity achieved by the planet over the 1 million year integration. Additionally, the eccentricity is an indication of the habitability of a terrestrial planet in the HZ; any value of e > 0.2 produces a significant temperature difference on a planet’s surface between apoapse and periapse. The results for possible stable orbits for terrestrial planets in habitable zones for the five systems are summarized as follows: for Gl 777 A nearly the entire HZ is stable, for 47 Uma, HD 72659 and HD 4208 terrestrial planets can survive for a sufficiently long time, while for Gl 614 our results exclude terrestrial planets moving in stable orbits within the HZ. Studies such as this one are of primary interest to future space missions dedicated to finding habitable terrestrial planets in other stellar systems. Assessing the likelihood of other habitable planets, and more generally the possibility of other life, is the central question of astrobiology today. Our investigation indicates that, from the dynamical point of view, habitable terrestrial planets seem to be quite compatible with many of the currently discovered extrasolar systems. Cited By (since 1996): 36, Export Date: 15 September 2011, Source: Scopus
- BookFull-text available
- Jan 1993

- Article
- Mar 1989
- NATURE

The results are described of a numerical integration, extending backwards over 200 million years, of an extensive analytic system of averaged differential equations containing the secular evolution of the orbits of the eight main planets. The solution is chaotic, with a maximum Lyapunov exponent that reaches the surprisingly large value of about 1.5/Myr. The motion of the solar system is thus shown to be chaotic, not quasi-periodic. In particular, predictability of the orbits of the inner planets, including the earth, is lost within a few tens of millions of years. - Article
- Nov 1991
- ASTRON J

The present study generalizes the mapping method of Wisdom (1982) to encompass all gravitational n-body problems with a dominant central mass. The rationale for the generalized mapping method is discussed as well as details for the mapping for the n-body problem. Some refinements of the method are considered, and the relationship of the mapping method to other symplectic integration methods is shown. The method is used to compute the evolution of the outer planets for a billion years. The resulting evolution is compared to the 845 million year evolution of the outer planets performed on the Digital Orerry using standard numerical integration techniques. This calculation provides independent numerical confirmation of the result of Sussman and Wisdom (1988) that the motion of the planet Pluto is chaotic. - Article
- Feb 1999
- Science

The GRAPE-4, the world's fastest computer in 1995-1997, has produced some major scientific results through a wide diversity of large-scale simulations in astrophysics. Applications have included planetary formation, the evolution of star clusters and galactic nuclei, and the formation of galaxies and clusters of galaxies. - Article
- Sep 2005
- ASTROBIOLOGY

Assuming that asteroidal and cometary impacts onto Earth can liberate material containing viable microorganisms, we studied the subsequent distribution of the escaping impact ejecta throughout the inner Solar System on time scales of 30,000 years. Our calculations of the delivery rates of this terrestrial material to Mars and Venus, as well as back to Earth, indicate that transport to great heliocentric distances may occur in just a few years and that the departure speed is significant. This material would have been efficiently and quickly dispersed throughout the Solar System. Our study considers the fate of all the ejected mass (not just the slowly moving material), and tabulates impact rates onto Venus and Mars in addition to Earth itself. Expressed as a fraction of the ejected particles, roughly 0.1% and 0.001% of the ejecta particles would have reached Venus and Mars, respectively, in 30,000 years, making the biological seeding of those planets viable if the target planet supported a receptive environment at the time. In terms of possibly safeguarding terrestrial life by allowing its survival in space while our planet cools after a major killing thermal pulse, we show via our 30,000- year integrations that efficient return to Earth continues for this duration. Our calculations indicate that roughly 1% of the launched mass returns to Earth after a major impact regardless of the impactor speed; although a larger mass is ejected following impacts at higher speeds, a smaller fraction of these ejecta is returned. Early bacterial life on Earth could have been safeguarded from any purported impact-induced extinction by temporary refuge in space. - The Digital Orrery has been used to perform an integration of the motion of the outer planets for 845 million years. This integration indicates that the long-term motion of the planet Pluto is chaotic. Nearby trajectories diverge exponentially with an e-folding time of only about 20 million years.
- The evolution of the entire planetary system has been numerically integrated for a time span of nearly 100 million years. This calculation confirms that the evolution of the solar system as a whole is chaotic, with a time scale of exponential divergence of about 4 million years. Additional numerical experiments indicate that the Jovian planet subsystem is chaotic, although some small variations in the model can yield quasiperiodic motion. The motion of Pluto is independently and robustly chaotic.
- We have designed and built the Orrery, a special computer for high-speed high-precision orbital mechanics computations. On the problems the Orrery was designed to solve, it achieves approximately 10 Mflops in about 1 ft3of space while consuming 150 W of power. The specialized parallelarchitecture of the Orrery, which is well matched to orbital mechanics problems, is the key to obtaining such high performance. In this paper we discuss the design, construction, and programming of the Orrery. Copyright © 1985 by The Institute of Electrical and Electronics Engineers, Inc.
- Article
- Dec 1998
- MON NOT R ASTRON SOC

Mixed-variable symplectic integrators exhibit no long-term accumulation of energy error, beyond that owing to round-off, and they are substantially faster than conventional N-body algorithms. This makes them the integrator of choice for many problems in Solar system astronomy. However, in their original formulation, they become inaccurate whenever two bodies approach one another closely. This occurs because the potential energy term for the pair undergoing the encounter becomes comparable to the terms representing the unperturbed motion in the Hamiltonian. The problem can be overcome using a hybrid method, in which the close encounter term is integrated using a conventional integrator, whilst the remaining terms are solved symplectically. In addition, using a simple separable potential technique, the hybrid scheme can be made symplectic even though it incorporates a non-symplectic component.