Title: Isolated pulsar population synthesis with simulation-based inference

Abstract: We combine pulsar population synthesis with simulation-based inference (SBI) to constrain the magneto-rotational properties of isolated Galactic radio pulsars. We first develop a framework to model neutron-star birth properties and their dynamical and magneto-rotational evolution. We specifically sample initial magnetic-field strengths, $B$, and spin periods, $P$, from log-normal distributions and capture the late-time magnetic-field decay with a power law. Each log-normal is described by a mean, $\mu_{\log B}, \mu_{\log P}$, and standard deviation, $\sigma_{\log B}, \sigma_{\log P}$, while the power law is characterized by the index, $a_{\rm late}$. We subsequently model the stars' radio emission and observational biases to mimic detections with three radio surveys, and produce a large database of synthetic $P$-$\dot{P}$ diagrams by varying our five magneto-rotational input parameters. We then follow an SBI approach that focuses on neural posterior estimation and train deep neural networks to infer the parameters' posterior distributions. After successfully validating these individual neural density estimators on simulated data, we use an ensemble of networks to infer the posterior distributions for the observed pulsar population. We obtain $\mu_{\log B} = 13.10^{+0.08}_{-0.10}$, $\sigma_{\log B} = 0.45^{+0.05}_{-0.05}$ and $\mu_{\log P} = -1.00^{+0.26}_{-0.21}$, $\sigma_{\log P} = 0.38^{+0.33}_{-0.18}$ for the log-normal distributions, and $a_{\rm late} = -1.80^{+0.65}_{-0.61}$ for the power law at $95\%$ credible interval. We contrast our results with previous studies and highlight uncertainties of the inferred $a_{\rm late}$ value. Our approach represents a crucial step towards robust statistical inference for complex population-synthesis frameworks and forms the basis for future multi-wavelength analyses of Galactic pulsars.