Effect on Test Masses

How gravitational waves physically distort spacetime, move freely falling particles, and how laser interferometers convert that motion into a measurable signal — from the geodesic deviation equation through matched filtering and the quantum noise limit.

1. Introduction — How Gravitational Waves Interact with Matter

A gravitational wave is a propagating perturbation of the spacetime metric. Unlike electromagnetic waves, which exert forces on charged particles, gravitational waves act on all matter through the geometry of spacetime itself. A single freely falling particle cannot detect a gravitational wave — it simply follows its geodesic. The observable effect requires at least two separated test masses whose relative displacement changes as the wave passes.

In the transverse-traceless (TT) gauge, the coordinates of freely falling particles do not change. Yet the proper distance between them oscillates. This is the fundamental observable: the gravitational wave strain $h(t)$ modulates the proper separation between test masses. For a wave propagating in the $z$-direction, the metric takes the form:

$ds^2 = -c^2 dt^2 + (1 + h_+)\,dx^2 + (1 - h_+)\,dy^2 + 2h_\times\,dx\,dy + dz^2$

where $h_+(t)$ and $h_\times(t)$ are the two independent polarizations. The plus polarization stretches one transverse axis while compressing the other; the cross polarization does the same but rotated by $45°$. A ring of test masses in the $(x,y)$ plane deforms into an ellipse that oscillates between prolate and oblate orientations at the wave frequency.

Key Physical Points

● A gravitational wave produces a tidal effect: it stretches space along one axis while compressing along the perpendicular axis.
● The fractional change in proper distance (the strain) is $\Delta L / L = \frac{1}{2} h$.
● For the strongest astrophysical sources (binary black hole mergers), $h \sim 10^{-21}$, corresponding to a length change of $\sim 10^{-18}$ m over a 4 km baseline.
● Detection requires extraordinary sensitivity: the displacement is a thousand times smaller than a proton.

The mathematical framework for quantifying these effects is the geodesic deviation equation, which relates the relative acceleration of nearby geodesics to the Riemann curvature tensor. In the sections that follow, we derive this equation in the linearized regime, then show how a Michelson interferometer converts the resulting displacement into a phase shift, how noise limits the sensitivity, and how matched filtering extracts the signal optimally.

2. Derivation: Geodesic Deviation Equation

Consider a one-parameter family of geodesics $x^\mu(\tau, s)$, where $\tau$ is proper time along each geodesic and $s$ labels neighbouring geodesics. The tangent vector is $u^\mu = \partial x^\mu / \partial \tau$ and the deviation vector is $\xi^\mu = \partial x^\mu / \partial s$. The geodesic deviation equation (also called the Jacobi equation) reads:

$\frac{D^2 \xi^\mu}{D\tau^2} = R^\mu_{\ \alpha\beta\gamma}\, u^\alpha\, u^\beta\, \xi^\gamma$

where $D/D\tau$ denotes the covariant derivative along the geodesic. This equation is exact in full general relativity. It tells us that spacetime curvature (encoded in the Riemann tensor $R^\mu_{\ \alpha\beta\gamma}$) causes nearby geodesics to accelerate toward or away from each other.

Step 1: Linearization for Slow-Moving Particles

For a gravitational wave detector on Earth, the test masses (mirrors) move at speeds$v \ll c$. In the detector's rest frame, the four-velocity reduces to:

$u^\mu \approx (c, 0, 0, 0)$

Substituting this into the geodesic deviation equation, the only surviving components have$\alpha = 0$ and $\beta = 0$:

$\frac{D^2 \xi^i}{D\tau^2} = c^2\, R^i_{\ 0j0}\, \xi^j$

where $i, j$ run over the spatial indices 1, 2, 3. Furthermore, in the linearized regime the covariant derivative reduces to an ordinary time derivative (the Christoffel symbols contribute only at second order in $h$), so:

$\ddot{\xi}^i = c^2\, R^i_{\ 0j0}\, \xi^j$

Step 2: Riemann Tensor in TT Gauge

