Dark Energy

A comprehensive graduate-level course on dark energy — from the 1998 supernova observations through the cosmological constant, dynamical dark energy models, observational probes, and the ultimate fate of the universe — with full MathJax derivations and Python simulations.

1. Introduction & Observational Evidence

In 1998, two independent teams — the Supernova Cosmology Project (Perlmutter et al.) and the High-z Supernova Search Team (Riess et al., Schmidt et al.) — discovered that the expansion of the universe is accelerating. This was awarded the 2011 Nobel Prize in Physics. Type Ia supernovae serve as standardizable candles because their peak luminosity correlates tightly with the width of their light curves (the Phillips relation). By measuring the apparent brightness of distant Type Ia supernovae and comparing to their known intrinsic luminosity, one infers their luminosity distances. The data showed that distant supernovae are fainter than expected in a decelerating universe — the expansion is speeding up.

The Acceleration Condition

From general relativity, the second Friedmann equation (the Raychaudhuri equation) gives the acceleration of the scale factor:

$$\frac{\ddot{a}}{a} = -\frac{4\pi G}{3}\left(\rho + 3p\right)$$

For acceleration ($\ddot{a} > 0$), we need the right-hand side to be positive:

$$\rho + 3p < 0 \quad \Longrightarrow \quad p < -\frac{\rho}{3}$$

This requires a component with strongly negative pressure. For a perfect fluid with equation of state $w = p/\rho$, the condition becomes $w < -1/3$. Ordinary matter ($w = 0$) and radiation ($w = 1/3$) both give deceleration. Something else is driving the acceleration — we call it dark energy.

Key Observational Milestones

●1998: Riess et al. and Perlmutter et al. publish evidence for cosmic acceleration from Type Ia supernovae
●2003: WMAP first-year data confirms a flat universe dominated by dark energy (~73%)
●2005: Baryon Acoustic Oscillation (BAO) detection in SDSS galaxy survey provides independent confirmation
●2011: Nobel Prize in Physics to Perlmutter, Schmidt, and Riess
●2018: Planck final results: $\Omega_\Lambda = 0.6847 \pm 0.0073$, $H_0 = 67.36 \pm 0.54\;\text{km/s/Mpc}$
●2024: DESI BAO results hint at possible time-varying dark energy ($w_0 w_a$CDM tension with $\Lambda$CDM at ~2.5σ)

The Cosmic Energy Budget

Modern observations consistently point to a universe composed of approximately:

68.5%

Dark Energy

26.5%

Dark Matter

Baryonic Matter

2. The Friedmann Equations

The Friedmann equations govern the expansion dynamics of a homogeneous, isotropic universe. They follow from the Einstein field equations applied to the Friedmann-Lemaître-Robertson-Walker (FLRW) metric.

The FLRW Metric

The most general metric for a homogeneous, isotropic spacetime is:

$$ds^2 = -c^2 dt^2 + a(t)^2\left[\frac{dr^2}{1-kr^2} + r^2\left(d\theta^2 + \sin^2\theta\;d\phi^2\right)\right]$$

where $a(t)$ is the scale factor, $k$ is the spatial curvature ($k = 0, +1, -1$ for flat, closed, or open geometry), and the Hubble parameter is defined as $H \equiv \dot{a}/a$.

Derivation from Einstein’s Equations

The Einstein field equations with cosmological constant are:

$$G_{\mu\nu} + \Lambda g_{\mu\nu} = \frac{8\pi G}{c^4}\,T_{\mu\nu}$$

For a perfect fluid at rest in comoving coordinates, the energy-momentum tensor is$T^{\mu\nu} = \text{diag}(\rho c^2, p, p, p)$. Computing the Ricci tensor and scalar curvature from the FLRW metric and substituting into Einstein’s equations, the $00$-component gives the first Friedmann equation:

$$H^2 = \frac{8\pi G}{3}\rho - \frac{kc^2}{a^2} + \frac{\Lambda}{3}$$

The $ij$-components, combined with the first equation, yield the second Friedmann equation (acceleration equation):

$$\dot{H} + H^2 = \frac{\ddot{a}}{a} = -\frac{4\pi G}{3}\left(\rho + \frac{3p}{c^2}\right) + \frac{\Lambda}{3}$$

In natural units ($c = 1$), these simplify to:

$$\boxed{H^2 = \frac{8\pi G}{3}\rho - \frac{k}{a^2} + \frac{\Lambda}{3}}$$

$$\boxed{\frac{\ddot{a}}{a} = -\frac{4\pi G}{3}\left(\rho + 3p\right) + \frac{\Lambda}{3}}$$

The Fluid Equation and Density Parameters

Combining the two Friedmann equations gives the continuity (fluid) equation:

$$\dot{\rho} + 3H\left(\rho + p\right) = 0$$

For a component with constant equation of state $w = p/\rho$, this integrates to:

$$\rho(a) = \rho_0\,a^{-3(1+w)}$$

Key cases: matter ($w = 0$, $\rho \propto a^{-3}$), radiation ($w = 1/3$, $\rho \propto a^{-4}$), cosmological constant ($w = -1$, $\rho = \text{const}$).

We define the dimensionless density parameters:

$$\Omega_i = \frac{\rho_i}{\rho_\text{crit}}, \quad \rho_\text{crit} = \frac{3H^2}{8\pi G}$$

so the first Friedmann equation becomes $\Omega_m + \Omega_r + \Omega_\Lambda + \Omega_k = 1$, where $\Omega_k = -k/(aH)^2$. Current best values (Planck 2018): $\Omega_m \approx 0.315$, $\Omega_\Lambda \approx 0.685$,$|\Omega_k| < 0.002$.

