Introduction

Observational cosmology provides the empirical foundation for our understanding of the universe. Through careful measurements of cosmic phenomena, we determine the fundamental parameters of cosmology and test theoretical predictions.

Modern cosmology is remarkably well-described by the ΛCDM model with just six parameters, determined to unprecedented precision by multiple independent observations.

1. Cosmic Microwave Background (CMB)

Discovery and Temperature

The CMB is relic radiation from the Big Bang, discovered by Penzias and Wilson in 1965. It has a nearly perfect blackbody spectrum at temperature:

$$T_{\text{CMB}} = 2.7255 \pm 0.0006 \text{ K}$$

Temperature Anisotropies

The CMB temperature varies across the sky at the level of $\Delta T/T \sim 10^{-5}$. These anisotropies are expanded in spherical harmonics:

$$\frac{\Delta T}{T}(\theta, \phi) = \sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell} a_{\ell m} Y_{\ell m}(\theta, \phi)$$

Angular Power Spectrum

The angular power spectrum characterizes the temperature fluctuations:

$$C_\ell = \frac{1}{2\ell+1}\sum_{m=-\ell}^{\ell}|a_{\ell m}|^2$$

Key features of the power spectrum:

• Sachs-Wolfe plateau ($\ell < 100$): Gravitational potential fluctuations
• Acoustic peaks ($\ell \sim 200-1000$): Sound waves in photon-baryon fluid
• Damping tail ($\ell > 1000$): Silk damping from diffusion

CMB Angular Power Spectrum Simulation

Python

Analytic model of the CMB power spectrum showing acoustic peaks, Sachs-Wolfe plateau, and Silk damping.

script.py60 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Simplified CMB angular power spectrum model
# Based on analytic approximations to the acoustic peaks

ell = np.arange(2, 2500)
# Sachs-Wolfe plateau + acoustic peaks model
# Parameters
A_s = 2.1e-9   # scalar amplitude
n_s = 0.965     # spectral index
theta_s = 0.0104  # sound horizon angle (radians)
ell_eq = 150    # matter-radiation equality

# Transfer function approximation
def cmb_power_spectrum(ell, A_s, n_s, theta_s):
    # Primordial spectrum
    P_prim = A_s * (ell / 200.0) ** (n_s - 1)
    # Acoustic oscillations
    k_s = ell * theta_s
    acoustic = np.cos(k_s) ** 2
    # Silk damping
    ell_silk = 1500
    damping = np.exp(-(ell / ell_silk) ** 1.4)
    # Combine
    Cl = P_prim * (1.0 + 5.0 * acoustic) * damping
    Dl = ell * (ell + 1) * Cl / (2 * np.pi) * 1e12
    return Dl

Dl = cmb_power_spectrum(ell, A_s, n_s, theta_s)

fig, ax = plt.subplots(figsize=(12, 6))
fig.patch.set_facecolor('#0a0a1a')
ax.set_facecolor('#0a0a1a')
for sp in ax.spines.values(): sp.set_edgecolor('#334155')
ax.tick_params(colors='#94a3b8')

ax.plot(ell, Dl, color='#60a5fa', lw=1.5, label='Model (analytic approx.)')
ax.set_xlabel('Multipole moment l', color='#94a3b8', fontsize=12)
ax.set_ylabel('D_l  [uK^2]', color='#94a3b8', fontsize=12)
ax.set_title('CMB Angular Power Spectrum', color='white', fontsize=14)
ax.set_xlim(2, 2500)
ax.set_xscale('log')
ax.legend(facecolor='#1e293b', labelcolor='white')
ax.grid(alpha=0.2, color='#334155')

# Mark acoustic peaks
peak_positions = [220, 540, 810, 1120, 1450]
ymax = np.max(Dl) * 0.95
for i, lp in enumerate(peak_positions):
    ax.axvline(lp, color='#f472b6', alpha=0.3, ls='--')
    ax.text(lp, ymax, f'Peak {i+1}', color='#f472b6', fontsize=8, ha='center', va='top')

plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')
plt.close()
print('CMB Power Spectrum computed')
print(f'First acoustic peak at l ~ {peak_positions[0]}')
print(f'Sound horizon angle: {theta_s*180/np.pi:.2f} degrees')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Physical Interpretation of Peaks

