Regularizing the Pulsar Timing Array likelihood: A path towards Fourier Space

Rutger van Haasteren; Serena Valtolina

REVIEW 3 major objections 4 minor 1 cited by

Regularizing the Pulsar Timing Array likelihood: A path towards Fourier Space

T0 review · 3 major / 4 minor · reviewed 2026-08-11 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read The paper shows that a pulsar timing array likelihood can be regularized into per-pulsar Gaussian Fourier summaries, so the array-level gravitational-wave background search becomes an analytic reweighting of those summaries.

desk verdict A useful two-step PTA workflow with correct Gaussian algebra and real-data validation; the main caveat is the unquantified Gaussian-mixture summary, whose limits the paper itself demonstrates in Appendix C. read the letter →

arxiv 2412.11894 v4 pith:HEL4QVKG submitted 2024-12-16 astro-ph.IM astro-ph.HEgr-qc

Serena Valtolina , Rutger van Haasteren This is my paper · ORCID

classification astro-ph.IMastro-ph.HEgr-qc

keywords pulsartimingarraysgravitationalwavebackgroundFourierdomainlikelihoodregularizationBayesianinferenceHellings-Downscorrelationgamma-raypulsarstwo-stepanalysis

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

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

The reading

This paper tries to make pulsar timing array searches modular: instead of fitting one giant joint model over every pulsar's noise and the gravitational-wave background (GWB), it proposes a two-step pipeline. In Step 1 each pulsar is analyzed alone, and all its timing information is compressed into a Gaussian distribution over Fourier coefficients, with a ridge-regression prior used to make the likelihood normalizable. In Step 2 the full array combines those per-pulsar summaries and reweights them analytically to search for the common, spatially correlated signal. On a 25-pulsar dataset the two-step results match the standard time-domain posteriors for both single-pulsar noise and the GWB amplitude and spectral index. If the method is correct, array-level searches become reusable artifacts that no longer need to re-run each pulsar's detailed noise model.

What carries the argument

The load-bearing object is the regularized per-pulsar likelihood: a ridge-regression prior $p(a|\rho_0)$ with fixed reference hyperparameters is added after marginalizing over timing-model parameters, turning the per-pulsar likelihood into an exact Gaussian $\mathcal{N}(a|\hat{a}_0,\Sigma_0)$ in the Fourier coefficients. Step 2 divides out that same Gaussian through the analytic ratio $p(a|\rho)/p(a|\rho_0)$, so the dummy prior cancels and the true prior, including inter-pulsar correlations, takes over. Assembling $\hat{a}_0$ and $\Sigma_0$ from MCMC samples of the noise parameters uses the covariance update formula of Eq. (19), and the flat-tail power-law reference spectrum keeps the covariance matrices well conditioned.

What would settle it

Simulate a PTA in which a strong deterministic signal overlaps the lowest Fourier frequencies, run both the two-step method and the full time-domain likelihood, and compare the GWB amplitude posteriors; if they differ by more than sampling noise, the Gaussian summary has discarded information the true posterior retains.

Watch

Extended reading notes

Core claim

The central result is Eq. (21): the full-array posterior factorizes as the product of per-pulsar Gaussian summaries $\mathcal{N}(a_k|\hat{a}_{0,k},\Sigma_{0,k})$, a reweighting ratio $p(a|\rho)/p(a|\rho_0)$ that removes the regularization prior, and the array-level prior $p(a|\rho)$ that carries the Hellings-Downs correlation, the quadrupolar angular correlation expected from a gravitational-wave background. This factorization is an analytical identity when the white-noise and deterministic-signal parameters are held fixed, and an approximation when they are marginalized over. The paper's concrete demonstration is that this Fourier-space formulation recovers the same single-pulsar noise and GWB posteriors as the standard time-domain likelihood on the 25-pulsar dataset used for validation.

Load-bearing premise

The claim stands or falls on treating the per-pulsar posterior of the Fourier coefficients, after white-noise and deterministic parameters are marginalized, as a single Gaussian; if any of those parameters is strongly correlated with the Fourier coefficients, the summary is lossy and the reweighting is biased.

Editorial extensions

If this is right

A GWB search can be split into per-pulsar noise analyses and a global Fourier-domain analysis that never touches the raw arrival times.
White-noise and deterministic signals, such as exponential DM dips, can be marginalized over in the per-pulsar step instead of being fixed to best-fit values, which is statistically more complete.
The Fourier summaries stay reusable: changing the array-level model only changes the reweighting factor and the prior, not the per-pulsar processing.
Radio and gamma-ray pulsar data, which currently use different likelihoods, could feed the same Fourier-space summary representation.
On the 25-pulsar dataset, the two-step search reproduces the standard GWB amplitude and spectral-index posteriors.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

