Back to Starlab home

# The Kira Integrator

Kira is Starlab's N-body integration program. Its basic function is to take an input (dyn or hdyn) N-body snapshot and evolve it forward for a specified period of time, producing snapshot and diagnostic (``log'') output at regular intervals. In addition to strictly dynamical evolution of stars and multiple stellar systems, kira also incorporates stellar and binary evolution (via the SeBa subpackage), and the possible influence of an external (``tidal'') gravitational field. The following discussion applies equally well to the GRAPE version of the program, kira_harp3.

Basic Operation

Particle motion is followed using a fourth-order, individual-timestep ``Hermite'' predictor-corrector scheme (Makino and Aarseth 1992). Briefly, during a time step, particle positions and velocities are first predicted to fourth order using the known acceleration and ``jerk'' (time derivative of the acceleration), the new acceleration and jerk are computed, and the motion is then corrected using the additional derivative information thereby obtained. Schematically, if pos, vel, acc, and jerk represent position, velocity, acceleration, and jerk at some time t, and the desired time step is dt, prediction consists of

```        pred_pos = pos + dt * (vel + (dt/2) * (acc + (dt/3) * jerk));
pred_vel = vel + dt * (acc + (dt/2) * jerk);
```
The acceleration and jerk are then recalculated at the predicted time, using pred_pos and pred_vel. Finally, a correction is based on the estimated higher-order derivatives:
```        a3 =  2 * (acc - pred_acc) + dt * (jerk + pred_jerk);
a2 = -3 * (acc - pred_acc) - dt * (2 * jerk + pred_jerk);
```
where
```        a2  =  a'' * dt*dt / 2        (a = acc, a' = jerk)
a3  =  a''' * dt*dt*dt / 6
```
so
```        pos = pred_pos + (a3/20 + a2/12) * dt * dt;
vel = pred_vel + (a3/4  + a2/3) * dt;
```

A single integration step in kira thus proceeds as follows:

• Determine which stars are to be updated next. Each star has associated with it an individual ``base time'' time, representing the time at which it was last advanced, and an individual timestep. The list of stars to be integrated consists of those with the least value of time+timestep. Timesteps are constrained to be powers of 2, allowing ``blocks'' of many stars to be advanced simultaneously.

• Before the step is taken, check for
• system reinitialization
• log output
• escaper removal
• termination of the run
• snapshot output

• Perform low-order prediction of all particles to the new time. This operation may be performed on the GRAPE, if present.

• Recompute the acceleration and jerk on all stars in the current block (using the GRAPE, if available), and correct their postions and velocities for fourth-order accuracy.

• Check for and initiate unperturbed motion.

• Check for collisions and mergers.

• Check for tree reorganization (see below).

• Check for and apply stellar and/or binary evolution, and correct the dynamics as necessary. A more detailed discussion of Starlab's stellar and binary evolution packages may be found in the SeBa section.

Log Output

The format of a snapshot is described in the internals section. The standard log output presently consists of

• information on bulk parameters of the system: total mass, energy, momentum, anisotropy, etc.
• technical information on CPU time, timestep distribution, etc.
• detailed information on the cluster mass distribution: core properties, Lagrangian radii, etc.
• stellar mass distribution and anisotropy by Lagrangian zone
• luminosity profile, and mass and luminosity functions
• cluster stellar content (by spectral type and luminosity class)
• detailed dynamical and physical data on all binary systems.
Additional data are constantly being added as circumstances warrant. The Starlab tool sys_stats may be used to generate similar information for any snapshot. As with all Starlab programs, log output is sent to stderr, allowing stdout to be used exclusively for snapshot information. Mechanisms exist within kira to allow the user to add her own analysis tools to the log output, either simply by modifying the Starlab code (found in .../hdyn/evolve/kira_log.C), or by linking user-written modules with kira.

Tree Structure

As discussed in more detail in the structure section, an N-body system in Starlab is represented as a linked-list structure, in the form of a mainly ``flat'' tree, individual stars being the ``leaves.'' The tree is flat in the sense that single stars (i.e. stars that are not members of any multiple system) are all represented by top-level nodes, having the root node as parent. Binary, triple, and more complex multiple systems are represented as binary trees below their top-level center of mass nodes. The tree structure determines both how node dynamics is implemented and how the long-range gravitational force is computed.

