  Spherical harmonic modes of gravitational waveforms for inspiraling compact
binaries in eccentric orbits from post-Newtonian (PN) theory accurate to third
post-Newtonian order, and those extracted from numerical relativity (NR)
simulations for binary black holes (BBHs) are compared. We combine results from
the two approaches (PN and NR) to construct time-domain hybrid waveforms that
describe the complete evolution of BBH mergers through inspiral-merger-ringdown
(IMR) stages. These hybrids are then used in constructing a fully analytical
dominant mode ($\ell$=2, $|m|$=2) eccentric IMR model. A simple extension to a
multi-mode model based on this dominant mode model is also presented. Overlaps
with quasi-circular IMR waveform models including the effect of higher modes,
maximized over a time- and phase-shift, hint at the importance (mismatches
$>1\%$) of including eccentricity in gravitational waveforms when analysing
BBHs lighter than $\sim 80 M_{\odot}$, irrespective of the binary's
eccentricity (as it enters the LIGO bands), or mass-ratio. Combined impact of
eccentricity and higher modes seems to become more apparent through smaller
overlaps with increasing inclination angles and mass ratios. Additionally, we
show that the state-of-the-art quasi-circular models including the effect of
higher modes will not be adequate in extracting source properties for signals
with initial eccentricities $e_0$ $\gtrsim0.1$.
