We have simulated the evolution of a galactic nucleus containing triple massive black holes. As we stated in the introduction, currently available numerical simulations suggest that when two galaxies merge, their central black holes will be left as a binary system. Thus, when the remnant of a galaxy merger again merges with a third galaxy, the final merged object is likely to host three massive black holes.
The evolution of a triple black hole system is entirely different from that of a binary system, since it is very unlikely that three black holes form a dynamically stable system. There are basically two possibilities for the final outcome. One is that one or more black holes will be ejected from the galaxy through three-body interaction of black holes. The other possibility is that during the three-body interaction two black holes come close enough that they merge through gravitational wave radiation. Thus, even after black holes have kicked out all nearby stars, they can still evolve in a complex way. On the other hand, a bound pair of two black holes is always dynamically stable. Thus, once a binary has kicked out all nearby stars, the evolution would be completely halted.
In our simulation, the galaxy is modeled with 786,432 equal-mass stars. The black holes are modeled as three point-mass particles with a mass of 1% of the total mass of the system. The relativistic effects were negligible in our simulation.
We performed a simulation for 2.5 dynamical time units, for which the
number of individual steps was 
. The whole
simulation, including on-the-fly analysis of the orbital elements
of black holes and file operations, took 19.45 hours. The total number
of floating point operations is 
, since one particle-particle
interaction amounts to 57 floating point operations. The resulting average computing speed is 992 Gflops.
As a short benchmark run, we tried the same simulation with 1,000,000
stars for 0.5 time units. For this run, the total number of timesteps
is 
, and calculation took 3.43 hours. The average
calculation speed was 1.11 Tflops.