Abstract
The recent development of new gravitational wave models with higher harmonic modes now allows us to perform, for the first time, fully Bayesian inference studies that can include the effects of these subdominant modes. With the use of NRHybSur3dq8, a recently-developed aligned-spin surrogate model built using numerical relativity data, along with the highly paralellizable rapid inference algorithm known as RIFT, we demonstrate the importance of higher modes on the parameter inference of coalescing massive binary black holes. We consider cases which are relevant to the current three-detector network of observatories, as well as probable candidates for future observing runs. We show that the exclusion of higher modes produces significant parameter biases for asymmetric binaries with mass ratios q > 1,and can even influence the posterior probability distributions for cases with comparable mass binaries and at low signal amplitude. Lastly, we discuss the impact of our results on individual spin measurability, source population inference, and self-consistency tests of general relativity.