In the TT gauge, the linearized Riemann tensor has the remarkably simple form:

$R^i_{\ 0j0} = -\frac{1}{2c^2}\, \ddot{h}^{TT}_{ij}$

To see this, recall that in linearized gravity the Riemann tensor is:

$R_{\mu\nu\rho\sigma} = \frac{1}{2}\left(\partial_\nu \partial_\rho h_{\mu\sigma} + \partial_\mu \partial_\sigma h_{\nu\rho} - \partial_\mu \partial_\rho h_{\nu\sigma} - \partial_\nu \partial_\sigma h_{\mu\rho}\right)$

Setting $\mu = i$, $\nu = 0$, $\rho = j$,$\sigma = 0$ and using the TT gauge conditions ($h^{TT}_{0\mu} = 0$), only the term $\partial_0 \partial_0 h_{ij}$ survives, giving us the result above.

Step 3: The Master Equation

Combining the two results, we obtain the fundamental equation describing how a gravitational wave affects the separation of test masses:

$\ddot{\xi}^i = \frac{1}{2}\, \ddot{h}^{TT}_{ij}\, \xi^j$

This can be integrated directly. Writing $\xi^i(t) = \xi^i_0 + \delta\xi^i(t)$where $\xi^i_0$ is the unperturbed separation, and working to first order in $h$:

$\delta\xi^i(t) = \frac{1}{2}\, h^{TT}_{ij}(t)\, \xi^j_0$

For the plus polarization with a wave propagating along $z$, the TT metric perturbation is:

$h^{TT}_{ij} = h_+(t)\begin{pmatrix} 1 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & 0 \end{pmatrix} + h_\times(t)\begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}$

For a particle initially displaced along the $x$-axis by $L$, the pure plus-polarization gives:

$\delta\xi^x = \frac{1}{2}\, h_+(t)\, L, \qquad \delta\xi^y = 0$

Similarly, for a particle initially along $y$:

$\delta\xi^x = 0, \qquad \delta\xi^y = -\frac{1}{2}\, h_+(t)\, L$

This is the key result: the $x$-arm stretches while the $y$-arm compresses (and vice versa half a cycle later). The fractional length change is:

$\frac{\Delta L}{L} = \frac{1}{2}\, h(t)$

This is the strain formula that underlies all interferometric gravitational wave detection. For LIGO, with $L = 4$ km and $h \sim 10^{-21}$, the displacement is $\Delta L \sim 2 \times 10^{-18}$ m — about one-thousandth the diameter of a proton.

3. Derivation: Interferometer Response

A laser interferometric detector such as LIGO is a Michelson interferometer with arms along the $x$ and $y$ directions. A laser beam is split at the beam splitter, travels down each arm, reflects off end mirrors, and recombines. The output intensity depends on the phase difference accumulated along the two arms.

Step 1: Round-Trip Phase in the x-Arm

A photon traveling along the $x$-axis in the presence of a plus-polarized gravitational wave follows a null geodesic with $ds^2 = 0$:

$0 = -c^2 dt^2 + \left[1 + h_+(t)\right] dx^2$

This gives $c\,dt = \sqrt{1 + h_+}\, dx \approx \left(1 + \frac{1}{2}h_+\right) dx$to first order in $h$. The round-trip light travel time from $x = 0$to $x = L$ and back is:

$c\,\Delta t_x = 2L + \int_0^{2L/c} \frac{1}{2}\, h_+(t)\, c\, dt' \approx 2L + \frac{1}{2}\, h_+(t)\cdot 2L$

where in the last step we used the fact that the gravitational wave period is much longer than the light travel time ($f_{GW} \ll c/(2L)$), so $h_+$ is approximately constant during one round trip. The accumulated phase for the x-arm is:

$\Phi_x = \frac{2\pi c}{\lambda_\text{laser}} \Delta t_x = \frac{4\pi L}{\lambda_\text{laser}}\left(1 + \frac{1}{2}\, h_+(t)\right)$

Step 2: Round-Trip Phase in the y-Arm

