pith. sign in

arxiv: 2601.07592 · v1 · pith:IJX2YTSSnew · submitted 2026-01-12 · 🌌 astro-ph.HE

Timing Gamma-ray Pulsars using Gibbs Sampling

Colin J. Clark , Serena Valtolina , Lars Nieder , Rutger van Haasteren This is my paper

Pith reviewed 2026-05-21 16:01 UTC · model grok-4.3

classification 🌌 astro-ph.HE
keywords gamma-ray pulsarspulsar timingGibbs samplingtiming noiseFermi-LATbinary pulsarsgravitational wave background
0
0 comments X

The pith

A Gibbs sampling method times gamma-ray pulsars by marginalizing over random photon assignments to pulse profile components while fitting timing and noise models.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper introduces a technique for analyzing arrival times of gamma-ray photons from pulsars detected by the Fermi Large Area Telescope. It converts the timing fit into a weighted least-squares problem by randomly assigning each photon to one Gaussian component in a template pulse profile, then uses Gibbs sampling to marginalize over those assignments. This approach estimates timing parameters and stochastic noise processes while incorporating uncertainties in the shape of the pulse profile. Tests on simulated data sets with power-law timing noise show that the recovered parameters remain robust. The same framework is used to model orbital period changes in binary systems and demonstrated on the pulsar B1957+20.

Core claim

By randomly assigning each detected gamma-ray photon to one of the Gaussian components in a pulse profile template and numerically marginalizing those assignments through a Gibbs sampling scheme, the timing analysis of discrete photon data is transformed into an iterative weighted least-squares problem that jointly fits timing model parameters and stochastic noise processes while propagating uncertainties from the unknown pulse shape.

What carries the argument

Gibbs sampling scheme that marginalizes over discrete photon-to-Gaussian-component assignments, converting the timing problem into repeated weighted least-squares fits.

If this is right

  • Enables fitting of power-law timing noise models directly to gamma-ray pulsar data.
  • Supports a Gaussian-process description of orbital period variations in black-widow and redback binaries.
  • Allows gamma-ray pulsars to contribute to searches for the stochastic gravitational wave background.
  • Yields unbiased timing solutions even when the pulse profile shape carries significant uncertainty.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The approach could be adapted to photon-counting data from other instruments beyond the Fermi telescope.
  • It may improve constraints on the timing noise properties of young gamma-ray pulsars.
  • Multi-wavelength timing that combines this gamma-ray method with radio observations could tighten overall pulsar timing array sensitivity.

Load-bearing premise

The pulse profile can be represented as a linear combination of a small number of Gaussian components whose widths and amplitudes are fixed or jointly fitted, and that the marginalization over assignments introduces no bias in the recovered timing or noise parameters.

What would settle it

In simulated Fermi-LAT data sets generated with known power-law timing noise parameters, the method should recover those input noise parameters within the expected posterior uncertainties.

Figures

Figures reproduced from arXiv: 2601.07592 by Colin J. Clark, Lars Nieder, Rutger van Haasteren, Serena Valtolina.

