Chapter 9: Weyl Tensor

The Weyl tensor is the trace-free part of the Riemann tensor. It describes the purely gravitational degrees of freedom—tidal forces and gravitational waves—independent of local matter content.

Definition

The Weyl tensor (conformal tensor) is defined as:

$C_{\rho\sigma\mu\nu} = R_{\rho\sigma\mu\nu} - \frac{2}{n-2}\left(g_{\rho[\mu}R_{\nu]\sigma} - g_{\sigma[\mu}R_{\nu]\rho}\right) + \frac{2}{(n-1)(n-2)} R \, g_{\rho[\mu}g_{\nu]\sigma}$

In 4 dimensions, this simplifies to:

$C_{\rho\sigma\mu\nu} = R_{\rho\sigma\mu\nu} - \left(g_{\rho[\mu}R_{\nu]\sigma} - g_{\sigma[\mu}R_{\nu]\rho}\right) + \frac{1}{3} R \, g_{\rho[\mu}g_{\nu]\sigma}$

Key Properties

Trace-Free

$C^\rho_{\;\sigma\rho\nu} = 0$ — All contractions vanish

Same Symmetries as Riemann

Antisymmetric in first pair, antisymmetric in second pair, symmetric under pair exchange

Conformally Invariant

$Under (\tilde{g}_{\mu\nu} = \Omega^2 g_{\mu\nu}): (\tilde{C}^\rho_{\;\sigma\mu\nu} = C^\rho_{\;\sigma\mu\nu})$

Independent Components

10 components in 4D (same as Ricci tensor) — total 20 = 10 (Weyl) + 10 (Ricci)

Physical Interpretation

Gravitational Waves

In vacuum (R_μν = 0), all curvature is in Weyl tensor. Gravitational waves are propagating Weyl curvature.

Tidal Effects

The Weyl tensor describes shape distortion without volume change (pure shear).

Free Gravitational Field

Represents curvature not directly caused by local matter—propagating from distant sources.

Petrov Classification

Spacetimes classified by algebraic type of Weyl tensor (Type I, II, D, III, N, O).

Python: Weyl Tensor Computation

Python

weyl_tensor.py212 lines

#!/usr/bin/env python3
"""
╔══════════════════════════════════════════════════════════════════════════════╗
║                        Weyl Tensor Computation                               ║
║                     Schwarzschild Spacetime (Numerical)                      ║
╚══════════════════════════════════════════════════════════════════════════════╝

Weyl tensor in 4D:
  C_rsmn = R_rsmn - (g_r[m R_n]s - g_s[m R_n]r) + (1/3) R g_r[m g_n]s

In vacuum (R_mn = 0):  C_rsmn = R_rsmn
Kretschmann scalar: C_abcd C^abcd = 48 M^2 / r^6
"""
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import io, base64

# ═══════════════════════════════════════════════════════════════
# Numerical helpers
# ═══════════════════════════════════════════════════════════════

def schwarzschild_metric(r_val, theta_val, M=1.0):
    f = 1.0 - 2.0*M/r_val
    g = np.diag([-f, 1.0/f, r_val**2, r_val**2 * np.sin(theta_val)**2])
    g_inv = np.diag([-1.0/f, f, 1.0/r_val**2, 1.0/(r_val**2 * np.sin(theta_val)**2)])
    return g, g_inv

def metric_derivs(r_val, theta_val, M=1.0, h=1e-6):
    n = 4
    dg = [np.zeros((n, n)) for _ in range(n)]
    gp, _ = schwarzschild_metric(r_val + h, theta_val, M)
    gm, _ = schwarzschild_metric(r_val - h, theta_val, M)
    dg[1] = (gp - gm) / (2*h)
    gp, _ = schwarzschild_metric(r_val, theta_val + h, M)
    gm, _ = schwarzschild_metric(r_val, theta_val - h, M)
    dg[2] = (gp - gm) / (2*h)
    return dg

def compute_christoffel(r_val, theta_val, M=1.0, h=1e-6):
    n = 4
    g, g_inv = schwarzschild_metric(r_val, theta_val, M)
    dg = metric_derivs(r_val, theta_val, M, h)
    Gamma = np.zeros((n, n, n))
    for rho in range(n):
        for mu in range(n):
            for nu in range(n):
                s = 0.0
                for sigma in range(n):
                    s += g_inv[rho, sigma] * (dg[mu][sigma, nu] + dg[nu][sigma, mu] - dg[sigma][mu, nu])
                Gamma[rho, mu, nu] = 0.5 * s
    return Gamma

