⭐

courseshub.world›Astrophysics & Computational Science›Quantum Gravity: Lattice & Information

Quantum Gravity: Lattice and Information Approaches

A comprehensive 8-chapter module exploring canonical quantum gravity, loop quantum gravity, spin foam models, the AdS/CFT correspondence, holographic entanglement entropy, tensor networks, the black hole information paradox, and quantum computing applications for quantum gravity.

8 Chapters24 Sections~45 Hours8 Code Labs

01Canonical Quantum GravityMini-Superspace WDW

02Loop Quantum GravityArea Spectrum

03Spin Foam Models6j Symbols & Foam Graph

04AdS/CFT CorrespondenceBTZ Entropy vs Cardy

05Holographic Entanglement EntropyRT Surface in AdS3

06Tensor NetworksMPS Entanglement

07Black Hole Information ParadoxPage Curve

08Quantum Computing for QGSpin Network Circuit

CH 01

Canonical Quantum Gravity

The canonical approach to quantum gravity begins by casting general relativity as a Hamiltonian system. The Arnowitt-Deser-Misner (ADM) formalism decomposes the four-dimensional spacetime metric into a family of spatial slices threaded together by a lapse function N and shift vector Nⁱ. The canonical momenta conjugate to the spatial metric γ_ij are related to the extrinsic curvature of each slice.

A remarkable feature of general relativity is that the Hamiltonian is a sum of constraints. There are no true dynamics in the conventional sense — instead, time evolution is generated by a constraint that must vanish on physical states. This is the essence of the “problem of time” in quantum gravity: the Hamiltonian annihilates physical states, so how do we recover a notion of temporal evolution?

The ADM Hamiltonian

In the ADM formalism, the total Hamiltonian is purely a sum of constraints weighted by the lapse and shift:

$$H = \int d^3x \left(N\mathcal{H} + N^i\mathcal{H}_i\right)$$

Here H is the Hamiltonian constraint (generating time reparametrizations) and H_i is the diffeomorphism constraint (generating spatial coordinate changes). The lapse N and shift Nⁱ are Lagrange multipliers enforcing these constraints. Physical states must satisfy both constraints simultaneously.

The Wheeler-DeWitt Equation

Quantization proceeds by promoting the constraints to operators acting on a wave functional Ψ[γ_ij] defined on the space of all three-metrics (superspace). The Hamiltonian constraint becomes the celebrated Wheeler-DeWitt equation:

$$\hat{\mathcal{H}}\Psi[\gamma_{ij}] = 0$$

This equation is the quantum gravity analogue of the Schrodinger equation, but with a crucial difference: there is no time derivative. The wave function of the universe does not evolve — it simply is. Time must be recovered internally, for example by using the scale factor of the universe as a clock variable in the mini-superspace approximation.

In the mini-superspace model, we freeze all degrees of freedom except the scale factor a and possibly a scalar field φ. This reduces the Wheeler-DeWitt equation to an ordinary differential equation that can be solved numerically, giving us quantum cosmological wave functions such as the Hartle-Hawking no-boundary state and the tunneling proposal of Vilenkin.

Lab: Toy Mini-Superspace WDW Equation

Mini-Superspace Wheeler-DeWitt Equation (1D Quantum Cosmology)

Python

script.py106 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

print("=" * 60)
print("  Mini-Superspace Wheeler-DeWitt Equation")
print("=" * 60)
print()

# Solve the 1D mini-superspace WDW equation:
#   -d^2 Psi/da^2 + V(a) Psi = 0
# where V(a) = a^2 - a^4 (closed FRW with cosmological constant)

N = 2000
a_max = 3.0
a = np.linspace(0.01, a_max, N)
da = a[1] - a[0]

# Potential: V(a) = Lambda * a^4 - k * a^2
# For a closed universe (k=1) with cosmological constant
Lambda_cosmo = 0.5
k = 1.0
V = Lambda_cosmo * a**4 - k * a**2

# Build the Hamiltonian matrix (finite differences)
H = np.zeros((N, N))
for i in range(1, N - 1):
    H[i, i - 1] = -1.0 / da**2
    H[i, i] = 2.0 / da**2 + V[i]
    H[i, i + 1] = -1.0 / da**2
H[0, 0] = 2.0 / da**2 + V[0]
H[0, 1] = -1.0 / da**2
H[-1, -2] = -1.0 / da**2
H[-1, -1] = 2.0 / da**2 + V[-1]

# Find the lowest eigenvalues/eigenvectors (WDW solutions)
from numpy.linalg import eigh
eigenvalues, eigenvectors = eigh(H)

# The WDW equation wants E=0 solutions; find states closest to E=0
idx_near_zero = np.argsort(np.abs(eigenvalues))[:4]
print(f"Eigenvalues closest to E=0:")
for i, idx in enumerate(idx_near_zero):
    print(f"  State {i}: E = {eigenvalues[idx]:.6f}")

# Plot
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.patch.set_facecolor('#0f172a')

# Panel 1: Potential
ax1 = axes[0, 0]
ax1.fill_between(a, V, alpha=0.3, color='#ec4899')
ax1.plot(a, V, color='#f472b6', linewidth=2)
ax1.axhline(y=0, color='#6b7280', linewidth=0.5, linestyle='--')
ax1.set_xlabel('Scale factor a', color='#d1d5db')
ax1.set_ylabel('V(a)', color='#d1d5db')
ax1.set_title('Mini-Superspace Potential', color='#f9a8d4', fontsize=13, fontweight='bold')
ax1.set_facecolor('#0f172a')
ax1.tick_params(colors='#9ca3af')
for spine in ax1.spines.values():
    spine.set_color('#334155')