Figure 1
Figure 1. Figure 1: Illustration of the Gibbs sampling procedure applied to the binary MSP J1526−2744. A: photon phases are determined by the current timing model parameters, with probability weights indicated by the greyscale. B: a short Metropolis-Hastings MCMC chain is produced over the template pulse profiles (faint black curves), given these photon phases, and the final step is chosen as the next sample for τ (solid red … view at source ↗
Figure 2
Figure 2. Figure 2: Comparison between the posterior samples obtained via Gibbs sampling (red) and those obtained using the existing emcee method (black), for the binary millisecond pulsar PSR J1526−2744. Template pulse profile parameters are shown in the upper right corner, while timing model parameters are shown in the lower left. There are no strong correlations between pairs of parameters across these two blocks. Contour … view at source ↗
Figure 3
Figure 3. Figure 3: PP-plot for the recovery of the TN amplitude and slope, as defined in Eq. 38, from 200 simulated pulsars. We also show the case of the amplitude recovery from datasets where the slope was fixed to 13/3 (predicted value for a GWB signal). The gray areas are the 1σ and 3σ confidence intervals from the predicted distribution (black diagonal). We simulated 200 pulsars with TN described by a smoothly broken pow… view at source ↗
Figure 4
Figure 4. Figure 4: Results of Gibbs-sampling analysis on PSR B1957+20. Left panels show the photon phases (bottom) and integrated pulse profile (top) according to the best-fitting timing model. Faint black curves on the pulse profile plot show individual posterior samples for the template pulse profile, with the best-fitting model shown in orange. Middle panels illustrate the posterior uncertainty on the photon phases (botto… view at source ↗
Figure 5
Figure 5. Figure 5: Posterior distribution for the GWB and OPV hyperparameters for PSR B1957+20. Contours on the 2-D distributions are at the 1σ and 2σ levels, while dashed vertical lines on the 1-D marginal distributions indicate the 5%, 50% and 95% quantiles. eters to their best-fitting values, but smaller than the AGWB < 6.15 × 10−14 limit obtained with a similarly free timing model. This method also allows us to make quan… view at source ↗
read the original abstract

Timing analyses of gamma-ray pulsars in the Fermi Large Area Telescope data set can provide sensitive probes of many astrophysical processes, including timing noise in young pulsars, orbital period variations in redback binaries, and the stochastic gravitational wave background (GWB). These goals can require careful accounting of stochastic noise processes, but existing methods developed to achieve this in radio pulsar timing analyses cannot be immediately applied to the discrete gamma-ray arrival time data. To address this, we have developed a new method for timing gamma-ray pulsars, in which the timing model fit is transformed into a weighted least squares problem by randomly assigning each photon to an individual Gaussian component of a template pulse profile. These random assignments are then numerically marginalised over through a Gibbs sampling scheme. This method allows for efficient estimation of timing and noise model parameters, while taking into account uncertainties in the pulse profile shape. We simulated Fermi-LAT data sets for gamma-ray pulsars with power-law timing noise processes, showing that this method provides robust estimates of timing noise parameters. We also describe a Gaussian-process model for orbital period variations in black-widow and redback binary systems that can be fit using this new timing method. We demonstrate this method on the black-widow binary millisecond pulsar B1957+20, where the orbital period varies significantly over the LAT data, but which provides one of the most stringent gamma-ray upper limits on the GWB.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 2 minor

Summary. The paper develops a Gibbs sampling procedure to marginalize over discrete assignments of gamma-ray photons to Gaussian components of a pulse-profile template. This converts the timing-model fit to a weighted least-squares problem that jointly estimates timing parameters, power-law timing noise, and profile uncertainties. Simulations of Fermi-LAT data with injected power-law timing noise are used to argue that the method recovers noise parameters robustly; a Gaussian-process model for orbital-period variations is also introduced and demonstrated on the black-widow pulsar B1957+20.

Significance. If the marginalization is shown to be unbiased even under realistic profile mismatch, the method would constitute a useful advance for gamma-ray pulsar timing, particularly for young pulsars with strong timing noise and for binary systems where orbital variations must be separated from a stochastic gravitational-wave background. The explicit handling of profile uncertainty via Gibbs sampling and the provision of a ready-to-use Gaussian-process orbital model are concrete strengths.

major comments (2)
  1. [§4] §4 (Simulation tests): the data sets are generated from the identical multi-Gaussian template model that is later used in the fit. This matched simulation design does not probe bias that would arise when the true pulse profile deviates from the assumed Gaussian-component representation, which is the central assumption required for the claim of robust recovery of timing-noise parameters.
  2. [Abstract and §4] Abstract and §4: the statement that the method 'provides robust estimates' is not accompanied by quantitative diagnostics (bias, coverage of credible intervals, or comparison against a known ground-truth estimator) that would allow the reader to judge the magnitude of any residual bias introduced by the weighted-least-squares approximation or by finite Gibbs mixing.
minor comments (2)
  1. [Methods] The number of Gaussian components and the criteria used to decide when the Gibbs chain has converged are not stated explicitly; these choices should be documented so that the marginalization procedure can be reproduced.
  2. [Methods] Notation for the weights that result from the photon-to-Gaussian assignment step is introduced without a clear equation reference; adding an explicit expression for the effective weight w_i would improve clarity.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for their constructive and detailed comments on our manuscript. These have highlighted important aspects of our simulation design and the presentation of robustness claims. We address each major comment point by point below, indicating the revisions we will incorporate.

read point-by-point responses

  1. Referee: [§4] §4 (Simulation tests): the data sets are generated from the identical multi-Gaussian template model that is later used in the fit. This matched simulation design does not probe bias that would arise when the true pulse profile deviates from the assumed Gaussian-component representation, which is the central assumption required for the claim of robust recovery of timing-noise parameters.

    Authors: We agree that the simulations employ the same multi-Gaussian template for both data generation and model fitting, and therefore do not directly test for bias arising from profile mismatch. This is a fair observation regarding the scope of the robustness claim. In the revised manuscript we will expand the discussion in §4 to explicitly note this modeling assumption and its implications for the method's applicability. We will also add a short supplementary test using a deliberately mismatched template (e.g., a single-Gaussian fit to multi-Gaussian data) to illustrate sensitivity to profile error. These changes will clarify the conditions under which the reported recovery of timing-noise parameters holds. revision: yes

  2. Referee: [Abstract and §4] Abstract and §4: the statement that the method 'provides robust estimates' is not accompanied by quantitative diagnostics (bias, coverage of credible intervals, or comparison against a known ground-truth estimator) that would allow the reader to judge the magnitude of any residual bias introduced by the weighted-least-squares approximation or by finite Gibbs mixing.

    Authors: We accept that the current text lacks explicit quantitative diagnostics to support the 'robust estimates' phrasing. In the revised version we will augment §4 with the following: (i) measured bias and root-mean-square error in the recovered power-law noise amplitude and spectral index across an ensemble of realizations; (ii) empirical coverage fractions for the 68 % and 95 % credible intervals on timing parameters; and (iii) effective sample sizes and autocorrelation times from the Gibbs chains to quantify mixing. These additions will provide readers with concrete metrics for assessing any residual bias from the weighted-least-squares step or finite sampling. revision: yes

Circularity Check

0 steps flagged

No significant circularity in the methodological derivation

full rationale

The paper presents a new numerical procedure that converts gamma-ray pulsar timing into a weighted least-squares problem via random photon-to-Gaussian assignments, then marginalizes the discrete assignments with Gibbs sampling. The central claims concern efficient parameter estimation and robust recovery of timing-noise parameters, supported by simulations generated from the same model class. No derivation step reduces by the paper's own equations to a quantity defined in terms of itself, no fitted input is relabeled as a prediction, and no load-bearing premise rests on a self-citation chain or imported uniqueness theorem. The method is therefore self-contained against the external simulation benchmarks described.

Axiom & Free-Parameter Ledger

0 free parameters · 2 axioms · 0 invented entities

The method rests on the assumption that the pulse profile is well approximated by a sum of Gaussians and that Gibbs sampling correctly marginalizes the discrete assignment variables; no new physical entities are introduced.

axioms (2)
  • domain assumption The observed gamma-ray photons can be treated as independent draws from a time-dependent intensity profile that is a linear combination of fixed Gaussian components.
    This modeling choice is required to turn the timing problem into a weighted least-squares form after random assignment.
  • standard math Gibbs sampling will converge to the correct marginal posterior over timing and noise parameters for the simulated and real data sets considered.
    Standard MCMC convergence assumption invoked to justify the numerical marginalization.

pith-pipeline@v0.9.0 · 5789 in / 1589 out tokens · 61524 ms · 2026-05-21T16:01:11.191154+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

63 extracted references · 63 canonical work pages · 1 internal anchor

  1. [1]

    A., Ackermann, M., Ajello, M., et al

    Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, Science, 325, 840

    work page 2009

  2. [2]

    A., Ajello, M., Allafort, A., et al

    Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS, 208, 17

    work page 2013

  3. [3]

    B., Axelsson, M., et al

    Ajello, M., Atwood, W. B., Axelsson, M., et al. 2021, ApJS, 256, 12

    work page 2021

  4. [4]

    2013, ApJ, 777, L2

    Allafort, A., Baldini, L., Ballet, J., et al. 2013, ApJ, 777, L2

    work page 2013

  5. [5]

    & Romano, J

    Allen, B. & Romano, J. D. 2025, Phys. Rev. Lett., 134, 031401

    work page 2025

  6. [6]

    Applegate, J. H. 1992, ApJ, 385, 621

    work page 1992

  7. [7]

    Applegate, J. H. & Patterson, J. 1987, ApJ, 322, L99

    work page 1987

  8. [8]

    S., & Taylor, J

    Arzoumanian, Z., Fruchter, A. S., & Taylor, J. H. 1994, ApJ, 426, L85

    work page 1994

  9. [9]

    B., Abdo, A

    Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071

    work page 2009

  10. [10]

    Balakrishnan, V ., Freire, P. C. C., Ransom, S. M., et al. 2023, ApJ, 942, L35

    work page 2023

  11. [11]

    2008, ApJ, 685, 384

    Bickel, P., Kleijn, B., & Rice, J. 2008, ApJ, 685, 384

    work page 2008

  12. [12]

    2015, ApJ, 813, L4

    Bochenek, C., Ransom, S., & Demorest, P. 2015, ApJ, 813, L4

    work page 2015

  13. [13]

    J., et al

    Burgay, M., Nieder, L., Clark, C. J., et al. 2024, A&A, 691, A315

    work page 2024

  14. [14]

    J., Di Mauro, M., Wu, J., et al

    Clark, C. J., Di Mauro, M., Wu, J., et al. 2025, ApJ, 994, 149

    work page 2025

  15. [15]

    J., Nieder, L., V oisin, G., et al

    Clark, C. J., Nieder, L., V oisin, G., et al. 2021, MNRAS, 502, 915

    work page 2021

  16. [16]

    J., Pletsch, H

    Clark, C. J., Pletsch, H. J., Wu, J., et al. 2015, ApJ, 809, L2

    work page 2015

  17. [17]

    J., Manchester, R

    Coles, W., Hobbs, G., Champion, D. J., Manchester, R. N., & Verbiest, J. P. W. 2011, MNRAS, 418, 561

    work page 2011

  18. [18]

    R., Gelman, A., & Rubin, D

    Cook, S. R., Gelman, A., & Rubin, D. B. 2006, Journal of Computational and Graphical Statistics, 15, 675

    work page 2006

  19. [19]

    A., Ransom, S

    Corcoran, K. A., Ransom, S. M., Rosenthal, A. C., et al. 2024, arXiv e-prints, arXiv:2412.08688

    work page arXiv 2024

  20. [20]

    Crisostomi, R

    Crisostomi, M., van Haasteren, R., Meyers, P. M., & Vallisneri, M. 2025, arXiv e-prints, arXiv:2506.13866 Csató, L. & Opper, M. 2002, Neural Computation, 14, 641

    work page arXiv 2025

  21. [21]

    2019, Journal of Process Control, 81, 209

    Daemi, A., Kodamana, H., & Huang, B. 2019, Journal of Process Control, 81, 209

    work page 2019

  22. [22]

    S., Ray, P

    Deneva, J. S., Ray, P. S., Camilo, F., et al. 2016, ApJ, 823, 105 Díaz, S. B., Thongmeearkom, T., Phosrisom, A., et al. 2025, MN- RAS[arXiv:2509.09605]

    work page arXiv 2016

  23. [23]

    T., Hobbs, G

    Edwards, R. T., Hobbs, G. B., & Manchester, R. N. 2006, MNRAS, 372, 1549

    work page 2006

  24. [24]

    A., Vallisneri, M., Taylor, S

    Ellis, J. A., Vallisneri, M., Taylor, S. R., & Baker, P. T. 2020, ENTERPRISE: En- hanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE, Zen- odo

    work page 2020

  25. [25]

    K., et al

    Fiori, A., Razzano, M., Harding, A. K., et al. 2024, A&A, 685, A70

    work page 2024

  26. [26]

    W., Lang, D., & Goodman, J

    Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

    work page 2013

  27. [27]

    & Geman, D

    Geman, S. & Geman, D. 1984, IEEE Transactions on Pattern Analysis and Ma- chine Intelligence, PAMI-6, 721

    work page 1984

  28. [28]

    J., Venter, C., et al

    Guillemot, L., Johnson, T. J., Venter, C., et al. 2012, ApJ, 744, 33

    work page 2012

  29. [29]

    & Cornish, N

    Gundersen, A. & Cornish, N. J. 2025, arXiv e-prints, arXiv:2510.12917

    work page arXiv 2025

  30. [30]

    2024, A&A, 692, A170

    Iraci, F., Chalumeau, A., Tiburzi, C., et al. 2024, A&A, 692, A170

    work page 2024

  31. [31]

    S., Johnston, S., Shannon, R

    Kerr, M., Ray, P. S., Johnston, S., Shannon, R. M., & Camilo, F. 2015, ApJ, 814, 128

    work page 2015

  32. [32]

    Larsen, B., Mingarelli, C. M. F., Hazboun, J. S., et al. 2024, ApJ, 972, 49

    work page 2024

  33. [33]

    P., et al

    Lentati, L., Alexander, P., Hobson, M. P., et al. 2014, MNRAS, 437, 3004

    work page 2014

  34. [34]

    P., et al

    Lentati, L., Alexander, P., Hobson, M. P., et al. 2013, Phys. Rev. D, 87, 104021

    work page 2013

  35. [35]

    E., Karastergiou, A., Johnston, S., et al

    Lower, M. E., Karastergiou, A., Johnston, S., et al. 2025, MNRAS, 538, 3104

    work page 2025

  36. [36]

    2021, ApJ, 911, 45

    Luo, J., Ransom, S., Demorest, P., et al. 2021, ApJ, 911, 45

    work page 2021

  37. [37]

    2010, Science, 329, 408

    Lyne, A., Hobbs, G., Kramer, M., Stairs, I., & Stappers, B. 2010, Science, 329, 408

    work page 2010

  38. [38]

    Neal, R. M. 2000, arXiv e-prints, physics/0009028

    work page internal anchor Pith review Pith/arXiv arXiv 2000

  39. [39]

    D., et al

    Ng, C., Bailes, M., Bates, S. D., et al. 2014, MNRAS, 439, 1865

    work page 2014

  40. [40]

    J., Bassa, C

    Nieder, L., Clark, C. J., Bassa, C. G., et al. 2019, ApJ, 883, 42

    work page 2019

  41. [41]

    J., Kandel, D., et al

    Nieder, L., Clark, C. J., Kandel, D., et al. 2020, ApJ, 902, L46

    work page 2020

  42. [42]

    J., et al

    Nieder, L., Kerr, M., Clark, C. J., et al. 2022, ApJ, 931, L3

    work page 2022

  43. [43]

    V ., Ransom, S

    Padmanabh, P. V ., Ransom, S. M., Freire, P. C. C., et al. 2024, A&A, 686, A166

    work page 2024

  44. [44]

    Panin, A. G. & Sokolova, E. V . 2025, A&A, 697, A178

    work page 2025

  45. [45]

    Pletsch, H. J. & Clark, C. J. 2015, ApJ, 807, 18

    work page 2015

  46. [46]

    Rasmussen, C. E. & Williams, C. K. I. 2005, Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning) (The MIT Press)

    work page 2005

  47. [47]

    S., Kerr, M., Parent, D., et al

    Ray, P. S., Kerr, M., Parent, D., et al. 2011, ApJS, 194, 17

    work page 2011

  48. [48]

    S., Polisensky, E., Parkinson, P

    Ray, P. S., Polisensky, E., Parkinson, P. S., et al. 2020, Research Notes of the American Astronomical Society, 4, 37

    work page 2020

  49. [49]

    S., Ransom, S

    Ray, P. S., Ransom, S. M., Cheung, C. C., et al. 2013, ApJ, 763, L13

    work page 2013

  50. [50]

    Ridolfi, A., Freire, P. C. C., Torne, P., et al. 2016, MNRAS, 462, 2918

    work page 2016

  51. [51]

    Roberts, M. S. E. 2013, in IAU Symposium, V ol. 291, Neutron Stars and Pulsars: Challenges and Opportunities after 80 years, ed. J. van Leeuwen, 127–132

    work page 2013

  52. [52]

    C., Ransom, S

    Rosenthal, A. C., Ransom, S. M., Corcoran, K. A., et al. 2025, ApJ, 982, 170

    work page 2025

  53. [53]

    2017, in 2017 Amer- ican Control Conference (ACC), 3134–3140

    Seiferth, D., Chowdhary, G., Mühlegg, M., & Holzapfel, F. 2017, in 2017 Amer- ican Control Conference (ACC), 3134–3140

    work page 2017

  54. [54]

    Shaifullah, G., Verbiest, J. P. W., Freire, P. C. C., et al. 2016, MNRAS, 462, 1029

    work page 2016

  55. [55]

    A., Abdollahi, S., Ajello, M., et al

    Smith, D. A., Abdollahi, S., Ajello, M., et al. 2023, ApJ, 958, 191

    work page 2023

  56. [56]

    C., Chalumeau, A., Tiburzi, C., et al

    Susarla, S. C., Chalumeau, A., Tiburzi, C., et al. 2024, A&A, 692, A18

    work page 2024

  57. [57]

    2025, ApJ, 980, 165

    Susobhanan, A. 2025, ApJ, 980, 165

    work page 2025

  58. [58]

    L., Archibald, A

    Susobhanan, A., Kaplan, D. L., Archibald, A. M., et al. 2024, ApJ, 971, 150

    work page 2024

  59. [59]

    H., Lin, L

    Takata, J., Wang, H. H., Lin, L. C. C., et al. 2020, ApJ, 890, 16

    work page 2020

  60. [60]

    2020, Validat- ing Bayesian Inference Algorithms with Simulation-Based Calibration The Fermi-LAT Collaboration, Ajello, M., Atwood, W

    Talts, S., Betancourt, M., Simpson, D., Vehtari, A., & Gelman, A. 2020, Validat- ing Bayesian Inference Algorithms with Simulation-Based Calibration The Fermi-LAT Collaboration, Ajello, M., Atwood, W. B., et al. 2022, Science, 376, 521

    work page 2020

  61. [61]

    J., Breton, R

    Thongmeearkom, T., Clark, C. J., Breton, R. P., et al. 2024, MNRAS, 530, 4676

    work page 2024

  62. [62]

    & van Haasteren, R

    Valtolina, S. & van Haasteren, R. 2025, Phys. Rev. D, 112, 043046 van Haasteren, R., Levin, Y ., McDonald, P., & Lu, T. 2009, MNRAS, 395, 1005 van Haasteren, R. & Vallisneri, M. 2014, Phys. Rev. D, 90, 104012 van Haasteren, R. & Vallisneri, M. 2015, MNRAS, 446, 1170 V oisin, G., Breton, R. P., & Summers, C. 2020a, MNRAS, 492, 1550 V oisin, G., Clark, C. J...

    work page 2025

  63. [63]

    Wilk, M. B. & Gnanadesikan, R. 1968, Biometrika, 55, 1 Article number, page 13 of 14 A&A proofs:manuscript no. ms Appendix A: Glossary Table A.1.Glossary of mathematical symbols used in this paper. Description Symbol Times-of-arrivalt Timing model parametersρ Pre-fit timing model parametersρ 0 Pre-fit rotational phases according to initial modelϕ 0 Timing...

    work page 1968