The first acoustic peak location determines the spatial curvature:

$$\ell_1 \approx \frac{\pi}{\theta_*} = \frac{\pi D_A(z_*)}{{r_s(z_*)}}$$

where $D_A$ is the angular diameter distance to recombination and $r_s$ is the sound horizon:

$$r_s(z) = \int_z^{\infty} \frac{c_s(z')}{H(z')} dz'$$

The sound speed in the photon-baryon plasma:

$$c_s = \frac{c}{\sqrt{3(1 + R)}}$$

where $R = 3\rho_b/4\rho_\gamma$ is the baryon-to-photon density ratio.

Polarization

Thomson scattering produces linear polarization, decomposed into E-modes and B-modes:

E-modes: Gradient-type patterns from density perturbations
B-modes: Curl-type patterns from gravitational waves (primordial) or lensing (secondary)

2. Cosmic Distance Ladder

Standard Candles

Objects with known intrinsic luminosity $L$ allow distance determination via:

$$d_L = \sqrt{\frac{L}{4\pi F}}$$

where $F$ is the observed flux and $d_L$ is the luminosity distance.

Type Ia Supernovae

SNe Ia result from white dwarf explosions near the Chandrasekhar limit, providing standardizable candles. After empirical corrections (stretch and color):

$$m_B = M_B + 5\log_{10}(d_L/10\text{ pc}) + \alpha X - \beta C$$

where $X$ is the light curve stretch factor and $C$ is the color excess.

Key discovery (1998): Distant SNe Ia are dimmer than expected in a decelerating universe, revealing cosmic acceleration!

Cepheid Variables

Period-luminosity relation:

$$M_V = -2.43(\log_{10} P - 1) - 4.05$$

where $P$ is the pulsation period in days.

Cosmic Distance Measures and Hubble Parameter

Python

Cosmological distance measures (comoving, luminosity, angular diameter) and Hubble parameter evolution for different universe models.

script.py65 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
def quad(f, a, b, limit=200, **kw):
    x = np.linspace(a if a != 0 else 1e-12, b, max(limit, 1000))
    y = np.array([f(xi) for xi in x]); return np.trapezoid(y, x), 0.0

c = 299792.458  # km/s
H0 = 70.0       # km/s/Mpc

def E(z, Om=0.3, OL=0.7):
    return np.sqrt(Om*(1+z)**3 + OL)

def comoving_distance(z, Om=0.3, OL=0.7):
    result, _ = quad(lambda zp: 1.0/E(zp, Om, OL), 0, z)
    return (c/H0) * result

def luminosity_distance(z, Om=0.3, OL=0.7):
    return (1+z) * comoving_distance(z, Om, OL)

def angular_diameter_distance(z, Om=0.3, OL=0.7):
    return comoving_distance(z, Om, OL) / (1+z)

z_arr = np.linspace(0.01, 3.0, 300)
d_C = np.array([comoving_distance(z) for z in z_arr])
d_L = np.array([luminosity_distance(z) for z in z_arr])
d_A = np.array([angular_diameter_distance(z) for z in z_arr])

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
fig.patch.set_facecolor('#0a0a1a')
for ax in [ax1, ax2]:
    ax.set_facecolor('#0a0a1a')
    for sp in ax.spines.values(): sp.set_edgecolor('#334155')
    ax.tick_params(colors='#94a3b8')

ax1.plot(z_arr, d_C, color='#60a5fa', lw=2, label='Comoving distance')
ax1.plot(z_arr, d_L, color='#f472b6', lw=2, label='Luminosity distance')
ax1.plot(z_arr, d_A, color='#34d399', lw=2, label='Angular diameter distance')
ax1.set_xlabel('Redshift z', color='#94a3b8')
ax1.set_ylabel('Distance [Mpc]', color='#94a3b8')
ax1.set_title('Cosmological Distances (LCDM)', color='white')
ax1.legend(facecolor='#1e293b', labelcolor='white')
ax1.grid(alpha=0.2, color='#334155')

