← Back to Part III: Einstein Field Equations

Chapter 14: Stress-Energy Tensor

The stress-energy tensor T_μν describes the distribution of energy, momentum, and stress in spacetime. It's the source of gravity in Einstein's equations, replacing Newtonian mass density with a full description of matter and fields.

Physical Interpretation

T⁰⁰ = Energy Density

$\rho c^2$ — mass-energy per unit volume

T⁰ⁱ = Momentum Density

Energy flux / momentum per unit volume

Tⁱ⁰ = Energy Flux

Energy flow per unit area per unit time

T^ij = Stress Tensor

Pressure (diagonal), shear stress (off-diagonal)

Key property: T_μν is symmetric (T_μν = T_νμ), ensuring angular momentum conservation.

Perfect Fluid

The most common matter model in cosmology and astrophysics:

$T_{\mu\nu} = (\rho + p/c^2) u_\mu u_\nu + p \, g_{\mu\nu}$

$\rho$ = energy density, $p$ = pressure, $u^\mu$ = 4-velocity

Dust (p = 0)

$T_{\mu\nu} = \rho u_\mu u_\nu$

Non-relativistic matter, galaxies

Radiation (p = ρc²/3)

Equation of state w = 1/3

Photons, early universe

Dark Energy (p = -ρc²)

$T_{\mu\nu} = -\rho c^2 g_{\mu\nu}$ (cosmological constant)

Accelerating expansion

Stiff Matter (p = ρc²)

Equation of state w = 1

Hypothetical, very early universe

Electromagnetic Field

$T_{\mu\nu} = \frac{1}{\mu_0}\left( F_{\mu\alpha}F_\nu^{\;\alpha} - \frac{1}{4}g_{\mu\nu}F_{\alpha\beta}F^{\alpha\beta} \right)$

Properties

• Trace-free: T^μ_μ = 0 (radiation-like)
• Energy density: $u = \frac{1}{2}(\epsilon_0 E^2 + B^2/\mu_0)$
• Poynting vector: $\vec{S} = \vec{E} \times \vec{B}/\mu_0$
• Maxwell stress tensor for T^ij

Conservation Law

$\nabla_\mu T^{\mu\nu} = 0$

Covariant conservation of energy-momentum

This equation generalizes energy and momentum conservation to curved spacetime. It contains:

ν = 0: Energy Conservation

Rate of change of energy density + divergence of energy flux = 0

ν = i: Momentum Conservation

Euler/Navier-Stokes equations emerge from spatial components

Python: Perfect Fluid Evolution

Python

stress_energy.py103 lines

#!/usr/bin/env python3
"""
stress_energy.py - Stress-energy tensor for different matter types
"""
import numpy as np

def perfect_fluid_stress_energy(rho, p, u_mu, g_mu_nu):
    """
    Compute perfect fluid stress-energy tensor
    T_mu_nu = (rho + p) u_mu u_nu + p g_mu_nu
    """
    T = np.zeros((4, 4))
    for mu in range(4):
        for nu in range(4):
            T[mu, nu] = (rho + p) * u_mu[mu] * u_mu[nu] + p * g_mu_nu[mu, nu]
    return T

def em_stress_energy(E, B, c=1):
    """
    Electromagnetic stress-energy tensor (3+1 form)
    """
    eps0 = 1  # Natural units
    mu0 = 1

T = np.zeros((4, 4))

# Energy density
    u = 0.5 * (eps0 * np.dot(E, E) + np.dot(B, B) / mu0)
    T[0, 0] = u

# Poynting vector (momentum density / energy flux)
    S = np.cross(E, B) / mu0
    T[0, 1:] = S / c
    T[1:, 0] = S / c

# Maxwell stress tensor
    for i in range(3):
        for j in range(3):
            T[i+1, j+1] = (eps0 * (E[i]*E[j] - 0.5*(i==j)*np.dot(E,E)) +
                          (B[i]*B[j] - 0.5*(i==j)*np.dot(B,B)) / mu0)

return T

# Flat spacetime metric
eta = np.diag([-1, 1, 1, 1])

# Fluid at rest: u^mu = (1, 0, 0, 0)
# Lower index: u_mu = g_mu_nu u^nu = (-1, 0, 0, 0)
u_mu = np.array([-1, 0, 0, 0])

print("Stress-Energy Tensors in General Relativity")
print("="*60)

# Dust (non-relativistic matter)
print("\n1. DUST (pressureless matter, p = 0)")
rho_dust = 1.0  # Energy density
p_dust = 0.0
T_dust = perfect_fluid_stress_energy(rho_dust, p_dust, u_mu, eta)
print("T_mu_nu =")
print(T_dust)
print(f"Trace T = {np.einsum('ii', np.einsum('ij,jk->ik', np.linalg.inv(eta), T_dust)):.4f}")

# Radiation
print("\n2. RADIATION (p = rho/3)")
rho_rad = 1.0
p_rad = rho_rad / 3
T_rad = perfect_fluid_stress_energy(rho_rad, p_rad, u_mu, eta)
print("T_mu_nu =")
print(np.round(T_rad, 4))
print(f"Trace T = {np.einsum('ii', np.einsum('ij,jk->ik', np.linalg.inv(eta), T_rad)):.4f}")
print("(Trace = 0 for radiation!)")