The motion of every node relative to its parent node is followed using the Hermite predictor-corrector scheme described above. The use of relative coordinates at every level ensures that high numerical precision is maintained at all times, even during very close encounters.

The tree evolves dynamically according to simple heuristic rules: particles that approach ``too close'' to one another are combined into a center of mass and binary node; and when a node becomes ``too large'' it is split into its binary components. These rules apply at all levels of the tree structure, allowing arbitrarily complex systems to be followed. The terms ``too close'' and ``too large'' are defined by command-line variables, as described below. The tree rearrangement corresponding to a simple three-body ``exchange'' encounter is:

A four-body exchange/disruption might look like:
and so on.

How the acceleration (and jerk) on a particle or node is computed depends on its location in the tree. Top-level nodes feel the force due to all other top-level nodes in the system. Forces are computed using direct summation over all other particles in the system; no tree or neighbor-list constructs are used. (The integrator is designed specifically to allow efficient computation of these forces using GRAPE hardware, if available.) Nearby binary and multiple systems may be resolved into their components, as necessary.

The internal motion of a binary component is naturally decomposed into two parts: (1) the dominant contribution due to its companion, and (2) the perturbative influence of the rest of the system. (Again, this decomposition is applied recursively, at all levels in a multiple system.) Since the perturbation (gamma) drops off rapidly with distance from the binary center of mass, in typical cases only a few near neighbors are significant perturbers of even a moderately hard binary. These neighbors are most efficiently handled by maintaining lists of perturbers for each binary, recomputed at each center of mass step, thereby greatly reducing the computational cost of the perturbation calculation:

A further efficiency measure is the imposition of unperturbed motion for binaries whose perturbation falls below some specified value. Unperturbed binaries may be followed analytically for many orbits as strictly two-body motion; they are also treated as point masses, from the point of view of their influence on other stars. Because unperturbed binaries are followed in time steps that are integer multiples of the orbit period, we can relax the perturbation threshold for unperturbed motion, relative to that for a perturbed step, as illustrated in the diagram. Perturbed binaries are resolved into their components, both for purposes of determining their center of mass motion and for determining their effect on other stars. Unperturbed treatments of multiple systems also are used, based on empirical studies of the stability of their internal motion.

Escaper Removal

Stars are removed (``stripped'') from the system when they exceed a specified distance from the cluster center of mass (or density center). For systems without an imposed Galactic tidal field, this stripping radius is arbitrary. For systems with a tidal field, the stripping radius is more usually tied to the Jacobi radius of the cluster. The Jacobi radius (governed by the imposed tidal parameters) may or may not have anything to do with the ``King'' radius -- the cutoff radius of an initial King model. Typical relationships between these radii are shown below:

Command-line Parameters

The behavior of kira is controlled by its command-line parameters. As with all Starlab tools, a list of parameters may be obtained by typing