The per-pulsar Gaussian summaries could be archived as a standard data product, letting future array-level searches test new common-signal models without redoing single-pulsar noise inference.
A practical safety diagnostic would compare the Gaussian summary's covariance with the sample covariance of Fourier coefficients drawn in the full model; where they diverge, Step 1 should be rerun rather than trusted.
The reference spectrum $\rho_0$ is a free knob, and its choice affects conditioning and summary fidelity, so the particular flat-tail choice used here should be validated on each pulsar rather than assumed universal.
A direct testable extension is to pass gamma-ray photon-to-photon Fourier samples through Eq. (21) and compare a combined radio+gamma-ray GWB posterior with the radio-only one.

Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

This paper gives PTA analysts a modular two-step workflow: per-pulsar Fourier summaries as Gaussian means and covariances, then an array-level reweighted Gaussian likelihood for the GWB search. The central derivation, Eqs. (21)-(23), is clean and correct, and the validation on the 25-pulsar EPTA DR2new dataset against ENTERPRISE is genuine evidence that the pipeline works. The authors also deserve credit for reporting the failed PCA variant and for stating plainly that there is no speedup over existing codes.

The genuinely new piece is the reweighting trick: regularize with a ridge prior at rho0, summarize each pulsar's posterior as a Gaussian, then multiply by the prior ratio p(a|rho)/p(a|rho0) to undo the regularization. That specific derivation is not in the cited literature, and the two-step separation is a real contribution even though all ingredients are standard.

The load-bearing approximation is Eq. (18)-(19): after marginalizing over white noise and deterministic parameters, the per-pulsar posterior of the Fourier coefficients is replaced by a single Gaussian matched to the conditional means and covariances. There is no error estimate or diagnostic for this step. The heuristic that theta are "not strongly correlated" with a is not the right criterion; a mixture can be strongly non-Gaussian even when those correlations vanish. So the abstract's "retains all information" is not established for the marginalized case. The analytic equivalence claim is proved only when theta is held fixed, which the conclusions do state, but the abstract overstates it. The Fig. 2 match is real but visual. I do not agree with the stress-test claim that both compared analyses use the same Gaussian-summary route: the blue curve is the standard time-domain ENTERPRISE likelihood with fixed white noise, so the comparison does provide some validation. It just does not bound the approximation error on noisier or more complex data. The Appendix C PCA failure is the honest warning: the same kind of compression can silently narrow posteriors.

This paper is for PTA methodologists, especially anyone wanting to combine radio and gamma-ray data or modularize noise and GWB inference. It deserves serious peer review. A referee should ask for diagnostics or injection tests on the summary approximation, and for a qualified abstract. I would accept it with revision, not desk reject.

Referee Report

3 major / 4 minor

Summary. The paper proposes a two-step, Fourier-domain reformulation of the pulsar timing array likelihood. Step 1 runs per-pulsar inference over white-noise and deterministic-signal parameters while holding the red-noise/DM/GWB hyperparameters fixed at reference values rho0, and compresses the result into a Gaussian summary N(a|ahat0,Sigma0) for the Fourier coefficients. Step 2 multiplies these per-pulsar summaries by the reweighting ratio p(a|rho)/p(a|rho0) and by the hyperparameter prior, yielding a full-array posterior for the GWB and red-noise parameters; the Fourier coefficients are then marginalized analytically in Eqs. (22)-(23). The authors show that when the white-noise/deterministic parameters theta are held fixed, the construction is algebraically identical to the standard time-domain marginalized likelihood, and they demonstrate the method on the EPTA DR2new dataset: a single-pulsar noise analysis for J1738+0333 (Fig. 1) and a 25-pulsar GWB search (Fig. 2). The main approximation is introduced in Eqs. (18)-(19), where the theta-marginalized conditional posterior of the Fourier coefficients is replaced by a single Gaussian with matched first two moments.

Significance. If the two-step construction is truly information-preserving, this is a useful contribution to PTA methodology. The Gaussian algebra in Eqs. (12)-(23) is clean and standard; the reweighting formula is derived inside the paper rather than assumed; and the comparison against an independent codebase (ENTERPRISE) on real EPTA DR2new data is a concrete implementation check. The authors also make code and tutorials available. The proposed split into per-pulsar summaries and an array-level search is attractive for combining radio and gamma-ray datasets and for reusing expensive per-pulsar noise analyses. However, the 'retains all information' claim currently rests on the unsupported Gaussian-mixture approximation in Eq. (19). Appendix C shows that a closely related moment/PCA compression can silently narrow posteriors on the same class of data, and no diagnostic or injection test is given for Eq. (19). The practical significance of the method therefore depends on closing this gap.

major comments (3)

