Pulsar Timing Array and
Nanohertz Gravitational Wave Detection Console
Live timing residual synthesis with Hellings-Downs angular correlation imposed via Cholesky decomposition of the array covariance, NANOGrav-style power-law characteristic strain, Mollweide-projected pulsar sky map, multi-pulsar pair cross-correlation analysis, and stochastic gravitational wave background detection statistics. Real NANOGrav 15-year array sample with literature timing precisions and astrometric positions in the J2000 ICRS frame.
Overview
Pulsar Timing Arrays use the rotation of millisecond pulsars as galactic-scale clocks. By precisely measuring the time-of-arrival of pulses from many pulsars across the sky over years to decades, we detect tiny correlated perturbations caused by gravitational waves passing through the galaxy. The signal regime is nanohertz - corresponding to GW wavelengths comparable to a light-year and source masses on the order of billions of solar masses (supermassive black hole binaries).
In June 2023, NANOGrav, EPTA, PPTA, and CPTA jointly announced strong evidence for a stochastic gravitational wave background. The smoking gun was the detection of a quadrupolar Hellings-Downs angular correlation between pulsar pairs - a signature unique to gravitational waves and impossible to fake with clock or ephemeris errors.
PTA-01 implements the standard analysis pipeline: synthesize timing residuals across the array with the correct Hellings-Downs covariance imposed via Cholesky decomposition, scatter the resulting pulsar-pair correlations against the theoretical HD curve, and characterise the underlying stochastic background by its power-law amplitude and spectral index. All numerical work runs in Float64. Pulsar positions are given in the J2000 ICRS frame.
Quick Start
Guided Workflows
Workflow A :: Recover the Hellings-Downs Curve
The clean demonstration. Inject a strong GW background and watch correlations align with the theoretical curve.
Workflow B :: NANOGrav 15-year Detection
Recreate the 2023 announcement. Published amplitude is A around 2.4e-15 with gamma near 3.2.
Workflow C :: Single Pulsar Deep Dive
Look at one pulsar's residuals and decompose noise budget.
Tab Guide
Pulsars
Mollweide-projected sky map of 18 NANOGrav-style millisecond pulsars at J2000 positions. Click any marker to select. Selection drives the Residuals tab. Pulsar list below the map gives spin period, period derivative, distance, and timing RMS.
Residuals
Time series of timing residuals for the selected pulsar over a 15-year baseline. Synthesised by summing Fourier modes with HD-correlated amplitudes (Cholesky-decomposed) plus per-pulsar white noise. Sliders control GW amplitude A and observation length.
Correlation
Cross-correlations of all 153 pulsar pairs scattered against angular separation, with the Hellings-Downs theoretical curve overlaid. Detection statistic and binned mean correlation reported.
Spectrum
Log-log timing residual power spectral density. The GW background appears as a steep power-law at low frequencies; pulsar white noise dominates at high frequencies. Sliders for amplitude and spectral index.
Stack / Method
Runtime capability probe and live memory accounting. Method tab gives mathematical and algorithmic documentation with citations.
Controls Reference
| Control | Tab | Action |
|---|---|---|
| Sky map click | Pulsars | Select pulsar for Residuals view. |
| List click | Pulsars | Same as sky map click. |
| A_GW slider | Residuals | Characteristic strain at f = 1/yr. Range 1e-16 to 1e-13. |
| T_obs slider | Residuals | Observation baseline in years. |
| GW ON / OFF | Residuals | Toggle GW component. White noise persists either way. |
| REGENERATE | Residuals/Correlation | Resample stochastic process; redraw scatter. |
| A slider (1e-15) | Spectrum | Strain amplitude in 1e-15 units. |
| gamma slider | Spectrum | Spectral index of timing residual PSD. Range 1 to 7. |
| SET 13/3 | Spectrum | Set gamma to 13/3 (canonical SMBHB). |
| NANOGRAV 15 | Spectrum | Set A=2.4e-15, gamma=3.2 (NANOGrav 15-year posterior centre). |
Reading the Plots
Sky Map
Mollweide equal-area projection. Right ascension increases from 0h on the right edge to 24h on the left edge; central meridian is 12h. Declination from -90 deg (bottom) to +90 deg (top). The dashed horizontal line is the celestial equator. The selected pulsar appears in cyan with a halo; unselected pulsars in dim grey.
Residuals
Vertical axis is timing residual in nanoseconds; horizontal axis is observation epoch in years. Amber dashed lines mark the white noise floor at +/- the literature timing RMS. Slow trends and waves through the band are the GW-induced red noise. For a single pulsar these are visually indistinguishable from instrumental red noise - which is why PTA detection requires the array.
Hellings-Downs Plot
Vertical axis is pair correlation coefficient; horizontal axis is angular separation in degrees. The amber curve is the theoretical HD function. Cyan dots are the 153 measured pair correlations. The dim teal horizontal line is the monopole reference (clock errors). The dim pink curve is the dipole reference (ephemeris errors). Detection appears as the cloud of points tracing the HD shape.
Spectrum
Log-log timing residual PSD. The GW power-law is straight on log-log with slope minus gamma. The amber dashed white-noise floor is flat. The dim pink line is the canonical SMBHB inspiral reference (gamma = 13/3, A = 2.4e-15).
Troubleshooting
- Pair correlations look like a flat cloud, not a curve
- GW amplitude is too low or T_obs is too short. Increase A or T_obs and regenerate. The HD shape emerges only when the GW signal is strong enough to dominate the white noise floor in the cross-correlation.
- Sky map markers overlap and are hard to click
- Use the pulsar list below the map. Each row is clickable. The list and the map are kept in sync.
- Residuals look identical with GW ON or OFF
- Correct behaviour for any single pulsar. The GW background is detectable only by cross-correlating across the array. Switch to the Correlation tab to see the signal.
- Spectrum slope does not visually match the slider value
- Power-law slopes can look misleading if the plot range is short. The slider directly controls gamma; the slope is exact.
- Stack shows WebGPU FALLBACK
- Browser does not advertise WebGPU yet, or the host page lacks COOP/COEP headers. CPU path is fully functional.
Scientific References
- Hellings and Downs 1983 - Upper limits on the isotropic gravitational radiation background from pulsar timing analysis. ApJ 265 L39. Original derivation of the angular correlation function.
- Detweiler 1979 - Pulsar timing measurements and the search for gravitational waves. ApJ 234 1100. Foundational paper on PTA principle.
- Sazhin 1978 - Opportunities for detecting ultralong gravitational waves. Soviet Astronomy 22 36. Earliest PTA proposal.
- NANOGrav 15-year Collaboration 2023 - The NANOGrav 15-year Data Set: Evidence for a Gravitational-Wave Background. ApJ Letters 951 L8.
- EPTA / InPTA Collaboration 2023 - The second data release from the European Pulsar Timing Array. A&A 678 A50.
- PPTA Collaboration 2023 - Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array. ApJ Letters 951 L6.
- Phinney 2001 - A Practical Theorem on Gravitational Wave Backgrounds. arXiv astro-ph 0108028.
- Sesana, Vecchio, Colacino 2008 - The stochastic gravitational-wave background from massive black hole binary systems. MNRAS 390 192.
- Romano and Cornish 2017 - Detection methods for stochastic gravitational-wave backgrounds: a unified treatment. Living Reviews in Relativity 20 2.
- Lentati et al. 2013 - Hyper-efficient model-independent Bayesian method for the analysis of pulsar timing data. Physical Review D 87 104021.
Glossary
| PTA | Pulsar Timing Array. Network of millisecond pulsars used as a galactic-scale gravitational wave detector. |
| MSP | Millisecond pulsar. Recycled neutron star with rotation period under about 30 ms and exceptional spin stability. |
| TOA | Time of arrival. Measured pulse arrival timestamp, referenced to the Solar System barycentre. |
| Timing residual | Difference between observed TOA and timing-model-predicted TOA. |
| HD curve | Hellings-Downs curve. Theoretical pulsar pair correlation versus angular separation under an isotropic stochastic GW background. |
| Stochastic GW background | Incoherent sum of GW signals from many unresolved sources; in PTA band primarily SMBHB inspirals. |
| SMBHB | Supermassive black hole binary. Inspirals at sub-parsec separations radiate nanohertz GWs. |
| Characteristic strain | h_c(f). Dimensionless GW amplitude per logarithmic frequency interval. |
| gamma | Power-law index of the timing residual PSD. gamma = 3 - 2 alpha where alpha is the strain index. SMBHB inspiral predicts gamma = 13/3. |
| NANOGrav | North American Nanohertz Observatory for Gravitational Waves. |
| Cholesky | Matrix factorisation C = L L^T used to impose a target covariance on iid samples. |
| nHz | Nanohertz. 1 nHz corresponds to a wave period of about 31.7 years. |
All numerical results carry a provenance hash recording the array, antenna pattern, integration scheme, frequency basis, and parameter set used to produce them. Persisted to IndexedDB for cross-session reproducibility.
- Array: NANOGrav 15-year sample, 18 millisecond pulsars
- Frame: J2000 ICRS for position; BCRS for time
- Antenna pattern: Hellings and Downs 1983
- Covariance imposition: Cholesky decomposition with 1e-10 ridge regularization
- Frequency basis: Fourier modes m/T_obs for m = 1..30
- Reference detection: NANOGrav 15-year, A around 2.4e-15
Hellings-Downs Antenna Pattern
For an isotropic stochastic gravitational wave background and two pulsars separated by angle theta_ab on the sky, the expected normalised cross-correlation of timing residuals is
HD(theta) = (3/2) x ln(x) - x/4 + 1/2
where x = (1 - cos theta) / 2 = sin^2(theta/2). Boundary values: HD(0) = 1/2 without the pulsar term and 1 with it; HD(180 deg) = 1/4; minimum near theta = 82.5 deg with HD around -0.15. The shape is unique to gravitational waves - clock errors give a monopole correlation, ephemeris errors give a dipole.
Power-Law Stochastic Background
The characteristic strain of an SMBHB-driven background follows h_c(f) = A (f/f_yr)^alpha with f_yr = 1/year. The timing residual power spectral density is then
S_r(f) = A^2 / (12 pi^2) * (f/f_yr)^(-gamma) * f_yr^(-3)
with gamma = 3 - 2 alpha. Pure circular SMBHB inspiral predicts alpha = -2/3 and gamma = 13/3. The NANOGrav 15-year posterior centres at A around 2.4e-15 with gamma around 3.2.
Residual Synthesis with HD Covariance
The console synthesises residuals by Fourier-mode summation with the Hellings-Downs covariance imposed via Cholesky decomposition.
- Build N x N pairwise HD covariance matrix C with C_ii = 1 (pulsar term) and C_ij = HD(theta_ij) for i not equal j
- Add ridge regularization 1e-10 I and Cholesky decompose C = L L^T
- For each Fourier mode m at frequency f_m = m/T_obs:
- Sample iid complex Gaussian z_m of size N (real and imaginary parts independent)
- Apply Cholesky: a_m = L z_m. The vector a_m has covariance C across pulsars.
- Scale by sqrt(S_r(f_m) df) where df = 1/T_obs
- Add cosine and sine projections to each pulsar residual time series
- Add per-pulsar white noise with sigma equal to the literature timing RMS
The resulting time series have the correct cross-correlation structure to leading order. With 30 modes and 200 time samples per pulsar the synthesis cost is below 100 ms on a typical mobile CPU.
Pair Correlation Statistic
For each pair (i, j) the empirical Pearson correlation is computed across the time samples. The detection statistic is the correlation between the binned mean of empirical pair correlations and the HD curve evaluated at the bin centres.
Mollweide Projection
The Mollweide equal-area projection maps spherical (lambda, phi) to (x, y) via the auxiliary angle theta solving 2 theta + sin(2 theta) = pi sin(phi). The console solves this Newton-Raphson per pulsar to convergence tolerance 1e-12.
Compute Pipeline
Cholesky factorization, mode synthesis, and pair correlation all run on Float64 buffers. State is held in Float64Array views over SharedArrayBuffer when cross-origin isolation is available, falling back to plain ArrayBuffer otherwise. Where the runtime advertises WebGPU and SharedArrayBuffer the mode-synthesis kernel can offload to a WGSL compute shader operating on Float64-emulated buffers.