def riemann_tensor(r_val, theta_val, M=1.0, h=1e-5):
    """Compute R^rho_sigma_mu_nu numerically"""
    n = 4
    G0 = compute_christoffel(r_val, theta_val, M)
    dGamma = [np.zeros((n, n, n)) for _ in range(n)]
    Gp = compute_christoffel(r_val + h, theta_val, M)
    Gm = compute_christoffel(r_val - h, theta_val, M)
    dGamma[1] = (Gp - Gm) / (2*h)
    Gp = compute_christoffel(r_val, theta_val + h, M)
    Gm = compute_christoffel(r_val, theta_val - h, M)
    dGamma[2] = (Gp - Gm) / (2*h)

R = np.zeros((n, n, n, n))
    for rho in range(n):
        for sigma in range(n):
            for mu in range(n):
                for nu in range(n):
                    val = dGamma[mu][rho, nu, sigma] - dGamma[nu][rho, mu, sigma]
                    for lam in range(n):
                        val += G0[rho, mu, lam]*G0[lam, nu, sigma] - G0[rho, nu, lam]*G0[lam, mu, sigma]
                    R[rho, sigma, mu, nu] = val
    return R, G0

def lower_riemann(R_up, g):
    return np.einsum('ra,asmn->rsmn', g, R_up)

def ricci_from_riemann(R_up):
    return np.einsum('rmrn->mn', R_up)

def weyl_tensor(R_lower, Ricci, R_scalar, g):
    n = 4
    C = np.zeros((n, n, n, n))
    for r in range(n):
        for s in range(n):
            for m in range(n):
                for nu in range(n):
                    t1 = R_lower[r, s, m, nu]
                    t2 = -(g[r, m]*Ricci[nu, s] - g[r, nu]*Ricci[m, s]
                           - g[s, m]*Ricci[nu, r] + g[s, nu]*Ricci[m, r]) / 2.0
                    t3 = R_scalar * (g[r, m]*g[nu, s] - g[r, nu]*g[m, s]) / 6.0
                    C[r, s, m, nu] = t1 + t2 + t3
    return C

# ═══════════════════════════════════════════════════════════════
# Evaluate at r = 6M, theta = pi/2
# ═══════════════════════════════════════════════════════════════
M = 1.0
r0, th0 = 6.0*M, np.pi/2
g, g_inv = schwarzschild_metric(r0, th0, M)
R_up, Gamma0 = riemann_tensor(r0, th0, M)
R_low = lower_riemann(R_up, g)
Ricci = ricci_from_riemann(R_up)
R_scalar = np.einsum('mn,mn->', g_inv, Ricci)
C = weyl_tensor(R_low, Ricci, R_scalar, g)

print("Weyl Tensor for Schwarzschild Spacetime (Numerical)")
print("=" * 60)
print(f"\nEvaluated at r = {r0:.1f}M, theta = pi/2")
print(f"\nRicci scalar R = {R_scalar:.2e}  (zero for vacuum)")

coord_names = ['t', 'r', 'theta', 'phi']
print("\nNon-trivial Weyl components C_rsmn (|val| > 1e-10):")
count = 0
for r in range(4):
    for s in range(r+1, 4):
        for m in range(4):
            for nu in range(m+1, 4):
                val = C[r, s, m, nu]
                if abs(val) > 1e-10:
                    print(f"  C_{coord_names[r]}{coord_names[s]}{coord_names[m]}{coord_names[nu]} = {val:+.8e}")
                    count += 1
print(f"  ({count} non-zero independent components)")

# Verify Weyl = Riemann in vacuum
diff_CR = np.max(np.abs(C - R_low))
print(f"\nMax |C_rsmn - R_rsmn| = {diff_CR:.2e}")
print("-> In vacuum, Weyl tensor equals the full Riemann tensor!")

