April 11, 2023

AR-CHAIN

The (more modern) and modified version of Seppo Mikkola’s AR-CHAIN

In particular, this guide is aimed to the modified version of AR-Chain which includes a recipe to merge black holes because of gravitational radiation following Lousto’s fitting equations and the role of an external (non-evolving) stellar potential which mimics the role of a stellar system affecting the system of study.

A manual

Firs of all, there is a detailed publication which explains how to use and run this modified version of AR-Chain, which includes some interesting features which can be switched off, so that we’re back to the default version of the programme (hopefully!).

Guide to the modified version.

The code with an example

This is the full code which corresponds to the pdf guide of the previous section, with an included example.

A step-by-step guide

(1) Choose code units / Conversion to physical units

For example, if we choose

  1. G=1
  2. as unit of mass the MBH,  M = 4.297e6 Msun
  3. as unit of length, D= 10 mpc

In such units,

  1. the speed of light in vacuum is c= 220.52482275418978.
  2. the corresponding unit of time is T = sqrt(D^3 / (G*M)) = 7.19255773861334 years.
  3. the unit of velocity V=D/T= 1359.4499442548772 km/s

Then, to convert from code units to physical units,

  1. Multiply length units by lu = 10 to get lengths in mpc
  2. Multiply time units by tu = 7.19255773861334 to get time years
  3. Multiply mass units by mu = 4.297e6 to get masses in Msun

(2) File PARAMETERS.dat

  1. Number of stellar bodies N, and how many of these are black holes (0 <= Nbh <= N)
  2. Type of units = 0 (our chosen units)
  3. Termination time TMAX over which you want to run the simulation in unit of T
  4. Time-step DeltaT for writing output: you will have an output file with positions and velocities of the bodies every DeltaT.
  5. Standard choice for the cmet parameter is: (1, 1.e-16, 0). See user manual for details.
  6. Speed of light in chosen units: in the previous example we gave is 220.52482275418978.
  7. Typically set the Ixc index to 1 (not very important).
  8. Orbital elements: set it to “yes” if you want the code to calculate the orbital elements for all the bodies, otherwise “no”.
  9. Tolerance. Set it in the numerical integration: Usually set it to 1.e-12
  10. Softening. Set “soft” to 0, unless you want softening of the Newtonian force.
  11. External potential. Unless you want to add an extended, smooth potential, following a Dehnen or Plummer profile, set both of them to “no no”.
  12. Spin values. Set ispinread = 1 and ispinrand = 0: read initial spin for each object from file “spin.dat”. If you set instead ispinread = 0 and ispinrand = 1, the code sets random initial spins for each body.

(3) File INPUT.dat

Each row should contain

  1. mass (in units of M, the mass of the central massive black hole),
  2. initial positions (x0,y0,z0) in units of D and
  3. initial velocities (vx0,vy0,vz0) in units of of each body.

First rows correspond to the black holes, then the other objects.

(4) File SPIN.txt

Each row should contain x, y, z components for the spin of each body (same order as initial conditions, first the black holes and then the stars).

The default version of Seppo Mikkola’s AR-CHAIN

Compilation

gfortran -O3 AR-CHAIN.f -o ARC.x

The executable ARC.x reads the input from an input file, which we call INP.ttt

Run the code

$ ./ARC.x < INP.ttt

Structure of the input file

  1.  a first line containing all the parameters of the run
  2.  a series of N lines containing the particles’ initial conditions in the format: mass, x, y, z, vx, vy, vz
  3.  a final line with as many 0s as the number of parameters in instruction

First line of input

The line must be as follows (this case is tailored to an EMRI-like binary):

0  2   1.d0    800.   0.e-3    0.0   1. 1.e-2 0.e0 456.99  O1  0 2  0.0 0.0
0.0 1.e-13  0  10000 /IWR,N,DELTAT,TMAX ,step,soft,cmet,Clight, out
file,Ixc,Nbh ,spin, EPS, ivelocity, KSMX !

with

IWR          id of the run -- nothing of interest
N            number of particles
DELTAT,TMAX  output frequency and integration time (in scaled units)
step         initial integration step (if set to 0, the code estimate this quantity on its own)
soft         softening parameter (in scaled units). Of course, this can be set to 0
cmet         a three-element vector that sets the integration method. Generally the choice (1 0.01 0) is good for EMRIs
Clight       speed of light in scaled units
outfile      name of the main output file, it contains time, (x, y, z) for all particles in scaled units, thus you have a total of 1 + 3N columns
Ixc          index that tells if exact output times are required,
               0 => no exact, output when t>tnext, (fast, but sometimes output interval too long)
               1 => exact time (not always accurate, but faster than the exact iteration)
               2 => exact iteration (time to tnext) (slower!)
Nbh          number of particles for which PN terms must be considered (these will be considered as the first Nbh lines in the input file)
spin         spin components of the heaviest BH in the sample, it must be placed in the first line (works well just with SMBHs and stellar objects, I'd say)
EPS          tolerance of the run, generally 1.E-12 is a good choice
ivelocity    index telling that there are v-dependent perturbations (=1) or not (=0). Set to 1 to include PN terms
KSMX         tells how many (maximum) integration steps the code is at most allowed to take between outputs

An example

This file contains an example of an EMRI

0  2   1.d1    8000.   0.e-3    0.0   1. 0.e0 0.e0 456.99  O1  0 2  0.0 0.0 0.0 1.e-13  0  10000 /IWR,N,DELTAT,TMAX ,step,soft,cmet,Clight, outfile,Ixc,Nbh ,spin, EPS, ivelocity, KSMX !
1.0 0.0 0.0 0.0 0.0 0.0 0.0
1.E-6 1.0 0.0 0.0 0.0 0.05 0.0
 0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

where we set the

mass unit = 10^6 Msun
length unit = 0.01 pc

This implies a

velocity unit = 656 km/s (thus Clight  = 456.99)
time unit ~ 15 yr

The orbit has an initial semiax of 0.01 pc and eccentricity = 0.997, which are typical of an EMRI.

Additional information

Besides the output file produced by the code, AR-CHAIN prints in the terminal some conservation quantities (not exactly the energy conservation when PN are on, but close to), and it also produces some files where there are the orbital elements of all bodies with respect to the first body in the file. This is particularly well-suited for EMRIs, less for nearly equal-mass bodies.

These files are:

  • axes.dat: time, semimajor axis in scaled units of the orbit between i-th body and the first one
  • eccs.dat : time, eccentricity of the orbit between the i-th body and the first one
  • incs.dat: time, inclination of the orbit between the i-th body and the first one
  • omes.dat: time, argument of periapsis
  • Omes.dat: time, longitude of the ascending node
  • spin.dat: time, spin components of the first body
  • merge.dat: time, indexes of merging bodies

Relativistic terms

The PN terms are implemented in the main code in the routine “Relativistic terms” at line 1773. They are called A2, A2p5, A3, A3p5, B2, B3 etc..