Hubble Parameter as a Function of Redshift

Using $a = 1/(1+z)$ and the density scaling relations, the Friedmann equation in terms of redshift becomes:

$$H(z) = H_0\sqrt{\Omega_r(1+z)^4 + \Omega_m(1+z)^3 + \Omega_k(1+z)^2 + \Omega_\text{DE}(1+z)^{3(1+w)}}$$

For $\Lambda$CDM with $w = -1$ and flat geometry ($\Omega_k = 0$), this simplifies to:

$$H(z) = H_0\sqrt{\Omega_m(1+z)^3 + \Omega_\Lambda}$$

3. The Cosmological Constant

The simplest and observationally most successful model of dark energy is the cosmological constant $\Lambda$, originally introduced by Einstein in 1917 to maintain a static universe, then abandoned after Hubble’s 1929 discovery of expansion, and resurrected after the 1998 supernova discovery. In the $\Lambda$CDM model, it is the dominant energy component today.

$\Lambda$ as Vacuum Energy

The cosmological constant can be interpreted as the energy density of the vacuum. The$\Lambda$ term in the Einstein equations acts as a perfect fluid with:

$$\rho_\Lambda = \frac{\Lambda}{8\pi G}, \quad p_\Lambda = -\rho_\Lambda c^2$$

The equation of state is exactly:

$$\boxed{w_\Lambda = \frac{p_\Lambda}{\rho_\Lambda c^2} = -1}$$

Proof that $\rho_\Lambda$ is constant: Substituting$w = -1$ into the fluid equation:

$$\dot{\rho}_\Lambda + 3H(\rho_\Lambda - \rho_\Lambda) = \dot{\rho}_\Lambda = 0$$

Hence $\rho_\Lambda$ remains constant as the universe expands, in stark contrast to matter and radiation which dilute as $a^{-3}$ and $a^{-4}$. Using the observed value:

$$\rho_\Lambda = \frac{\Lambda}{8\pi G} \approx 6.9 \times 10^{-27}\;\text{kg/m}^3 \approx 3.8\;\text{GeV/m}^3$$

The Cosmological Constant Problem

This is arguably the worst prediction in all of physics. In quantum field theory, every field contributes zero-point energy to the vacuum. For a field of mass $m$, the vacuum energy density up to a UV cutoff $\Lambda_\text{UV}$ is:

$$\rho_\text{vac}^\text{QFT} \sim \frac{\Lambda_\text{UV}^4}{\hbar^3 c^5}$$

Taking the Planck scale as the natural cutoff ($\Lambda_\text{UV} = E_\text{Pl} \sim 10^{19}\;\text{GeV}$):

$$\rho_\text{vac}^\text{QFT} \sim 10^{76}\;\text{GeV}^4 \sim 10^{112}\;\text{erg/cm}^3$$

The observed dark energy density is:

$$\rho_\Lambda^\text{obs} \sim 10^{-47}\;\text{GeV}^4 \sim 10^{-8}\;\text{erg/cm}^3$$

The ratio is staggering:

$$\boxed{\frac{\rho_\text{vac}^\text{QFT}}{\rho_\Lambda^\text{obs}} \sim 10^{120}}$$

This 120-order-of-magnitude discrepancy is the cosmological constant problem. Even using the electroweak scale ($\Lambda_\text{UV} \sim 100\;\text{GeV}$) as a cutoff, the discrepancy is still $\sim 10^{55}$. No known symmetry or mechanism can explain why the vacuum energy is so exquisitely fine-tuned to nearly zero but not exactly zero.

The Coincidence Problem

Why are the dark energy and matter densities comparable today? Since$\rho_m \propto a^{-3}$ while $\rho_\Lambda = \text{const}$, their ratio changes dramatically with time:

$$\frac{\rho_\Lambda}{\rho_m} = \frac{\Omega_\Lambda}{\Omega_m}(1+z)^{3} \approx 2.2\,(1+z)^3$$

At $z = 10$ this ratio was $\sim 10^{-3}$; at $z = 1000$(CMB) it was $\sim 10^{-9}$. We happen to live in the cosmologically brief epoch where $\rho_\Lambda \sim \rho_m$ — this curious timing is the coincidence problem, one motivation for dynamical dark energy models.

4. Dark Energy Models Beyond $\Lambda$

While $\Lambda$CDM fits observations remarkably well, the cosmological constant problem and the coincidence problem motivate the exploration of dynamical dark energy. Here we derive the key models.

4.1 Quintessence: Scalar Field Dark Energy

Quintessence models dark energy as a slowly-rolling scalar field $\phi$, analogous to the inflaton in early-universe inflation. The Lagrangian density is:

$$\mathcal{L} = \frac{1}{2}g^{\mu\nu}\partial_\mu\phi\,\partial_\nu\phi - V(\phi) = \frac{1}{2}\dot{\phi}^2 - V(\phi)$$

where the second equality holds for a spatially homogeneous field. The energy density and pressure are:

$$\rho_\phi = \frac{1}{2}\dot{\phi}^2 + V(\phi), \quad p_\phi = \frac{1}{2}\dot{\phi}^2 - V(\phi)$$

The equation of state is:

$$\boxed{w_\phi = \frac{p_\phi}{\rho_\phi} = \frac{\dot{\phi}^2/2 - V(\phi)}{\dot{\phi}^2/2 + V(\phi)}}$$

Key properties:

