Quantum Simulator

Real-time quantum physics simulations rendered with OpenGL. Explore tunneling, interference, scattering, and other quantum phenomena through interactive visualizations.

Build

1. Generate GLAD (one-time setup)

GLAD is an OpenGL function loader that must be generated before building.

1. Go to https://glad.dav1d.de/
2. Set the following options:

• Language: C
• Specification: OpenGL
• API: gl 4.5
• Profile: Core
• Generate a loader: yes

3. Click Generate and download the zip
4. Place the output files into the repo:

   external/glad/src/glad.c
   external/glad/include/glad/glad.h
   external/glad/include/KHR/khrplatform.h
   

2. Configure and build

cmake -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build --parallel

3. Run

./build/quantum_sim

4. Run tests

cd build && ctest --output-on-failure

---

Simulations

#	Name	Description
1	Quantum Tunneling	Gaussian wave packet vs. rectangular barrier
2	Double Slit	2D interference pattern from wave-particle duality
3	Harmonic Oscillator	Coherent state oscillating in a parabolic well
4	Particle Scattering	Two charged packets interacting via Coulomb potential
5	Photon Interference	Superposition and beat patterns from two photon packets
6	3D Harmonic Oscillator	3D coherent state in an isotropic harmonic trap (xy projection)
7	Spin-1/2	Two-component spinor in a magnetic field — Rabi oscillations
8	Entanglement	Two-particle joint density in configuration space with linear entropy
9	3D Oscillator (Volume)	Coherent state in a 3D harmonic trap — volumetric ray marching

Simulation Details

All simulations use the split-step Fourier method: each time step alternates between a half-step in position space (applying the potential phase exp(-i V(x) dt/2ℏ)) and a full step in momentum space (applying the kinetic phase exp(-i ℏk²dt/2m)). This keeps the evolution unitary and conserves probability automatically.

Every slider change resets the simulation to a fresh initial state.

---

1 — Quantum Tunneling

A Gaussian wave packet starts to the left of a rectangular potential barrier and propagates toward it. In classical mechanics a particle with energy below the barrier height cannot cross. Quantum mechanically, the wave function penetrates the barrier and a fraction emerges on the far side — this is quantum tunneling.

The display shows the probability density |Ψ(x)|² as a blue line. You will see the packet split into a reflected part (moving back left) and a transmitted part (continuing right). The ratio depends exponentially on barrier width and the gap between barrier height and packet energy.

Slider	Range	Effect
Height	0.1 – 20	Barrier potential V₀. Higher → exponentially less transmission. Near V₀ ≈ k₀²/2 is the classical threshold.
Width	0.1 – 4	Barrier thickness. Tunneling probability falls off as exp(-2κL) where κ = √(2m(V₀−E))/ℏ.
Momentum	0.5 – 15	Initial wave-number k₀ of the packet (energy ∝ k₀²). Push it above the barrier height to watch full classical transmission with no reflection.
Sigma	0.1 – 2	Spatial width of the Gaussian. Wider σ → narrower momentum spread → cleaner single-energy behavior. Narrow σ → broad momentum spread, partial energy components both above and below barrier.
dt	0.0001 – 0.01	Time step. Smaller = more accurate; larger = faster but may introduce drift in total probability.

What to look for: set Height = 8, Momentum = 4 to see a clear split. Increase Momentum past ~4 to watch the transmitted fraction grow. Increase Width to see tunneling suppressed even for the same energy.

---

2 — Double Slit

A 2D Gaussian wave packet with a large transverse spread (quasi-plane-wave) propagates in the +x direction and strikes a wall containing two narrow slits. Each slit acts as an independent coherent source. The two transmitted waves overlap and produce the canonical interference fringe pattern: bright bands where path lengths differ by a whole wavelength, dark bands where they differ by a half wavelength.

The display is a 2D heatmap (viridis colormap). A gamma curve (γ = 0.4) is applied to the density so that faint outer fringes — at roughly 1% of peak intensity — are lifted to ~16% on the color scale and remain visible alongside the bright central beam.

