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.

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 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 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.

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 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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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 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.

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.

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.

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.