# Hubble parameter vs redshift for different models
z_h = np.linspace(0, 3, 200)
H_lcdm = H0 * np.array([E(z, 0.3, 0.7) for z in z_h])
H_eds = H0 * np.array([E(z, 1.0, 0.0) for z in z_h])
H_deS = H0 * np.array([E(z, 0.0, 1.0) for z in z_h])

ax2.plot(z_h, H_lcdm, color='#60a5fa', lw=2, label='LCDM (0.3, 0.7)')
ax2.plot(z_h, H_eds, color='#fbbf24', lw=2, label='EdS (1.0, 0.0)')
ax2.plot(z_h, H_deS, color='#a78bfa', lw=2, label='de Sitter (0.0, 1.0)')
ax2.set_xlabel('Redshift z', color='#94a3b8')
ax2.set_ylabel('H(z) [km/s/Mpc]', color='#94a3b8')
ax2.set_title('Hubble Parameter Evolution', color='white')
ax2.legend(facecolor='#1e293b', labelcolor='white')
ax2.grid(alpha=0.2, color='#334155')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')
plt.close()
print(f'At z=1: d_L = {luminosity_distance(1.0):.1f} Mpc, d_A = {angular_diameter_distance(1.0):.1f} Mpc')
print(f'At z=1100 (CMB): d_A = {angular_diameter_distance(1100):.1f} Mpc')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

3. Baryon Acoustic Oscillations (BAO)

Physical Origin

Sound waves in the pre-recombination plasma create a characteristic scale in the matter distribution. The sound horizon at drag epoch provides a "standard ruler":

$$r_d = \int_0^{z_d} \frac{c_s(z)}{H(z)} dz \approx 147 \text{ Mpc}$$

Two-Point Correlation Function

The galaxy correlation function $\xi(r)$ shows a bump at $r \approx r_d$:

$$\xi(r) = \frac{1}{2\pi^2}\int_0^{\infty} P(k) \frac{\sin(kr)}{kr} k^2 dk$$

where $P(k)$ is the matter power spectrum.

Cosmological Constraints

BAO measurements constrain the combination:

$$D_V(z) = \left[(1+z)^2 D_A^2(z) \frac{cz}{H(z)}\right]^{1/3}$$

This provides geometric constraints complementary to SNe Ia and CMB.

BAO Scale Calculator (Fortran)

Fortran

Computes BAO-related distance measures at multiple redshifts using Planck 2018 cosmological parameters.

program.f9065 lines

program bao_calculator
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  real(dp), parameter :: c_light = 299792.458_dp  ! km/s
  real(dp), parameter :: H0 = 67.4_dp             ! km/s/Mpc (Planck 2018)
  real(dp) :: Om, Ob, OL, z, r_s, d_V, d_A, d_H
  real(dp) :: z_drag, theta_s
  integer :: i

Om = 0.315_dp
  Ob = 0.0493_dp
  OL = 1.0_dp - Om

! Sound horizon at drag epoch (Eisenstein & Hu 1998 fitting formula)
  z_drag = 1059.94_dp  ! Planck 2018
  r_s = 147.09_dp      ! Mpc (Planck 2018 sound horizon)

print *, '=== BAO Scale Calculator ==='
  print *, ''
  print '(A,F8.2,A)', ' Sound horizon r_s = ', r_s, ' Mpc'
  print '(A,F8.2)',   ' z_drag = ', z_drag
  print '(A,F6.3)',   ' Omega_m = ', Om
  print '(A,F6.4)',   ' Omega_b = ', Ob
  print '(A,F6.2,A)', ' H0 = ', H0, ' km/s/Mpc'
  print *, ''
  print *, 'BAO measurements at different redshifts:'
  print *, '  z       d_A(z)/r_s   H(z)*r_s     d_V(z)/r_s   theta_BAO [deg]'
  print *, '  -----   ----------   ----------   ----------   ---------------'

do i = 1, 8
    z = 0.1_dp + 0.3_dp * (i - 1)
    call compute_distances(z, Om, OL, H0, c_light, r_s, d_A, d_H, d_V, theta_s)
    print '(F7.2, 4F13.4)', z, d_A/r_s, H0*E_func(z,Om,OL)*r_s, d_V/r_s, theta_s
  end do

