  Gravitational waves from the collision of binary neutron stars provide a
unique opportunity to study the behaviour of supranuclear matter, the
fundamental properties of gravity, and the cosmic history of our Universe.
However, given the complexity of Einstein's Field Equations, theoretical models
that enable source-property inference suffer from systematic uncertainties due
to simplifying assumptions. We develop a hypermodel approach to compare and
measure the uncertainty gravitational-wave approximants. Using state-of-the-art
models, we apply this new technique to the binary neutron star observations
GW170817 and GW190425 and the sub-threshold candidate GW200311_103121. Our
analysis reveals subtle systematic differences between waveform models, and a
frequency-dependence study suggests that this is due to the treatment of the
tidal sector. This new technique provides a proving ground for model
development, and a means to identify waveform-systematics in future observing
runs where detector improvements will increase the number and clarity of binary
neutron star collisions we observe.
