Symbolic Tensor Algebra

from itensorpy import Metric # Schwarzschild metric (non-rotating black hole) schwarzschild = Metric.schwarzschild(M=1) # M is mass parameter # Kerr metric (rotating black hole) kerr = Metric.kerr(M=1, a=0.9) # M is mass, a is spin parameter # FLRW metric (expanding universe) flrw = Metric.friedmann_lemaitre_robertson_walker(k=0) # k is curvature parameter

from itensorpy import Metric, symbols from sympy import sin, exp # Define coordinate symbols t, r, theta, phi = symbols('t r theta phi') # Define metric components for a simple wormhole metric alpha = symbols('alpha', positive=True) # Throat parameter g_tt = -exp(2*alpha) g_rr = 1 g_thth = r**2 + alpha**2 g_phph = (r**2 + alpha**2) * sin(theta)**2 # Create metric from components (specify coordinate order and components) wormhole = Metric( coords=[t, r, theta, phi], components=[ [g_tt, 0, 0, 0], [0, g_rr, 0, 0], [0, 0, g_thth, 0], [0, 0, 0, g_phph] ], name="Ellis-Bronnikov Wormhole" )

# Computing Christoffel symbols for the Schwarzschild metric christoffel = schwarzschild.compute_christoffel_symbols() # Display all non-zero components print(christoffel.display_components(only_nonzero=True)) # Display a specific component (using index notation) print(christoffel[t, r, t]) # Γ^t_rt component # Show derivation steps christoffel_with_steps = schwarzschild.compute_christoffel_symbols(show_steps=True) print(christoffel_with_steps.derivation)

# Computing the Riemann curvature tensor riemann = schwarzschild.compute_riemann_tensor() # Display all non-zero components print(riemann.display_components(only_nonzero=True)) # Computing a specific component print("R^t_rtr = ", riemann[t, r, t, r].simplify())

# Computing the Ricci tensor (contraction of Riemann tensor) ricci = schwarzschild.compute_ricci_tensor() print(ricci.display_components()) # Computing the Ricci scalar (contraction of Ricci tensor) ricci_scalar = schwarzschild.compute_scalar_curvature() print("Ricci scalar: ", ricci_scalar.simplify()) # For Schwarzschild (a vacuum solution), these should all be zero # This can be verified numerically at a specific point print("Max component magnitude: ", ricci.max_abs_component())

# Computing the Einstein tensor einstein = schwarzschild.compute_einstein_tensor() print(einstein.display_components()) # For non-vacuum spacetimes, the Einstein tensor equals 8πG times # the stress-energy tensor, per Einstein's field equations: # G_μν = 8πG T_μν # Example with a simple dust matter distribution from itensorpy import StressEnergyTensor rho = symbols('rho') # Matter density dust = StressEnergyTensor.dust(density=rho) G = symbols('G') # Gravitational constant # Einstein's field equations field_equations = einstein - 8*pi*G*dust print(field_equations.display_components())

from itensorpy import Tensor, symbols, indices # Define coordinates and base tensors x, y, z = symbols('x y z') # Define arbitrary tensors A = Tensor('A', rank=(1, 1)) # Type (1,1) tensor - one upper, one lower index B = Tensor('B', rank=(1, 2)) # Type (1,2) tensor # Define indices for operations i, j, k, l = indices('i j k l') # Tensor addition (tensors must have same index structure) C1 = A + A print(C1.display()) # 2*A^i_j # Tensor multiplication (creates tensor with combined indices) C2 = A * B print(C2.display()) # A^i_j * B^k_lm # Tensor contraction (sum over repeated indices) C3 = A[i, j] * B[j, k, l] print(C3.display()) # A^i_j * B^j_kl # Lower/raise indices with the metric g = Metric.minkowski() # Flat spacetime metric A_lower = A.lower_index(0, g) # Lower first index using metric g print(A_lower.display()) # g_ij * A^j_k # Symmetrization C4 = A[i, j] * B[k, i, l] C4_sym = C4.symmetrize([0, 2]) # Symmetrize over 1st and 3rd indices print(C4_sym.display())

from itensorpy import symbols from sympy import solve # Create Schwarzschild metric schwarzschild = Metric.schwarzschild(M=1) # Get geodesic equations geodesic_eqs = schwarzschild.geodesic_equations() # Equations are parameterized by an affine parameter λ # For each coordinate x^μ, we get a second-order ODE: d²x^μ/dλ² + Γ^μ_ρσ * (dx^ρ/dλ) * (dx^σ/dλ) = 0 # For circular orbits at r=constant, we need dr/dλ = 0 and d²r/dλ² = 0 # Let's solve for the orbital velocity (dφ/dt) for a circular orbit: # First, get the r component of the geodesic equation r_eq = geodesic_eqs[1] # 0=t, 1=r, 2=θ, 3=φ # For a circular orbit in the equatorial plane: # θ = π/2, dθ/dλ = 0, r = constant, dr/dλ = 0 # Substitute these conditions from sympy import symbols, pi r_val = symbols('r_val', positive=True) dt_dlambda = symbols('dt_dlambda', positive=True) dphi_dlambda = symbols('dphi_dlambda') substitutions = { 'r': r_val, 'theta': pi/2, 'd_r_dlambda': 0, 'd_theta_dlambda': 0, 'd_t_dlambda': dt_dlambda, 'd_phi_dlambda': dphi_dlambda } # Apply substitutions to the r geodesic equation r_eq_circular = r_eq.subs(substitutions) # Solve for the ratio (dphi_dlambda/dt_dlambda)² which is (orbital angular velocity)² omega_squared = solve(r_eq_circular, (dphi_dlambda/dt_dlambda)**2)[0] print("Orbital angular velocity squared: ", omega_squared) # Result should be: M/r³ (Kepler's third law in GR for small r) # The exact formula is: Ω² = M/r³/(1 - 3M/r)

from itensorpy import Metric, symbols from sympy import sin, cos # Define coordinates t, x, y, z = symbols('t x y z') # Define wave parameters omega, k = symbols('omega k', positive=True) h_plus = symbols('h_plus', positive=True) # Wave amplitude for + polarization # Define a linearized metric with a gravitational wave perturbation # We use η_μν + h_μν where η is Minkowski metric and h is the small perturbation # For a wave traveling in z direction with + polarization: # h_xx = -h_yy = h_plus * cos(omega*t - k*z) # Create the perturbed metric eta = Metric.minkowski() # Flat spacetime wave_phase = omega*t - k*z h_xx = h_plus * cos(wave_phase) h_yy = -h_plus * cos(wave_phase) # Start with Minkowski components g_components = [ [-1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1] ] # Add the perturbation g_components[1][1] += h_xx # Add h_xx to g_xx g_components[2][2] += h_yy # Add h_yy to g_yy # Create the metric g_wave = Metric([t, x, y, z], g_components, name="Gravitational Wave") # Compute the Ricci tensor to first order in h_plus ricci = g_wave.compute_ricci_tensor(simplify=True) # In vacuum, Einstein's equations require R_μν = 0 # This gives us the wave equation: □h_μν = 0 # where □ is the d'Alembertian operator # Extract the wave equation from the Ricci tensor wave_eq_xx = ricci[1, 1].simplify() print("Wave equation (xx component): ", wave_eq_xx) # Verify it matches the expected form: (∂²/∂t² - ∂²/∂z²)h_xx = 0 # This is the wave equation with propagation speed c=1 (speed of light) # Should yield: -omega² + k² = 0, which means omega = k (dispersion relation)