contains

function E_func(z, Om, OL) result(E)
    real(dp), intent(in) :: z, Om, OL
    real(dp) :: E
    E = sqrt(Om*(1+z)**3 + OL)
  end function

subroutine compute_distances(z, Om, OL, H0, c, rs, dA, dH, dV, theta)
    real(dp), intent(in) :: z, Om, OL, H0, c, rs
    real(dp), intent(out) :: dA, dH, dV, theta
    real(dp) :: dc, dz, E_val, sum_val
    integer :: j, nsteps

! Comoving distance via trapezoidal integration
    nsteps = 10000
    dz = z / nsteps
    sum_val = 0.0_dp
    do j = 1, nsteps
      sum_val = sum_val + dz / E_func((j-0.5_dp)*dz, Om, OL)
    end do
    dc = (c / H0) * sum_val

dA = dc / (1.0_dp + z)
    dH = c / (H0 * E_func(z, Om, OL))
    dV = (dA**2 * (1+z)**2 * dH * z)**(1.0_dp/3.0_dp)
    theta = rs / dA * 180.0_dp / 3.14159265358979_dp
  end subroutine

end program bao_calculator

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

4. Hubble Constant and Tensions

Local Measurements

SH0ES collaboration (Cepheids + SNe Ia):

$$H_0 = 73.04 \pm 1.04 \text{ km s}^{-1}\text{ Mpc}^{-1}$$

CMB-based Inference

Planck 2018 (assuming ΛCDM):

$$H_0 = 67.4 \pm 0.5 \text{ km s}^{-1}\text{ Mpc}^{-1}$$

Hubble Tension

The $\sim 5\sigma$ discrepancy between early-universe and late-universe measurements presents a major challenge. Possible explanations:

• Systematic errors in distance ladder or CMB analysis
• New physics: early dark energy, modified gravity, extra relativistic species
• Late-time modifications: evolving dark energy $w(z)$

Hubble Tension: Measurement Comparison

Python

Compilation of H0 measurements from early-universe (CMB, BAO) and late-universe (Cepheids, lensing, TRGB) probes showing the ~5-sigma tension.

script.py48 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Hubble constant measurements compilation
measurements = [
    ('Planck 2018 (CMB)', 67.4, 0.5, '#60a5fa', 'Early Universe'),
    ('ACT+WMAP', 67.6, 1.1, '#818cf8', 'Early Universe'),
    ('SPT-3G', 67.5, 1.2, '#a78bfa', 'Early Universe'),
    ('BAO + BBN', 67.6, 0.8, '#c084fc', 'Early Universe'),
    ('DES Y3 + BAO + BBN', 67.9, 0.8, '#d946ef', 'Early Universe'),
    ('SH0ES 2022 (Cepheids)', 73.0, 1.0, '#f472b6', 'Late Universe'),
    ('H0LiCOW (Lensing)', 73.3, 1.8, '#fb7185', 'Late Universe'),
    ('CCHP (TRGB)', 69.8, 1.7, '#f59e0b', 'Late Universe'),
    ('Megamasers', 73.9, 3.0, '#fbbf24', 'Late Universe'),
    ('SBF', 73.3, 2.5, '#34d399', 'Late Universe'),
    ('TDCOSMO', 74.2, 1.6, '#fb923c', 'Late Universe'),
]

fig, ax = plt.subplots(figsize=(12, 7))
fig.patch.set_facecolor('#0a0a1a')
ax.set_facecolor('#0a0a1a')
for sp in ax.spines.values(): sp.set_edgecolor('#334155')
ax.tick_params(colors='#94a3b8')

y_pos = np.arange(len(measurements))
for i, (name, val, err, color, epoch) in enumerate(measurements):
    ax.errorbar(val, i, xerr=err, fmt='o', color=color, markersize=8, capsize=5, capthick=2, lw=2)
    ax.text(val + err + 0.5, i, f'{val:.1f} +/- {err:.1f}', color=color, va='center', fontsize=9)

