Moving Planets Around: An Introduction to N-Body Simulations Applied to Exoplanetary Systems

Moving Planets Around

An Introduction to N-Body Simulations Applied to Exoplanetary Systems

Javier Roa, Adrian S. Hamers, Maxwell X. Cai, and Nathan W. C. Leigh

The MIT Press

Cambridge, Massachusetts

London, England

© 2020 Massachusetts Institute of Technology

This work is subject to a Creative Commons CC BY-NC 4.0 license. Subject to such license, all rights are reserved.

Library of Congress Cataloging-in-Publication Data

Names: Roa, Javier (Javier Roa Vicens), author. | Hamers, Adrian S. (Adrian Sven), 1988- author. | Cai, Maxwell X. (Maxwell Xu), author. | Leigh, Nathan W. C., author.

Title: Moving planets around: an introduction to n-body simulations applied to exoplanetary systems / Javier Roa, Adrian S. Hamers, Maxwell X. Cai, and Nathan W. C. Leigh.

Description: Cambridge, Massachusetts: The MIT Press, [2019] | Includes bibliographical references and index.

Identifiers: LCCN 2019057050 | ISBN 9780262539340 (paperback)

Subjects: LCSH: Extrasolar planets–Orbits–Computer simulation. | Many-body problem.

Classification: LCC QB820.R63 2019 | DDC 523.2/4–dc23

LC record available at https://lccn.loc.gov/2019057050

d_r0

Dr. Javier Roa

Jet Propulsion Laboratory

California Institute of Technology¹

Pasadena, CA 91109, USA

Dr. Maxwell X. Cai

Leiden Observatory

Leiden University

2300 RA Leiden, The Netherlands

SURF Cooperative

1098 XG Amsterdam

The Netherlands

Dr. Adrian S. Hamers

School of Natural Sciences

Institute for Advanced Study

Princeton, NJ 08540, USA

Max Planck Institute for Astrophysics

Karl-Schwarzschild-Str. 1

Garching, D-85741, Germany

Dr. Nathan W. C. Leigh

Departamento de Astronomía

Facultad de Ciencias Físicas y Matemáticas

Universidad de Concepción

Concepción, Chile

Department of Astrophysics

American Museum of Natural History

New York, NY 10024, USA

¹ This work was done as a private venture and not in the author’s capacity as an employee of the Jet Propulsion Laboratory, California Institute of Technology.

Contents

List of Code Snippets

Acknowledgments

Foreword

Prologue

1 A Bestiary of Planets

1.1 Many Moons

1.2 Multiple Star Systems

1.3 Dead Stars

1.4 Supermassive Black Holes

1.5 The Cast of Characters

2 Physics Background

2.1 Newton’s Laws

2.2 Coordinate Systems

2.3 The Two-Body Problem

2.4 The General N-Body Problem

2.5 Multiplanet Systems

3 First Two-Body Code

3.1 Choosing a Programming Language

3.2 Forward Euler Integrator

4 Accuracy and Performance of the Integration

4.1 Structuring the Code

4.2 Speeding Up the Code: Allocating Arrays

4.3 Checking the Conservation Laws

5 Fixed Step-Size Integration

5.1 Truncation and Round-Off Errors

5.2 Runge–Kutta Methods

5.3 Multistep Methods

5.4 Leapfrog Integrator

5.5 Symplectic Integrators

5.6 Numerical Performance

6 Variable Step-Size Integration

6.1 Estimating the Local Truncation Error

6.2 Step-Size Control

6.3 Initial Step Size

6.4 Implementing Flexibility

6.5 Numerical Performance

6.6 Thoughts on Choosing the Integration Method

6.7 Code Review

7 The Three- and N-Body Problems

7.1 From the Two-Body to the Three-Body Problem

7.2 The N-Body Problem

8 Gauss–Radau Integrator of the Fifteenth Order

8.1 Preparing the Second-Order System

8.2 Numerical Quadratures: Choosing the Right Sequence

8.3 Approximating the Integrand Using a Polynomial

8.4 Approximating the Force Function

8.5 Computing the Coefficients

8.6 Radau Sequence and Integration Order

8.7 Coding the Main Components of the Integrator

8.8 Advancing One Integration Step

8.9 The Complete Integrator

8.10 Testing the Integrator

8.11 Code Review

9 Symplectic Map for Long-Term Integration

9.1 Understanding Time Scales

