Analog Gravity: Lensing by Density Gradients

Visualizing light bending in a variable density medium.

Objective

To demonstrate that gravitational lensing—the bending of light by massive objects—can be reproduced as a hydrodynamic phenomenon. In the Dynamic Background Hypothesis, gravity is not a geometric curvature of empty space, but a refraction effect caused by variations in the vacuum density.

Methodology

We simulate the Gross-Pitaevskii equation with a high-momentum wave packet (representing a photon) passing near a region of low density.

1. The Mass: We introduce a static repulsive potential \(V_{ext}\) at the center. This creates a “hole” in the superfluid density (\(\rho < \rho_0\)).
2. The Metric: Since the speed of sound in the medium depends on density (\(c_s \propto \sqrt{\rho}\)), the low-density region acts as a medium with a high refractive index.
3. The Lensing: According to Huygens’ principle, the part of the wavefront closest to the mass travels slower, causing the trajectory to bend towards the center

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML

# 1. CONFIGURATION
N = 256             # High resolution for wave packet detail
L = 40.0            # Large domain
dt = 0.01           # Time step
g = 1.0             # Interaction strength
rho_0 = 1.0         # Base density

x = np.linspace(-L/2, L/2, N)
y = np.linspace(-L/2, L/2, N)
X, Y = np.meshgrid(x, y)
dx = x[1] - x[0]

# k-space
k = 2 * np.pi * np.fft.fftfreq(N, d=dx)
KX, KY = np.meshgrid(k, k)
K2 = KX**2 + KY**2

# 2. THE "MASSIVE OBJECT" (EXTERNAL POTENTIAL)
# A repulsive potential creates a low-density region (The "Star")
R_center = np.sqrt(X**2 + Y**2)
V_star_amp = 1.5 * g * rho_0
V_star_width = 4.0
V_ext = V_star_amp * np.exp(-(R_center / V_star_width)**2)

# Absorbing boundary conditions (soft trap)
V_trap = 0.1 * (R_center / (0.45*L))**20
V_total = V_ext + V_trap

# 3. BACKGROUND PREPARATION (THE METRIC)
# Imaginary time evolution to create the static curved background
Psi_bg = np.sqrt(rho_0) * np.ones((N, N), dtype=complex)
dt_im = 0.01

# Relaxation loop (silent for web output)
for i in range(300):
    rho = np.abs(Psi_bg)**2
    Psi_bg *= np.exp(-(V_total + g*rho) * dt_im)
    Psi_bg *= np.sqrt(rho_0) / np.max(np.abs(Psi_bg[-1,:])) # Local renormalization
    Psi_k = np.fft.fft2(Psi_bg)
    Psi_k *= np.exp(-(K2/2)*dt_im)
    Psi_bg = np.fft.ifft2(Psi_k)

# 4. PHOTON INJECTION
# High momentum Gaussian packet
k_photon = 3.0      
y_impact = 5.0      # Impact parameter
x_start = -15.0     
width_photon = 1.5  

Packet = 0.3 * np.sqrt(rho_0) * \
         np.exp(-((X - x_start)**2 + (Y - y_impact)**2) / (2*width_photon**2)) * \
         np.exp(1j * k_photon * X)

Psi = Psi_bg + Packet

# 5. EVOLUTION & VISUALIZATION
U_kin = np.exp(-1j * (K2/2) * dt)

def step(psi_in):
    psi_mod = psi_in * np.exp(-1j * (V_total + g*np.abs(psi_in)**2) * (dt/2))
    psi_k = np.fft.fft2(psi_mod)
    psi_k *= U_kin
    psi_mod = np.fft.ifft2(psi_k)
    psi_out = psi_mod * np.exp(-1j * (V_total + g*np.abs(psi_mod)**2) * (dt/2))
    return psi_out

# Plot setup
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))
plt.close() # Prevent static plot output

# Panel 1: The Metric (Background Density)
im1 = ax1.imshow(np.abs(Psi_bg)**2, extent=[-L/2, L/2, -L/2, L/2], 
                 cmap='gray', origin='lower')
ax1.set_title("Background Metric (Density)")
ax1.set_xlabel("Dark region = Slow light speed")

# Panel 2: The Photon (Perturbation)
diff = np.abs(Psi)**2 - np.abs(Psi_bg)**2
im2 = ax2.imshow(diff, extent=[-L/2, L/2, -L/2, L/2], 
                 cmap='plasma', origin='lower', vmin=0, vmax=0.1)
ax2.set_title("Photon Trajectory")
ax2.axhline(y=y_impact, color='white', linestyle='--', alpha=0.3, label="Newtonian Path")
ax2.legend(loc='upper right')

def update(frame):
    global Psi
    for _ in range(4):
        Psi = step(Psi)
    
    # Visualize only the photon (subtract background)
    photon_density = np.abs(Psi)**2 - np.abs(Psi_bg)**2
    im2.set_data(photon_density)
    return im1, im2

ani = animation.FuncAnimation(fig, update, frames=150, interval=30, blit=True)
display(HTML(ani.to_html5_video()))