  With an increasing number of expected gravitational-wave detections of binary
neutron star mergers, it is essential that gravitational-wave models employed
for the analysis of observational data are able to describe generic compact
binary systems. This includes systems in which the individual neutron stars are
millisecond pulsars for which spin effects become essential. In this work, we
perform numerical-relativity simulations of binary neutron stars with aligned
and anti-aligned spins within a range of dimensionless spins of $\chi \sim
[-0.28,0.58]$. The simulations are performed with multiple resolutions, show a
clear convergence order and, consequently, can be used to test existing
waveform approximants. We find that for very high spins gravitational-wave
models that have been employed for the interpretation of GW170817 and GW190425
are not capable of describing our numerical-relativity dataset. We verify
through a full parameter estimation study in which clear biases in the estimate
of the tidal deformability and effective spin are present. We hope that in
preparation of the next gravitational-wave observing run of the Advanced LIGO
and Advanced Virgo detectors our new set of numerical-relativity data can be
used to support future developments of new gravitational-wave models.