For the $y$-arm, the metric factor is $(1 - h_+)$, so the sign of the perturbation flips:

$\Phi_y = \frac{4\pi L}{\lambda_\text{laser}}\left(1 - \frac{1}{2}\, h_+(t)\right)$

Step 3: Differential Phase

The interferometer output depends on the phase difference:

$\Delta\Phi = \Phi_x - \Phi_y = \frac{4\pi L}{\lambda_\text{laser}}\, h_+(t) = \frac{2\pi L}{\lambda_\text{laser}}\cdot 2h_+(t)$

The output intensity at the antisymmetric port is $I \propto \sin^2(\Delta\Phi / 2)$. For LIGO, with $L = 4$ km, $\lambda_\text{laser} = 1064$ nm, and $h \sim 10^{-21}$, the phase shift is $\Delta\Phi \sim 10^{-11}$ radians. Fabry-Pérot cavities in the arms effectively increase $L$ by a factor of about 300, boosting the signal proportionally.

Step 4: Antenna Pattern Functions

The derivation above assumed the wave arrives along the $z$-axis (perpendicular to the detector plane). For a wave arriving from an arbitrary sky direction $(\theta, \phi)$with polarization angle $\psi$, the detector response is:

$h(t) = F_+(\theta, \phi, \psi)\, h_+(t) + F_\times(\theta, \phi, \psi)\, h_\times(t)$

where the antenna pattern functions for an L-shaped interferometer with arms along $\hat{x}$and $\hat{y}$ are:

$F_+(\theta, \phi, \psi) = \frac{1}{2}(1 + \cos^2\theta)\cos 2\phi\,\cos 2\psi - \cos\theta\,\sin 2\phi\,\sin 2\psi$

$F_\times(\theta, \phi, \psi) = \frac{1}{2}(1 + \cos^2\theta)\cos 2\phi\,\sin 2\psi + \cos\theta\,\sin 2\phi\,\cos 2\psi$

These functions encode the directional sensitivity of the detector. Key properties:

● Maximum response ($|F| = 1$) when the wave arrives perpendicular to the detector plane ($\theta = 0$ or $\pi$) and is aligned with the arms ($\phi = 0$, $\psi = 0$).
● Zero response (blind spots) for waves propagating along the bisector of the arms ($\phi = \pi/4$, $\theta = \pi/2$).
● The combined response $F_+^2 + F_\times^2$ is independent of $\psi$, so the total sensitivity to unpolarized radiation depends only on the source direction.
● Averaged over all sky directions and polarizations: $\langle F_+^2 + F_\times^2 \rangle = 2/5$.

The derivation proceeds by constructing the detector tensor $D^{ij} = \frac{1}{2}(\hat{x}^i\hat{x}^j - \hat{y}^i\hat{y}^j)$and contracting it with the wave-frame polarization tensors rotated into the detector frame via the Euler angles $(\theta, \phi, \psi)$. The polarization tensors in the wave frame are:

$e^+_{ij} = \hat{m}_i \hat{m}_j - \hat{n}_i \hat{n}_j$

$e^\times_{ij} = \hat{m}_i \hat{n}_j + \hat{n}_i \hat{m}_j$

where $\hat{m}$ and $\hat{n}$ are unit vectors perpendicular to the propagation direction $\hat{\Omega}$ and to each other, rotated by $\psi$. The antenna pattern functions are then $F_+ = D^{ij} e^+_{ij}$ and$F_\times = D^{ij} e^\times_{ij}$, which yield the expressions above after explicit computation.

4. Derivation: Strain Sensitivity and Noise

The detectability of a gravitational wave signal depends on its amplitude relative to the detector noise. This competition is quantified by the signal-to-noise ratio (SNR). To derive it, we work in the frequency domain.

Step 1: Frequency-Domain Representation

The detector output is $s(t) = h(t) + n(t)$, where $h(t)$ is the gravitational wave signal and $n(t)$ is noise. We define the Fourier transform convention:

