8. Gravitational Wave Cosmology

Gravitational waves open an entirely new window on the cosmos — enabling direct measurements of cosmic distances without any reliance on the electromagnetic distance ladder, testing general relativity on cosmological scales, and probing the earliest instants of the universe through primordial gravitational wave backgrounds.

1. Introduction: GWs as a New Window on the Cosmos

Electromagnetic observations have built our understanding of cosmology through a chain of distance indicators: Cepheids calibrate Type Ia supernovae, which in turn reach cosmological distances. Each rung of this "cosmic distance ladder" introduces systematic uncertainties that compound. Gravitational waves offer a fundamentally different approach.

The key insight, first articulated by Bernard Schutz in 1986, is that the gravitational wave signal from a compact binary inspiral encodes the absolute luminosity distance to the source directly in the waveform amplitude. No calibration is needed — the distance measurement follows from general relativity alone. Such sources are called standard sirens, the gravitational analog of standard candles.

Why GW Cosmology Matters

● No distance ladder: luminosity distance is measured directly from the waveform
● Independent H₀ measurement: can arbitrate the Hubble tension between CMB and local probes
● Tests of gravity: GW propagation probes modified gravity on cosmological baselines
● Primordial universe: stochastic GW backgrounds encode physics from inflation to phase transitions
● Dark energy: multi-messenger events constrain the dark energy equation of state w(z)

The 2017 detection of the binary neutron star merger GW170817, accompanied by electromagnetic counterparts across the spectrum, demonstrated the power of this approach: a single event provided the first standard siren measurement of the Hubble constant. With next-generation detectors, hundreds to thousands of such measurements will transform precision cosmology.

2. Derivation: Standard Sirens

The gravitational wave strain from a compact binary inspiral in the quadrupole approximation takes the form:

$h(t) = \frac{4}{d_L} \left(\frac{G \mathcal{M}_c}{c^2}\right)^{5/3} \left(\frac{\pi f(t)}{c}\right)^{2/3} \cos\Phi(t)$

where $\mathcal{M}_c = (m_1 m_2)^{3/5}/(m_1+m_2)^{1/5}$ is the chirp mass,$f(t)$ is the instantaneous GW frequency, $d_L$ is the luminosity distance, and $\Phi(t)$ is the orbital phase. The crucial point is that the chirp mass is independently determined from the frequency evolution:

$\dot{f} = \frac{96}{5}\pi^{8/3} \left(\frac{G\mathcal{M}_c}{c^3}\right)^{5/3} f^{11/3}$

By measuring $f$ and $\dot{f}$ from the data, we extract $\mathcal{M}_c$. With $\mathcal{M}_c$ known, the amplitude of the waveform gives $d_L$ directly:

Distance Extraction

$d_L = \frac{4}{h_0} \left(\frac{G\mathcal{M}_c}{c^2}\right)^{5/3} \left(\frac{\pi f}{c}\right)^{2/3}$

where $h_0$ is the measured strain amplitude. This is a direct, absolute distance measurement — no calibration against other distance indicators is required.

From Distance to Hubble Constant

If an electromagnetic counterpart identifies the host galaxy and provides a redshift $z$, then in the low-redshift limit the Hubble law gives:

$H_0 = \frac{cz}{d_L}$

More generally, for sources at cosmological distances, the luminosity distance depends on the full cosmological model:

$d_L(z) = \frac{c(1+z)}{H_0} \int_0^z \frac{dz'}{E(z')}$

where $E(z) = \sqrt{\Omega_m(1+z)^3 + \Omega_\Lambda}$ for a flat $\Lambda$CDM universe.

Precision Scaling

For a single event, the fractional uncertainty on $d_L$ is dominated by the signal-to-noise ratio (SNR) and the degeneracy between distance and inclination angle $\iota$:

$\frac{\sigma(d_L)}{d_L} \sim \frac{1}{\text{SNR}} \times \text{geometric factor}$

The peculiar velocity of the host galaxy contributes an additional uncertainty$\sigma_v \sim 200\text{ km/s}$, which dominates for nearby events. For $N$independent standard siren measurements, the combined uncertainty on $H_0$ scales as:

