Integration Algorithms

Introduction

Particles can be divided into three basic groups:

1. HWPs (Heavy Weight Particles), which always includes the dominant body, whose self- and external-interactions are always accounted for.
2. LWPs (Light Weight Particles), whose interactions with HWPs (and, optionally, ZWPs) is computed, but whose self-potential is always ignored.
3. ZWPs (Zero Weight Particles), whose self- and external-interactions are always neglected—i.e., the particles are treated as massless. They are always perturbed by HWPs, but may or may not be perturbed by LWPs (user-selectable).

This logical division of particles is useful whenever M_H >> M_L >> M_Z, where M_k is the total mass of particles in group k. If used to follow close encounters between bodies, the package is able to shift particles between the LWP and ZWP classes to maintain a balance between accuracy and efficiency (note: this feature is currently disabled).

For background information on the symplectic algorithms employed, see:

• Wisdom & Holman (1991), AJ, 102, 1520
• Touma & Wisdom (1994), AJ, 107, 1189
• Saha & Tremaine (1994), AJ, 108, 1962
• Wisdom, Holman, & Touma (1996), Fields Institute Comm., 10, 217
• Duncan, Levison, & Lee (1998), AJ, 116, 2067
• Rauch & Holman (1999), AJ, 117, 1087
• Rauch & Hamilton (2005), AJ, in preparation

Runge-Kutta and Bulirsh-Stoer ODE integrators are also available; they have been included mainly for the purposes of testing and generating comparison runs.

Integration Coordinate Systems

The HWP/LWP system integration coordinates can be either Jacobi (barycentric) or canonical bodycentric (a.k.a. heliocentric, in the Solar System case). ZWPs are integrated as perturbed two-body orbits, and perturb nothing themselves; their integration coordinates can be either barycentric (specifically, the barycenter of the HWP/LWP subsystem) or canonical bodycentric, and their coordinates (as well as the timestep used) need not agree with those used for the other two groups. HWPs and LWPs always share the same coordinate choice and timestep.

Bodycentric, barycentric, and (time) regularized barycentric coordinates are available for ZWPs; bodycentric and Jacobi coordinates are available for the HWP/LWPs. An optional Modified bodycentric coordinate choice (see IntegCoord) implements Touma & Wisdom's form of the Hamiltonian splitting, which is exact for the two-body problem (unlike canonical bodycentric coordinates).

The default ordering of Jacobi indices i in terms of particle array indices j is (nH=# HWPs, nL=# LWPs):

      i: 0  1       ... nL,  1+nL ... nH-1+nL
      j: 0, nL+nH-1 ... nH,  1    ... nH-1

Thus, array indices of LWPs need not correspond to semi-major axes in any way since all LWPs will be integrated using near-bodycentric coordinates; HWPs should normally appear in order of increasing semi-major axis (this maximizes accuracy). The ordering of Jacobi indices can be changed by specifying a JacIndex field in InputOrder.

To the base coordinates a number of predefined modifiers can be added specifying things such as a 4th order integration or use of a symplectic corrector. Current options:

• base coordinates: Jacobi (HWP/LWPs only) Barycentric (ZWPs only), Bodycentric, RegJacobi, RegBarycentric, ModBodycentric
• splitting order: Kick_Drift (default), Drift_Kick, Shift_Kick
• integrator order: Order2 (default), Order4, Order4mk
• corrector order: Uncorrected (default), Corrected2, Corrected4, Corrected6, Corrected (use highest available order),
• gravitation type: Newtonian (default), PostNewtonian

Example: coord=Corrected+Order4+Jacobi

Let K = `kick', S = `shift', and D = `drift'. By default integrations proceed as K/2 - S/2 - D - S/2 - K/2 (Kick_Drift); option Drift_Kick implies D/2 - S/2 - K - S/2 - D/2 and option Shift_Drift implies S/2 - K/2 - D - K/2 - S/2. In practice hnb_step() proceeds as S/2 - K - S/2 - D, S/2 - K - S/2 - D, or K/2 - S - K/2 - D, respectively; hnb_halfstep() performs the realignment to "half-step" output phase. Note that the shift term is absent unless either bodycentric coordinates or post-Newtonian corrections are used.

Currently, resolving close encounters between particles requires use of one of the ODE integration schemes (see Integrator). Support for close encounters in symplectic integrations will be available in a future release (see Problems and Limitations).