Isolated pulsar population synthesis with simulation-based inference

TL;DR

This work tackles the problem of inferring the birth properties and magnetorotational evolution of isolated Galactic radio pulsars using simulation-based inference (SBI). It builds a forward pulsar population synthesis model that couples Milky Way dynamics, magnetorotational evolution with a five-parameter magnetic-field model, and radio emission/beaming physics, producing synthetic $P$--$\dot{P}$ diagrams. Neural density estimators are trained in an SBI framework to recover the posterior distributions of $μ_{\log B}$, $σ_{\log B}$, $μ_{\log P}$, $σ_{\log P}$, and $a_{\rm late}$, with an ensemble providing robust posteriors for the observed pulsar population. The inferred results show $μ_{\log B}=13.10^{+0.08}_{-0.10}$, $σ_{\log B}=0.45^{+0.05}_{-0.05}$, $μ_{\log P}=-1.00^{+0.26}_{-0.21}$, $σ_{\log P}=0.38^{+0.33}_{-0.18}$, and $a_{\rm late}=-1.80^{+0.65}_{-0.61}$ (95% credibility), demonstrating SBI’s effectiveness for complex population studies and enabling future multiwavelength pulsar analyses.

Abstract

We combine pulsar population synthesis with simulation-based inference (SBI) to constrain the magnetorotational properties of isolated Galactic radio pulsars. We first develop a framework to model neutron star birth properties and their dynamical and magnetorotational evolution. We specifically sample initial magnetic field strengths, $B$, and spin periods, $P$, from lognormal distributions and capture the late-time magnetic field decay with a power law. Each lognormal is described by a mean, $μ_{\log B}, μ_{\log P}$, and standard deviation, $σ_{\log B}, σ_{\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 we produce a large database of synthetic $P$--$\dot{P}$ diagrams by varying our five magnetorotational 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 $μ_{\log B} = 13.10^{+0.08}_{-0.10}$, $σ_{\log B} = 0.45^{+0.05}_{-0.05}$ and $μ_{\log P} = -1.00^{+0.26}_{-0.21}$, $σ_{\log P} = 0.38^{+0.33}_{-0.18}$ for the lognormal distributions and $a_{\rm late} = -1.80^{+0.65}_{-0.61}$ for the power law at the $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 toward robust statistical inference for complex population synthesis frameworks and forms the basis for future multiwavelength analyses of Galactic pulsars.

Figures (2)

• Figure 1: The key ingredients for pulsar population synthesis. Starting from the bottom left, this approach relies on modeling the neutron stars' dynamical evolution, as well as their magnetorotational properties. For a given beaming geometry and luminosity model, we then determine the pulsars' radio emission and its propagation across the Galaxy toward Earth. For the neutron stars pointing toward us, we subsequently invoke survey limitations and sensitivity thresholds to determine those objects that are detectable. The resulting synthetic populations are compared to the observed ones to constrain input physics.
• Figure 2: Illustration of the $B$-field parameterization used for this study. The five solid curves represent realistic two-dimensional simulations of magnetothermal evolution in the neutron star crust Vigano2021. We fit these together with the late-time power-law evolution of the magnetic field with several broken power laws. The dashed curves shown here are determined for $a_{\rm late} = -3.0$. The colors represent the initial magnetic field strength, $B_0$. To avoid the field decaying to unrealistically small numbers at very late times, we sample the final fields from a Gaussian distribution. The procedure, which allows us to easily extract the dipolar field strength, $B$, at different times, $t$, to study the magnetorotational evolution of our synthetic pulsars, is described in detail in Appendix \ref{['app:B-field']}.