$\frac{\sigma(H_0)}{H_0} \sim \frac{1}{\sqrt{N}} \sqrt{\left\langle\left(\frac{\sigma(d_L)}{d_L}\right)^2 + \left(\frac{\sigma_v}{cz}\right)^2\right\rangle}$

For $N \sim 50$ BNS events with SNR $\sim 20$ at $z \sim 0.05$, this yields $\sigma(H_0)/H_0 \sim 2\%$, competitive with current best measurements.

GW170817: The First Standard Siren

The binary neutron star merger GW170817 (August 17, 2017) was accompanied by a short gamma-ray burst GRB 170817A and an optical/infrared kilonova in the galaxy NGC 4993 at redshift $z = 0.0098$. The LIGO-Virgo measurement yielded:

$H_0 = 70.0^{+12.0}_{-8.0} \text{ km s}^{-1}\text{ Mpc}^{-1}$

Abbott et al. (2017), Nature 551, 85

This single-event measurement is consistent with both the CMB value ($67.4 \pm 0.5$ km/s/Mpc) and the local distance ladder value ($73.0 \pm 1.0$ km/s/Mpc). While the uncertainty is large from one event, it demonstrates the method works and is free of systematic errors that plague other approaches. The $\sim 14\%$ uncertainty is dominated by the distance-inclination degeneracy.

3. Derivation: Dark Sirens (Statistical Method)

Most gravitational wave events — particularly binary black hole (BBH) mergers — lack electromagnetic counterparts. Without a direct redshift measurement, the standard siren method seems inapplicable. However, the statistical dark siren method recovers cosmological information by cross-correlating GW sky localizations with galaxy catalogs.

The Bayesian Framework

For a single GW event with measured luminosity distance $\hat{d}_L$ and sky localization$\Delta\Omega$, the likelihood of the cosmological parameters (here just $H_0$) is obtained by marginalizing over all potential host galaxies:

$p(\hat{d}_L | H_0) = \sum_{i \in \text{galaxies}} w_i \, p(\hat{d}_L | d_L(z_i, H_0))$

where $w_i$ is the weight of galaxy $i$ (based on luminosity, mass, star formation rate), and $d_L(z_i, H_0)$ is the luminosity distance at the galaxy's redshift for a given $H_0$.

The individual galaxy likelihood is typically a Gaussian centered on the GW-measured distance:

$p(\hat{d}_L | d_L(z_i, H_0)) = \frac{1}{\sqrt{2\pi}\sigma_{d_L}} \exp\left[-\frac{(\hat{d}_L - d_L(z_i, H_0))^2}{2\sigma_{d_L}^2}\right]$

Combining Multiple Events

For $N$ independent GW events, the posterior on $H_0$ is obtained by multiplying the individual likelihoods:

Combined Dark Siren Posterior

$p(H_0 | \{\text{data}\}) \propto \pi(H_0) \prod_{k=1}^{N} \left[\sum_{i \in \text{gal}_k} w_i \, p(\hat{d}_{L,k} | d_L(z_i, H_0))\right]$

where $\pi(H_0)$ is the prior on the Hubble constant and the product runs over all events.

Each event contributes a broad, multi-peaked likelihood (one peak per candidate galaxy). As events accumulate, the true value of $H_0$ is consistently supported while false peaks are suppressed. The effective scaling is:

$\frac{\sigma(H_0)}{H_0} \sim \frac{1}{\sqrt{N_{\text{eff}}}}$

where $N_{\text{eff}} \leq N$ accounts for the dilution from multiple candidate hosts per event.

Projected Constraints

●O3 (LIGO/Virgo): ~50 BBH events yield $\sigma(H_0)/H_0 \sim 10\text{-}20\%$
●O5 (design sensitivity): ~100 well-localized BBH events yield $\sigma(H_0)/H_0 \sim 2\text{-}5\%$
●Next-generation (ET, CE): thousands of events with sub-degree localization yield sub-percent precision

The dark siren method has already been applied to the BBH event GW190814 and to the full GWTC-3 catalog, yielding $H_0 = 68^{+8}_{-6}$ km/s/Mpc from ~47 events — demonstrating the method's viability even before reaching percent-level precision.