ax.set_yticks(y_pos)
ax.set_yticklabels([m[0] for m in measurements], color='#94a3b8', fontsize=10)
ax.axvline(67.4, color='#60a5fa', ls='--', alpha=0.5, label='Planck (67.4)')
ax.axvline(73.0, color='#f472b6', ls='--', alpha=0.5, label='SH0ES (73.0)')
ax.axvspan(66.9, 67.9, alpha=0.1, color='#60a5fa')
ax.axvspan(72.0, 74.0, alpha=0.1, color='#f472b6')
ax.set_xlabel(r'$H_0$ [km/s/Mpc]', color='#94a3b8', fontsize=12)
ax.set_title('The Hubble Tension', color='white', fontsize=14)
ax.legend(facecolor='#1e293b', labelcolor='white', fontsize=9)
ax.set_xlim(63, 80)
ax.grid(alpha=0.2, color='#334155', axis='x')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')
plt.close()
print('Hubble Tension visualization computed')
print(f'Tension: (73.0 - 67.4) / sqrt(1.0^2 + 0.5^2) = {(73.0-67.4)/np.sqrt(1.0**2+0.5**2):.1f} sigma')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

5. Gravitational Lensing

Weak Lensing

Cosmic shear from large-scale structure distorts background galaxy images. The convergence and shear are:

$$\kappa(\theta) = \int_0^{\chi_H} d\chi \, W(\chi) \delta(\chi\theta, \chi)$$

where the lensing weight function is:

$$W(\chi) = \frac{3H_0^2\Omega_m}{2c^2}\frac{\chi}{a(\chi)}\int_{\chi}^{\chi_H}d\chi' \, n(\chi')\frac{\chi' - \chi}{\chi'}$$

Strong Lensing

The lens equation relates source and image positions:

$$\vec{\beta} = \vec{\theta} - \vec{\alpha}(\vec{\theta})$$

where $\vec{\alpha}$ is the deflection angle. For a point mass:

$$\alpha = \frac{4GM}{c^2 b}$$

where $b$ is the impact parameter.

Gravitational Lensing: SIS Model

Python

Singular Isothermal Sphere lens model showing image positions, Einstein ring, and magnification as a function of source position.

script.py72 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Gravitational lensing simulation
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
fig.patch.set_facecolor('#0a0a1a')
for ax in [ax1, ax2]:
    ax.set_facecolor('#0a0a1a')
    for sp in ax.spines.values(): sp.set_edgecolor('#334155')
    ax.tick_params(colors='#94a3b8')

# Einstein ring and critical curves
theta = np.linspace(0, 2*np.pi, 500)
# SIS lens model
theta_E = 1.0  # Einstein radius in arcsec

# Source positions
beta_values = [0.0, 0.3, 0.6, 0.9, 1.2]
colors = ['#f472b6', '#a78bfa', '#60a5fa', '#34d399', '#fbbf24']

for beta, color in zip(beta_values, colors):
    # Image positions for SIS: theta = beta +/- theta_E
    if beta == 0:
        # Einstein ring
        x_ring = theta_E * np.cos(theta)
        y_ring = theta_E * np.sin(theta)
        ax1.plot(x_ring, y_ring, color=color, lw=2, label=f'beta={beta:.1f} (ring)')
    else:
        theta_plus = beta + theta_E
        theta_minus = beta - theta_E
        ax1.plot(theta_plus, 0, 'o', color=color, ms=8)
        if theta_minus > 0:
            ax1.plot(-theta_minus, 0, 's', color=color, ms=6)
        ax1.plot(beta, 0, 'x', color=color, ms=10, mew=2)

ax1.plot(0, 0, '*', color='white', ms=15, label='Lens center')
circle = plt.Circle((0, 0), theta_E, fill=False, color='white', ls='--', alpha=0.3)
ax1.add_patch(circle)
ax1.set_xlabel('theta_x [arcsec]', color='#94a3b8')
ax1.set_ylabel('theta_y [arcsec]', color='#94a3b8')
ax1.set_title('SIS Lens: Image Positions', color='white')
ax1.legend(facecolor='#1e293b', labelcolor='white', fontsize=8)
ax1.set_xlim(-2.5, 2.5)
ax1.set_ylim(-2.5, 2.5)
ax1.set_aspect('equal')
ax1.grid(alpha=0.2, color='#334155')

