Constraining axion quadratic couplings with the Hulse-Taylor binary system

The orbital evolution of the Hulse-Taylor binary neutron star system is described to high precision by general relativity, in which gravity is the only long-range force and gravitational waves provide the dominant energy-loss channel. We use this precision test of relativistic binary dynamics to derive new constraints on axion couplings to stable neutron star constituents: neutrons, electrons, and muons. Quadratic shift symmetry breaking axion-fermion couplings allow binary systems to lose energy through dipole and quadrupole emission of axion waves. These couplings also mediate long range, spin independent forces in two different regimes: in an ambient DM background, and when a tachyonic phase transition is triggered inside the neutron stars. For light QCD axions our constraints can be recast as limits on the axion decay constant, which are complementary to other probes for $m_a\lesssim 10^{-12}\text{ eV}$ and $f_a\lesssim M_\textrm{pl}$. We also place the strongest constraint available on the axion-muon quadratic coupling, which is otherwise only constrained by supernova cooling.