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
- G=1
- as unit of mass the MBH, M = 4.297e6 Msun
- as unit of length, D= 10 mpc
In such units,
- the speed of light in vacuum is c= 220.52482275418978.
- the corresponding unit of time is T = sqrt(D^3 / (G*M)) = 7.19255773861334 years.
- the unit of velocity V=D/T= 1359.4499442548772 km/s
Then, to convert from code units to physical units,
- Multiply length units by lu = 10 to get lengths in mpc
- Multiply time units by tu = 7.19255773861334 to get time years
- Multiply mass units by mu = 4.297e6 to get masses in Msun
(2) File PARAMETERS.dat
- Number of stellar bodies N, and how many of these are black holes (0 <= Nbh <= N)
- Type of units = 0 (our chosen units)
- Termination time TMAX over which you want to run the simulation in unit of T
- Time-step DeltaT for writing output: you will have an output file with positions and velocities of the bodies every DeltaT.
- Standard choice for the cmet parameter is: (1, 1.e-16, 0). See user manual for details.
- Speed of light in chosen units: in the previous example we gave is 220.52482275418978.
- Typically set the Ixc index to 1 (not very important).
- Orbital elements: set it to “yes” if you want the code to calculate the orbital elements for all the bodies, otherwise “no”.
- Tolerance. Set it in the numerical integration: Usually set it to 1.e-12
- Softening. Set “soft” to 0, unless you want softening of the Newtonian force.
- External potential. Unless you want to add an extended, smooth potential, following a Dehnen or Plummer profile, set both of them to “no no”.
- 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
- mass (in units of M, the mass of the central massive black hole),
- initial positions (x0,y0,z0) in units of D and
- 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
- a first line containing all the parameters of the run
- a series of N lines containing the particles’ initial conditions in the format: mass, x, y, z, vx, vy, vz
- 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..