[Sec. II D, Eqs. (18)-(19)] The final step of Eq. (18)-(19) replaces the true marginal p(a|delta_t,rho0)=integral dtheta N(a|ahat0(theta),Sigma0(theta)) p(theta|delta_t,rho0), which is generally a Gaussian mixture, with the single Gaussian N(a|ahat0,Sigma0) obtained by matching first and second moments. This replacement is exact only under restrictive conditions (e.g., ahat0(theta) linear in theta and Sigma0(theta) independent of theta); the stated heuristic that the white-noise parameters are 'not strongly correlated' with the Fourier coefficients is not the correct criterion, since a symmetric two-component mixture can have zero correlation and still be strongly non-Gaussian. The approximation error propagates directly into Eq. (21) through the reweighting ratio p(a|rho)/p(a|rho0), so it biases the Step-2 GWB posterior whenever the summary is inaccurate. The paper provides no error estimate, no diagnostic, and no injection test for this step. This is load-bearing because the abstract's claim that the method 'retains all information' goes beyond the Conclusions' equivalence statement, which is explicitly restricted to the case where theta is held fixed.
[Appendix C, Fig. 4] The admitted failure of the PCA-based moment compression in Appendix C is directly relevant to the soundness of Eq. (19). On the same EPTA pulsar J1738+0333, the moment-summary approach of Eqs. (C10)-(C11) produces RN posteriors that are visibly narrower than the full time-domain result (Fig. 4). The authors attribute this to poor frequency coverage and non-Gaussianity, but the same pulsar is used in Fig. 1 to illustrate Eq. (19). This does not by itself invalidate Eq. (19), but it removes the default assumption that Gaussian moment summaries are benign in this setting. The manuscript needs either a diagnostic that distinguishes the regime where Eq. (19) is safe from the regime where the Appendix C compression fails, or a restriction of the information-preservation claim to cases satisfying that diagnostic.
[Sec. III B, Fig. 2] The validation in Fig. 2 is visual only and its scope is narrower than the abstract suggests. The blue reference curve is a standard time-domain analysis with the white-noise parameters fixed to single-pulsar maximum-likelihood values, while the orange curve uses the Fourier-domain summary with those parameters marginalized. Agreement between the two therefore demonstrates the fixed-theta identity of Eqs. (20)-(23) together with the practical insignificance of the WN marginalization for this dataset, but it does not validate the information-preservation property of Eq. (19). A quantitative comparison (e.g., posterior overlap or a Bayes-factor check between the two formulations) and, more importantly, an injection-recovery test with a known GWB signal would be needed to support the 'retains all information' wording in the abstract.

minor comments (4)

[Eq. (18)] Equation (18) appears to contain a typo: the line after the proportionality sign repeats the integral twice ('integral dtheta N(...)p(theta) x integral dtheta N(...)p(theta)'), whereas the intended expression is a single integral. Please correct.
[Eqs. (16)-(17) and Table I] The notation for the fixed-theta conditional quantities in Eqs. (16)-(17) is not kept distinct from the marginalized quantities in Eq. (19); the primes/superscripts are introduced and then dropped, and Table I lists only ahat0 and Sigma0. Please introduce consistent symbols such as ahat0(theta) and Sigma0(theta) for the fixed-theta quantities and use them throughout.
[Sec. III A] The text says the regularized-likelihood posteriors are 'equivalent' to a full single-pulsar noise analysis while also noting that they are slightly wider because the white-noise parameters are marginalized. Please rephrase to distinguish 'equivalent in sampling target' from 'numerically identical'.
[Appendix A 2] The sentence 'the distribution at the denominator has to be wider than the condition number at the numerator' presumably means that the covariance matrix of the denominator distribution must be larger in the relevant sense; the reference to 'condition number' is confusing and should be rewritten.

Circularity Check

0 steps flagged · score 1.0 of 10

No significant circularity: the reweighted Gaussian-summary derivation is internal analytic algebra with an explicit approximation, and the Fig. 2 validation is independent of the summary construction.

full rationale

The derivation chain in Eqs. (13)-(23) is self-contained analytic algebra. Eq. (13) is a quadratic-completion identity; Eqs. (15)-(17) follow by marginalizing the timing-model parameters; and Eq. (21) is obtained by inserting the factor p(a|rho0)/p(a|rho0) and then using Eq. (18), which the paper explicitly treats as an approximation rather than an equality. The only lossy step is the moment-matched Gaussian summary of the theta-marginalized distribution in Eqs. (18)-(19), and the paper does not claim it is exact for correlated theta; it says it is 'quite accurate when the parameters theta are not strongly correlated with the Fourier coefficients a' (Sec. II D). That is a correctness and robustness caveat, not a circular input: the summary is not fitted to the GWB posterior, and the reweighting ratio is computed in closed form rather than tuned to force agreement. The validation in Fig. 2 compares the Fourier-domain posterior with the time-domain likelihood implemented in ENTERPRISE (blue curves), which is an independent route, so the agreement is not manufactured by construction. Appendix C's failed PCA experiment is explicitly disclosed as unsatisfactory and is a limitation of a related compression, not evidence that the main reweighting is defined in terms of its own output. The self-citations to earlier likelihood derivations [38-41] are standard references for the time-domain likelihood and are not load-bearing here, because the paper re-derives its own identities and checks against an external codebase. No predicted quantity reduces by construction to a fitted parameter or to a self-citation chain.

