Open quantum system simulation

An atom in a lossy cavity does not evolve by the Schrödinger equation alone, because it keeps leaking photons into an environment nobody tracks. The Lindblad master equation is the standard tool for such open systems, evolving the density matrix under both coherent dynamics and dissipation. This project, the final project for Columbia’s graduate quantum optics course (a small class of Ph.D. students, taken in my first year), simulates N excited two-level atoms decaying into a shared cavity mode and into free space.

Background

The starting point is a single excited atom sitting at an anti-node of the cavity field, on resonance with the cavity, so no coherent Hamiltonian term survives and the dynamics is pure dissipation. The atom decays into the cavity at a rate $\Gamma_c$ and into free space at a rate $\Gamma_0$, and its master equation reads

$\dot{\rho} = \frac{\Gamma_c + \Gamma_0}{2}\left(2\sigma_-\rho\,\sigma_+ - \sigma_+\sigma_-\rho - \rho\,\sigma_+\sigma_-\right)$

where $\rho$ is the density matrix and $\sigma_\pm$ are the raising and lowering operators of the two-level atom. The excited-state population then obeys $\partial_t \rho_{ee} = -(\Gamma_c + \Gamma_0)\,\rho_{ee}$, plain exponential decay, which serves as the analytic benchmark. The simulation works in units of $\Gamma_0 t$ with $\Gamma_c = 5\Gamma_0$ unless otherwise specified. Direct integration with an ODE solver (scipy’s solve_ivp) reproduces this analytic spontaneous-emission decay exactly for one atom.

For N atoms sharing the cavity mode the master equation generalizes to

$\dot{\rho} = \frac{1}{2}\sum_{i,j=1}^{N}\left(\Gamma_c + \Gamma_0\,\delta_{ij}\right)\left(2\sigma_-^{j}\rho\,\sigma_+^{i} - \sigma_+^{i}\sigma_-^{j}\rho - \rho\,\sigma_+^{i}\sigma_-^{j}\right)$

where $\sigma_\pm^{i}$ acts on atom $i$ alone, built as a tensor product of identities with one Pauli operator inserted at slot $i$. The cross terms with $i \neq j$ come from the shared cavity and carry $\Gamma_c$; the Kronecker delta adds $\Gamma_0$ only on the diagonal, since each atom leaks into free space independently.

Quantum trajectories

The second numerical route unravels the master equation into stochastic wavefunction evolutions. The equation is rewritten as

$\dot{\rho} = -i\left(H_{\text{eff}}\,\rho - \rho\,H_{\text{eff}}^{\dagger}\right) + \sum_m c_m\,\rho\,c_m^{\dagger}$

with the effective non-Hermitian Hamiltonian

$H_{\text{eff}} = H - \frac{i}{2}\sum_m c_m^{\dagger} c_m$

where the $c_m$ are jump operators. Each trajectory propagates a state under $H_{\text{eff}}$, whose anti-Hermitian part shrinks the norm, and at each step the norm deficit is compared against a random number to decide whether a quantum jump fires. A single trajectory is nearly meaningless, a step function of one emission event, but averaging a few hundred trajectories (over 200 in practice) converges to the master-equation result while never storing more than a wavefunction.

For one atom the jump operator is simply $c = \sigma_-$. For N atoms the operators have to be chosen so that the unraveling reproduces the N-atom master equation above. A collective operator $c_1 = \sqrt{\Gamma_c}\,\sum_i \sigma_-^{i}$ generates all the cavity cross terms, and one operator $\sqrt{\Gamma_0}\,\sigma_-^{i}$ per atom supplies the independent free-space channel, giving $N+1$ jump operators in total. With several decay channels open, a second random number selects which operator fires at each jump. The multi-atom trajectories agree well with the master equation at short times, though they drift from it at intermediate times as N grows.

Scaling

The many-atom problem is where it gets expensive. The density matrix grows as $2^N \times 2^N$, and the simulation reaches seven atoms in about fifty seconds on a laptop; beyond seven atoms a run takes more than two minutes. Trajectories should win here in principle, since a wavefunction is exponentially smaller than a density matrix, yet the implementation stored the jump operators and intermediate quantities as full tensor-product arrays, so at seven atoms the trajectory method actually ran about ten seconds slower than direct integration. The writeup proposes the fix. Each basis state of N two-level atoms is a bit string, which maps one to one onto an integer, and the jump operators just flip bits, so the whole trajectory procedure could run on integers instead of exponentially large arrays.

Superradiance

The physics payoff sits in the emission rate, the negative time derivative of the total excited-state population $R = -\,d\langle \sum_i \sigma_{ee}^{i}\rangle/dt$. With free-space decay switched off ($\Gamma_0 = 0$, $\Gamma_c \neq 0$), the rate of energy loss is not monotonic; it rises to an early peak that grows steeply with atom number, which is the onset of the Dicke superradiant burst, the atoms synchronizing their emission through the shared cavity mode. Superradiance lives entirely in the interaction of the N excited atoms with a common light field, so the burst needs that shared mode to dominate. Letting free-space decay dominate instead ($\Gamma_0 \gg \Gamma_c$) makes the peak vanish, because the atoms then decay independently and no collective emission builds up.