9.2 Long-Term Evolution

9.3 Wisdom–Holman Integrator

9.4 Propagating the Solar System for a Million Years

9.5 Code Review

10 Building a Production Code

10.1 An Object-Oriented Discussion

10.2 Requirement Analysis: What Do We Want Our Code to Be?

10.3 Constructing the Software Framework

10.4 Accelerating the Code Using Native C

11 Defining the Project

11.1 What Project to Tackle?

11.2 Circumbinary Planets

12 Setting Up the Project

12.1 Initial Conditions

12.2 The First Simulation

12.3 Running on a Cluster

13 Running and Analyzing the Simulations

13.1 Problems with the Simulations

13.2 A Suite of Integrations for Kepler-16

13.3 The Other Kepler Systems

14 How to Write a Publishable Research Paper

14.1 Abstract

14.2 Introduction

14.3 Methods

14.4 Results

14.5 Discussion

14.6 Summary

15 Conclusions

15.1 Physics of the N-Body Problem

15.2 Programming Languages

15.3 Numerical Integrators

15.4 Research Project

15.5 Publishing the Results

APPENDICES

A Derivation of Kepler’s Third Law and the Kepler Equation

A.1 The Kepler Equation

A.2 Conserved Quantities in the N-Body Problem

B Keplerian Propagator in Universal Variables

B.1 Stumpff Functions

B.2 Solving the Universal Kepler Equation

B.3 Two-Body Propagator

C Introduction to Matrices

D Derivations in the Lidov–Kozai Problem

D.1 The Three-Body Hamiltonian

D.2 Averaging the Hamiltonian

D.3 The Equations of Motion

D.4 Analytic Solutions in the Test Particle Quadrupole-Order Limit

Bibliography

Index

List of Figures

Figure 2.1

Several vectors, R, R1, and R2, in a two-dimensional plane. R is a vector pointing from the origin to an arbitrary point P, whereas R1 and R2 point to two bodies labeled 1 and 2, respectively. The difference vector r12 points from the endpoint of R2 to R1 where bodies 1 and 2 are located, and its magnitude is the distance between the bodies.

Figure 2.2

The relative position vector r in the orbital plane spanned by the unit vectors ê and . The vector r can be constructed from its projections along ê and .

Figure 2.3

Relative orbits in the two-body problem for different values of the eccentricity e.

Figure 2.4

Illustration of the orbital elements in the (x, y, z)-plane.

Figure 2.5

The evolution of the eccentricities and inclinations of the planets in the Solar System according to Laplace–Lagrange theory.

Figure 3.1

Plot generated by the PYTHON script representing f(x) = sin(x)+ sin(2x).

Figure 3.2

First attempt of the forward Euler integrator, with Δt = 0.1.

Figure 3.3

Second attempt of the forward Euler integrator, with Δt = 0.01.

Figure 3.4

Third and fourth attempts of the forward Euler integrator.

Figure 3.5

Two equal-mass particles (P1 and P2) orbiting their common center of mass.

Figure 4.1

Time evolution of the total energy in the two-body problem.

Figure 4.2

Time evolution of the total angular momentum in the two-body problem.

Figure 4.3

Conservation of the total angular momentum using different time steps.

Figure 4.4

Time evolution of the eccentricity compared to the orbital period. Note the clearly periodic oscillations.

Figure 5.1

Error in energy conservation for different integrators (Δt = 0.01 and tf = 500). The dash-dot line at the bottom represents the machine zero, ε.

Figure 5.2

Performance of different integration schemes. Ef denotes the value of the energy at the end of the integration interval.

Figure 6.1

An example of the trajectories followed by two gravitationally bound equal-mass particles in an eccentric orbit.

Figure 6.2

A comparison of the performances of different variable step-size integrators.

Figure 7.1

The figure-8 orbital configuration.

Figure 7.2

Solution to the Pythagorean three-body problem.

Figure 7.3

Double figure-8 orbits with five bodies.

Figure 7.4

Integration of the Solar System for 1,000 years with RKF4(5) and ε = 10−10.

Figure 8.1

Area under the curve F(x) and position of the five Legendre nodes in table 8.1.

Figure 8.2

Area under the polynomial P(x) and position of five Legendre nodes.

Figure 8.3

Steps for computing the coefficients Bi at every step i.

Figure 8.4

Solution to the elliptic two-body problem using GR15.

Figure 8.5

Solution to the figure-8 problem with N = 5.

