Discovering pulsars in compact binaries with a hidden Markov model

TL;DR

The paper tackles the computational challenge of finding pulsars in compact binaries with Doppler-modulated $f_{ m p}(t)$ by introducing a fast semi-coherent detection scheme that merges a hidden Markov model with a Schuster periodogram and solves it with the Viterbi algorithm. By discretizing the pulse frequency into $N_Q$ bins and segmenting data into $N_T$ coherent blocks across $N_B$ subbands, the method tracks the frequency evolution efficiently without exhaustive orbital searches. Validation on synthetic Parkes-like data shows detection down to $S \approx 0.50$ mJy for $P_{ m b} \gtrsim 0.012$ day, with runtimes scaling approximately as $T_{ m run} \propto N_B (N_T)(N_Q \ln N_Q)$, enabling practical deployment alongside traditional searches. The approach complements existing techniques, offering a fast, flexible tool for blind searches of millisecond pulsars in binaries and providing groundwork for applying to real data with extensions to colored noise and RFI handling.

Abstract

Discovering radio pulsars in compact binaries, whose orbital periods $P_{\rm b}$ satisfy $P_{\rm b} \lesssim 1 \, \rm{day}$, is computationally challenging, because the time-dependent pulse frequency $f_{\rm p}(t)$ is strongly Doppler modulated by the binary motion. Here we present a new, fast, semi-coherent detection scheme based on a hidden Markov model (HMM) combined with a maximum likelihood matched filter, the Schuster periodogram. The HMM scheme complements traditional acceleration searches by dividing $f_{\rm p}(t)$ into piecewise-constant blocks and tracking the block-to-block evolution efficiently using dynamic programming. Monte Carlo simulations show that the new method can detect compact binaries with flux densities $S \geq 0.50 \, \rm{mJy}$ and orbital periods $P_{\rm b} \geq 0.012 \, \rm{day}$ under observing conditions (e.g.\ cadence) typical of radio pulsar surveys, with and without impulsive, narrowband radio frequency interference. The new method is fast; it employs the classic Viterbi algorithm to solve the HMM recursively. The central processing unit run time scales nominally as $T_{\rm run} \approx 2.8 \, N_B (N_T/10^2) (N_Q \ln N_Q/10^4 \ln 10^4) \, {\rm s}$ for $N_B$ subbands, $N_T$ coherent segments, and $N_Q$ frequency bins.

Figures (11)

