  The "gravitational million-body problem," to model the dynamical evolution of
a self-gravitating, collisional N-body system with ~10^6 particles over many
relaxation times, remains a major challenge in computational astrophysics.
Unfortunately, current techniques to model such systems suffer from severe
limitations. A direct N-body simulation with more than 10^5 particles can
require months or even years to complete, while an orbit-sampling Monte Carlo
approach cannot adequately model the dynamics in a dense cluster core,
particularly in the presence of many black holes. We have developed a new
technique combining the precision of a direct N-body integration with the speed
of a Monte Carlo approach. Our Rapid And Precisely Integrated Dynamics code,
the RAPID code, statistically models interactions between neighboring stars and
stellar binaries while integrating directly the orbits of stars or black holes
in the cluster core. This allows us to accurately simulate the dynamics of the
black holes in a realistic globular cluster environment without the burdensome
N^2 scaling of a full N-body integration. We compare RAPID models of idealized
globular clusters to identical models from the direct N-body and Monte Carlo
methods. Our tests show that RAPID can reproduce the half-mass radii, core
radii, black hole ejection rates, and binary properties of the direct N-body
models far more accurately than a standard Monte Carlo integration while
remaining significantly faster than a full N-body integration. With this
technique, it will be possible to create more realistic models of Milky Way
globular clusters with sufficient rapidity to explore the full parameter space
of dense stellar clusters.