Assumptions & free parameters 3 free parameters · 6 assumptions · 0 invented entities

No new physical entity is introduced: the 'regularized likelihood' is a statistical device, not a new signal or force. What the method does add are tuned statistical inputs: the regularization prior rho0 (chosen for numerical stability), the flat-tail spectral floor k (partly motivated by the same stability), the bin grid, and the moment-matched Gaussian summary. The item that carries the most risk is the Gaussian-summary approximation itself, because it is the difference between the abstract's 'retains all information' and what is actually guaranteed, which is retention of the first two moments.

free parameters (3)

Regularization prior rho0 = (log10 A0, gamma0, log10 k0) = log10 A0 = -12, gamma0 = 5, log10 k0 = -5 (hand-chosen)
Fixed in Step 1 for all pulsars; chosen for numerical stability of the reweighting, not fit to the science result. The paper argues the final posterior is insensitive, but rho0 conditions the positivity of Sigma^{-1} in Eq. (23) and the safe range is data-dependent (Appendix A2).
Flat-tail floor k per Gaussian process = Sampled in Step 2 under log10 k in [-9,-4]
The flat-tail power law Eq. (24) is adopted partly to guarantee the numerical stability of the regularization (Appendix A1), so the spectral model is partially load-bearing for the method rather than purely physical.
Number of Fourier bins (30 for RN, 100 for DM) = 30 (RN), 100 (DM) per pulsar
Hand-chosen analysis grid applied uniformly to all 25 pulsars; it sets the dimension of the Gaussian summaries and the cost of Step 2, and is standard PTA practice.

assumptions (6)

domain assumption Timing residuals are linear-Gaussian with white noise: p(delta_t|b,theta) = N(delta_t|Tb, N)
Sec. II B, Eqs. (1)-(8). The entire Fourier-domain derivation in Eqs. (12)-(23) presumes this Gaussian structure; non-Gaussian processes would invalidate the reweighting identities.
domain assumption Improper infinite prior on timing-model offsets xi, B = diag(infinity, phi)
Eq. (9). Standard pulsar-timing practice per [41]; the non-normalizability that motivates the ridge regularization is a direct consequence. If a proper prior were used instead, the likelihood would be normalizable and the regularization unnecessary.
ad hoc to paper The theta-marginalized Fourier-coefficient posterior is a single Gaussian with matched first two moments
Eqs. (18)-(19), Sec. II D. The load-bearing approximation of the method: exact only when theta and a are uncorrelated. The paper gives no error bound, normality diagnostic, or injection test for the EPTA applications.
domain assumption Fourier prior covariance phi is diagonal for single-pulsar noise and HD-correlated across pulsars for the GWB
Eq. (6), Sec. II A. The diagonal form cites [44]; the paper itself notes first-principles corrections in [45,46]. The cross-pulsar term assumes a Gaussian, stationary, isotropic, unpolarized GWB with Hellings-Downs correlations.
domain assumption White noise and deterministic signals (EFAC, EQUAD, DM dips) are not covariant with the GWB modes
Secs. II D and III A. This justifies confining them to Step 1 and is the stated reason the Gaussian summary is accurate. The J1713+0747 DM dip is marginalized in Step 1; residual covariance with low-frequency Fourier modes would bias Step 2.
ad hoc to paper The reweighting ratio p(a|rho)/p(a|rho0) is numerically well-conditioned over the sampled rho
Appendix A2. Requires phi0 wider than phi(rho) so Sigma^{-1} = Sigma0^{-1} + phi^{-1} - phi0^{-1} stays positive definite; the safe range of log10 k0 is 'pulsar specific and depends on the observation cadence', with -5 adopted for EPTA DR2new.

how reviews work

0 comments

Cite this review

Pith. "Pith review of Regularizing the Pulsar Timing Array likelihood: A path towards Fourier Space." pith.science (2026). https://pith.science/paper/HEL4QVKG

@misc{pith2026241211894,
  author       = {Pith},
  title        = {Pith review of: Regularizing the Pulsar Timing Array likelihood: A path towards Fourier Space},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/HEL4QVKG}},
  note         = {Machine review of arXiv:2412.11894}
}

read the original abstract

The recent announcement of evidence for a stochastic background of gravitational waves (GWB) in pulsar timing array (PTA) data has piqued interest across the scientific community. A combined analysis of all currently available data holds the promise of confirming the announced evidence as a solid detection of a GWB. However, the complexity of individual pulsar noise models and the variety of modeling tools used for different types of pulsars present significant challenges for a truly unified analysis. In this work we propose a novel approach to the analysis of PTA data: first a posterior distribution over Fourier modes is produced for each pulsar individually. Then, in a global analysis of all pulsars these posterior distributions can be re-used for a GWB search, which retains all information regarding the signals of interest without the added complexity of the underlying noise models or implementation differences. This approach facilitates combining radio and gamma-ray pulsar data, while reducing the complexity of the model and of its implementations when carrying out a GWB search with PTA data.

