Transdimensional Bayesian approach to pulsar timing noise analysis

The modeling of intrinsic noise in pulsar timing residual data is of crucial importance for gravitational wave detection and pulsar timing (astro)physics in general. The noise budget in pulsars is a collection of several well-studied effects including radiometer noise, pulse-phase jitter noise, dispersion measure variations, and low-frequency spin noise. However, as pulsar timing data continue to improve, nonstationary and non-power-law noise terms are beginning to manifest which are not well modeled by current noise analysis techniques. In this work, we use a transdimensional approach to model these nonstationary and non-power-law effects through the use of a wavelet basis and an interpolation-based adaptive spectral modeling. In both cases, the number of wavelets and the number of control points in the interpolated spectrum are free parameters that are constrained by the data and then marginalized over in the final inferences, thus fully incorporating our ignorance of the noise model. We show that these new methods outperform standard techniques when nonstationary and non-power-law noise is present. We also show that these methods return results consistent with the standard analyses when no such signals are present.

Figure 1

Summary of RJMCMC results from the uniform SNR prior run. In the top panels, the daily averaged residuals are plotted in blue, the injected waveform in green, the MAP waveform (chosen from the model with the highest Bayes factor) in red, and the 90% uncertainty on the waveform (marginalized over all models) in gray. In the bottom-left plot we show the utility as a function of the number of wavelets. In the bottom right, we plot the histogram of the normalized residuals (normalized by the TOA uncertainties) when including (blue) and not including (red) the MAP wavelet model. The green curve is a zero-mean unit variance Gaussian distribution.

Figure 2

Same as Fig. 1, but for the log-uniform wavelet amplitude RJMCMC run.

Figure 3

Waveform reconstruction using a power-law (top) and free-spectrum (bottom) parameterization for the red noise with no additional nonstationary terms. The green is the injected waveform, the red is the MAP recovered waveform, and the gray is the 90% uncertainty on the waveform.

Figure 4

Match between the injected white noise burst waveform and waveform reconstructions for the wavelet and standard noise models. We see that the wavelet models have a higher match overall, thus indicating better fitting, and the power-law and free spectral models have similar matches that are both lower than the wavelet model.

Figure 5

Top: Marginalized posterior distribution of the dimensionless strain amplitude of a stochastic GW background using the transdimensional wavelet model (blue) and standard power-law red noise model (green). Bottom: probability-probability plot for the GW amplitude parameter. On the $x$ axis is the probability $p$ contained in a credible interval, and on the $y$ axis is the fraction of true values that lie within that interval. The black diagonal line is the ideal distribution, where the credible interval is perfectly consistent with the fraction of recovered injections. Again, the blue and green curves represent the wavelet and standard noise models, respectively, with the shaded area accounting for the uncertainty due to the limited number (500) of injections.

Figure 6

Results from simulation with injected GWB and no white noise burst. Top: Marginalized 2D posterior of GWB amplitude and spectral index for a standard model with no wavelets (green), wavelet model marginalized over number of wavelets (blue), and wavelet model using only samples from the RJMCMC that use zero wavelets. The black lines and cross are the injected values. We see that the posteriors are nearly identical in all three cases, showing that there is no evidence of GW absorption by the wavelet model. Bottom: Utility vs number of wavelets for this injection. The data clearly favor zero wavelets.

Figure 7

Results of BayesSpecPTA and the two standard analyses (power-law and free spectral coefficients) on the null case. The top two panels show the recovered spectrum (middle panel) and the favored control points for the spectral interpolation (top panel). The bottom panel is the recovered spectrum from the powerlaw analysis (density) and the free spectral coefficients (error bars). In both density plots, the dashed line is the injected PSD, and the solid black lines represent the median and 90% confidence region. In the bottom panel, the gray error bars represent the median and 90% confidence interval on each spectral coefficient.

Figure 8

Favored number of spectral components from the BayesSpecPTA analysis of the null case. In this case, where the injected spectrum is a simple power law, the spectrum is best described by two parameters as one would expect.

Figure 9

Same as Fig. 7, but for the intermediate case.

Figure 10

Favored number of spectral components from the BayesSpecPTA analysis of the intermediate case. In this case, where the injected spectrum has a singe turnover, the spectrum is best described by three parameters as one would expect, and two parameters (power law) is strongly disfavored while more than three components is allowed.

Figure 11

Same content as Fig. 7 for the extreme case.

Figure 12

Favored number of spectral components from the BayesSpecPTA analysis of the extreme case. In this case, the injected spectrum has two peaks. The data can either be described by a power law or a more complicated spectra that requires $\ge 4$ control points.

Figure 13

Summary plots for combined run. In the top panel, we plot the daily averaged residuals in blue, the MAP red noise and wavelet realization in red and the injected value in green. The gray region is the 90% confidence region of the signal. In the middle panel, we plot the probability density of the recovered PSD and the injected PSD in black. The bottom panels show the utility for the number of wavelets and number of spectral components.