• Figure 1: Workflow of the binary pulsar search pipeline for a single subband. The start and end points of the pipeline are in gray ovals. Processes are reported in green rectangles. Inputs and outputs are reported as red and blue parallelograms, respectively. Decision points are drawn as yellow diamonds. The acronym DFT stands for discrete Fourier transforms. For the tests in Section \ref{['Sec:Validation']}, the workflow is repeated for $N_B = 100$ subbands from 50 Hz to 1050 Hz.
• Figure 2: Synthetic Parkes "Murriyang" multibeam system search-mode data, generated using the simulatesearch software package. The data without RFI are 1-bit digitized and constructed from Gaussian radiometer data $n(t)$ added to an injected pulsar signal $f_{\rm p}(t)$, simulated using the simulateSystemNoise, simulateComplexPsr, and createSearchFile subroutines, details of which are given in Section \ref{['SubSec:SyntheticData']}. The radio telescope and binary input parameters are reported in Table \ref{['tab:SourceParameters']}. The observing frequency $f_0$ (MHz) is reported on the vertical axis, centered at 1374 MHz Edwards_2001Manchester_2001. Time (units: seconds) is reported on the horizontal axis. A version of the diagram with RFI included appears in Figure \ref{['Fig:specRFI']}.
• Figure 3: Tail of the noise-only PDF $p(\mathcal{L}$) for $\mathcal{L}> \mathcal{L}_{\rm tail}$, used to set the detection threshold ${\cal L}_{\rm th}$ as a function of the per-subband false alarm probability $\alpha'$ . The gray histogram corresponds to the empirical PDF $p(\mathcal{L}$) for $\mathcal{L}> \mathcal{L}_{\rm tail}$ computed by analyzing $N_{\rm real} = 10^4$ Monte-Carlo realizations of a 10 Hz subband with $N_T = 32$ time bins and $N_Q = 3120$ frequency bins using the Viterbi algorithm in Appendix \ref{['App:Viterbi']}. Overplotted as a red, dashed line is Equation (\ref{['Eq:TailDist']}) for $k = 10^{-4}$ and $\Tilde{\lambda} = 0.51$. The likelihood threshold $\mathcal{L}_{\rm th} = 96.2\pm 0.052$, calculated in Section \ref{['SubSec:NoiseOnly']}, is overplotted as a blue, shaded, vertical region.
• Figure 4: Frequency tracking results for the representative test source in Section \ref{['Sec:Validation']} across the full search band $50 \leq f_0/ (1\, {\rm Hz}) \leq 1050$ (top panel) and in a single subband $220 \leq f_0 / (1\, {\rm Hz}) \leq 230$ which contains an above-threshold outlier with ${\cal L} > {\cal L}_{\rm th}$ (middle and bottom panels). (Top panel.) Histogram of the $100$$\mathcal{L}$ values returned by the Viterbi algorithm for the maximal paths in the $N_B = 100$ subbands for the synthetic noise plus signal validation test in Section \ref{['SubSec:FrequencyTracking']}. Six out of the 100 maximal paths are above-threshold outliers. The black, dashed, vertical line indicates the log-likelihood threshold $\mathcal{L}_{\rm th}$ set in Section \ref{['SubSec:Threshold']}. (Middle panel.) Log-likelihood $\mathcal{L} = \ln P(Q^*|O)$ of the Viterbi paths ending in 3120 frequency bins versus the terminating bin frequency $f_0$ (units: Hz), plotted as a black curve. The blue, dotted, vertical line indicates the injected pulse frequency $f_{\rm p, inj} = 225.02 \, \rm{Hz}$, reported in the top section of Table \ref{['tab:SourceParameters']}. The gray, dotted, horizontal line corresponds to the likelihood threshold $\mathcal{L}_{\rm th} = 96.2 \pm 0.052$, calculated in Section \ref{['SubSec:NoiseOnly']}. (Bottom panel.) Magnified subset of the frequency-time spectrogram with 32$\times$12 pixels whose coloring indicates in a heat map the value of the normalized Fourier power, calculated according to Equation (\ref{['Eq:Power']}). The red, dashed curve is the optimal hidden state sequence $f_{\rm p}(t)$ output by the Viterbi algorithm.
• Figure 5: Minimum detectable flux density $S_{\rm min}$ (units: mJy) for the representative test source in Section \ref{['Sec:Validation']} in a single subband $220 \leq f_0/(1 \, \rm{Hz}) \leq 230$. (Top panel.) Log-likelihood $\mathcal{L} = \ln P(Q^*|O)$ of the optimal Viterbi paths versus the injected flux density $10^{-1} \leq S/(1 \, \rm{mJy}) \leq 2.0$. For every injected $S$, the experiment is repeated using 50 Monte-Carlo realizations of $n(t)$. The median value of $\mathcal{L}$ as a black curve. The dark blue to light blue shaded regions correspond to the 68%, 95%, and 99% credible intervals, respectively. The black, dashed, horizontal line indicates the log-likelihood threshold set in Sections \ref{['SubSec:Threshold']} and \ref{['SubSec:NoiseOnly']}. The injected flux density $S = 1.0 \, \rm{mJy}$ (gray, solid line) as well as the inferred $\mathcal{L} = 345.18$ (gray, dotted line) from the validation tests in Section \ref{['Sec:Validation']} are overplotted by way of comparison. The red, dashed, vertical line indicates the inferred minimum detectable flux density $S_{\rm min} = 0.50 \, \rm{mJy}$. (Bottom panel.) Probability of detection $P_{\rm d}$ versus the injected flux density $10^{-1} \leq S/(1 \, \rm{mJy}) \leq 2.0$ using 50 Monte-Carlo realizations of $n(t)$. Every blue point corresponds to the number of realizations that exceed $\mathcal{L}_{\rm th}$ divided by the total number of Monte-Carlo realizations of $n(t)$. The red, dashed, vertical line indicates $S_{\rm min} = 0.50 \, \rm{mJy}$.