Figure 8.6

Relative error of the integration of the Solar System for 1,000 years.

Figure 9.1

Integration step size when propagating the Solar System with Gauss–Radau 15.

Figure 9.2

Time scales associated with the motion of Mercury and Jupiter.

Figure 9.3

Integration step size when propagating the Solar System without Mercury and Venus.

Figure 9.4

Schematic representation of the Wisdom–Holman integrator for two consecutive steps.

Figure 9.5

Graphical representation of the Jacobi coordinates.

Figure 9.6

Solution to the elliptic two-body problem using the WH integrator.

Figure 9.7

Energy error during the integration of the Solar System.

Figure 9.8

Evolution of the Solar System over a million years.

Figure 10.1

The software architecture of ABIE. The top layer is the presentation layer, consisting of a command-line interface (CLI), a simple graphical user interface (GUI) for plotting, and a PYTHON application programming interface (API); the middle layer is the logic layer, consisting of an abstract integrator API and the concrete implementation of various integrator algorithms; the bottom layer is the data layer, consisting of the input/output facility, the in-memory data structure for particles and particle sets, the utilities of orbital element conversions, and so on.

Figure 10.2

The CUDA memory model. Processors are grouped in units of blocks. A GPU processor has exclusive access to its own registers and local memory; the processors in the same block have shared access to the block-wide shared memory, but the shared memory of one block is not accessible from another block. All blocks have read/write access to the global memory, but the constant memory and the texture memory are read-only. Multiple blocks are grouped into a grid (not shown).

Figure 12.1

Derivation of the Hill radius.

Figure 12.2

Orbital elements as a function of time for the Kepler-16 system with am = 2rH ≃ 0.068 AU and ipm = 0°. Black, blue, and red lines correspond to the binary, planet, and moon, respectively. The moon quickly becomes unbound.

Figure 12.3

The stars (black lines), planet (blue line), and the moon (red line) in the (x,z)-plane for Kepler-16, using the same parameters as in Figure 12.2. Left (right)-hand panel: after 10 (1, 000) orbits of the planet around the binary.

Figure 12.4

Orbital elements as a function of time for the Kepler-16 system with am = rH ≃ 0.034 AU and ipm = 0°. The moon still becomes unbound, although it takes a little longer.

Figure 12.5

Orbital elements as a function of time for the Kepler-16 system with am = 0.01 AU and ipm = 0°. The moon now remains bound.

Figure 12.6

The stars (black lines), planet (blue line), and the moon (red line) in the (x,z)-plane for Kepler-16 and the same parameters as in Figure 12.5, that is, am = 0.01 AU and ipm = 0°. Left (right)-hand panel: after 10 (1, 000) orbits of the planet around the binary.

Figure 13.1

Orbital elements as a function of time for the Kepler-16 system (am = 0.005 AU and ipm = 81°). Black, blue, and red lines correspond to the binary, planet, and moon, respectively. The moon becomes highly eccentric and changes its inclination due to Lidov–Kozai oscillations.

Figure 13.2

Stability map for the Kepler-16 system. In each panel, each point represents ten realizations with different initial lunar mean anomaly, which are used to define stable (green filled circles), completely unstable (red crosses), and marginally stable (yellow open circles) systems. Collisions of the moon with the planet and stars are indicated with − and ⋆, respectively, where the color encodes the frequency with respect to the ten realizations of the initial lunar mean anomaly (yellow to red for one to ten). The frequencies of the most common collision type among the mean anomalies are shown. For collisions of the moon with the stars, the large and small ⋆ symbols correspond to collisions with the primary and secondary star, respectively. The three panels in each row correspond to different integration times: 400, 700, and 1,000 Pp, where Pp is the orbital period of the planet around the stellar binary. The mass ratio of the moon to the planet is qm = Mm/Mp = 0.001.

Figure 13.3

Stability maps for the Kepler-16 system. Top row: qm = 0.001; middle row: qm = 0.1; bottom row: qm = 0.001, with the stellar binary replaced by a point mass. The vertical black solid lines show the location of the 1:1 MMC of the moon with the binary (equation 13.8); the vertical black dashed lines show the locations of the 2:1 and 1:2 MMCs. The red dashed lines show the boundary for collisions of the moon with the planet due to LK evolution, computed according to equation (13.7). In the third row, the black dotted lines show the Hill stability boundary for the single-star case according to Grishin et al. (2017).