Figures

Figures reproduced from arXiv: 2412.11894 by the authors.

**Figure 2.** FIG. 2. Posteriors for the GWB hyperparameters obtained from [PITH_FULL_IMAGE:figures/full_fig_p010_2.png] view at source ↗

**Figure 3.** FIG. 3. Visualization of the PCA approach. The squares are schematic representations of the covariance matrices of [PITH_FULL_IMAGE:figures/full_fig_p015_3.png] view at source ↗

Figures from the paper (1 more)

**Figure 4.** Figure 4: FIG. 4. Posteriors for J [PITH_FULL_IMAGE:figures/full_fig_p016_4.png]

Discussion (0). Continue with ORCID to comment.

Forward citations

Cited by 1 Pith paper

Reviewed papers in the Pith corpus that reference this work. Sorted by Pith novelty score. Full citation record

Bayesian pulsar timing and noise analysis with Vela.jl: an overview
astro-ph.IM 2024-12 conditional novelty 6.0 of 10

A new Julia package, Vela.jl, provides an independent, parallelized Bayesian pulsar timing and noise analysis implementation with a Python interface, validated against PINT and tempo2.

Reference graph

Works this paper leans on

78 extracted references · 68 canonical work pages · cited by 1 Pith paper

[1]

Timing model The timing residuals δt are obtained from the observed Tobs as δt ¼ Tobs − fðt; β0Þ [Eq. (1)], where the term fðt; β0Þ corresponds to the TOAs predicted by the timing model evaluated at the reference values β0 (usually best- guesses from the timing analysis) for the timing parameters β.F r o mE q .(1), we see that δt can be rewritten as a sum...

work page
[2]

calibra- tion parameters,

White noise The TOAs are calculated from averaged pulse profiles by comparing the observed pulse profile with the template profile. From this comparison, we can obtain the TOAs with a certain measurement error σ TOA. If this process were perfect, the root-mean-square of δtWN would be σTOA.D u e to inaccuracies in this process, and other potential instru- ...

work page
[3]

chromatic

Red noise and chromatic noise The turbulent ionized interstellar medium (ISM), the solar wind, and similar effects all influence the pulse propagation differently at different wavelengths, causing a frequency-dependent delay. This effect is modeled as an observing-frequency dependent time-correlated delay δt DM, which is proportional to the square of the ...

work page 2025
[4]

Sesana, A

A. Sesana, A. V ecchio, and C. N. Colacino, Mon. Not. R. Astron. Soc. 390, 192 (2008)

work page 2008
[5]

The end result is an approximation of the posterior distribution as a multivariate Gaussian distribution in the Fourier coefficients a. Step 2 : inference on the hyperparameters of the GWB and all signals that are typically modeled in terms of the Fourier coefficients (RN, DM variations, chro- matic noise, and similar processes). Here, the analysis runs o...

work page
[6]

Let us start with the marginalization over ξ only 2Note that in this section we will first marginalize over ξ, so the quantities b, T will be replaced by a, F

