We combine analytical understanding of resonant dynamics in two-planet systems with machine learning techniques to train a model capable of robustly classifying stability in compact multi-planet systems over long timescales of $10^9$ orbits. Our Stability of Planetary Orbital Configurations Klassifier (SPOCK) predicts stability using physically motivated summary statistics measured in integrations of the first $10^4$ orbits, thus achieving speed-ups of up to $10^5$ over full simulations. This computationally opens up the stability constrained characterization of multi-planet systems. Our model, trained on $\approx 100,000$ three-planet systems sampled at discrete resonances, generalizes both to a sample spanning a continuous period-ratio range, as well as to a large five-planet sample with qualitatively different configurations to our training dataset. Our approach significantly outperforms previous methods based on systems' angular momentum deficit, chaos indicators, and parametrized fits to numerical integrations. We use SPOCK to constrain the free eccentricities between the inner and outer pairs of planets in the Kepler-431 system of three approximately Earth-sized planets to both be below 0.05. Our stability analysis provides significantly stronger eccentricity constraints than currently achievable through either radial velocity or transit duration measurements for small planets, and within a factor of a few of systems that exhibit transit timing variations (TTVs). Given that current exoplanet detection strategies now rarely allow for strong TTV constraints (Hadden et al., 2019), SPOCK enables a powerful complementary method for precisely characterizing compact multi-planet systems. We publicly release SPOCK for community use.

It is common in classical mechanics to encounter systems whose Hamiltonian $H$ is the sum of an often exactly integrable Hamiltonian $H_0$ and a small perturbation $\epsilon H_1$ with $\epsilon\ll1$. Such near-integrability can be exploited to construct particularly accurate operator splitting methods to solve the equations of motion of $H$. However, in many cases, for example in problems related to planetary motion, it is computationally expensive to obtain the exact solution to $H_0$. In this paper we present a new family of embedded operator splitting (EOS) methods which do not use the exact solution to $H_0$, but rather approximate it with yet another, embedded operator splitting method. Our new methods have all the desirable properties of classical methods which solve $H_0$ directly. But in addition they are very easy to implement and in some cases faster. When applied to the problem of planetary motion, our EOS methods have error scalings identical to that of the often used Wisdom-Holman method but do not require a Kepler solver, nor any coordinate transformations, or the allocation of memory. The only two problem specific functions that need to be implemented are the straight-forward kick and drift steps typically used in the standard second order leap-frog method.

Symplectic integrators have made it possible to study the long-term evolution of planetary systems with direct N-body simulations. In this paper we reassess the accuracy of such simulations by running a convergence test on 20Myr integrations of the Solar System using various symplectic integrators. We find that the specific choice of metric for determining a simulation's accuracy is important. Only looking at metrics related to integrals of motions such as the energy error can overestimate the accuracy of a method. As one specific example, we show that symplectic correctors do not improve the accuracy of secular frequencies compared to the standard Wisdom-Holman method without symplectic correctors, despite the fact that the energy error is three orders of magnitudes smaller. We present a framework to trace the origin of this apparent paradox to one term in the shadow Hamiltonian. Specifically, we find a term that leads to negligible contributions to the energy error but introduces non-oscillatory errors that result in artificial periastron precession. This term is the dominant error when determining secular frequencies of the system. We show that higher order symplectic methods such as the Wisdom-Holman method with a modified kernel or the SABAC family of integrators perform significantly better in secularly evolving systems because they remove this specific term.

Symplectic methods, in particular the Wisdom-Holman map, have revolutionized our ability to model the long-term, conservative dynamics of planetary systems. However, many astrophysically important effects are dissipative. The consequences of incorporating such forces into otherwise symplectic schemes is not always clear. We show that moving to a general framework of non-commutative operators (dissipative or not) clarifies many of these questions, and that several important properties of symplectic schemes carry over to the general case. In particular, we show that explicit splitting schemes generically exploit symmetries in the applied external forces which often strongly suppress integration errors. Furthermore, we demonstrate that so-called `symplectic correctors' (which reduce energy errors by orders of magnitude at fixed computational cost) apply equally well to weakly dissipative systems and can thus be more generally thought of as `weak splitting correctors.' Finally, we show that previously advocated approaches of incorporating additional forces into symplectic methods work well for dissipative forces, but give qualitatively wrong answers for conservative but velocity-dependent forces like post-Newtonian corrections. We release REBOUNDx, an open-source C library for incorporating additional effects into REBOUND N-body integrations, together with a convenient Python wrapper. All effects are machine-independent and we provide a binary format that interfaces with the SimulationArchive class in REBOUND to enable the sharing and reproducibility of results. Users can add effects from a list of pre-implemented astrophysical forces, or contribute new ones.

GJ4276 is an M4.0 dwarf star with an inferred Neptune mass planet from radial velocity (RV) observations. We reanalyse the RV data for this system and focus on the possibility of a second, super earth mass, planet. We compute the timescale for fast resonant librations in the eccentricity to be $\sim 2\,000\,\mathrm{days}$. Given that the observations were taken over $700\,\mathrm{days}$, we expect to see the effect of these librations in the observations. We perform a fully dynamical fit to test this hypothesis. Similar to previous results, we determine that the data could be fit by two planets in a 2:1 mean motion resonance. However, we also find solutions near the 5:4 mean motion resonance which are not present when planet-planet interactions are ignored. Using the MEGNO indicator, we analyze the stability of the system and find that our solutions lie in a stable region of parameter space. We also find that though out of resonance solutions are possible, the system favours a configuration which is in a first order mean motion resonance. The existence of mean motion resonances has important implications in many planet formation theories. Although we do not attempt to distinguish between the one and two planet models in this work, in either case, the predicted orbital parameters are interesting enough to merit further study. Future observations should be able to distinguish between the different scenarios within the next $5\,\mathrm{years}$.