Slider	Range	Effect
Slit Width	0.05 – 2	Width of each individual opening. Narrower → more diffraction → fringes spread wider. Wider → the single-slit envelope narrows, concentrating intensity in the center.
Slit Spacing	0.5 – 6	Center-to-center distance d. Fringe spacing ≈ λ/d, so increasing d packs more fringes into the same angular region.
Wall Height	10 – 500	Potential of the opaque wall. Higher → more complete blocking between and outside the slits → cleaner two-source pattern.
dt	0.0001 – 0.01	Time step.

What to look for: watch the wave reach the wall, pass through the slits, and fan out into stripes. Slide Slit Spacing up to compress the fringes, or reduce Slit Width to widen the diffraction envelope and reveal more orders.

---

3 — Harmonic Oscillator

A Gaussian wave packet (a coherent state) is placed off-center inside a parabolic potential well V(x) = ½ω²x². A coherent state is the quantum state that most closely resembles a classical oscillator: it oscillates back and forth without spreading, with period T = 2π/ω.

The display shows |Ψ(x)|² as a red line. The background draws the parabolic potential well to scale.

Slider	Range	Effect
Omega	0.5 – 10	Angular frequency ω of the well. Higher → steeper walls → faster oscillation. Period T = 2π/ω, so ω = 2 gives T ≈ 3.14 time units.
x0	−5 – 5	Initial displacement from center. Sets the classical amplitude. x0 = 0 places the packet at the ground-state center — it stays still.
dt	0.0001 – 0.01	Time step.

What to look for: set x0 to ±3 and watch the packet oscillate like a pendulum. Increase Omega to speed up the oscillation. The packet shape stays approximately Gaussian throughout — that is the defining property of a coherent state.

---

4 — Particle Scattering

Two same-charge wave packets (charge +1 each) start on opposite sides of the grid and move toward each other with equal and opposite momenta. They interact through a mean-field Coulomb potential: each packet evolves in the electric field generated by the other packet's probability density distribution. The repulsive interaction deflects and scatters both packets.

The display draws each particle's |Ψ(x)|² as a separate colored line strip (cyan = particle 1, orange = particle 2).

Slider	Range	Effect
Momentum	0.5 – 10	Initial momentum magnitude for both packets (they approach each other). Low momentum → Coulomb repulsion dominates, packets bounce back. High momentum → packets tunnel through each other.
Separation	2 – 12	Initial center-to-center distance. Smaller → interaction begins sooner with higher Coulomb energy; larger → packets propagate freely longer before repulsion kicks in.
dt	0.0001 – 0.01	Time step.

What to look for: at low momentum the packets slow down, deform, and reverse direction — quantum reflection from the Coulomb barrier. At high momentum they pass through each other, which has no classical analogue. Intermediate values show partial transmission and broad spread-out density after the collision.

---

5 — Photon Interference

Two plane-wave components with different wave numbers k₁ and k₂ are superposed inside a shared Gaussian envelope and then evolved freely (V = 0). The result is a beat pattern: rapid carrier oscillations at frequency (k₁ + k₂)/2 amplitude-modulated by a slow envelope at |k₁ − k₂|/2. Both components travel at their own phase velocities so the pattern evolves as they drift apart.

The display draws the total density |Ψ|² (yellow) plus the real and imaginary parts of Ψ as separate overlaid curves.

Slider	Range	Effect
k1	1 – 15	Wave number of the first component (momentum p₁ = ℏk₁). Sets its phase velocity and wavelength.
k2	1 – 15	Wave number of the second component. When k₁ = k₂ there are no beats — single pure frequency. As |k₁ − k₂| grows, beat nodes become more pronounced and tightly spaced.
dt	0.0001 – 0.01	Time step.

What to look for: set k1 = 5, k2 = 6 for gentle beating with a slowly drifting envelope. Push them to k1 = 3, k2 = 12 to see rapid high-frequency oscillations inside a broad envelope. The Re(Ψ) and Im(Ψ) traces reveal the carrier oscillation directly.

---

6 — 3D Harmonic Oscillator

