  The recent discovery of GW150914, the binary black hole merger detected by
Advanced LIGO, has the potential to revolutionize observational astrophysics.
But to fully utilize this new window into the universe, we must compare these
new observations to detailed models of binary black hole formation throughout
cosmic time. Expanding upon our previous work (Rodriguez et al., 2015), we
study merging binary black holes formed in globular clusters using our Monte
Carlo approach to stellar dynamics. We have created a new set of 52 cluster
models with different masses, metallicities, and radii to fully characterize
the binary black hole merger rate. These models include all the relevant
dynamical processes (such as two-body relaxation, strong encounters, and
three-body binary formation) and agree well with detailed direct N-body
simulations. In addition, we have enhanced our stellar evolution algorithms
with updated metallicity-dependent stellar wind and supernova prescriptions,
allowing us to compare our results directly to the most recent population
synthesis predictions for merger rates from isolated binary evolution. We
explore the relationship between a cluster's global properties and the
population of binary black holes that it produces. In particular, we derive a
numerically calibrated relationship between the merger times of ejected black
hole binaries and a cluster's mass and radius. We explore the masses and mass
ratios of these binaries as a function of redshift, and find a merger rate of
~5 Gpc$^{-3}$ yr$^{-1}$ in the local universe, with 80% of sources having total
masses from $32M_{\odot}$ to $64M_{\odot}$. Under standard assumptions,
approximately 1 out of every 7 binary black hole mergers in the local universe
will have originated in a globular cluster, but we also explore the sensitivity
of this result to different assumptions for binary stellar evolution.
[Abridged]