Direct N-body simulations and symplectic integrators are effective tools to study the long-term evolution of planetary systems. The Wisdom-Holman (WH) integrator in particular has been used extensively in planetary dynamics as it allows for large timesteps at good accuracy. One can extend the WH method to achieve even higher accuracy using several different approaches. In this paper we survey integrators developed by Wisdom et al. (1996), Laskar & Robutel (2001, and Blanes et al. (2013). Since some of these methods are harder to implement and not as readily available to astronomers compared to the standard WH method, they are not used as often. This is somewhat unfortunate given that in typical simulations it is possible to improve the accuracy by up to six orders of magnitude (!) compared to the standard WH method without the need for any additional force evaluations. To change this, we implement a variety of high order symplectic methods in the freely available N-body integrator REBOUND. In this paper we catalogue these methods, discuss their differences, describe their error scalings, and benchmark their speed using our implementations.

Hybrid symplectic integrators such as MERCURY are widely used to simulate complex dynamical phenomena in planetary dynamics that could otherwise not be investigated. A hybrid integrator achieves high accuracy during close encounters by using a high order integration scheme for the duration of the encounter while otherwise using a standard 2nd order Wisdom-Holman scheme, thereby optimizing both speed and accuracy. In this paper we reassess the criteria for choosing the switching function that determines which parts of the Hamiltonian are integrated with the high order integrator. We show that the original motivation for choosing a polynomial switching function in MERCURY is not correct. We explain the nevertheless excellent performance of the MERCURY integrator and then explore a wide range of different switching functions including an infinitely differentiable function and a Heaviside function. We find that using a Heaviside function leads to a significantly simpler scheme compared to MERCURY, while maintaining the same accuracy in short term simulations.

The number of isolating integrals of motion of the Trappist-1 system - a late M-dwarf orbited by seven Earth-sized planets - was determined numerically, using an adapted version of the correlation dimension method. It was found that over the investigated time-scales of up to 20 000 years the number of isolating integrals of motion is the same as one would find for a system of seven non-interacting planets - despite the fact that the planets in the Trappist-1 system are strongly interacting. Considering perturbed versions of the Trappist-1 system shows that the system may occupy an atypical part of phase-space with high stability. These findings are consistent with earlier studies.

In single star systems like our own Solar system, comets dominate the mass budget of bodies that are ejected into interstellar space, since they form further away and are less tightly bound. However 1I/`Oumuamua, the first interstellar object detected, appears asteroidal in its spectra and in its lack of detectable activity. We argue that the galactic budget of interstellar objects like 1I/`Oumuamua should be dominated by planetesimal material ejected during planet formation in circumbinary systems, rather than in single star systems or widely separated binaries. We further show that in circumbinary systems, rocky bodies should be ejected in comparable numbers to icy ones. This suggests that a substantial fraction of additional interstellar objects discovered in the future should display an active coma. We find that the rocky population, of which 1I/`Oumuamua seems to be a member, should be predominantly sourced from A-type and late B-star binaries.

On February 6th, 2018 SpaceX launched a Tesla Roadster on a Mars-crossing orbit. We perform N-body simulations to determine the fate of the object over the next 15 Myr. The orbital evolution is initially dominated by close encounters with the Earth. While a precise orbit can not be predicted beyond the next several centuries due to these repeated chaotic scatterings, one can reliably predict the long-term outcomes by statistically analyzing a large suite of possible trajectories with slightly perturbed initial conditions. Repeated gravitational scatterings with Earth lead to a random walk. Collisions with the Earth, Venus and the Sun represent primary sinks for the Roadster's orbital evolution. Collisions with Mercury and Mars, or ejections from the Solar System by Jupiter, are highly unlikely. We calculate a dynamical half-life of the Tesla of approximately 15 Myr, with some 22%, 12% and 12% of Roadster orbit realizations impacting the Earth, Venus, and the Sun within one half-life, respectively. Because the eccentricities and inclinations in our ensemble increase over time due to mean-motion and secular resonances, the impact rates with the terrestrial planets decrease beyond a few million years, whereas the impact rate on the Sun remains roughly constant.

Hamiltonian systems such as the gravitational N-body problem have time-reversal symmetry. However, all numerical N-body integration schemes, including symplectic ones, respect this property only approximately. In this paper, we present the new N-body integrator JANUS, for which we achieve exact time-reversal symmetry by combining integer and floating point arithmetic. JANUS is explicit, formally symplectic and satisfies Liouville's theorem exactly. Its order is even and can be adjusted between two and ten. We discuss the implementation ofJANUS and present tests of its accuracy and speed by performing and analyzing long-term integrations of the Solar System. We show that JANUS is fast and accurate enough to tackle a broad class of dynamical problems. We also discuss the practical and philosophical implications of running exactly time-reversible simulations.

Two Jovian-sized planets are orbiting the star HD155358 near exact mean motion resonance (MMR) commensurability. In this work we re-analyze the radial velocity (RV) data previously collected by Robertson et al. (2012). Using a Bayesian framework we construct two models - one that includes and one that excludes gravitational planet-planet interactions (PPI). We find that the orbital parameters from our PPI and noPPI models differ by up to $2\sigma$, with our noPPI model being statistically consistent with previous results. In addition, our new PPI model strongly favours the planets being in MMR while our noPPI model strongly disfavours MMR. We conduct a stability analysis by drawing samples from our PPI model's posterior distribution and simulating them for $10^9$ years, finding that our best-fit values land firmly in a stable region of parameter space. We explore a series of formation models that migrate the planets into their observed MMR. We then use these models to directly fit to the observed RV data, where each model is uniquely parameterized by only three constants describing its migration history. Using a Bayesian framework we find that a number of migration models fit the RV data surprisingly well, with some migration parameters being ruled out. Our analysis shows that planet-planet interactions are important to take into account when modelling observations of multi-planetary systems. The additional information that one can gain from interacting models can help constrain planet migration parameters.

TRAPPIST-1 is a late M-dwarf orbited by seven Earth-sized planets with orbital period ratios near a chain of mean motion resonances. Due to uncertain system parameters, most orbital configurations drawn from the inferred posterior distribution are unstable on short timescales of $\sim$ 0.5 Myr, even when including the eccentricity damping effect of tides. By contrast, we show that most physically plausible resonant configurations generated through disk migration are stable even without tidal dissipation on timescales of at least 50 Myr ($10^{10}$ orbits), an increase of at least two orders of magnitude. This result, together with the remarkable chain of period ratios in the system, provide strong evidence for convergent migration naturally emplacing the system near an equilibrium configuration forced by the resonant chain. We provide an openly available database of physically plausible initial conditions for TRAPPIST-1 generated through parametrized planet-disk interactions, as well as bit-by-bit reproducible N-body integrations over $10^9-10^{10}$ orbits.

The Atacama Large Millimeter/submilimeter Array (ALMA) recently revealed a set of nearly concentric gaps in the protoplanetary disk surrounding the young star HL Tau. If these are carved by forming gas giants, this provides the first set of orbital initial conditions for planets as they emerge from their birth disks. Using N-body integrations, we have followed the evolution of the system for 5 Gyr to explore the possible outcomes. We find that HL Tau initial conditions scaled down to the size of typically observed exoplanet orbits naturally produce several populations in the observed exoplanet sample. First, for a plausible range of planetary masses, we can match the observed eccentricity distribution of dynamically excited radial velocity giant planets with eccentricities $>$ 0.2. Second, we roughly obtain the observed rate of hot Jupiters around FGK stars. Finally, we obtain a large efficiency of planetary ejections of $\approx 2$ per HL Tau-like system, but the small fraction of stars observed to host giant planets makes it hard to match the rate of free-floating planets inferred from microlensing observations. In view of upcoming GAIA results, we also provide predictions for the expected mutual inclination distribution, which is significantly broader than the absolute inclination distributions typically considered by previous studies.

The reproducibility of experiments is one of the main principles of the scientific method. However, numerical N-body experiments, especially those of planetary systems, are currently not reproducible. In the most optimistic scenario, they can only be replicated in an approximate or statistical sense. Even if authors share their full source code and initial conditions, differences in compilers, libraries, operating systems or hardware often lead to qualitatively different results. We provide a new set of easy-to-use, open-source tools that address the above issues, allowing for exact (bit-by-bit) reproducibility of N-body experiments. In addition to generating completely reproducible integrations, we show that our framework also offers novel and innovative ways to analyze these simulations. As an example, we present a high-accuracy integration of the Solar System spanning 10Gyrs, requiring several weeks to run on a modern CPU. In our framework we can not only easily access simulation data at predefined intervals for which we save snapshots, but at any time during the integration. We achieve this by integrating an on-demand reconstructed simulation forward in time from the nearest snapshot. This allows us to extract arbitrary quantities at any point in the saved simulation exactly (bit-by-bit), and within seconds rather than weeks. We believe that the tools we present in this paper offer a new paradigm for how N-body simulations are run, analyzed, and shared across the community.

This chapter explores the physics shared by planetary rings and the various disks that populate the Universe. It begins with an observational overview, ranging from protoplanetary disks to spiral galaxies, and then compares and contrasts these astrophysical disks with the rings of the Solar System. Emphasis is placed on fundamental physics and dynamics, and how research into the two classes of object connects. Topics covered include disk formation, accretion, collisional processes, waves, instabilities, and satellite-disk interactions.

The requirement that planetary systems be dynamically stable is often used to vet new discoveries or set limits on unconstrained masses or orbital elements. This is typically carried out via computationally expensive N-body simulations. We show that characterizing the complicated and multi-dimensional stability boundary of tightly packed systems is amenable to machine learning methods. We find that training an XGBoost machine learning algorithm on physically motivated features yields an accurate classifier of stability in packed systems. On the stability timescale investigated ($10^7$ orbits), it is 3 orders of magnitude faster than direct N-body simulations. Optimized machine learning classifiers for dynamical stability may thus prove useful across the discipline, e.g., to characterize the exoplanet sample discovered by the upcoming Transiting Exoplanet Survey Satellite (TESS). This proof of concept motivates investing computational resources to train algorithms capable of predicting stability over longer timescales and over broader regions of phase space.

The Astropy Project (http://astropy.org) is, in its own words, "a community effort to develop a single core package for Astronomy in Python and foster interoperability between Python astronomy packages." For five years this project has been managed, written, and operated as a grassroots, self-organized, almost entirely volunteer effort while the software is used by the majority of the astronomical community. Despite this, the project has always been and remains to this day effectively unfunded. Further, contributors receive little or no formal recognition for creating and supporting what is now critical software. This paper explores the problem in detail, outlines possible solutions to correct this, and presents a few suggestions on how to address the sustainability of general purpose astronomical software.

First-order variational equations are widely used in N-body simulations to study how nearby trajectories diverge from one another. These allow for efficient and reliable determinations of chaos indicators such as the Maximal Lyapunov characteristic Exponent (MLE) and the Mean Exponential Growth factor of Nearby Orbits (MEGNO). In this paper we lay out the theoretical framework to extend the idea of variational equations to higher order. We explicitly derive the differential equations that govern the evolution of second-order variations in the N-body problem. Going to second order opens the door to new applications, including optimization algorithms that require the first and second derivatives of the solution, like the classical Newton's method. Typically, these methods have faster convergence rates than derivative-free methods. Derivatives are also required for Riemann manifold Langevin and Hamiltonian Monte Carlo methods which provide significantly shorter correlation times than standard methods. Such improved optimization methods can be applied to anything from radial-velocity/transit-timing-variation fitting to spacecraft trajectory optimization to asteroid deflection. We provide an implementation of first and second-order variational equations for the publicly available REBOUND integrator package. Our implementation allows the simultaneous integration of any number of first and second-order variational equations with the high-accuracy IAS15 integrator. We also provide routines to generate consistent and accurate initial conditions without the need for finite differencing.

A number of Kepler planet pairs lie just wide of first-order mean motion resonances (MMRs). Tides have been frequently proposed to explain these pileups, but it is still an ongoing discussion. We contribute to this discussion by calculating an optimistic theoretical estimate on the minimum initial eccentricity required by Kepler planets to explain the current observed spacing, and compliment these calculations with N-body simulations. In particular, we investigate 27 Kepler systems having planets within 6% of the 2:1 MMR, and find that the initial eccentricities required to explain the observed spacings are unreasonable from simple dynamical arguments. Furthermore, our numerical simulations reveal resonant tugging, an effect which conspires against the migration of resonant planets away from the 2:1 MMR, requiring even higher initial eccentricities in order to explain the current Kepler distribution. Overall, we find that tides alone cannot explain planets close to 2:1 MMR, and additional mechanisms are required to explain these systems.

The Kepler mission was launched in 2009 and has discovered thousands of planet candidates. In a recent paper, Esteves et al. (2013) found a periodic signal in the light curves of KOI-13 and HAT-P-7, with a frequency triple the orbital frequency of a transiting planet. We found similar harmonics in many systems with a high occurrence rate. At this time, the origins of the signal are not entirely certain. We look carefully at the possibility of errors being introduced through our data processing routines but conclude that the signal is real. The harmonics on multiples of the orbital frequency are a result of non-sinusoidal periodic signals. We speculate on their origin and generally caution that these harmonics could lead to wrong estimates of planet albedos, beaming mass estimates, and ellipsoidal variations.

We present WHFast, a fast and accurate implementation of a Wisdom-Holman symplectic integrator for long-term orbit integrations of planetary systems. WHFast is significantly faster and conserves energy better than all other Wisdom-Holman integrators tested. We achieve this by significantly improving the Kepler-solver and ensuring numerical stability of coordinate transformations to and from Jacobi coordinates. These refinements allow us to remove the linear secular trend in the energy error that is present in other implementations. For small enough timesteps we achieve Brouwer's law, i.e. the energy error is dominated by an unbiased random walk due to floating-point round-off errors. We implement symplectic correctors up to order eleven that significantly reduce the energy error. We also implement a symplectic tangent map for the variational equations. This allows us to efficiently calculate two widely used chaos indicators the Lyapunov characteristic number (LCN) and the Mean Exponential Growth factor of Nearby Orbits (MEGNO). WHFast is freely available as a flexible C package, as a shared library, and as an easy-to-use python module.

A recent ALMA image revealed several concentric gaps in the protoplanetary disk surrounding the young star HL Tau. We consider the hypothesis that these gaps are carved by planets, and present a general framework for understanding the dynamical stability of such systems over typical disk lifetimes, providing estimates for the maximum planetary masses. We collect these easily evaluated constraints into a workflow that can help guide the design and interpretation of new observational campaigns and numerical simulations of gap opening in such systems. We argue that the locations of resonances should be significantly shifted in massive disks like HL Tau, and that theoretical uncertainties in the exact offset, together with observational errors, imply a large uncertainty in the dynamical state and stability in such disks. This presents an important barrier to using systems like HL Tau as a proxy for the initial conditions following planet formation. An important observational avenue to breaking this degeneracy is to search for eccentric gaps, which could implicate resonantly interacting planets. Unfortunately, massive disks like HL Tau should induce swift pericenter precession that would smear out any such eccentric features of planetary origin. This motivates pushing toward more typical, less massive disks. For a nominal non-resonant model of the HL Tau system with five planets, we find a maximum mass for the outer three bodies of approximately 2 Neptune masses. In a resonant configuration, these planets can reach at least the mass of Saturn. The inner two planets' masses are unconstrained by dynamical stability arguments.

The multi-planetary system HD128311 hosts at least two planets. Its dynamical formation history has been studied extensively in the literature. We reanalyse the latest radial velocity data for this system with the affine-invariant Markov chain Monte Carlo sampler EMCEE. Using the high order integrator IAS15, we perform a fully dynamical fit, allowing the planets to interact during the sampling process. A stability analysis using the MEGNO indicator reveals that the system is located in a stable island of the parameter space. In contrast to a previous study, we find that the system is locked in a 2:1 mean motion resonance. The resonant angle $\varphi_1$ is librating with a libration amplitude of approximately 37{\deg}. The existence of mean motion resonances has important implication for planet formation theories. Our results confirm predictions of models involving planet migration and stochastic forces.

We present IAS15, a 15th-order integrator to simulate gravitational dynamics. The integrator is based on a Gau\ss-Radau quadrature and can handle conservative as well as non-conservative forces. We develop a step-size control that can automatically choose an optimal timestep. The algorithm can handle close encounters and high-eccentricity orbits. The systematic errors are kept well below machine precision and long-term orbit integrations over $10^9$ orbits show that IAS15 is optimal in the sense that it follows Brouwer's law, i.e. the energy error behaves like a random walk. Our tests show that IAS15 is superior to a mixed-variable symplectic integrator (MVS) and other popular integrators, including high-order ones, in both speed and accuracy. In fact, IAS15 preserves the symplecticity of Hamiltonian systems better than the commonly-used nominally symplectic integrators to which we compared it. We provide an open-source implementation of IAS15. The package comes with several easy-to-extend examples involving resonant planetary systems, Kozai-Lidov cycles, close encounters, radiation pressure, quadrupole moment, and generic damping functions that can, among other things, be used to simulate planet-disc interactions. Other non-conservative forces can be added easily.

The detection of strong thermochemical disequilibrium in the atmosphere of an extrasolar planet is thought to be a potential biosignature. In this article we present a new kind of false positive that can mimic a disequilibrium or any other biosignature that involves two chemical species. We consider a scenario where the exoplanet hosts a moon that has its own atmosphere and neither of the atmospheres is in chemical disequilibrium. Our results show that the integrated spectrum of the planet and the moon closely resembles that of a single object in strong chemical disequilibrium. We derive a firm limit on the maximum spectral resolution that can be obtained for both directly-imaged and transiting planets. The spectral resolution of even idealized space-based spectrographs that might be achievable in the next several decades is in general insufficient to break the degeneracy. Both chemical species can only be definitively confirmed in the same object if absorption features of both chemicals can be unambiguously identified and their combined depth exceeds 100%.

We present the Superparticle Model/Algorithm for Collisions in Kuiper belts and debris disks (SMACK), a new method for simultaneously modeling, in 3-D, the collisional and dynamical evolution of planetesimals in a debris disk with planets. SMACK can simulate azimuthal asymmetries and how these asymmetries evolve over time. We show that SMACK is stable to numerical viscosity and numerical heating over 10^7 yr, and that it can reproduce analytic models of disk evolution. We use SMACK to model the evolution of a debris ring containing a planet on an eccentric orbit. Differential precession creates a spiral structure as the ring evolves, but collisions subsequently break up the spiral, leaving a narrower eccentric ring.

The Kepler-36 system consists of two planets that are spaced unusually close together, near the 7:6 mean motion resonance. While it is known that mean motion resonances can easily form by convergent migration, Kepler-36 is an extreme case due to the close spacing and the relatively high planet masses of 4 and 8 times that of the Earth. In this paper, we investigate whether such a system can be obtained by interactions with the protoplanetary disc. These discs are thought to be turbulent and exhibit density fluctuations which might originate from the magneto-rotational instability. We adopt a realistic description for stochastic forces due to these density fluctuations and perform both long term hydrodynamical and N-body simulations. Our results show that planets in the Kepler-36 mass range can be naturally assembled into a closely spaced planetary system for a wide range of migration parameters in a turbulent disc similar to the minimum mass solar nebula. The final orbits of our formation scenarios tend to be Lagrange stable, even though large parts of the parameter space are chaotic and unstable.

We present results from large-scale particle simulations of the viscous overstability in Saturn's rings. The overstability generates a variety of structure on scales covering a few hundred metres to several kilometres, including axisymmetric wavetrains and larger-scale modulations. Such patterns have been observed in Saturn's rings by the Cassini spacecraft. Our simulations model the collisional evolution of particles in a co-rotating patch of the disk. These are the largest N-body simulations of the viscous overstability yet performed. The radial box size is five orders of magnitude larger than a typical particle radius, and so describes a 20-50 km radial portion of the rings. Its evolution is tracked for more than 10,000 orbits. In agreement with hydrodynamics, our N-body simulations reveal that the viscous overstability exhibits a rich set of dynamics characterised by nonlinear travelling waves with wavelengths of a few hundred meters. In addition, wave defects, such as sources and shocks, punctuate this bed of waves and break them up into large-scale divisions of radial width ~5 km. We find that the wavelength of the travelling waves is positively correlated with the mean optical depth. In order to assess the role of the numerical boundary conditions and also background ring structure, we include simulations of broad spreading rings and simulations with a gradient in the background surface density. Overall, our numerical results and approach provide a tool with which to interpret Cassini occultation observations of microstructure in Saturn's rings. We present an example of such a synthetic occultation observation and discuss what features to expect. We also make the entire source code freely available.

I present a new kind of astronomical database based on small text files and a distributed version control system. This encourages the community to work collaboratively. It creates a decentralized, completely open and democratic way of managing small to medium sized heterogeneous astronomical databases and catalogues. The use of the XML file format allows an easy to parse and read, yet dynamic and extendable database structure. The Open Exoplanet Catalogue is based on these principles and presented as an example. It is a catalogue of all discovered extra-solar planets. It is the only catalogue that can correctly represent the orbital structure of planets in arbitrary binary, triple and quadruple star systems, as well as orphan planets.

The Kepler planet candidates are an interesting testbed for planet formation scenarios. We present results from N-body simulations of multi-planetary systems that resemble those observed by Kepler. We add both smooth (Type I/II) and stochastic migration forces. The observed period ratio distribution is inconsistent with either of those two scenarios on its own. However, applying both stochastic and smooth migration forces to the planets simultaneously results in a period ratio distribution that is similar to the observed one. This is a natural scenario if planets form in a turbulent proto-planetary disk where these forces are always present. We show how the observed period ratio and eccentricity distribution can constrain the relative strength of these forces, a parameter which has been notoriously hard to predict for decades. We make the source code of our simulations and the initial conditions freely available to enable the community to expand this study and include effect other than planetary migration.

Kilometer-sized moonlets in Saturn's A ring create S-shaped wakes called "propellers" in surrounding material. The Cassini spacecraft has tracked the motions of propellers for several years and finds that they deviate from Keplerian orbits having constant semimajor axes. The inferred orbital migration is known to switch sign. We show using a statistical test that the time series of orbital longitudes of the propeller Bl\'eriot is consistent with that of a time-integrated Gaussian random walk. That is, Bl\'eriot's observed migration pattern is consistent with being stochastic. We further show, using a combination of analytic estimates and collisional N-body simulations, that stochastic migration of the right magnitude to explain the Cassini observations can be driven by encounters with ring particles 10-20 m in radius. That the local ring mass is concentrated in decameter-sized particles is supported on independent grounds by occultation analyses.

At least two multi-planetary systems in a 4:3 mean motion resonance have been found by radial velocity surveys. These planets are gas giants and the systems are only stable when protected by a resonance. Additionally the Kepler mission has detected at least 4 strong candidate planetary systems with a period ratio close to 4:3. This paper investigates traditional dynamical scenarios for the formation of these systems. We systematically study migration scenarios with both N-body and hydro-dynamic simulations. We investigate scenarios involving the in-situ formation of two planets in resonance. We look at the results from finely tuned planet-planet scattering simulations with gas disk damping. Finally, we investigate a formation scenario involving isolation-mass embryos. Although the combined planet-planet scattering and damping scenario seems promising, none of the above scenarios is successful in forming enough systems in 4:3 resonance with planetary masses similar to the observed ones. This is a negative result but it has important implications for planet formation. Previous studies were successful in forming 2:1 and 3:2 resonances. This is generally believed to be evidence of planet migration. We highlight the main differences between those studies and our failure in forming a 4:3 resonance. We also speculate on more exotic and complicated ideas. These results will guide future investigators toward exploring the above scenarios and alternative mechanisms in a more general framework.

We describe the dynamics of a stream of equally spaced macroscopic particles in orbit around a central body (e.g. a planet or star). A co-orbital configuration of small bodies may be subject to gravitational instability, which takes the system to a spreading, disordered and collisional state. We detail the linear instability's mathematical and physical features using the shearing sheet model and subsequently track its nonlinear evolution with local N-body simulations. This model provides a convenient tool with which to understand the gravitational and collisional dynamics of narrow belts, such as Saturn's F-ring and the streams of material wrenched from tidally disrupted bodies. In particular, we study the tendency of these systems to form long-lived particle aggregates. Finally, we uncover an unexpected connection between the linear dynamics of the gravitational instability and the magnetorotational instability.

We study the interaction of a proto-planetary disk and a planet on a highly inclined orbit in the linear regime. The evolution of the planet is dominated by dynamical friction for planet masses above several Earth-masses. Smaller planets are dominated by aerodynamic drag, especially for very high inclinations and retrograde orbits. The time-scales associated with migration and inclination damping are calculated. For certain values of the inclination, the inclination damping time-scale is longer than the migration time-scale and the disk lifetime. This result shows that highly inclined planets can not (re-)align with the proto-planetary disk. We discuss the dependence of numerical simulations on the gravitational softening parameter. We find only a logarithmic dependence, making global three dimensional simulations of this process computationally feasible. A large fraction of Hot Jupiters is on highly inclined orbits with respect to the rotation axis of the star. On the other hand small-mass planetary systems discovered by the Kepler mission have low mutual inclinations. This shows that there are two distinct formation mechanisms at work. The process that creates inclined Hot Jupiters does not operate on small mass planets because the damping timescales are so long that these systems would still be inclined today.

REBOUND is a new multi-purpose N-body code which is freely available under an open-source license. It was designed for collisional dynamics such as planetary rings but can also solve the classical N-body problem. It is highly modular and can be customized easily to work on a wide variety of different problems in astrophysics and beyond. REBOUND comes with three symplectic integrators: leap-frog, the symplectic epicycle integrator (SEI) and a Wisdom-Holman mapping (WH). It supports open, periodic and shearing-sheet boundary conditions. REBOUND can use a Barnes-Hut tree to calculate both self-gravity and collisions. These modules are fully parallelized with MPI as well as OpenMP. The former makes use of a static domain decomposition and a distributed essential tree. Two new collision detection modules based on a plane-sweep algorithm are also implemented. The performance of the plane-sweep algorithm is superior to a tree code for simulations in which one dimension is much longer than the other two and in simulations which are quasi-two dimensional with less than one million particles. In this work, we discuss the different algorithms implemented in REBOUND, the philosophy behind the code's structure as well as implementation specific details of the different modules. We present results of accuracy and scaling tests which show that the code can run efficiently on both desktop machines and large computing clusters.

The shearing sheet is a model dynamical system that is used to study the small-scale dynamics of astrophysical disks. Numerical simulations of particle trajectories in the shearing sheet usually employ the leapfrog integrator, but this integrator performs poorly because of velocity-dependent (Coriolis) forces. We describe two new integrators for this purpose; both are symplectic, time-reversible and second-order accurate, and can easily be generalized to higher orders. Moreover, both integrators are exact when there are no small-scale forces such as mutual gravitational forces between disk particles. In numerical experiments these integrators have errors that are often several orders of magnitude smaller than competing methods. The first of our new integrators (SEI) is well-suited for disks in which the typical inter-particle separation is large compared to the particles' Hill radii (e.g., planetary rings), and the second (SEKI) is designed for disks in which the particles are on bound orbits or the separation is smaller than the Hill radius (e.g., irregular satellites of the giant planets).

The increasing number of extra-solar planets opens a new opportunity for studies of the formation of planetary systems. Resonant systems are of particular interest because their dynamical configuration provides constraints on the unobservable formation and migration phase. In this thesis, formation scenarios for the planetary systems HD128311 and HD45364 are presented. N-body simulations of two planets and two dimensional hydrodynamical simulations of proto-planetary discs are used to model the migration phase and the capture into resonance. The results indicate that the proto-planetary disc has a larger surface density than previously thought. Proto-planets are exposed to stochastic forces, generated by density fluctuations in the disc. A generic model of both a single planet, and two planets in resonance, being stochastically forced is presented. The system GJ876, for example, is stable for reasonable strengths of the stochastic forces, but systems with lighter planets can get disrupted. Even if they are not disrupted completely, stochastic forces create characteristic, observable libration patterns. Turbulence plays also an important role during the early phases of the planet formation process. Due to the large separation of scales, the gravitational collapse in the core accretion model is very hard to model numerically. A scaled method is presented, that allows for the correct treatment of self-gravity in a marginally collisional system by taking into account the relevant small scale processes. Interestingly, this system is dynamically very similar to Saturn's rings. The stochastic migration of small bodies in Saturn's rings is also studied. Analytic predictions of the interactions of a moonlet/propeller with ring particles are compared to collisional N-body simulations with up to a million particles. The random walk is fast enough to be directly observable by the Cassini spacecraft.

Many small moonlets, creating propeller structures, have been found in Saturn's rings by the Cassini spacecraft. We study the dynamical evolution of such 20-50m sized bodies which are embedded in Saturn's rings. We estimate the importance of various interaction processes with the ring particles on the moonlet's eccentricity and semi-major axis analytically. For low ring surface densities, the main effects on the evolution of the eccentricity and the semi-major axis are found to be due to collisions and the gravitational interaction with particles in the vicinity of the moonlet. For large surface densities, the gravitational interaction with self-gravitating wakes becomes important. We also perform realistic three dimensional, collisional N-body simulations with up to a quarter of a million particles. A new set of pseudo shear periodic boundary conditions is used which reduces the computational costs by an order of magnitude compared to previous studies. Our analytic estimates are confirmed to within a factor of two. On short timescales the evolution is always dominated by stochastic effects caused by collisions and gravitational interaction with self-gravitating ring particles. These result in a random walk of the moonlet's semi-major axis. The eccentricity of the moonlet quickly reaches an equilibrium value due to collisional damping. The average change in semi-major axis of the moonlet after 100 orbital periods is 10-100m. This translates to an offset in the azimuthal direction of several hundred kilometres. We expect that such a shift is easily observable.

Hundred meter sized objects have been identified by the Cassini spacecraft in Saturn's A ring through the so-called "propeller" features they create in the ring. These moonlets should migrate, due to their gravitational interaction with the ring ; in fact, some orbital variation have been detected. The standard theory of type I migration of planets in protoplanetary disks can't be applied to the ring system, as it is pressureless. Thus, we compute the differential torque felt by a moonlet embedded in a two-dimensional disk of solid particles, with flat surface density profile, both analytically and numerically. We find that the corresponding migration rate is too small to explain the observed variations of the propeller's orbit in Saturn's A-ring. However, local density fluctuations (due to gravity wakes in the marginally gravitationally stable A-ring) may exert a stochastic torque on a moonlet. Our simulations show that this torque can be large enough to account for the observations, depending on the parameters of the rings. We find that on time scales of several years the migration of propellers is likely to be dominated by stochastic effects (while the former, non-stochastic migration dominates after ~ 10^{4-5} years). In that case, the migration rates provided by observations so far suggests that the surface density of the A ring should be of the order of 700 kg/m^2. The age of the propellers shouldn't exceed 1 to 100 million years, depending on the dominant migration regime.

The formation mechanism of planetesimals in protoplanetary discs is hotly debated. Currently, the favoured model involves the accumulation of meter-sized objects within a turbulent disc, followed by a phase of gravitational instability. At best one can simulate a few million particles numerically as opposed to the several trillion meter-sized particles expected in a real protoplanetary disc. Therefore, single particles are often used as super-particles to represent a distribution of many smaller particles. It is assumed that small scale phenomena do not play a role and particle collisions are not modeled. The super-particle approximation can only be valid in a collisionless or strongly collisional system, however, in many recent numerical simulations this is not the case. In this work we present new results from numerical simulations of planetesimal formation via gravitational instability. A scaled system is studied that does not require the use of super-particles. We find that the scaled particles can be used to model the initial phases of clumping if the properties of the scaled particles are chosen such that all important timescales in the system are equivalent to what is expected in a real protoplanetary disc. Constraints are given for the number of particles needed in order to achieve numerical convergence. We compare this new method to the standard super-particle approach. We find that the super-particle approach produces unreliable results that depend on artifacts such as the gravitational softening in both the requirement for gravitational collapse and the resulting clump statistics. Our results show that short range interactions (collisions) have to be modelled properly.

The recently discovered planetary system HD45364 which consists of a Jupiter and Saturn mass planet is very likely in a 3:2 mean motion resonance. The standard scenario to form planetary commensurabilities is convergent migration of two planets embedded in a protoplanetary disc. When the planets are initially separated by a period ratio larger than two, convergent migration will most likely lead to a very stable 2:1 resonance for moderate migration rates. To avoid this fate, formation of the planets close enough to prevent this resonance may be proposed. However, such a simultaneous formation of the planets within a small annulus, seems to be very unlikely. Rapid type III migration of the outer planet crossing the 2:1 resonance is one possible way around this problem. In this paper, we investigate this idea in detail. We present an estimate for the required convergent migration rate and confirm this with N-body and hydrodynamical simulations. If the dynamical history of the planetary system had a phase of rapid inward migration that forms a resonant configuration, we predict that the orbital parameters of the two planets are always very similar and hence should show evidence of that. We use the orbital parameters from our simulation to calculate a radial velocity curve and compare it to observations. Our model can explain the observational data as good as the previously reported fit. The eccentricities of both planets are considerably smaller and the libration pattern is different. Within a few years, it will be possible to observe the planet-planet interaction directly and thus distinguish between these different dynamical states.

Many multi-planet systems have been discovered in recent years. Some of them are in mean-motion resonances (MMR). Planet formation theory was successful in explaining the formation of 2:1, 3:1 and other low resonances as a result of convergent migration. However, higher order resonances require high initial orbital eccentricities in order to be formed by this process and these are in general unexpected in a dissipative disk. We present a way of generating large initial eccentricities using additional planets. This procedure allows us to form high order MMRs and predict new planets using a genetic N-body code.

We summarize the analytic model and numerical simulations of stochastically forced planets in a turbulent disk presented in a recent paper by Rein and Papaloizou. We identify two modes of libration in systems with planets in mean motion resonance which react differently to random forces. The slow mode, which mostly corresponds to motion of the angle between the apsidal lines of the two planets, is converted to circulation more readily than the fast mode which is associated with oscillations of the semi-major axes. We therefore conclude that stochastic forcing due to disk turbulence may have played an important role in shaping the configuration of observed systems in mean motion resonance. For example, it naturally provides a mechanism for accounting for the HD128311 system for which the fast mode librates and the slow mode does not.

Aims. We clarify the response of extrasolar planetary systems in a 2:1 mean motion commensurability with masses ranging from the super Jovian range to the terrestrial range to stochastic forcing that could result from protoplanetary disk turbulence. The behaviour of the different libration modes for a wide range of system parameters and stochastic forcing magnitudes is investigated. The growth of libration amplitudes is parameterized as a function of the relevant physical parameters. The results are applied to provide an explanation of the configuration of the HD128311 system. Methods. We first develop an analytic model from first principles without making the assumption that both eccentricities are small. We also perform numerical N-body simulations with additional stochastic forcing terms to represent the effects of putative disk turbulence. Results. Systems are quickly destabilized by large magnitudes of stochastic forcing but some stability is imparted should systems undergo a net orbital migration. The slow mode, which mostly corresponds to motion of the angle between the apsidal lines of the two planets, is converted to circulation more readily than the fast mode which is associated with oscillations of the semi-major axes. This mode is also vulnerable to the attainment of small eccentricities which causes oscillations between periods of libration and circulation. Conclusions. Stochastic forcing due to disk turbulence may have played a role in shaping the configurations of observed systems in mean motion resonance. It naturally provides a mechanism for accounting for the HD128311 system.