4. Derivation: Modified Gravity and GW Propagation

Gravitational waves that travel cosmological distances serve as powerful probes of deviations from general relativity. In GR, GWs propagate at exactly the speed of light, are non-dispersive, carry only two tensor polarizations, and suffer no anomalous damping. Modified gravity theories can violate any or all of these properties.

(a) Anomalous Speed: $c_T \neq c$

In general scalar-tensor and Horndeski theories, the GW propagation speed can differ from the speed of light. The arrival time difference between GWs and photons from the same event constrains:

$\frac{\Delta t}{t_{\text{travel}}} = \frac{c_T - c}{c} \equiv \alpha_T$

GW170817 and GRB 170817A arrived within $\Delta t = 1.7$ seconds of each other after traveling$\sim 40$ Mpc ($\sim 1.3 \times 10^8$ years), yielding the extraordinary constraint:

$-3 \times 10^{-15} \leq \alpha_T \leq +7 \times 10^{-16}$

This single measurement eliminated large classes of modified gravity theories in one stroke.

(b) Massive Graviton: Dispersion

If the graviton has a non-zero mass $m_g$, gravitational waves obey a modified dispersion relation:

$E^2 = p^2 c^2 + m_g^2 c^4$

$v_g = c\sqrt{1 - \frac{m_g^2 c^4}{E^2}} \approx c\left(1 - \frac{1}{2}\frac{m_g^2 c^4}{\hbar^2 \omega^2}\right)$

This frequency-dependent speed causes dispersion: lower-frequency components arrive later than higher-frequency components. The accumulated phase shift over a distance $D$ is:

$\delta\Psi = -\frac{\pi^2 D \mathcal{M}_c}{\lambda_g^2 (1+z)} \frac{1}{(\pi \mathcal{M}_c f)^{1/3}}$

where $\lambda_g = h/(m_g c)$ is the graviton Compton wavelength.

Current LIGO-Virgo constraints give $\lambda_g > 1.6 \times 10^{13}$ km, corresponding to $m_g < 7.7 \times 10^{-23}$ eV/$c^2$.

(c) Extra Polarizations

In GR, gravitational waves have exactly two tensor polarizations ($h_+$ and $h_\times$). Alternative theories of gravity can predict up to six polarization modes:

●Tensor modes ($h_+, h_\times$): present in GR and all metric theories
●Scalar breathing mode ($h_b$): isotropic expansion/contraction; appears in Brans-Dicke theory
●Scalar longitudinal mode ($h_L$): deformation along propagation direction; massive gravity
●Vector modes ($h_x, h_y$): appear in Einstein-Aether and TeVeS theories

Detecting or constraining extra polarizations requires a network of non-co-aligned detectors. With LIGO Hanford, LIGO Livingston, Virgo, KAGRA, and LIGO-India, the network can distinguish tensor from non-tensor modes for sufficiently loud signals.

(d) Amplitude Damping: GW Friction

In many modified gravity theories, GW propagation on a cosmological (FRW) background acquires an additional friction term beyond the standard Hubble damping. The modified propagation equation reads:

$\ddot{h}_A + (3 + \alpha_M)H\dot{h}_A + c_T^2 \frac{k^2}{a^2}h_A = 0$

where $\alpha_M = \frac{d\ln M_*^2}{d\ln a}$ is the Planck-mass run rate and $A = +, \times$ labels the polarization.

The extra friction $\alpha_M H \dot{h}$ modifies the GW amplitude relative to GR. This means the luminosity distance inferred from GWs ($d_L^{\text{GW}}$) differs from the electromagnetic luminosity distance ($d_L^{\text{EM}}$):

