Over the past year, a handful of new gravitational wave models have been developed to include multiple harmonic modes thereby enabling for the first time fully Bayesian inference studies including higher modes to be performed. Using one recently developed numerical relativity surrogate model, NRHybSur3dq8, we investigate the importance of higher modes on parameter inference of coalescing massive binary black holes. We focus on examples relevant to the current three-detector network of observatories, with a detector-frame mass set to
120 M⊙ and with signal amplitude values that are consistent with plausible candidates for the next few observing runs. We show that for such systems the higher mode content will be important for interpreting coalescing binary black holes, reducing systematic bias, and computing properties of the remnant object. Even for comparable-mass binaries and at low signal amplitude, the omission of higher modes can influence posterior probability distributions. We discuss the impact of our results on source population inference and self-consistency tests of general relativity. Our work can be used to better understand asymmetric binary black hole merger events, such as GW190412. Higher modes are critical for such systems, and their omission usually produces substantial parameter biases
