Simulation engine reference · v0.8
The simulation works in the tensor-product space spin ⊗ motion, with basis ordering:
The total Hilbert space dimension is 2 × Nmax. All operators are represented as complex matrices of this dimension. Complex numbers are stored as pairs of Float64 values [re, im] in a flat array of length 2 × (2Nmax)².
The truncation at Nmax is the primary source of systematic error. After each simulation run, the population in the top 3 Fock states (both spin components) is measured:
Green (OK): boundary population < 0.01% — results converged.
Amber (Marginal): 0.01–1% — possible truncation artefacts.
Red (Unreliable): > 1% — results not converged. Recommended Nmax ≥ ⟨n⟩ + 5√⟨n⟩ + 10.
In the rotating frame at the laser frequency, the effective Hamiltonian for one stroboscopic pulse at detuning δ₀ is:
where:
The operator C is computed exactly in the Fock basis via matrix exponentiation of the position quadrature operator X = a + a†. No Lamb–Dicke approximation or sideband truncation is applied. The matrix exponential uses a Taylor(12) series with scaling and squaring.
In the 2Nmax × 2Nmax matrix:
| Block | Rows | Cols | Content |
|---|---|---|---|
| ↓↓ | 0…N−1 | 0…N−1 | −δ₀/2 on diagonal |
| ↑↑ | N…2N−1 | N…2N−1 | +δ₀/2 on diagonal |
| ↓↑ (σ₋) | 0…N−1 | N…2N−1 | (Ω/2) · C |
| ↑↓ (σ₊) | N…2N−1 | 0…N−1 | (Ω/2) · C† |
The protocol applies Np pulses, each separated by one motional period Tm = 2π/ωm. The per-pulse unitary is:
The pulse duration δt is derived from the requirement that the total accumulated rotation equals π/2:
where Ωeff = Ω · exp(−η²/2) is the Debye–Waller suppressed carrier Rabi rate.
Changing Np does not change the total rotation — it changes the per-pulse rotation angle θ = π/(2Np) and the duty cycle δt/Tm. The total sequence duration is Np × Tm.
Between pulses, the ion evolves freely for one motional period. In the Fock basis:
This is exact for integer Fock states: exp(−i2πn) = 1 for all n ∈ {0, 1, 2, …}. The stroboscopic condition is therefore built into the Fock-basis representation — the inter-pulse free evolution is the identity operator. The full Np-pulse propagator is simply UNp.
For the default parameters (η = 0.397, Ω/(2π) = 0.300 MHz, ωm/(2π) = 1.30 MHz):
| Np | θ/pulse | δt [ns] | Duty cycle | Tseq [μs] |
|---|---|---|---|---|
| 10 | 9.0° | 90 | 11.7% | 7.7 |
| 22 | 4.1° | 41 | 5.3% | 16.9 |
| 50 | 1.8° | 18 | 2.3% | 38.5 |
The displaced coherent state D(α)|0⟩ is computed in the Fock basis as:
For real α (convention used throughout), all coefficients are real.
After the pulse sequence, the reduced spin density matrix is obtained by tracing over the motional degrees of freedom:
From the 2×2 spin density matrix:
| Observable | Formula |
|---|---|
| ⟨σx⟩ | 2 Re(ρ↓↑) |
| ⟨σy⟩ | −2 Im(ρ↓↑) |
| ⟨σz⟩ | ρ↑↑ − ρ↓↓ |
| Coherence | √(⟨σx⟩² + ⟨σy⟩² + ⟨σz⟩²) |
| Entropy S | −λ₊ log₂ λ₊ − λ₋ log₂ λ₋ |
where λ± = (1 ± r)/2 and r is the Bloch vector length.
When any decoherence parameter is non-zero, the simulation switches from pure unitary evolution to the quantum trajectory (Monte Carlo wavefunction) method. Three Lindblad channels act during the inter-pulse free evolution periods:
| Channel | Collapse operator | Rate |
|---|---|---|
| Spin decay | L₁ = √γ₁ · σ₋ ⊗ I | γ₁ = 1/T₁ |
| Pure dephasing | L₂ = √(γφ/2) · σz ⊗ I | γφ = 1/T₂ − 1/(2T₁) |
| Motional heating | L₃ = √γh · I ⊗ a† | γh = (dn/dt) / 1000 μs⁻¹ |
For each inter-pulse interval Tm, the jump probabilities are:
A random number selects which (if any) jump occurs. If no jump: the non-Hermitian correction (I − Σ L†L/2 · Tm)|ψ⟩ is applied and renormalised. Results are averaged over Ntraj trajectories.
The pure dephasing rate γφ = 1/T₂ − 1/(2T₁) is clamped to ≥ 0. If T₂ > 2T₁, the pure dephasing contribution is zero — all decoherence comes from population decay.
When Nrep > 0, each spin expectation value is sampled as a Bernoulli experiment:
The ±1σ error bar is 2√(p(1−p)/Nrep). This models the quantum projection noise of a real single-ion experiment.
On browsers with WebGPU support (Chrome 113+, Edge 113+), the complex matrix multiplication — the bottleneck of the Taylor-series matrix exponential — is offloaded to a GPU compute shader. The WGSL kernel dispatches 8×8 workgroups for tiled complex GEMM.
The GPU path converts Float64 → Float32 for the shader and back. This provides ~7 significant digits, which is adequate for η ≈ 0.4 physics (where the coupling matrix elements span ~3 orders of magnitude). The CPU fallback uses Float64 throughout. The provenance manifest records which engine was used.
Every simulation result carries a SHA-256 hash computed over a JSON payload containing:
| Field | Content |
|---|---|
| code_version | Semantic version of the simulation engine |
| repo | GitHub repository URL |
| params | All input parameters (including decoherence and noise settings) |
| timestamp | ISO 8601 timestamp of the run |
| data_fingerprint | Array sums for σx,y,z, entropy, coherence; point count; max Fock leakage |
Manifest schema v2.0 supports three modes: single_run, sweep_1d, state_comparison. See ARCHITECTURE.md and schemas/manifest_v2.schema.json.
Structure: { manifest: {...}, data: {...} } (v1.2) or { schema_version: "2.0", payload: {...} } (v2.0)
The downloadable stroboscopic_sweep.py is a self-contained Python/scipy engine that reproduces the browser simulation with identical physics. It supports three modes: single run, 1D parameter sweep, and state comparison. See Getting Started for usage.