# Panel 2-4: Wave functions
for panel_idx, state_idx in enumerate(idx_near_zero[:3]):
    ax = axes[(panel_idx + 1) // 2, (panel_idx + 1) % 2]
    psi = eigenvectors[:, state_idx]
    psi = psi / np.max(np.abs(psi))
    ax.plot(a, psi, color='#f472b6', linewidth=1.5, label=f'Psi (E={eigenvalues[state_idx]:.4f})')
    ax.fill_between(a, psi, alpha=0.2, color='#ec4899')
    ax.axhline(y=0, color='#6b7280', linewidth=0.5, linestyle='--')
    ax.set_xlabel('Scale factor a', color='#d1d5db')
    ax.set_ylabel('Psi(a)', color='#d1d5db')
    ax.set_title(f'WDW State {panel_idx} (E={eigenvalues[state_idx]:.4f})', color='#f9a8d4', fontsize=12, fontweight='bold')
    ax.set_facecolor('#0f172a')
    ax.tick_params(colors='#9ca3af')
    ax.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#d1d5db', fontsize=9)
    for spine in ax.spines.values():
        spine.set_color('#334155')

# Panel 4: Probability density
ax4 = axes[1, 1]
for i, idx in enumerate(idx_near_zero[:3]):
    psi = eigenvectors[:, idx]
    prob = np.abs(psi)**2
    prob = prob / np.max(prob)
    ax4.plot(a, prob, linewidth=1.5, label=f'|Psi_{i}|^2')
ax4.set_xlabel('Scale factor a', color='#d1d5db')
ax4.set_ylabel('|Psi(a)|^2', color='#d1d5db')
ax4.set_title('Probability Densities', color='#f9a8d4', fontsize=12, fontweight='bold')
ax4.set_facecolor('#0f172a')
ax4.tick_params(colors='#9ca3af')
ax4.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#d1d5db', fontsize=9)
for spine in ax4.spines.values():
    spine.set_color('#334155')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()

print()
print("The WDW wave function peaks near the classical turning points")
print("of the potential, analogous to a quantum particle in a well.")
print("The tunneling tail beyond V(a)=0 corresponds to quantum")
print("creation of the universe from 'nothing'.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

CH 02

Loop Quantum Gravity

Loop Quantum Gravity (LQG) reformulates general relativity using connection variables instead of the metric. The key insight, due to Abhay Ashtekar, is to use a complexified connection as the fundamental configuration variable, analogous to the vector potential in Yang-Mills theory. This enables the application of non-perturbative quantization techniques borrowed from gauge theory.

The Ashtekar-Barbero connection combines the spin connection Γ with the extrinsic curvature K, weighted by the Barbero-Immirzi parameter γ:

$$A_a^i = \Gamma_a^i + \gamma K_a^i$$

The Barbero-Immirzi parameter γ is a free parameter of the theory that does not affect classical physics but has profound quantum consequences — it sets the scale of the discrete area and volume spectra. Its value is typically fixed by requiring agreement with the Bekenstein-Hawking black hole entropy formula, giving γ ≈ 0.2375.

Holonomies and Wilson Loops

Rather than quantizing the connection directly, LQG works with holonomies — parallel transport operators along edges of a graph embedded in the spatial manifold. The holonomy of the connection along an edge e is:

$$h_e[A] = \mathcal{P}\exp\left(\int_e A\right)$$

The “P” denotes path ordering, necessary because the connection is a matrix-valued one-form. Holonomies are the basic building blocks of the quantum theory: states in the kinematical Hilbert space are functionals of holonomies, known as cylindrical functions. The space of all such states, after taking a projective limit over all graphs, forms the kinematical Hilbert space of LQG.

Area Quantization

The most celebrated prediction of LQG is the discreteness of geometric operators. The area operator, when acting on a spin network state pierced by edges carrying spin labels j_p, has a discrete spectrum:

$$\hat{A} = 8\pi\gamma\ell_P^2 \sum_p \sqrt{j_p(j_p+1)}$$

This means that area comes in discrete quanta, with a minimum non-zero area eigenvalue (the area gap) occurring at j = 1/2. This discreteness has profound implications: it suggests that spacetime has an atomic structure at the Planck scale, resolving the classical singularities of general relativity.

Lab: LQG Area Spectrum

LQG: Discrete Area Spectrum for j = 1/2 to j = 5

Python

script.py113 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

print("=" * 60)
print("  Loop Quantum Gravity: Discrete Area Spectrum")
print("=" * 60)
print()

# Barbero-Immirzi parameter (fixed by BH entropy)
gamma = 0.2375
# Planck length squared (set to 1 in Planck units)
lP2 = 1.0

# Compute area eigenvalues for j = 1/2, 1, 3/2, ..., 5
j_values = np.arange(0.5, 5.5, 0.5)
area_eigenvalues = 8 * np.pi * gamma * lP2 * np.sqrt(j_values * (j_values + 1))

print(f"Barbero-Immirzi parameter gamma = {gamma}")
print(f"Area eigenvalues (in Planck units):")
print(f"{'j':>6s}  {'sqrt(j(j+1))':>14s}  {'A / (8*pi*gamma*lP^2)':>22s}  {'A (Planck)':>12s}")
print("-" * 60)
for j, A in zip(j_values, area_eigenvalues):
    sq = np.sqrt(j * (j + 1))
    print(f"{j:6.1f}  {sq:14.6f}  {sq:22.6f}  {A:12.6f}")

# Multi-puncture area spectra: surfaces pierced by N punctures
print()
print("Area spectrum for surfaces with 1-3 punctures (j=1/2 each):")
for n_punctures in range(1, 4):
    A_min = n_punctures * area_eigenvalues[0]
    print(f"  {n_punctures} puncture(s): A_min = {A_min:.6f} lP^2")

# Plot
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
fig.patch.set_facecolor('#0f172a')

# Panel 1: Discrete area eigenvalues (stem plot)
ax1 = axes[0]
markerline, stemlines, baseline = ax1.stem(j_values, area_eigenvalues, linefmt='-', markerfmt='o', basefmt='')
plt.setp(markerline, color='#f472b6', markersize=8, zorder=5)
plt.setp(stemlines, color='#ec4899', linewidth=1.5, alpha=0.7)
plt.setp(baseline, visible=False)
ax1.set_xlabel('Spin quantum number j', color='#d1d5db', fontsize=12)
ax1.set_ylabel('Area eigenvalue (Planck units)', color='#d1d5db', fontsize=12)
ax1.set_title('LQG Area Spectrum', color='#f9a8d4', fontsize=14, fontweight='bold')
ax1.set_facecolor('#0f172a')
ax1.tick_params(colors='#9ca3af')
ax1.grid(True, alpha=0.15, color='#6b7280')
for spine in ax1.spines.values():
    spine.set_color('#334155')

# Annotate the area gap
ax1.annotate(f'Area gap\nA_min = {area_eigenvalues[0]:.4f}',
             xy=(0.5, area_eigenvalues[0]),
             xytext=(2.0, area_eigenvalues[0] + 2),
             arrowprops=dict(arrowstyle='->', color='#f9a8d4', lw=1.5),
             color='#f9a8d4', fontsize=10,
             bbox=dict(boxstyle='round,pad=0.3', facecolor='#1e1b4b', edgecolor='#ec4899', alpha=0.8))

# Panel 2: Area vs j with classical comparison
ax2 = axes[1]
j_cont = np.linspace(0.5, 5.0, 200)
A_cont = 8 * np.pi * gamma * lP2 * np.sqrt(j_cont * (j_cont + 1))
A_classical = 8 * np.pi * gamma * lP2 * j_cont  # classical: A ~ j for large j
ax2.plot(j_cont, A_cont, color='#f472b6', linewidth=2, label='Quantum: sqrt(j(j+1))')
ax2.plot(j_cont, A_classical, color='#818cf8', linewidth=1.5, linestyle='--', label='Classical: j (large j limit)')
ax2.scatter(j_values, area_eigenvalues, color='#ec4899', s=60, zorder=5, edgecolors='white', linewidth=0.5)
ax2.set_xlabel('j', color='#d1d5db', fontsize=12)
ax2.set_ylabel('Area', color='#d1d5db', fontsize=12)
ax2.set_title('Quantum vs Classical Scaling', color='#f9a8d4', fontsize=14, fontweight='bold')
ax2.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#d1d5db', fontsize=10)
ax2.set_facecolor('#0f172a')
ax2.tick_params(colors='#9ca3af')
ax2.grid(True, alpha=0.15, color='#6b7280')
for spine in ax2.spines.values():
    spine.set_color('#334155')

# Panel 3: Degeneracy (number of ways to reach a given total area)
ax3 = axes[2]
# For a surface with up to 4 punctures, each carrying j=1/2 or j=1
j_options = [0.5, 1.0, 1.5, 2.0]
single_areas = [8 * np.pi * gamma * lP2 * np.sqrt(j * (j + 1)) for j in j_options]
# Simple combinatorial: count multi-puncture states up to 4 punctures
from itertools import combinations_with_replacement
total_areas = []
for n in range(1, 5):
    for combo in combinations_with_replacement(j_options, n):
        A_total = sum(8 * np.pi * gamma * lP2 * np.sqrt(j * (j + 1)) for j in combo)
        total_areas.append(A_total)
total_areas = np.array(total_areas)
bins = np.linspace(0, max(total_areas) + 1, 50)
ax3.hist(total_areas, bins=bins, color='#ec4899', alpha=0.7, edgecolor='#be185d')
ax3.set_xlabel('Total area (Planck units)', color='#d1d5db', fontsize=12)
ax3.set_ylabel('Number of states', color='#d1d5db', fontsize=12)
ax3.set_title('Area Spectrum Degeneracy', color='#f9a8d4', fontsize=14, fontweight='bold')
ax3.set_facecolor('#0f172a')
ax3.tick_params(colors='#9ca3af')
ax3.grid(True, alpha=0.15, color='#6b7280')
for spine in ax3.spines.values():
    spine.set_color('#334155')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()

print()
print("The area spectrum is discrete with a minimum 'area gap'.")
print("At large j, the quantum spectrum approaches the classical linear scaling.")
print("The degeneracy of states at a given total area grows exponentially,")
print("which is the origin of black hole entropy in LQG.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

CH 03

Spin Foam Models

While loop quantum gravity gives a kinematical picture of quantum geometry through spin networks, spin foam models provide the dynamics. A spin foam is a two-complex (a combinatorial structure of vertices, edges, and faces) whose faces are labeled by spins and whose edges carry intertwiners. The spin foam can be thought of as the “history” of a spin network — a path integral over quantum geometries.

The transition amplitude between an initial and final spin network state is given by summing over all spin foams interpolating between them. The EPRL (Engle-Pereira-Rovelli-Livine) model is the most developed spin foam model, giving amplitudes of the form:

$$Z = \sum_{\text{colorings}} \prod_f A_f \prod_e A_e \prod_v A_v$$

Here the sum is over all possible spin labelings (colorings) of the two-complex. The face amplitudes A_f are simple dimension factors (2j+1), the edge amplitudes A_e enforce the simplicity constraints (which ensure we are quantizing gravity rather than a topological theory), and the vertex amplitudes A_v encode the dynamics.

Vertex Amplitudes and 15j-Symbols

The vertex amplitude in the EPRL model involves 15j-symbols — a higher-order generalization of the familiar Clebsch-Gordan coefficients. Each vertex corresponds to a 4-simplex (the four-dimensional analogue of a tetrahedron), and the 15j-symbol encodes how the ten triangular faces of this 4-simplex (each carrying a spin label) are combined consistently. In the large-spin (semiclassical) limit, the vertex amplitude is dominated by the Regge action, providing evidence that spin foams reproduce classical general relativity.

6j-Symbols: Building Blocks

The 6j-symbol (or Wigner 6j-symbol) is the fundamental recoupling coefficient for three angular momenta. It appears as a building block in all spin foam vertex amplitudes. For the simpler Ponzano-Regge model (3D quantum gravity), the partition function is a product of 6j-symbols, one per tetrahedron.

Lab: 6j Symbols and Spin Foam Graph

Spin Foam Models: 6j Symbols & Foam Graph Visualization

Python

script.py165 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from math import factorial, sqrt

print("=" * 60)
print("  Spin Foam Models: 6j Symbols & Foam Visualization")
print("=" * 60)
print()

# ── Wigner 6j symbol computation ──
# {j1 j2 j3}
# {j4 j5 j6}
# Using Racah formula

def triangle_coeff(a, b, c):
    """Triangle coefficient Delta(a,b,c)"""
    a, b, c = int(2*a), int(2*b), int(2*c)
    num = factorial((a+b-c)//2) * factorial((a-b+c)//2) * factorial((-a+b+c)//2)
    den = factorial((a+b+c)//2 + 1)
    return num / den

def wigner_6j(j1, j2, j3, j4, j5, j6):
    """Compute Wigner 6j symbol using Racah formula"""
    # Check triangle conditions
    triads = [(j1,j2,j3), (j1,j5,j6), (j4,j2,j6), (j4,j5,j3)]
    for a, b, c in triads:
        if abs(a-b) > c or a+b < c:
            return 0.0
        if (2*(a+b+c)) % 2 != 0:
            return 0.0

prefactor = sqrt(triangle_coeff(j1,j2,j3) * triangle_coeff(j1,j5,j6) *
                     triangle_coeff(j4,j2,j6) * triangle_coeff(j4,j5,j3))

# Sum over t
    t_min = max(int(j1+j2+j3), int(j1+j5+j6), int(j4+j2+j6), int(j4+j5+j3))
    t_max = min(int(j1+j2+j4+j5), int(j2+j3+j5+j6), int(j3+j1+j6+j4))

result = 0.0
    for t in range(t_min, t_max + 1):
        num = factorial(t + 1)
        d1 = factorial(t - int(j1+j2+j3))
        d2 = factorial(t - int(j1+j5+j6))
        d3 = factorial(t - int(j4+j2+j6))
        d4 = factorial(t - int(j4+j5+j3))
        d5 = factorial(int(j1+j2+j4+j5) - t)
        d6 = factorial(int(j2+j3+j5+j6) - t)
        d7 = factorial(int(j3+j1+j6+j4) - t)
        result += ((-1)**t * num) / (d1*d2*d3*d4*d5*d6*d7)

return prefactor * result

# Compute 6j symbols for various spin configurations
print("Wigner 6j Symbols for selected configurations:")
print(f"{'j1':>4s} {'j2':>4s} {'j3':>4s} {'j4':>4s} {'j5':>4s} {'j6':>4s}  {'6j-symbol':>12s}")
print("-" * 52)

configs = [
    (0.5, 0.5, 1, 0.5, 0.5, 1),
    (1, 1, 1, 1, 1, 1),
    (0.5, 0.5, 0, 0.5, 0.5, 0),
    (1, 1, 2, 1, 1, 2),
    (1.5, 1, 0.5, 1, 1.5, 0.5),
    (1, 1, 0, 1, 1, 0),
    (2, 1, 1, 1, 2, 1),
    (1.5, 1.5, 1, 1.5, 1.5, 1),
]

results_6j = []
for c in configs:
    val = wigner_6j(*c)
    results_6j.append(val)
    print(f"{c[0]:4.1f} {c[1]:4.1f} {c[2]:4.1f} {c[3]:4.1f} {c[4]:4.1f} {c[5]:4.1f}  {val:12.6f}")

# ── Spin foam graph visualization ──
print()
print("Visualizing a simple spin foam (triangulation of a single 4-simplex)...")

fig, axes = plt.subplots(1, 2, figsize=(16, 7))
fig.patch.set_facecolor('#0f172a')

# Panel 1: 6j symbol values as bar chart
ax1 = axes[0]
labels = [f'({c[0]},{c[1]},{c[2]})\n({c[3]},{c[4]},{c[5]})' for c in configs]
colors = ['#f472b6' if v > 0 else '#818cf8' for v in results_6j]
bars = ax1.bar(range(len(results_6j)), results_6j, color=colors, edgecolor='#be185d', width=0.6)
ax1.set_xticks(range(len(results_6j)))
ax1.set_xticklabels(labels, fontsize=7, color='#d1d5db')
ax1.set_ylabel('6j-symbol value', color='#d1d5db', fontsize=12)
ax1.set_title('Wigner 6j-Symbols', color='#f9a8d4', fontsize=14, fontweight='bold')
ax1.axhline(y=0, color='#6b7280', linewidth=0.5)
ax1.set_facecolor('#0f172a')
ax1.tick_params(colors='#9ca3af')
ax1.grid(True, alpha=0.15, color='#6b7280', axis='y')
for spine in ax1.spines.values():
    spine.set_color('#334155')

# Panel 2: Spin foam diagram (2-complex)
ax2 = axes[1]
ax2.set_xlim(-1.5, 1.5)
ax2.set_ylim(-1.5, 1.5)
ax2.set_aspect('equal')
ax2.set_title('Spin Foam: 4-Simplex Dual', color='#f9a8d4', fontsize=14, fontweight='bold')
ax2.set_facecolor('#0f172a')
ax2.tick_params(colors='#9ca3af')
for spine in ax2.spines.values():
    spine.set_color('#334155')

# Vertices of the spin foam (5 vertices of a 4-simplex projected to 2D)
angles = np.linspace(0, 2*np.pi, 6)[:-1]
vx = np.cos(angles)
vy = np.sin(angles)
vertex_labels = ['v0', 'v1', 'v2', 'v3', 'v4']

# Draw faces (triangles) as colored patches
from matplotlib.patches import Polygon
from itertools import combinations
face_spins = {(0,1,2): 1.0, (0,1,3): 0.5, (0,2,3): 1.5,
              (1,2,3): 1.0, (0,1,4): 0.5, (1,3,4): 1.0,
              (0,3,4): 0.5, (2,3,4): 1.5, (0,2,4): 1.0, (1,2,4): 0.5}

cmap = plt.cm.magma
for (i, j, k), spin in face_spins.items():
    triangle = Polygon([(vx[i], vy[i]), (vx[j], vy[j]), (vx[k], vy[k])],
                       alpha=0.08 + 0.05*spin, facecolor=cmap(spin/2.0),
                       edgecolor='none')
    ax2.add_patch(triangle)

# Draw edges
edge_spins = {}
for (i, j, k), spin in face_spins.items():
    for a, b in [(i,j), (j,k), (i,k)]:
        key = (min(a,b), max(a,b))
        if key not in edge_spins:
            edge_spins[key] = spin

for (i, j), spin in edge_spins.items():
    lw = 0.5 + spin
    ax2.plot([vx[i], vx[j]], [vy[i], vy[j]], color='#f472b6',
             linewidth=lw, alpha=0.7, zorder=2)
    # Label edges with spin
    mx, my = (vx[i]+vx[j])/2, (vy[i]+vy[j])/2
    ax2.text(mx*1.05, my*1.05, f'j={spin}', fontsize=6, color='#fda4af',
             ha='center', va='center',
             bbox=dict(boxstyle='round,pad=0.15', facecolor='#1e1b4b', edgecolor='none', alpha=0.8))

# Draw vertices
for i in range(5):
    ax2.scatter(vx[i], vy[i], s=120, color='#ec4899', edgecolors='white',
                linewidth=1.5, zorder=5)
    offset = 0.18
    ax2.text(vx[i]*(1+offset), vy[i]*(1+offset), vertex_labels[i],
             fontsize=10, color='#fda4af', ha='center', va='center', fontweight='bold')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()

print()
print("The spin foam graph shows the dual of a 4-simplex triangulation.")
print("Each face (triangle) carries a spin label j, and the vertex amplitude")
print("is constructed from 6j-symbols evaluated at these spins.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

CH 04

AdS/CFT Correspondence

The AdS/CFT correspondence, proposed by Juan Maldacena in 1997, is arguably the most important theoretical development in quantum gravity since the discovery of Hawking radiation. It states that a theory of quantum gravity in (d+1)-dimensional anti-de Sitter space is exactly dual to a conformal field theory living on the d-dimensional boundary of that space. This is a concrete realization of the holographic principle: all the information about the bulk gravitational physics is encoded in a lower-dimensional non-gravitational theory.

The precise statement of the duality, often called the GKPW dictionary (Gubser-Klebanov-Polyakov-Witten), equates the gravitational partition function with specified boundary conditions to the generating functional of the boundary CFT:

$$Z_{\text{gravity}}[\phi_0] = \langle e^{\int \phi_0 \mathcal{O}} \rangle_{\text{CFT}}$$

Here φ₀ is the boundary value of a bulk field, and O is the dual CFT operator. The mass of the bulk field determines the conformal dimension of the boundary operator. This dictionary allows us to translate difficult quantum gravity questions into (sometimes tractable) field theory calculations and vice versa.

Central Charge and Brown-Henneaux

In three-dimensional gravity (AdS₃), the asymptotic symmetry algebra is two copies of the Virasoro algebra with central charge determined by the AdS radius and Newton’s constant:

$$c = \frac{3L}{2G_N}$$

This Brown-Henneaux central charge provides the bridge between the bulk gravitational theory and the boundary CFT. It determines the density of states at high energies via the Cardy formula, which in turn reproduces the Bekenstein-Hawking entropy of BTZ black holes — a remarkable check of the correspondence.

Lab: BTZ Black Hole Entropy vs Cardy Formula

BTZ Black Hole Entropy vs Cardy Formula

Python

script.py107 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

print("=" * 60)
print("  AdS/CFT: BTZ Black Hole Entropy vs Cardy Formula")
print("=" * 60)
print()

# The BTZ black hole in AdS3
# Metric: ds^2 = -(r^2/L^2 - M) dt^2 + (r^2/L^2 - M)^{-1} dr^2 + r^2 dphi^2
# Horizon radius: r_+ = L * sqrt(M)
# Bekenstein-Hawking entropy: S_BH = 2*pi*r_+ / (4*G_N) = pi*L*sqrt(M) / (2*G_N)

# Parameters
L = 1.0          # AdS radius
G_N = 0.01       # Newton's constant (small for semiclassical regime)

# Central charge
c = 3 * L / (2 * G_N)
print(f"AdS radius L = {L}")
print(f"Newton's constant G_N = {G_N}")
print(f"Brown-Henneaux central charge c = 3L/(2G_N) = {c:.1f}")
print()

# BTZ mass range
M_values = np.linspace(0.1, 20.0, 200)

# Bekenstein-Hawking entropy
r_plus = L * np.sqrt(M_values)
S_BH = 2 * np.pi * r_plus / (4 * G_N)

# Cardy formula: S = 2*pi*sqrt(c*Delta/6) where Delta = M*L/2 for non-rotating BTZ
# For the BTZ: Delta_L = Delta_R = M*L/2 (non-rotating)
# Total: S_Cardy = 2*pi*sqrt(c*M*L/12) + 2*pi*sqrt(c*M*L/12) = 4*pi*sqrt(c*M*L/12)
# Simplify: with c = 3L/(2G_N):
#   S_Cardy = 4*pi*sqrt(3L/(2G_N) * M*L / 12) = 4*pi*sqrt(M*L^2/(8*G_N))
#           = 4*pi*L*sqrt(M)/(2*sqrt(2*G_N)) ... let's just compute directly

Delta = M_values * L / 2  # conformal dimension
S_Cardy = 2 * np.pi * np.sqrt(c * Delta / 6) * 2  # factor of 2 for L+R movers

print("Comparison at selected masses:")
print(f"{'M':>8s}  {'r+':>8s}  {'S_BH':>10s}  {'S_Cardy':>10s}  {'Ratio':>8s}")
print("-" * 52)
for M in [0.5, 1.0, 2.0, 5.0, 10.0, 20.0]:
    rp = L * np.sqrt(M)
    s_bh = 2 * np.pi * rp / (4 * G_N)
    delta = M * L / 2
    s_c = 2 * 2 * np.pi * np.sqrt(c * delta / 6)
    print(f"{M:8.1f}  {rp:8.4f}  {s_bh:10.4f}  {s_c:10.4f}  {s_c/s_bh:8.6f}")

# Plot
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
fig.patch.set_facecolor('#0f172a')

# Panel 1: S_BH vs S_Cardy
ax1 = axes[0]
ax1.plot(M_values, S_BH, color='#f472b6', linewidth=2.5, label='Bekenstein-Hawking')
ax1.plot(M_values, S_Cardy, color='#818cf8', linewidth=2, linestyle='--', label='Cardy formula')
ax1.set_xlabel('BTZ mass M', color='#d1d5db', fontsize=12)
ax1.set_ylabel('Entropy S', color='#d1d5db', fontsize=12)
ax1.set_title('BTZ Entropy: BH vs Cardy', color='#f9a8d4', fontsize=14, fontweight='bold')
ax1.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#d1d5db', fontsize=10)
ax1.set_facecolor('#0f172a')
ax1.tick_params(colors='#9ca3af')
ax1.grid(True, alpha=0.15, color='#6b7280')
for spine in ax1.spines.values():
    spine.set_color('#334155')

# Panel 2: Ratio S_Cardy / S_BH
ax2 = axes[1]
ratio = S_Cardy / S_BH
ax2.plot(M_values, ratio, color='#f472b6', linewidth=2)
ax2.axhline(y=1.0, color='#6b7280', linewidth=1, linestyle='--')
ax2.set_xlabel('BTZ mass M', color='#d1d5db', fontsize=12)
ax2.set_ylabel('S_Cardy / S_BH', color='#d1d5db', fontsize=12)
ax2.set_title('Agreement Ratio', color='#f9a8d4', fontsize=14, fontweight='bold')
ax2.set_ylim(0.9, 1.1)
ax2.set_facecolor('#0f172a')
ax2.tick_params(colors='#9ca3af')
ax2.grid(True, alpha=0.15, color='#6b7280')
for spine in ax2.spines.values():
    spine.set_color('#334155')

# Panel 3: Hawking temperature
T_H = np.sqrt(M_values) / (2 * np.pi * L)
ax3 = axes[2]
ax3.plot(M_values, T_H, color='#fb923c', linewidth=2)
ax3.set_xlabel('BTZ mass M', color='#d1d5db', fontsize=12)
ax3.set_ylabel('Hawking Temperature T_H', color='#d1d5db', fontsize=12)
ax3.set_title('BTZ Hawking Temperature', color='#f9a8d4', fontsize=14, fontweight='bold')
ax3.set_facecolor('#0f172a')
ax3.tick_params(colors='#9ca3af')
ax3.grid(True, alpha=0.15, color='#6b7280')
for spine in ax3.spines.values():
    spine.set_color('#334155')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()

print()
print("The Cardy formula exactly reproduces the Bekenstein-Hawking entropy!")
print("This is a precision test of the AdS3/CFT2 correspondence.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

CH 05

Holographic Entanglement Entropy

One of the most profound consequences of AdS/CFT is the geometric realization of entanglement entropy. In 2006, Ryu and Takayanagi (RT) proposed that the entanglement entropy of a spatial region A in the boundary CFT is computed by the area of a certain minimal surface in the bulk:

$$S_A = \frac{\text{Area}(\gamma_A)}{4G_N}$$

Here γ_A is the minimal area surface in the bulk that is homologous to the boundary region A (i.e., shares the same boundary ∂A). This formula has the same structure as the Bekenstein-Hawking entropy formula, but now applied to surfaces that are not event horizons. The RT formula was later proved by Lewkowycz and Maldacena using the replica trick in gravity.

RT Surface in AdS₃

In AdS₃, the minimal surface anchored to an interval of length lon the boundary is a geodesic in the hyperbolic plane. The Poincare patch metric is ds² = L²(dz² + dx²)/z², and the geodesic connecting the endpoints of the interval at (x₁, 0) and (x₂, 0) is a semicircle. Its regulated length gives the famous logarithmic entanglement entropy of a 1+1 dimensional CFT: S = (c/3) log(l/ε).

Lab: RT Surface in AdS₃

RT Surface in AdS3: Geodesic Computation & Entanglement Entropy

Python

script.py142 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

print("=" * 60)
print("  Holographic Entanglement Entropy: RT Surface in AdS3")
print("=" * 60)
print()

# In AdS3 Poincare patch: ds^2 = L^2 (dz^2 + dx^2) / z^2
# The RT surface for an interval [x1, x2] on the boundary (z=0)
# is a semicircle: (x - x_center)^2 + z^2 = R^2
# where x_center = (x1+x2)/2 and R = (x2-x1)/2

L_AdS = 1.0   # AdS radius
G_N = 0.01    # Newton's constant
c = 3 * L_AdS / (2 * G_N)  # central charge
epsilon = 0.01  # UV cutoff

print(f"AdS radius L = {L_AdS}")
print(f"G_N = {G_N}, central charge c = {c:.1f}")
print(f"UV cutoff epsilon = {epsilon}")
print()

# Compute RT surface for various interval lengths
intervals = [(0, 1), (-1, 1), (-2, 2), (0, 3)]

print("Entanglement entropy for different intervals:")
print(f"{'Interval':>12s}  {'Length l':>8s}  {'S_RT':>10s}  {'S_CFT':>10s}")
print("-" * 48)

for x1, x2 in intervals:
    l = x2 - x1
    # RT formula: Area of geodesic (regularized) = 2*L*log(l/epsilon)
    # S = Area / (4*G_N) = L * log(l/epsilon) / (2*G_N) = (c/3) * log(l/epsilon)
    S_RT = 2 * L_AdS * np.log(l / epsilon) / (4 * G_N)
    S_CFT = (c / 3) * np.log(l / epsilon)
    print(f"  [{x1:+.0f}, {x2:+.0f}]  {l:8.1f}  {S_RT:10.4f}  {S_CFT:10.4f}")

# Plot
fig, axes = plt.subplots(1, 3, figsize=(18, 7))
fig.patch.set_facecolor('#0f172a')

# Panel 1: RT surfaces (geodesics) in the Poincare half-plane
ax1 = axes[0]
theta = np.linspace(0, np.pi, 300)

colors_rt = ['#f472b6', '#c084fc', '#818cf8', '#fb923c']
for idx, (x1, x2) in enumerate(intervals):
    R = (x2 - x1) / 2
    xc = (x1 + x2) / 2
    x_geo = xc + R * np.cos(theta)
    z_geo = R * np.sin(theta)
    ax1.plot(x_geo, z_geo, color=colors_rt[idx], linewidth=2,
             label=f'[{x1},{x2}], l={x2-x1}')
    # Mark the boundary endpoints
    ax1.scatter([x1, x2], [0, 0], color=colors_rt[idx], s=40, zorder=5)

# Draw the boundary
ax1.axhline(y=0, color='#6b7280', linewidth=1.5, linestyle='-')
ax1.text(3.5, 0.05, 'Boundary (z=0)', color='#9ca3af', fontsize=9)

# UV cutoff line
ax1.axhline(y=epsilon, color='#ef4444', linewidth=0.8, linestyle=':', alpha=0.5)
ax1.text(3.5, epsilon + 0.05, f'cutoff z={epsilon}', color='#ef4444', fontsize=8, alpha=0.7)

ax1.set_xlabel('x (boundary coordinate)', color='#d1d5db', fontsize=12)
ax1.set_ylabel('z (radial / bulk)', color='#d1d5db', fontsize=12)
ax1.set_title('RT Geodesics in AdS3 Poincare Patch', color='#f9a8d4', fontsize=13, fontweight='bold')
ax1.set_ylim(-0.1, 3.5)
ax1.set_xlim(-4, 5)
ax1.invert_yaxis()
ax1.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#d1d5db', fontsize=9, loc='lower right')
ax1.set_facecolor('#0f172a')
ax1.tick_params(colors='#9ca3af')
ax1.grid(True, alpha=0.1, color='#6b7280')
for spine in ax1.spines.values():
    spine.set_color('#334155')

# Panel 2: Entanglement entropy vs interval length
ax2 = axes[1]
l_range = np.linspace(0.1, 10.0, 300)
S_vs_l = (c / 3) * np.log(l_range / epsilon)
ax2.plot(l_range, S_vs_l, color='#f472b6', linewidth=2.5)
# Mark the computed intervals
for idx, (x1, x2) in enumerate(intervals):
    l = x2 - x1
    S = (c / 3) * np.log(l / epsilon)
    ax2.scatter([l], [S], color=colors_rt[idx], s=80, zorder=5, edgecolors='white', linewidth=1)

ax2.set_xlabel('Interval length l', color='#d1d5db', fontsize=12)
ax2.set_ylabel('Entanglement entropy S', color='#d1d5db', fontsize=12)
ax2.set_title('S = (c/3) log(l/epsilon)', color='#f9a8d4', fontsize=13, fontweight='bold')
ax2.set_facecolor('#0f172a')
ax2.tick_params(colors='#9ca3af')
ax2.grid(True, alpha=0.15, color='#6b7280')
for spine in ax2.spines.values():
    spine.set_color('#334155')

# Panel 3: Mutual information (two disjoint intervals)
ax3 = axes[2]
# I(A:B) = S(A) + S(B) - S(A union B)
# For two intervals of length l separated by distance d:
d_range = np.linspace(0.1, 5.0, 200)
l_A = 2.0
l_B = 2.0
S_A = (c / 3) * np.log(l_A / epsilon)
S_B = (c / 3) * np.log(l_B / epsilon)
# S(A union B) involves two possible RT surfaces:
# disconnected: S = S_A + S_B
# connected: S = (c/3) * [log((l_A+d)/eps) + log((l_B+d)/eps) - ...]
# Simplified: mutual information for adjacent intervals
S_AB = (c / 3) * np.log((l_A + l_B + d_range) / epsilon)
# Actually use the proper formula
S_union_disconnected = np.full_like(d_range, S_A + S_B)
S_cross = (c / 3) * (np.log((l_A + d_range) / epsilon) + np.log((l_B + d_range) / epsilon))
S_union_connected = np.minimum(S_union_disconnected, S_cross)

MI = S_A + S_B - S_union_connected
MI = np.maximum(MI, 0)  # mutual information is non-negative

ax3.plot(d_range, MI, color='#f472b6', linewidth=2.5)
ax3.fill_between(d_range, MI, alpha=0.15, color='#ec4899')
ax3.set_xlabel('Separation d', color='#d1d5db', fontsize=12)
ax3.set_ylabel('Mutual Information I(A:B)', color='#d1d5db', fontsize=12)
ax3.set_title(f'Holographic MI (l_A=l_B={l_A})', color='#f9a8d4', fontsize=13, fontweight='bold')
ax3.set_facecolor('#0f172a')
ax3.tick_params(colors='#9ca3af')
ax3.grid(True, alpha=0.15, color='#6b7280')
for spine in ax3.spines.values():
    spine.set_color('#334155')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()

print()
print("RT surfaces are semicircular geodesics in the Poincare half-plane.")
print("The entanglement entropy grows logarithmically with interval length,")
print("matching the universal CFT result S = (c/3) log(l/epsilon).")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

CH 06

Tensor Networks

Tensor networks provide a powerful language for describing many-body quantum states that lies at the intersection of quantum information theory, condensed matter physics, and quantum gravity. The key idea is to represent a quantum state not as a single exponentially large vector, but as a network of smaller tensors connected by contracted indices. The pattern of contractions encodes the entanglement structure of the state.

MERA and Holography

The Multi-scale Entanglement Renormalization Ansatz (MERA) is a tensor network that builds a quantum state through layers of disentanglers and isometries:

$$|\Psi\rangle = \prod_{\text{layers}} U_{\text{disentanglers}} \cdot U_{\text{isometries}} |0\rangle^{\otimes n}$$

Brian Swingle observed that the MERA network has the geometric structure of a discretized AdS space, with the layers corresponding to different radial (energy) scales in the bulk. This provides an explicit tensor network realization of the AdS/CFT correspondence, where the entanglement structure of the boundary state is literally “built up” layer by layer from the bulk.

Matrix Product States (MPS)

The simplest tensor network is the Matrix Product State (MPS), which represents a 1D quantum state as a chain of matrices. For a system of n qubits, each site carries a tensor A^[k]_s of dimensions χ × χ, where χ is the bond dimension controlling the entanglement capacity. The entanglement entropy across any bipartition is bounded by log(χ), making MPS efficient for gapped ground states that obey an area law.

Lab: MPS for 6 Qubits

Matrix Product State: 6-Qubit Entanglement Analysis

Python

script.py172 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

print("=" * 60)
print("  Tensor Networks: MPS for 6 Qubits")
print("=" * 60)
print()

# Build a Matrix Product State (MPS) for 6 qubits
# Physical dimension d=2 (qubit), bond dimension chi
n_sites = 6
d = 2  # physical dimension (qubit)

# We'll build MPS for various bond dimensions and compute
# entanglement entropy across the central cut

def build_random_mps(n, d, chi):
    """Build a random MPS with given bond dimension"""
    tensors = []
    for i in range(n):
        chi_l = 1 if i == 0 else min(chi, d**i, d**(n-i))
        chi_r = 1 if i == n-1 else min(chi, d**(i+1), d**(n-i-1))
        # Random tensor of shape (chi_l, d, chi_r)
        A = np.random.randn(chi_l, d, chi_r) + 1j * np.random.randn(chi_l, d, chi_r)
        tensors.append(A)
    return tensors

def contract_mps(tensors):
    """Contract MPS to get the full state vector"""
    n = len(tensors)
    d = tensors[0].shape[1]
    # Start from left
    result = tensors[0]  # shape (1, d, chi_r)
    for i in range(1, n):
        # result shape: (1, d^i, chi_current)
        # tensor shape: (chi_current, d, chi_next)
        chi_l = result.shape[2]
        result = np.tensordot(result, tensors[i], axes=([2], [0]))
        # result shape: (1, d^i, d, chi_next)
        # Reshape to merge physical indices
        shape = result.shape
        result = result.reshape(shape[0], -1, shape[-1])
    # Final shape: (1, d^n, 1) -> (d^n,)
    return result.reshape(-1)

def entanglement_entropy_mps(tensors, cut):
    """Compute entanglement entropy across a cut between site 'cut-1' and 'cut'"""
    n = len(tensors)
    # Contract left block
    left = tensors[0]
    for i in range(1, cut):
        left = np.tensordot(left, tensors[i], axes=([2], [0]))
        shape = left.shape
        left = left.reshape(shape[0], -1, shape[-1])
    # Contract right block
    right = tensors[cut]
    for i in range(cut + 1, n):
        right = np.tensordot(right, tensors[i], axes=([2], [0]))
        shape = right.shape
        right = right.reshape(shape[0], -1, shape[-1])
    # Form the state matrix: left has shape (1, d^cut, chi), right has (chi, d^(n-cut), 1)
    left_mat = left.reshape(-1, left.shape[-1])    # (d^cut, chi)
    right_mat = right.reshape(right.shape[0], -1)  # (chi, d^(n-cut))
    # Full state matrix
    psi_mat = left_mat @ right_mat  # (d^cut, d^(n-cut))
    # SVD
    U, s, Vh = np.linalg.svd(psi_mat, full_matrices=False)
    s = s[s > 1e-12]
    s = s / np.linalg.norm(s)
    probs = s**2
    entropy = -np.sum(probs * np.log2(probs + 1e-30))
    return entropy

# Compute entanglement profiles for different bond dimensions
chi_values = [1, 2, 4, 8, 16]
np.random.seed(42)

print(f"System: {n_sites} qubits, d = {d}")
print()
print("Entanglement entropy across central cut (site 2|3):")
print(f"{'chi':>6s}  {'S (bits)':>10s}  {'S_max = log2(chi)':>18s}")
print("-" * 40)

all_profiles = {}
for chi in chi_values:
    tensors = build_random_mps(n_sites, d, chi)
    psi = contract_mps(tensors)
    psi = psi / np.linalg.norm(psi)

# Compute entanglement profile (entropy at each cut)
    profile = []
    for cut in range(1, n_sites):
        S = entanglement_entropy_mps(tensors, cut)
        profile.append(S)
    all_profiles[chi] = profile

S_central = profile[n_sites // 2 - 1]
    S_max = np.log2(chi) if chi > 1 else 0
    print(f"{chi:6d}  {S_central:10.4f}  {S_max:18.4f}")

# Plot
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
fig.patch.set_facecolor('#0f172a')

# Panel 1: Entanglement profiles
ax1 = axes[0]
colors_mps = ['#6b7280', '#818cf8', '#c084fc', '#f472b6', '#fb923c']
cuts = list(range(1, n_sites))
for idx, chi in enumerate(chi_values):
    ax1.plot(cuts, all_profiles[chi], 'o-', color=colors_mps[idx],
             linewidth=2, markersize=6, label=f'chi={chi}')
ax1.set_xlabel('Cut position', color='#d1d5db', fontsize=12)
ax1.set_ylabel('Entanglement entropy (bits)', color='#d1d5db', fontsize=12)
ax1.set_title('MPS Entanglement Profiles', color='#f9a8d4', fontsize=14, fontweight='bold')
ax1.set_xticks(cuts)
ax1.set_xticklabels([f'{i}|{i+1}' for i in cuts])
ax1.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#d1d5db', fontsize=10)
ax1.set_facecolor('#0f172a')
ax1.tick_params(colors='#9ca3af')
ax1.grid(True, alpha=0.15, color='#6b7280')
for spine in ax1.spines.values():
    spine.set_color('#334155')

# Panel 2: Central entropy vs bond dimension
ax2 = axes[1]
central_entropies = [all_profiles[chi][n_sites // 2 - 1] for chi in chi_values]
max_entropies = [np.log2(chi) if chi > 1 else 0 for chi in chi_values]
ax2.plot(chi_values, central_entropies, 'o-', color='#f472b6', linewidth=2,
         markersize=8, label='Measured S')
ax2.plot(chi_values, max_entropies, 's--', color='#818cf8', linewidth=1.5,
         markersize=6, label='S_max = log2(chi)')
ax2.set_xlabel('Bond dimension chi', color='#d1d5db', fontsize=12)
ax2.set_ylabel('Central cut entropy (bits)', color='#d1d5db', fontsize=12)
ax2.set_title('Entropy vs Bond Dimension', color='#f9a8d4', fontsize=14, fontweight='bold')
ax2.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#d1d5db', fontsize=10)
ax2.set_facecolor('#0f172a')
ax2.tick_params(colors='#9ca3af')
ax2.grid(True, alpha=0.15, color='#6b7280')
for spine in ax2.spines.values():
    spine.set_color('#334155')

# Panel 3: Schmidt spectrum for chi=8
ax3 = axes[2]
chi_show = 8
tensors = build_random_mps(n_sites, d, chi_show)
psi = contract_mps(tensors)
psi = psi / np.linalg.norm(psi)
# Get Schmidt values at central cut
psi_mat = psi.reshape(d**(n_sites//2), d**(n_sites//2))
U, s, Vh = np.linalg.svd(psi_mat, full_matrices=False)
s = s / np.linalg.norm(s)
ax3.bar(range(len(s[:16])), s[:16]**2, color='#f472b6', edgecolor='#be185d', alpha=0.8)
ax3.set_xlabel('Schmidt index', color='#d1d5db', fontsize=12)
ax3.set_ylabel('Schmidt weight (lambda^2)', color='#d1d5db', fontsize=12)
ax3.set_title(f'Schmidt Spectrum (chi={chi_show})', color='#f9a8d4', fontsize=14, fontweight='bold')
ax3.set_facecolor('#0f172a')
ax3.tick_params(colors='#9ca3af')
ax3.grid(True, alpha=0.15, color='#6b7280')
for spine in ax3.spines.values():
    spine.set_color('#334155')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()

print()
print("The MPS entanglement entropy is bounded by log2(chi).")
print("Increasing chi allows the MPS to capture more entanglement.")
print("This area-law bound is the tensor network analogue of the RT formula.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

CH 07

Black Hole Information Paradox

The black hole information paradox, first articulated by Stephen Hawking in 1976, is arguably the sharpest puzzle in quantum gravity. Hawking showed that black holes emit thermal radiation at a temperature inversely proportional to their mass. As the black hole evaporates, the radiation appears to be in a mixed thermal state regardless of the initial pure state that formed the black hole. If the black hole evaporates completely, this represents a violation of unitarity — pure states evolve into mixed states, destroying quantum information.

The Page Curve

Don Page argued that if black hole evaporation is unitary, the entanglement entropy of the Hawking radiation must follow a specific curve. Initially, the entropy grows as more radiation is emitted. At the Page time (roughly when half the black hole has evaporated), the entropy must begin to decrease, eventually returning to zero when the black hole is gone:

$$S(R) = \min\left(S_{\text{no-island}}, S_{\text{island}}\right)$$

The Island Formula

The recent breakthrough by Penington, Almheiri, Mahajan, Maldacena, and Zhao (2019) showed how to derive the Page curve from semiclassical gravity using the “island formula” for the entropy of Hawking radiation:

$$S = \text{ext}\left[\frac{\text{Area}(\partial I)}{4G_N} + S_{\text{matter}}(R \cup I)\right]$$

Here I is the “island” — a region behind the black hole horizon that is secretly entangled with the radiation R. The formula says to extremize over all possible islands and take the minimum. Before the Page time, the island is empty and the entropy grows. After the Page time, a non-trivial island appears behind the horizon, and the entropy begins to decrease. The competition between these two saddle points generates the Page curve, resolving the information paradox within the framework of semiclassical gravity.

Lab: The Page Curve

The Page Curve: Island Phase Transition

Python

script.py141 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

print("=" * 60)
print("  Black Hole Information Paradox: The Page Curve")
print("=" * 60)
print()

# Model the Page curve for a black hole evaporation
# Units: time in units of the total evaporation time

# BH starts with entropy S_BH_0 (= initial Bekenstein-Hawking entropy)
S_BH_0 = 100.0  # in some units

# Time parameter: t from 0 to 1 (full evaporation)
t = np.linspace(0, 0.999, 1000)

# Black hole mass decreases as M(t) = M_0 * (1 - t)^(1/3) for Hawking evap
# BH entropy: S_BH(t) = S_0 * (1 - t)^(2/3)
S_BH = S_BH_0 * (1 - t)**(2/3)

# Hawking's (incorrect) prediction: S_rad grows monotonically
# Thermal radiation entropy: S_thermal ~ S_0 * (1 - (1-t)^{2/3})
S_hawking = S_BH_0 * (1 - (1 - t)**(2/3))

# No-island saddle: S grows linearly (simplified)
# S_no_island ~ k * t (matter entropy of radiation)
k_matter = 1.5 * S_BH_0
S_no_island = k_matter * t

# Island saddle: S = BH entropy (area of island boundary) + small correction
# S_island ~ S_BH(t) + small correction
S_island = S_BH + 2.0  # small correction for matter entropy on island

# Page curve: min of the two saddles
S_page = np.minimum(S_no_island, S_island)

# Find Page time
page_idx = np.argmin(np.abs(S_no_island - S_island))
t_page = t[page_idx]
S_page_val = S_page[page_idx]

print(f"Initial BH entropy: S_BH(0) = {S_BH_0:.1f}")
print(f"Page time: t_Page = {t_page:.3f}")
print(f"Entropy at Page time: S = {S_page_val:.2f}")
print()
print("Before Page time: no-island saddle dominates (S grows)")
print("After Page time:  island saddle dominates (S decreases)")
print("At t=1: BH fully evaporated, S -> 0 (unitarity preserved)")

# Plot
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
fig.patch.set_facecolor('#0f172a')

# Panel 1: The Page Curve
ax1 = axes[0, 0]
ax1.plot(t, S_hawking, color='#ef4444', linewidth=2, linestyle='--', label="Hawking (non-unitary)", alpha=0.7)
ax1.plot(t, S_no_island, color='#818cf8', linewidth=1.5, linestyle=':', label='No-island saddle', alpha=0.6)
ax1.plot(t, S_island, color='#c084fc', linewidth=1.5, linestyle=':', label='Island saddle', alpha=0.6)
ax1.plot(t, S_page, color='#f472b6', linewidth=3, label='Page curve (unitary)')
ax1.axvline(x=t_page, color='#fbbf24', linewidth=1, linestyle='--', alpha=0.5)
ax1.annotate(f'Page time\nt = {t_page:.2f}', xy=(t_page, S_page_val),
             xytext=(t_page + 0.15, S_page_val + 20),
             arrowprops=dict(arrowstyle='->', color='#fbbf24', lw=1.5),
             color='#fbbf24', fontsize=11,
             bbox=dict(boxstyle='round,pad=0.3', facecolor='#1e1b4b', edgecolor='#fbbf24', alpha=0.8))
ax1.set_xlabel('Time (units of evaporation time)', color='#d1d5db', fontsize=12)
ax1.set_ylabel('Radiation entropy S(R)', color='#d1d5db', fontsize=12)
ax1.set_title('The Page Curve', color='#f9a8d4', fontsize=14, fontweight='bold')
ax1.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#d1d5db', fontsize=10, loc='upper right')
ax1.set_facecolor('#0f172a')
ax1.tick_params(colors='#9ca3af')
ax1.grid(True, alpha=0.15, color='#6b7280')
for spine in ax1.spines.values():
    spine.set_color('#334155')

# Panel 2: BH entropy and radiation entropy
ax2 = axes[0, 1]
ax2.fill_between(t, S_BH, alpha=0.2, color='#ec4899')
ax2.plot(t, S_BH, color='#f472b6', linewidth=2, label='BH entropy S_BH')
ax2.plot(t, S_page, color='#818cf8', linewidth=2, label='Radiation entropy (Page)')
ax2.fill_between(t, S_page, alpha=0.15, color='#818cf8')
ax2.axvline(x=t_page, color='#fbbf24', linewidth=1, linestyle='--', alpha=0.5)
ax2.set_xlabel('Time', color='#d1d5db', fontsize=12)
ax2.set_ylabel('Entropy', color='#d1d5db', fontsize=12)
ax2.set_title('BH vs Radiation Entropy', color='#f9a8d4', fontsize=14, fontweight='bold')
ax2.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#d1d5db', fontsize=10)
ax2.set_facecolor('#0f172a')
ax2.tick_params(colors='#9ca3af')
ax2.grid(True, alpha=0.15, color='#6b7280')
for spine in ax2.spines.values():
    spine.set_color('#334155')

# Panel 3: Island size
ax3 = axes[1, 0]
# Before Page time: no island. After: island grows
island_area = np.zeros_like(t)
island_area[t > t_page] = S_BH_0 * (1 - t[t > t_page])**(2/3)
ax3.fill_between(t, island_area, alpha=0.3, color='#c084fc')
ax3.plot(t, island_area, color='#c084fc', linewidth=2)
ax3.axvline(x=t_page, color='#fbbf24', linewidth=1, linestyle='--', alpha=0.5)
ax3.set_xlabel('Time', color='#d1d5db', fontsize=12)
ax3.set_ylabel('Island area / 4G_N', color='#d1d5db', fontsize=12)
ax3.set_title('Island Contribution', color='#f9a8d4', fontsize=14, fontweight='bold')
ax3.set_facecolor('#0f172a')
ax3.tick_params(colors='#9ca3af')
ax3.grid(True, alpha=0.15, color='#6b7280')
ax3.text(0.2, max(island_area)*0.5, 'No island', color='#6b7280', fontsize=14, ha='center')
ax3.text(0.75, max(island_area)*0.5, 'Island present', color='#c084fc', fontsize=14, ha='center')
for spine in ax3.spines.values():
    spine.set_color('#334155')

# Panel 4: Which saddle dominates
ax4 = axes[1, 1]
dominance = np.where(S_no_island < S_island, 0, 1)  # 0 = no-island, 1 = island
ax4.fill_between(t, 0, 1, where=(dominance == 0), alpha=0.3, color='#818cf8', label='No-island phase')
ax4.fill_between(t, 0, 1, where=(dominance == 1), alpha=0.3, color='#f472b6', label='Island phase')
ax4.axvline(x=t_page, color='#fbbf24', linewidth=2, linestyle='-', label=f'Phase transition (t={t_page:.2f})')
ax4.set_xlabel('Time', color='#d1d5db', fontsize=12)
ax4.set_ylabel('Dominant saddle', color='#d1d5db', fontsize=12)
ax4.set_title('Island Phase Transition', color='#f9a8d4', fontsize=14, fontweight='bold')
ax4.set_yticks([0.25, 0.75])
ax4.set_yticklabels(['No-island', 'Island'], color='#d1d5db')
ax4.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#d1d5db', fontsize=10, loc='upper right')
ax4.set_facecolor('#0f172a')
ax4.tick_params(colors='#9ca3af')
for spine in ax4.spines.values():
    spine.set_color('#334155')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()

print()
print("The Page curve shows a phase transition between two gravitational saddles.")
print("Before the Page time, the entropy grows (Hawking's calculation is valid).")
print("After the Page time, the island saddle dominates and entropy decreases.")
print("This resolves the information paradox: unitarity is preserved!")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

CH 08

Quantum Computing for Quantum Gravity

Quantum computers offer a natural platform for simulating quantum gravitational systems. The fundamental objects of loop quantum gravity — spin networks — are essentially networks of entangled quantum states, which can be directly mapped onto qubit registers and entangling gates. Similarly, the path integrals of spin foam models involve sums over quantum amplitudes that are exponentially costly on classical computers but potentially efficient on quantum hardware.

A spin network is a graph whose edges carry spin labels (j = 0, 1/2, 1, 3/2, ...) and whose vertices carry intertwiners (coupling tensors). On a quantum computer, each edge can be represented by a set of qubits, and the intertwiner constraints can be enforced using entangling gates. The simplest non-trivial spin network is the tetrahedron graph: four nodes connected by six edges, dual to a single tetrahedron in 3D space.

Graph States and Spin Networks

A graph state is a quantum state associated with a graph G = (V, E). Starting from |+〉^⊗n, we apply a controlled-Z gate for each edge in the graph. The resulting state has entanglement structure that mirrors the graph topology. For the tetrahedron (complete graph K₄), this produces a highly entangled 4-qubit state that serves as a toy model for a spin network node.

Lab: Spin Network on 4 Qubits

Spin Network State on 4 Qubits: Tetrahedron Graph State

Python

script.py183 lines

from qiskit import QuantumCircuit
from qiskit.quantum_info import Statevector, entropy, partial_trace, DensityMatrix
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

print("=" * 60)
print("  Quantum Computing for Quantum Gravity:")
print("  Spin Network State on 4 Qubits (Tetrahedron)")
print("=" * 60)
print()

# Create a 4-qubit spin network: tetrahedron graph state
# K4 complete graph = 6 edges
qc = QuantumCircuit(4)
qc.h(range(4))  # All qubits to |+>

# Entangle along tetrahedron edges (all pairs = complete graph K4)
edges = [(0,1), (0,2), (0,3), (1,2), (1,3), (2,3)]
for i, j in edges:
    qc.cz(i, j)

print("Circuit (tetrahedron graph state on 4 qubits):")
print(qc.draw(output='text'))
print()

# Get the statevector
sv = Statevector.from_instruction(qc)
print("Statevector amplitudes (non-zero):")
for i, amp in enumerate(sv.data):
    if abs(amp) > 1e-10:
        phase = np.angle(amp) / np.pi
        print(f"  |{i:04b}> : amp = {amp.real:+.4f}{amp.imag:+.4f}j  (phase = {phase:+.3f}*pi)")
print()

# Compute entanglement entropy for all bipartitions
print("Entanglement entropy for each single-qubit subsystem:")
for qubit in range(4):
    # Trace out all other qubits
    others = [q for q in range(4) if q != qubit]
    rho_reduced = partial_trace(sv, others)
    S = entropy(rho_reduced, base=2)
    print(f"  Qubit {qubit}: S = {S:.4f} bits")

print()
print("Entanglement entropy for 2-qubit subsystems:")
bipartitions_2 = [(0,1), (0,2), (0,3), (1,2), (1,3), (2,3)]
S_2qubit = []
for q1, q2 in bipartitions_2:
    others = [q for q in range(4) if q not in (q1, q2)]
    rho = partial_trace(sv, others)
    S = entropy(rho, base=2)
    S_2qubit.append(S)
    print(f"  Qubits ({q1},{q2}): S = {S:.4f} bits")

# Compare with a product state (no entanglement)
qc_prod = QuantumCircuit(4)
qc_prod.h(range(4))
sv_prod = Statevector.from_instruction(qc_prod)

print()
print("For comparison - product state |++++>:")
for qubit in range(4):
    others = [q for q in range(4) if q != qubit]
    rho = partial_trace(sv_prod, others)
    S = entropy(rho, base=2)
    print(f"  Qubit {qubit}: S = {S:.4f} bits (no entanglement)")

# Compare with GHZ state
qc_ghz = QuantumCircuit(4)
qc_ghz.h(0)
for i in range(3):
    qc_ghz.cx(i, i+1)
sv_ghz = Statevector.from_instruction(qc_ghz)

print()
print("For comparison - GHZ state (|0000> + |1111>)/sqrt(2):")
for qubit in range(4):
    others = [q for q in range(4) if q != qubit]
    rho = partial_trace(sv_ghz, others)
    S = entropy(rho, base=2)
    print(f"  Qubit {qubit}: S = {S:.4f} bits")

# Plot
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
fig.patch.set_facecolor('#0f172a')

# Panel 1: Tetrahedron graph
ax1 = axes[0, 0]
ax1.set_xlim(-1.5, 1.5)
ax1.set_ylim(-1.5, 1.5)
ax1.set_aspect('equal')
positions = [(-0.8, -0.8), (0.8, -0.8), (0, 0.9), (0, -0.1)]
labels = ['q0', 'q1', 'q2', 'q3']
for (i, j) in edges:
    x1, y1 = positions[i]
    x2, y2 = positions[j]
    ax1.plot([x1, x2], [y1, y2], color='#f472b6', linewidth=2, alpha=0.6, zorder=1)
    mx, my = (x1+x2)/2, (y1+y2)/2
    ax1.text(mx, my, 'CZ', fontsize=7, color='#fda4af', ha='center', va='center',
             bbox=dict(boxstyle='round,pad=0.15', facecolor='#1e1b4b', edgecolor='none', alpha=0.8))
for idx, (x, y) in enumerate(positions):
    ax1.scatter(x, y, s=300, color='#ec4899', edgecolors='white', linewidth=2, zorder=5)
    ax1.text(x, y, labels[idx], fontsize=11, color='white', ha='center', va='center', fontweight='bold')
ax1.set_title('Tetrahedron Spin Network', color='#f9a8d4', fontsize=14, fontweight='bold')
ax1.set_facecolor('#0f172a')
ax1.tick_params(colors='#0f172a')
for spine in ax1.spines.values():
    spine.set_color('#334155')

# Panel 2: Single-qubit entanglement entropies comparison
ax2 = axes[0, 1]
states_names = ['Graph (K4)', 'Product |+>^4', 'GHZ']
S_graph_1q = []
S_prod_1q = []
S_ghz_1q = []
for qubit in range(4):
    others = [q for q in range(4) if q != qubit]
    S_graph_1q.append(entropy(partial_trace(sv, others), base=2))
    S_prod_1q.append(entropy(partial_trace(sv_prod, others), base=2))
    S_ghz_1q.append(entropy(partial_trace(sv_ghz, others), base=2))

x_pos = np.arange(4)
width = 0.25
ax2.bar(x_pos - width, S_graph_1q, width, color='#f472b6', label='Graph K4', edgecolor='#be185d')
ax2.bar(x_pos, S_prod_1q, width, color='#818cf8', label='Product', edgecolor='#4338ca')
ax2.bar(x_pos + width, S_ghz_1q, width, color='#fb923c', label='GHZ', edgecolor='#c2410c')
ax2.set_xlabel('Qubit', color='#d1d5db', fontsize=12)
ax2.set_ylabel('Entanglement entropy (bits)', color='#d1d5db', fontsize=12)
ax2.set_title('Single-Qubit Entanglement', color='#f9a8d4', fontsize=14, fontweight='bold')
ax2.set_xticks(x_pos)
ax2.set_xticklabels([f'q{i}' for i in range(4)])
ax2.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#d1d5db', fontsize=10)
ax2.set_facecolor('#0f172a')
ax2.tick_params(colors='#9ca3af')
ax2.grid(True, alpha=0.15, color='#6b7280', axis='y')
for spine in ax2.spines.values():
    spine.set_color('#334155')

# Panel 3: 2-qubit subsystem entropies
ax3 = axes[1, 0]
pair_labels = [f'({q1},{q2})' for q1, q2 in bipartitions_2]
ax3.bar(range(len(S_2qubit)), S_2qubit, color='#f472b6', edgecolor='#be185d', alpha=0.8)
ax3.set_xlabel('Qubit pair', color='#d1d5db', fontsize=12)
ax3.set_ylabel('Entanglement entropy (bits)', color='#d1d5db', fontsize=12)
ax3.set_title('2-Qubit Bipartition Entropies', color='#f9a8d4', fontsize=14, fontweight='bold')
ax3.set_xticks(range(len(pair_labels)))
ax3.set_xticklabels(pair_labels, fontsize=10)
ax3.set_facecolor('#0f172a')
ax3.tick_params(colors='#9ca3af')
ax3.grid(True, alpha=0.15, color='#6b7280', axis='y')
for spine in ax3.spines.values():
    spine.set_color('#334155')

# Panel 4: Probability distribution
ax4 = axes[1, 1]
probs = np.abs(sv.data)**2
nonzero_mask = probs > 1e-12
states_labels = [f'|{i:04b}>' for i in range(16)]
colors_bar = ['#f472b6' if p > 0.01 else '#334155' for p in probs]
ax4.bar(range(16), probs, color=colors_bar, edgecolor='#be185d', alpha=0.8)
ax4.set_xlabel('Computational basis state', color='#d1d5db', fontsize=12)
ax4.set_ylabel('Probability', color='#d1d5db', fontsize=12)
ax4.set_title('Measurement Probabilities', color='#f9a8d4', fontsize=14, fontweight='bold')
ax4.set_xticks(range(16))
ax4.set_xticklabels(states_labels, fontsize=7, rotation=45, ha='right')
ax4.set_facecolor('#0f172a')
ax4.tick_params(colors='#9ca3af')
ax4.grid(True, alpha=0.15, color='#6b7280', axis='y')
for spine in ax4.spines.values():
    spine.set_color('#334155')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()

print()
print("The tetrahedron graph state has uniform entanglement across all qubits,")
print("reflecting the symmetry of the complete graph K4.")
print("This is a toy model for a spin network vertex in loop quantum gravity.")
print("The entanglement structure encodes the quantum geometry of a tetrahedron.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Share:X Reddit LinkedIn

Quantum Gravity: Lattice and Information Approaches

Table of Contents

Canonical Quantum Gravity

The ADM Hamiltonian

The Wheeler-DeWitt Equation

Lab: Toy Mini-Superspace WDW Equation

Mini-Superspace Wheeler-DeWitt Equation (1D Quantum Cosmology)

Loop Quantum Gravity

Holonomies and Wilson Loops

Area Quantization

Lab: LQG Area Spectrum

LQG: Discrete Area Spectrum for j = 1/2 to j = 5

Spin Foam Models

Vertex Amplitudes and 15j-Symbols

6j-Symbols: Building Blocks

Lab: 6j Symbols and Spin Foam Graph

Spin Foam Models: 6j Symbols & Foam Graph Visualization

AdS/CFT Correspondence

Central Charge and Brown-Henneaux

Lab: BTZ Black Hole Entropy vs Cardy Formula

BTZ Black Hole Entropy vs Cardy Formula

Holographic Entanglement Entropy

RT Surface in AdS3

Lab: RT Surface in AdS3

RT Surface in AdS3: Geodesic Computation & Entanglement Entropy

Tensor Networks

MERA and Holography

Matrix Product States (MPS)

Lab: MPS for 6 Qubits

Matrix Product State: 6-Qubit Entanglement Analysis

Black Hole Information Paradox

The Page Curve

The Island Formula

Lab: The Page Curve

The Page Curve: Island Phase Transition

Quantum Computing for Quantum Gravity

Graph States and Spin Networks

Lab: Spin Network on 4 Qubits

Spin Network State on 4 Qubits: Tetrahedron Graph State

RT Surface in AdS₃

Lab: RT Surface in AdS₃