Bayesian Evidence in PTA Analysis

• Bayesian evidence in PTA data analysis is defined as the marginal likelihood integrating likelihood and prior to compare astrophysical models.
• Advanced methods like nested sampling, thermodynamic integration, and normalizing flows enhance convergence and reliability in high-dimensional parameter spaces.
• The application of these techniques in PTA pipelines leads to efficient model discrimination for gravitational-wave backgrounds and exotic signal detections.

Bayesian evidence, or marginal likelihood, is the central model-comparison statistic in the analysis of Pulsar Timing Array (PTA) data. In the Bayesian paradigm, evidence supports or refutes competing models for stochastic gravitational-wave backgrounds (GWB), continuous-wave sources, or exotic signatures such as ultralight axion dark matter, taking full account of noise models, prior choices, and complex signal parameterizations. The computational challenge arises from the high dimensionality (often $D \gtrsim 50$) of PTA likelihoods and highly structured prior spaces. Recent advances, particularly in nested sampling, normalizing flow-based proposals, and thermodynamic integration, have yielded orders-of-magnitude improvements in both convergence speed and evidence reliability. This article summarizes the mathematical framework, computational strategies, and operational performance of Bayesian evidence calculation in PTA inference.

1. Bayesian Evidence: Definition and Role in PTA Inference

The Bayesian evidence $Z$ is the marginal likelihood of the data under a model, integrating the likelihood $L(\theta) = p(d\,|\,\theta)$ weighted by the prior $\pi(\theta)$ over the parameter space: $Z = \int L(\theta)\,\pi(\theta)\,d\theta$ Here, $\theta$ is a vector aggregating timing-model, noise, and signal parameters, such as white and red noise amplitudes, GWB spectrum, and spectral indices. The posterior $p(\theta|d) \propto L(\theta)\,\pi(\theta)$ is used for parameter estimation, while $Z$ enables objective model selection.

