Total: 1
Numerical-relativity simulations offer a unique approach to investigating the dynamics of binary neutron star mergers and provide the most accurate predictions of waveforms in the late inspiral phase. However, the numerical predictions are prone to systematic biases originating from the construction of initial quasi-circular binary configurations, the numerical methods used to evolve them, and to extract gravitational signals. To assess uncertainties arising from these aspects, we analyze mergers of highly spinning neutron stars with dimensionless spin parameter χ=0.5. The initial data are prepared by two solvers, \textsc{FUKA} and \textsc{SGRID}, which are then evolved by two independent codes, \textsc{SACRA} and \textsc{BAM}. We assess the impact of numerical discretizations, finite extraction radii, and differences in numerical frameworks on the resulting gravitational waveforms. Our analysis reveals that the primary source of uncertainty in numerical waveforms is the evolution code, while the initial data solver has a smaller impact. We also compare our numerical-relativity waveforms with state-of-the-art analytical models, finding that the discrepancies between them exceed the estimated numerical uncertainties. Few suggestions are offered: (i) the analytic waveform becomes an inadequate approximation after the two neutron stars come into contact and the binary enters the essentially-one-body phase, (ii) the analytical models may not capture finite-size effects beyond quadrupole moment, and (iii) the inconsistent use of the binary black hole baseline in the analytical models may also be contributing to these discrepancies. The presented results benchmark the error budget for numerical waveforms of binary neutron star mergers, and provide information for the analytic models to explore further the high spin parameter space of binary neutron star mergers.