A chain regularization method is combined with special purpose computer hardware to study the evolution of massive black hole binaries at the centers of galaxies. Preliminary results with up to N = 0.26 × 10 6 particles are presented. The decay rate of the binary is shown to decrease with increasing N , as expected on the basis of theoretical arguments. The eccentricity of the binary remains small.Coalescence of binary supermassive black holes is potentially the strongest source of gravitational waves in the universe [1]. The coalescence rate is limited by the efficiency with which massive binaries can interact with stars and gas in a galaxy and reach the relativistic regime at separations of ∼ 10 −3 pc. Exchange of energy between a binary black hole and stars should also leave observable traces in the stellar distribution, perhaps allowing us to infer something about the merger history of galaxies from their nuclear structure [2]. Henry Kandrup worked on this problem shortly before his death. In "Supermassive Black Hole Binaries as Galactic Blenders" [3], Kandrup et al. investigated the effects of a massive binary on the stellar orbits near the center of spherical and nearly spherical galaxies. They showed that the periodically-varying potential due to the binary, coupled with the fixed potential from the galaxy, was effective at inducing chaos in the stellar orbits, leading to diffusion in both energy and configuration space and to ejection of stars from the nucleus. This study was a complement to earlier studies based on scattering experiments [4,5] in which the potential of the galaxy was ignored.Another approach to the binary black hole problem is via direct N -body techniques [6][7][8]. This approach is computationally challenging because of the need to handle close interactions between the star-and black hole particles with high precision. In addition, large particle numbers are required to avoid the effects of spurious relaxation [9,10]. Here, we present preliminary results of Nbody integrations of the binary black hole problem, in which close interactions between the black holes and stars are handled via the Mikkola-Aarseth chainregularization algorithm [11,12]. Recently Aarseth [13] described an application 1