# Magnification vs source position
beta_arr = np.linspace(0.01, 3.0, 500)
# SIS magnification
mu_plus = 1.0 + theta_E / beta_arr
mu_minus = np.where(beta_arr < theta_E, -1.0 + theta_E / beta_arr, 0)
mu_total = mu_plus + np.abs(mu_minus)

ax2.plot(beta_arr, mu_total, color='#60a5fa', lw=2, label='Total |mu|')
ax2.plot(beta_arr, mu_plus, color='#34d399', lw=1.5, ls='--', label='mu+ (primary)')
ax2.plot(beta_arr, np.abs(mu_minus), color='#f472b6', lw=1.5, ls='--', label='|mu-| (secondary)')
ax2.axvline(theta_E, color='white', ls=':', alpha=0.3, label='Einstein radius')
ax2.set_xlabel('Source position beta/theta_E', color='#94a3b8')
ax2.set_ylabel('Magnification |mu|', color='#94a3b8')
ax2.set_title('SIS Lens: Magnification', color='white')
ax2.legend(facecolor='#1e293b', labelcolor='white', fontsize=8)
ax2.set_ylim(0, 10)
ax2.grid(alpha=0.2, color='#334155')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')
plt.close()
print(f'Einstein radius: {theta_E:.1f} arcsec')
print(f'Magnification at beta=0.5*theta_E: {1+theta_E/0.5 + abs(-1+theta_E/0.5):.1f}x')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Time Delay Cosmography

Multiple images of lensed quasars arrive at different times. The time delay constrains $H_0$:

$$\Delta t = \frac{D_{\Delta t}}{c}[(1+z_L)\Delta\phi - \frac{1}{2}\Delta|\vec{\theta}|^2]$$

where $D_{\Delta t} = \frac{D_L D_S}{D_{LS}}$ is the time-delay distance.

6. Large-Scale Structure

Matter Power Spectrum

The power spectrum quantifies clustering as a function of scale:

$$P(k) = \langle|\delta_k|^2\rangle = A k^{n_s} T^2(k)$$

where $n_s \approx 0.96$ is the spectral index and $T(k)$ is the transfer function.

Matter Power Spectrum (Fortran - BBKS Transfer Function)

Fortran

Computes the matter power spectrum P(k) using the BBKS transfer function approximation with Planck 2018 cosmological parameters.

program.f9049 lines

program matter_power_spectrum
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  real(dp), parameter :: pi = 3.14159265358979_dp
  real(dp) :: k, Pk, T_k, delta_H
  real(dp) :: Om, Ob, h, n_s, sigma8
  real(dp) :: q, Gamma
  integer :: i

! Cosmological parameters (Planck 2018)
  Om = 0.315_dp
  Ob = 0.0493_dp
  h = 0.674_dp
  n_s = 0.965_dp
  sigma8 = 0.811_dp
  Gamma = Om * h * exp(-Ob * (1.0_dp + sqrt(2.0_dp*h)/Om))

print *, '=== Matter Power Spectrum (BBKS Transfer Function) ==='
  print *, ''
  print '(A,F6.3)', ' Omega_m = ', Om
  print '(A,F6.4)', ' Omega_b = ', Ob
  print '(A,F6.3)', ' h = ', h
  print '(A,F6.3)', ' n_s = ', n_s
  print '(A,F6.3)', ' sigma_8 = ', sigma8
  print '(A,F6.4)', ' Shape Gamma = ', Gamma
  print *, ''
  print *, '  k [h/Mpc]     T(k)        P(k) [Mpc/h]^3'
  print *, '  ----------     ----        ---------------'

do i = 1, 30
    k = 10.0_dp ** (-3.0_dp + 0.2_dp * (i-1))
    q = k / Gamma
    ! BBKS transfer function
    T_k = log(1.0_dp + 2.34_dp*q) / (2.34_dp*q) * &
          (1.0_dp + 3.89_dp*q + (16.1_dp*q)**2 + (5.46_dp*q)**3 + (6.71_dp*q)**4)**(-0.25_dp)
    ! Power spectrum P(k) = A * k^n_s * T(k)^2
    Pk = k**n_s * T_k**2
    ! Normalize (approximate)
    Pk = Pk * (sigma8/0.8_dp)**2 * 2.0e4_dp

