Chapter 14: Stress-Energy Tensor

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

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^μ = (1, 0, 0, 0)
# Lower index: u_μ = g_μν u^ν = (-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_μν =")
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 = ρ/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_μν =")
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_de = 1.0
p_de = -rho_de
T_de = perfect_fluid_stress_energy(rho_de, p_de, u_mu, eta)
print("T_μν =")
print(T_de)
print("(Proportional to metric: T_μν = -ρ g_μν)")

# 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_μν =")
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 ρ c²")
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)")