The numerical simulations of massive collisional stellar systems, such as globular clusters (GCs), are very time-consuming. Until now, only a few realistic million-body simulations of GCs with a small fraction of binaries (5%) have been performed by using the NBODY6++GPU code. Such models took half a year computational time on a GPU based super-computer. In this work, we develop a new N-body code, PeTar, by combining the methods of Barnes-Hut tree, Hermite integrator and slow-down algorithmic regularization (SDAR). The code can accurately handle an arbitrary fraction of multiple systems (e.g. binaries, triples) while keeping a high performance by using the hybrid parallelization methods with MPI, OpenMP, SIMD instructions and GPU. A few benchmarks indicate that PeTar and NBODY6++GPU have a very good agreement on the long-term evolution of the global structure, binary orbits and escapers. On a highly configured GPU desktop computer, the performance of a million-body simulation with all stars in binaries by using PeTar is 11 times faster than that of NBODY6++GPU. Moreover, on the Cray XC50 supercomputer, PeTar well scales when number of cores increase. The ten million-body problem, which covers the region of ultra compact dwarfs and nuclearstar clusters, becomes possible to be solved.

An accurate and efficient method dealing with the few-body dynamics is important for simulating collisional N-body systems like star clusters and to follow the formation and evolution of compact binaries. We describe such a method which combines the time-transformed explicit symplectic integrator (Preto & Tremaine 1999; Mikkola & Tanikawa 1999) and the slow-down method (Mikkola & Aarseth 1996). The former conserves the Hamiltonian and the angular momentum for a long-term evolution, while the latter significantly reduces the computational cost for a weakly perturbed binary. In this work, the Hamilton equations of this algorithm are analyzed in detail. We mathematically and numerically show that it can correctly reproduce the secular evolution like the orbit averaged method and also well conserve the angular momentum. For a weakly perturbed binary, the method is possible to provide a few order of magnitude faster performance than the classical algorithm. A publicly available code written in the c++ language, SDAR, is available on GitHub (https://github.com/lwang-astro/SDAR). It can be used either as a stand alone tool or a library to be plugged in other $N$-body codes. The high precision of the floating point to 62 digits is also supported.

$N$-body integrations are used to model a wide range of astrophysical dynamics, but they suffer from errors which make their orbits diverge exponentially in time from the correct orbits. Over long time-scales, their reliability needs to be established. We address this reliability by running a three-body planetary system over about $200$ e-folding times. Using nearby initial conditions, we can construct statistics of the long-term phase-space structure and compare to rough estimates of resonant widths of the system. We compared statistics for a wide range of numerical methods, including a Runge--Kutta method, Wisdom--Holman method, symplectic corrector methods, and a method by Laskar and Robutel. "Improving" an integrator did not increase the phase space accuracy, but simply increasing the number of initial conditions did. In fact, the statistics of a higher order symplectic corrector method were inconsistent with the other methods in one test.

In this paper, we report the implementation and measured performance of our extreme-scale global simulation code on Sunway TaihuLight and two PEZY-SC2 systems: Shoubu System B and Gyoukou. The numerical algorithm is the parallel Barnes-Hut tree algorithm, which has been used in many large-scale astrophysical particle-based simulations. Our implementation is based on our FDPS framework. However, the extremely large numbers of cores of the systems used (10M on TaihuLight and 16M on Gyoukou) and their relatively poor memory and network bandwidth pose new challenges. We describe the new algorithms introduced to achieve high efficiency on machines with low memory bandwidth. The measured performance is 47.9, 10.6 PF, and 1.01PF on TaihuLight, Gyoukou and Shoubu System B (efficiency 40\%, 23.5\% and 35.5\%). The current code is developed for the simulation of planetary rings, but most of the new algorithms are useful for other simulations, and are now available in the FDPS framework.

In this paper, we describe the algorithms we implemented in FDPS to make efficient use of accelerator hardware such as GPGPUs. We have developed FDPS to make it possible for many researchers to develop their own high-performance parallel particle-based simulation programs without spending large amount of time for parallelization and performance tuning. The basic idea of FDPS is to provide a high-performance implementation of parallel algorithms for particle-based simulations in a "generic" form, so that researchers can define their own particle data structure and interparticle interaction functions and supply them to FDPS. FDPS compiled with user-supplied data type and interaction function provides all necessary functions for parallelization, and using those functions researchers can write their programs as though they are writing simple non-parallel program. It has been possible to use accelerators with FDPS, by writing the interaction function that uses the accelerator. However, the efficiency was limited by the latency and bandwidth of communication between the CPU and the accelerator and also by the mismatch between the available degree of parallelism of the interaction function and that of the hardware parallelism. We have modified the interface of user-provided interaction function so that accelerators are more efficiently used. We also implemented new techniques which reduce the amount of work on the side of CPU and amount of communication between CPU and accelerators. We have measured the performance of N-body simulations on a systems with NVIDIA Volta GPGPU using FDPS and the achieved performance is around 27 \% of the theoretical peak limit. We have constructed a detailed performance model, and found that the current implementation can achieve good performance on systems with much smaller memory and communication bandwidth.

The new technique of measuring frequency by optical lattice clocks now approaches to the relative precision of $(\Delta f/f)=O(10^{-18})$. We propose to place such precise clocks in space and to use Doppler tracking method for detecting low-frequency gravitational wave below 1 Hz. Our idea is to locate three spacecrafts at one A.U. distance (say at L1, L4 & L5 of the Sun-Earth orbit), and apply the Doppler tracking method by communicating "the time" each other. Applying the current available technologies, we obtain the sensitivity for gravitational wave with three or four-order improvement ($h_{\rm n}\sim 10^{-17}$ or $10^{-18}$ level in $10^{-5}$Hz -- $1$ Hz) than that of Cassini spacecraft in 2001. This sensitivity enables us to observe black-hole mergers of their mass greater than $10^5 M_\odot$ in the cosmological scale. Based on the hierarchical growth model of black-holes in galaxies, we estimate the event rate of detection will be 20-50 a year. We nickname "INO" (Interplanetary Network of Optical Lattice Clocks) for this system, named after Tadataka Ino (1745--1818), a Japanese astronomer, cartographer, and geodesist.

In most of mesh-free methods, the calculation of interactions between sample points or particles is the most time consuming. When we use mesh-free methods with high spatial orders, the order of the time integration should also be high. If we use usual Runge-Kutta schemes, we need to perform the interaction calculation multiple times per one time step. One way to reduce the number of interaction calculations is to use Hermite schemes, which use the time derivatives of the right hand side of differential equations, since Hermite schemes require smaller number of interaction calculations than RK schemes do to achieve the same order. In this paper, we construct a Hermite scheme for a mesh-free method with high spatial orders. We performed several numerical tests with fourth-order Hermite schemes and Runge-Kutta schemes. We found that, for both of Hermite and Runge-Kutta schemes, the overall error is determined by the error of spatial derivatives, for timesteps smaller than the stability limit. The calculation cost at the timestep size of the stability limit is smaller for Hermite schemes. Therefore, we conclude that Hermite schemes are more efficient than Runge-Kutta schemes and thus useful for high-order mesh-free methods for Lagrangian Hydrodynamics.

Numerical simulations based on particle methods have been widely used in various fields including astrophysics. To date, simulation softwares have been developed by individual researchers or research groups in each field, with a huge amount of time and effort, even though numerical algorithms used are very similar. To improve the situation, we have developed a framework, called FDPS, which enables researchers to easily develop massively parallel particle simulation codes for arbitrary particle methods. Until version 3.0, FDPS have provided API only for C++ programing language. This limitation comes from the fact that FDPS is developed using the template feature in C++, which is essential to support arbitrary data types of particle. However, there are many researchers who use Fortran to develop their codes. Thus, the previous versions of FDPS require such people to invest much time to learn C++. This is inefficient. To cope with this problem, we newly developed a Fortran interface layer in FDPS, which provides API for Fortran. In order to support arbitrary data types of particle in Fortran, we design the Fortran interface layer as follows. Based on a given derived data type in Fortran representing particle, a Python script provided by us automatically generates a library that manipulates the C++ core part of FDPS. This library is seen as a Fortran module providing API of FDPS from the Fortran side and uses C programs internally to interoperate Fortran with C++. In this way, we have overcome several technical issues when emulating `template' in Fortran. By using the Fortran interface, users can develop all parts of their codes in Fortran. We show that the overhead of the Fortran interface part is sufficiently small and a code written in Fortran shows a performance practically identical to the one written in C++.

Advanced LIGO achieved the first detection of the gravitational wave, which was from a merging binary black hole (BBH). In the near future, more merger events will be observed, and the mass distribution of them will become available. The mass distribution of merger events reflects the evolutionary path of BBHs: dynamical formation in dense star clusters or common envelope evolution from primordial binaries. In this paper, we estimate the detection rate of merging BBHs which dynamically formed in dense star clusters by combining the results of $N$-body simulations, modeling of globular clusters, and cosmic star-cluster formation history. We estimate that the merger rate density in the local universe within the redshift of 0.1 is 13--57 Gpc$^{-3}$ yr$^{-1}$. We find that the detection rate is 0.23--4.6 per year for the current sensitivity limit and that it would increase to 5.1--99 per year for the designed sensitivity which will be achieved in 2019. The distribution of merger rate density in the local universe as a function of redshifted chirp mass has a peak close to the low-mass end.The chirp mass function of the detected mergers, on the other hand, has a peak at the high-mass end, but is almost flat. This difference is simply because the detection range is larger for more massive BBHs.

We investigated whether outward Planetesimal Driven Migration (PDM) takes place or not in simulations when the self gravity of planetesimals is included. We performed $N$-body simulations of planetesimal disks with large width (0.7 - 4AU) which ranges over the ice line. The simulations consisted of two stages. The first stage simulations were carried out to see the runaway growth phase using the planetesimals of initially the same mass. The runaway growth took place both at the inner edge of the disk and at the region just outside the ice line. This result was utilized for the initial setup of the second stage simulations in which the runaway bodies just outside the ice line were replaced by the protoplanets with about the isolation mass. In the second stage simulations, the outward migration of the protoplanet was followed by the stopping of the migration due to the increase of the random velocity of the planetesimals. Due to this increase of random velocities, one of the PDM criteria derived in Minton and Levison (2014) was broken. In the current simulations, the effect of the gas disk is not considered. It is likely that the gas disk plays an important role in planetesimal driven migration, and we plan to study its effect in future papers.

In fluid dynamical simulations in astrophysics, large deformations are common and surface tracking is sometimes necessary. Smoothed Particle Hydrodynamics (SPH) method has been used in many of such simulations. Recently, however, it has been shown that SPH cannot handle contact discontinuities or free surfaces accurately. There are several reasons for this problem. The first one is that SPH requires that the density is continuous and differentiable. The second one is that SPH does not have the consistency, and thus the accuracy is zeroth order in space. In addition, we cannot express accurate boundary conditions with SPH. In this paper, we propose a novel, high-order scheme for particle-based hydrodynamics of compress- ible fluid. Our method is based on kernel-weighted high-order fitting polynomial for intensive variables. With this approach, we can construct a scheme which solves all of the three prob- lems described above. For shock capturing, we use a tensor form of von-Neumann-Richtmyer artificial viscosity. We have applied our method to many test problems and obtained excel- lent result. Our method is not conservative, since particles do not have mass or energy, but only their densities. However, because of the Lagrangian nature of our scheme, the violation of the conservation laws turned out to be small. We name this method Consistent Particle Hydrodynamics in Strong Form (CPHSF).

The giant impact (GI) is one of the most important hypotheses both in planetary science and geoscience, since it is related to the origin of the Moon and also the initial condition of the Earth. A number of numerical simulations have been done using the smoothed particle hydrodynamics (SPH) method. However, GI hypothesis is currently in a crisis. The "canonical" GI scenario failed to explain the identical isotope ratio between the Earth and the Moon. On the other hand, little has been known about the reliability of the result of GI simulations. In this paper, we discuss the effect of the resolution on the results of the GI simulations by varying the number of particles from $3 \times10^3$ to $10^8$. We found that the results does not converge, but shows oscillatory behaviour. We discuss the origin of this oscillatory behaviour.

We present the basic idea, implementation, measured performance and performance model of FDPS (Framework for developing particle simulators). FDPS is an application-development framework which helps the researchers to develop particle-based simulation programs for large-scale distributed-memory parallel supercomputers. A particle-based simulation program for distributed-memory parallel computers needs to perform domain decomposition, redistribution of particles, and gathering of particle information for interaction calculation. Also, even if distributed-memory parallel computers are not used, in order to reduce the amount of computation, algorithms such as Barnes-Hut tree method should be used for long-range interactions. For short-range interactions, some methods to limit the calculation to neighbor particles are necessary. FDPS provides all of these necessary functions for efficient parallel execution of particle-based simulations as "templates", which are independent of the actual data structure of particles and the functional form of the interaction. By using FDPS, researchers can write their programs with the amount of work necessary to write a simple, sequential and unoptimized program of O(N^2) calculation cost, and yet the program, once compiled with FDPS, will run efficiently on large-scale parallel supercomputers. A simple gravitational N-body program can be written in around 120 lines. We report the actual performance of these programs and the performance model. The weak scaling performance is very good, and almost linear speedup was obtained for up to the full system of K computer. The minimum calculation time per timestep is in the range of 30 ms (N=10^7) to 300 ms (N=10^9). These are currently limited by the time for the calculation of the domain decomposition and communication necessary for the interaction calculation. We discuss how we can overcome these bottlenecks.

The Santa Barbara cluster comparison project (Frenk et al. Frenk+1999) revealed that there is a systematic difference between entropy profiles of clusters of galaxies obtained by Eulerian mesh and Lagrangian smoothed particle hydrodynamics (SPH) codes: Mesh codes gave a core with a constant entropy whereas SPH codes did not. One possible reason for this difference is that mesh codes are not Galilean invariant. Another possible reason is the problem of the SPH method, which might give too much "protection" to cold clumps because of the unphysical surface tension induced at contact discontinuities. In this paper, we apply the density independent formulation of SPH (DISPH), which can handle contact discontinuities accurately, to simulations of a cluster of galaxies, and compare the results with those with the standard SPH. We obtained the entropy core when we adopt DISPH. The size of the core is, however, significantly smaller than those obtained with mesh simulations, and is comparable to those obtained with quasi-Lagrangian schemes such as "moving mesh" and "mesh free" schemes. We conclude that both the standard SPH without artificial conductivity and Eulerian mesh codes have serious problems even such an idealized simulation, while DISPH, SPH with artificial conductivity, and quasi-Lagrangian schemes have sufficient capability to deal with it.

Hydrodynamical simulations of rotating disk play important roles in the field of astrophysical and planetary science. Smoothed Particle Hydrodynamics (SPH) has been widely used for such simulations. It, however, has been known that with SPH, a cold and thin Kepler disk breaks up due to the unwanted angular momentum transfer. Two possible reasons have been suggested for this breaking up of the disk; the artificial viscosity (AV) and the numerical error in the evaluation of pressure gradient in SPH. Which one is dominant has been still unclear. In this paper, we investigate the reason for this rapid breaking up of the disk. We implemented most of popular formulations of AV and switches and measured the angular momentum transfer due to both AV and the error of SPH estimate of pressure gradient. We found that the angular momentum transfer due to AV at the inner edge triggers the breaking up of the disk. We also found that the classical von-Neumann-Richtmyer-Landshoff type AV with a high order estimate for $\nabla \cdot \vec{v}$ can maintain the disk for $\sim 100$ orbits even when used with the standard formulation of SPH.

At present, the giant impact (GI) is the most widely accepted model for the origin of the Moon. Most of the numerical simulations of GI have been carried out with the smoothed particle hydrodynamics (SPH) method. Recently, however, it has been pointed out that standard formulation of SPH (SSPH) has difficulties in the treatment of a contact discontinuity such as a core-mantle boundary and a free surface such as a planetary surface. This difficulty comes from the assumption of differentiability of density in SSPH. We have developed an alternative formulation of SPH, density independent SPH (DISPH), which is based on differentiability of pressure instead of density to solve the problem of a contact discontinuity. In this paper, we report the results of the GI simulations with DISPH and compare them with those obtained with SSPH. We found that the disk properties, such as mass and angular momentum produced by DISPH is different from that of SSPH. In general, the disks formed by DISPH are more compact: while formation of a smaller mass moon for low-oblique impacts is expected with DISPH, inhibition of ejection would promote formation of a larger mass moon for high-oblique impacts. Since only the improvement of core-mantle boundary significantly affects the properties of circumplanetary disks generated by GI and DISPH has not been significantly improved from SSPH for a free surface, we should be very careful when some conclusions are drawn from the numerical simulations for GI. And it is necessary to develop the numerical hydrodynamical scheme for GI that can properly treat the free surface as well as the contact discontinuity.

We describe the implementation and performance of the ${\rm P^3T}$ (Particle-Particle Particle-Tree) scheme for simulating dense stellar systems. In ${\rm P^3T}$, the force experienced by a particle is split into short-range and long-range contributions. Short-range forces are evaluated by direct summation and integrated with the fourth order Hermite predictor-corrector method with the block timesteps. For long-range forces, we use a combination of the Barnes-Hut tree code and the leapfrog integrator. The tree part of our simulation environment is accelerated using graphical processing units (GPU), whereas the direct summation is carried out on the host CPU. Our code gives excellent performance and accuracy for star cluster simulations with a large number of particles even when the core size of the star cluster is small.

As an entry for the 2012 Gordon-Bell performance prize, we report performance results of astrophysical N-body simulations of one trillion particles performed on the full system of K computer. This is the first gravitational trillion-body simulation in the world. We describe the scientific motivation, the numerical algorithm, the parallelization strategy, and the performance analysis. Unlike many previous Gordon-Bell prize winners that used the tree algorithm for astrophysical N-body simulations, we used the hybrid TreePM method, for similar level of accuracy in which the short-range force is calculated by the tree algorithm, and the long-range force is solved by the particle-mesh algorithm. We developed a highly-tuned gravity kernel for short-range forces, and a novel communication algorithm for long-range forces. The average performance on 24576 and 82944 nodes of K computer are 1.53 and 4.45 Pflops, which correspond to 49% and 42% of the peak speed.

In this paper, we present a new formulation of smoothed particle hydrodynamics (SPH), which, unlike the standard SPH (SSPH), is well-behaved at the contact discontinuity. The SSPH scheme cannot handle discontinuities in density (e.g. the contact discontinuity and the free surface), because it requires that the density of fluid is positive and continuous everywhere. Thus there is inconsistency in the formulation of the SSPH scheme at discontinuities of the fluid density. To solve this problem, we introduce a new quantity associated with particles and "density" of that quantity. This "density" evolves through the usual continuity equation with an additional artificial diffusion term, in order to guarantee the continuity of "density". We use this "density" or pseudo density, instead of the mass density, to formulate our SPH scheme. We call our new method as SPH with smoothed pseudo-density (SPSPH). We show that our new scheme is physically consistent and can handle discontinuities quite well.

Many massive objects have been found in the outer region of the Solar system. How they were formed and evolved has not been well understood, although there have been intensive studies on accretion process of terrestrial planets. One of the mysteries is the existence of binary planetesimals with near-equal mass components and highly eccentric orbits. These binary planetesimals are quite different from the satellites observed in the asteroid belt region. The ratio of the Hill radius to the physical radius of the planetesimals is much larger for the outer region of the disk, compared to the inner region of the disk. The Hill radius increases with the semi major axis. Therefore, planetesimals in the outer region can form close and eccentric binaries, while those in the inner region would simply collide. In this paper, we carried out $N$-body simulations in different regions of the disk and studied if binaries form in the outer region of the disk. We found that large planetesimals tend to form binaries. A significant fraction of large planetesimals are components of the binaries. Planetesimals that become the components of binaries eventually collide with a third body, through three-body encounters. Thus, the existence of binaries can enhance the growth rate of planetesimals in the Trans-Neptunian Object (TNO) region.

At the moment of deepest core collapse, a star cluster core contains less than ten stars. This small number makes the traditional treatment of hard binary formation, assuming a homogeneous background density, suspect. In a previous paper, we have found that indeed the conventional wisdom of binary formation, based on three-body encounters, is incorrect. Here we refine that insight, by further dissecting the subsequent steps leading to hard binary formation. For this purpose, we add some analysis tools in order to make the study less subjective. We find that the conventional treatment does remain valid for direct three-body scattering, but fails for resonant three-body scattering. Especially democratic resonance scattering, which forms an important part of the analytical theory of three-body binary formation, takes too much space and time to be approximated as being isolated, in the context of a cluster core around core collapse. We conclude that, while three-body encounters can be analytically approximated as isolated, subsequent strong perturbations typically occur whenever those encounters give rise to democratic resonances. We present analytical estimates postdicting our numerical results. If we only had been a bit more clever, we could have predicted this qualitative behaviour.

We present a method to couple N-body star cluster simulations to a cosmological tidal field, using the Astrophysical Multipurpose Software Environment. We apply this method to star clusters embedded in the CosmoGrid dark matter-only LambdaCDM simulation. Our star clusters are born at z = 10 (corresponding to an age of the Universe of about 500 Myr) by selecting a dark matter particle and initializing a star cluster with 32,000 stars on its location. We then follow the dynamical evolution of the star cluster within the cosmological environment. We compare the evolution of star clusters in two Milky-Way size haloes with a different accretion history. The mass loss of the star clusters is continuous irrespective of the tidal history of the host halo, but major merger events tend to increase the rate of mass loss. From the selected two dark matter haloes, the halo that experienced the larger number of mergers tends to drive a smaller mass loss rate from the embedded star clusters, even though the final masses of both haloes are similar. We identify two families of star clusters: native clusters, which become part of the main halo before its final major merger event, and the immigrant clusters, which are accreted upon or after this event; native clusters tend to evaporate more quickly than immigrant clusters. Accounting for the evolution of the dark matter halo causes immigrant star clusters to retain more mass than when the z=0 tidal field is taken as a static potential. The reason for this is the weaker tidal field experienced by immigrant star clusters before merging with the larger dark matter halo.

A compact gas cloud G2 is predicted to reach the pericenter of its orbit around the super massive black hole (SMBH) of our galaxy, Sagittarius A* (Sgr A*). This event will give us a rare opportunity to observe the interaction between SMBH and gas around it. We report the result of the fully three-dimensional simulation of the evolution of G2 during the first pericenter passage. The strong tidal force by the SMBH stretches the cloud along its orbit, and compresses it strongly in the vertical direction, resulting in the heating up and flaring up of the cloud. The bolometric luminosity will reach the maximum of $\sim100 L_{\odot}$. This flare should be easily observed in the near infrared.

The smoothed particle hydrodynamics (SPH) method is a useful numerical tool for the study of a variety of astrophysical and planetlogical problems. However, it turned out that the standard SPH algorithm has problems in dealing with hydrodynamical instabilities. This problem is due to the assumption that the local density distribution is differentiable. In order to solve this problem, a new SPH formulation, which does not require the differentiability of the density, have been proposed. This new SPH method improved the treatment of hydrodynamical instabilities. This method, however, is applicable only to the equation of state (EOS) of the ideal gas. In this paper, we describe how to extend the new SPH method to non-ideal EOS. We present the results of various standard numerical tests for non-ideal EOS. Our new method works well for non-ideal EOS. We conclude that our new SPH can handle hydrodynamical instabilities for an arbitrary EOS and that it is an attractive alternative to the standard SPH.

Earth-mass dark matter microhalos with size of ~100 AUs are the first structures formed in the universe, if the dark matter of the Universe are made of neutralino. Here, we report the results of ultra-high-resolution simulations of the formation and evolution of these microhalos. We found that microhalos have the central density cusps of the form $\rho \propto r^{-1.5}$, much steeper than the cusps of larger dark halos. The central regions of these microhalos survive the encounters with stars except in very inner region of the galaxy down to the radius of a few hundreds pcs from the galactic center. The annihilation signals from nearest microhalos are observed as gamma-ray point-sources (radius less than 1'), with unusually large proper motions of ~0.2 degree per year. Their surface brightnesses are ~10% of that of the galactic center. Their S/N ratios might be better if they are far from the galactic plane. Luminosities of subhalos are determined only by their mass, and they are more than one order of magnitude luminous than the estimation by Springel et al. (2008): A boost factor can be larger than 1000. Perturbations to the millisecond pulsars by gravitational attractions of nearby earth-mass microhalos can be detected by the observations of Parkes Pulsar Timing Array (PPTA).

We present the results of the "Cosmogrid" cosmological N-body simulation suites based on the concordance LCDM model. The Cosmogrid simulation was performed in a 30Mpc box with 2048^3 particles. The mass of each particle is 1.28x10^5 Msun, which is sufficient to resolve ultra-faint dwarfs. We found that the halo mass function shows good agreement with the Sheth & Tormen fitting function down to ~10^7 Msun. We have analyzed the spherically averaged density profiles of the three most massive halos which are of galaxy group size and contain at least 170 million particles. The slopes of these density profiles become shallower than -1 at the inner most radius. We also find a clear correlation of halo concentration with mass. The mass dependence of the concentration parameter cannot be expressed by a single power law, however a simple model based on the Press-Schechter theory proposed by Navarro et al. gives reasonable agreement with this dependence. The spin parameter does not show a correlation with the halo mass. The probability distribution functions for both concentration and spin are well fitted by the log-normal distribution for halos with the masses larger than ~10^8 Msun. The subhalo abundance depends on the halo mass. Galaxy-sized halos have 50% more subhalos than ~10^{11} Msun halos have.

The standard formulation of the smoothed particle hydrodynamics (SPH) assumes that the local density distribution is differentiable. This assumption is used to derive the spatial derivatives of other quantities. However, this assumption breaks down at the contact discontinuity. At the contact discontinuity, the density of the low-density side is overestimated while that of the high-density side is underestimated. As a result, the pressure of the low (high) density side is over (under) estimated. Thus, unphysical repulsive force appears at the contact discontinuity, resulting in the effective surface tension. This tension suppresses fluid instabilities. In this paper, we present a new formulation of SPH, which does not require the differentiability of density. Instead of the mass density, we adopt the internal energy density (pressure), and its arbitrary function, which are smoothed quantities at the contact discontinuity, as the volume element used for the kernel integration. We call this new formulation density independent SPH (DISPH). It handles the contact discontinuity without numerical problems. The results of standard tests such as the shock tube, Kelvin-Helmholtz and Rayleigh-Taylor instabilities, point like explosion, and blob tests are all very favorable to DISPH. We conclude that DISPH solved most of known difficulties of the standard SPH, without introducing additional numerical diffusion or breaking the exact force symmetry or energy conservation. Our new SPH includes the formulation proposed by Ritchie & Thomas (2001) as a special case. Our formulation can be extended to handle a non-ideal gas easily.

We present a data format for the output of general N-body simulations, allowing the presence of individual time steps. By specifying a standard, different N-body integrators and different visualization and analysis programs can all share the simulation data, independent of the type of programs used to produce the data. Our Particle Stream Data Format, PSDF, is specified in YAML, based on the same approach as XML but with a simpler syntax. Together with a specification of PSDF, we provide background and motivation, as well as specific examples in a variety of computer languages. We also offer a web site from which these examples can be retrieved, in order to make it easy to augment existing codes in order to give them the option to produce PSDF output.

Ultraluminous infrared galaxies (ULIRGs) with multiple ($\ge 3$) nuclei are frequently observed. It has been suggested that these nuclei are produced by multiple major mergers of galaxies. The expected rate of such mergers is, however, too low to reproduce the observed number of ULIRGs with multiple nuclei. We have performed high-resolution simulations of the merging of two gas-rich disk galaxies. We found that extremely massive and compact star clusters form from the strongly disturbed gas disks after the first or second encounter between the galaxies. The mass of such clusters reaches $\sim 10^8 M_{\odot}$, and their half-mass radii are $20-30 \rm{pc}$. Since these clusters consist of young stars, they appear to be several bright cores in the galactic central region ($\sim \rm{kpc}$). The peak luminosity of these clusters reaches $\sim 10%$ of the total luminosity of the merging galaxy. These massive and compact clusters are consistent with the characteristics of the observed multiple nuclei in ULIRGs. Multiple mergers are not necessary to explain multiple nuclei in ULIRGs.

The conventional wisdom for the formation of the first hard binary in core collapse is that three-body interactions of single stars form many soft binaries, most of which are quickly destroyed, but eventually one of them survives. We report on direct N-body simulations to test these ideas, for the first time. We find that both assumptions are often incorrect: 1) quite a few three-body interactions produce a hard binary from scratch; 2) and in many cases there are more than three bodies directly and simultaneously involved in the production of the first binary. The main reason for the discrepancies is that the core of a star cluster, at the first deep collapse, contains typically only five or so stars. Therefore, the homogeneous background assumption, which still would be reasonable for, say, 25 stars, utterly breaks down. There have been some speculations in this direction, but we demonstrate this result here explicitly, for the first time.

It has been believed that spirals in pure stellar disks, especially the ones spontaneously formed, decay in several galactic rotations due to the increase of stellar velocity dispersions. Therefore, some cooling mechanism, for example dissipational effects of the interstellar medium, was assumed to be necessary to keep the spiral arms. Here we show that stellar disks can maintain spiral features for several tens of rotations without the help of cooling, using a series of high-resolution three-dimensional $N$-body simulations of pure stellar disks. We found that if the number of particles is sufficiently large, e.g., $3\times 10^6$, multi-arm spirals developed in an isolated disk can survive for more than 10 Gyrs. We confirmed that there is a self-regulating mechanism that maintains the amplitude of the spiral arms. Spiral arms increase Toomre's $Q$ of the disk, and the heating rate correlates with the squared amplitude of the spirals. Since the amplitude itself is limited by the value of $Q$, this makes the dynamical heating less effective in the later phase of evolution. A simple analytical argument suggests that the heating is caused by gravitational scattering of stars by spiral arms, and that the self-regulating mechanism in pure-stellar disks can effectively maintain spiral arms on a cosmological timescale. In the case of a smaller number of particles, e.g., $3\times 10^5$, spiral arms grow faster in the beginning of the simulation (while $Q$ is small) and they cause a rapid increase of $Q$. As a result, the spiral arms become faint in several Gyrs.

As to April 2010, 48 TNO (trans-Neptunian Object) binaries have been found. This is about 6% of known TNOs. However, in previous theoretical studies of planetary formation in the TNO region, the effect of binary formation has been neglected. TNO binaries can be formed through a variety of mechanisms, such as three-body process, dynamical friction on two massive bodies, inelastic collisions between two bodies etc. Most of these mechanisms become more effective as the distance from the Sun increases. In this paper, we studied three-body process using direct N-body simulations. We systematically changed the distance from the Sun, the number density of planetesimals, and the radius of the planetesimals and studied the effect of the binaries on the collision rate of planetesimals. In the TNO region, binaries are involved in 1/3 - 1/2 of collisions, and the collision rate is increased by about a factor of a few compared to the theoretical estimate for the direct two-body collisions. Thus, it is possible that the binaries formed through three-body process significantly enhance the collision rate and reduce the growth time scale. In the terrestrial planet region, binaries are less important, because the ratio between the Hill radius and physical size of the planetesimals is relatively small. Although the time scale of our simulations is short, they clearly demonstrated that the accretion process in the TNO region is quite different from that in the terrestrial planet region. Simulations which cover longer time scale are required to obtain more accurate estimate for the accretion enhancement.

In this paper, we present a new hybrid algorithm for the time integration of collisional N-body systems. In this algorithm, gravitational force between two particles is divided into short-range and long-range terms, using a distance-dependent cutoff function. The long-range interaction is calculated using the tree algorithm and integrated with the constant-timestep leapfrog integrator. The short-range term is calculated directly and integrated with the high-order Hermite scheme. We can reduce the calculation cost per orbital period from O(N^2) to O(N log N), without significantly increasing the long-term integration error. The results of our test simulations show that close encounters are integrated accurately. Long-term errors of the total energy shows random-walk behaviour, because it is dominated by the error caused by tree approximation.

We studied the formation process of star clusters using high-resolution N-body/smoothed particle hydrodynamcs simulations of colliding galaxies. The total number of particles is 1.2x10^8 for our high resolution run. The gravitational softening is 5 pc and we allow gas to cool down to \sim 10 K. During the first encounter of the collision, a giant filament consists of cold and dense gas found between the progenitors by shock compression. A vigorous starburst took place in the filament, resulting in the formation of star clusters. The mass of these star clusters ranges from 10^{5-8} Msun. These star clusters formed hierarchically: at first small star clusters formed, and then they merged via gravity, resulting in larger star clusters.

We propose a symmetrized form of the softened gravitational potential which is a natural extension of the Plummer potential. The gravitational potential at the position of particle i (x_i,y_i,z_i), induced by particle j at (x_j,y_j,z_j), is given by: phi_ij = -G m_j/|r_ij^2+e_i^2+e_j^2|^1/2, where G is the gravitational constant, m_j is the mass of particle j, r_ij = |(x_i-x_j)^2+(y_i-y_j)^2+(z_i-z_j)^2|^1/2 and e_i and e_j are the gravitational softening lengths of particles i and j, respectively. This form is formally an extension of the Newtonian potential to five dimensions. The derivative of this equation in the x,y, and z directions correspond to the gravitational accelerations in these directions and these are always symmetric between two particles. When one applies this potential to a group of particles with different softening lengths, as is the case with a tree code, an averaged gravitational softening length for the group can be used. We find that the most suitable averaged softening length for a group of particles is <e_j^2> = sum_j^N m_j e_j^2 / M, where M = sum_j^N m_j and N are the mass and number of all particles in the group, respectively. The leading error related to the softening length is O(sum_j r_j d(e_j^2)/r_ij^3), where r_j is the distance between particle j and the center of mass of the group and d(e_j^2) = e_j^2 - <e_j^2>. Using this averaged gravitational softening length with the tree method, one can use a single tree to evaluate the gravitational forces for a system of particles with a wide variety of gravitational softening lengths. Consequently, this will reduce the calculation cost of the gravitational force for such a system with different softenings without the need for complicated forms of softening. We present the result of simple numerical tests. We found that our modification of the Plummer potential works well.

We report on the performance of our cold-dark matter cosmological N-body simulation which was carried out concurrently using supercomputers across the globe. We ran simulations on 60 to 750 cores distributed over a variety of supercomputers in Amsterdam (the Netherlands, Europe), in Tokyo (Japan, Asia), Edinburgh (UK, Europe) and Espoo (Finland, Europe). Regardless the network latency of 0.32 seconds and the communication over 30.000 km of optical network cable we are able to achieve about 87% of the performance compared to an equal number of cores on a single supercomputer. We argue that using widely distributed supercomputers in order to acquire more compute power is technically feasible, and that the largest obstacle is introduced by local scheduling and reservation policies.

In recent numerical simulations \citep{matsubayashi07,lockmann08}, it has been found that the eccentricity of supermassive black hole(SMBH) - intermediate black hole(IMBH) binaries grows toward unity through interactions with stellar background. This increase of eccentricity reduces the merging timescale of the binary through the gravitational radiation to the value well below the Hubble Time. It also gives the theoretical explanation of the existence of eccentric binary such as that in OJ287 \citep{lehto96, valtonen08}. In self-consistent N-body simulations, this increase of eccentricity is always observed. On the other hand, the result of scattering experiment between SMBH binaries and field stars \citep{quinlan96} indicated no increase of eccentricity. This discrepancy leaves the high eccentricity of the SMBH binaries in $N$-body simulations unexplained. Here we present a stellar-dynamical mechanism that drives the increase of the eccentricity of an SMBH binary with large mass ratio. There are two key processes involved. The first one is the Kozai mechanism under non-axisymmetric potential, which effectively randomizes the angular momenta of surrounding stars. The other is the selective ejection of stars with prograde orbits. Through these two mechanisms, field stars extract the orbital angular momentum of the SMBH binary. Our proposed mechanism causes the increase in the eccentricity of most of SMBH binaries, resulting in the rapid merger through gravitational wave radiation. Our result has given a definite solution to the "last-parsec problem".

Within the distance of 1 pc from the Galactic center (GC), more than 100 young massive stars have been found. The massive stars at 0.1-1 pc from the GC are located in one or two disks, while those within 0.1 pc from the GC, S-stars, have an isotropic distribution. How these stars are formed is not well understood, especially for S-stars. Here we propose that a young star cluster with an intermediate-mass black hole (IMBH) can form both the disks and S-stars. We performed a fully self-consistent $N$-body simulation of a star cluster near the GC. Stars escaped from the tidally disrupted star cluster were carried to the GC due to an 1:1 mean motion resonance with the IMBH formed in the cluster. In the final phase of the evolution, the eccentricity of the IMBH becomes very high. In this phase, stars carried by the 1:1 resonance with the IMBH were dropped from the resonance and their orbits are randomized by a chaotic Kozai mechanism. The mass function of these carried stars is extremely top-heavy within 10". The surface density distributions of young massive stars has a slope of -1.5 within 10" from the GC. The distribution of stars in the most central region is isotropic. These characteristics agree well with those of stars observed within 10" from the GC.

We describe a new algorithm for the integration of self-gravitating fluid systems using SPH method. We split the Hamiltonian of a self-gravitating fluid system to the gravitational potential and others (kinetic and internal energies) and use different time-steps for their integrations. The time integration is done in the way similar to that used in the mixed variable or multiple stepsize symplectic schemes. We performed three test calculations. One was the spherical collapse and the other was an explosion. We also performed a realistic test, in which the initial model was taken from a simulation of merging galaxies. In all test calculations, we found that the number of time-steps for gravitational interaction were reduced by nearly an order of magnitude when we adopted our integration method. In the case of the realistic test, in which the dark matter potential dominates the total system, the total calculation time was significantly reduced. Simulation results were almost the same with those of simulations with the ordinary individual time-step method. Our new method achieves good performance without sacrificing the accuracy of the time integration.

Understanding the universe is hampered by the elusiveness of its most common constituent, cold dark matter. Almost impossible to observe, dark matter can be studied effectively by means of simulation and there is probably no other research field where simulation has led to so much progress in the last decade. Cosmological N-body simulations are an essential tool for evolving density perturbations in the nonlinear regime. Simulating the formation of large-scale structures in the universe, however, is still a challenge due to the enormous dynamic range in spatial and temporal coordinates, and due to the enormous computer resources required. The dynamic range is generally dealt with by the hybridization of numerical techniques. We deal with the computational requirements by connecting two supercomputers via an optical network and make them operate as a single machine. This is challenging, if only for the fact that the supercomputers of our choice are separated by half the planet, as one is located in Amsterdam and the other is in Tokyo. The co-scheduling of the two computers and the 'gridification' of the code enables us to achieve a 90% efficiency for this distributed intercontinental supercomputer.

In this paper, we describe the implementation and performance of GreeM, a massively parallel TreePM code for large-scale cosmological N-body simulations. GreeM uses a recursive multi-section algorithm for domain decomposition. The size of the domains are adjusted so that the total calculation time of the force becomes the same for all processes. The loss of performance due to non-optimal load balancing is around 4%, even for more than 10^3 CPU cores. GreeM runs efficiently on PC clusters and massively-parallel computers such as a Cray XT4. The measured calculation speed on Cray XT4 is 5 \times 10^4 particles per second per CPU core, for the case of an opening angle of \theta=0.5, if the number of particles per CPU core is larger than 10^6.

We introduce AstroSim, a Second Life based prototype application for synchronous collaborative visualization targeted at astronomers.

Supermassive black holes (SMBHs) found in the centers of many galaxies have been recognized to play a fundamental active role in the cosmological structure formation process. In hierarchical formation scenarios, SMBHs are expected to form binaries following the merger of their host galaxies. If these binaries do not coalesce before the merger with a third galaxy, the formation of a black hole triple system is possible. Numerical simulations of the dynamics of triples within galaxy cores exhibit phases of very high eccentricity (as high as $e \sim 0.99$). During these phases, intense bursts of gravitational radiation can be emitted at orbital periapsis. This produces a gravitational wave signal at frequencies substantially higher than the orbital frequency. The likelihood of detection of these bursts with pulsar timing and the Laser Interferometer Space Antenna ({\it LISA}) is estimated using several population models of SMBHs with masses $\gtrsim 10^7 {\rm M_\odot}$. Assuming a fraction of binaries $\ge 0.1$ in triple system, we find that few to few dozens of these bursts will produce residuals $>1$ ns, within the sensitivity range of forthcoming pulsar timing arrays (PTAs). However, most of such bursts will be washed out in the underlying confusion noise produced by all the other 'standard' SMBH binaries emitting in the same frequency window. A detailed data analysis study would be required to assess resolvability of such sources. Implementing a basic resolvability criterion, we find that the chance of catching a resolvable burst at a one nanosecond precision level is 2-50%, depending on the adopted SMBH evolution model. On the other hand, the probability of detecting bursts produced by massive binaries (masses $\gtrsim 10^7\msun$) with {\it LISA} is negligible.

Recent VLBI (Very Long Baseline Interferometer) observations determined the distances and proper motions of star-forming regions in spiral arms directly. They showed that star-forming regions and young stars have large peculiar motions, as large as 30 km/s with complex structures. Such a large peculiar motion is incompatible with the prediction of the standard theory of quasi-stationary spiral arms. We use a high-resolution, self-consistent N-body+hydrodynamical simulation to explore how the spiral arms are formed and maintained, and how star-forming regions move. We found that arms are not quasi-stationary but transient and recurrent, as suggested in alternative theories of spiral structures. Because of this transient nature of the spiral arms, star-forming regions exhibit a trend of large and complex non-circular motions, which is qualitatively consistent with the VLBI observations. Owing to this large non-circular motion, a kinematically estimated gas map of our Galaxy has a large systematic errors of ~2-3 kpc in the distance from the Sun.

We analyzed the statistics of subhalo abundance of galaxy-sized and giant-galaxy-sized halos formed in a high-resolution cosmological simulation of a 46.5Mpc cube with the uniform mass resolution of $10^6 M_{\odot}$. We analyzed all halos with mass more than $1.5 \times 10^{12}M_{\odot}$ formed in this simulation box. The total number of halos was 125. We found that the subhalo abundance, measured by the number of subhalos with maximum rotation velocity larger than 10% of that of the parent halo, shows large halo-to-halo variations. The results of recent ultra-high-resolution runs fall within the variation of our samples. We found that the concentration parameter and the radius at the moment of the maximum expansion shows fairly tight correlation with the subhalo abundance. This correlation suggests that the variation of the subhalo abundance is at least partly due to the difference in the formation history. Halos formed earlier have smaller number of subhalos at present.

We show that the smoothed particle hydrodynamics (SPH) method, used with individual time-steps in the way described in the literature, cannot handle strong explosion problems correctly. In the individual time-step scheme, particles determine their time-steps essentially from a local Courant condition. Thus they cannot respond to a strong shock, if the pre-shock timescale is too long compared to the shock timescale. This problem is not severe in SPH simulations of galaxy formation with a temperature cutoff in the cooling function at $10^4 {\rm K}$, while it is very dangerous for simulations in which the multiphase nature of the interstellar medium under $10^4 {\rm K}$ is taken into account. A solution for this problem is to introduce a time-step limiter which reduces the time-step of a particle if it is too long compared to the time-steps of its neighbor particles. Thus this kind of time-step constraint is essential for the correct treatment of explosions in high-resolution SPH simulations with individual time-steps.

We performed, for the first time, the simulation of spiral-in of a star cluster formed close to the Galactic center (GC) using a fully self-consistent $N$-body model. In our model, the central super-massive black hole (SMBH) is surrounded by stars and the star cluster. Not only are the orbits of stars and the cluster stars integrated self-consistently, but the stellar evolution, collisions and merging of the cluster stars are also included. We found that an intermediate-mass black hole (IMBH) is formed in the star cluster and stars escaped from the cluster are captured into a 1:1 mean motion resonance with the IMBH. These "Trojan" stars are brought close to the SMBH by the IMBH, which spirals into the GC due to the dynamical friction. Our results show that, once the IMBH is formed, it brings the massive stars to the vicinity of the central SMBH even after the star cluster itself is disrupted. Stars carried by the IMBH form a disk similar to the observed disks and the core of the cluster including the IMBH has properties similar to those of IRS13E, which is a compact assembly of several young stars.

We investigated the evolution of interacting disk galaxies using high-resolution $N$-body/SPH simulations, taking into account the multiphase nature of the interstellar medium (ISM). In our high-resolution simulations, a large-scale starburst occurred naturally at the collision interface between two gas disks at the first encounter, resulting in the formation of star clusters. This is consistent with observations of interacting galaxies. The probability distribution function (PDF) of gas density showed clear change during the galaxy-galaxy encounter. The compression of gas at the collision interface between the gas disks first appears as an excess at $n_{\rm H} \sim 10{\rm cm^{-3}}$ in the PDF, and then the excess moves to higher densities ($n_{\rm H} \gtrsim 100{\rm cm^{-3}}$) in a few times $10^7$ years where starburst takes place. After the starburst, the PDF goes back to the quasi-steady state. These results give a simple picture of starburst phenomena in galaxy-galaxy encounters.

We have performed fully self-consistent $N$-body simulations of star clusters near the Galactic center (GC). Such simulations have not been performed because it is difficult to perform fast and accurate simulations of such systems using conventional methods. We used the Bridge code, which integrates the parent galaxy using the tree algorithm and the star cluster using the fourth-order Hermite scheme with individual timestep. The interaction between the parent galaxy and the star cluster is calculate with the tree algorithm. Therefore, the Bridge code can handle both the orbital and internal evolutions of star clusters correctly at the same time. We investigated the evolution of star clusters using the Bridge code and compared the results with previous studies. We found that 1) the inspiral timescale of the star clusters is shorter than that obtained with "traditional" simulations, in which the orbital evolution of star clusters is calculated analytically using the dynamical friction formula and 2) the core collapse of the star cluster increases the core density and help the cluster survive. The initial conditions of star clusters is not so severe as previously suggested.

Recent high-resolution simulations of the formation of dark-matter halos have shown that the distribution of subhalos is scale-free, in the sense that if scaled by the velocity dispersion of the parent halo, the velocity distribution function of galaxy-sized and cluster-sized halos are identical. For cluster-sized halos, simulation results agreed well with observations. Simulations, however, predicted far too many subhalos for galaxy-sized halos. Our galaxy has several tens of known dwarf galaxies. On the other hands, simulated dark-matter halos contain thousands of subhalos. We have performed simulation of a single large volume and measured the abundance of subhalos in all massive halos. We found that the variation of the subhalo abundance is very large, and those with largest number of subhalos correspond to simulated halos in previous studies. The subhalo abundance depends strongly on the local density of the background. Halos in high-density regions contain large number of subhalos. Our galaxy is in the low-density region. For our simulated halos in low-density regions, the number of subhalos is within a factor of three to that of our galaxy. We argue that the ``missing dwarf problem'' is not a real problem but caused by the biased selection of the initial conditions in previous studies, which were not appropriate for field galaxies.

We performed 3-dimensional N-body/SPH simulations to study how mass resolution and other model parameters such as the star formation efficiency parameter, C* and the threshold density, nth affect structures of the galactic gaseous/stellar disk in a static galactic potential. We employ 10^6 - 10^7 particles to resolve a cold and dense (T < 100 K & n_H > 100 cm^{-3}) phase. We found that structures of the ISM and the distribution of young stars are sensitive to the assumed nth. High-nth models with nth = 100 cm^{-3} yield clumpy multi-phase features in the ISM. Young stars are distributed in a thin disk of which half-mass scale height is 10 - 30 pc. In low-nth models with nth = 0.1 cm^{-3}, the stellar disk is found to be several times thicker, and the gas disk appears smoother than the high-nth models. A high-resolution simulation with high-nth is necessary to reproduce the complex structure of the gas disk. The global properties of the model galaxies in low-nth models, such as star formation histories, are similar to those in the high-nth models when we tune the value of C* so that they reproduce the observed relation between surface gas density and surface star formation rate density. We however emphasize that high-nth models automatically reproduce the relation, regardless of the values of C*. The ISM structure, phase distribution, and distributions of young star forming region are quite similar between two runs with values of C* which differ by a factor of 15. We also found that the timescale of the flow from n_H ~1 cm^{-3} to n_H > 100 cm^{-3} is about 5 times as long as the local dynamical time and is independent of the value of C*. The use of a high-nth criterion for star formation in high-resolution simulations makes numerical models fairy insensitive to the modelling of star formation. (Abridged)

We present sixth- and eighth-order Hermite integrators for astrophysical $N$-body simulations, which use the derivatives of accelerations up to second order ({\it snap}) and third order ({\it crackle}). These schemes do not require previous values for the corrector, and require only one previous value to construct the predictor. Thus, they are fairly easy to implemente. The additional cost of the calculation of the higher order derivatives is not very high. Even for the eighth-order scheme, the number of floating-point operations for force calculation is only about two times larger than that for traditional fourth-order Hermite scheme. The sixth order scheme is better than the traditional fourth order scheme for most cases. When the required accuracy is very high, the eighth-order one is the best. These high-order schemes have several practical advantages. For example, they allow a larger number of particles to be integrated in parallel than the fourth-order scheme does, resulting in higher execution efficiency in both general-purpose parallel computers and GRAPE systems.

In this paper, we investigate the structures of galaxies which either have or have had three BHs using $N$-body simulations, and compare them with those of galaxies with binary BHs. We found that the cusp region of a galaxy which have (or had) triple BHs is significantly larger and less dense than that of a galaxy with binary BHs of the same mass. Moreover, the size of the cusp region depends strongly on the evolution history of triple BHs, while in the case of binary BHs, the size of the cusp is determined by the mass of the BHs. In galaxies which have (or had) three BHs, there is a region with significant radial velocity anisotropy, while such a region is not observed in galaxies with binary BH. These differences come from the fact that with triple BHs the energy deposit to the central region of the galaxy can be much larger due to multiple binary-single BH scatterings. Our result suggests that we can discriminate between galaxies which experienced triple BH interactions with those which did not, through the observable signatures such as the cusp size and velocity anisotropy.

We present performance measurements of direct gravitational N -body simulation on the grid, with and without specialized (GRAPE-6) hardware. Our inter-continental virtual organization consists of three sites, one in Tokyo, one in Philadelphia and one in Amsterdam. We run simulations with up to 196608 particles for a variety of topologies. In many cases, high performance simulations over the entire planet are dominated by network bandwidth rather than latency. With this global grid of GRAPEs our calculation time remains dominated by communication over the entire range of N, which was limited due to the use of three sites. Increasing the number of particles will result in a more efficient execution. Based on these timings we construct and calibrate a model to predict the performance of our simulation on any grid infrastructure with or without GRAPE. We apply this model to predict the simulation performance on the Netherlands DAS-3 wide area computer. Equipping the DAS-3 with GRAPE-6Af hardware would achieve break-even between calculation and communication at a few million particles, resulting in a compute time of just over ten hours for 1 N -body time unit. Key words: high-performance computing, grid, N-body simulation, performance modelling

Black holes and neutron stars present extreme forms of matter that cannot be created as such in a laboratory on Earth. Instead, we have to observe and analyze the experiments that are ongoing in the Universe. The most telling observations of black holes and neutron stars come from dense stellar systems, where stars are crowded close enough to each other to undergo frequent interactions. It is the interplay between black holes, neutron stars and other objects in a dense environment that allows us to use observations to draw firm conclusions about the properties of these extreme forms of matter, through comparisons with simulations. The art of modeling dense stellar systems through computer simulations forms the main topic of this review.

We developed a new direct-tree hybrid N-body algorithm for fully self-consistent N-body simulations of star clusters in their parent galaxies. In such simulations, star clusters need high accuracy, while galaxies need a fast scheme because of the large number of the particles required to model it. In our new algorithm, the internal motion of the star cluster is calculated accurately using the direct Hermite scheme with individual timesteps and all other motions are calculated using the tree code with second-order leapfrog integrator. The direct and tree schemes are combined using an extension of the mixed variable symplectic (MVS) scheme. Thus, the Hamiltonian corresponding to everything other than the internal motion of the star cluster is integrated with the leapfrog, which is symplectic. Using this algorithm, we performed fully self-consistent N-body simulations of star clusters in their parent galaxy. The internal and orbital evolutions of the star cluster agreed well with those obtained using the direct scheme. We also performed fully self-consistent N-body simulation for large-N models ($N=2\times 10^6$). In this case, the calculation speed was seven times faster than what would be if the direct scheme was used.

More than half of stars reside in binary or multiple star systems and many planets have been found in binary systems. From theoretical point of view, however, whether or not the planetary formation proceeds in a binary system is a very complex problem, because secular perturbation from the companion star can easily stir up the eccentricity of the planetesimals and cause high-velocity, destructive collisions between planetesimals. Early stage of planetary formation process in binary systems has been studied by restricted three-body approach with gas drag and it is commonly accepted that accretion of planetesimals can proceed due to orbital phasing by gas drag. However, the gas drag becomes less effective as the planetesimals become massive. Therefore it is still uncertain whether the collision velocity remains small and planetary accretion can proceed, once the planetesimals become massive. We performed {\it N}-body simulations of planetary formation in binary systems starting from massive planetesimals whose size is about 100-500 km. We found that the eccentricity vectors of planetesimals quickly converge to the forced eccentricity due to the coupling of the perturbation of the companion and the mutual interaction of planetesimals if the initial disk model is sufficiently wide in radial distribution. This convergence decreases the collision velocity and as a result accretion can proceed much in the same way as in isolated systems. The basic processes of the planetary formation, such as runaway growth and oligarchic growth and final configuration of the protoplanets are essentially the same in binary systems and single star systems, at least in the late stage where the effect of gas drag is small.

We investigate the evolution of collisionally merged stars with mass of ~100 MSun which might be formed in dense star clusters. We assumed that massive stars with several tens Msun collide typically after ~1Myr of the formation of the cluster and performed hydrodynamical simulations of several collision events. Our simulations show that after the collisions, merged stars have extended envelopes and their radii are larger than those in the thermal equilibrium states and that their interiors are He-rich because of the stellar evolution of the progenitor stars. We also found that if the mass-ratio of merging stars is far from unity, the interior of the merger product is not well mixed and the elemental abundance is not homogeneous. We then followed the evolution of these collision products by a one dimensional stellar evolution code. After an initial contraction on the Kelvin-Helmholtz (thermal adjustment) timescale (~1000-10000 yr), the evolution of the merged stars traces that of single homogeneous stars with corresponding masses and abundances, while the initial contraction phase shows variations which depend on the mass ratio of the merged stars. We infer that, once runaway collisions have set in, subsequent collisions of the merged stars take place before mass loss by stellar winds becomes significant. Hence, stellar mass loss does not inhibit the formation of massive stars with mass of ~1000Msun.

(Simplified) We present a theoretical framework which establishes how the core radius of a star cluster varies with the mass of an assumed central black hole. Our result is that the ratio of core to half-mass radius varies as the 3/4 power of the ratio of the black hole to cluster mass, when the system is well relaxed. The theory compares favourably with a number of simulations of this problem, which extend to black hole masses of order 10% of the cluster mass. Though strictly limited as yet to clusters with stars of equal mass, our conclusion strengthens the view that clusters with large core radii are the most promising candidates in which to find a massive black hole.

We have developed PGPG (Pipeline Generator for Programmable GRAPE), a software which generates the low-level design of the pipeline processor and communication software for FPGA-based computing engines (FBCEs). An FBCE typically consists of one or multiple FPGA (Field-Programmable Gate Array) chips and local memory. Here, the term "Field-Programmable" means that one can rewrite the logic implemented to the chip after the hardware is completed, and therefore a single FBCE can be used for calculation of various functions, for example pipeline processors for gravity, SPH interaction, or image processing. The main problem with FBCEs is that the user need to develop the detailed hardware design for the processor to be implemented to FPGA chips. In addition, she or he has to write the control logic for the processor, communication and data conversion library on the host processor, and application program which uses the developed processor. These require detailed knowledge of hardware design, a hardware description language such as VHDL, the operating system and the application, and amount of human work is huge. A relatively simple design would require 1 person-year or more. The PGPG software generates all necessary design descriptions, except for the application software itself, from a high-level design description of the pipeline processor in the PGPG language. The PGPG language is a simple language, specialized to the description of pipeline processors. Thus, the design of pipeline processor in PGPG language is much easier than the traditional design. For real applications such as the pipeline for gravitational interaction, the pipeline processor generated by PGPG achieved the performance similar to that of hand-written code. In this paper we present a detailed description of PGPG version 1.0.

We study the first 100Myr of the evolution of isolated star clusters initially containing 144179 stars, including 13107 (10%) primordial hard binaries. Our calculations include the effects of both stellar and binary evolution. Gravitational interactions among the stars are computed by direct N-body integration using high precision GRAPE-6 hardware. The evolution of the core radii and central concentrations of our simulated clusters are compared with the observed sample of young (about 100Myr) star clusters in the large Magellanic cloud. Even though our simulations start with a rich population of primordial binaries, core collapse during the early phase of the cluster evolution is not prevented. Throughout the simulations, the fraction of binaries remains roughly constant (about 10%). Due to the effects of mass segregation the mass function of intermediate-mass main-sequence stars becomes as flat as $\alpha=-1.8$ in the central part of the cluster (where the initial Salpeter mass function had $\alpha=-2.35$). About 6--12% of the neutron stars were retained in our simulations; the fraction of retained black holes is 40--70%. In each simulation about three neutron stars become members of close binaries with a main-sequence companion. Such a binary will eventually become an x-ray binary, when the main-sequence star starts to fill its Roche lobe. Black holes are found more frequently in binaries; in each simulated cluster we find about 11 potential x-ray binaries containing a black hole. Abstract abbreviated....

For a rigid model satellite, Chandrasekhar's dynamical friction formula describes the orbital evolution quite accurately, when the Coulomb logarithm is chosen appropriately. However, it is not known if the orbital evolution of a real satellite with the internal degree of freedom can be described by the dynamical friction formula. We performed N-body simulation of the orbital evolution of a self-consistent satellite galaxy within a self-consistent parent galaxy. We found that the orbital decay of the simulated satellite is significantly faster than the estimate from the dynamical friction formula. The main cause of this discrepancy is that the stars stripped out of the satellite are still close to the satellite, and increase the drag force on the satellite through two mechanisms. One is the direct drag force from particles in the trailing tidal arm, a non-axisymmetric force that slows the satellite down. The other is the indirect effect that is caused by the particles remaining close to the satellite after escape. The force from them enhances the wake caused in the parent galaxy by dynamical friction, and this larger wake in turn slows the satellite down more than expected from the contribution of its bound mass. We found these two have comparable effects, and the combined effect can be as large as 20% of the total drag force on the satellite.

In this paper, we describe the performance of an $N$-body simulation of star cluster with 64k stars on a Cray XD1 system with 400 dual-core Opteron processors. A number of astrophysical $N$-body simulations were reported in SCxy conferences. All previous entries for Gordon-Bell prizes used at least 700k particles. The reason for this preference of large numbers of particles is the parallel efficiency. It is very difficult to achieve high performance on large parallel machines, if the number of particles is small. However, for many scientifically important problems the calculation cost scales as $O(N^{3.3})$, and it is very important to use large machines for relatively small number of particles. We achieved 2.03 Tflops, or 57.7% of the theoretical peak performance, using a direct $O(N^2)$ calculation with the individual timestep algorithm, on 64k particles. The best efficiency previously reported on similar calculation with 64K or smaller number of particles is 12% (9 Gflops) on Cray T3E-600 with 128 processors. Our implementation is based on highly scalable two-dimensional parallelization scheme, and low-latency communication network of Cray XD1 turned out to be essential to achieve this level of performance.

The method of choice for integrating the equations of motion of the general N-body problem has been to use an individual time step scheme. For the sake of efficiency, block time steps have been the most popular, where all time step sizes are smaller than a maximum time step size by an integer power of two. We present the first successful attempt to construct a time-symmetric integration scheme, based on block time steps. We demonstrate how our scheme shows a vastly better long-time behavior of energy errors, in the form of a random walk rather than a linear drift. Increasing the number of particles makes the improvement even more pronounced.

We investigated the formation and evolution of CDM subhalos in galaxy-sized and cluster-sized halos by means of N-body simulations. Our aim is to make clear what the ``dwarf galaxy problem'' is. It has been argued that the number of subhalos in simulated galaxy-sized halos is too large compared to the observed number of dwarfs in the local group, while that in cluster-sized halos is consistent with observed number of galaxies in clusters such as the Virgo cluster. We simulated nine halos with several different mass resolutions and physical scales. We found that the dependence of the cumulative number of subhalos N_c on their maximum circular velocity V_c is given by N_c \propto V_c^-3, down to the reliability limit, independent of the mass of the main halo. This implies that simulations for cluster-sized halos give too many halos with V_c ~ 140km/s or less. Previous comparisons of cluster-sized halos gave much smaller number of subhalos in this regime simply because of their limited resolution. Our result implies that any theory which attempts to resolve the missing dwarf problem should also explain the discrepancy of the simulation and observation in cluster-sized halos.

In order to examine the ``giant impact hypothesis'' for the Moon formation, we run the first grid-based, high-resolution hydrodynamic simulations for an impact between proto-Earth and a proto-planet. The spatial resolution for the impact-generated disk is greatly improved from previous particle-based simulations. This allows us to explore fine structures of a circumterrestrial debris disk and its long-term evolution. We find that in order to form a debris disk from which a lunar-sized satellite can be accumulated, the impact must result in a disk of mostly liquid or solid debris, where pressure is not effective, well before the accumulation process starts. If the debris is dominated by vapor gas, strong spiral shocks are generated, and therefore the circumterrestrial disk cannot survive more than several days. This suggests that there could be an appropriate mass range for terrestrial planets to harbor a large moon as a result of giant impacts, since vaporization during an impact depends on the impact energy.

We simulate the inner 100pc of the Milky-Way Galaxy to study the formation and evolution of the population of star clusters and intermediate mass black holes. For this study we perform extensive direct N-body simulations of the star clusters which reside in the bulge, and of the inner few tenth of parsecs of the super massive black hole in the Galactic center. In our N-body simulations the dynamical friction of the star cluster in the tidal field of the bulge are taken into account via (semi)analytic soluations. The N-body calculations are used to calibrate a (semi)analytic model of the formation and evolution of the bulge. We find that about 10% of the clusters born within 100pc of the Galactic center undergo core collapse during their inward migration and form intermediate-mass black holes (IMBHs) via runaway stellar merging. After the clusters dissolve, these IMBHs continue their inward drift, carrying a few of the most massive stars with them. We predict that region within about 10 parsec of the SMBH is populated by about 50IMBHs of some 1000Msun. Several of these are expected to be accompanied still by some of the most massive stars from the star cluster. We also find that within a few milliparsec of the SMBH there is a steady population of several IMBHs. This population drives the merger rate between IMBHs and the SMBH at a rate of about one per 10Myr, sufficient to build the accumulate majority of mass of the SMBH. Mergers of IMBHs with SMBHs throughout the universe are detectable by LISA, at a rate of about two per week.

Recent observations (Burnes2002,Veillet2002,Margot2002a) have revealed an unexpectedly high binary fraction among the Trans-Neptunian Objects (TNOs) that populate the Kuiper Belt. The TNO binaries are strikingly different from asteroid binaries in four respects (Veillet2002): their frequency is an order of magnitude larger, the mass ratio of their components is closer to unity, and their orbits are wider and highly eccentric. Two explanations have been proposed for their formation, one assuming large numbers of massive bodies (Weidenschilling2002), and one assuming large numbers of light bodies (Goldreich2002). We argue that both assumptions are unwarranted, and we show how TNO binaries can be produced from a modest number of intermediate-mass bodies of the type predicted by the gravitational instability theory for the formation of planetesimals (Goldreich and Ward1973). We start with a TNO binary population similar to the asteroid binary population, but subsequently modified by three-body exchange reactions, a process that is far more efficient in the Kuiper belt, because of the much smaller tidal perturbations by the Sun. Our mechanism can naturally account for all four characteristics that distinguish TNO binaries from main-belt asteroid binaries.

Resent observations and theoretical interpretations suggest that IMBHs (intermediate-mass black hole) are formed in the centers of young and compact star clusters born close to the center of their parent galaxy. Such a star cluster would sink toward the center of the galaxy, and at the same time stars are stripped out of the cluster by the tidal field of the parent galaxy. We investigated the orbital evolution of the IMBH, after its parent cluster is completely disrupted by the tidal field of the parent galaxy, by means of large-scale N-body simulations. We constructed a model of the central region of our galaxy, with an SMBH (supermassive black hole) and Bahcall-Wolf stellar cusp, and placed an IMBH in a circular orbit of radius 0.086pc. The IMBH sinks toward the SMBH through dynamical friction, but dynamical friction becomes ineffective when the IMBH reached the radius inside which the initial stellar mass is comparable to the IMBH mass. This is because the IMBH kicks out the stars. This behavior is essentially the same as the loss-cone depletion observed in simulations of massive SMBH binaries. After the evolution through dynamical friction stalled, the eccentricity of the orbit of the IMBH goes up, resulting in the strong reduction in the merging timescale through gravitational wave radiation. Our result indicates that the IMBHs formed close to the galactic center can merge with the central SMBH in short time. The number of merging events detectable with DECIGO is estimated to be around 50 per year. Event rate for LISA would be similar or less, depending on the growth mode of IMBHs.

Recent X-ray observations and theoretical modelling have made it plausible that some ultraluminous X-ray sources (ULX) are powered by intermediate-mass black holes (IMBHs). N-body simulations have also shown that runaway merging of stars in dense star clusters is a way to form IMBHs. In the present paper we have performed N-body simulations of young clusters such as MGG-11 of M82 in which IMBHs form through runaway merging. We took into account the effect of tidal heating of stars by the IMBH to study the tidal capture and disruption of stars by IMBHs. Our results show that the IMBHs have a high chance of capturing stars through tidal heating within a few core relaxation times and we find that 1/3 of all runs contain a ULX within the age limits of MGG-11, a result consistent with the fact that a ULX is found in this galaxy. Our results strengthen the case for some ULX being powered by intermediate-mass black holes.

We present the result of $N$-body simulations of dynamical evolution of triple massive blackhole (BH) systems in galactic nuclei. We found that in most cases two of the three BHs merge through gravitational wave (GW) radiation in the timescale much shorter than the Hubble time, before ejecting one BH through a slingshot. In order for a binary BH to merge before ejecting out the third one, it has to become highly eccentric since the gravitational wave timescale would be much longer than the Hubble time unless the eccentricity is very high. We found that two mechanisms drive the increase of the eccentricity of the binary. One is the strong binary-single BH interaction resulting in the thermalization of the eccentricity. The second is the Kozai mechanism which drives the cyclic change of the inclination and eccentricity of the inner binary of a stable hierarchical triple system. Our result implies that many of supermassive blackholes are binaries.

Time-symmetric integration schemes share with symplectic schemes the property that their energy errors show a much better behavior than is the case for generic integration schemes. Allowing adaptive time steps typically leads to a loss of symplecticity. In contrast, time symmetry can be easily maintained, at least for a continuous choice of time step size. In large-scale N-body simulations, however, one often uses block time steps, where all time steps are forced to take on values as powers of two. This greatly facilitates parallelization, and hence code efficiency. Straightforward implementation of time-symmetry, translated to block time steps, faces significant hurdles. For example, iteration can lead to oscillatory behavior, and even when such behavior is suppressed, energy errors show a linear drift in time. We present an approach that circumvents these problems.

The main performance bottleneck of gravitational N-body codes is the force calculation between two particles. We have succeeded in speeding up this pair-wise force calculation by factors between two and ten, depending on the code and the processor on which the code is run. These speedups were obtained by writing highly fine-tuned code for x86_64 microprocessors. Any existing N-body code, running on these chips, can easily incorporate our assembly code programs. In the current paper, we present an outline of our overall approach, which we illustrate with one specific example: the use of a Hermite scheme for a direct N^2 type integration on a single 2.0 GHz Athlon 64 processor, for which we obtain an effective performance of 4.05 Gflops, for double precision accuracy. In subsequent papers, we will discuss other variations, including the combinations of N log N codes, single precision implementations, and performance on other microprocessors.

We describe a modified SIMD architecture suitable for single-chip integration of a large number of processing elements, such as 1,000 or more. Important differences from traditional SIMD designs are: a) The size of the memory per processing elements is kept small. b) The processors are organized into groups, each with a small buffer memory. Reduction operation over the groups is done in hardware. The first change allows us to integrate a very large number of processing elements into a single chip. The second change allows us to achieve a close-to-peak performance for many scientific applications like particle-based simulations and dense-matrix operations.

In this paper, we describe the design and performance of GRAPE-6A, a special-purpose computer for gravitational many-body simulations. It was designed to be used with a PC cluster, in which each node has one GRAPE-6A. Such configuration is particularly effective in running parallel tree algorithm. Though the use of parallel tree algorithm was possible with the original GRAPE-6 hardware, it was not very cost-effective since a single GRAPE-6 board was still too fast and too expensive. Therefore, we designed GRAPE-6A as a single PCI card to minimize the reproduction cost and optimize the computing speed. The peak performance is 130 Gflops for one GRAPE-6A board and 3.1 Tflops for our 24 node cluster. We describe the implementation of the tree, TreePM and individual timestep algorithms on both a single GRAPE-6A system and GRAPE-6A cluster. Using the tree algorithm on our 16-node GRAPE-6A system, we can complete a collisionless simulation with 100 million particles (8000 steps) within 10 days.

Close encounters and physical collisions between stars in young dense clusters can result in new channels for stellar evolution, and may lead to the formation of very massive stars and black holes via runaway merging. We present some details of this process, using the results of N-body simulations and simple analytical estimates to place limits on the cluster parameters for which it expected to occur. For small clusters, the mass of the runaway is effectively limited by the total number of high-mass stars in the system. For larger clusters, the runaway mass is determined by the fraction of stars that can mass-segregate to the cluster core while still on the main sequence. In typical cases, the result is in the range commonly cited for intermediate-mass black holes. This mechanism may therefore have important implications for the formation of massive black holes and black-hole binaries in dense cluster cores.

It has been assumed that intermediate-mass black holes (IMBHs) in globular clusters can only reside in the most centrally concentrated clusters, with a so-called `core-collapsed' density profile. While this would be a natural guess, it is in fact wrong. We have followed the evolution of star clusters containing IMBHs with masses between 125 \le M_{BH} \le 1000 M_{\odot} through detailed N-body simulations, and we find that a cluster with an IMBH, in projection, appears to have a relatively large `core' with surface brightness only slightly rising toward the center. This makes it highly unlikely that any of `core-collapsed' clusters will harbor an IMBH. On the contrary, the places to look for an IMBH are those clusters that can be fitted well by medium-concentration King models. The velocity dispersion of the visible stars in a globular cluster with an IMBH is nearly constant well inside the apparent core radius. For a cluster of mass M_C containing an IMBH of mass M_{BH}, the influence of the IMBH becomes significant only at a fraction 2.5 M_{BH}/M_C of the half-mass radius, deep within the core, where it will affect only a small number of stars. In conclusion, observational detection of an IMBH may be possible, but will be challenging.

We have investigated whether or not a tidal stripping scenario can reproduce the observed surface-brightness profile of omega Centauri using N-body simulations. Assuming that the progenitor of omega Centauri is a dwarf elliptical galaxy, we model it with a King model with a core radius being the same as that of omega Centauri. A dark matter halo of the dwarf is not taken into account. We consider two different models of the Milky Way potential: a singular isothermal sphere and a three-component model. The progenitor dwarf is expressed as an N-body system, which orbits in the fixed Galactic potential. The dwarf lost more than 90 per cent of its mass during the first few pericenter passages. Thereafter, the mass remains practically constant. The final surface-density profile is in good agreement with the observational data on omega Centauri, if the pericenter distance of the orbit of the progenitor dwarf is around 500 pc. This value is within the error bar of the current proper motion data on omega Centauri and Galactic parameters. Although our simulation is limited to a King-like progenitor dwarf without a dark matter halo, it strongly suggests that the current density profile of omega Centauri is nicely reproduced by a tidal stripping scenario, in other words, that omega Centauri can plausibly be identified with a stripped dwarf elliptical.

We have followed the evolution of multi-mass star clusters containing massive central black holes through collisional N-body simulations done on GRAPE6. Each cluster is composed of between 16,384 to 131,072 stars together with a black hole with an initial mass of M_BH=1000 Msun. We follow the evolution of the clusters under the combined influence of two-body relaxation, stellar mass-loss and tidal disruption of stars. The (3D) mass density profile follows a power-law distribution \rho \sim r^{-\alpha} with slope \alpha=1.55. This leads to a constant density profile of bright stars in projection, which makes it highly unlikely that core collapse clusters contain intermediate-mass black holes (IMBHs). Instead globular clusters containing IMBHs can be fitted with standard King profiles. The disruption rate of stars is too small to form an IMBH out of a M_BH \approx 50 Msun progenitor black hole, unless a cluster starts with a central density significantly higher than what is seen in globular clusters. Kinematical studies can reveal 1000 Msun IMBHs in the closest clusters. IMBHs in globular clusters are only weak X-ray sources since the tidal disruption rate of stars is low and the star closest to the IMBH is normally another black hole. For globular clusters, dynamical evolution can push compact stars near the IMBH to distances small enough that they become detectable through gravitational radiation. If 10% of all globular clusters contain IMBHs, extragalactic globular clusters could be one of the major sources for {\it LISA}. (abridged)

In this paper we report results of collisional N-body simulations of the dynamical evolution of equal-mass star clusters containing a massive central black hole. Each cluster is composed of between 5,000 to 180,000 stars together with a central black hole which contains between 0.2% to 10% of the total cluster mass. We find that for large enough black hole masses, the central density follows a power-law distribution with slope \rho \sim r^{-1.75} inside the radius of influence of the black hole, in agreement with predictions from earlier Fokker Planck and Monte Carlo models. The tidal disruption rate of stars is within a factor of two of that derived in previous studies. It seems impossible to grow an intermediate-mass black hole (IMBH) from a M \le 100 Msun progenitor in a globular cluster by the tidal disruption of stars, although M = 10^3 Msun IMBHs can double their mass within a Hubble time in dense globular clusters. The same is true for the supermassive black hole at the centre of the Milky Way. Black holes in star clusters will feed mainly on stars tightly bound to them and the re-population of these stars causes the clusters to expand, reversing core-collapse without the need for dynamically active binaries. Close encounters of stars in the central cusp also lead to an increased mass loss rate in the form of high-velocity stars escaping from the cluster. A companion paper will extend these results to the multi-mass case.

We have followed the evolution of multi-mass star clusters containing massive central black holes by N-body simulations on the GRAPE6 computers of Tokyo University. We find a strong cluster expansion and significant structural changes of the clusters. Star clusters with IMBHs have power-law density profiles $rho ~ r^{-alpha}$ with slopes $alpha=1.55$ inside the influence sphere of the central black hole. This leads to a constant density profile of bright stars in projection, which rules out the presence of intermediate mass black holes in core collapse clusters. If the star clusters are surrounded by a tidal field, a central IMBH speeds up the destruction of the cluster until a remnant of a few hundred stars remains, which stays bound to the IMBH for a long time. We also discuss the efficiency of different detection mechanisms for finding IMBHs in star clusters.

A luminous X-ray source is associated with a cluster (MGG-11) of young stars \~200pc from the center of the starburst galaxy M82. The properties of the X-ray source are best explained by a black hole with a mass of at least 350Msun, which is intermediate between stellar-mass and supermassive black holes. A nearby but somewhat more massive star cluster (MGG-9) shows no evidence of such an intermediate mass black hole, raising the issue of just what physical characteristics of the clusters can account for this difference. Here we report numerical simulations of the evolution and the motions of stars within the clusters, where stars are allowed to mergers with each other. We find that for MGG-11 dynamical friction leads to the massive stars sinking rapidly to the center of the cluster to participate in a runaway collision, thereby producing a star of 800-3000Msun, which ultimately collapses to an black hole of intermediate mass. No such runaway occurs in the cluster MGG-9 because the larger cluster radius leads to a mass-segregation timescale a factor of five longer than for MGG-11.

Recent observations have revealed an unexpectedly high binary fraction among the Trans-Neptunian Objects (TNOs) that populate the Kuiper-belt. The discovered binaries have four characteristics they comprise a few percent of the TNOs, the mass ratio of their components is close to unity, their internal orbits are highly eccentric, and the orbits are more than 100 times wider than the primary's radius. In contrast, theories of binary asteroid formation tend to produce close, circular binaries. Therefore, a new approach is required to explain the unique characteristics of the TNO binaries. Two models have been proposed. Both, however, require extreme assumptions on the size distribution of TNOs. Here we show a mechanism which is guaranteed to produces binaries of the required type during the early TNO growth phase, based on only one plausible assumption, namely that initially TNOs were formed through gravitational instabilities of the protoplanetary dust layer.

In this paper, we describe the architecture and performance of the GRAPE-6 system, a massively-parallel special-purpose computer for astrophysical $N$-body simulations. GRAPE-6 is the successor of GRAPE-4, which was completed in 1995 and achieved the theoretical peak speed of 1.08 Tflops. As was the case with GRAPE-4, the primary application of GRAPE-6 is simulation of collisional systems, though it can be used for collisionless systems. The main differences between GRAPE-4 and GRAPE-6 are (a) The processor chip of GRAPE-6 integrates 6 force-calculation pipelines, compared to one pipeline of GRAPE-4 (which needed 3 clock cycles to calculate one interaction), (b) the clock speed is increased from 32 to 90 MHz, and (c) the total number of processor chips is increased from 1728 to 2048. These improvements resulted in the peak speed of 64 Tflops. We also discuss the design of the successor of GRAPE-6.

We present the result of large-scale N-body simulations of the stellar-dynamical evolution of a massive black-hole binary at the center of a spherical galaxy. We focus on the dependence of the hardening rate on the relaxation timescale of the parent galaxy. A simple theoretical argument predicts that a binary black hole creates the ``loss cone'' around it. Once the loss cone is formed, the hardening rate is determined by the rate at which field stars diffuse into the loss cone. Therefore the hardening timescale becomes proportional to the relaxation timescale. Recent N-body simulations, however, have failed to confirm this theory and various explanations have been proposed. By performing simulations with sufficiently large N (up to $10^6$) for sufficiently long time, we found that the hardening rate does depend on N. Our result is consistent with the simple theoretical prediction that the hardening timescale is proportional to the relaxation timescale. This dependence implies that most massive black hole binaries are unlikely to merge within the Hubble time through interaction with field stars and gravitational wave radiation alone.

I overview the current understanding of the evolution of massive black hole (MBH) binaries in the center of the host stellar system. One of the main questions is whether the stellar dynamical effect can make the MBH binary hard enough that they can merge through gravitational wave radiation. So far, theories and numerical simulations suggested otherwize, since the hardening timescale becomes very long once "loss cone" is depleted. I'll present the result of recent simulations on this hardening timescale, and discuss its implication on the formation history of massive black holes.

We investigate the structure of the dark matter halo formed in the cold dark matter scenarios by N-body simulations with parallel treecode on GRAPE cluster systems. We simulated 8 halos with the mass of $4.4\times 10^{14}M_{\odot}$ to $1.6\times 10^{15}M_{\odot}$ in the SCDM and LCDM model using up to 30 million particles. With the resolution of our simulations, the density profile is reliable down to 0.2 percent of the virial radius. Our results show that the slope of inner cusp within 1 percent virial radius is shallower than -1.5, and the radius where the shallowing starts exhibits run-to-run variation, which means the innermost profile is not universal.

We present a comparison between the observational data on the kinematical structure of G1 in M31, obtained with the Hubble WFPC2 and STIS instruments, and the results of dynamical simulations carried out using the special-purpose computer GRAPE-6. We have obtained good fits for models starting with single cluster King-model initial conditions and even better fits when starting our simulations with a dynamically constructed merger remnant of two star clusters. In the latter case, the results from our simulations are in excellent agreement with the observed profiles of luminosity, velocity dispersion, rotation, and ellipticity. We obtain a mass-to-light ratio of $M/L = 4.0 \pm 0.4$ and a total cluster mass of $M=(8\pm 1)\times 10^6 M_\odot$. Given that our dynamical model can fit all available observational data very well, there seems to be no need to invoke the presence of an intermediate-mass black hole in the center of G1.

Recent observations (Burnes et al 2002,Veillet et al 2002, Margot et al 2002a) have revealed an unexpectedly high binary fraction among the Trans-Neptunian Objects (TNOs) that populate the Kuiper Belt. The TNO binaries are strikingly different from asteroid binaries in four respects (Veillet et al 2002): their frequency is an order of magnitude larger, the mass ratio of their components is closer to unity, and their orbits are wider and highly eccentric. Two explanations have been proposed for their formation, one assuming large numbers of massive bodies (Weidenschilling 2002), and one assuming large numbers of light bodies (Goldreich et al 2002). We argue that both assumptions are unwarranted, and we show how TNO binaries can be produced from a modest number of intermediate-mass bodies of the type predicted by the gravitational instability theory for the formation of planetesimals (Goldreich and Ward 1973). We start with a TNO binary population similar to the asteroid binary population, but subsequently modified by three-body exchange reactions, a process that is far more efficient in the Kuiper belt, because of the much smaller tidal perturbations by the Sun. Our mechanism can naturally account for all four characteristics that distinguish TNO binaries from main-belt asteroid binaries.

This paper considers the phenomenon of deep core collapse in collisional stellar systems, with stars of equal mass. The collapse takes place on some multiple, $\xi ^{-1}$, of the central relaxation time, and produces a density profile in which $\rho \propto r^{-\alpha}$, where $\alpha$ is a constant. The parameters $\alpha$ and $\xi$ have usually been determined from simplified models, such as gas and Fokker-Planck models, often with the simplification of isotropy. Here we determine the parameters directly from N-body simulations carried out using the newly completed GRAPE-6.

We have modeled in detail the evolution of rich open star clusters such as the Pleiades, Praesepe and Hyades, using simulations that include stellar dynamics as well as the effects of stellar evolution. The dynamics is modeled via direct N-body integration, while the evolution of single stars and binaries is followed through the use of fitting formulae and recipes. The feedback of stellar and binary evolution on the dynamical evolution of the stellar system is taken into account self-consistently. Our model clusters dissolve in the tidal field of the Galaxy in a time span on the order of a billion years. The rate of mass loss is rather constant, about 1Msun per million years. The binary fraction at first is nearly constant in time, then increases slowly near the end of a cluster's lifetime. For clusters which are more than about 100 million years old the fractions of stars in the form of binaries, giants and collision products in the inner few core radii are considerably higher than in the outer regions, beyond the cluster's half mass radius. When stars with masses >2Msun escape from the cluster, they tend to do so with velocities higher than average.... <Abstract abbreviated, read manuscript for rest>

We present a detailed comparison between the latest observational data on the kinematical structure of the core of M15, obtained with the Hubble STIS and WFPC2 instruments, and the results of dynamical simulations carried out using the special-purpose GRAPE-6 computer. The observations imply the presence of a significant amount of dark matter in the cluster core. In our dynamical simulations, neutron stars and/or massive white dwarfs concentrate to the center through mass segregation, resulting in a sharp increase in $M/L$ toward the center. While consistent with the presence of a central black hole, the Hubble data can also be explained by this central concentration of stellar-mass compact objects. The latter interpretation is more conservative, since such remnants result naturally from stellar evolution, although runaway merging leading to the formation of a black hole may also occur for some range of initial conditions. We conclude that no central massive object is required to explain the observational data, although we cannot conclusively exclude such an object at the level of $\sim500-1000$ solar masses. Our findings are unchanged when we reduce the assumed neutron-star retention fraction in our simulations from 100% to 0%.

We report results of a large set of N-body calculations aimed to study the evolution of multi-mass star clusters in external tidal fields. Our clusters start with the same initial mass-functions, but varying particle numbers, orbital types and density profiles. Our main focus is to study how the stellar mass-function and other cluster parameters change under the combined influence of stellar evolution, two-body relaxation and the external tidal field. We find that the lifetimes of star clusters moving on similar orbits scale as T sim T_RH^x, where T_RH is the relaxation time, and the exponent x depends on the initial concentration of the cluster and is around x approx 0.75. From the results for the lifetimes, we predict that between 53% to 67% of all galactic globular clusters will be destroyed within the next Hubble time. Low-mass stars are preferentially lost and the depletion is strong enough to turn initially increasing mass-functions into mass-functions which decrease towards the low-mass end. The details of this depletion are insensitive to the starting condition of the cluster. The preferential depletion of low-mass stars from star clusters leads to a decrease of their mass-to-light ratios except for a short period close to final dissolution. The fraction of compact remnants is increasing throughout the evolution and they are more strongly concentrated towards the cluster cores than main-sequence stars. For a sample of galactic globular clusters with well observed parameters, we find a correlation between the observed slope of the mass-function and the lifetimes predicted by us. It seems possible that galactic globular clusters started with a mass-function similar to what one observes for the average mass-function of the galactic disc and bulge. (Abridged)

We summarize the main results from MODEST-1, the first workshop on MOdeling DEnse STellar systems. Our goal is to go beyond traditional population synthesis models, by introducing dynamical interactions between single stars, binaries, and multiple systems. The challenge is to define and develop a software framework to enable us to combine in one simulation existing computer codes in stellar evolution, stellar dynamics, and stellar hydrodynamics. With this objective, the workshop brought together experts in these three fields, as well as other interested astrophysicists and computer scientists. We report here our main conclusions, questions and suggestions for further steps toward integrating stellar evolution and stellar (hydro)dynamics.

We investigated the orbital evolution of satellite galaxies using numerical simulations. It has been long believed that the orbit suffers circularization due to the dynamical friction from the galactic halo during orbital decay. This circularization was confirmed by numerous simulations where dynamical friction is added as external force. However, some of the resent N-body simulations demonstrated that circularization is much slower than expected from approximate calculations. We found that the dominant reason for this discrepancy is the assumption that Coulomb logarithm log Lambda is constant, which has been used in practically all recent calculations. Since the size of the satellite is relatively large, accurate determination of the outer cutoff radius is crucial to obtain good estimate for the dynamical friction. An excellent agreement between $N$-body simulations and approximate calculations was observed when the outer cutoff radius is taken to be the distance of the satellite to the center of the galaxy. When satellite is at the perigalacticon, the distance to the center is smaller and therefore log Lambda becomes smaller. As a result, the dynamical friction becomes less effective. We apply our result to the Large Magellanic Cloud. We found that the expected lifetime of the LMC is twice as long as that would be predicted with previous calculations. Previous study predicts that the LMC will merge into the Milky Way after 7 G years, while we found that the merging will take place after 14 G years from now. Our result suggests that generally satellites formed around a galaxy have longer lifetime than previous estimates.

Challenging stellar dynamical problems, such as the study of gravothermal oscillations in star clusters, have in the past initiated the very successful building of GRAPE special purpose computers. It is discussed, that present day tasks such as the formation and evolution of galactic nuclei with one or more massive black holes and the coupled stellar and gas dynamical processes in the formation of nuclei and star clusters, demand a new kind of hybrid architecture, using both GRAPE and a reconfigurable logics board called RACE. For such a system we have developed first implementations and floating point performance studies in the case of the SPH algorithm (smoothed particle hydrodynamics), which will be of great advantage for SPH modelling and also for direct $N$-body simulations using the more efficient Ahmad-Cohen neighbour schemes.

We investigated the dynamical evolution of clusters of galaxies in virial equilibrium using Fokker-Planck models and self-consistent N-body models. In particular, we focused on the growth of a common halo, which is a cluster-wide halo formed by matter stripped from galaxies, and the development of a central density cusp. The Fokker-Planck models include the effects of two-body gravitational encounters both between galaxies and between galaxies and common halo particles. The effects of tidal mass stripping from the galaxies due to close galaxy-galaxy encounters and accompanying dissipation of the orbital kinetic energies of the galaxies were also taken into account in the Fokker-Planck models. We find that the results of the Fokker-Planck models are in excellent agreement with those of the N-body models regarding the growth of the common halo mass and the evolution of the cluster density profiles. In the central region of the cluster, a shallow density cusp, approximated by $\rho (r) \propto r^{-\alpha}$ ($\alpha \sim$ 1), develops. This shallow cusp results from the combined effects of two-body relaxation and tidal stripping. The cusp steepness, $\alpha$, weakly depends on the relative importance of the tidal stripping. When the effect of stripping is important, the central velocity dispersion decreases as the central density increases and, consequently, a shallow ($\alpha <2$) cusp is formed. In the limit of no stripping, usual gravothermal core collapse occurs, i.e. the central velocity dispersion increases as the central density increases with a steep ($\alpha >2$) cusp left. We conclude from our consideration of the origin of the cusp demonstrated here that shallow cusps should develop in real galaxy clusters.

During the last couple of years, observers have started to make plans for a Virtual Observatory, as a federation of existing data bases, connected through levels of software that enable rapid searches, correlations, and various forms of data mining. We propose to extend the notion of a Virtual Observatory by adding archives of simulations, together with interactive query and visualization capabilities, as well as ways to simulate observations of simulations in order to compare them with observations. For this purpose, we have already organized two small workshops, earlier in 2001, in Tucson and Aspen. We have also provided concrete examples of theory data, designed to be federated with a Virtual Observatory. These data stem from a project to construct an archive for our large-scale simulations using the GRAPE-6 (a 32-Teraflops special purpose computer for stellar dynamics). We are constructing interfaces by which remote observers can observe these simulations. In addition, these data will enable detailed comparisons between different simulations.

We study the evolution and observability of young, compact star clusters near the Galactic center, such as the Arches and Quintuplet systems. The clusters are modeled by integrating the equations of motion of all stars while accounting for the internal evolution of stars and binaries, as well as the effect of the Galactic tidal field. We find that clusters within 150 pc of the Galactic center dissolve within ~55 Myr, but their projected densities drop below the background density in the direction of the Galactic center within only a few Myr, effectively making these clusters undetectable after that time. Detailed observations of the Arches cluster, taken at face value, suggest that its mass function is unusually flat and that the cluster contains an overabundance of stars more massive than 20 Msun. Our dynamical analysis, however, shows that the observed characteristics of the Arches cluster are consistent with a perfectly normal initial mass function. The observed anomalies are then caused by a combination of observational selection effects and the dynamical evolution of the cluster. We calibrate the current parameters of the Arches cluster using a normal initial mass function and conclude that the cluster is more massive than 40000 Msun, has a half mass radius of about 0.35 pc and is located between 50 and 90 pc from the Galactic center.

We present a novel, highly efficient algorithm to parallelize O(N^2) direct summation method for N-body problems with individual timesteps on distributed-memory parallel machines such as Beowulf clusters. Previously known algorithms, in which all processors have complete copies of the N-body system, has the serious problem that the communication-computation ratio increases as we increase the number of processors, since the communication cost is independent of the number of processors. In the new algorithm, p processors are organized as a $\sqrt{p}\times \sqrt{p}$ two-dimensional array. Each processor has $N/\sqrt{p}$ particles, but the data are distributed in such a way that complete system is presented if we look at any row or column consisting of $\sqrt{p}$ processors. In this algorithm, the communication cost scales as $N /\sqrt{p}$, while the calculation cost scales as $N^2/p$. Thus, we can use a much larger number of processors without losing efficiency compared to what was practical with previously known algorithms.

We investigate the structure of the dark matter halo formed in three different cold dark matter scenarios. We performed N-body simulations of formation of 13 cluster-sized halos. In all runs, density cusps proportional to $r^{-1.5}$ developed at the center. This result was independent of the cosmological models we simulated. We could not reproduce the cusp shallower than $r^{-1.5}$, which was obtained in some of previous studies. We also found that in all runs the density structure evolves in a self-similar way, even in $\Omega\ne 1$ universes. These results show that the formation of structural form is a process decoupled from a background cosmology.

Observations of stellar kinematics, gas dynamics and masers around galactic nuclei have now firmly established that many galaxies host central supermassive black holes (SMBHs) with masses in the range $10^6 \sim 10^9$M$_{\odot}$. However, how these SMBHs formed is not well understood. One reason for this situation is the lack of observations of intermediate-mass BHs (IMBHs), which could bridge the gap between stellar-mass BHs and SMBHs. Recently, this missing link (i.e., an IMBH) has been found in observations made by the ASCA and the Chandra of the central region of the starburst galaxy M82 \citep{MT99, PG99, MT01, Ka01}. Subsequent observations by SUBARU have revealed that this IMBH apparently coincides with a young compact star cluster. Based on these findings, we suggest a new formation scenario for SMBHs. In this scenario, IMBHs first form in young compact star clusters through runaway merging of massive stars. While these IMBHs are forming, the host star clusters sink toward the galactic nucleus through dynamical friction, and upon evaporation deposit their IMBHs near the galactic center. The IMBHs then form binaries and eventually merge via gravitational radiation, forming an SMBH.

We investigate the structure of the dark matter halo formed in the cold dark matter scenario using $N$-body simulations. We simulated 12 halos with the mass of $6.6\times 10^{11}M_{\odot}$ to $8.0\times 10^{14}M_{\odot}$. In almost all runs, the halos have density cusps proportional to $r^{-1.5}$ developed at the center, which is consistent with the results of recent high-resolution calculations. The density structure evolves in a self-similar way, and is universal in the sense that it is independent of the halo mass and initial random realization of density fluctuation. The density profile is in good agreement with the profile proposed by Moore et al. (1999), which has central slope proportional to $r^{-1.5}$ and outer slope proportional to $r^{-3}$. The halo grows through repeated accretion of diffuse smaller halos. We argue that the cusp is understood as a convergence slope for the accretion of tidally disrupted matter.

In dense clusters a bewildering variety of interactions between stars can be observed, ranging from simple encounters to collisions and other mass-transfer encounters. With faster and special-purpose computers like GRAPE, the amount of data per simulation is now exceeding 1TB. Visualization of such data has now become a complex 4D data-mining problem, combining space and time, and finding interesting events in these large datasets. We have recently starting using the virtual reality simulator, installed in the Hayden Planetarium in the American Museum for Natural History, to tackle some of these problem. This work (http://www.astro.umd.edu/nemo/amnh/) reports on our first ``observations'', modifications needed for our specific experiments, and perhaps field ideas for other fields in science which can benefit from such immersion. We also discuss how our normal analysis programs can be interfaced with this kind of visualization.

We study the occurrence of physical collisions between stars in young and compact star cluster. The calculations are performed on the GRAPE-4 with the starlab software environment which include the dynamical evolution and the nuclear evolution of all stars and binaries. The selection of the initial conditions is based on existing and well observed star clusters, such as R136 in the 30 Doradus region in the Large Magellanic Cloud and the Arches and Quintuplet star clusters in the vicinity of the Galactic center. Collisions between stars occurred rather frequently in our models. At any time a single star dominates the collision history of the system. The collision rate of this runaway merger scales with the initial relaxation time of the cluster and is independent on other cluster parameters, such as the initial mass function or the initial density profile of the cluster. Subsequent encounters result in a steady grow in mass of the coagulating star, until it escapes or explodes in a supernova. The collision rate in these models is about 0.00022 collisions per star per Myr for a cluster with an initial relaxation time of 1 Myr.

We investigated the evolution of clusters of galaxies using self-consistent N-body simulations. We varied the initial model of galaxies and clusters, and studied the dependence of evolution on initial conditions. We found that the growth rate of the common halo depends only weakly on galaxy models. On the other hand, the growth rate depends strongly on cluster models. Initially the growth rate is higher for cluster models with higher central density. However, this high growth rate drops in a few crossing times of the cluster, and after several crossing times, roughly half of the total mass is in the common halo in all models we considered. In the central region of clusters density cusps with the profile \rho ~ r^{-1.2} develops regardless of the models of clusters and galaxies. We also found that the galaxies evolved so as to satisfy the relation between the masses of galaxies m_{gx} and their velocity dispersion \sigma_{gx} expressed as m_{gx} \propto \sigma_{gx}^{3~4} for all galaxy models as a consequence of their dynamical evolution through galaxy-galaxy interactions. We discuss the relation between our result and the observed Faber-Jackson relation.

In this letter we describe the pseudoparticle multipole method (P2M2), a new method to express multipole expansion by a distribution of pseudoparticles. We can use this distribution of particles to calculate high order terms in both the Barnes-Hut treecode and FMM. The primary advantage of P2M2 is that it works on GRAPE. GRAPE is a special-purpose hardware for the calculation of gravitational force between particles. Although the treecode has been implemented on GRAPE, we could handle terms only up to dipole, since GRAPE can calculate forces from point-mass particles only. Thus the calculation cost grows quickly when high accuracy is required. With P2M2, the multipole expansion is expressed by particles, and thus GRAPE can calculate high order terms. Using P2M2, we implemented an arbitrary-order treecode on GRAPE-4. Timing result shows GRAPE-4 accelerates the calculation by a factor between 10 (for low accuracy) to 150 (for high accuracy). Even on general-purpose programmable computers, our method offers the advantage that the mathematical formulae and therefore the actual program is much simpler than that of the direct implementation of multipole expansion.

We study the evolution and observability of young compact star clusters within about 200pc of the Galactic center. Calculations are performed using direct N-body integration on the GRAPE-4, including the effects of both stellar and binary evolution and the external influence of the Galaxy. The results of these detailed calculations are used to calibrate a simplified model applicable over a wider range of cluster initial conditions. We find that clusters within 200 pc from the Galactic center dissolve within about 70 Myr. However, their projected densities drop below the background density in the direction of the Galactic center within 20 Myr, effectively making these clusters undetectable after that time. Clusters farther from the Galactic center but at the same projected distance are more strongly affected by this selection effect, and may go undetected for their entire lifetimes. Based on these findings, we conclude that the region within 200 pc of the Galactic center could easily harbor some 50 clusters with properties similar to those of the Arches or the Quintuplet systems.

In this paper we describe the current status of the GRAPE-6 project to develop a special-purpose computer with a peak speed exceeding 100 Tflops for the simulation of astrophysical N-body problems. One of the main targets of the GRAPE-6 project is the simulation of dense stellar systems. In this paper, therefore, we overview the basic algorithms we use for the simulation of dense stellar systems and their characteristics. We then describe how we designed GRAPE hardwares to meet the requirements of these algorithms. GRAPE-6 will be completed by the year 2001. As an example of what science can be done on GRAPE-6, we describe our work on the galactic center with massive black holes performed on GRAPE-4, the predecessor of GRAPE-6.

The evolution of star clusters is studied using N-body simulations in which the evolution of single stars and binaries are taken self-consistently into account. Initial conditions are chosen to represent relatively young Galactic open clusters, such as the Pleiades, Praesepe and the Hyades. The calculations include a realistic mass function, primordial binaries and the external potential of the parent Galaxy. Our model clusters are generally significantly flattened in the Galactic tidal field, and dissolve before deep core collapse occurs. The binary fraction decreases initially due to the destruction of soft binaries, but increases later because lower mass single stars escape more easily than the more massive binaries. At late times, the cluster core is quite rich in giants and white dwarfs. There is no evidence for preferential evaporation of old white dwarfs, on the contrary the formed white dwarfs are likely to remain in the cluster. Stars tend to escape from the cluster through the first and second Lagrange points, in the direction of and away from the Galactic center. Mass segregation manifests itself in our models well within an initial relaxation time. As expected, giants and white dwarfs are much more strongly affected by mass segregation than main-sequence stars. Open clusters are dynamically rather inactive. However, the combined effect of stellar mass loss and evaporation of stars from the cluster potential drives its dissolution on a much shorter timescale than if these effects are neglected. The often-used argument that a star cluster is barely older than its relaxation time and therefore cannot be dynamically evolved is clearly in error for the majority of star clusters.

As an entry for the 1999 Gordon Bell price/performance prize, we report an astrophysical N-body simulation performed with a treecode on GRAPE-5 (Gravity Pipe 5) system, a special-purpose computer for astrophysical N-body simulations. The GRAPE-5 system has 32 pipeline processors specialized for the gravitational force calculation. Other operations, such as tree construction, tree traverse and time integration, are performed on a general purpose workstation. The total cost for the GRAPE-5 system is 40,900 dollars. We performed a cosmological N-body simulation with 2.1 million particles, which sustained a performance of 5.92 Gflops averaged over 8.37 hours. The price per performance obtained is 7.0 dollars per Mflops.

We investigated the evolution of clusters of galaxies using self-consistent $N$-body simulations in which each galaxy was modeled by many particles. We carried out simulations for about 20 cases using different initial conditions. In all simulations, clusters were initially in virial equilibrium. We found that more than half of the total mass escaped from individual galaxies within a few crossing times of the cluster, and that a diffuse halo was formed. The growth rate of the common halo depended on the size of individual galaxies only weakly. The stripping of the mass from galaxies was mainly due to the interaction of galaxies, not due to the effect of the tidal field of the cluster potential. The amount of stripped mass was larger for galaxies in the central region than for those in the outer region, since the interactions were more frequent in the central region. As a result, a positive correlation between the distance from the center and the mass of the galaxy developed. The volume-density profile of the common halo is expressed as $\rho\propto r^{-1}$ in the central region. This mass distribution is consistent with the mass distribution in clusters estimated using X-ray observations.

In this letter, we show that a massive black hole (MBH) which falls into the center of a galaxy in dynamical timescale leaves a weak cusp ($\rho \propto r^{-1/2}$) around it, which is in good agreement with the recent observations of luminous ellipticals by Hubble Space Telescope. Such event is a natural outcome of merging of two galaxies which have central MBHs. This is the only known mechanism to form weak cusps in luminous ellipticals. Therefore, the existence of the weak cusps indicates the central BHs of luminous ellipticals have fallen to the center from outside, most likely during a major merger event.

We have developed a special-purpose computer for gravitational many-body simulations, GRAPE-5. GRAPE-5 is the successor of GRAPE-3. Both consist of eight custom pipeline chips (G5 chip and GRAPE chip). The difference between GRAPE-5 and GRAPE-3 are: (1) The G5 chip contains two pipelines operating at 80 MHz, while the GRAPE chip had one at 20 MHz. Thus, the calculation speed of the G5 chip and that of GRAPE-5 board are 8 times faster than that of GRAPE chip and GRAPE-3 board. (2) The GRAPE-5 board adopted PCI bus as the interface to the host computer instead of VME of GRAPE-3, resulting in the communication speed one order of magnitude faster. (3) In addition to the pure 1/r potential, the G5 chip can calculate forces with arbitrary cutoff functions, so that it can be applied to Ewald or P^3M methods. (4) The pairwise force calculated on GRAPE-5 is about 10 times more accurate than that on GRAPE-3. On one GRAPE-5 board, one timestep of 128k-body simulation with direct summation algorithm takes 14 seconds. With Barnes-Hut tree algorithm (theta = 0.75), one timestep of 10^6-body simulation can be done in 16 seconds.

We have developed PROGRAPE-1 (PROgrammable GRAPE-1), a programmable multi-purpose computer for many-body simulations. The main difference between PROGRAPE-1 and "traditional" GRAPE systems is that the former uses FPGA (Field Programmable Gate Array) chips as the processing elements, while the latter rely on the hardwired pipeline processor specialized to gravitational interactions. Since the logic implemented in FPGA chips can be reconfigured, we can use PROGRAPE-1 to calculate not only gravitational interactions but also other forms of interactions such as van der Waals force, hydrodynamical interactions in SPH calculation and so on. PROGRAPE-1 comprises two Altera EPF10K100 FPGA chips, each of which contains nominally 100,000 gates. To evaluate the programmability and performance of PROGRAPE-1, we implemented a pipeline for gravitational interaction similar to that of GRAPE-3. One pipeline fitted into a single FPGA chip, which operated at 16 MHz clock. Thus, for gravitational interaction, PROGRAPE-1 provided the speed of 0.96 Gflops-equivalent. PROGRAPE will prove to be useful for wide-range of particle-based simulations in which the calculation cost of interactions other than gravity is high, such as the evaluation of SPH interactions.

In this paper, we study the thermal relaxation in the one-dimensional self-gravitating system, or the so-called sheet model. According to the standard argument, the thermal relaxation time of the system is around $Nt_c$, where $N$ is the number of sheets and $t_c$ is the crossing time. It has been claimed that the system does not reach the thermal equilibrium in this thermal relaxation timescale, and that it takes much longer time for the system to reach true thermal equilibrium. We demonstrate that this behavior is explained simply by the fact that the relaxation time is long. The relaxation time of sheets with average binding energy is $\sim 20 Nt_c$, and that of sheets with high energy can exceed $1000 Nt_c$. Thus, one needs to take the average over the relaxation timescale of high-energy sheets, if one wants to look at the thermal characteristic of these high energy sheets.

The evolution of young compact star clusters is studied using N-body simulations in which both stellar evolution and physical collisions between stars are taken into account. The initial conditions are chosen to represent R136, a compact star cluster in the 30 Doradus region of the Large Magellanic Cloud. The present runs do not include the effects of primordial binaries. We find that physical collisions between stars in these models are frequent, and that the evolution of the most massive stars and the dynamical evolution of the cluster are closely coupled. In all cases, a single star grows steadily in mass through mergers with other stars, forming a very massive (>100 Msun) star in less than 3-4 Myr. The growth rate of this runaway merger is much larger than estimates based on simple cross-section arguments, mainly because the star is typically found in the core and tends to form binaries with other massive stars there. The runaway is ``rejuvenated'' by each new collision, and its lifetime is extended considerably as a consequence. Observationally, such a star will appear in the Hertzsprung-Russell diagram as a blue straggler. When the runaway forms a black hole, the binary in which it is found is usually dissociated.

We present the pseudo-particle multipole method (P2M2), a new method to handle multipole expansion in fast multipole method and treecode. This method uses a small number of pseudo-particles to express multipole expansion. With this method, the implementation of FMM and treecode with high-order multipole terms is greatly simplified. We applied P2M2 to treecode and combined it with special-purpose computer GRAPE. Extensive tests on the accuracy and calculation cost demonstrate that the new method is quite attractive.

The next decade will be an exciting time for computational physicists. After 50 years of being forced to use standardized commercial equipment, it will finally become relatively straightforward to adapt one's computing tools to one's own needs. The breakthrough that opens this new era is the now wide-spread availability of programmable chips that allow virtually every computational scientist to design his or her own special-purpose computer.

The GRAPE-4, the world's fastest computer in 1995-1997, has produced some major scientific results, through a wide diversity of large-scale simulations in astrophysics. Applications have ranged from planetary formation, through the evolution of star clusters and galactic nuclei, to the formation of galaxies and clusters of galaxies.

We investigated the dynamical reaction of the central region of galaxies to a falling massive black hole by N-body simulations. As the initial galaxy model, we used an isothermal King model and placed a massive black hole at around the half-mass radius of the galaxy. We found that the central core of the galaxy is destroyed by the heating due to the black hole and that a very weak density cusp ($\rho \propto r^{-\alpha}$, with $\alpha \sim 0.5$) is formed around the black hole. This result is consistent with recent observations of large elliptical galaxies with Hubble Space Telescope. The velocity of the stars becomes tangentially anisotropic in the inner region, while in the outer region the stars have radially anisotropic velocity dispersion. The radius of the weak cusp region is larger for larger black hole mass. Our result naturally explains the formation of the weak cusp found in the previous simulations of galaxy merging, and implies that the weak cusp observed in large elliptical galaxies may be formed by the heating process by sinking black holes during merging events.

In this paper we describe a new approach to implement the O(N) fast multipole method and $O(N\log N)$ tree method, which uses pseudoparticles to express the potential field. The new method is similar to Anderson's method, which uses the values of potential at discrete points to represent the potential field. However, for the same expansion order the new method is more accurate and computationally efficient.

Using direct N-body simulations which include both the evolution of single stars and the tidal field of the parent galaxy, we study the dynamical evolution of globular clusters and rich open clusters. We compare our results with other N-body simulations and Fokker-Planck calculations. Our simulations, performed on the GRAPE-4, employ up to 32,768 stars. The results are not in agreement with Fokker-Planck models, in the sense that the lifetimes of stellar systems derived using the latter are an order of magnitude smaller than those obtained in our simulations. For our standard run, Fokker-Plank calculations obtained a lifetime of 0.28 Gyr, while our equivalent N-body calculations find about 4 Gyr. The principal reason for the discrepancy is that a basic assumption of the Fokker-Plank approach is not valid for typical cluster parameters. The stellar evolution timescale is comparable to the dynamical timescale, and therefore the assumption of dynamical equilibrium leads to an overestimate of the dynamical effects of mass loss. Our results suggest that the region in parameter space for which Fokker-Planck studies of globular cluster evolution, including the effects of both stellar evolution and the galactic tidal field, are valid is limited. The discrepancy is largest for clusters with short lifetimes.

We present the results from a two-day study in which we discussed various implementations of Smooth Particle Hydrodynamics (SPH), one of the leading methods used across a variety of areas of large-scale astrophysical simulations. In particular, we evaluated the suitability of designing special hardware extensions, to further boost the performance of the high-end general purpose computers currently used for those simulations. We considered a range of hybrid architectures, consisting of a mix of custom LSI and reconfigurable logic, combining the extremely high throughput of Special-Purpose Devices (SPDs) with the flexibility of reconfigurable structures, based on Field Programmable Gate Arrays (FPGAs). The main findings of our workshop consist of a clarification of the decomposition of the computational requirements, together with specific estimates for cost/performance improvements that can be obtained at each stage in this decomposition, by using enabling hardware technology to accelerate the performance of general purpose computers.

The performance and accuracy of a GRAPE-3 system for collisionless N-body simulations is discussed. After a description of the hardware configurations available to us at Marseille, and the usefulness of on-line analysis, we concentrate on the performance and accuracy of the direct summation and tree code software. For the former we discuss the sources of round-off errors. The standard Barnes-Hut tree code has been modified for GRAPE-3, by dividing particles into blocks and performing the tree traversal for the whole block. The forces are then calculated by direct summation over the whole interaction list. The performance of the tree code depends on the number of particles in the block, the optimum number depending on the speed of the front end and the number of boards. The time per step scales with the number of particles, and decreases only weakly as the tolerance increases. We discuss the effect of the front end and compare performances with those of general purpose computers. The accuracy of both direct summation and the tree code is discussed as function of the number of particles and the softening. The increased role of direct summation in the force calculation enhances the accuracy compared to the standard tree code. Finally, we follow the evolution of an isolated barred galaxy using different hardware and software in order to assess the reliability of our results. We find excellent agreement between the pattern speed of the bar in direct summation simulations run on the high precision GRAPE-4 machines and on our GRAPE-3 system. The agreement with the tree code is very good for tolerance values smaller than about 1.0. We conclude that GRAPE-3 systems are well suited for collisionless simulations, in particular those of galaxies.

We developed a PCI interface for GRAPE systems. GRAPE(GRAvity piPE) is a special-purpose computer for gravitational N-body simulations. A GRAPE system consists of GRAPE processor boards and a host computer. GRAPE processors perform the calculation of gravitational forces between particles. The host computer performs the rest of calculations. The newest of GRAPE machines, the GRAPE-4, achieved the peak performance of 1.08 Tflops. The GRAPE-4 system uses TURBOChannel for the interface to the host, which limits the selection of the host computer. The TURBOChannel bus is not supported by any of recent workstations. We developed a new host interface board which adopts the PCI bus instead of the TURBOChannel. PCI is an I/O bus standard developed by Intel. It has fairly high peak transfer speed, and is available on wide range of computers, from PCs to supercomputers. Thus, the new interface allows us to connect GRAPE-4 to a wide variety of host computers. In test runs with a Barnes-Hut treecode, we found that the performance of new system with PCI interface is 40% better than that of the original system.

We present two types of meta-algorithm that can greatly improve the accuracy of existing algorithms for integrating the equations of motion of dynamical systems. The first meta-algorithm takes an integrator that is time-symmetric only for constant time steps, and ensures time-symmetry even in the case of varying time steps. The second meta-algorithm can be applied to any self-starting integration scheme to create time-symmetry, both for constant and for variable time steps, even if the original scheme was not time-symmetric. Our meta-algorithms are most effective for Hamilton systems or systems with periodic solutions. If the system is not Hamiltonian (for example, if some dissipative force exists), our methods are still useful so long as the dissipation is small.

We present the outline of a research project aimed at designing and constructing a hybrid computing system that can be easily scaled up to petaflops speeds. As a first step, we envision building a prototype which will consist of three main components: a general-purpose, programmable front end, a special-purpose, fully hardwired computing engine, and a multi-purpose, reconfigurable system. The driving application will be a suite of particle-based large-scale simulations in various areas of physics. The prototype system will achieve performance in the $\sim 50 - 100$ teraflops range for a broad class of applications in this area. The combination of a hardwired petaflops-class computational engine and a front end with sustained speed on the order of 10 gigaflops can produce extremely high performance, but only for the limited class of problems in which there exists a single bottleneck with computing cost dominating the total. While the calculation for which the Grape-4 (our system's immediate predecessor) was designed is a prime example of such a problem, in many other applications the primary computational bottleneck, while still related to an inverse-square (gravitational, Coulomb, etc.) force, requires less than 99% of the computing power. Although the remainder of the CPU time is typically dominated by just one secondary bottleneck, its nature varies greatly from problem to problem. It is not cost-effective to attempt to design custom chips for each new problem that arises. FPGA-based systems can restore the balance, guaranteeing scalability from the teraflops to the petaflops domain, while still retaining significant flexibility. (abbreviated abstract)

We present cross sections and reaction rates for merging to occur during encounters of equal-mass spherical galaxies. As an application, we determine the rate of galaxy merging in clusters of galaxies. We present results for two types of Plummer models (a full and a truncated one), two King models and the Hernquist model. Cross sections are determined on the basis of a large number ($\sim 500$) of simulations of galaxy encounters, using the 10-Gigaflops GRAPE 3A special-purpose computer. We characterize the overall merger rate of galaxies in a galaxy cluster by a single number, derived from our cross sections by an integration over galaxy encounter velocities in the limit of a constant density in velocity space. For small clusters, where the cluster velocity dispersion may not significantly exceed the internal velocity dispersion of the individual galaxies, this constant-density approximation may not be valid. For those cases, we present separate results, based on integrations of our cross sections over Maxwellian velocity distributions. Finally, tidal effects from the cluster potential as well as from neighboring galaxies may prevent a barely bound galaxy pair from spiraling in after their first encounter. We give a quantitative estimate of the resulting reduction in the actual merger rate, due to these tidal interactions.

We discuss the merging rates in compact groups of 5 identical elliptical galaxies. All groups have the same mass and binding energy. We consider both cases with individual halos and cases where the halo is common to all galaxies and enveloping the whole group. In the latter situation the merging rate is slower if the halo is more massive. The mass of individual halos has little influence on the merging rates, due to the fact that all galaxies in our simulations have the same mass, and so the more extended ones have a smaller velocity dispersion. Groups with individual halos merge faster than groups with common halos if the configuration is centrally concentrated, like a King distribution of index 10. On the other hand for less concentrated configurations the merging is initially faster for individual halo cases, and slower after part of the group has merged. In cases with common halo, centrally concentrated configurations merge faster for high halo-to-total mass ratios and slower for low halo-to-total mass ratios. Groups whose virial ratio is initially less than one merge faster, while groups that have initially cylindrical rotation merge slower than groups starting in virial equilibrium. In order to test how long a virialised group can survive before merging we followed the evolution of a group with a high halo-to-total mass ratio and a density distribution with very little central concentration. We find that the first merging occurred only after a large number of crossing times, which with areasonable calibration should be larger than a Hubble time. Hence, at least for appropriate initial conditions, the longevity of compact groups is not necessarily a problem, which is an alternative explanation to why we observe so many compact groups despite the fact that their lifetimes seem short.

Observed cusps with density profiles $\rho\propto r^{-1}$ or shallower, in the central regions of galaxies, cannot be reproduced in the standard Cold Dark Matter (CDM) picture of hierarchical clustering. Previous claims to the contrary were based on simulations with relatively few particles, and substantial softening. We present simulations with particle numbers an order of magnitude higher, and essentially no softening, and show that typical central density profiles are clearly steeper than $\rho\propto r^{-1}$. The observed shallower profiles may have formed through the smoothing effect of the spiral-in of central black holes in previous merger phases. In addition, we confirm the presence of a temperature inversion in the inner 5 kpc of massive galactic halos, and illustrate its formation as a natural result of the merging of unequal progenitors.

We report on resent N-body simulations of galaxy formation performed on the GRAPE-4 (GRAvity PipE) system, a special-purpose computer for astrophysical N-body simulations. We review the astrophysical motivation, the algorithm, the actual performance, and the price per performance. The performance obtained is 332 Gflops averaged over 185 hours for a simulation of a galaxy formation with 786,400 particles. The price per performance obtained is 4,600 dollars per Gflops. The configuration used for the simulation consists of 1,269 pipeline processors and has a peak speed of 663 Gflops.

We investigated the evolution of the black hole binary formed by the merging of two galaxies each containing a central massive black hole. Our main goal here is to determine if the black hole binary can merge through the hardening by dynamical friction and the gravitational wave radiation. We performed $N$-body simulations of merging of two galaxies with wide range of total number of particles to investigate the effect of the number of particles on the evolution of the black hole binary. We found that the evolution timescale was independent of the number of particles in the galaxy $N$ until the separation reaches a critical value. After the separation became smaller than this critical value, the evolution timescale was longer for larger number of particles. Qualitatively, this behavior is understood naturally the result the ``loss-cone'' effect. However, the dependence of the timescale on $N$ is noticeably weaker than the theoretical prediction. In addition, the critical separation is smaller than the theoretical prediction. The timescale of evolution through gravitational radiation at this critical separation is longer than the Hubble time. We discuss the reason of these discrepancy and the implication of the present result on the structure of the ellipticals and QSO activities.

A number of globular clusters appear to have undergone core collapse, in the sense that their predicted collapse time is much shorter than their current age. Simulations using gas models and Fokker-Planck approximation have shown that the central density of a globular cluster after the collapse undergoes nonlinear oscillation with large amplitude (gravothermal oscillation). However, whether such an oscillation actually takes place in a real $N$-body system has remained unsolved, because an $N$-body simulation with a sufficiently high resolution would have required the computing resource of the order of several Gflops$\cdot$years. In the present paper, we report the result of such a simulation, performed on a dedicated special-purpose computer GRAPE-4. We simulated the evolution of isolated point-mass systems with up to 32,768 particles. The largest number of particles reported previously is 10,000. We confirmed that gravothermal oscillation takes place in an $N$-body system. The expansion phase shows all signatures that are considered as the evidences of the gravothermal nature of the oscillation. At the maximum expansion, the core radius is $\sim 1$\% of the half-mass radius for the run with 32,768 particles. The maximum core size $r_c$ depends on $N$, as $<r_c> \propto N^{-1/3}$.

In this paper we describe a new algorithm for the long-term numerical integration of the two-body problem, in which two particles interact under a Newtonian gravitational potential. Although analytical solutions exist in the unperturbed and weakly perturbed cases, numerical integration is necessary in situations where the perturbation is relatively strong. Kustaanheimo--Stiefel (KS) regularization is widely used to remove the singularity in the equations of motion, making it possible to integrate orbits having very high eccentricity. However, even with KS regularization, long-term integration is difficult, simply because the required accuracy is usually very high. We present a new time-integration algorithm which has no secular error in either the binding energy or the eccentricity, while allowing variable stepsize. The basic approach is to take a time-symmetric algorithm, then apply an implicit criterion for the stepsize to ensure strict time reversibility. We describe the algorithm in detail and present the results of numerical tests involving long-term integration of binaries and hierarchical triples. In all cases studied, we found no systematic error in either the energy or the angular momentum. We also found that its calculation cost does not become higher than those of existing algorithms. By contrast, the stabilization technique, which has been widely used in the field of collisional stellar dynamics, conserves energy very well but does not conserve angular momentum.

The homogeneity of the cosmic microwave background radiation (CBR) is one of the most severe constraint for theories of the structure formation in the universe. We investigated the effect of the gravitational scattering (lensing) of galaxies, clusters of galaxies, and superclusters on the anisotropy of the CBR by numerical simulations. Although this effect was thought to be unimportant, we found that the gravitational scatterings by superclusters can significantly reduce the anisotropy of the CBR. We took into account the exponential growth of the distance between two rays due to multiple scatterings. The bending angle of each ray grows through the random walk process. On the other hand, difference between two rays grows exponentially while it is small. This exponential growth is caused by coherent scatterings that two rays suffer, and was neglected in the previous studies. The gravitational scattering by superclusters reduces the observed temperature anisotropy of the CBR at present time approximately by 40--60 $\%$ from that at the recombination time for angular scale up to a few degrees if the supercluster were formed at $z=$2--4.

We give an estimate of the relative error in the angular measurement of observations for high redshift objects induced by gravitational scatterings (lensing). Gunn (1967) concluded that the gravitational scatterings by galaxies induce the relative error of a few percent in the observations for objects at $z=1$. This estimate has been considered as a fundamental limitation of accuracy of the angular measurements in the observational cosmology. In multiple graviational scatterings, bending angle of single ray grows through the random work process. Gunn (1967) assumed that the difference of nearby rays also grows through the random walk process. However, distance between nearby photons grows exponentially because the two rays suffer coherent scatterings by the same scattering object. This exponential growth continues as long as the scattering is coherent. In the case of scattering by individual galaxies, the exponential growth continues until the angular distance reaches an arcminute or so. The relative error of the angular measurements under an arcminute due to the exponential growth is $\sim 30\%$ at $z=1$ and exceeds $100\%$ at $z=3$, in the case that the density parameter of galaxies is 0.2. The effects of clusters of galaxies or superclusters are more difficult to estimate accurately, but might be significant. In the case of supercluster the angular measurements up to a few degrees could be affected.