  Gravitational wave observations of large mass ratio compact binary mergers
like GW190814 highlight the need for reliable, high-accuracy waveform templates
for such systems. We present NRHybSur2dq15, a new surrogate model trained on
hybridized numerical relativity (NR) waveforms with mass ratios $q\leq15$, and
aligned spins $|\chi_{1z}|\leq0.5$ and $\chi_{2z}=0$. We target the parameter
space of GW190814-like events as large mass ratio NR simulations are very
expensive. The model includes the (2,2), (2,1), (3,3), (4,4), and (5,5)
spin-weighted spherical harmonic modes, and spans the entire LIGO bandwidth
(with $f_{\mathrm{low}}=20$ Hz) for total masses $M \gtrsim 9.5 \, M_{\odot}$.
NRHybSur2dq15 accurately reproduces the hybrid waveforms, with mismatches below
$\sim 2 \times 10^{-3}$ for total masses $10 \, M_{\odot} \leq M \leq 300 \,
M_{\odot}$. This is at least an order of magnitude improvement over existing
semi-analytical models for GW190814-like systems. Finally, we reanalyze
GW190814 with the new model and obtain source parameter constraints consistent
with previous work.
