MZIMesh.jl

Simulator for programmable Mach–Zehnder interferometer meshes. Pure Julia, no external dependencies beyond the standard library.

See docs/mesh-simulator-plan.md for the full planning document. Current status: Phases 0–4 complete (scaffolding, mesh core, Clements decomposition, Fock backend, Makie visualization), plus photon-counting statistics (photon_counting, photon_moments) and an upgraded interactive dashboard (N up to 8, pattern-based inputs, automatic fallback from the exact joint distribution to per-mode moments). Phase 5 (Keifer topology) pending on the open question of its decomposition algorithm. Gaussian boson sampling (the hafnian/covariance-matrix regime) is future work — the current Fock backend implements Fock-input (Aaronson–Arkhipov) sampling via permanents.

Architecture invariant

The mesh layer is a pure function producing U(θ⃗, φ⃗) ∈ U(N). Backends consume U with their own state representations. Nothing about photons or Fock spaces lives in the mesh types (src/mesh.jl, src/decompose.jl); the Fock backend (src/fock.jl, src/permanent.jl) sees only a unitary matrix. A future Gaussian/CV backend attaches at the same interface.

Conventions

• MZI block: B(θ, φ) = [e^{iφ}cosθ -sinθ; e^{iφ}sinθ cosθ], acting on adjacent modes (m, m+1). Reflectivity is cos²θ, so a 50:50 splitter is θ = π/4 (papers with the half-angle convention T ∝ cos(θ/2) call the same device θ = π/2). Full transmission range is θ ∈ [0, π/2].
• Mesh function: compile(mesh, phases) = D · T_K ⋯ T_1 with placements in application order (first entry = first crossing light encounters) and D = Diagonal(exp.(iδ)) the output phase screen.
• Fock amplitudes: ⟨T|φ(U)|S⟩ = Perm(U_{T,S}) / √(∏ sᵢ! ∏ tⱼ!), permanent via Ryser's formula.

Quick start

using MZIMesh

# Forward: phases → unitary
mesh = clements(4)                       # or reck(4)
ps = PhaseSettings(fill(π/4, nmzi(mesh)), zeros(nmzi(mesh)), zeros(4))
U = compile(mesh, ps)

# Inverse: unitary → phases (Clements 2016 nulling)
mesh2, ps2 = decompose(U)                # compile(mesh2, ps2) ≈ U to ~1e-13

# Fock backend: Hong–Ou–Mandel
bs = compile(clements(2), PhaseSettings([π/4], [0.0], zeros(2)))
states, probs = output_distribution(bs, [1, 1])   # P(1,1) = 0, P(2,0) = P(0,2) = 1/2

# Photon counting: exact output-state distribution + full counting statistics
# (per-mode means/variances/marginals, plus the per-source decomposition of
# the means: sources[j,i] = mean photons at output j that came from input mode i)
stats = photon_counting(bs, [1, 1])
stats.mean          # [1.0, 1.0]
stats.var           # [1.0, 1.0] — bunching: twice the distinguishable value
stats.sources       # [0.5 0.5; 0.5 0.5]
mean_photons(bs, [1, 1])   # means alone, closed form (no permanents)

# Exact per-mode moments for any photon number, O(N³), no permanents —
# scales where photon_counting's joint distribution would explode
pm = photon_moments(bs, [10, 10])
pm.mean   # [10.0, 10.0]
pm.var    # [55.0, 55.0] — twin-Fock bunching m(m+1)/2, exact for any photon number

# Mesh bookkeeping: dof(mesh) = 2K+N phase parameters (= N² = dim U(N) for
# Clements/Reck), depth(mesh) = scheduled layer count
dof(clements(4)), depth(clements(4)), depth(reck(4))   # (16, 4, 5)

# Visualization (loads the MZIMeshMakieExt package extension)
using GLMakie            # interactive; or CairoMakie for static output
mesh_figure(mesh, ps)                 # annotated mesh diagram
distribution_figure(bs, [1, 1])       # output-distribution bar chart
dashboard()                           # reactive: topology/N(2–8)/pattern+line+photons-per-line/preset controls
                                       # + per-MZI θ/φ tuning sliders (selected MZI highlighted in the
                                       # diagram) + stats panel; switches between an exact output-distribution
                                       # bar chart and per-mode mean±σ moments once the Fock space is too big

Tests

julia --project=. -e 'using Pkg; Pkg.test()'

Covers: unitarity for random phases on both topologies (N ≤ 10), MZI counts, hand-checked N = 2, 3 cases, Haar-random round-trip decomposition to 1e-12, permanent vs. naive reference, HOM dip, balanced-tritter statistics, classical-routing limit, distribution normalization, and photon-counting statistics (means vs. the interference-free closed form, source-contribution sum rules, marginals, HOM super-classical variance, dof/depth). Also photon_moments: cross-checked against the exact joint distribution (mean, variance, and full covariance matrix), photon-number conservation sum rules (Σⱼ⟨n̂ⱼ⟩ = n, each cov row sums to 0) for a 40-photon case no joint distribution could handle, and the twin-Fock [10, 10] variance through a 50:50 beamsplitter.