# Dark energy / Cosmological constant
print("\n3. DARK ENERGY / COSMOLOGICAL CONSTANT (p = -rho)")
rho_de = 1.0
p_de = -rho_de
T_de = perfect_fluid_stress_energy(rho_de, p_de, u_mu, eta)
print("T_mu_nu =")
print(T_de)
print("(Proportional to metric: T_mu_nu = -rho g_mu_nu)")

# Electromagnetic field
print("\n4. ELECTROMAGNETIC FIELD")
E = np.array([1.0, 0, 0])  # Electric field in x-direction
B = np.array([0, 1.0, 0])  # Magnetic field in y-direction
T_em = em_stress_energy(E, B)
print(f"E = {E}")
print(f"B = {B}")
print("T_mu_nu =")
print(np.round(T_em, 4))
print(f"Trace T = {np.trace(np.einsum('ij,jk->ik', np.linalg.inv(eta), T_em)):.4f}")
print("(Trace = 0 for EM field!)")

# Equation of state
print("\n" + "="*60)
print("Equation of State: p = w rho c^2")
print("-"*60)
print("w = 0    : Non-relativistic matter (dust)")
print("w = 1/3  : Radiation (photons, relativistic particles)")
print("w = -1   : Cosmological constant (dark energy)")
print("w = -1/3 : Cosmic strings")
print("w = 1    : Stiff matter (speculative)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran: Conservation Check

Python

conservation_plot.py97 lines

import subprocess, tempfile, os
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

fortran_source = r"""
program conservation_law
    implicit none
    integer, parameter :: dp = selected_real_kind(15)
    real(dp) :: a, rho0, w
    integer :: i, n_pts, j
    real(dp) :: w_values(4)

w_values(1) = 0.0_dp
    w_values(2) = 1.0_dp/3.0_dp
    w_values(3) = -1.0_dp
    w_values(4) = -1.0_dp/3.0_dp

rho0 = 1.0_dp
    n_pts = 200

do j = 1, 4
        w = w_values(j)
        write(*,'(A,F8.4)') 'DATA_START w=', w
        do i = 1, n_pts
            a = 0.1_dp + real(i-1, dp) * 9.9_dp / real(n_pts-1, dp)
            write(*,'(F14.8,A,ES16.8)') a, ',', rho0 * a**(-3.0_dp*(1.0_dp + w))
        end do
        print '(A)', 'DATA_END'
    end do
end program conservation_law
"""

workdir = tempfile.mkdtemp()
src = os.path.join(workdir, 'cons.f90')
exe = os.path.join(workdir, 'cons')

with open(src, 'w') as f:
    f.write(fortran_source)

comp = subprocess.run(['gfortran', '-O2', '-o', exe, src],
                      capture_output=True, text=True, timeout=30)
if comp.returncode != 0:
    print("Compilation error:", comp.stderr)
else:
    result = subprocess.run([exe], capture_output=True, text=True, timeout=30, cwd=workdir)
    output = result.stdout

sections = {}
    current = None
    for line in output.split('\n'):
        if 'DATA_START w=' in line:
            current = line.split('w=')[1].strip()
            sections[current] = []
            continue
        if 'DATA_END' in line:
            current = None
            continue
        if current and line.strip():
            sections[current].append(line.strip())

fig, ax = plt.subplots(1, 1, figsize=(10, 7))
    fig.patch.set_facecolor('#0f172a')
    ax.set_facecolor('#1e293b')

labels_map = {
        '0.0000': ('Dust (w=0): $\\rho \\propto a^{-3}$', '#22d3ee'),
        '0.3333': ('Radiation (w=1/3): $\\rho \\propto a^{-4}$', '#f472b6'),
        '-1.0000': ('Dark Energy (w=-1): $\\rho = const$', '#4ade80'),
        '-0.3333': ('Curvature (w=-1/3): $\\rho \\propto a^{-2}$', '#fbbf24'),
    }

for w_key, data in sections.items():
        a_arr, rho_arr = [], []
        for line in data:
            vals = [float(x) for x in line.split(',')]
            a_arr.append(vals[0])
            rho_arr.append(vals[1])
        label, color = labels_map.get(w_key, (f'w={w_key}', '#94a3b8'))
        ax.semilogy(a_arr, rho_arr, color=color, linewidth=2.5, label=label)

ax.axvline(x=1.0, color='#94a3b8', linestyle=':', alpha=0.5, label='Present (a=1)')
    ax.set_xlabel('Scale factor a', color='#cbd5e1', fontsize=12)
    ax.set_ylabel('Energy density $\\rho / \\rho_0$', color='#cbd5e1', fontsize=12)
    ax.set_title('$\\nabla_\\mu T^{\\mu\\nu} = 0$: Energy Density vs Scale Factor', color='white', fontsize=14, fontweight='bold')
    ax.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1', fontsize=10, loc='upper right')
    ax.tick_params(colors='#cbd5e1')
    for spine in ax.spines.values():
        spine.set_color('#334155')
    ax.grid(True, color='#334155', alpha=0.4)

plt.tight_layout()
    plt.savefig('conservation_plot.png', dpi=150, bbox_inches='tight', facecolor=fig.get_facecolor())
    plt.close()
    print("Conservation: rho ~ a^(-3(1+w)) for each equation of state")