  All neutron star progenitors in neutron-star High-Mass X-ray Binaries (NS
HMXBs) undergo a supernova event that may lead to a significant natal kick
impacting the motion of the whole binary system. The space observatory Gaia
performs a deep optical survey with exquisite astrometric accuracy, for both
position and proper motions, that can be used to study natal kicks in NS HMXBs.
  We aim to survey the observed Galactic NS HMXB population and to quantify the
magnitude of the kick imparted onto their NSs, and to highlight any possible
differences arising in between the various HMXB types.
  We perform a census of Galactic NS HMXBs and cross-match existing detections
in X-rays, optical and infrared with the Gaia Early Data Release 3 database. We
retrieve their parallaxes, proper motions, and radial velocities (when
available), and perform a selection based on the quality of the parallax
measurement. We then compute their peculiar velocities with respect to the
rotating reference frame of the Milky Way, and including their respective
masses and periods, we estimate their kick velocities through Markov Chain
Monte Carlo simulations of the orbit undergoing a supernova event.
  We infer the posterior kick distributions of 35 NS HMXBs. After an
inconclusive attempt at characterising the kick distributions with Maxwellian
statistics, we find that the observed NS kicks are best reproduced by a Gamma
distribution of mean $116^{+18}_{-15}$km.s$^{-1}$. We note that supergiant
systems tend to have higher kick velocities than Be High-Mass X-ray Binaries.
The peculiar velocity versus non-degenerate companion mass plane hints at a
similar trend, supergiant systems having a higher peculiar velocity
independently of their companion mass.