$\tilde{h}(f) = \int_{-\infty}^{\infty} h(t)\, e^{-2\pi i f t}\, dt$

The noise is characterized by its one-sided power spectral density (PSD) $S_n(f)$, defined through:

$\langle \tilde{n}(f)\, \tilde{n}^*(f') \rangle = \frac{1}{2}\, S_n(f)\, \delta(f - f')$

The one-sided PSD has units of $\text{Hz}^{-1}$ (or equivalently,$\text{strain}^2/\text{Hz}$ when referring to strain noise). The characteristic strain noise is $h_n(f) = \sqrt{f\, S_n(f)}$.

Step 2: Noise-Weighted Inner Product

We define a noise-weighted inner product between two time series:

$(a | b) = 4 \,\text{Re} \int_0^{\infty} \frac{\tilde{a}(f)\, \tilde{b}^*(f)}{S_n(f)}\, df$

This inner product has the property that $\langle (n|a)(n|b) \rangle = (a|b)$when $n$ is stationary Gaussian noise with PSD $S_n(f)$.

Step 3: Optimal SNR

The optimal signal-to-noise ratio for a known signal $h(t)$ embedded in stationary Gaussian noise is:

$\text{SNR}^2 = (h | h) = 4 \int_0^{\infty} \frac{|\tilde{h}(f)|^2}{S_n(f)}\, df$

This result follows from the Wiener filter (or equivalently from the matched filter, as we derive in the next section). The integrand $|\tilde{h}(f)|^2 / S_n(f)$ tells us which frequency bands contribute most to the SNR: it is the ratio of signal power to noise power at each frequency.

Step 4: Fundamental Noise Sources

The noise PSD of a laser interferometer has contributions from several fundamental sources:

Shot Noise (High Frequencies)

The quantum nature of light creates photon counting noise. The strain PSD due to shot noise is:

$S_n^{\text{shot}}(f) = \frac{\hbar c \lambda_\text{laser}}{4\pi P_\text{laser} L^2}$

This is frequency-independent (white noise) and decreases with increasing laser power$P_\text{laser}$. More photons means better counting statistics.

Radiation Pressure Noise (Low Frequencies)

The photon momentum fluctuations push the mirrors randomly:

$S_n^{\text{rad}}(f) = \frac{4\hbar P_\text{laser}}{M^2 c \lambda_\text{laser} (2\pi f)^4 L^2}$

This increases with laser power and dominates at low frequencies (scales as $f^{-4}$). The mirror mass $M$ enters because heavier mirrors are harder to push.

Standard Quantum Limit

Shot noise and radiation pressure noise trade off against each other. Minimizing the total quantum noise $S_n^{\text{shot}} + S_n^{\text{rad}}$ with respect to$P_\text{laser}$ gives the standard quantum limit (SQL):

$S_n^{\text{SQL}}(f) = \frac{2\hbar}{M (2\pi f)^2 L^2}$

The SQL is a manifestation of the Heisenberg uncertainty principle applied to the mirror position and momentum. It can be circumvented using squeezed light (as done in Advanced LIGO+).

Seismic Noise (Very Low Frequencies)

Ground vibrations couple to the test masses despite sophisticated seismic isolation. The seismic noise PSD falls steeply with frequency but creates a "noise wall" below about 10 Hz for ground-based detectors:

$S_n^{\text{seis}}(f) \propto f^{-10} \quad \text{(after multi-stage isolation)}$

This is why LIGO is effectively blind below $\sim 10$ Hz. Accessing lower frequencies requires space-based detectors (LISA) or different techniques (pulsar timing arrays).

Thermal Noise (Mid Frequencies)

Brownian motion in the mirror coatings and suspensions contributes noise in the$\sim 10\text{-}100$ Hz band. The fluctuation-dissipation theorem relates this to the mechanical loss angle $\phi_\text{loss}$:

$S_n^{\text{thermal}}(f) \propto \frac{k_B T \, \phi_\text{loss}}{f}$