Figure 13.4

Stability maps for nine Kepler CBP systems, similar to Figure 13.3. The name of the system is indicated above each panel. In all cases, qm = 0.001, and the integration time is 1,000 Pp. The vertical black solid lines show the location of the 1:1 MMC of the moon with the binary (equation 13.8); the vertical black dashed lines show the locations of the 2:1 and 1:2 MMCs. The red dashed lines show the boundary for collisions of the moon with the planet due to LK evolution, computed according to equation (13.7). The vertical blue dashed lines (not visible in all panels) show the Hill radius, equation (12.5).

Figure 13.5

Similar to Figure 13.4, now for the case when the stellar binary is replaced by a point mass. The black dotted lines show the Hill stability boundary according to Grishin et al. (2017).

Figure 13.6

Left panel: largest stable value of am, am, crit, normalized to abin and plotted as a function of Mp/Mbin, where Mbin = M1+M2, for each system in the binary case with qm = 0.001 and for a given value of ipm (refer to the legend). The solid black line shows equation (13.8), the location of the 1:1 MMC. Right panel: am, crit normalized to ap. The black dotted line shows the Hill radius (equation 12.5), and the black solid line shows the estimate in equation (13.12) assuming the CBP is close to the critical boundary for stability.

Figure A.1

An illustration of the unit vectors and in polar coordinates.

Figure A.2

The relation between the true anomaly θ and the eccentric anomaly ℰ for a Kepler orbit. The focal point closest to periapsis is F, and P is a point on the ellipse. P′ is the vertical extension of P on a greater circle with radius a centered at the (geometric) center of the ellipse, C. Q is the projection of P on the major axis of the ellipse.

Figure A.3

Newton–Raphson method.

Figure D.1

Schematic depiction of a three-body system.

Figure D.2

A depiction of the vectors and and their sum . All the vectors lie in the same plane, the plane of the drawing. The invariable plane is perpendicular to this plane, that is, perpendicular to , and to the plane of the drawing.

Figure D.3

Top: the maximum inner orbit eccentricity, e1, max, as a function of θ0 ≡ cos (itot, 0), where itot, 0 is the initial mutual inclination angle, for several values of e1, 0 (0.01, 0.1, 0.2, and 0.5) and ω1, 0 (0°, 45°, and 90°). The black dashed lines show

, the maximum value of θ0 for eccentricity excitation in the limit that ω1, 0 = 90°. Bottom: the (cos (ω1), e1) space for e1, 0 = 0.5. These curves are obtained from conservation of energy and angular momentum, that is, equation (D.109) with θ given by equation (D.108). Solid blue curves correspond to ω1, 0 = 0° and red dashed curves to ω1, 0 = 90°.

List of Tables

Table 5.1

Butcher tableau for explicit Runge–Kutta methods.

Table 8.1

Nodes and weights for the 5-point Gauss–Legendre quadrature.

Table 8.2

Position of the Gauss–Radau nodes over the interval [0, 1] for integration orders 3 to 19 (the integration order is 2n − 1).

Table 8.3

Final position and velocity of the particles in the two-body test problem.

List of Code Snippets

4.1 Function ode_two_body_ first_order: Evaluate the first-order two-body acceleration.

4.2 Function integrate_euler: Forward Euler integrator.

4.3 Function initialize_problem: Define initial conditions for the integration script.

4.4 Function plot_trajectory: Plot the two-body orbits.

4.5 Function compute_energy: Compute the two-body energy.

4.6 Function compute_angular_momentum: Compute the magnitude of the two-body angular momentum.

4.7 Function compute_eccentricity: Compute the eccentricity of body 1 with respect to body 2 in the two-body problem.

5.1 Function integrate_runge_kutta: Runge–Kutta integrator.

5.2 Function integrate_adams_bashforth: Adams–Bashforth integrator.

5.3 Function integrate_leapfrog: Leapfrog integrator.

6.1 Estimating the initial step for variable step-size integrators.

6.2 Function integrate_rk_embedded: Integrate using variable step-size RK integrator.

6.3 Main integration loop for adaptive Runge–Kutta methods.

6.4 Function butcher_tableaus_rk: Butcher tableaus for embedded Runge–Kutta methods.

6.5 Function integrate_rk_embedded: Integrate using an embedded (variable step-size) Runge–Kutta method.

7.1 Function ode_three_body_ first_ order: Differential equations of the three-body problem.

