QuantumPulse.jl

QuantumPulse.jl is a high-performance computational physics library in Julia for the simulation of open quantum systems, density matrix dynamics, and dissipative quantum state evolution governed by the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) master equation.

Designed for quantum optics, cavity quantum electrodynamics (QED), and qubit decoherence modeling, QuantumPulse.jl implements both exact analytical matrix-exponential propagators on vectorized Liouville space and adaptive high-order numerical ODE integrators, alongside exact steady-state solvers and quantum information metrics.

---

Key Features

• Exact Liouvillian Superoperator Vectorization: Transforms the Lindblad master equation into an isomorphic linear dynamical system $\frac{d}{dt}|\rho\rangle!\rangle = \mathcal{K} |\rho\rangle!\rangle$ on $\mathbb{C}^{N^2 \times N^2}$, enabling exact analytical matrix exponential evolution and spectral decomposition.
• Multiple Time-Evolution Engines:
  • mesolve_analytic: Analytical matrix exponential propagator $\exp(\mathcal{K} t)$ with eigensystem caching ($O(N^4)$ per time step).
  • mesolve_rk4: Classical 4th-order Runge-Kutta integrator for arbitrary time-dependent driving Hamiltonians $H(t)$.
  • mesolve_adaptive: Embedded 5(4) Dormand-Prince (DP54) adaptive step-size integrator with continuous local truncation error control.
• Exact Null-Space Steady-State Solver: Directly determines the stationary density matrix $\rho_{\mathrm{ss}}$ satisfying $\mathcal{K}|\rho_{\mathrm{ss}}\rangle!\rangle = 0$ and $\mathrm{Tr}(\rho_{\mathrm{ss}}) = 1$ via singular value decomposition (SVD) and QR factorizations, achieving machine-precision agreement with analytical optical Bloch solutions ($< 10^{-15}$ trace distance).
• Comprehensive Quantum Information Diagnostics: State purity $\gamma = \mathrm{Tr}(\rho^2)$, von Neumann entropy $S(\rho) = -\mathrm{Tr}(\rho \ln \rho)$, Uhlmann quantum fidelity $F(\rho, \sigma)$, trace distance $D(\rho, \sigma)$, expectation values $\langle A \rangle = \mathrm{Tr}(\rho A)$, and geometric Bloch sphere vector mapping $(u, v, w)$ with radius $r \le 1$.
• Zero External Dependencies: Built purely on Julia's standard library (LinearAlgebra, Printf, Random, Test), guaranteeing minimal compile latency and maximum cross-platform reproducibility.

---

Architecture Overview

graph TD
    subgraph HilbertSpace ["Hilbert Space Operators"]
        H["Hamiltonian H(t)"]
        C["Collapse Ops {L_k, γ_k}"]
        State["State |ψ⟩ / ρ0"]
    end

    subgraph Superoperators ["Liouville Space Superoperators"]
        Unitary["-i (I ⊗ H - H^T ⊗ I)"]
        Dissipator["∑_k γ_k [ L_k* ⊗ L_k - 0.5 (I ⊗ L_k†L_k + (L_k†L_k)^T ⊗ I) ]"]
        K["Liouvillian Generator Matrix K ∈ ℂ^(N² × N²)"]
    end

    subgraph Solvers ["Evolution Engines"]
        Exp["mesolve_analytic: exp(K t)"]
        RK4["mesolve_rk4: 4th-Order Runge-Kutta"]
        DP54["mesolve_adaptive: Dormand-Prince 5(4)"]
        SS["steadystate: Augmented SVD / QR Kernel"]
    end

    subgraph Observables ["Information-Theoretic Diagnostics"]
        Purity["Purity Tr(ρ²)"]
        Entropy["von Neumann Entropy S(ρ)"]
        Fid["Uhlmann Fidelity F(ρ, σ)"]
        Dist["Trace Distance D(ρ, σ)"]
        Bloch["Bloch Sphere Trajectory (u, v, w)"]
    end

    H --> Unitary
    C --> Dissipator
    Unitary --> K
    Dissipator --> K
    State --> Exp
    State --> RK4
    State --> DP54
    K --> Exp
    K --> SS
    Exp --> Observables
    RK4 --> Observables
    DP54 --> Observables
    SS --> Observables

Loading

---

Mathematical Formulation

1. The Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) Equation

For an open quantum system coupled to a Markovian thermal reservoir, the reduced density operator $\rho(t) \in \mathbb{C}^{N \times N}$ evolves according to the GKSL master equation ($\hbar = 1$):

$$\frac{d\rho(t)}{dt} = -i [H(t), \rho(t)] + \sum_{k} \gamma_k \mathcal{D}[L_k]\rho(t)$$

