  Gravitomagnetic quasi-normal modes of neutron stars are resonantly excited by
tidal effects during a binary inspiral, leading to a potentially measurable
effect in the gravitational-wave signal. We take an important step towards
incorporating these effects in waveform models by developing a relativistic
effective action for the gravitomagnetic dynamics that clarifies a number of
subtleties. Working in the slow-rotation limit, we first consider the
post-Newtonian approximation and explicitly derive the effective action from
the equations of motion. We demonstrate that this formulation opens a way to
compute mode frequencies, yields insights into the relevant matter variables,
and elucidates the role of a shift symmetry of the fluid properties under a
displacement of the gravitomagnetic mode amplitudes. We then construct a fully
relativistic action based on the symmetries and a power counting scheme. This
action involves four coupling coefficients that depend on the internal
structure of the neutron star and characterize the key matter parameters
imprinted in the gravitational waves. We show that, after fixing one of the
coefficients by normalization, the other three directly involve the two kinds
of gravitomagnetic Love numbers (static and irrotational), and the mode
frequencies. We discuss several interesting features and dynamical consequences
of this action, and analyze the frequency-domain response function (the
frequency-dependent ratio between the induced flux quadrupole and the external
gravitomagnetic field), and a corresponding Love operator representing the
time-domain response. Our results provide the foundation for deriving precision
predictions of gravitomagnetic effects, and the nuclear physics they encode,
for gravitational-wave astronomy.