A 3D Gaussian wave packet (coherent state) placed off-center in an isotropic harmonic potential V(x,y,z) = ½ω²(x²+y²+z²). In 3D the packet oscillates simultaneously along all three axes with the same period T = 2π/ω.

The display is a heatmap of the xy-plane marginal density: ρ(x,y) = ∫|Ψ(x,y,z)|² dz — the probability density projected onto the horizontal plane by integrating out the z direction.

Slider	Range	Effect
Omega	0.5 – 5	Angular frequency ω. Period T = 2π/ω.
x0, y0, z0	−6 – 6	Initial displacement along each axis. Setting only x0 gives purely horizontal oscillation visible in the projection; nonzero y0 or z0 adds perpendicular components.
dt	0.0001 – 0.02	Time step.

What to look for: set x0 = 3, y0 = 0, z0 = 0 and watch the projected blob oscillate left-right. Set x0 = 2, y0 = 2 to see the packet trace an ellipse in the xy projection. Because the potential is isotropic, all three periods are equal and the motion is periodic.

---

7 — Spin-1/2

A Gaussian wave packet representing a spin-1/2 particle propagating in a uniform magnetic field B = (Bx, 0, Bz). The field has two distinct effects:

• Bz (ωZ — Larmor frequency): splits the energy of ↑ and ↓ by ℏωZ. Each component acquires a different phase at the same rate; no population transfer.
• Bx (ωX — Rabi frequency): drives coherent spin-flip transitions. Population oscillates between ↑ and ↓ at the Rabi frequency ωX with period T = 2π/ωX.

The blue line shows |Ψ↑(x)|² and the red line shows |Ψ↓(x)|². The ImGui panel displays P(↑), P(↓) and the spin-z expectation value ⟨Sz⟩ in real time.

Slider	Range	Effect
x0	−8 – 8	Initial packet center.
k0	−10 – 10	Initial wave number (group velocity).
Sigma	0.2 – 2	Packet width.
omegaZ	−20 – 20	Larmor frequency from Bz. Adds a relative phase between ↑ and ↓; non-zero omegaZ combined with omegaX shifts the oscillation axis.
omegaX	0 – 20	Rabi frequency from Bx. Drives spin-flip oscillations. Period T = 2π/ωX.
Spin θ	0 – π	Initial spin polar angle. 0 = pure ↑, π = pure ↓, π/2 = balanced superposition.
dt	0.0001 – 0.02	Time step.

What to look for: set omegaZ = 0, omegaX = 5 and watch probability oscillate between the blue and red traces at period T ≈ 1.26. Increase omegaZ to tilt the rotation axis — the oscillation no longer reaches full inversion. Set θ = π/2 for a balanced starting superposition.

---

8 — Entanglement

Two distinguishable particles confined in a shared harmonic trap with a repulsive soft-core Coulomb interaction between them. The state is the two-particle wave function Ψ(x₁, x₂) in configuration space — a 2D grid where the horizontal axis is particle-1 position and the vertical axis is particle-2 position.

An unentangled (product) state Ψ = φ₁(x₁)·φ₂(x₂) appears as a simple blob with no correlations. As the interaction builds entanglement, the density develops structure: repulsion pushes the particles apart so high density appears away from the x₁ = x₂ diagonal.

The linear entropy S_L = 1 − Tr[ρ₁²] ∈ [0, 1) quantifies entanglement:

• S_L = 0: product state (no entanglement)
• S_L > 0: entangled — the two particles cannot be described independently

Slider	Range	Effect
Omega	0.2 – 4	Harmonic trap frequency ω. Tighter trap increases interaction energy.
Coupling	0 – 10	Repulsion strength g. Higher → stronger Coulomb-like repulsion → faster entanglement buildup.
Separation	0.5 – 4	Initial distance
dt	0.0001 – 0.02	Time step.

What to look for: start with Coupling = 0 — the blob oscillates without distortion. Increase Coupling gradually and watch S_L grow as the density develops an anti-diagonal tilt. With high coupling the probability concentrates in two separated blobs (anti-bunching), reflecting particle repulsion.