Step 1: Per-pulsar analysis In Step 1 , we analyze each pulsar individually with a ridge regression regularization prior. Let us start with the marginalization over ξ only 2Note that in this section we will first marginalize over ξ, so the quantities b, T will be replaced by a, F. SERENA V ALTOLINA and RUTGER V AN HAASTEREN PHYS. REV . D 112, 043046 (2025...

work page 2025
[7]

In reality, we are actually interested in pða; ρjδtÞ

Step 2: Full array analysis The analysis of Step 1 was carried out with a regulari- zation prior where ρ ¼ ρ0 was held fixed. In reality, we are actually interested in pða; ρjδtÞ. Following Eq. (18),w e can write pðδt; a; ρÞ¼ Z dθpðδtja; θÞpðajρÞpðθÞpðρÞ ≈ Z dθpðδtja; θÞpðajρ0ÞpðθÞpðρÞ pðajρÞ pðajρ0Þ pða; ρjδtÞ ≈ N ðaj ˆa0; Σ0Þ pðajρÞ pðajρ0Þ pðρÞ: ð20Þ W...

work page 2025
[8]

standard

In most contemporary PTA analyses, the white noise parameters θ are typically held fixed for computational efficiency. Although it is formally not correct, because the θ are assumed to not be covariant with ρ and a, holding θ fixed to their maximum likelihood estimators of a single- pulsar analysis is deemed sufficiently accurate for a full PTA analysis. ...

work page 2025

Show all 78 references

[9]

Regularization In classical time-series analysis, the data can be con- verted from the time-domain to the Fourier domain by the discrete Fourier transform (DFT): an invertible linear transformation of the data which we can write as δt ¼ F ˜δt; ðA1Þ where ˜δt is the DFT of the ...
[10]

Numerical resolution The ρ 0 define the regularization of the Fourier coeffi- cients that we use in Sec. II C. Care needs to be taken in choosing the values of ρ0 so that numerical issues like under/overflow and roundoff errors are avoided. The dis- tribution of Fourier coeffi...

2025
[11]

Notebook tutorials and the codes used to produce the figures in this appendix are available at [52]

Implementation: Sampling over the Fourier coefficients We describe here in detail the implementation of a PCA to marginalize over the DM variation hyperparameters when carrying out an inference run on PTA data. Notebook tutorials and the codes used to produce the figures in th...
[12]

percentage of covariance infor- mation

These two quantities not only describe the distribution of the Fourier coefficients aRN when the RN hyperparameters are set to ρRN0 , but also store information about their covariance with the DM hyperparameters ρDM. We describe here how to obtain the best estimates for these ...

2025
[13]

As we mentioned before, we did not get satisfactory results with this approach for some pulsars

Results We present the results obtained with the method pre- sented in Appendix C1 in the case of SPNA on the EPTA DR2new data [49]. As we mentioned before, we did not get satisfactory results with this approach for some pulsars. As an example, we show the results obtained wit...

2025
[14]

Rajagopal and R

M. Rajagopal and R. W. Romani, Astrophys. J. 446, 543 (1995)

1995
[15]

A. H. Jaffe and D. C. Backer, Astrophys. J. 583, 616 (2003)

2003
[16]

J. S. B. Wyithe and A. Loeb, Astrophys. J. 590, 691 (2003)

2003
[17]

Tomita, Prog

K. Tomita, Prog. Theor. Phys. 37, 831 (1967)

1967
[18]

P . A. Rosado, A. Sesana, and J. Gair, Mon. Not. R. Astron. Soc. 451, 2417 (2015)

2015
[19]

D. R. Lorimer and M. Kramer, Handbook of Pulsar Astronomy (Cambridge University Press, Cambridge, UK, 2004), V ol. 4

2004
[20]

R. W. Hellings and G. S. Downs, Astrophys. J. Lett. 265, L39 (1983)

1983
[21]

H. Xu, S. Chen et al. (The CPTA Collaboration), Res. Astron. Astrophys. 23, 075024 (2023)

2023
[22]

Antoniadis et al

J. Antoniadis et al. (EPTA and InPTA Collaborations), Astron. Astrophys. 678, A50 (2023)

2023
[23]

Agazie et al

G. Agazie et al. (The NANOGrav Collaboration),Astrophys. J. Lett. 951, L8 (2023)

2023
[24]

D. J. Reardon, A. Zic et al. (The PPTA Collaboration), Astrophys. J. Lett. 951, L6 (2023)

2023
[25]

M. T. Miles et al. , Mon. Not. R. Astron. Soc. 536, 1489 (2024)

2024
[26]

Agazie et al

G. Agazie et al. (The International Pulsar Timing Array Collaboration), arXiv:2309.00693

arXiv
[27]

J. P . W. V erbiest, S. J. Vigeland, N. K. Porayko, S. Chen, and D. J. Reardon, Results Phys. 61, 107719 (2024)

2024
[28]

Damour and A

T. Damour and A. Vilenkin, Phys. Rev. Lett. 85, 3761 (2000)

2000
[29]

Quelquejay Leclere et al

H. Quelquejay Leclere et al. , Phys. Rev. D 108, 123527 (2023)

2023
[30]

Kerr, Astrophys

M. Kerr, Astrophys. J. 885, 92 (2019)

2019
[31]

Matarrese, O

S. Matarrese, O. Pantano, and D. Saez, Phys. Rev. D 47, 1311 (1993)

1993
[32]

Kosowsky, M

A. Kosowsky, M. S. Turner, and R. Watkins, Phys. Rev. D 45, 4514 (1992)

1992
[33]

Hindmarsh, S

M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Phys. Rev. Lett. 112, 041301 (2014)

2014
[34]

Bartolo, S

N. Bartolo, S. Matarrese, A. Riotto, and A. Väihkönen, Phys. Rev. D 76, 061302 (2007)

2007
[35]

L. A. Boyle and A. Buonanno, Phys. Rev. D 78, 043531 (2008)

2008
[36]

regularization

or FORTY-TWO [37] are designed specifically for radio TOAs, whereas the method of [30] is only able to analyze gamma-ray photon data. Modifying these packages for a combined radio and gamma-ray dataset is highly non- trivial, and requires significant development work. Our work...

2025
[37]

Sorbo, J

L. Sorbo, J. Cosmol. Astropart. Phys. 06 (2011) 003

2011
[38]

Antoniadis et al

J. Antoniadis et al. (EPTA and InPTA Collaborations), arXiv:2306.16227

arXiv
[39]

Afzal et al

A. Afzal et al. (The NANOGrav Collaboration), Astrophys. J. Lett. 951, L11 (2023)

2023
[40]

Fermi-LA T Collaboration, Science 376, 521 (2022)

2022
[41]

Bickel, B

P . Bickel, B. Kleijn, and J. Rice, Astrophys. J. 685, 384 (2008)

2008
[42]

Kerr, Astrophys

M. Kerr, Astrophys. J. 732, 38 (2011)

2011
[43]

Bruel, Astron

P . Bruel, Astron. Astrophys. 622, A108 (2019)

2019
[44]

M. Kerr, A. Parthasarathy, and T. Cromartie (Large Area Telescope Collaboration), Proc. Sci. ICRC2023 ( 2023) 1595

2023
[45]

J. Luo, S. Ransom, P . Demorest, P . S. Ray, A. Archibald, M. Kerr, R. J. Jennings, M. Bachetti, R. van Haasteren, C. A. Champagne, J. Colen, C. Phillips, J. Zimmerman, K. Stovall, M. T. Lam, and F. A. Jenet, Astrophys. J. 911,4 5 (2021)

2021
[46]

Susobhanan et al

A. Susobhanan et al. , arXiv:2405.01977

arXiv
[47]

G. B. Hobbs, R. T. Edwards, and R. N. Manchester, Mon. Not. R. Astron. Soc. 369, 655 (2006)

2006
[48]

R. T. Edwards, G. B. Hobbs, and R. N. Manchester, Mon. Not. R. Astron. Soc. 372, 1549 (2006)

2006
[49]

J. A. Ellis, M. V allisneri, S. R. Taylor, and P . T. Baker, ENTERPRISE: Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE, Astrophysics Source Code Library, record ascl:1912.015 (2019), ascl:1912.015

2019
[50]

R. N. Caballero et al., Mon. Not. R. Astron. Soc. 457, 4421 (2016)

2016
[51]

Lentati, P

L. Lentati, P . Alexander, M. P . Hobson, S. Taylor, J. Gair, S. T. Balan, and R. van Haasteren, Phys. Rev. D 87, 104021 (2013)

2013
[52]

van Haasteren, Y

R. van Haasteren, Y . Levin, P . McDonald, and T. Lu, Mon. Not. R. Astron. Soc. 395, 1005 (2009)

2009
[53]

van Haasteren and Y

R. van Haasteren and Y . Levin, Mon. Not. R. Astron. Soc. 428, 1147 (2012)

2012
[54]

van Haasteren and M

R. van Haasteren and M. V allisneri, Phys. Rev. D 90, 104012 (2014)

2014
[55]

Antoniadis et al

J. Antoniadis et al. (EPTA and InPTA Collaborations), 10.5281/zenodo.8164425 (2023)

2023 doi
[56]

Virtanen et al

P . Virtanen et al. (SciPy 1.0 Contributors), Nat. Methods 17, 261 (2020)

2020
[57]

van Haasteren and M

R. van Haasteren and M. V allisneri, Mon. Not. R. Astron. Soc. 446, 1170 (2015)

2015
[58]

Allen and J

B. Allen and J. D. Romano, Phys. Rev. Lett. 134, 031401 (2025)

2025
[59]

Crisostomi, R

M. Crisostomi, R. van Haasteren, P . M. Meyers, and M. V allisneri,arXiv:2506.13866. REGULARIZING THE PULSAR TIMING ARRA Y LIKELIHOOD: … PHYS. REV . D 112, 043046 (2025) 043046-17

2025 arXiv
[60]

V allisneri and R

M. V allisneri and R. van Haasteren, Mon. Not. R. Astron. Soc. 466, stx069 (2017)

2017
[61]

A. E. Hoerl and R. W. Kennard, Technometrics 12,5 5 (1970)

1970
[62]

Antoniadis et al

J. Antoniadis et al. (EPTA and InPTA Collaborations), Astron. Astrophys. 678, A48 (2023)

2023
[63]

Antoniadis et al

J. Antoniadis et al. (EPTA and InPTA Collaborations), Astron. Astrophys. 678, A49 (2023)

2023
[64]

Ellis and R

J. Ellis and R. van Haasteren, jellis18/ PTMCMCSampler: Official release, 10.5281/zenodo.1037579 (2017)

2017 doi
[65]

https://github.com/serevaltolina/PTA_FourierLikelihood
[66]

Goncharov and S

B. Goncharov and S. Sardana, Mon. Not. R. Astron. Soc. 537, 3470 (2025)

2025
[67]

Thongmeearkom et al

T. Thongmeearkom et al. , Mon. Not. R. Astron. Soc. 530, 4676 (2024)

2024
[68]

C. R. Harris et al. , Nature (London) 585, 357 (2020)

2020
[69]

J. D. Hunter, Comput. Sci. Eng. 9, 90 (2007)

2007
[70]

Foreman-Mackey, J

D. Foreman-Mackey, J. Open Source Software 1,2 4 (2016)

2016
[71]

V allisneri, LIBSTEMPO: Python wrapper for TEMPO2, Astrophysics Source Code Library, record ascl:2002.017 (2020), ascl:2002.017

M. V allisneri, LIBSTEMPO: Python wrapper for TEMPO2, Astrophysics Source Code Library, record ascl:2002.017 (2020), ascl:2002.017

2020
[72]

J. S. Hazboun, LA_FORGE, 10.5281/zenodo.4152550 (2020)

2020 doi
[73]

W. G. Lamb, S. R. Taylor, and R. van Haasteren, Phys. Rev. D 108, 103019 (2023)

2023
[74]

J. W. Cooley and J. W. Tukey, Math. Comput. 19, 297 (1965)

1965
[75]

C. M. Bishop and H. Bishop, Deep Learning: Foundations and Concepts (Springer International Publishing, Cham, 2024)

2024
[76]

Coles, G

W. Coles, G. Hobbs, D. J. Champion, R. N. Manchester, and J. P . W. V erbiest,Mon. Not. R. Astron. Soc. 418, 561 (2011)

2011
[77]

R. V ershynin, High-Dimensional Probability: An Introduc- tion with Applications in Data Science , Cambridge Series in Statistical and Probabilistic Mathematics (Cambridge University Press, Cambridge, England, 2018)

2018
[78]

Jolliffe and J

I. Jolliffe and J. Cadima, Phil. Trans. R. Soc. A 374, 20150202 (2016) . SERENA V ALTOLINA and RUTGER V AN HAASTEREN PHYS. REV . D 112, 043046 (2025) 043046-18

2016

Pith tools

Reviewed August 11, 2026 · model on record in the stance chip above.

[1] [1]

Timing model The timing residuals δt are obtained from the observed Tobs as δt ¼ Tobs − fðt; β0Þ [Eq. (1)], where the term fðt; β0Þ corresponds to the TOAs predicted by the timing model evaluated at the reference values β0 (usually best- guesses from the timing analysis) for the timing parameters β.F r o mE q .(1), we see that δt can be rewritten as a sum...

work page

[2] [2]

calibra- tion parameters,

White noise The TOAs are calculated from averaged pulse profiles by comparing the observed pulse profile with the template profile. From this comparison, we can obtain the TOAs with a certain measurement error σ TOA. If this process were perfect, the root-mean-square of δtWN would be σTOA.D u e to inaccuracies in this process, and other potential instru- ...

work page

[3] [3]

chromatic

Red noise and chromatic noise The turbulent ionized interstellar medium (ISM), the solar wind, and similar effects all influence the pulse propagation differently at different wavelengths, causing a frequency-dependent delay. This effect is modeled as an observing-frequency dependent time-correlated delay δt DM, which is proportional to the square of the ...

work page 2025

[4] [4]

Sesana, A

A. Sesana, A. V ecchio, and C. N. Colacino, Mon. Not. R. Astron. Soc. 390, 192 (2008)

work page 2008

[5] [5]

The end result is an approximation of the posterior distribution as a multivariate Gaussian distribution in the Fourier coefficients a. Step 2 : inference on the hyperparameters of the GWB and all signals that are typically modeled in terms of the Fourier coefficients (RN, DM variations, chro- matic noise, and similar processes). Here, the analysis runs o...

work page

[6] [6]

Let us start with the marginalization over ξ only 2Note that in this section we will first marginalize over ξ, so the quantities b, T will be replaced by a, F

Step 1: Per-pulsar analysis In Step 1 , we analyze each pulsar individually with a ridge regression regularization prior. Let us start with the marginalization over ξ only 2Note that in this section we will first marginalize over ξ, so the quantities b, T will be replaced by a, F. SERENA V ALTOLINA and RUTGER V AN HAASTEREN PHYS. REV . D 112, 043046 (2025...

work page 2025

[7] [7]

In reality, we are actually interested in pða; ρjδtÞ

Step 2: Full array analysis The analysis of Step 1 was carried out with a regulari- zation prior where ρ ¼ ρ0 was held fixed. In reality, we are actually interested in pða; ρjδtÞ. Following Eq. (18),w e can write pðδt; a; ρÞ¼ Z dθpðδtja; θÞpðajρÞpðθÞpðρÞ ≈ Z dθpðδtja; θÞpðajρ0ÞpðθÞpðρÞ pðajρÞ pðajρ0Þ pða; ρjδtÞ ≈ N ðaj ˆa0; Σ0Þ pðajρÞ pðajρ0Þ pðρÞ: ð20Þ W...

work page 2025

[8] [8]

standard

In most contemporary PTA analyses, the white noise parameters θ are typically held fixed for computational efficiency. Although it is formally not correct, because the θ are assumed to not be covariant with ρ and a, holding θ fixed to their maximum likelihood estimators of a single- pulsar analysis is deemed sufficiently accurate for a full PTA analysis. ...

work page 2025