  Gravitational-wave signals from binary neutron star coalescences carry
information about the star's equation of state in their tidal signatures. A
major issue in the inference of the tidal parameters (or directly of the
equation of state) is the systematic error introduced by the waveform
approximants. We use a bottom-up approach based on gauge-invariant phase
analysis and the Fisher information matrix to investigate waveform systematics
and help identifying biases in parameter estimation. A mock analysis of 15
different binaries indicates that systematics in current waveform models
dominate over statistical errors at signal-to-noise ratio (SNR) ${\gtrsim} 80$.
This implies biases in the inference of the reduced tidal parameter that are
are larger than the statistical $90\%$ credible-intervals. For example, while
the neutron-star radius could be constrained at ${\sim} 5\%$ level at SNR 80,
systematics can be at the ${\sim} 10\%$ level. We apply our approach to
GW170817 (SNR ${\sim}30$) and confirm that no significant systematic effects
are present. Using an optimal frequency range for the analysis, we estimate a
neutron-star radius of $12.5^{+1.1}_{-1.8}\,$km. The latter is consistent with
an electromagnetic-informed prior and the recent NICER measurement. Exploring
SNR ${\gtrsim}100$ in view of third-generation detectors, we find that all the
current waveform models lead to differences of at least 1-sigma in the
inference of the reduced tidal parameter (for any value of the latter). We
conclude that current waveform models, including those from numerical
relativity, are insufficient to infer the equation of state in the loudest (and
potentially most informative) events that will be observed by advanced and
third generation detectors.
