Density Matrices in Code

A density matrix is a complex array together with a basis convention and physical constraints. Storing a square matrix is easy; preserving its interpretation as a state requires checking Hermiticity, unit trace, and positivity separately. This page develops those checks and subsystem index contractions. Circuit evolution and master-equation integration use this representation.

Helpful background. Reduced density operators explains the physical meaning of a mixed state; tensor products in code fixes array ordering.

Density operators as complex arrays

Use an orthonormal basis and store $\rho_{ij}=\langle i|\rho|j\rangle$. For an ensemble of normalized columns $\psi_a$ with nonnegative weights $p_a$ summing to one,

$$\rho=\sum_a p_a\psi_a\psi_a^\dagger.$$

The conjugation is essential. An outer product without complex conjugation need not be Hermitian. For any vector $v$, $v^\dagger\rho v=\sum_a p_a|\psi_a^\dagger v|^2\geq0$, which establishes positivity without diagonalization. A low-rank factor $X$ with columns $\sqrt{p_a}\psi_a$ gives $\rho=XX^\dagger$ and can reduce storage when rank is small.

For a dense $d$-dimensional state, complex double precision costs approximately $16d^2$ bytes before work arrays. Hermiticity does not make a generic matrix sparse. Sparse storage helps only while the physical representation has few nonzero entries.

Diagnostics with scale-aware tolerances

Compute the anti-Hermitian residual before calling a Hermitian eigensolver:

$$r_H=\frac{\|\rho-\rho^\dagger\|_F}{\max(1,\|\rho\|_F)}, \qquad r_T=|\operatorname{Tr}\rho-1|.$$

Then diagonalize the Hermitian part $R=(\rho+\rho^\dagger)/2$ for diagnostic eigenvalues. Record $\lambda_{\min}(R)$ and the negative weight $\nu=\sum_{\lambda_j<0}|\lambda_j|$. Small roundoff-scale negative eigenvalues are different from a negative weight that grows with time or persists under refinement.

For a physical state, $1/d\leq\operatorname{Tr}\rho^2\leq1$. Purity within these bounds does not prove positivity: the eigenvalues $(0.6,0.5,-0.1)$ sum to one and have purity $0.62$. Likewise, trace normalization cannot repair a wrong channel.

A useful tolerance combines machine precision, matrix size, the scale of the computation, and the measured solver error. Do not impose one absolute tolerance on all state sizes. If projecting negative eigenvalues to zero is part of an estimator, report the projection and its magnitude; do not hide an unstable evolution behind it.

The negative weight gives a quantitative repair diagnostic. For a Hermitian, trace-one matrix $R$, write $R=R_+-R_-$ with positive parts on orthogonal eigenspaces and $\nu=\operatorname{Tr}R_-$. Clipping and renormalizing gives $\sigma=R_+/(1+\nu)$, because $\operatorname{Tr}R_+=1+\nu$. Consequently

$$\frac12\|\sigma-R\|_1=\nu.$$

The positive and negative spectral sectors each contribute $\nu$ to the trace norm. For the example $(0.6,0.5,-0.1)$, the repair changes the matrix by $0.1$ in this measure; it cannot be described as negligible roundoff. This calculation quantifies one specified repair, not an error bound relative to the unknown physical state.

Partial trace as an index contraction

With composite ordering $i=a\,d_B+b$, reshape $\rho$ into $R_{a b a' b'}$. The reduced state is

$$(\rho_A)_{a a'}=\sum_b R_{a b a' b}.$$

For a pure state reshape $\psi$ into $C_{ab}$ and use $\rho_A=CC^\dagger$, avoiding the full $d_A d_B$ square matrix. Nonadjacent subsystems require an explicit axis permutation before reshaping. Tracing arbitrary axes of an array without documenting their labels is a common source of plausible wrong answers.

import numpy as np

def reduce_pure_state(psi, da, db):
    c = np.asarray(psi, dtype=complex).reshape(da, db)
    return c @ c.conj().T

psi = np.array([1, 0, 0, 1j], dtype=complex) / np.sqrt(2)
rho_a = reduce_pure_state(psi, 2, 2)
assert np.allclose(rho_a, np.eye(2) / 2)

The phase $i$ makes this a stronger conjugation test than an entirely real Bell state. Test also a product state: its reduction must have rank one. These tests exercise different failure modes.

An equal-dimension Bell pair does not test which subsystem was retained. Use a $2\times3$ fixture,

$$|\psi\rangle=\frac{|0,0\rangle+i|1,2\rangle}{\sqrt2}.$$

The two reductions must have different shapes: $\rho_A=I_2/2$ and $\rho_B=\operatorname{diag}(1/2,0,1/2)$. In the same coefficient convention, $\rho_B=C^{\mathsf T}C^*$; $C^\dagger C$ is its transpose and is generally different for complex states. Compare a general complex fixture against the explicit index sum as well as the simple analytic example.

Expectations and truncation

Evaluate $\langle A\rangle=\operatorname{Tr}(\rho A)$, equivalently $\sum_{ij}\rho_{ij}A_{ji}$. A Hermitian observable should have a real expectation to numerical accuracy. An unexpected imaginary part is a useful signal of mismatched conjugations or index order.

In a truncated oscillator basis, a positive unit-trace matrix is only a valid state within the truncated model. Monitor boundary occupation and increase the cutoff. These physical truncation errors are independent of matrix roundoff and time-step error.

Exercise

For $\rho=p|0\rangle\langle0|+(1-p)|+\rangle\langle+|$, derive the determinant and purity and identify the pure cases.

Solution

In the ordered basis $(|0\rangle,|1\rangle)$, the diagonal entries are $(1+p)/2$ and $(1-p)/2$, and the off-diagonal entries are $(1-p)/2$. Hence $\det\rho=p(1-p)/2$ and $\operatorname{Tr}\rho^2=1-p+p^2$. Positivity holds for $0\leq p\leq1$; the state is pure only at $p=0,1$. Sampling intermediate values tests mixtures of nonorthogonal states.

References

• Johansson, J. R., Nation, P. D., and Nori, F., “QuTiP: An open-source Python framework for the dynamics of open quantum systems,” Computer Physics Communications 183, 1760–1772 (2012), doi:10.1016/j.cpc.2012.02.021 — density-operator computation and solver context.
• QuTiP developers, “Lindblad Master Equation Solver”, QuTiP documentation — state representations and expectation-value calculations.
• Watrous, J., The Theory of Quantum Information, Cambridge University Press (2018), author’s text — reduced states and trace-norm diagnostics.