---

9 — 3D Oscillator (Volume)

The same 3D Gaussian coherent state as simulation 6, but rendered by orthographic ray marching through the full volumetric probability density |Ψ(x,y,z)|². Each screen pixel casts a parallel ray from a fixed camera position (looking from direction (1, 0.7, 1)) into the grid, sampling the 3D density texture along the ray and accumulating color and opacity front-to-back using the viridis colormap with gamma compression.

This shows the true 3D shape of the wave packet — a blob whose center traces an elliptical path in space — rather than a planar projection.

Slider	Range	Effect
Omega	0.5 – 5	Angular frequency ω. Period T = 2π/ω.
x0, y0, z0	−6 – 6	Initial displacement along each axis.
dt	0.0001 – 0.02	Time step.
Extinction	0.5 – 30	Opacity accumulation coefficient. Low → transparent cloud; high → opaque solid core.
View Scale	0.5 – 3	Camera frustum half-extent. Increase to zoom out; decrease to zoom in on the packet.
Steps	32 – 256	Number of ray march samples per pixel. More → smoother but slower.

What to look for: set x0 = 3, y0 = 0, z0 = 0 and watch the glowing blob oscillate left-right. Set x0 = 2, y0 = 2 to see it trace a circular path. Use high Extinction to see the dense core; reduce it to reveal the faint tails. Compare with simulation 6 to see what information the marginal projection loses.

---

Architecture

The project follows a layered architecture with strict dependency direction:

┌──────────────────────────────────────────┐
│  app/         SimulationRunner, UI,      │
│               concrete simulations       │
├──────────────────────────────────────────┤
│  infra/       OpenGL rendering (RAII),   │
│               GLFW window, input, ImGui  │
├──────────────────────────────────────────┤
│  domain/      Pure physics: wave funcs,  │
│               potentials, particles,     │
│               interactions (no GL deps)  │
└──────────────────────────────────────────┘

The domain layer has zero infrastructure dependencies. It can be tested, benchmarked, and reused independently of any rendering framework.

GLSL shaders live in assets/shaders/ and are loaded from disk at runtime via ShaderProgram::fromFiles(). The SHADER_DIR compile definition (set in CMakeLists.txt) resolves to the absolute path of that directory at build time.

File	Used by	Purpose
line.vert	WaveFunctionRenderer, ParticleRenderer	Vertex stage for all 1D line-strip draws (density, Re/Im components, potential overlay, particle densities)
line.frag	WaveFunctionRenderer, ParticleRenderer	Solid-color fill for line strips (uColor uniform)
heatmap.vert	WaveFunctionRenderer	Fullscreen-quad vertex stage for the 2D heatmap pass
heatmap.frag	WaveFunctionRenderer	Viridis colormap with γ = 0.4 compression for 2D probability density

Prerequisites

• C++23 compiler (GCC 13+, Clang 17+, or MSVC 17.8+)
• CMake 3.25+
• OpenGL 4.5 capable GPU + drivers
• GLAD (see external/glad/README.md for generation instructions)

All other dependencies (GLFW, GLM, Dear ImGui, Google Test) are fetched automatically via CMake FetchContent.

Controls

Key	Action
1–9	Switch simulation
Space	Pause / Resume
R	Reset current simulation
Esc	Quit

Parameters are adjustable in real time via the ImGui panel on the right.

Coding Standards

See AGENTS.md for the full set of rules. Key points:

• C++23 standard, CMake-only build system
• Smart pointers only — no raw new / malloc
• RAII for all resources (GPU buffers, shaders, windows)
• Domain logic free of infrastructure concerns
• Google Test for all public APIs
• clang-format + clang-tidy enforced in CI

Future Work

• Replace naive DFT with FFTW / cuFFT for real-time performance
• GPU-accelerated evolution via compute shaders (.comp) or CUDA (.cu)
• 3D wave function rendering with ray marching
• Spin-1/2 particle simulations (Pauli equation)
• Entanglement visualization (two-particle joint probability)
• Recording / export of simulation frames