  For a brief moment, a binary black hole (BBH) merger can be the most powerful
astrophysical event in the visible universe. Here we present a model fit for
this gravitational-wave peak luminosity of nonprecessing quasicircular BBH
systems as a function of the masses and spins of the component black holes,
based on numerical relativity (NR) simulations and the hierarchical fitting
approach introduced by X. Jim\'enez-Forteza et al. [Phys. Rev. D 95, 064024
(2017), arXiv:1611.00332]. This fit improves over previous results in accuracy
and parameter-space coverage and can be used to infer posterior distributions
for the peak luminosity of future astrophysical signals like GW150914 and
GW151226. The model is calibrated to the l<=6 modes of 378 nonprecessing NR
simulations up to mass ratios of 18 and dimensionless spin magnitudes up to
0.995, and includes unequal-spin effects. We also constrain the fit to
perturbative numerical results for large mass ratios. Studies of key
contributions to the uncertainty in NR peak luminosities, such as (i) mode
selection, (ii) finite resolution, (iii) finite extraction radius, and (iv)
different methods for converting NR waveforms to luminosity, allow us to use NR
simulations from four different codes as a homogeneous calibration set. This
study of systematic fits to combined NR and large-mass-ratio data, including
higher modes, also paves the way for improved inspiral-merger-ringdown waveform
models.