```        kira  --help
```
The (excessively long) result is
```    Starlab version 3.4
program created ...
source file     ...

kira:  Hermite N-body integrator with evolving hierarchical tree
structure, stellar and binary evolution, and an external
tidal field.

Basic options are as follows.  Default values for a new simulation are
indicated in square brackets.  For restart (continuation of a previous
kira calculation), the defaults for many options [*] are determined from
the input snapshot, making it possible to continue a run without having
to re-specify all the command-line parameters previously used.  (These
new defaults may still be overridden on the command line, of course.)
Some options typically having to do with initial setup may be overridden
by data from the input snapshot (if present), as noted below.  Kira may
also turn *on* some options (the B, G, Q,  S, and u settings) if they
were turned on in the previous run, but are not specified on the current
command line.  To prevent this, use the "-o" switch.

The first page of log output gives complete information on parameter
of doubt, read the log file!

Options:     -a    specify accuracy parameter [0.1][*]
-B    turn on binary evolution [off][*]
-c    include comment [none]
-C    specify GRAPE release interval, in seconds [15]
-d    specify log output interval [0.25][*]
-D    specify snapshot interval [end of run]
-e    specify softening length [0][*]
-E    use exact calculation [false]
-f    specify close-encounter distance [1 --> 1/N][*]
-F    specify tidal field type (0 = none, 1 = point-mass,
2 = isothermal, 3 = disk)
[0 if -Q not set; 1 otherwise][*]
-g    specify hysteresis factor [2.5][*]
-G    specify initial stripping radius [none][*]
-h    specify stellar-evolution time step [0.015625 = 1/64][*]
-I    specify (re)initialization timescale [1][*]
-k    specify perturbation factor [1.e-7][*]
-L    specify CPU time limit, in seconds [none]
-M    specify initial total mass, in solar masses
[none; take from input snap if present]
-n    stop at specified number of particles [10]
-N    specify frequency of CPU check output [25000]
-O    save (and overwrite) extra snapshot at each output [no]
-q    specify initial virial ratio [0.5]
-o    prevent kira from overriding some settings (BGQSu)
based on input snapshot data [allow]
-Q    use tidal field [none][*]
-r    specify initial virial radius in code units
[1; take from input snap if present]
-R    specify initial virial radius, in pc
[none; take from input snap if present]
-s    specify random seed [take from system clock]
-S    turn on stellar evolution [off][*]
-t    specify time span of calculation [10]
-T    specify initial virial time scale, in Myr
[none; take from input snap if present]
-u    toggle unperturbed multiple motion [disabled][*]
-U    toggle all unperturbed motion [enabled][*]
-v    toggle "verbose" mode [on]
-X    specify escaper removal timescale [reinit][*]
-y    specify stellar encounter criterion
[0 N-body units or solar radii][*]
-z    specify stellar merger criterion [0 stellar radii][*]
-Z    specify stellar tidal dissipation criterion

As a convenient shorthand, any "dt" interval specified less than zero
is interpreted as a power of 2, i.e. "-d -3" sets dt_log = 0.125.
```

Many of the above options are fairly self-explanatory. For reference (in the choice of default units) ``standard'' N-body units (see Heggie and Mathieu, 1986) are typically defined so that the length and time scales of the system are unity, and the energy is -1/4.

A few specifics:

• The accuracy parameter specified by -a is an overall multiplier to the standard Aarseth time scale (Aarseth 1985), defined in terms of the acceleration and its derivatives at the end of each particle time step. The per-step integration error scales as the fifth power of this parameter. In order to ensure efficient use of the GRAPE hardware, particle time steps are constrained to be powers of 2, allowing many stars to be advanced simultaneously.

• The softening parameter is commonly used to avoid problems associated with the singularity in the point-mass gravitational field. This option is included mainly for historical reasons. It is not used in kira applications.

• The ``exact'' (-E) option means that perturber lists are not used, and that all binaries are resolved into their components.

• The -f (close encounter) and -g (hysteresis) options control the creation and dissolution of ``binary'' pairs in kira. The number specified by -f specifies the threshhold for binary creation; it is scaled internally by 1/N, where N is the total number of stars in the system. (Elementary dynamical considerations indicate that the separation of a 1 kT binary is O(1/N).) The number specified by -g multiplies the binary creation threshhold to determine the maximum separation at which a binary reverts to two separate stars. Internal reorganization of multiple systems is based on simple nearest-neighbor considerations.

• If an external Galactic field is imposed (-Q option), -J specifies the system Jabobi radius, in units of the system virial radius. Tidal parameters may also be passed to kira via the input snapshot using the Starlab story mechanism. Whether or not a tidal field is used, stars may be stripped when they cross the radius defined by -G (again specified in units of the virial radius).

• The perturbation factor specified by -k sets the (dimensionless) threshhold at which a star is considered a perturber of a given binary. Because the perturbative terms are periodic and so do not contribute over a complete binary orbit, the threshhold for fully unperturbed motion is taken substantially higher, usually by a factor of 10 - 100. More technical details on unperturbed binary and multiple motion may be found here.

• The -M, -T, -R, and -q options allow the user to specify explicitly mass-, length-, and time-scales, for use in defining the interface between the dynamics and the tellar evolution portions of the program. To a large extent, these options reproduce functionality already provided by the Starlab tool addstar; they are included in kira to allow the user to override earlier choices, if necessary.

• The -y option specifies the threshhold (in absolute units: dynamical or physical, depending on whether or not stellar evolution is enabled) for diagnostic output on stellar close encounters. The stellar merger criterion is set by -z; the threshhold for tidal dissipation (applied at periastron) is set by -Z (both ``z'' options specify distances in terms of stellar radii).