Typical PTA analyses define $L(\theta)$ as a multivariate Gaussian in the residuals after subtraction of deterministic timing models, with the covariance dependent on both intrinsic red noise and correlated GWB (or other signal) contributions: $$p(d|\theta) = \frac{\exp\left[ -\frac{1}{2} (d-s(\theta))^{T} C^{-1}(\theta) (d-s(\theta)) \right] }{ \sqrt{(2\pi)^{N_{\mathrm{TOA}}} \det C(\theta)} }$$ where $Z$0 includes noise and GWB terms, and $Z$1 is the model signal. The evidence is typically evaluated over a high-dimensional box prior determined by domain constraints on GW amplitude, spectral index, and noise components (Villa et al., 3 Nov 2025, Kaiser et al., 2022, Li et al., 5 Jun 2025, D'Amico et al., 5 Nov 2025).

Bayes factors, the ratio $Z$2 between models, provide formal model comparison metrics; thresholds for "detection" are set by empirically calibrated values, e.g., $Z$3 for continuous-wave signals (Ellis, 2013).

2. Evidence Estimation Algorithms

Efficient and accurate computation of $Z$4 underlies credible Bayesian detection claims. The following methodologies dominate advanced PTA pipelines:

a. Nested Sampling and Importance Nested Sampling

Standard nested sampling replaces brute-force quadrature by a sequence of "live" points exploring constrained priors with ascending likelihood thresholds. At each iteration, the lowest-likelihood point is removed, replaced by a new prior-constrained sample, and the cumulative volume and evidence estimate are updated: $Z$5 where $Z$6 estimates the enclosed prior volume. The final $Z$7 converges with internal uncertainty $Z$8 in $Z$9 even for $L(\theta) = p(d\,|\,\theta)$0 (Villa et al., 3 Nov 2025).

Importance nested sampling (i-nessai), enhances efficiency by allowing all proposal points to contribute via importance weights, and leverages advanced proposal distributions parameterized by normalizing flows. The live set evolution, flow training on constrainted points, and meta-proposal construction enable adaptive focusing on posterior regions with significant evidence (Villa et al., 3 Nov 2025).

b. Thermodynamic Integration and Steppingstone Methods

Thermodynamic integration (TI) computes

$L(\theta) = p(d\,|\,\theta)$1

with $L(\theta) = p(d\,|\,\theta)$2 parameterizing a "power posterior." Parallel tempering MCMC with a ladder of $L(\theta) = p(d\,|\,\theta)$3 values computes averages across interpolating distributions (Ellis, 2013).

Generalized steppingstone sampling (GSS) improves on TI by introducing a reference density $L(\theta) = p(d\,|\,\theta)$4 approximating the posterior: $L(\theta) = p(d\,|\,\theta)$5 GSS enables lower-temperature chains to sample effectively and requires only $L(\theta) = p(d\,|\,\theta)$6–$L(\theta) = p(d\,|\,\theta)$7 steps for stable results, reducing MCMC cost and bias (Zahraoui et al., 2024).

3. Advanced Proposal Distributions: Normalizing Flows

Normalizing flows (NF) are invertible deep neural networks trained to map simple base distributions to complex target densities. PTA analyses utilize NFs to model:

• The posterior or constrained-prior surfaces required for nested sampling;
• Hierarchical priors (e.g., in hyperparameterized noise models), enabling factorization or reparameterization of the joint prior (Villa et al., 3 Nov 2025, D'Amico et al., 5 Nov 2025).

Flows provide explicit Jacobians, permitting proper weight updates under transformations. In i-nessai, new live points $L(\theta) = p(d\,|\,\theta)$8 are drawn from the learned flow $L(\theta) = p(d\,|\,\theta)$9, subject to $\pi(\theta)$0, significantly boosting the efficiency in high-dimensional spaces (up to $\pi(\theta)$1–$\pi(\theta)$2 higher effective sample rate versus ellipsoidal or rejection-based proposals) (Villa et al., 3 Nov 2025). In hierarchical Bayesian workflows, NFs are trained for both marginal and conditional prior densities, decorrelating physical and nuisance parameter sets, and further improving evidence robustness (D'Amico et al., 5 Nov 2025).

4. Pipeline Architectures and Empirical Performance

PTA Bayesian inference pipelines integrate these evidence estimation tools within modular frameworks such as Enterprise and PTArcade:

• Enterprise + i-nessai: Integrates NF-based nested sampling, extracting parameter bounds and evidence on timescales of minutes-to-hours compared to days-to-weeks for PTMCMC (Villa et al., 3 Nov 2025).
• PTArcade: Implements classical and cross-correlated PTA likelihoods, with nested sampling (e.g., dynesty) or MCMC-based evidence routines for induced SGWB cosmological models (Domènech et al., 2024).
• Standard pipelines (PTMCMC, AM+PT): Employ parallel-tempered Metropolis samplers, adaptive covariance proposals, and TI for evidence, as in continuous-wave analyses (Ellis, 2013).

Typical performance improvements for i-nessai versus PTMCMC are:

• Effective sample rate (ESS/s) increases by $\pi(\theta)$3–$\pi(\theta)$4 (see Table below).
• Evidence uncertainty $\pi(\theta)$5 in $\pi(\theta)$6 is routinely achieved in $\pi(\theta)$7–$\pi(\theta)$8.
• Stability diagnostics (e.g., $\pi(\theta)$9 convergence, corner plots of log-likelihood and weights) confirm robustness.

Pulsars	Sampler	Wall-time (s)	# Likelihood Evals	ESS/s
1	i-nessai	61	42,000	74.8
1	PTMCMC	61	29,762	0.75
2	i-nessai	146	102,000	52.1
2	PTMCMC	146	41,487	0.19
3	i-nessai	360	184,000	25.1
3	PTMCMC	360	75,268	0.08

Further, GSS delivers evidence precision $Z = \int L(\theta)\,\pi(\theta)\,d\theta$0 in log-evidence and enables reliable Bayes-factor discrimination among models, at sharply reduced computational cost versus TI or nested sampling (Zahraoui et al., 2024).

5. Hierarchical and Hypermodel Extensions

PTA datasets demand joint inference over both astrophysical (signal) and complex noise (including red noise, intrinsic and correlated) parameters. Hierarchical Bayesian formalism treats population-level (hyper-)parameters as part of the inference: $Z = \int L(\theta)\,\pi(\theta)\,d\theta$1 where $Z = \int L(\theta)\,\pi(\theta)\,d\theta$2 are hyperparameters. Normalizing flows trained on this full hierarchy allow reparameterizations that reduce posterior correlations, leading to tighter, better-calibrated $Z = \int L(\theta)\,\pi(\theta)\,d\theta$3 estimates and lower prior dependence (D'Amico et al., 5 Nov 2025).

Hypermodel (product-space) methods embed multiple competing models into a single joint inference over an expanded parameter space, allowing rapid Bayes-factor extraction by tracking model-occupancy fractions in posterior samples (Kaiser et al., 2022).

6. Accuracy, Diagnostics, and Recommended Practice

Robust evidence calculation in high-dimensional PTA inference requires:

• Diagnostic convergence tools: live point state plots, evidence evolution, effective sample size per temperature/layer.
• Empirical accuracy validation: replicated runs yield $Z = \int L(\theta)\,\pi(\theta)\,d\theta$4 values consistent with Monte Carlo expectations; parameter recovery biases are negligible compared to prior–likelihood structure (Villa et al., 3 Nov 2025, D'Amico et al., 5 Nov 2025).
• Best practices for GSS: select a reference density from posterior MCMC draws, choose $Z = \int L(\theta)\,\pi(\theta)\,d\theta$5–$Z = \int L(\theta)\,\pi(\theta)\,d\theta$6 $Z = \int L(\theta)\,\pi(\theta)\,d\theta$7-steps, and run multiple independent estimates to ensure $Z = \int L(\theta)\,\pi(\theta)\,d\theta$8 typical inter-model log Bayes factor (Zahraoui et al., 2024).
• Hierarchical workflows: use NF-driven reparameterization to mitigate the impact of prior–hyperparameter coupling, lowering variance in evidence estimates (D'Amico et al., 5 Nov 2025).

7. Application Domains and Example Results

Bayesian evidence estimation underlies detection and characterization of GWBs (supermassive black-hole binaries, primordial GW), continuous sources, and exotic signatures (e.g., ultralight axion dark matter):

• In two-component GWB searches, Bayes factors based on $Z = \int L(\theta)\,\pi(\theta)\,d\theta$9 extracted from nested sampling or hypermodel sampling reveal separability of astrophysical vs. cosmological backgrounds as a function of dataset length and SNR. At 20 years, median BF$\theta$016 for a secondary cosmological component with $\theta$1 (Kaiser et al., 2022).
• Induced-GW primordial cosmology with a free equation of state $\theta$2 exploits nested sampling via PTArcade, marginalizing over the full GW model, with $\theta$3 constrained to $\theta$4 at 95% HPDI in the monochromatic scenario, and constraints on PBH formation mapped directly via the same Bayesian pipeline (Domènech et al., 2024).
• Full PTA runs for stochastic backgrounds (10 pulsars, 52 parameters) using i-nessai reduce inference times from $\theta$51 week to $\theta$613 hours, with stable evidence convergence and robust posteriors (Villa et al., 3 Nov 2025).
• In PTA-PPA synergistic probes of axion dark matter, the Gaussian approximation for the formal likelihood is justified at currently accessible signal strengths, enabling evidence computation via standard nested sampling (Li et al., 5 Jun 2025).

Bayesian evidence techniques reviewed herein are thus foundational to PTA science, providing rigorous and scalable model discrimination across a spectrum of signal hypotheses and noise complexities.