Hierarchical Bayesian estimation of population-level torque law parameters from anomalous pulsar braking indices

Abridged. Stochastic fluctuations in the spin frequency $\nu$ of a rotation-powered pulsar affect how accurately one measures the power-law braking index, $n_{\rm pl}$, defined through $\dot{\nu}=K\nu^{n_{\rm pl}}$, and can lead to measurements of anomalous braking indices, with $\vert n \vert = \vert \nu \ddot{\nu}/ \dot{\nu}^{2} \vert \gg1$, where the overdot symbolizes a derivative with respect to time. Previous studies show that the variance of the measured $n$ obeys the predictive, falsifiable formula $\langle n^{2} \rangle = n_{\rm pl}^{2}+\sigma^{2}_{\ddot{\nu}}\nu^{2}\gamma_{\ddot{\nu}}^{-2}\dot{\nu}^{-4}T_{\rm obs}^{-1}$ for $\dot{K}=0$, where $\sigma_{\ddot{\nu}}$ is the timing noise amplitude, $\gamma_{\ddot{\nu}}^{-1}$ is a stellar damping time-scale, and $T_{\rm obs}$ is the total observing time. Here we combine this formula with a hierarchical Bayesian scheme to infer the population-level distribution of $n_{\rm pl}$ for a pulsar population of size $M$. The scheme is validated using synthetic data. For a plausible test population with $M=100$ and injected $n_{\rm pl}$ values drawn from a population-level Gaussian with mean $\mu_{\rm pl}=4$ and standard deviation $\sigma_{\rm pl}=0.5$, intermediate between electromagnetic braking and mass quadrupole gravitational radiation reaction, the Bayesian scheme infers $\mu_{\rm pl}=3.89^{+0.24}_{-0.23}$ and $\sigma_{\rm pl}=0.43^{+0.21}_{-0.14}$. The $M=100$ per-pulsar posteriors for $n_{\rm pl}$ and $\sigma^{2}_{\ddot{\nu}}\gamma_{\ddot{\nu}}^{-2}$ contain $87\%$ and $69\%$, respectively, of the injected values within their $90\%$ credible intervals. Comparable accuracy is achieved for (i) population sizes spanning the range $50 \leq M \leq 300$, and (ii) wide priors satisfying $\mu_{\rm pl} \leq 10^{3}$ and $\sigma_{\rm pl} \leq 10^{2}$, which accommodate plausible spin-down mechanisms with $\dot{K}\neq0$.