  Coalescing binaries of neutron stars (NS) and black holes (BH) are one of the
most important sources of gravitational waves for the upcoming network of
ground based detectors. Detection and extraction of astrophysical information
from gravitational-wave signals requires accurate waveform models. The
Effective-One-Body and other phenomenological models interpolate between
analytic results and $10-30$ orbit numerical relativity (NR) merger
simulations. In this paper we study the accuracy of these models using new NR
simulations that span $36-88$ orbits, with mass-ratios and black hole spins
$(q,\chi_{BH}) = (7, \pm 0.4), (7, \pm 0.6)$, and $(5, -0.9)$. We find that:
(i) the recently published SEOBNRv1 and SEOBNRv2 models of the
Effective-One-Body family disagree with each other (mismatches of a few
percent) for black hole spins $\geq 0.5$ or $\leq -0.3$, with waveform mismatch
accumulating during early inspiral; (ii) comparison with numerical waveforms
indicate that this disagreement is due to phasing errors of SEOBNRv1, with
SEOBNRv2 in good agreement with all of our simulations; (iii) Phenomenological
waveforms disagree with SEOBNRv2 over most of the NSBH binary parameter space;
(iv) comparison with NR waveforms shows that most of the model's dephasing
accumulates near the frequency interval where it switches to a phenomenological
phasing prescription; and finally (v) both SEOBNR and post-Newtonian (PN)
models are effectual for NSBH systems, but PN waveforms will give a significant
bias in parameter recovery. Our results suggest that future gravitational-wave
detection searches and parameter estimation efforts targeted at NSBH systems
with $q\lesssim 7$ and $\chi_\mathrm{BH} \approx [-0.9, +0.6]$ will benefit
from using SEOBNRv2 templates. For larger black hole spins and/or binary
mass-ratios, we recommend the models be further investigated as suitable NR
simulations become available.