print '(ES13.3, F12.6, ES18.6)', k, T_k, Pk
  end do

print *, ''
  print *, 'Key scales:'
  print '(A,F8.4,A)', ' k_eq ~ ', Om*h**2 / 3000.0_dp, ' h/Mpc (matter-radiation equality)'
  print '(A,F8.4,A)', ' k_silk ~ ', 0.1_dp, ' h/Mpc (Silk damping scale)'

end program matter_power_spectrum

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Growth of Structure

Linear growth is characterized by the growth function $D(a)$:

$$\delta(k, a) = D(a) \delta(k, a_i)$$

The growth rate parameter:

$$f(a) = \frac{d\ln D}{d\ln a} \approx \Omega_m(a)^{0.55}$$

Redshift-Space Distortions

Peculiar velocities distort clustering measurements in redshift space:

$$P^s(k, \mu) = \left(1 + \beta\mu^2\right)^2 P^r(k)$$

where $\beta = f/b$, with $b$ the galaxy bias.

7. Modern Galaxy Surveys

Major Surveys

SDSS (Sloan Digital Sky Survey): Mapped millions of galaxies, foundational for BAO
DES (Dark Energy Survey): Weak lensing, galaxy clusters, SNe Ia
DESI (Dark Energy Spectroscopic Instrument): 40 million galaxy redshifts
Euclid: Space-based weak lensing and galaxy clustering
LSST/Vera Rubin: 20 billion galaxies, time-domain astronomy

21cm Cosmology

Neutral hydrogen 21cm line traces matter at high redshift:

$$T_b(z) \approx 27 x_{HI}(1+\delta)\left(\frac{H(z)}{dv_{\parallel}/dr_{\parallel}}\right)\left(\frac{1+z}{10}\right)^{1/2} \text{ mK}$$

This probes the Dark Ages and Epoch of Reionization.

8. ΛCDM Best-Fit Parameters

Planck 2018 + BAO + Pantheon SNe (TT,TE,EE+lowE+lensing):

$\Omega_b h^2 = 0.02237 \pm 0.00015$

$\Omega_c h^2 = 0.1200 \pm 0.0012$

$H_0 = 67.4 \pm 0.5$ km/s/Mpc

$\tau = 0.054 \pm 0.007$

$n_s = 0.9649 \pm 0.0042$

$\ln(10^{10}A_s) = 3.044 \pm 0.014$

Derived Parameters

$\Omega_m = 0.315 \pm 0.007$

$\Omega_\Lambda = 0.685 \pm 0.007$

$\sigma_8 = 0.811 \pm 0.006$

$t_0 = 13.80 \pm 0.02$ Gyr

These parameters describe a spatially flat universe dominated by dark energy (68.5%), with dark matter (26.8%) and ordinary matter (4.9%). The universe is 13.8 billion years old.

Observational Cosmology

Introduction

1. Cosmic Microwave Background (CMB)

Discovery and Temperature

Temperature Anisotropies

Angular Power Spectrum

CMB Angular Power Spectrum Simulation

Physical Interpretation of Peaks

Polarization

2. Cosmic Distance Ladder

Standard Candles

Type Ia Supernovae

Cepheid Variables

Cosmic Distance Measures and Hubble Parameter

3. Baryon Acoustic Oscillations (BAO)

Physical Origin

Two-Point Correlation Function

Cosmological Constraints

BAO Scale Calculator (Fortran)

4. Hubble Constant and Tensions

Local Measurements

CMB-based Inference

Hubble Tension

Hubble Tension: Measurement Comparison

5. Gravitational Lensing

Weak Lensing

Strong Lensing

Gravitational Lensing: SIS Model

Time Delay Cosmography

6. Large-Scale Structure

Matter Power Spectrum

Matter Power Spectrum (Fortran - BBKS Transfer Function)

Growth of Structure

Redshift-Space Distortions

7. Modern Galaxy Surveys

Major Surveys

21cm Cosmology

8. ΛCDM Best-Fit Parameters

Derived Parameters

Related Courses

Cosmology

Dark Matter

Dark Energy

Astrophysics