  Aims. Neutron stars (NSs) produced in the Milky Way are supposedly ten to the
eighth - ten to the ninth, of which only $\sim 2 \times 10^{3}$ are observed.
Constraining the phase space distribution of NSs may help to characterize the
yet undetected population of stellar remnants. Methods. We perform Monte Carlo
simulations of NS orbits, under different assumptions concerning the Galactic
potential and the distribution of progenitors and birth velocities. We study
the resulting phase space distributions, focusing on the statistical properties
of the NS populations in the disk and in the solar neighbourhood. Results. It
is shown that $\sim 80$ percent of NSs are in bound orbits. The fraction of NSs
located in a disk of radius 20 kpc and width 0.4 kpc is $\lesssim 20$ percent.
Therefore the majority of NSs populate the halo. Fits for the surface density
of the disk, the distribution of heights on the Galactic plane and the velocity
distribution of the disk, are given. We also provide sky maps of the projected
number density in heliocentric Galactic coordinates (l, b). Our results are
compared with previous ones reported in the literature. Conclusions. Obvious
applications of our modelling are in the revisiting of accretion luminosities
of old isolated NSs, the issue of the observability of the nearest NS and the
NS optical depth for microlensing events. These will be the scope of further
studies.
