Formation and evolution of binary black holes in $N$-body simulations of star clusters with up to two million stars

Understanding binary black hole (BBH) dynamics in dense star clusters is key to interpreting the gravitational wave detections by LIGO and Virgo. Here, we perform $N$-body simulations of star clusters, focusing on BBH formation mechanisms, dynamical evolution and merging properties. We explore a wide parameter space of initial conditions, with cluster masses ranging from $10^{4}$ to $10^{6}~\mathrm{M_{\odot}}$, densities from $10^{3}$ to $10^{5}~\rm M_{\odot}pc^{-3}$, and up to $100\%$ of massive stars in binaries. We show that most BBH mergers originate from the primordial binary population rather than being dynamically assembled, and that the evolution towards merger for most of these binaries is not significantly altered by dynamical encounters. As a result, the overall number of BBH mergers from the $N$-body simulations is nearly identical to that obtained when the same stellar population is evolved in isolation. Contrary to theoretical expectations, nearly all dynamically formed BBH mergers occur when the binary is still bound to its host cluster, with $\simeq 90\%$ of all dynamical mergers occurring within the cluster core region. In about half of these mergers the binary is part of a stable black hole-triple system. In one model, stellar mergers lead to the formation of a $\simeq 200\,\mathrm{M_\odot}$ black hole, which then grows to $\simeq 300\,\mathrm{M_\odot}$ through black hole mergers. Our study highlights the importance of detailed $N$-body simulations in capturing the evolution of black hole populations in dense clusters and challenges conclusions based on semi-analytical and Monte Carlo methods.