The total noise PSD is the sum of all contributions: $S_n(f) = S_n^{\text{seis}} + S_n^{\text{thermal}} + S_n^{\text{shot}} + S_n^{\text{rad}} + \ldots$. The detector achieves best sensitivity (lowest $S_n$) in a "sweet spot" around 100–300 Hz for LIGO, where thermal and shot noise are both relatively low.

5. Derivation: Matched Filtering

Matched filtering is the optimal technique for detecting a known signal buried in stationary Gaussian noise. The theory is grounded in the Neyman-Pearson lemma, which states that the likelihood ratio test is the most powerful test at any given false alarm rate.

Step 1: The Likelihood Ratio

We want to decide between two hypotheses: $\mathcal{H}_0$ (noise only,$s = n$) versus $\mathcal{H}_1$ (signal plus noise,$s = h + n$). For stationary Gaussian noise, the log-likelihood ratio is:

$\ln \Lambda = (s | h) - \frac{1}{2}(h | h)$

where we use the noise-weighted inner product defined earlier. The first term is a filter: we correlate the data $s$ with the template $h$. The second term is a normalization constant. The optimal test statistic is thus:

$\hat{\rho} = \frac{(s | h)}{\sqrt{(h | h)}}$

Step 2: Deriving the Optimal Filter

Expanding the inner product in the frequency domain:

$(s | h) = 4 \,\text{Re} \int_0^{\infty} \frac{\tilde{s}(f)\, \tilde{h}^*(f)}{S_n(f)}\, df$

The optimal filter applied to the data in the frequency domain is therefore:

$\tilde{K}(f) = \frac{\tilde{h}^*(f)}{S_n(f)}$

This is the matched filter: it whitens the noise (divides by $S_n(f)$) and then correlates with the expected signal shape ($\tilde{h}^*(f)$). Frequencies where the signal is strong relative to the noise are automatically up-weighted.

Step 3: Expected SNR

When a signal $h$ is present, the expected value of $\hat{\rho}$ is:

$\langle \hat{\rho} \rangle = \sqrt{(h|h)}$

and the variance is $\text{Var}(\hat{\rho}) = 1$ (by construction of the normalization). Therefore the expected SNR is:

$\text{SNR}^2 = (h | h) = 4 \int_0^{\infty} \frac{|\tilde{h}(f)|^2}{S_n(f)}\, df$

This confirms the result from the previous section and shows it is optimal: no other linear filter can achieve a higher SNR.

Step 4: Template Banks and Computational Costs

In practice, the signal waveform $h(t; \vec{\lambda})$ depends on unknown parameters $\vec{\lambda}$ (masses, spins, sky location, etc.). The search must be repeated for a discrete grid of templates covering the parameter space.

The match between a signal with parameters $\vec{\lambda}$ and a template with parameters $\vec{\lambda}'$ is:

$\mathcal{M}(\vec{\lambda}, \vec{\lambda}') = \frac{(h(\vec{\lambda}) | h(\vec{\lambda}'))}{\sqrt{(h(\vec{\lambda}) | h(\vec{\lambda}))(h(\vec{\lambda}') | h(\vec{\lambda}'))}}$

The template bank is constructed so that for any possible signal, at least one template in the bank has a match $\mathcal{M} \geq 1 - \epsilon_{\text{MM}}$, where$\epsilon_{\text{MM}} \sim 0.03$ is the maximum allowed mismatch (corresponding to a 10% loss in detection rate). The number of templates scales as:

$N_{\text{templates}} \sim \frac{V_{\text{param}}}{\epsilon_{\text{MM}}^{d/2}}$

where $d$ is the intrinsic dimensionality of the parameter space and$V_{\text{param}}$ is its volume measured with the Fisher information metric (the "parameter space metric"). For LIGO binary searches:

● Non-spinning BBH: $d = 2$ (two mass parameters), $N_{\text{templates}} \sim 10^5$
● Aligned-spin BBH: $d = 4$, $N_{\text{templates}} \sim 10^6$
● Precessing BBH: $d \geq 7$, $N_{\text{templates}} \sim 10^{10+}$ — computationally prohibitive with brute force