where the Lindblad dissipator $\mathcal{D}[L_k]\rho$ is defined as:

$$\mathcal{D}[L_k]\rho = L_k \rho L_k^\dagger - \frac{1}{2} \left{ L_k^\dagger L_k, \rho \right} = L_k \rho L_k^\dagger - \frac{1}{2} L_k^\dagger L_k \rho - \frac{1}{2} \rho L_k^\dagger L_k$$

Here:

• $H(t)$ is the system Hamiltonian.
• $L_k$ are system collapse/jump operators describing environmental interaction channels (e.g. spontaneous emission $\sigma_-$, pure dephasing $\sigma_z$, cavity photon annihilation $a$).
• $\gamma_k \ge 0$ are the corresponding dissipation rates.
• $\mathcal{D}[L_k]$ is completely positive and trace-preserving (CPTP), ensuring $\mathrm{Tr}(\rho(t)) = 1$ and $\rho(t) \ge 0$ for all $t \ge 0$.

2. Liouville Space Vectorization & Matrix Superoperators

Using the column-major vectorization isomorphism $\mathrm{vec}(A X B) = (B^T \otimes A) \mathrm{vec}(X)$, the master equation is recast into a linear ordinary differential equation in $N^2$-dimensional Liouville space:

$$\frac{d}{dt} |\rho(t)\rangle!\rangle = \mathcal{K} |\rho(t)\rangle!\rangle$$

The Liouvillian generator matrix $\mathcal{K} \in \mathbb{C}^{N^2 \times N^2}$ is given by:

$$\mathcal{K} = -i \left( I_N \otimes H - H^T \otimes I_N \right) + \sum_k \gamma_k \left[ L_k^* \otimes L_k - \frac{1}{2} I_N \otimes (L_k^\dagger L_k) - \frac{1}{2} (L_k^\dagger L_k)^T \otimes I_N \right]$$

For time-independent $H$, the formal solution is computed via the matrix exponential:

$$|\rho(t)\rangle!\rangle = \exp(\mathcal{K} t) |\rho(0)\rangle!\rangle$$

When $\mathcal{K}$ is non-defective ($\mathcal{K} = V \Lambda V^{-1}$), arbitrary time evaluation simplifies to:

$$|\rho(t)\rangle!\rangle = V \exp(\Lambda t) V^{-1} |\rho(0)\rangle!\rangle$$

3. Exact Stationary State Solution

The steady-state density operator $\rho_{\mathrm{ss}}$ satisfies $\frac{d\rho}{dt} = 0$, equivalent to finding the null-space of $\mathcal{K}$ subject to $\mathrm{Tr}(\rho_{\mathrm{ss}}) = 1$:

$$\begin{pmatrix} \mathcal{K} \ \mathbf{t}^T \end{pmatrix} |\rho_{\mathrm{ss}}\rangle!\rangle = \begin{pmatrix} \mathbf{0} \ 1 \end{pmatrix}$$

where $\mathbf{t} \in \mathbb{R}^{N^2}$ has $t_{(i-1)N + i} = 1$ for $i = 1, \dots, N$ and zeros elsewhere.

4. Quantum Information Diagnostics

Observable	Mathematical Definition	Physical Interpretation
Expectation Value	$\langle A \rangle = \mathrm{Tr}(\rho A)$	Macroscopic measurement expectation
Purity	$\gamma = \mathrm{Tr}(\rho^2) \in [1/N, 1]$	$\gamma = 1 \implies$ pure state; $\gamma = 1/N \implies$ maximally mixed
von Neumann Entropy	$S(\rho) = -\sum_i \lambda_i \ln \lambda_i$	Information entropy in natural units (nats)
Uhlmann Fidelity	$F(\rho, \sigma) = \left( \mathrm{Tr}\sqrt{\sqrt{\rho}\sigma\sqrt{\rho}} \right)^2$	State overlap metric ($F \in [0, 1]$)
Trace Distance	$D(\rho, \sigma) = \frac{1}{2} \mathrm{Tr}\sqrt{(\rho-\sigma)^\dagger(\rho-\sigma)}$	Distinguishability metric ($D \in [0, 1]$)
Bloch Vector (Qubit)	$\vec{r} = (\langle\sigma_x\rangle, \langle\sigma_y\rangle, \langle\sigma_z\rangle)$	Coordinates on/inside the unit Bloch sphere ($r \le 1$)

---

Performance Benchmarks

Benchmarks executed on AMD x86_64 CPU under WSL Ubuntu (Julia 1.13.3):

=================================================================
  QuantumPulse.jl — High-Performance Lindblad Solver Benchmark   
=================================================================