$d_L^{\text{GW}}(z) = d_L^{\text{EM}}(z) \exp\left[\frac{1}{2}\int_0^z \frac{\alpha_M(z')}{1+z'} dz'\right]$

A discrepancy between GW and EM distance measurements to the same source would be a smoking gun for modified gravity. Current data constrain $|\alpha_M| \lesssim \mathcal{O}(1)$.

5. Derivation: Primordial Gravitational Wave Background

Inflation generically predicts a stochastic background of gravitational waves produced by quantum vacuum fluctuations of the tensor metric perturbations. This primordial GW background carries information about the energy scale of inflation and subsequent cosmological phase transitions.

Inflationary Tensor Spectrum

During inflation, tensor perturbations $h_{ij}$ are amplified from quantum fluctuations. The primordial tensor power spectrum at horizon crossing is:

$\mathcal{P}_T(k) = \frac{2}{\pi^2}\frac{H_{\text{inf}}^2}{M_{\text{Pl}}^2} \left(\frac{k}{k_*}\right)^{n_T}$

where $n_T = -r/8$ is the tensor spectral index (consistency relation) and $r = \mathcal{P}_T/\mathcal{P}_S$ is the tensor-to-scalar ratio.

Energy Density Spectrum Today

The present-day energy density of the primordial GW background, normalized to the critical density, is:

Primordial GW Energy Density

$\Omega_{\text{GW}}(f) = \frac{1}{\rho_c}\frac{d\rho_{\text{GW}}}{d\ln f} = \Omega_{\text{rad}} \times \frac{r}{24} \times T^2(f) \times \left(\frac{f}{f_*}\right)^{n_T}$

where $\Omega_{\text{rad}} \approx 9 \times 10^{-5}$ is the radiation density parameter and $T(f)$ is the transfer function encoding the suppression of modes that entered the horizon during the matter-dominated era.

The transfer function has two key regimes:

$T^2(f) \approx \begin{cases} 1 & f \gg f_{\text{eq}} \approx 10^{-16}\text{ Hz (radiation era entry)} \\ (f_{\text{eq}}/f)^2 & f \ll f_{\text{eq}}\text{ (matter era entry)} \end{cases}$

For the standard inflationary prediction with $r = 0.01$ (near current upper limits):

$\Omega_{\text{GW}}^{\text{inf}} \sim \Omega_{\text{rad}} \times \frac{r}{24} \sim 9 \times 10^{-5} \times 4 \times 10^{-4} \sim 4 \times 10^{-8} \times T^2(f)$

In the LIGO/LISA band ($f \gg f_{\text{eq}}$), $T^2 \approx 1$ but the spectrum is further reduced by the number of relativistic species, giving $\Omega_{\text{GW}} \sim 10^{-16}$for $r = 0.01$ — far below current detector sensitivity.

Phase Transitions: Electroweak and QCD

First-order cosmological phase transitions produce GW backgrounds through three mechanisms: bubble collisions, sound waves in the plasma, and turbulence. The peak frequency and amplitude depend on the transition temperature $T_*$ and strength $\alpha$:

$f_{\text{peak}} \sim 10^{-3}\text{ Hz} \times \left(\frac{T_*}{100\text{ GeV}}\right) \times \left(\frac{\beta}{H_*}\right)$

$\Omega_{\text{GW}}^{\text{PT}} \sim 10^{-6} \times \left(\frac{\alpha}{1+\alpha}\right)^2 \times \left(\frac{H_*}{\beta}\right)^2$

where $\beta$ is the inverse duration of the transition and $H_*$ is the Hubble rate at the transition.

Key Phase Transitions

●Electroweak ($T_* \sim 100$ GeV): peak at $f \sim 10^{-3}$ Hz (LISA band). Only produces GWs if the transition is first-order (requires BSM physics).
●QCD ($T_* \sim 150$ MeV): peak at $f \sim 10^{-8}$ Hz (PTA band). Standard Model predicts a crossover, but BSM scenarios allow first-order transition.

Cosmic Strings

Cosmic strings — topological defects from symmetry-breaking phase transitions — produce a nearly scale-invariant GW background:

$\Omega_{\text{GW}}^{\text{CS}}(f) \sim 10^{-8} \times \left(\frac{G\mu}{10^{-11}}\right)$

where $G\mu$ is the string tension parameter. Current PTA data constrain $G\mu \lesssim 10^{-11}$.

6. Applications

Resolving the Hubble Tension

The $\sim 5\sigma$ discrepancy between the CMB-inferred value ($H_0 \approx 67.4$ km/s/Mpc) and the local distance-ladder value ($H_0 \approx 73.0$ km/s/Mpc) is one of the most pressing problems in cosmology. Standard sirens offer a completely independent measurement with different systematics. With $\sim 50$ bright sirens or $\sim 200$ dark sirens, GW observations can definitively distinguish between the two values.

Dark Energy from GW+EM

Multi-messenger standard siren observations at higher redshifts ($z \sim 0.1\text{-}2$) constrain the dark energy equation of state parameter $w(z)$. By measuring $d_L(z)$at multiple redshifts, one can reconstruct the expansion history $H(z)$ and test whether dark energy is truly a cosmological constant ($w = -1$) or evolves with time.

Testing Lorentz Invariance

The coincident detection of GWs and gamma rays from GW170817 constrains violations of Lorentz invariance in the gravitational sector. Combined with the non-observation of graviton-mass-induced dispersion in the waveform, these measurements constrain Lorentz-violating terms in the graviton dispersion relation at energy scales up to $\sim 10^{16}$ GeV.

Probing the Early Universe

Unlike electromagnetic radiation, which decouples at the CMB ($z \sim 1100$), gravitational waves propagate freely from the Planck epoch. The primordial GW spectrum encodes the energy scale of inflation ($H_{\text{inf}}$), the number of relativistic species $g_*(T)$, and any exotic phase transitions in the early universe.

7. Historical Context

1986

Schutz Proposes Standard Sirens

Bernard Schutz publishes "Determining the Hubble constant from gravitational wave observations" in Nature, showing that compact binary inspirals are self-calibrating distance indicators. This foundational paper establishes the entire field of GW cosmology.

2005

Dalal et al.: Dark Siren Statistics

The statistical galaxy catalog method for measuring $H_0$ without electromagnetic counterparts is developed, extending the standard siren concept to binary black hole mergers.

2012

Nissanke et al.: Realistic Forecasts

Detailed Fisher-matrix and Bayesian forecasts show that a network of advanced detectors observing$\sim 50$ BNS events can measure $H_0$ to $\sim 2\%$.

2017

GW170817: First Standard Siren Measurement

The binary neutron star merger GW170817, with its electromagnetic counterpart, delivers the first standard siren measurement of $H_0 = 70^{+12}_{-8}$ km/s/Mpc. The same event constrains$|c_T/c - 1| < 10^{-15}$, ruling out vast swaths of modified gravity parameter space.

2021

GWTC-3 Dark Siren Analysis

The full third gravitational wave transient catalog enables a dark siren measurement of $H_0$using ~47 BBH events combined with galaxy catalog cross-correlation, achieving $\sim 12\%$ precision.

2023

NANOGrav Stochastic Background

NANOGrav, EPTA, PPTA, and CPTA report evidence for a stochastic gravitational wave background in pulsar timing data at nanohertz frequencies. The signal is consistent with supermassive black hole binaries but also admits interpretations involving primordial sources (cosmic strings, phase transitions).

8. Python Simulation: Standard Sirens & Primordial GW Spectrum

This simulation demonstrates two key aspects of GW cosmology: (1) measuring $H_0$ from simulated standard siren events using both bright and dark siren methods, and (2) computing the primordial GW energy density spectrum for different inflationary and phase transition models.

Python

script.py233 lines

import numpy as np

# ==============================================================
# GW Cosmology Simulation
# Part 1: H0 posterior from standard sirens (bright + dark)
# Part 2: Primordial GW background spectrum
# ==============================================================

np.random.seed(42)

# --- Physical constants and fiducial cosmology ---
c_km = 3.0e5          # speed of light in km/s
H0_true = 70.0        # true H0 in km/s/Mpc
sigma_v = 200.0       # peculiar velocity dispersion km/s

# ==============================================================
# PART 1: H0 from Standard Sirens
# ==============================================================

# --- Simulate bright siren events (BNS with EM counterpart) ---
N_bright = 15
z_bright = np.random.uniform(0.01, 0.08, N_bright)
dL_true_bright = c_km * z_bright / H0_true  # Mpc (low-z approx)
sigma_dL_frac = 0.15  # 15% distance uncertainty per event
sigma_dL_bright = sigma_dL_frac * dL_true_bright
dL_obs_bright = dL_true_bright + np.random.normal(0, sigma_dL_bright)
# Redshift uncertainty from peculiar velocity
sigma_z_bright = sigma_v / c_km

# --- Simulate dark siren events (BBH, no direct redshift) ---
N_dark = 80
z_dark_true = np.random.uniform(0.02, 0.15, N_dark)
dL_true_dark = c_km * z_dark_true / H0_true
sigma_dL_dark = 0.25 * dL_true_dark  # 25% distance uncertainty (poorer localization)
dL_obs_dark = dL_true_dark + np.random.normal(0, sigma_dL_dark)

# For each dark event, generate candidate galaxies in the error volume
N_gal_per_event = 20
z_gal_catalog = []
for i in range(N_dark):
    # True host + random field galaxies
    z_cands = np.random.uniform(
        max(0.005, z_dark_true[i] - 0.04),
        z_dark_true[i] + 0.04,
        N_gal_per_event - 1
    )
    z_cands = np.append(z_cands, z_dark_true[i] + np.random.normal(0, 0.002))
    np.random.shuffle(z_cands)
    z_gal_catalog.append(z_cands)

# --- Compute H0 posterior ---
H0_grid = np.linspace(50.0, 100.0, 500)
log_posterior = np.zeros_like(H0_grid)

# Bright sirens contribution
for i in range(N_bright):
    z_i = z_bright[i]
    sigma_z = sigma_z_bright
    for j, H0 in enumerate(H0_grid):
        dL_model = c_km * z_i / H0
        sigma_total = np.sqrt(sigma_dL_bright[i]**2 + (c_km * sigma_z / H0)**2)
        log_posterior[j] += -0.5 * ((dL_obs_bright[i] - dL_model) / sigma_total)**2 - np.log(sigma_total)

# Dark sirens contribution (marginalize over candidate galaxies)
for i in range(N_dark):
    for j, H0 in enumerate(H0_grid):
        likelihoods = np.zeros(N_gal_per_event)
        for k, z_g in enumerate(z_gal_catalog[i]):
            if z_g <= 0:
                continue
            dL_model = c_km * z_g / H0
            sigma_total = sigma_dL_dark[i]
            likelihoods[k] = np.exp(-0.5 * ((dL_obs_dark[i] - dL_model) / sigma_total)**2) / sigma_total
        marg_like = np.mean(likelihoods)
        if marg_like > 0:
            log_posterior[j] += np.log(marg_like)

# Normalize
log_posterior -= np.max(log_posterior)
posterior = np.exp(log_posterior)
posterior /= np.trapz(posterior, H0_grid)

# Also compute bright-only posterior
log_post_bright = np.zeros_like(H0_grid)
for i in range(N_bright):
    z_i = z_bright[i]
    sigma_z = sigma_z_bright
    for j, H0 in enumerate(H0_grid):
        dL_model = c_km * z_i / H0
        sigma_total = np.sqrt(sigma_dL_bright[i]**2 + (c_km * sigma_z / H0)**2)
        log_post_bright[j] += -0.5 * ((dL_obs_bright[i] - dL_model) / sigma_total)**2 - np.log(sigma_total)

log_post_bright -= np.max(log_post_bright)
post_bright = np.exp(log_post_bright)
post_bright /= np.trapz(post_bright, H0_grid)

# Summary statistics
H0_median = H0_grid[np.searchsorted(np.cumsum(posterior) / np.sum(posterior), 0.5)]
cumul = np.cumsum(posterior) / np.sum(posterior)
H0_lo = H0_grid[np.searchsorted(cumul, 0.16)]
H0_hi = H0_grid[np.searchsorted(cumul, 0.84)]

H0_med_bright = H0_grid[np.searchsorted(np.cumsum(post_bright) / np.sum(post_bright), 0.5)]
cumul_b = np.cumsum(post_bright) / np.sum(post_bright)
H0_lo_b = H0_grid[np.searchsorted(cumul_b, 0.16)]
H0_hi_b = H0_grid[np.searchsorted(cumul_b, 0.84)]

print("=== H0 Posterior from Standard Sirens ===")
print(f"True H0 = {H0_true:.1f} km/s/Mpc")
print(f"Bright only ({N_bright} BNS): H0 = {H0_med_bright:.1f} +{H0_hi_b - H0_med_bright:.1f} / -{H0_med_bright - H0_lo_b:.1f}")
print(f"Combined ({N_bright} bright + {N_dark} dark): H0 = {H0_median:.1f} +{H0_hi - H0_median:.1f} / -{H0_median - H0_lo:.1f}")
print()

# Print posterior for plotting
print("H0_posterior_plot")
print("H0(km/s/Mpc)  Bright_Only  Combined  CMB_value  SHoES_value")
step = max(1, len(H0_grid) // 120)
for i in range(0, len(H0_grid), step):
    H0_val = H0_grid[i]
    # Gaussian peaks for CMB and SHoES reference
    cmb = np.exp(-0.5 * ((H0_val - 67.4) / 0.5)**2)
    shoes = np.exp(-0.5 * ((H0_val - 73.0) / 1.0)**2)
    # Scale reference curves
    cmb_scaled = cmb * np.max(posterior) * 0.6
    shoes_scaled = shoes * np.max(posterior) * 0.6
    print(f"{H0_val:.2f}  {post_bright[i]:.6f}  {posterior[i]:.6f}  {cmb_scaled:.6f}  {shoes_scaled:.6f}")

# ==============================================================
# PART 2: Primordial GW Background Spectrum
# ==============================================================
print()
print("=== Primordial GW Background Spectrum ===")

# Frequency range (Hz)
f = np.logspace(-18, 4, 300)

# Physical parameters
Omega_rad = 9.0e-5
f_eq = 2.0e-17        # matter-radiation equality frequency (Hz)
g_star_0 = 3.36       # effective relativistic dof today
g_star_inf = 106.75    # dof during inflation

# Transfer function
def transfer_sq(freq, f_equality):
    T2 = np.ones_like(freq)
    mask = freq < f_equality
    T2[mask] = (freq[mask] / f_equality)**2
    # Additional suppression from change in g_star
    g_factor = (g_star_0 / g_star_inf)**(1.0/3.0)
    return T2 * g_factor

# --- Model 1: Inflation with r = 0.06 ---
r1 = 0.06
nT1 = -r1 / 8.0
f_star = 7.7e-17  # pivot scale in Hz
Omega_inf_1 = Omega_rad * (r1 / 24.0) * transfer_sq(f, f_eq) * (f / f_star)**nT1

# --- Model 2: Inflation with r = 0.01 ---
r2 = 0.01
nT2 = -r2 / 8.0
Omega_inf_2 = Omega_rad * (r2 / 24.0) * transfer_sq(f, f_eq) * (f / f_star)**nT2

# --- Model 3: Inflation with r = 0.001 ---
r3 = 0.001
nT3 = -r3 / 8.0
Omega_inf_3 = Omega_rad * (r3 / 24.0) * transfer_sq(f, f_eq) * (f / f_star)**nT3

# --- Electroweak phase transition (peaked at ~mHz) ---
f_ew = 3.0e-3       # peak frequency Hz
alpha_ew = 0.5       # transition strength
beta_H_ew = 10.0     # beta/H*
Omega_ew_peak = 1.0e-6 * (alpha_ew / (1 + alpha_ew))**2 * (1.0 / beta_H_ew)**2
# Broken power law shape
Omega_ew = Omega_ew_peak * (f / f_ew)**3 / (1 + (f / f_ew)**2)**2.5

# --- QCD phase transition (peaked at ~nHz) ---
f_qcd = 3.0e-8      # peak frequency Hz
alpha_qcd = 0.3
beta_H_qcd = 20.0
Omega_qcd_peak = 1.0e-6 * (alpha_qcd / (1 + alpha_qcd))**2 * (1.0 / beta_H_qcd)**2
Omega_qcd = Omega_qcd_peak * (f / f_qcd)**3 / (1 + (f / f_qcd)**2)**2.5

# --- Cosmic strings ---
Gmu = 1.0e-11
Omega_cs = 1.0e-8 * (Gmu / 1.0e-11) * np.ones_like(f)
# Suppress at very low and very high frequencies
Omega_cs *= 1.0 / (1.0 + (1.0e-15 / f)**2) * 1.0 / (1.0 + (f / 1.0e3)**2)

# --- Detector sensitivity curves (approximate) ---
# LIGO design sensitivity
f_ligo = np.logspace(0, 3.5, 100)
Omega_ligo = 1.0e-9 * ((f_ligo / 25.0)**(-4) + 1.0 + (f_ligo / 300.0)**2)

# LISA sensitivity
f_lisa = np.logspace(-4, -0.5, 100)
Omega_lisa = 1.0e-12 * ((f_lisa / 3.0e-3)**(-2) + 1.0 + (f_lisa / 0.1)**2)

# PTA sensitivity
f_pta = np.logspace(-9, -6.5, 50)
Omega_pta = 1.0e-9 * (f_pta / 1.0e-8)**(-4.0/3.0)

# NANOGrav signal region
f_nano = np.logspace(-9, -7.5, 30)
Omega_nano = 6.0e-10 * (f_nano / 3.17e-8)**(-2.0/3.0)

print("GW_spectrum_plot")
print("f(Hz)  Inflation_r=0.06  Inflation_r=0.01  Inflation_r=0.001  EW_PhaseTransition  QCD_PhaseTransition  CosmicStrings")
step2 = max(1, len(f) // 150)
for i in range(0, len(f), step2):
    print(f"{f[i]:.4e}  {Omega_inf_1[i]:.4e}  {Omega_inf_2[i]:.4e}  {Omega_inf_3[i]:.4e}  {Omega_ew[i]:.4e}  {Omega_qcd[i]:.4e}  {Omega_cs[i]:.4e}")

# Print detector curves separately
print()
print("Detector_curves")
print("f(Hz)  Sensitivity")
print("# LIGO")
for i in range(0, len(f_ligo), 5):
    print(f"{f_ligo[i]:.4e}  {Omega_ligo[i]:.4e}")
print("# LISA")
for i in range(0, len(f_lisa), 5):
    print(f"{f_lisa[i]:.4e}  {Omega_lisa[i]:.4e}")
print("# PTA")
for i in range(0, len(f_pta), 3):
    print(f"{f_pta[i]:.4e}  {Omega_pta[i]:.4e}")

print()
print("=== Key Results ===")
print(f"Inflation r=0.06: Omega_GW(f=100Hz) ~ {Omega_inf_1[np.argmin(np.abs(f - 100))]:.2e}")
print(f"Inflation r=0.01: Omega_GW(f=100Hz) ~ {Omega_inf_2[np.argmin(np.abs(f - 100))]:.2e}")
print(f"EW phase transition peak: Omega_GW ~ {Omega_ew_peak:.2e} at f ~ {f_ew:.0e} Hz")
print(f"QCD phase transition peak: Omega_GW ~ {Omega_qcd_peak:.2e} at f ~ {f_qcd:.0e} Hz")
print(f"Cosmic strings (Gmu=1e-11): Omega_GW ~ {Omega_cs[len(f)//2]:.2e}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Key Takeaways

Distance Measurements

●GW inspiral waveforms encode $d_L$ directly — no distance ladder needed
●Bright sirens (BNS+EM) give $H_0$ from single events; precision scales as $1/\sqrt{N}$
●Dark sirens (BBH+galaxy catalogs) use Bayesian marginalization over candidate hosts
●GW170817 gave $H_0 = 70 \pm 12$ km/s/Mpc from a single event

Fundamental Physics

●GW speed constrained to $|c_T/c - 1| < 10^{-15}$
●Graviton mass bounded: $m_g < 7.7 \times 10^{-23}$ eV/$c^2$
●GW friction ($\alpha_M$) probes modified gravity via $d_L^{\text{GW}} \neq d_L^{\text{EM}}$
●Primordial backgrounds encode inflation ($r$), phase transitions, and cosmic strings

←7. Major Detections 9. Future Prospects→