  Tidal interactions in coalescing binary neutron stars modify the dynamics of
the inspiral and hence imprint a signature on their gravitational wave (GW)
signals in the form of an extra phase shift. We need accurate models for the
tidal phase shift in order to constrain the supranuclear equation of state from
observations. In previous studies, GW waveform models were typically
constructed by treating the tide as a linear response to a perturbing tidal
field. In this work, we incorporate non-linear corrections due to hydrodynamic
three- and four-mode interactions and show how they can improve the accuracy
and explanatory power of waveform models. We set up and numerically solve the
coupled differential equations for the orbit and the modes and analytically
derive solutions of the system's equilibrium configuration. Our analytical
solutions agree well with the numerical ones up to the merger and involve only
algebraic relations, allowing for fast phase shift and waveform evaluations for
different equations of state over a large parameter space. We find that, at
Newtonian order, non-linear fluid effects can enhance the tidal phase shift by
$\gtrsim 1\,{\rm radian}$ at a GW frequency of 1000 Hz, corresponding to a
$10-20\%$ correction to the linear theory. The scale of the additional phase
shift near the merger is consistent with the difference between numerical
relativity and theoretical predictions that account only for the linear tide.
Non-linear fluid effects are thus important when interpreting the results of
numerical relativity and in the construction of waveform models for current and
future GW detectors.