[Benchmark 1] Two-Level Damped Rabi System (N = 2, Liouvillian N² = 4)
  Analytic (matrix exp):     2.443 ms (201 steps)
  Numerical RK4:             5.653 ms (201 steps)
  Max trajectory error:   1.47e-06 (Trace distance)

[Benchmark 2] Jaynes-Cummings Cavity QED (N = 20, Liouvillian N² = 400)
  Analytic (matrix exp):   272.588 ms (101 steps)
  Final state purity:       0.5014
  Final state entropy:      0.6917 nats

[Benchmark 3] Steady-State Solver Accuracy & Convergence
  Steady-State Trace Distance: 4.06e-16 (vs Analytical Optical Bloch)
  Steady-State Quantum Fidelity: 1.0000000000

---

Installation & Environment Setup

Prerequisites

• Julia 1.9 or newer (tested on Julia 1.13.3)

Installation

Clone the repository and instantiate the package environment:

git clone https://github.com/ZengLLQ/quantum-pulse.git
cd quantum-pulse
julia --project=. -e "using Pkg; Pkg.instantiate()"

To run the full automated test suite:

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

---

Quick Start & Code Examples

1. Damped Rabi Oscillations of a Qubit

using QuantumPulse

# Define a driven two-level system with relaxation (γ1) and dephasing (γ_phi)
model = TwoLevelRabi(Δ = 0.0, Ω = 1.0, γ1 = 0.1, γ_phi = 0.05)
H = hamiltonian(model)
c_ops = collapse_operators(model)
r = rates(model)

# Initial pure state: Excited state |e⟩
rho0 = basis_state(2, 0)
tlist = collect(range(0.0, stop=25.0, length=51))

# Track Pauli observables
e_ops = Dict("sx" => sigma_x(), "sy" => sigma_y(), "sz" => sigma_z())

# Evolve system analytically via matrix exponential
res = mesolve_analytic(H, rho0, tlist; c_ops=c_ops, rates=r, e_ops=e_ops)

# Print final purity and Bloch coordinates
bvec, r_bloch = bloch_vector(res.states[end])
println("Final Bloch vector: ", bvec)
println("Final Purity: ", purity(res.states[end]))

2. Cavity QED Vacuum Rabi Oscillations (Jaynes-Cummings)

using QuantumPulse
using LinearAlgebra

dim_cavity = 6
jc = JaynesCummings(ω_a = 1.0, ω_c = 1.0, g = 0.15, γ_a = 0.01, κ_c = 0.02, dim_cavity = dim_cavity)

H = hamiltonian(jc)
c_ops = collapse_operators(jc)
r = rates(jc)

# Initial state: Excited atom |e⟩ ⊗ Cavity vacuum |0⟩
atom_e = basis_state(2, 0)
cav_0 = basis_state(dim_cavity, 0)
rho0 = QuantumState(Matrix{ComplexF64}(kron(atom_e.data, cav_0.data)))

tlist = collect(range(0.0, stop=30.0, length=61))
e_ops = Dict(
    "sz" => kron(sigma_z(), identity_op(dim_cavity)),
    "n"  => kron(identity_op(2), number_op(dim_cavity))
)

result = mesolve_analytic(H, rho0, tlist; c_ops=c_ops, rates=r, e_ops=e_ops)

---

CLI Commands

The package includes a built-in command-line interface for running benchmarks and simulations:

# Run performance and accuracy benchmarks
julia --project=. src/cli.jl benchmark

# Run Rabi oscillation trajectory demo with ASCII table
julia --project=. src/cli.jl rabi

# Execute standalone example scripts
julia --project=. examples/rabi_decay.jl
julia --project=. examples/cavity_decay.jl
julia --project=. examples/jaynes_cummings_dynamics.jl

---

Design Trade-Offs & Known Limitations

1. Dense vs Sparse Liouvillian Scaling: The vectorized Liouvillian superoperator matrix $\mathcal{K}$ scales as $N^2 \times N^2$ for an $N$-dimensional Hilbert space. While exact matrix exponential propagation is optimal for $N \le 40$ (e.g. qubit arrays, cavity QED with moderate Fock truncation), macroscopic Hilbert spaces ($N > 100$) require sparse Krylov subspace or quantum trajectory Monte Carlo wave-function (MCWF) methods.
2. Markovian Approximation: The GKSL master equation assumes weak system-bath coupling and Markovian memoryless reservoirs. Non-Markovian structured environments with spectral memory require hierarchical equations of motion (HEOM) or pseudo-mode expansions.
3. Truncated Fock Spaces: Infinite-dimensional bosonic operators ($a, a^\dagger$) are represented in finite truncated Fock bases. The highest Fock state $|N_{\mathrm{fock}}-1\rangle$ experiences artificial truncation effects if populations reach the boundary.

---

License

This project is licensed under the MIT License.