●When $\dot{\phi}^2 \ll V(\phi)$ (slow-roll): $w_\phi \to -1$ (mimics $\Lambda$)
●When $\dot{\phi}^2 \gg V(\phi)$ (kinetic dominated): $w_\phi \to +1$ (stiff fluid)
●The range $-1 \leq w_\phi \leq +1$ is always satisfied for canonical quintessence

The field obeys the Klein-Gordon equation in an expanding background:

$$\ddot{\phi} + 3H\dot{\phi} + V'(\phi) = 0$$

The $3H\dot{\phi}$ term acts as Hubble friction. Common potentials include the exponential ($V = V_0 e^{-\lambda\phi/M_\text{Pl}}$, which admits scaling solutions where $\Omega_\phi$ is constant), the inverse power-law ($V = M^{4+\alpha}\phi^{-\alpha}$), and the cosine ($V = \Lambda^4[1 + \cos(\phi/f)]$, pseudo-Nambu-Goldstone boson).

4.2 Phantom Energy and the Big Rip

What if $w < -1$? This is phantom dark energy, modeled by a scalar field with a wrong-sign kinetic term:

$$\mathcal{L}_\text{phantom} = -\frac{1}{2}(\partial\phi)^2 - V(\phi)$$

For phantom energy with constant $w < -1$, the energy density grows with expansion:

$$\rho_\text{ph} \propto a^{-3(1+w)} = a^{|3(1+w)|} \quad (\text{since } 1+w < 0)$$

The scale factor diverges in finite time (Caldwell, Kamionkowski, Weinberg, 2003):

$$a(t) \propto (t_\text{rip} - t)^{2/[3(1+w)]}$$

The time to the Big Rip from today:

$$\boxed{t_\text{rip} - t_0 = \frac{2}{3|1+w|H_0\sqrt{1-\Omega_m}}}$$

For $w = -1.1$, this gives $t_\text{rip} - t_0 \sim 100\;\text{Gyr}$. In the approach to the Big Rip, gravitationally bound structures are torn apart: galaxy clusters, then galaxies, solar systems, planets, and ultimately atoms are ripped apart by the ever-increasing expansion rate. However, phantom fields violate the null energy condition and generically suffer from quantum instabilities, making this scenario theoretically problematic.

4.3 k-essence

k-essence generalizes quintessence by allowing a non-standard kinetic term. The Lagrangian is a general function of the field and its kinetic energy $X = \frac{1}{2}(\partial\phi)^2$:

$$\mathcal{L} = K(X, \phi)$$

The energy density and pressure are:

$$\rho = 2X K_X - K, \quad p = K$$

where $K_X = \partial K/\partial X$. k-essence can produce acceleration through purely kinetic effects, without a potential, and naturally achieves $w < -1/3$through the non-linear kinetic terms. It also produces a variable speed of sound$c_s^2 = K_X/(K_X + 2X K_{XX})$, which distinguishes it observationally from quintessence.

4.4 The CPL Parametrization

Rather than specifying a Lagrangian, one can phenomenologically parametrize the equation of state. The Chevallier-Polarski-Linder (CPL) form is the most widely used:

$$\boxed{w(a) = w_0 + w_a(1-a) = w_0 + w_a\frac{z}{1+z}}$$

Here $w_0$ is the present-day value and $w_a$ captures time evolution. For $\Lambda$CDM: $w_0 = -1$, $w_a = 0$. The dark energy density evolves as:

$$\rho_\text{DE}(a) = \rho_{\text{DE},0}\;a^{-3(1+w_0+w_a)}\;\exp\!\left[-3w_a(1-a)\right]$$

The DESI 2024 BAO results, combined with CMB and supernova data, found$w_0 = -0.55^{+0.39}_{-0.21}$ and $w_a = -1.32^{+0.58}_{-0.82}$ (at 68% CL), showing $\sim 2.5\sigma$ tension with $\Lambda$CDM in the$w_0$-$w_a$ plane.

4.5 Modified Gravity: f(R) Theories

Instead of adding a new energy component, one can modify gravity itself. In f(R) theories, the Einstein-Hilbert action is generalized:

$$S = \frac{1}{16\pi G}\int d^4x\,\sqrt{-g}\;f(R) + S_\text{matter}$$

where $f(R)$ is a general function of the Ricci scalar. The modified field equations are:

$$f'(R)R_{\mu\nu} - \frac{1}{2}f(R)g_{\mu\nu} + \left(g_{\mu\nu}\Box - \nabla_\mu\nabla_\nu\right)f'(R) = 8\pi G\,T_{\mu\nu}$$

Setting $f(R) = R - 2\Lambda$ recovers GR with a cosmological constant. Popular models include Starobinsky ($f(R) = R + \alpha R^2$) and Hu-Sawicki ($f(R) = R - \mu R_c\frac{(R/R_c)^n}{(R/R_c)^n + 1}$).

Crucially, f(R) theories are equivalent to scalar-tensor (Brans-Dicke) theories with parameter $\omega_\text{BD} = 0$, via the conformal transformation $\tilde{g}_{\mu\nu} = f'(R)g_{\mu\nu}$. The extra scalar degree of freedom $\varphi = f'(R)$ mediates a fifth force that must be screened on solar system scales (chameleon mechanism).

5. Observational Probes of Dark Energy

Constraining dark energy requires precise measurements of the expansion history and growth of structure. Multiple independent probes are combined to break parameter degeneracies.

5.1 Distance-Redshift Relations

The comoving distance to an object at redshift $z$ is:

$$\chi(z) = \int_0^z \frac{c\,dz'}{H(z')}$$

The luminosity distance (measured via standard candles like Type Ia SNe):

