Density Matrix Simulation

Density-matrix circuit simulation propagates the ensemble state exactly within a finite noise model, up to floating-point error. It removes trajectory sampling variance from expectation values at the cost of $4^n$ complex entries for $n$ qubits. It is particularly useful as a small-system reference for testing approximate simulators.

Helpful background. Density-matrix storage supplies representation checks; quantum channels defines Kraus maps.

Local gates on a density tensor

With $n$ binary ket indices and $n$ bra indices, store

$$\rho_{i_0\cdots i_{n-1};j_0\cdots j_{n-1}}.$$

A gate $U$ on a chosen $k$-qubit target set contracts $U$ into the ket axes and $U^*$ into the corresponding bra axes. For one target,

$$\rho'_{\ldots i\ldots;\ldots j\ldots} =\sum_{ab}U_{ia}\rho_{\ldots a\ldots;\ldots b\ldots}U^*_{jb}.$$

Apply the two contractions sequentially. This avoids constructing a $2^n$ square full-register gate or a $4^n$ square superoperator. For fixed local gate size the work is proportional to the density-array size, with constants set by target dimension.

For several targets, keep their ordered tuple explicit. Permute ket axes so that those targets precede the untouched axes, and apply the matching permutation to bra axes. Reshape to $\rho_{ar;bs}$ with target dimension $d_t=2^k$ and remainder dimension $d_r=2^{n-k}$. Then

$$\rho'_{ir;js}=\sum_{ab}U_{ia}\rho_{ar;bs}U^*_{jb}.$$

Undo both permutations after the contractions. A gate whose local basis is ordered by targets $(2,0)$ generally changes meaning if that tuple is sorted to $(0,2)$ without also transforming the gate. Test an asymmetric complex gate on nonadjacent targets, since real gates or symmetric states can hide conjugation and ordering errors. The two dense contractions cost $O(d_t\,4^n)$ operations for a fixed gate; applying $r_K$ Kraus operators costs $O(r_Kd_t\,4^n)$.

A channel sums the same contractions over its local Kraus operators. Use separate input and accumulation buffers: overwriting the input after the first Kraus term changes later terms and can still produce a normalized-looking result.

Bell pair through amplitude damping

Prepare $|\Phi^+\rangle=(|00\rangle+|11\rangle)/\sqrt2$ and damp the second qubit with probability $p$. The output is

$$\rho'=\frac12 \begin{pmatrix} 1&0&0&\sqrt{1-p}\\ 0&0&0&0\\ 0&0&p&0\\ \sqrt{1-p}&0&0&1-p \end{pmatrix}.$$

In ordering $(00,01,10,11)$ the decay event sends $|11\rangle$ to $|10\rangle$. This distinguishes which qubit was targeted; damping the first qubit puts the population in $|01\rangle$ instead.

The exact measurement probabilities are $(1/2,0,p/2,(1-p)/2)$, and Bell-state fidelity is $(2-p+2\sqrt{1-p})/4$. At $p=1$ the result is a mixture of $|00\rangle$ and $|10\rangle$, not a pure reset of the entire register.

Measurements and classical branches

For outcome $m$ with measurement operator $M_m$, the unnormalized branch is $\tilde\rho_m=M_m\rho M_m^\dagger$, with probability $p_m=\operatorname{Tr}\tilde\rho_m$. If the outcome is ignored, sum the branches. If later gates depend on $m$, keep the classical branch label until branches can validly be merged.

A density matrix with all branches averaged together is insufficient for arbitrary outcome-dependent control. Two branches can have the same quantum state yet different future classical instructions. Branch weights should sum to one; never normalize a branch with zero probability.

One exact representation retains unnormalized blocks $\rho_{\rm CQ}=\sum_m|m\rangle\langle m|\otimes\tilde\rho_m$. If the next operation is the unitary $U_m$ selected by the classical outcome, the final quantum state is

$$\rho_{\rm out}=\sum_m U_m\tilde\rho_m U_m^\dagger.$$

For example, measure $|+\rangle$ in the $Z$ basis, then apply $X$ only after outcome 1. The two blocks are $|0\rangle\langle0|/2$ and $|1\rangle\langle1|/2$, and the controlled continuation produces $|0\rangle\langle0|$. Averaging first gives $I/2$; no outcome-independent unitary acting on that average can reproduce the result. Store branch traces as weights, or store normalized states with separate weights, but do not count the probability twice.

Sampling outcomes introduces shot noise even when state propagation is deterministic. Exact probabilities and finite-shot histograms are different outputs; use measurement sampling for the latter.

Resource and correctness checks

Complex double storage is $16\,4^n$ bytes: $n=15$ already needs 16 GiB for one state array. Kraus buffers and contraction copies can multiply peak memory. Plan memory from actual workspace requirements.

Test local contractions against full Kronecker matrices for two or three qubits, including nonadjacent targets and complex gates. Check each channel’s Choi matrix separately and state diagnostics after circuit segments. Compare a noiseless run to the projector from a statevector simulator.

Do not truncate tiny off-diagonal entries without an observable-error study. Their collective contribution can matter, and entrywise sparsification is not a positivity-preserving operation.

Exercise

Why can all terminal computational-basis shots be sampled from one evolved density matrix, but mid-circuit adaptive measurements may require repeated branching?

Solution

Terminal shots are independent categorical draws from the unchanged diagonal distribution. A mid-circuit outcome changes the conditional state and can select subsequent gates, so each branch has a different continuation. Exact branch enumeration or correctly sampled conditional trajectories is required.

References

• QuTiP developers, “Basic Operations on Quantum Objects”, QuTiP documentation — tensor dimensions and operator representations.
• Watrous, J., The Theory of Quantum Information, Cambridge University Press (2018), author’s text — channels, instruments, and classical outcomes.