Each template requires an FFT-based correlation with the data, costing $\mathcal{O}(N \log N)$operations where $N$ is the number of data samples. For Advanced LIGO with a sample rate of 16,384 Hz and year-long observing runs, the computational budget can reach$\sim 10^{18}$ floating-point operations. Hierarchical search strategies, reduced-order models, and GPU acceleration are essential for practical searches.

6. Applications

LIGO Noise Budget

The Advanced LIGO design sensitivity curve is a composite of the noise sources discussed above. At design sensitivity, the noise amplitude spectral density reaches a minimum of$\sqrt{S_n} \approx 4 \times 10^{-24}\, \text{Hz}^{-1/2}$ near 100 Hz. The noise budget is typically presented as a plot of $\sqrt{S_n(f)}$ versus frequency, showing the individual contributions:

● Below 10 Hz: Seismic noise dominates completely, creating the low-frequency "wall."
● 10–50 Hz: Suspension thermal noise (pendulum modes of the mirror suspensions).
● 50–200 Hz: Coating Brownian noise (the dominant thermal noise source in the sweet spot).
● Above 200 Hz: Shot noise, which rises as a flat floor in strain PSD.

Quantum Noise Limit and Squeezed Light

The standard quantum limit (SQL) arises from the trade-off between shot noise and radiation pressure noise. In the language of quantum optics, the vacuum fluctuations entering the antisymmetric port of the interferometer set this limit. The SQL in displacement is:

$x_{\text{SQL}} = \sqrt{\frac{2\hbar}{M(2\pi f)^2}} \approx 2 \times 10^{-20}\,\text{m}\,\sqrt{\frac{40\,\text{kg}}{M}\left(\frac{100\,\text{Hz}}{f}\right)^2}$

Squeezed light injection circumvents the SQL by replacing the vacuum fluctuations with a squeezed vacuum state, reducing quantum noise in the measured quadrature. Advanced LIGO has operated with frequency-independent squeezing since O3, achieving about 3 dB of broadband quantum noise reduction. The A+ upgrade uses frequency-dependent squeezing (via a 300-m filter cavity) to reduce shot noise at high frequencies and radiation pressure noise at low frequencies simultaneously.

Signal Processing in Practice

In a real LIGO search pipeline, matched filtering involves several additional steps beyond the basic theory:

● Data conditioning: The raw strain data is whitened, band-passed, and cleaned of instrumental artifacts (glitches) before filtering.
● PSD estimation: The noise PSD $S_n(f)$ is estimated from the data itself using Welch's method (averaged, overlapping periodograms), updated every few seconds to track non-stationarity.
● Trigger generation: The matched filter output is thresholded (typically at SNR $> 4$ per detector) to generate candidate triggers.
● Coincidence: Triggers must appear in multiple detectors within the light travel time (10 ms for LIGO Hanford-Livingston) with consistent parameters.
● Background estimation: The false alarm rate is computed by time-sliding the detector data streams relative to each other, destroying any real coincidences.
● Significance: A detection is claimed when the false alarm rate is below $\sim 1$ per $10^5$ years (corresponding to $\sim 5\sigma$).

7. Historical Context

Weber Bars (1960s)

