  We present a frequency domain reduced order model (ROM) for the aligned-spin
effective-one-body (EOB) model for binary black holes (BBHs) SEOBNRv4HM that
includes the spherical harmonics modes $(\ell, |m|) = (2,1),(3,3),(4,4),(5,5)$
beyond the dominant $(\ell, |m|) = (2,2)$ mode. These higher modes are crucial
to accurately represent the waveform emitted from asymmetric BBHs. We discuss a
decomposition of the waveform, extending other methods in the literature, that
allows us to accurately and efficiently capture the morphology of higher mode
waveforms. We show that the ROM is very accurate with median (maximum) values
of the unfaithfulness against SEOBNRv4HM lower than $0.001\% (0.03\%)$ for
total masses in $[2.8,100] M_\odot$. For a total mass of $M = 300 M_\odot$ the
median (maximum) value of the unfaithfulness increases up to $0.004\%
(0.17\%)$. This is still at least an order of magnitude lower than the
estimated accuracy of SEOBNRv4HM compared to numerical relativity simulations.
The ROM is two orders of magnitude faster in generating a waveform compared to
SEOBNRv4HM. Data analysis applications typically require
$\mathcal{O}(10^6-10^8)$ waveform evaluations for which SEOBNRv4HM is in
general too slow. The ROM is therefore crucial to allow the SEOBNRv4HM waveform
to be used in searches and Bayesian parameter inference. We present a targeted
parameter estimation study that shows the improvements in measuring binary
parameters when using waveforms that includes higher modes and compare against
three other waveform models.