7.2 Initial conditions of the figure-8 three-body orbit.

7.3 Function plot_trajectory: Updated version to handle N bodies.

7.4 Function ode_n_body_ first_order: Differential equations for the N-body problem.

7.5 Initial conditions of the double figure-8 problem with five bodies.

7.6 Initial conditions of the planets in the Solar System relative to the Solar System Barycenter.

8.1 Function ode_n_body_ second_order: Evaluate the gravitational acceleration from N bodies.

8.2 Procedure compute_nodes: Compute the nodes of the Radau sequence.

8.3 Functions approx_pos and approx_vel: Approximate the position and velocity of the particles using equation (8.56).

8.4 Predictor-corrector scheme.

8.5 Function initial_time_step: Estimate the initial size of the time step.

8.6 Initial conditions for an example of the elliptic two-body problem.

8.7 Function radau_spacing: Return the Gauss–Radau sequence for integration order 15.

8.8 Function compute_bs_from_gs: Compute Bi from Gj up to ih using equation (8.41).

8.9 Function compute_cs: Precompute the coefficients cij in equation (8.47).

8.10 Function compute_gs: Compute Gi up to ih using equation (8.43). The acceleration Fi is stored in ddys[i,:].

8.11 Function compute_rs: Return the values of the coefficients rij defined in equation (8.44).

8.12 Function refine_bs: Refine the coefficients Bi using equation (8.60).

8.13 Function gauss_radau15_step: Advance one integration step using the Gauss–Radau integrator.

8.14 Function integrate_gauss_radau15: Gauss–Radau integrator of the fifteenth order.

9.1 Function wh_advance_step: Advance one step using the Wisdom–Holman integrator.

9.2 Function wh_kick: Compute the velocity kick.

9.3 Function wh_drift: Drift the state of all bodies (Keplerian propagation).

9.4 Function jacobi2cart: Convert from Jacobi to Cartesian coordinates following Rein and Tamayo (2015).

9.5 Function cart2jacobi: Convert from Cartesian to Jacobi coordinates using the transformation from Rein and Tamayo (2015).

9.6 Function compute_accel: Compute perturbing acceleration in Jacobi coordinates.

9.7 Function propagate_wh: Propagate N-body system using the Wisdom– Holman integrator.

10.1 First version of the Particle class.

10.2 Second version of the Particle class.

10.3 First version of the adding/removing particle methods.

10.4 Second version of the adding/removing particle methods.

10.5 Third version of the adding/removing particles method.

10.6 Minimum code for writing an array into a HDF5 file.

10.7 Minimum code for reading an array into a HDF5 file.

10.8 Output buffer.

10.9 First version of the Integrator API.

10.10 integrator_euler.py (abridged).

10.11 integrator_leap_frog.py (abridged).

10.12 Second version of the Integrator API.

10.13 Third version of the Integrator API.

10.14 Obtaining the command-line arguments from the main function.

10.15 Command-line argument parsing with PYTHON.

10.16 A minimum example of running ABIE programmatically.

10.17 Function fibonacci.c: Compute the Fibonacci sequence.

10.18 test_ctypes.py: Call a C function from PYTHON using ctypes.

10.19 Minimum performance benchmark code abie_benchmark.py.

10.20 Integrate a simplified Solar System with 1,000 Oort cloud particles.

10.21 A simple CUDA kernel of adding two vectors.

10.22 A naive implementation of gravity calculation on GPU.

10.23 The driver of the GPU kernel.

B.1 Function stumpff_functions: Evaluate the first four Stumpff functions.

B.2 Function propagate_kepler: Analytic propagation of Kepler’s problem in universal variables.

Acknowledgments

This book would not have been possible without the inspiration and support by Prof. Piet Hut (IAS) and Prof. Jun Makino (RIKEN). We are extremely grateful to them for getting us involved in this endeavor and for supporting us along the way by providing feedback on important decisions and even inviting us to RIKEN in Kobe, Japan, to meet in person and to work on the book. We thank Yuko Wakamatsu, Yoshie Yamaguchi, and the other support staff at RIKEN for taking care of many practical details of our unforgettable trip to Japan.

We also thank Piet Hut for writing a foreword in which he puts our new work in perspective of the history of the development of astrophysical computer software and its educational aspects.

Lastly, we thank Prof. Scott Tremaine (IAS) and Prof. Simon Portegies Zwart (Leiden University) for carefully reading our manuscript and for writing statements of recommendation.