$$\boxed{d_L(z) = (1+z)\int_0^z \frac{c\,dz'}{H(z')}}$$

The angular diameter distance (measured via standard rulers like BAO):

$$d_A(z) = \frac{1}{1+z}\int_0^z \frac{c\,dz'}{H(z')} = \frac{d_L(z)}{(1+z)^2}$$

These relations encode the full expansion history $H(z)$, which depends on all the density parameters and the dark energy equation of state.

5.2 Type Ia Supernovae

The distance modulus relates the apparent magnitude $m$ to the absolute magnitude $M$:

$$\mu = m - M = 5\log_{10}\!\left(\frac{d_L}{\text{Mpc}}\right) + 25$$

Modern compilations (Pantheon+, DES-SN5YR, Union3) contain $\sim 1500$ supernovae spanning $0.01 < z < 2.3$. They constrain $\Omega_m$and $w$ primarily through the shape of $d_L(z)$.

5.3 Baryon Acoustic Oscillations (BAO)

BAO are the frozen imprint of sound waves in the baryon-photon plasma before recombination. The sound horizon at the drag epoch is:

$$r_d = \int_0^{t_d} \frac{c_s\,dt}{a} = \int_{z_d}^\infty \frac{c_s\,dz}{H(z)} \approx 147\;\text{Mpc}$$

This is a standard ruler. Measuring the BAO scale at multiple redshifts in the galaxy power spectrum gives both $d_A(z)/r_d$ (transverse) and $c/(H(z)r_d)$(line-of-sight). DESI, SDSS-IV, and Euclid provide BAO measurements at $z = 0.1$ to $z = 4.2$, tightly constraining the expansion history.

5.4 CMB Power Spectrum

The CMB constrains dark energy mainly through two effects: (1) the angular diameter distance to the last scattering surface at $z_* \approx 1090$, which determines the angular scale of the acoustic peaks via $\theta_* = r_*/d_A(z_*)$; and (2) the integrated Sachs-Wolfe (ISW) effect — late-time evolution of gravitational potentials due to dark energy produces excess power at large angular scales ($\ell \lesssim 30$).

The ISW contribution to the temperature anisotropy is:

$$\frac{\Delta T}{T}\bigg|_\text{ISW} = -2\int_0^{\chi_*}\dot{\Phi}\;d\chi$$

5.5 Weak Lensing and Cluster Counts

Weak gravitational lensing measures the distortion of background galaxy shapes by foreground mass distributions. The lensing convergence power spectrum depends on both the geometry (distance ratios) and the growth of structure:

$$C_\ell^{\kappa\kappa} = \int_0^{\chi_*} d\chi\;\frac{W(\chi)^2}{\chi^2}\;P_\delta\!\left(\frac{\ell}{\chi}, z(\chi)\right)$$

Galaxy cluster counts probe dark energy through the growth rate of structure. The number of clusters above mass $M$ per unit redshift:

$$\frac{dN}{dz} = \frac{dV}{dz}\int_M^\infty \frac{dn}{dM'}\,dM'$$

where $dn/dM$ is the halo mass function. Both probes are sensitive to the combination $S_8 = \sigma_8\sqrt{\Omega_m/0.3}$, providing complementary constraints to geometric probes.

6. The Fate of the Universe

The ultimate fate of the universe depends critically on the nature of dark energy, specifically on the value of $w$ and whether it evolves with time.

Scenario 1: Big Freeze / Heat Death ($w = -1$)

In the $\Lambda$CDM model, the universe asymptotes to de Sitter space. The scale factor grows exponentially:

$$a(t) \to e^{H_\infty t}, \quad H_\infty = H_0\sqrt{\Omega_\Lambda} \approx 56\;\text{km/s/Mpc}$$

The cosmic event horizon forms at a comoving distance $\chi_\text{EH} = c/H_\infty \approx 17\;\text{Gpc}$. Objects beyond this horizon can never send signals to us. Galaxies outside the Local Group recede beyond the horizon within $\sim 100\;\text{Gyr}$. The universe approaches maximum entropy — the heat death.

Scenario 2: Big Rip ($w < -1$)

For phantom dark energy, the scale factor diverges in a finite time $t_\text{rip}$. As $t \to t_\text{rip}$, the Hubble parameter diverges $H \to \infty$, and all bound structures are progressively disrupted. The approximate timeline before the rip (for $w = -1.5$):

●$\sim 60$ Myr before: Galaxy clusters fly apart
●$\sim 60$ Myr before: Milky Way disrupted
●$\sim 3$ months before: Solar system unbound
●$\sim 30$ minutes before: Earth destroyed
●$\sim 10^{-19}$ seconds before: Atoms dissociate

Scenario 3: Big Crunch ($w > -1/3$ or Dark Energy Decaying)

If dark energy decays or has $w > -1/3$, it eventually ceases to drive acceleration. In a matter-dominated closed universe ($k = +1$), the expansion halts and reverses:

$$a(t) = a_\text{max}\sin^{2/3}\!\left(\frac{3H_0\sqrt{\Omega_m}}{2}t\right)$$

The universe recollapses to a singularity. Some cyclic models (Steinhardt-Turok) envision a series of big bangs and crunches. In practice, current observations strongly disfavor$w > -1/3$ for the present dark energy component.

The Deceleration Parameter

The deceleration parameter quantifies whether expansion is accelerating or decelerating:

$$q = -\frac{\ddot{a}a}{\dot{a}^2} = -1 - \frac{\dot{H}}{H^2}$$

For a flat universe with matter and dark energy ($w = \text{const}$):

$$q(z) = \frac{1}{2}\frac{\Omega_m(1+z)^3 + (1+3w)\,\Omega_\text{DE}(1+z)^{3(1+w)}}{\Omega_m(1+z)^3 + \Omega_\text{DE}(1+z)^{3(1+w)}}$$

The transition from deceleration ($q > 0$) to acceleration ($q < 0$) occurred at the transition redshift:

$$z_t = \left(\frac{-(1+3w)\Omega_\text{DE}}{\Omega_m}\right)^{-1/3w} - 1 \approx 0.67 \quad (\text{for } \Lambda\text{CDM})$$

7. Python Simulations

These simulations numerically solve the Friedmann equations and compute key cosmological observables for different dark energy equations of state.

7.1 Scale Factor Evolution for Different Dark Energy Models

Numerically integrate the Friedmann equation to find how the scale factor evolves with cosmic time for different values of the dark energy equation of state parameter w.

Scale Factor a(t) for Different w Values

Python

script.py79 lines

import numpy as np

# Friedmann equation: (da/dt)^2 = H0^2 * [Om * a^(-1) + Ode * a^(-(1+3w))]
# We solve da/dt = H0 * a * sqrt(Om/a^3 + Ode * a^(-3(1+w)))

H0 = 67.4  # km/s/Mpc
Om = 0.315
Ode = 1.0 - Om
c_km = 3e5  # km/s

# Convert H0 to 1/Gyr
H0_Gyr = H0 * 3.156e16 / 3.086e22  # ~ 0.0688 / Gyr

w_values = [-1.0, -0.8, -1.2, -1.5, -0.5]
labels = ["w = -1.0 (Lambda)", "w = -0.8 (quint.)", "w = -1.2 (phantom)", "w = -1.5 (phantom)", "w = -0.5 (quint.)"]

# Integrate da/dt = H0_Gyr * a * E(a)
# E(a) = sqrt(Om * a^-3 + Ode * a^(-3(1+w)))

dt = 0.005  # Gyr
t_max = 50.0
n_steps = int(t_max / dt)

print("=== Scale Factor Evolution a(t) ===")
print("Columns: t(Gyr), a_w-1.0, a_w-0.8, a_w-1.2, a_w-1.5, a_w-0.5")

# Store results
results = []
a_arr = [0.001] * len(w_values)  # start from small a

# Find t offset so that a=1 corresponds to ~13.8 Gyr
# First pass: evolve and record when a crosses 1
t_cross = [None] * len(w_values)

for step in range(n_steps):
    t = step * dt
    row = [t]
    for i, w in enumerate(w_values):
        a = a_arr[i]
        if a > 0 and a < 1e6:
            E = np.sqrt(max(Om * a**(-3) + Ode * a**(-3*(1+w)), 1e-30))
            da = H0_Gyr * a * E * dt
            a_arr[i] = a + da
            if t_cross[i] is None and a_arr[i] >= 1.0:
                t_cross[i] = t
        row.append(a_arr[i])
    results.append(row)

# Now re-do with time shifted so a=1 at t~13.8 Gyr
print()
for i, w in enumerate(w_values):
    tc = t_cross[i] if t_cross[i] else 13.8
    print(f"{labels[i]}: a=1 at t={tc:.1f} Gyr")

# Output sampled points for chart (shifted time)
a_arr2 = [0.001] * len(w_values)
t_offsets = [tc if tc else 13.8 for tc in t_cross]

print()
print("t_shifted a_LambdaCDM a_w-0.8 a_w-1.2 a_w-1.5 a_w-0.5")

for step in range(n_steps):
    t = step * dt
    for i, w in enumerate(w_values):
        a = a_arr2[i]
        if a > 0 and a < 1e8:
            E = np.sqrt(max(Om * a**(-3) + Ode * a**(-3*(1+w)), 1e-30))
            da = H0_Gyr * a * E * dt
            a_arr2[i] = a + da

if step % 100 == 0:
        t_shifted = t - t_offsets[0] + 13.8
        vals = " ".join(f"{a_arr2[i]:.4f}" for i in range(len(w_values)))
        print(f"{t_shifted:.2f} {vals}")

print()
print("Key: t_shifted is cosmic time in Gyr (today = 13.8 Gyr)")
print("Phantom (w<-1) models show faster acceleration; quintessence (w>-1) slower")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

7.2 Distance-Redshift Relation

Compute the luminosity distance $d_L(z)$ and distance modulus$\mu(z)$ for different dark energy models, showing how the supernova Hubble diagram constrains $w$.

Luminosity Distance and Hubble Diagram

Python

script.py63 lines

import numpy as np

H0 = 67.4   # km/s/Mpc
Om = 0.315
Ode = 1.0 - Om
c = 3.0e5   # km/s

def E_inv(z, w):
    """1/E(z) = 1/sqrt(Om*(1+z)^3 + Ode*(1+z)^(3(1+w)))"""
    zp1 = 1.0 + z
    val = Om * zp1**3 + Ode * zp1**(3*(1+w))
    return 1.0 / np.sqrt(max(val, 1e-30))

def luminosity_distance(z, w, n=500):
    """d_L in Mpc via trapezoidal integration"""
    zs = np.linspace(0, z, n+1)
    dz = z / n
    integral = 0.0
    for i in range(n):
        integral += 0.5 * (E_inv(zs[i], w) + E_inv(zs[i+1], w)) * dz
    return (c / H0) * (1 + z) * integral

def dist_modulus(dL_Mpc):
    return 5.0 * np.log10(max(dL_Mpc, 1e-10)) + 25.0

w_values = [-1.0, -0.8, -1.2, -0.5]
labels = ["w=-1.0", "w=-0.8", "w=-1.2", "w=-0.5"]

print("=== Luminosity Distance d_L(z) and Distance Modulus mu(z) ===")
print()

# Output table for chart
print("z dL_w-1.0 dL_w-0.8 dL_w-1.2 dL_w-0.5")
z_values = np.linspace(0.01, 2.0, 40)
for z in z_values:
    dLs = [luminosity_distance(z, w) for w in w_values]
    print(f"{z:.3f} " + " ".join(f"{d:.1f}" for d in dLs))

print()
print("--- Distance Modulus mu(z) ---")
print("z mu_w-1.0 mu_w-0.8 mu_w-1.2 mu_w-0.5")
for z in z_values:
    mus = [dist_modulus(luminosity_distance(z, w)) for w in w_values]
    print(f"{z:.3f} " + " ".join(f"{m:.3f}" for m in mus))

print()
print("Phantom (w=-1.2) gives LARGER distances => SNe appear FAINTER")
print("Quintessence (w=-0.8) gives SMALLER distances => SNe appear BRIGHTER")
print("This is how Type Ia supernovae constrain the dark energy equation of state")

# Compute differences from LambdaCDM
print()
print("--- Fractional difference in d_L from LambdaCDM ---")
print("z delta_w-0.8(%) delta_w-1.2(%) delta_w-0.5(%)")
for z in [0.1, 0.3, 0.5, 0.7, 1.0, 1.5, 2.0]:
    dL_ref = luminosity_distance(z, -1.0)
    deltas = [(luminosity_distance(z, w) - dL_ref) / dL_ref * 100 for w in [-0.8, -1.2, -0.5]]
    print(f"{z:.1f}  {deltas[0]:+.2f}  {deltas[1]:+.2f}  {deltas[2]:+.2f}")

print()
print("Percent-level differences require sub-percent distance precision")
print("=> need thousands of well-calibrated supernovae")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

7.3 Deceleration Parameter and Cosmic Acceleration

Compute the deceleration parameter $q(z)$ and identify the transition redshift where the universe switched from deceleration to acceleration.

Deceleration Parameter q(z)

Python

script.py61 lines

import numpy as np

Om = 0.315
Ode = 1.0 - Om

def q_of_z(z, w):
    """Deceleration parameter q(z) for flat universe with matter + DE"""
    zp1 = 1.0 + z
    rho_m = Om * zp1**3
    rho_de = Ode * zp1**(3*(1+w))
    E2 = rho_m + rho_de
    q = (0.5 * rho_m + 0.5 * (1 + 3*w) * rho_de) / E2
    return q

def find_zt(w, z_min=0.0, z_max=5.0, tol=1e-6):
    """Find transition redshift where q(z_t) = 0 by bisection"""
    if q_of_z(z_min, w) > 0:
        return None  # already accelerating at z_min, shouldn't happen
    # q should be negative at z=0 for DE-dominated, positive at high z
    if q_of_z(z_max, w) < 0:
        return None  # always accelerating (phantom)
    lo, hi = z_min, z_max
    for _ in range(100):
        mid = (lo + hi) / 2
        if q_of_z(mid, w) < 0:
            lo = mid
        else:
            hi = mid
        if hi - lo < tol:
            break
    return (lo + hi) / 2

w_values = [-1.0, -0.8, -1.2, -0.5, -1.5]
labels = ["w=-1.0", "w=-0.8", "w=-1.2", "w=-0.5", "w=-1.5"]

print("=== Deceleration Parameter q(z) ===")
print()

# Transition redshifts
print("--- Transition Redshifts (q = 0) ---")
for w, label in zip(w_values, labels):
    zt = find_zt(w)
    if zt is not None:
        print(f"{label}: z_t = {zt:.4f} (acceleration began at z < z_t)")
    else:
        print(f"{label}: always accelerating (phantom)")
print()

# Output table for chart
print("z q_w-1.0 q_w-0.8 q_w-1.2 q_w-0.5 q_w-1.5")
z_values = np.linspace(0.0, 3.0, 60)
for z in z_values:
    qs = [q_of_z(z, w) for w in w_values]
    print(f"{z:.3f} " + " ".join(f"{q:.4f}" for q in qs))

print()
print("q < 0 means acceleration, q > 0 means deceleration")
print("The universe transitioned from deceleration to acceleration at z ~ 0.67")
print("More negative w => earlier and stronger acceleration")
print("Quintessence (w > -1) => later transition, weaker acceleration")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

7.4 Energy Density Evolution and Dark Energy Domination

Track how the fractional energy densities $\Omega_i(z)$ evolve, showing the transitions between radiation, matter, and dark energy domination.

Cosmic Energy Density Evolution

Python

script.py54 lines

import numpy as np

Om0 = 0.315
Or0 = 9.15e-5   # radiation (photons + 3 massless neutrinos)
Ode0 = 1.0 - Om0 - Or0

print("=== Cosmic Energy Density Evolution ===")
print()
print("Current values:")
print(f"  Omega_matter  = {Om0:.4f}")
print(f"  Omega_rad     = {Or0:.6f}")
print(f"  Omega_DE      = {Ode0:.4f}")
print()

# For w=-1 (LambdaCDM)
def omegas(z, w=-1.0):
    zp1 = 1.0 + z
    rho_m = Om0 * zp1**3
    rho_r = Or0 * zp1**4
    rho_de = Ode0 * zp1**(3*(1+w))
    E2 = rho_m + rho_r + rho_de
    return rho_m/E2, rho_r/E2, rho_de/E2

# Find equality epochs
# Matter-radiation equality
z_eq_mr = Om0 / Or0 - 1
print(f"Matter-radiation equality: z_eq = {z_eq_mr:.0f}")

# Matter-DE equality (LambdaCDM)
z_eq_mde = (Ode0 / Om0)**(1.0/3) - 1
print(f"Matter-DE equality: z_eq = {z_eq_mde:.2f}")
print()

# Output for chart
print("z Omega_m Omega_r Omega_DE")
# Log-spaced redshifts from z=0 to z=10000
z_vals = np.concatenate([
    np.linspace(0, 2, 20),
    np.logspace(np.log10(3), 4, 30)
])

for z in z_vals:
    om, or_, ode = omegas(z)
    print(f"{z:.2f} {om:.6f} {or_:.6f} {ode:.6f}")

print()
print("=== Epoch Summary ===")
print("z > 3400: Radiation dominated (Omega_r > 0.5)")
print("0.3 < z < 3400: Matter dominated (Omega_m > 0.5)")
print("z < 0.3: Dark energy dominated (Omega_DE > 0.5)")
print()
print("The coincidence problem: we live near z ~ 0.3 where Omega_m ~ Omega_DE")
print("This is a very brief epoch in cosmic history")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

7.5 Big Rip Timescale for Phantom Dark Energy

For phantom dark energy ($w < -1$), compute the time remaining until the Big Rip and show how different structures are disrupted.

Big Rip Timeline

Python

script.py79 lines

import numpy as np

H0_Gyr = 67.4 * 3.156e16 / 3.086e22  # H0 in 1/Gyr
Om = 0.315
Ode = 1.0 - Om
t_now = 13.8  # Gyr

print("=== Big Rip Analysis for Phantom Dark Energy ===")
print()

# Time to Big Rip: t_rip - t_0 = 2 / (3|1+w| * H0 * sqrt(Ode))
print("--- Time to Big Rip ---")
print("w  t_rip-t_now(Gyr)  t_rip(Gyr)")
w_phantom = [-1.05, -1.1, -1.2, -1.3, -1.5, -2.0]
for w in w_phantom:
    dt_rip = 2.0 / (3.0 * abs(1+w) * H0_Gyr * np.sqrt(Ode))
    print(f"{w:.2f}   {dt_rip:.1f}            {t_now + dt_rip:.1f}")

print()

# Detailed timeline for w = -1.5
w = -1.5
dt_rip = 2.0 / (3.0 * abs(1+w) * H0_Gyr * np.sqrt(Ode))
t_rip = t_now + dt_rip
print(f"Detailed Big Rip timeline for w = {w}")
print(f"Big Rip occurs at t = {t_rip:.2f} Gyr ({dt_rip:.2f} Gyr from now)")
print()

# Scale factor evolution toward the rip
# a(t) ~ (t_rip - t)^(2/(3(1+w)))
exp = 2.0 / (3.0 * (1 + w))  # negative for w < -1

print("t_before_rip(Gyr) a/a_now H(t)/H0")
times_before = [dt_rip, dt_rip*0.8, dt_rip*0.5, dt_rip*0.2, dt_rip*0.1,
                1.0, 0.1, 0.01, 0.001, 1e-4, 1e-6, 1e-9, 1e-12]
for dt_b in times_before:
    if dt_b > 0:
        a_ratio = (dt_b / dt_rip)**exp
        H_ratio = abs(exp) / dt_b / H0_Gyr if dt_b > 1e-20 else float('inf')
        if a_ratio < 1e15:
            print(f"{dt_b:.2e}  {a_ratio:.4e}  {H_ratio:.4e}")
        else:
            print(f"{dt_b:.2e}  -> infinity  -> infinity")

print()

# Structure disruption times
print("=== Structure Disruption Timeline (w = -1.5) ===")
print()
# A bound structure of size R and binding force F is disrupted when")
# the "rip acceleration" exceeds the binding: H^2 * R > F/m
# This happens at t_rip - t ~ (binding timescale)

structures = [
    ("Galaxy clusters", 60e6, "yr"),
    ("Milky Way-size galaxy", 60e6, "yr"),
    ("Solar system (Sun-Earth)", 0.25, "yr"),
    ("Earth destroyed", 0.5/12, "yr"),
    ("Atoms dissociated", 1e-19/(3.156e7), "yr"),
]

print(f"Time before Big Rip | Structure disrupted")
print("-" * 50)
for name, dt_yr, _ in structures:
    if dt_yr > 1e6:
        print(f"~{dt_yr/1e6:.0f} Myr before     | {name}")
    elif dt_yr > 1:
        print(f"~{dt_yr:.1f} yr before      | {name}")
    elif dt_yr > 1/12:
        print(f"~{dt_yr*12:.0f} months before  | {name}")
    else:
        print(f"~{dt_yr*3.156e7:.0e} s before  | {name}")

print()
print("w t_rip_from_now_Gyr")
for w in [-1.05, -1.1, -1.15, -1.2, -1.3, -1.4, -1.5, -1.6, -1.8, -2.0]:
    dt_rip = 2.0 / (3.0 * abs(1+w) * H0_Gyr * np.sqrt(Ode))
    print(f"{w:.2f} {dt_rip:.2f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

7.6 CPL Dark Energy: $w_0$-$w_a$ Parametrization

Explore the Chevallier-Polarski-Linder parametrization $w(a) = w_0 + w_a(1-a)$and see how time-varying dark energy affects the expansion history.

CPL w(z) and H(z) for Time-Varying Dark Energy

Python

script.py78 lines

import numpy as np

H0 = 67.4
Om = 0.315
Ode = 1.0 - Om
c = 3.0e5

def w_cpl(z, w0, wa):
    """CPL equation of state"""
    a = 1.0 / (1.0 + z)
    return w0 + wa * (1.0 - a)

def rho_de_cpl(z, w0, wa):
    """DE density ratio rho_DE(z)/rho_DE(0) for CPL"""
    a = 1.0 / (1.0 + z)
    return a**(-3*(1+w0+wa)) * np.exp(-3*wa*(1-a))

def H_cpl(z, w0, wa):
    """H(z)/H0 for flat FLRW with CPL dark energy"""
    return np.sqrt(Om*(1+z)**3 + Ode * rho_de_cpl(z, w0, wa))

# Parameter sets to compare
params = [
    (-1.0, 0.0, "LambdaCDM"),
    (-0.55, -1.32, "DESI best-fit"),
    (-0.8, -0.5, "Thawing"),
    (-1.1, 0.3, "Freezing"),
]

print("=== CPL Dark Energy: w(z) and H(z) ===")
print()

# w(z) evolution
print("z w_LCDM w_DESI w_Thaw w_Freeze")
z_vals = np.linspace(0, 3, 50)
for z in z_vals:
    ws = [w_cpl(z, w0, wa) for w0, wa, _ in params]
    print(f"{z:.3f} " + " ".join(f"{w:.4f}" for w in ws))

print()

# H(z)/H0
print("z H_LCDM H_DESI H_Thaw H_Freeze")
for z in z_vals:
    Hs = [H_cpl(z, w0, wa) for w0, wa, _ in params]
    print(f"{z:.3f} " + " ".join(f"{h:.4f}" for h in Hs))

print()

# Luminosity distance comparison
def dL_cpl(z, w0, wa, n=300):
    zs = np.linspace(0, z, n+1)
    dz = z / n
    integral = 0.0
    for i in range(n):
        integral += 0.5 * (1.0/H_cpl(zs[i], w0, wa) + 1.0/H_cpl(zs[i+1], w0, wa)) * dz
    return (c/H0) * (1+z) * integral

print("--- Distance Modulus Differences from LambdaCDM ---")
print("z delta_mu_DESI delta_mu_Thaw delta_mu_Freeze")
for z in [0.1, 0.3, 0.5, 0.7, 1.0, 1.5, 2.0]:
    dL_ref = dL_cpl(z, -1.0, 0.0)
    mu_ref = 5*np.log10(max(dL_ref,1)) + 25
    for w0, wa, label in params[1:]:
        dL_val = dL_cpl(z, w0, wa)
        mu_val = 5*np.log10(max(dL_val,1)) + 25
    deltas = []
    for w0, wa, label in params[1:]:
        dL_val = dL_cpl(z, w0, wa)
        mu_val = 5*np.log10(max(dL_val,1)) + 25
        deltas.append(mu_val - mu_ref)
    print(f"{z:.1f}  " + "  ".join(f"{d:+.4f}" for d in deltas))

print()
print("The DESI best-fit (w0=-0.55, wa=-1.32) implies dark energy was")
print("weaker in the past (w closer to 0) and is strengthening over time.")
print("This is a 'thawing' behavior: DE was frozen and is now evolving.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Dark Energy

1. Introduction & Observational Evidence

The Acceleration Condition

Key Observational Milestones

The Cosmic Energy Budget

2. The Friedmann Equations

The FLRW Metric

Derivation from Einstein’s Equations

The Fluid Equation and Density Parameters

Hubble Parameter as a Function of Redshift

3. The Cosmological Constant

$\Lambda$ as Vacuum Energy

The Cosmological Constant Problem

The Coincidence Problem

4. Dark Energy Models Beyond $\Lambda$

4.1 Quintessence: Scalar Field Dark Energy

4.2 Phantom Energy and the Big Rip

4.3 k-essence

4.4 The CPL Parametrization

4.5 Modified Gravity: f(R) Theories

5. Observational Probes of Dark Energy

5.1 Distance-Redshift Relations

5.2 Type Ia Supernovae

5.3 Baryon Acoustic Oscillations (BAO)

5.4 CMB Power Spectrum

5.5 Weak Lensing and Cluster Counts

6. The Fate of the Universe

Scenario 1: Big Freeze / Heat Death ($w = -1$)

Scenario 2: Big Rip ($w < -1$)

Scenario 3: Big Crunch ($w > -1/3$ or Dark Energy Decaying)

The Deceleration Parameter

7. Python Simulations

7.1 Scale Factor Evolution for Different Dark Energy Models

Scale Factor a(t) for Different w Values

7.2 Distance-Redshift Relation

Luminosity Distance and Hubble Diagram

7.3 Deceleration Parameter and Cosmic Acceleration

Deceleration Parameter q(z)

7.4 Energy Density Evolution and Dark Energy Domination

Cosmic Energy Density Evolution

7.5 Big Rip Timescale for Phantom Dark Energy

Big Rip Timeline

7.6 CPL Dark Energy: $w_0$-$w_a$ Parametrization

CPL w(z) and H(z) for Time-Varying Dark Energy

Prerequisites

General Relativity

Cosmology

Quantum Field Theory

Mathematics

Further Reading

Textbooks

Key Review Articles

Landmark Observational Papers

Related Courses

Dark Matter

Cosmic Inflation

Observational Cosmology

Gravitational Waves

Black Holes

Quantum Gravity