Implementation of MVS integrators to cosmological N-body codes

1. Implementation of MVS integrators to cosmological N-body codes Ryuji Morishima (UCLA/JPL)

2. N-body codes for planet formation • N-body code: Gravity solver + Integrator • Gravity solver must be fast and handle close encounters • Special hardware (N2): Grape, GPS • Tree (N log N): PKDGRAV, Gadget • Integrator must take a large time step with good accuracy • Bulirsch-Stoer • Hermite: often used with Grape • Mixed Variable Symplectic (MVS) integrators: SYMBA, Mercury • This work: Implementation of SYMBA to PKDGRAV

3. Cosmological N-bodycode PKDGRAV • Developed by Stadel (2001) • Source is open in The astro-code wiki • Tree gravity: 4th order multiple moments • Adaptive to various parallel environments (shared memory, mpi) • Different functions and integrators • Collisions (Richardson et al. 2000) • SPH (Wadsleyet al. 2004) • Fragmentation (Leinhardt & Richardson 2005) • SYMBA integrator (Morishima et al. 2010)

4. Multi-pole expansion Up to 4th order (Hexadecapole) Error estimation from cosmological simulations

5. Tree build k-D Tree Spatial binary tree Spatial binary tree can reduce the higher order multi-pole moments It is also as efficient as k-D tree in neighboring search.

6. Parallel architecture • Mater layer • Controls overall flows of program call • Machine dependent layers • Interface to parallel primitives (e.g. MPI) • Processor Set Tree (PST) layer • Assigns tasks to processors call • Parallel KD layers • Executes tasks in each core One needs to understand PST format but not parallel primitives such as MPI

7. Mixed variable symplectic integrator • Specialized for systems with a massive central body • Mixed variables: Cartesian and Keplerian co-ordinates • A large time step along Keplerian orbit • Time-reversible (no secular error) • Handling close encounters: • SYMBA (Duncan et al. 1998) • Mercury (Chambers 1999) • Most of N-body simulations for planet formation have been performed by these two codes in last decade • But both codes use N2 gravity calculations

8. SYMBA (Symplectic Massive Body Algorithm) • Democratic co-ordinate • Heliocentric position +barycentric velocity Hkep>>Hint (if there is no close encounter) Hkep>>Hsun

9. Multiple time step in close encounters • Potential (or Force) decomposition based on mutual distance normalized by the Hill radius • A higher order potential component is calculated with a small block-sized time step Kepler Drift with F1 Kick (F0) Kick (F0)

10. For time reversibility… • The time step size needs to be determined by the minimum mutual distance during particle drift. • This distance must be estimated by using particle co-ordinates at the beginning and ending of particle drift symmetrically (e.g. Hut et al. 1995).

11. Loop for PKDGRAV-SYMBA • Half kick (t0/2) due to Sun’s motion • Half Kick (t0/2) due to force F0 from other particles • Kepler drift (t0) for all particles • Tree build and neighboring search (after drift) • Particles in close encounters • Sent back topre-drift positions and velocities • Put into a single core (domain decomposition) • SYMBA multiple time stepping • Collisions are also handled here • Tree build and gravity calculation and neighboring search (before drift) • 2 and 1

12. Example (Morishima et al. 2013, EPSL 366) With gas Without gas

13. Evolution of De Martian meteorites Accretion truncated at 14 My Accretion Extrapolated

14. Energy conservation (Mtot =3 MEarth, N = 100, a = 1.0-1.3 AU, Dt = 6 days)

15. Scaling with particle number

16. Multi-body encounter and chain of encounters • In 3-body encounter, time stepping of these 3 bodies is synchronized (for time symmetry) • If the system’s number density is high, all particles share the time step….

17. Summary • PKDGRAV-SYMBA works as desired unless the system’s number density is so high that most of bodies are in multi-body encounters. • For such systems, time symmetry needs to be sacrificed?