  Binary pulsars are a powerful tool for probing strong gravity that still
outperforms direct gravitational wave observations in a number of directions
due to the remarkable accuracy of the pulsar timing. They can constrain very
precisely the presence of additional charges of the orbiting neutron stars
leading to new channels of energy and angular momentum loss, such as the scalar
dipole radiation. In the present paper, we explore in detail the possibility of
constraining different classes of scalar-Gauss-Bonnet gravity with binary
pulsars. Additionally, the existing constraints related to the observed maximum
mass of neutron stars are also updated. Interestingly, depending on the
equation of state, the resulting limits on the theory coupling parameters can
outperform the constraints coming from binary merger observations by up to a
factor of 2 even for the so-called Einstein-dilaton-Gauss-Bonnet gravity where
neutron stars are often underestimated as relevant theory probes. As an
additional merit, precise Bayesian methods are compared with approximate
approaches with the latter showing very good performance despite their
simplicity.
