 |
|
N-body
Simulation Testbed
NBodyLab.org
A
laboratory for experimenting with small astrophysical N-body
simulations,
using a desktop GRAPE-6a supercomputer.
|
Some small-N experiments
(Prepared by Vicki Johnson, last modified February 2007)
Contents:
Introduction
Small-N simulations
illustrate complex interactions that must be handled
accurately and efficiently in codes like NBODY4 and NBODY6 for large N. These small-N simulations do not use the GRAPE card.
The simulations for 3-body and 4-body simulations below use the fast and accurate codes triple2 and chain from Sverre Aarseth's
download site,
listed under the "multireg" category. Refer also to the accompanying manual manreg.ps.
Augmenting these codes for iterative experimentation is a straightforward process.
www.sverre.com presents Java animations of
triple and chain and discussion of interesting initial conditions. Chapter 18 of
Aarseth's book
Gravitational
N-body Simulations Tools and Algorithms discusses small-N experiments with his codes and surveys research by others.
The NBodyLab adaptations of triple2 and chain below
support a (greatly simplified) perturbation option for initial conditions, run the standard Fortran distributions of
triple2 and chain on the server side to produce full diagnostic output, and can display animations in 3D.
[Added February 2007] The NBodyLab adaptations of triple2 now includes three body simulations with relativistic effects.
Simulation events
Small-N simulation events have been classified in many ways.
For example, in a study by Giersz and Spurzem 2006,
descriptions of three and four body encounters include the following:
hardening: the remaining two stars are
bound with well defined semi-major axis and eccentricity.
softening: the binary loses binding energy.
three bound: a bound three-body subsystem is formed.
dissolution: the kinetic energy of the incoming star
is large enough to disrupt the binary.
stable: two binaries are well defined and bound after
the four-body encounter.
hierarchical: one star escapes after four-body encounter and a hierarchical triple
remains.
four bound: all four stars form a bound subsystem.
The codes triple2 and chain quantify simulation events such as binary formation.
Perturbations
In the NBodyLab adaptation of triple2 and chain, a very simple perturbation capability is supported:
a single input variable (e.q. x1, vx2, m1, zmp) may be varied linearly, without any rescaling. While this often
produces a wide range of visually interesting outcomes, much more rigor is needed for scientific study. For example,
Heggie 2000 describes a perturbation scheme where
"the coordinates and velocity components were chosen randomly in (0, 1) and then converted to the barycentric frame. By
adding the same suitably chosen angular velocity to all three particles the angular momentum was then set to zero. Finally
velocities, coordinates and masses were scaled to virial equilibrium in standard units".
3-body simulations
Figure 8 periodic orbit
This remarkable periodic orbit was discovered by
Moore in 1993 and
rediscovered by Chenciner and Montgomery in 2001.
Montgomery 2001 states
"numerical experiments done by Heggie suggest that the probability of an eight is somewhere between one per galaxy and one per universe."
Heggie 2000 describes a series of scattering experiments and states
if the coordinates of the initial central
particle are changed by as much as about 0.13 (in standard units),
the 8-like behavior does not persist after about T = 100.
|
|
Run the simulation:
figure 8 orbit
Run the simulation:
figure 8 orbit, perturbing x3
Note binary formations after perturbation > .12
|
|
Earth Jupiter Sun idealized triple system
Pythagorean Problem
These initial conditions for the Pythagorean Problem are given on www.sverre.com
and discussed in Aarseth's book.
The masses are in the ratio 3:4:5, on the apexes of a right triangle.
In 1893 the mathematician
Meissel conjectured these initial conditions would produce a periodic orbit.
Burrau in 1913 and Szebehely and Peters in 1967 investigated this intriguing configuration numerically.
|
|
Run the simulation:
Pythagorean Problem
Run the simulation:
Pythagorean Problem, perturbing m1
Run the simulation:
Pythagorean Problem, perturbing y3
The different pairings are close encounters rather
than binaries, although the two heavies are often bound
|
|
Criss-Cross Orbit
This periodic orbit is described in a
2005 paper by Moore and Nauenberg.
It was first obtained by Henon in 1976, rediscovered by Moore in 1993, and more recently rigorously proved to be a member
of a family of retrograde orbits. The authors state similar orbits exist for different mass ratios - can you find one?
|
|
|
|
4-body simulations
Great Circle Unstable Orbit
Symmetrical 4-body exchange example
3-body simulations with relativistic effects
Triple2 with post-Newtonian terms can be used to simulate pure two-body
motion influenced by general relativity (GR), using a dominant binary and distant
perturber of low mass.
Such runs can be useful for studies of massive black holes (BH) in order
to explore the parameter space where GR effect may be important.
The investigations can also act as guides for more realistic simulations
of large stellar systems. Moreover, the simplicity of the formulation is
beneficial for educational purposes.
Triple2 with post-Newtonian terms can also be used to study the general
three-body problem of massive members (say 3 BHs). In this case only the
most dominant two-body motion will be subject to the GR effect but this
will be OK in practice.
In considering a binary, only the sum of masses is relevant. This also
applies for BH's and GR. Real systems might have one massive BH + a star or
two BH's. Comparable masses were chosen for convenience in the examples below.
The GR effect itself becomes relevant when the orbital velocity reaches
an appreciable fraction of speed of light; say V^2=G*M/a ~1D-04 c^2.
GR effects are controlled by choosing initial c (CLIGHT) instead of
working with physical units for M and a.
GR coalescence is defined as the distance of three Schwarschild radii
(a relativistic concept) which is 6*M/c^2.
A design objective of triple2 is to demonstrate that reasonable
relativistic coalescence times can be obtained without using the full GR effect all the
time. Experiments can be made with different values of IPN and CLIGHT.
For example, IPN=0 specifies beginning the simulation with relativistic radiation
but
other post-Newtonian terms are activated if the time-scale (a/(da/dt)) becomes
sufficiently
short.
For more information, see the section "Post-Newtonian terms" in the
triple2 documentation (manreg.ps),
and
Post-Newtonian N-body simulations (Aarseth 2007).
Circular binary and distant small body perturber
These initial conditions produce a GR (general relativity) inspiral to coalescence. Compare the GR coalescence times,
which indicate computation can be saved by using lower order terms.
|
|
Run the simulation, beginning with only the relativistic radiation:
IPN=0
Run the simulation, beginning with the full relativistic effect:
IPN=3
|
|
Eccentric binary and distant small body perturber
As above, compare the GR coalescence times.
|
|
Run the simulation, beginning with only the relativistic radiation:
IPN=0
Run the simulation, beginning with full relativistic effect:
IPN=3
|
|
For the same semi-major axis the GR coalescence occurs much sooner for large
eccentricity. Starting with the circular binary example,
choosing e=0.60 gives a simple relation: velocities smaller by factor 2
and apocentre distance 1.6 instead of 1 for e=0. Hence velocities 1.25 and
-0.75, and X1=1.6. This makes the total energy same as for circular orbits
as consistency check. Compare the GR coalescence times.
|
|
Run the simulation, beginning with only the relativistic radiation:
IPN=0
Run the simulation, beginning with full relativistic effect:
IPN=3
Run the simulation, IPN=0,
perturbing CLIGHT from 20 to 25
Run the simulation, CLIGHT=25,
perturbing IPN from 0 to 3
|
|
3D Animation
The simulation evolution can be viewed with a 3D Java applet. It may take 30-60 seconds to load on the first use, but will load quickly thereafter.
Try this example.
Start/stop animation | Double click in the plot display area to start/stop
|
Rotate | Left-click and drag
|
Zoom in/out | Left-click while holding the shift key and drag vertically
|
Control animation speed/direction: | Right-click and drag horizontally.
(Mac: hold down Apple key while dragging horizontally)
|
Reset view |
Press Home key to set view to x-y plane
|
Reload | Reload browser window to restart
|
|
Note that time in the animation does not progress linearly, but reflects sampling driven by the integration algorithm.
The so-called regularized time advances more slowly than
real time during the close encounters.
Triple2 (NBodyLab adaptation) input parameters and output
TOL | Tolerance for Bulirsch-Stoer integrator.
|
TCRIT0 | Maximum time (terminates on escape or TCRIT).
|
ECRIT | Energy error criterion (time reversal only).
|
IREV | Time reversal option (IREV = 0: standard).
|
PLOTS | Plot display option. 0 = output only; 1 = add orbit trace plot; 2 = add 3D applet; 3 = orbit trace plots only.
|
AXESLEN | Axes length for plots.
|
unused |
|
unused |
|
ZMP | Mass of Plummer model (zero for isolated system).
|
EPS | Plummer softening (square is saved).
|
CLIGHT | Speed of light in scaled units (3.0D+05/SQRT(GM/R), if GM/R yields km/sec).
No GR effects are included if CLIGHT > 100,000.
GR effects are extremely weak for large values like ~1000, depending on other parameters.
If CLIGHT = 20: strong GR.
Small CLIGHT produces strong relativistic effects.
|
IPN | IPN is used for relativistic effects, and unused for non-GR runs, as indicated by CLIGHT.
IPN options: =1: IPN1; =2: IPN1 + IPN2; =3: IPN1 -> IPN3. (IPN=-1: precession without radiation.)
IPN=3 means full GR is included initially.
IPN=0 can change to 1, 2, 3 if the instantaneous time-scale becomes sufficiently short.
For IPN=0 only the relativistic radiation is included initially; however, the effect
may be extremely small. This continues until the time-scale for the energy
loss ( a/(da/dt) ) becomes sufficiently short (by some reasonable definition),
in which case the first-order precession terms are added. For shorter time-scales
the second-order terms are included (IPN=2) and finally at the shortest times
the third-order terms are also included.
|
M,R,V | Mass, coordinates & velocities for three bodies, each on one row
|
PERTURB PSTART PEND PDELTA | Variable to perturb (e.q. x1, vx1, m1, zmp), and perturbation range and delta.
|
Some of the parameters in
triple2
have been modified or renamed for the NBodyLab adaptation. (Refer to the original triple.f code
listing for the original parameters.)
RGRAV | Gravitational radius (characteristic size=2a for binary)
|
R1 | Distance between the body currently located in position 1
|
R2 | Distance between the body currently located in position 2
|
R | The third 2-body separation
|
ENERGY | Twice the initial total energy
|
RMIN | Minimum two body separation
|
RCOLL | Minimum two-body separation (osculating pericentre)
|
RMAX | MAX(R1,R2,R)
|
NREG | Counter for switching dominant body m3
|
NSTEPS | Number of calls to integrator
|
RIJ | Minimum pairwise separations (osculating pericentre).
|
EBE | Binding energy of binary scaled by total energy
|
RZ | GR coalescence distance (6*M/C^2)
|
V/C | ratio of mean orbital velocity (V^2=M/a) over speed of light
|
TZ | GR time-scale for coalescence. The actual time for a given
two-body orbit is shorter - around a factor of 4 in some cases.
|
Chain/Plummer (NBodyLab adaptation) input parameters and output
N | Particle number.
|
TOL | Tolerance for Bulirsch-Stoer integrator.
|
TCRIT | Maximum time (in units of crossing time).
|
RESC | Escape distance (scaled by gravitational radius).
|
ECRIT | Error criterion (rms for time reversal check).
|
IREV | Time reversal option (IREV = 0: standard).
|
JCOLL | Close encounter (=1) and collision (>1) indicator.
|
unused |
|
PLOTS | Plot display option. 0 = output only; 1 = add orbit trace plot; 2 = add 3D applet; 3 = orbit trace plots only.
|
ZMP | Mass of Plummer model (zero for isolated system).
|
EPS | Plummer softening (square is saved).
|
unused |
|
AXESLEN | Axes length for plots
|
unused |
|
M,R,V | Mass, coordinates & velocities for each of the bodies, one per row
|
SIZE | Individual radii (JCOLL > 0; = 0 without collision).
|
PERTURB PSTART PEND PDELTA | Variable to perturb (e.q. x1, vx1, m1, zmp), and perturbation range and delta.
For example, use x1 0 1.0 .1 for a sequence of runs varying x1 from 0 to 1.0 by .1.
|
Some of the parameters in
chain
have been modified or renamed for the NBodyLab adaptation. (Refer to the original chain.f code
listing for the original parameters.)
TCR | Crossing time.
|
RG | Gravitational radius.
|
MB | Combined mass of two bodies.
|
More information
Aarseth's book, small-N chapters in Hut and Heggie, and on-line N-body
demonstrations by Bob Jenkins.
If you'd like demonstration versions Windows or Mac OS binaries of triple2 or chain, please write to nbodylab@interconnect.com.
To sign up for news, get help or request permission for longer simulations,
write to nbodylab@interconnect.com.
|
