Quantum Tunneling in Your Browser: Solving the 2D Schrödinger Equation in WebGL

After simulating classical continuum mechanics with incompressible Navier–Stokes fluids and Turing reaction–diffusion morphogenesis, we can push the same ping-pong WebGL Framebuffer architecture into the quantum regime: evolving a complex probability amplitude $\psi(\vec x, t) \in \mathbb{C}$ governed by the 2D time-dependent Schrödinger equation.

Where a reaction–diffusion PDE is purely dissipative and damps out oscillations, the Schrödinger equation is unitary: total probability $\int |\psi(\vec x, t)|^2\,d^2x = 1$ must be conserved forever. Naive forward-Euler time stepping fails catastrophically here – its amplification factor $|\lambda| = \sqrt{1 + (\Delta t E)^2} \gt 1$ pumps artificial energy into every mode until the wavefunction explodes to infinity in a dozen frames. Using Visscher's (1991) staggered symplectic leapfrog integrator, we can evolve the real and imaginary parts $(\psi_R, \psi_I)$ in alternating half-steps with exact unitarity ($|\lambda| \equiv 1$), running 32 shader passes per frame at 60 fps.

0. Anatomy of a GPU Quantum Wavepacket Solver

Working in natural grid units $\hbar = m = \Delta x = 1$, each texel of our floating-point RGBA16F FBO stores the real and imaginary components of the matter wave $\psi(\vec x, t) = \psi_R(\vec x, t) + i\,\psi_I(\vec x, t + \tfrac{1}{2}\Delta t)$:

Governing PDE

$i\partial_t \psi = \hat H\psi = \big(-\tfrac{1}{2}\nabla^2 + V(\vec x)\big)\psi$, coupling the curvature of $\psi_R$ to the time derivative of $\psi_I$ and vice versa.

Wavepacket State

Coherent minimum-uncertainty Gaussian $\psi(\vec x, 0) = e^{-|\vec x - \vec x_0|^2 / 2\sigma^2} e^{i\vec k_0 \cdot (\vec x - \vec x_0)}$ with carrier kinetic energy $E_0 = \tfrac{1}{2}|\vec k_0|^2$.

Symplectic Leapfrog

Alternating explicit updates to $\psi_R^n$ and $\psi_I^{n+1/2}$ keep the time-evolution eigenvalues strictly on the complex unit circle $|\lambda| = 1$ for all $\Delta t \lt 2 / E_{\max}$.

Domain Coloring

Maps complex phase $\arg\psi = \mathrm{atan2}(\psi_I, \psi_R) \in [-\pi, \pi]$ to the HSV hue wheel and probability density $|\psi|^2$ to luminance.

Evanescent Tunneling

Inside a classically forbidden barrier ($V_0 \gt E_0$), the wavenumber becomes imaginary ($k = i\kappa = i\sqrt{2(V_0 - E_0)}$), decaying exponentially as $e^{-\kappa x}$ before emerging on the far side.

1. Why Forward Euler Explodes – and How Visscher's Symplectic Leapfrog Fixes It

In natural units ($\hbar = m = 1$), the 2D time-dependent Schrödinger equation is a first-order complex PDE that splits cleanly into two coupled real equations for $\psi(\vec x, t) = \psi_R(\vec x, t) + i\,\psi_I(\vec x, t)$:

Because $\dot\psi_R = +\hat H\psi_I$ and $\dot\psi_I = -\hat H\psi_R$ have the exact structure of Hamilton's canonical equations, applying forward Euler to both components at the same timestamp $t_n$ spirals outward in phase space ($|\lambda|^2 = 1 + E^2\Delta t^2 \gt 1$). Instead, Visscher (1991) staggered $\psi_R$ at integer steps $t_n$ and $\psi_I$ at half-integer steps $t_{n+1/2}$, updating them in two alternating ping-pong passes per substep:

// Exact unit-grid isotropic 9-point Laplacian: (10/3) * (0.20*cardinal + 0.05*diagonal - center)
vec2 lap = (10.0 / 3.0) * (0.20 * (N + S + E + W) + 0.05 * (NE + NW + SE + SW) - C);
float V  = evalPotential(vUv);

float psiR = C.r;
float psiI = C.g;

if (uPassType == 0) {
  // Pass A: advance Real part psi_R^(n+1) using H * psi_I^(n+1/2)
  float hPsiI = -0.5 * lap.g + V * psiI;
  psiR = (psiR + uDt * hPsiI) * damp;
} else {
  // Pass B: advance Imag part psi_I^(n+3/2) using H * psi_R^(n+1)
  float hPsiR = -0.5 * lap.r + V * psiR;
  psiI = (psiI - uDt * hPsiR) * damp;
}
gl_FragColor = vec4(psiR, psiI, 0.0, 1.0);

2. Coherent Wavepackets, Group Velocity, and Evanescent Tunneling

When you click or drag on the canvas, the shader injects a minimum-uncertainty 2D Gaussian wavepacket centered at $\vec x_0$ with spatial width $\sigma$ and carrier wavevector $\vec k_0 = (k_x, k_y)$:

What happens when this wavepacket hits a potential barrier of height $V_0$ and thickness $L$ (click the Tunneling Wall preset in Controls)? Classical mechanics says that if $V_0 \gt E_0$, 100% of the particles bounce back. Quantum mechanics disagrees: inside the barrier, the effective kinetic energy $E_0 - V_0$ is negative, turning the oscillatory wave $e^{ikx}$ into a real exponential decay $e^{-\kappa x}$ with $\kappa = \sqrt{2(V_0 - E_0)}$. If the barrier thickness $L$ is comparable to $1/\kappa$, the evanescent tail reaches the far face with non-zero amplitude and reconverts into a propagating wave:

3. Visualizing a Complex Field: Phase Wheels vs. Probability Density

A quantum wavefunction has two degrees of freedom at every point in space: amplitude $|\psi(\vec x, t)| = \sqrt{\psi_R^2 + \psi_I^2}$ and phase angle $\varphi(\vec x, t) = \mathrm{atan2}(\psi_I, \psi_R) \in [-\pi, \pi]$. Plotting $|\psi|^2$ alone hides the momentum $\vec k = \nabla\varphi$ that tells you which way the wave is moving. Using complex domain coloring, the display pass encodes both quantities at once:

vec2 psi  = texture2D(uState, vUv).rg;
float prob = dot(psi, psi);                  // Born probability density |psi|^2 = psi_R^2 + psi_I^2
float amp  = sqrt(probSmooth);

// Map complex phase angle arg(psi) in [-pi, pi] to [0, 1] HSV hue cycle
float phase = atan(psi.g, psi.r);
float hue   = phase / 6.2831853 + 0.5;
float val   = smoothstep(0.0, 0.85, pow(amp * 1.45, 0.78));
vec3 waveCol = hsv2rgb(vec3(hue, 0.82, val));

4. The Staggered Leapfrog GPU Pipeline

Each animation frame runs 16 symplectic leapfrog cycles (32 alternating Pass A / Pass B draw calls) on the RGBA16F wavefunction FBO before rendering the domain-colored phase wheel and potential barrier contours:

5. Four Quantum Scattering Experiments in One Shader

Each preset button in the control panel swaps the analytical potential $V(\vec x)$ evaluated inside the leapfrog stencil:

Switch to Resonant Ring or Bragg Lattice in the controls above, then click-and-drag at different angles on the canvas to fire aimed wavepackets into the target – watching a matter wave diffract around a crystal lattice in real time makes wave-particle duality tangible at 60 frames per second.