  Natal supernova kicks, the linear momentum compact remnants receive during
their formation, are an essential part of binary population synthesis (BPS)
models. Although these kicks are well-supported by evidence, their underlying
distributions and incorporation into BPS models is uncertain. In this work, we
investigate the nature of natal kicks using a previously proposed analytical
prescription where the strength of the kick is given by
$v_\text{k}=\alpha\frac{m_\text{ejecta}}{m_\text{remnant}}+\beta~\text{km
s}^{-1}$, for free parameters $\alpha$ and $\beta$. We vary the free parameters
over large ranges of possible values, comparing these synthetic populations
simultaneously against four constraints: the merger rate of compact binary
neutron star (BNS) systems, the period-eccentricity distribution of galactic
BNSs, the velocity distribution of single-star pulsars, and the likelihood for
low-ejecta mass supernovae to produce low-velocity kicks. We find that
different samples of the parameter space satisfy each tests, and only 1 per
cent of the models satisfy all four constraints simultaneously. Although we
cannot identify a single best kick model, we report
$\alpha=115^{+40}_{-55}~\text{km s}^{-1}, \beta=15^{+10}_{-15}~\text{km
s}^{-1}$ as the center of the region of the parameter space that fulfils all of
our constraints, and expect $\beta\geq0~\text{km s}^{-1}$ as a further
constraint. We also suggest further observations that will enable future
refinement of the kick model. A sensitive test for the kick model will be the
redshift evolution of the BNS merger rate since this is effectively a direct
measure of the delay-time distribution for mergers. For our best fitting
values, we find that the peak of the BNS merger rate is the present-day.