# Kretschmann scalar
R_up4 = np.einsum('ra,sb,rsmn->abmn', g_inv, g_inv, R_low)
R_fully_up = np.einsum('mc,nd,abmn->abcd', g_inv, g_inv, R_up4)
K_numerical = np.einsum('abcd,abcd->', R_low, R_fully_up)
K_analytic = 48.0 * M**2 / r0**6
print(f"\nKretschmann scalar K = {K_numerical:.8e}")
print(f"Analytic 48M^2/r^6   = {K_analytic:.8e}")
print(f"Relative error        = {abs(K_numerical - K_analytic)/K_analytic:.2e}")

# ═══════════════════════════════════════════════════════════════
# Plot Weyl components vs radius
# ═══════════════════════════════════════════════════════════════
r_arr = np.linspace(2.5*M, 20*M, 150)
C_trtr = np.zeros_like(r_arr)
C_thph_thph = np.zeros_like(r_arr)
K_arr = np.zeros_like(r_arr)
K_exact = 48.0 * M**2 / r_arr**6

for i, rv in enumerate(r_arr):
    gi, gi_inv = schwarzschild_metric(rv, th0, M)
    Ri_up, _ = riemann_tensor(rv, th0, M)
    Ri_low = lower_riemann(Ri_up, gi)
    Ric_i = ricci_from_riemann(Ri_up)
    Rs_i = np.einsum('mn,mn->', gi_inv, Ric_i)
    Ci = weyl_tensor(Ri_low, Ric_i, Rs_i, gi)
    C_trtr[i] = Ci[0, 1, 0, 1]
    C_thph_thph[i] = Ci[2, 3, 2, 3]
    Ri4 = np.einsum('ra,sb,rsmn->abmn', gi_inv, gi_inv, Ri_low)
    Ri_fu = np.einsum('mc,nd,abmn->abcd', gi_inv, gi_inv, Ri4)
    K_arr[i] = np.einsum('abcd,abcd->', Ri_low, Ri_fu)

fig, axes = plt.subplots(1, 3, figsize=(18, 5))
fig.patch.set_facecolor('#0a0a1a')
fig.suptitle('Weyl Tensor & Kretschmann Scalar for Schwarzschild',
             color='white', fontsize=14, fontweight='bold', y=1.02)

ax = axes[0]
ax.set_facecolor('#0a0a1a')
ax.plot(r_arr/M, C_trtr, color='#f472b6', linewidth=2, label=r'$C_{trtr}$')
ax.plot(r_arr/M, C_thph_thph, color='#38bdf8', linewidth=2, label=r'$C_{\theta\phi\theta\phi}$')
ax.axhline(0, color='#94a3b8', linewidth=0.5, linestyle='--')
ax.set_xlabel('r / M', color='#94a3b8')
ax.set_ylabel('Weyl component', color='#94a3b8')
ax.set_title('Weyl Tensor Components', color='white')
ax.legend(facecolor='#0a0a1a', edgecolor='#334155', labelcolor='#94a3b8')
ax.tick_params(colors='#94a3b8')
for spine in ax.spines.values():
    spine.set_color('#334155')

ax = axes[1]
ax.set_facecolor('#0a0a1a')
ax.semilogy(r_arr/M, K_arr, color='#c084fc', linewidth=2, label='Numerical')
ax.semilogy(r_arr/M, K_exact, color='#fbbf24', linewidth=1.5, linestyle='--', label=r'$48M^2/r^6$')
ax.axvline(x=2.0, color='#ef4444', linestyle=':', alpha=0.7, label='Horizon')
ax.set_xlabel('r / M', color='#94a3b8')
ax.set_ylabel('K', color='#94a3b8')
ax.set_title('Kretschmann Scalar (log)', color='white')
ax.legend(facecolor='#0a0a1a', edgecolor='#334155', labelcolor='#94a3b8', fontsize=9)
ax.tick_params(colors='#94a3b8')
for spine in ax.spines.values():
    spine.set_color('#334155')

ax = axes[2]
ax.set_facecolor('#0a0a1a')
rel_err = np.abs(K_arr - K_exact) / K_exact
ax.semilogy(r_arr/M, rel_err, color='#4ade80', linewidth=2)
ax.set_xlabel('r / M', color='#94a3b8')
ax.set_ylabel('Relative error', color='#94a3b8')
ax.set_title('Kretschmann: Numerical vs Analytic', color='white')
ax.tick_params(colors='#94a3b8')
for spine in ax.spines.values():
    spine.set_color('#334155')

