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 modelling the spin up and down of pulsars self-consistently. We analyse the merger rate R-MW, orbital period P-orb, eccentricity e, spin period P, and spin period derivative (P)over dot 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 (R-MW approximate to 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-orb, e, P, (P)over dot) 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 similar to 20 new BNSs in the Milky Way.