Foreword

This book is unique in providing all that you need to start writing your own software for simulating the dynamical evolution of planets, stars, or other heavenly bodies. For decades, the knowledge gathered here between two covers could only be found scattered across the literature, at best, while partly existing only in the form of oral traditions, passed on from teacher to student or just from student to student.

The authors have asked me to write this foreword in order to provide a bit of prehistory, spanning more than three decades, leading up to the writing of this book. It started in 1985, at the Institute for Advanced Study in Princeton, when Joshua Barnes and I were developing new algorithms for stellar dynamics, such as the tree code that could efficiently speed up gravitational interactions from 𝒪(N²) to 𝒪(N log N).¹

At that time, both of us already had years of experience writing, running, and analyzing the results of N-body codes. In conversations, we both deplored the fact that there were no good books or even review articles about developing such codes. As a result, many users (then and even now!) relied on codes that were at least partly treated like black boxes. This carries a double risk: at the side of the code writer or modifier, there may be latent bugs in the code that had not shown up before and cannot easily be traced, and at the side of the code user, it is all too easy to use a black-boxed code outside the range of applicability that the writer had foreseen.

In addition, Josh and I were frustrated by the clumsy way in which N-body codes were used back then. There was no standard, well-documented software environment in which to develop and test new codes or modify and improve old ones. As a result, just about any astrophysicist working in stellar dynamics had cobbled together a grab bag of programs and scripts to run their own codes, with very little compatibility among them.

In order to address both problems, we developed our own toolbox of small codes that could be combined, using the modularity of the Unix operating system. We used the Unix pipes system to let the user choose a module for setting up initial conditions, from which the data would flow to an integrator that would be used for the actual run, and from there on to one or more analysis programs. A popular article about our approach was published soon afterward in Scientific American by Gerald Sussman and me.²

We quickly attracted other users for our toolbox-based environment. During the next few years, our software forked into a more focused branch and a more general branch with wrappers containing programs provided by others. Josh continued to develop the first branch, as the original version, Zeno.³ Peter Teuben, who had later joined our efforts, maintained a publicly supported version, Nemo.⁴

Meanwhile, I started to design a third version, Starlab, with a different underlying data structure. While Zeno and Nemo were very efficient for conducting experiments in galactic dynamics, through the use of homogeneous arrays to store particle data, they lacked the ability to simulate small-N systems such as star clusters. In order to model the role of binary stars, for example, a more flexible data format was called for, allowing individual stars to be treated separately, each in a different way.

I also used this extra freedom to develop a flexible way of reporting interesting events, through what I called a story mechanism, in which stars undergoing unusual interactions could scribble their own reports. These were implemented as ascii strings, with data structures that could be designed on the fly, as a kind of forerunner of scripting languages, all written in C. In this way, I laid the groundwork for Starlab, with inspiration from helpful conversations with Jun Makino, whom I was visiting during my first sabbatical, at Tokyo University.

After returning to Princeton, Steve McMillan joined me in the development, and soon afterward we met Simon Portegies Zwart, who offered his stellar evolution code to be integrated with the dynamics part of the new environment. By that time, we switched to C++, mainly to use the convenience of data encapsulation, from vector notation, to carrying stellar evolution data around. With Jun Makino joining us as well, the four of us spent a decade to extend Starlab into a workhorse for simulating star cluster evolution.

In the spring of 2001, Jun Makino and I began to develop a vision for a simpler and more educational version of Starlab while hiking in the desert near Tucson before attending one of the first astrophysics conferences on virtual observatories. We worked on two parallel approaches: first to write an introduction to stellar dynamics, from scratch, and second to write an extension in which we introduced very simple, almost cartoon-like versions of stellar evolution.

In the fall of 2001, we started to work on the introduction, which we finished in 2003 as an online book Moving Stars Around, written in C++, followed a few years later by a version with the same title, written in RUBY, a scripting language much like PYTHON.⁵

In the summer of 2002, we produced our first cartoon version of stellar evolution, during the inaugural MODEST meeting in the Museum for Natural History in New York City.⁶ MODEST stands for MOdeling DEnse STellar systems.⁷

A central theme of the MODEST initiative was to develop a multipurpose software environment for astrophysical applications in which different existing numerical codes are incorporated into a single framework. Using STARLAB, compiled as a single code, as our taking-off point, the first steps toward such an environment were taken in Amsterdam in 2006, at the MODEST-6d meeting. The environment was called MUSE, for MUltiscale MUltiphysics Simulation Environment.