Joseph Weber pioneered experimental gravitational wave detection with resonant bar detectors in the 1960s. A Weber bar is a large aluminum cylinder (typically 1–2 tons) instrumented with piezoelectric sensors. A passing gravitational wave would excite the fundamental longitudinal mode of the bar (at $\sim 1660$ Hz for Weber's original bars). Weber claimed detections in 1969, but other groups could not reproduce his results. The sensitivity of bar detectors ($h \sim 10^{-16}$) was later shown to be far insufficient for astrophysical sources. Nevertheless, Weber's work launched the field of experimental gravitational wave physics.

Weiss's Interferometer Concept (1972)

Rainer Weiss published a landmark MIT technical report in 1972 analyzing the use of a laser Michelson interferometer for gravitational wave detection. He systematically identified all major noise sources — seismic, thermal, shot noise, radiation pressure — and showed that a kilometer-scale interferometer could reach the sensitivity needed to detect astrophysical gravitational waves. This analysis became the blueprint for LIGO. The key insight was that the interferometer's sensitivity scales with arm length ($\Delta\Phi \propto L \cdot h$), while most noise sources grow more slowly, making longer arms advantageous.

LIGO Proposal (1989) and Construction

Kip Thorne (Caltech) and Rainer Weiss (MIT) led the proposal for LIGO (Laser Interferometer Gravitational-Wave Observatory) to the National Science Foundation in 1989. The project was approved in 1990 — the largest single investment the NSF had ever made. Two identical 4-km interferometers were built: one in Hanford, Washington, and one in Livingston, Louisiana. The separation of $\sim 3000$ km enables coincidence analysis and rudimentary source localization. Initial LIGO operated from 2002–2010 without a detection, reaching its design sensitivity of $h \sim 10^{-22}$.

Advanced LIGO and GW150914

Advanced LIGO, with a factor of $\sim 10$ improvement in strain sensitivity (and $\sim 1000\times$ in volume surveyed), began observations in September 2015. On September 14, 2015 — just two days into the first observing run — the two detectors simultaneously recorded the inspiral, merger, and ringdown of a$36 + 29 \, M_\odot$ binary black hole system at a distance of$\sim 410$ Mpc. The signal, designated GW150914, matched the predictions of general relativity with extraordinary precision. The peak strain was$h \sim 10^{-21}$, and the event was detected with a combined SNR of 24. The 2017 Nobel Prize in Physics was awarded to Weiss, Thorne, and Barish for this discovery.

8. Python Simulation

The simulation below computes two visualizations. First, it shows the deformation of a ring of test particles over one full gravitational wave period for both the plus ($+$) and cross ($\times$) polarizations, plotting snapshots at eight equally spaced phases. Second, it computes the antenna pattern functions $F_+$ and$F_\times$ as sky maps, displaying the directional sensitivity of an L-shaped interferometer. All computation uses NumPy only.

Test Mass Ring Deformation and Antenna Pattern Sky Maps

Python

script.py127 lines

import numpy as np

# ============================================================
# Part 1: Ring of Test Particles Deformation
# ============================================================
# A ring of N particles in the (x, y) plane is deformed by
# a gravitational wave propagating along z.

N_particles = 60
theta_ring = np.linspace(0, 2 * np.pi, N_particles, endpoint=False)
x0 = np.cos(theta_ring)
y0 = np.sin(theta_ring)

# Strain amplitude (exaggerated for visibility)
h0 = 0.35

# Phase snapshots over one full period
N_phases = 8
phases = np.linspace(0, 2 * np.pi, N_phases, endpoint=False)

print("RING_DEFORMATION_DATA")
print("=" * 60)

# Plus polarization: h_ij = h0 cos(phase) * diag(1, -1)
print("PLUS_POLARIZATION")
for k, phi in enumerate(phases):
    h_plus = h0 * np.cos(phi)
    # delta_x = 0.5 * h_plus * x0, delta_y = -0.5 * h_plus * y0
    x_def = x0 * (1 + 0.5 * h_plus)
    y_def = y0 * (1 - 0.5 * h_plus)
    phase_deg = phi * 180.0 / np.pi
    coords = " ".join(f"{x_def[i]:.6f},{y_def[i]:.6f}" for i in range(N_particles))
    print(f"phase={phase_deg:.1f} {coords}")

# Cross polarization: h_ij = h0 cos(phase) * [[0,1],[1,0]]
print("CROSS_POLARIZATION")
for k, phi in enumerate(phases):
    h_cross = h0 * np.cos(phi)
    # delta_x = 0.5 * h_cross * y0, delta_y = 0.5 * h_cross * x0
    x_def = x0 + 0.5 * h_cross * y0
    y_def = y0 + 0.5 * h_cross * x0
    phase_deg = phi * 180.0 / np.pi
    coords = " ".join(f"{x_def[i]:.6f},{y_def[i]:.6f}" for i in range(N_particles))
    print(f"phase={phase_deg:.1f} {coords}")

# ============================================================
# Part 2: Antenna Pattern Functions
# ============================================================
# F_+(theta, phi, psi) and F_x(theta, phi, psi) for an L-shaped
# interferometer with arms along x and y.
# We set psi = 0 and compute over the sky.

N_theta = 91
N_phi = 181

theta_arr = np.linspace(0, np.pi, N_theta)
phi_arr = np.linspace(0, 2 * np.pi, N_phi)
THETA, PHI = np.meshgrid(theta_arr, phi_arr, indexing='ij')

psi = 0.0
cos_psi = np.cos(2 * psi)
sin_psi = np.sin(2 * psi)

# Antenna pattern functions
F_plus = (0.5 * (1 + np.cos(THETA)**2) * np.cos(2 * PHI) * cos_psi
          - np.cos(THETA) * np.sin(2 * PHI) * sin_psi)

F_cross = (0.5 * (1 + np.cos(THETA)**2) * np.cos(2 * PHI) * sin_psi
           + np.cos(THETA) * np.sin(2 * PHI) * cos_psi)

# Combined response (independent of psi)
F_combined = np.sqrt(F_plus**2 + F_cross**2)

print()
print("ANTENNA_PATTERN_DATA")
print("=" * 60)
print(f"N_theta={N_theta} N_phi={N_phi}")

# Output F_plus sky map
print("F_PLUS_MAP")
for i in range(N_theta):
    row = " ".join(f"{F_plus[i, j]:.4f}" for j in range(N_phi))
    print(row)

print("F_CROSS_MAP")
for i in range(N_theta):
    row = " ".join(f"{F_cross[i, j]:.4f}" for j in range(N_phi))
    print(row)

print("F_COMBINED_MAP")
for i in range(N_theta):
    row = " ".join(f"{F_combined[i, j]:.4f}" for j in range(N_phi))
    print(row)

# Summary statistics
print()
print("SUMMARY_STATISTICS")
print(f"max |F_+| = {np.max(np.abs(F_plus)):.4f}")
print(f"max |F_x| = {np.max(np.abs(F_cross)):.4f}")
print(f"max |F_combined| = {np.max(F_combined):.4f}")

# Sky-averaged response
dOmega = np.sin(THETA) * (theta_arr[1] - theta_arr[0]) * (phi_arr[1] - phi_arr[0])
avg_F2 = np.sum((F_plus**2 + F_cross**2) * dOmega) / np.sum(dOmega)
print(f"<F_+^2 + F_x^2> = {avg_F2:.4f}  (analytic: 0.4000)")
print(f"RMS response = {np.sqrt(avg_F2):.4f}")

# Also compute for a few specific directions
directions = [
    ("Overhead (theta=0)", 0.0, 0.0),
    ("In-plane, along arm (theta=pi/2, phi=0)", np.pi/2, 0.0),
    ("In-plane, bisector (theta=pi/2, phi=pi/4)", np.pi/2, np.pi/4),
    ("45 deg elevation (theta=pi/4, phi=0)", np.pi/4, 0.0),
]

print()
print("SPECIFIC_DIRECTIONS (psi=0)")
for name, th, ph in directions:
    fp = (0.5 * (1 + np.cos(th)**2) * np.cos(2*ph) * cos_psi
          - np.cos(th) * np.sin(2*ph) * sin_psi)
    fc = (0.5 * (1 + np.cos(th)**2) * np.cos(2*ph) * sin_psi
          + np.cos(th) * np.sin(2*ph) * cos_psi)
    print(f"  {name}: F_+ = {fp:.4f}, F_x = {fc:.4f}, |F| = {np.sqrt(fp**2+fc**2):.4f}")

print()
print("Simulation complete.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Prev: Plane Waves & Polarizations Next: Quadrupole Formula →