  Black hole binaries formed dynamically in globular clusters are believed to
be one of the main sources of gravitational waves in the Universe. Here, we use
our new population synthesis code, cBHBd, to determine the redshift evolution
of the merger rate density and masses of black hole binaries formed in globular
clusters. We simulate $\sim 2$ million models to explore the parameter space
that is relevant to real clusters and over all mass scales. We show that when
uncertainties on the initial cluster mass function and density are properly
taken into account, they become the two dominant factors in setting the
theoretical error bars on merger rates. Other model parameters (e.g., natal
kicks, black hole masses, metallicity) have virtually no effect on the local
merger rate density, although they affect the masses of the merging black
holes. Modelling the merger rate density as a function of redshift as
$R(z)=R_0(1+z)^\kappa$ at $z<2$, and marginalizing over uncertainties, we find:
$R_0=7.2^{+21.5}_{-5.5}{\rm Gpc^{-3}yr^{-1}}$ and $\kappa=1.6^{+0.4}_{-0.6}$
($90\%$ credibility). The rate parameters for binaries that merge inside the
clusters are ${R}_{\rm 0,in}=1.6^{+1.9}_{-1.0}{\rm Gpc^{-3}yr^{-1}}$ and
$\kappa_{\rm in}=2.3^{+1.3}_{-1.0}$; $\sim 20\%$ of these form as the result of
a gravitational-wave capture, implying that eccentric mergers from globular
clusters contribute $\lesssim 0.4 \rm Gpc^{-3}yr^{-1}$ to the local rate. A
comparison to the merger rate reported by LIGO-Virgo shows that a scenario in
which most of the detected black hole mergers are formed in globular clusters
is consistent with current constraints, and requires initial cluster half-mass
densities $\gtrsim 10^4 M_\odot \rm pc^{-3}$. Such models also reproduce the
inferred primary black hole mass distribution for masses $13-30 M_\odot$, but
under-predict the data outside this range.