MUSE, in turn, led to AMUSE, for Astrophysical Multipurpose Software Environment, developed at Leiden Observatory by Simon Portegies Zwart and coworkers, supported by a multimillion Euro grant from the Dutch government. The AMUSE kickoff workshop was held there in the fall of 2009, followed by various other meetings. Further development remains ongoing, and the package is well maintained, with a friendly staff happy to help users who run into problems.⁸

A very clear and helpful introduction to AMUSE can be found in a recent book.⁹ In many ways, their book is complementary to the current book. Where this book teaches you to write and understand individual codes in detail, the AMUSE book teaches you to set up larger experiments by using individual codes as gray or black boxes that the user can combine by using the tools in the AMUSE environment.

Fast-forwarding to the fall of 2016, I had some conversations about few-body dynamics with Adrian Hamers and Javier Roa at IAS in Princeton. When we talked about various aspects of software development, I mentioned my work with Jun Makino on an ambitious project called the Art of Computational Science, ACS for short, that Jun and I had started in 2003.¹⁰ I added that we had lost momentum after a few years, for various reasons, although our website was still being used by many visitors, either for the codes presented there or the educational texts that we provided.

I especially mentioned the Moving Stars Around books, in their C++ and RUBY versions, mentioned above. I also made a casual remark that it would be nice if someone would write something similar, in a more modern way, and actually carry it to completion in the form of a published book. To my happy surprise, they showed great interest in doing so and immediately started scheming as to how that could be done.

When I next talked with Jun Makino, he invited Adrian and Javier and me to visit him in Kobe, in the summer of 2017, to discuss together what such a newer version of our project might look like. By the time we all met, Maxwell Cai and Nathan Leigh had joined as well, and the six of us had a great time discussing the options for a new book.

Jun and I strongly felt that it would be best for them to start from scratch and only use the few ingredients of our original online book that they found useful. We suggested the use of a dialogue as a pedagogical device and also the inclusion of a discussion about debugging, something that is typically left out in books about software. But other than that, all the ideas presented in the current book are those of the authors.

Indeed, the current book has incorporated many great improvements over the original text. Instead of stars, the particles now represent exoplanets, a very timely and important topic of both observational and theoretical studies. Instead of RUBY, the scripting language used is PYTHON, by now the obvious choice for a widespread and well-supported language, suitable for numerical work. The actual integrators presented and developed in the book from scratch also go well beyond the much more elementary algorithms used in the original version and are among the best that are currently used at the cutting edge of research using simulations of exoplanet dynamics.

Most important, the current book includes a realistic student research project, leading to an actual publication in the scientific literature; see the last five chapters of the book. This was a great idea, since there is no better way to teach something than to show a full example of how to actually do something, from beginning to end. I was stunned when they handed me the draft of their manuscript, since I had never seen anything like that in a textbook.

When I asked them whether they were thinking of publishing their research article and, if so, whether they would do so under the names of the protagonists in the book, Alice and Bob, they smiled. They told me they had already submitted the paper¹¹ but after some discussion had decided to use their own names, in order not to create confusion about the identity of the mysterious Alice and Bob authors.

So in this book, you will find an abundance of original work, from a how-to manual of designing and implementing state-of-the-art simulation codes, up to and including a how-to manual for publishing your first article. The proof of the pudding may be in the eating, but the proof of a piece of research is in the publishing—for real students reading this book as well as for the fictional students who set the example in the last few chapters of this book.

When Jun Makino and I invited the authors to visit us in Kobe two years ago, we never expected to see such a great book growing so quickly out of the small seed that we handed them to get started. Congratulations, Javier, Adrian, Maxwell, and Nathan!

Piet Hut

Institute for Advanced Study

Princeton, NJ

March 2019

¹ Barnes, J. & Hut, P., 1986, Nature 324, 446–449.

² Hut, P. & Sussman, G. J., 1987, Scientific American 255, 145–153.

³http://www.ifa.hawaii.edu/~barnes/zeno/index.html

⁴https://bima.astro.umd.edu/nemo/index_old.html

⁵http://www.artcompsci.org/#msa

⁶ Hut, P., Shara, M. M., Aarseth, S. J., Klessen, R. S., Lombardi, J. C., Makino, J., McMillan, S.,