  Galactic binary neutron stars (BNSs) are a unique laboratory to probe the
evolution of BNSs and their progenitors. Here, we use a new version of the
population synthesis code SEVN to evolve the population of Galactic BNSs, by
modeling the spin up and down of pulsars self-consistently. We analyze the
merger rate $\mathcal{R}_{\rm MW}$, orbital period $P_{\rm orb}$, eccentricity
$e$, spin period $P$, and spin period derivative $\dot{P}$ of the BNS
population. Values of the common envelope parameter $\alpha=1 - 3$ and an
accurate model of the Milky Way star formation history best reproduce the BNS
merger rate in our Galaxy ($\mathcal{R}_{\rm MW}\approx{}30$ Myr$^{-1}$). We
apply radio-selection effects to our simulated BNSs and compare them to the
observed population. Using a Dirichlet process Gaussian mixture method, we
evaluate the four-dimensional likelihood in the $(P_{\rm orb}, e, P, \dot{P})$
space, by comparing our radio-selected simulated pulsars against Galactic BNSs.
Our analysis favours an uniform initial distribution for both the magnetic
field ($10^{10-13}$ G) and the spin period ($10-100$ ms). The implementation of
radio selection effects is critical to match not only the spin period and
period derivative, but also the orbital period and eccentricity of Galactic
BNSs. According to our fiducial model, the Square Kilometre Array will detect
$\sim 20$ new BNSs in the Milky Way.
