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 linearly proportional to the ejecta-remnant mass ratio. 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 πΌ = 80 Β± 30 km s β1 , π½ = 0 Β± 20 km s β1 as the center of the region of the parameter space that fulfils all of our constraints, and expect π½ β₯ 0 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.