plt.tight_layout()
buf = io.BytesIO()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a'))
# Image saved as output.png
buf.close()
plt.close()

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran: Weyl Scalar Calculation

Python

weyl_plot.py105 lines

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

fortran_source = r"""
program weyl_scalar
    implicit none
    integer, parameter :: dp = selected_real_kind(15)
    real(dp), parameter :: M = 1.0_dp
    real(dp) :: r, K_weyl, tidal_accel
    integer :: i, n_pts

n_pts = 200
    print '(A)', 'DATA_START'
    do i = 1, n_pts
        r = 2.01_dp + real(i-1, dp) * 47.99_dp / real(n_pts-1, dp)
        K_weyl = 48.0_dp * M**2 / r**6
        tidal_accel = sqrt(K_weyl)  ! ~ M/r^3
        write(*,'(F14.8,A,ES16.8,A,ES16.8)') r, ',', K_weyl, ',', tidal_accel
    end do
    print '(A)', 'DATA_END'
end program weyl_scalar
"""

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

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

lines = []
    in_data = False
    for line in output.split('\n'):
        if 'DATA_START' in line:
            in_data = True
            continue
        if 'DATA_END' in line:
            break
        if in_data and line.strip():
            lines.append(line.strip())

r_arr, K_arr, tidal = [], [], []
    for line in lines:
        vals = [float(x) for x in line.split(',')]
        r_arr.append(vals[0])
        K_arr.append(vals[1])
        tidal.append(vals[2])

r_arr = np.array(r_arr)
    K_arr = np.array(K_arr)
    tidal = np.array(tidal)

fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    fig.patch.set_facecolor('#0f172a')
    fig.suptitle('Weyl/Kretschmann Scalar: Tidal Forces in Schwarzschild', color='white', fontsize=14, fontweight='bold')

ax = axes[0]
    ax.set_facecolor('#1e293b')
    ax.semilogy(r_arr, K_arr, color='#f472b6', linewidth=2)
    ax.axvline(x=2.0, color='#ef4444', linestyle='--', alpha=0.7, label='Horizon (r=2M)')
    ax.axvline(x=6.0, color='#fbbf24', linestyle=':', alpha=0.7, label='ISCO (r=6M)')
    K_horizon = 48.0 / 2.0**6
    ax.axhline(y=K_horizon, color='#94a3b8', linestyle=':', alpha=0.5)
    ax.annotate(f'K(2M) = {K_horizon:.2f}\n(finite!)', xy=(2.5, K_horizon),
                color='#94a3b8', fontsize=9)
    ax.set_xlabel('r / M', color='#cbd5e1')
    ax.set_ylabel('$K = C_{\\alpha\\beta\\gamma\\delta}C^{\\alpha\\beta\\gamma\\delta}$', color='#cbd5e1')
    ax.set_title('Kretschmann Scalar (log scale)', color='white')
    ax.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')
    ax.tick_params(colors='#cbd5e1')
    for spine in ax.spines.values():
        spine.set_color('#334155')
    ax.grid(True, color='#334155', alpha=0.4)

ax2 = axes[1]
    ax2.set_facecolor('#1e293b')
    ax2.semilogy(r_arr, tidal, color='#4ade80', linewidth=2, label='$a_{tidal} \\sim \\sqrt{K}$')
    ax2.axvline(x=2.0, color='#ef4444', linestyle='--', alpha=0.7, label='Horizon')
    ax2.set_xlabel('r / M', color='#cbd5e1')
    ax2.set_ylabel('Tidal acceleration ($\\sim M/r^3$)', color='#cbd5e1')
    ax2.set_title('Tidal Force vs Radius', color='white')
    ax2.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')
    ax2.tick_params(colors='#cbd5e1')
    for spine in ax2.spines.values():
        spine.set_color('#334155')
    ax2.grid(True, color='#334155', alpha=0.4)

plt.tight_layout()
    buf = io.BytesIO()
    plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a'))
    buf.close()
    plt.close()
    # Image saved as output.png

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Ricci Tensor and Scalar Next: Bianchi Identities →