Gravitational-Wave Data Analysis

A laser-interferometric gravitational-wave detector such as LIGO measures the strain h(t), the fractional change in the differential arm length caused by a passing wave. For a compact binary coalescence the dominant quadrupole radiation produces a characteristic “chirp”: rising frequency and amplitude as the two bodies spiral inward. The waveform in the restricted post-Newtonian approximation is governed by a single mass combination, the chirp mass.

Here D_L is the luminosity distance and τ = t_c − t is the time to coalescence. The chirp mass M encapsulates the leading-order dependence on the component masses:

For GW150914 the component masses were approximately m₁ ≈ 36 M_⊙ and m₂ ≈ 29 M_⊙, yielding a chirp mass of roughly 28.3 M_⊙. The orbital frequency sweeps from about 35 Hz (entering the LIGO band) to merger near 150 Hz within the final fraction of a second. The code below generates this chirp waveform using the leading-order frequency evolution.

No detector is silent: seismic motion dominates below ~10 Hz, thermal noise in the mirror suspensions peaks near 40 Hz, and photon shot noise rises at high frequencies. The one-sided power spectral density S_n(f) captures the frequency-dependent noise floor. It is defined as the ensemble-averaged squared modulus of the Fourier-transformed noise:

The Advanced LIGO design sensitivity is often parameterised analytically. A common fit uses a combination of a low-frequency seismic wall, a thermal bump, and a high-frequency shot-noise rise. The amplitude spectral density (ASD) is simply the square root of the PSD and has units of strain / Hz^1/2.

To generate coloured noise with the correct spectral shape, one draws white Gaussian noise in the frequency domain and multiplies by the square root of the PSD before inverse-Fourier transforming back to the time domain. The code below implements the aLIGO design PSD analytically and plots the resulting noise curve.

Matched filtering is the optimal linear filter for detecting a known signal buried in stationary Gaussian noise. It cross-correlates the data stream with a bank of template waveforms, weighting each frequency bin by the inverse noise PSD. The output is the signal-to-noise ratio (SNR) time series. A detection is claimed when the SNR exceeds a preset threshold (typically 8 for single-detector searches).

The expected SNR for a signal with template perfectly matching the true waveform is:

In practice one computes the matched-filter integral via an inverse FFT, giving the SNR at every possible time offset simultaneously. The peak of the SNR time series marks the best-fit arrival time. Below we inject a chirp signal into synthetic coloured noise, perform matched filtering, and recover the SNR peak.

The Fisher information matrix provides a quick estimate of the best achievable parameter uncertainties (the Cramer-Rao bound) without running a full Bayesian analysis. Each element of the matrix is the noise-weighted inner product of waveform partial derivatives with respect to the source parameters θ_i:

The noise-weighted inner product, central to gravitational-wave data analysis, is defined as:

The inverse of the Fisher matrix gives the parameter covariance matrix: Σ_ij = (Γ^-1)_ij. The diagonal elements Σ_ii give the variance on each parameter (the 1-σ uncertainty is √Σ_ii), while off-diagonal elements encode correlations. Below we compute the Fisher matrix for a simplified 4-parameter waveform (chirp mass, symmetric mass ratio, coalescence time, phase) and visualise the resulting error ellipses.

While the Fisher matrix gives a Gaussian approximation to the posterior, real gravitational-wave signals often have multi-modal, correlated posteriors that require full Bayesian inference. By Bayes’ theorem, the posterior probability of the parameters given the data is:

The likelihood function for stationary Gaussian noise is determined entirely by the noise-weighted residual between data and template:

Modern analyses use nested sampling (via dynesty in Bilby) or parallel-tempered MCMC. For pedagogy we implement the simplest Markov Chain Monte Carlo algorithm — the Metropolis-Hastings sampler — on a toy 2-parameter problem: amplitude A and frequency f of a sinusoidal signal in Gaussian noise. The proposal distribution is a simple Gaussian random walk.

Modern gravitational-wave data analysis rests on two major open-source Python stacks. PyCBC provides a complete pipeline for template-bank generation, matched filtering, coincidence testing, and significance estimation. Bilby is the standard Bayesian inference library, wrapping dynestyand emcee samplers with modular likelihood, prior, and waveform interfaces. Both use the LALSuite waveform approximants (IMRPhenom, SEOBNR families) under the hood.

Beyond classical computing, an emerging research direction explores quantum algorithmsfor gravitational-wave analysis. One concrete entry point is amplitude encoding: mapping a discretised waveform into the amplitudes of a quantum state, enabling subsequent quantum inner-product estimation (swap test) or quantum machine learning classifiers. The circuit below encodes a 4-sample waveform snippet into 2 qubits using Qiskit’s state-preparation routine.

The Laser Interferometer Space Antenna (LISA), an ESA-led mission with NASA participation scheduled for launch in the mid-2030s, will open the millihertz gravitational-wave window (0.1 mHz – 0.1 Hz). With 2.5-million-km arms in a heliocentric orbit, LISA will achieve strain sensitivities of order:

LISA’s science targets span three major source classes. Extreme mass-ratio inspirals (EMRIs) — stellar-mass compact objects spiralling into massive black holes — will map the Kerr metric to exquisite precision. Supermassive black hole binaries (SMBHBs) from galaxy mergers at z > 10 will trace the assembly history of cosmic structure. A stochastic background from unresolved galactic binaries and potentially cosmological phase transitions will fill the mHz band.

The code below plots the LISA and aLIGO sensitivity curves together in characteristic strain, overlaid with approximate source population tracks for stellar-mass BBH (LIGO band), EMRIs, and SMBHBs (LISA band). This “waterfall” plot is the standard way to visualise the complementary frequency coverage of ground- and space-based detectors.