Density Matrices and Mixed States in Qiskit

If you have only worked with pure quantum states (state vectors like |ψ⟩ = α|0⟩ + β|1⟩), you have been working with an idealisation. Real quantum hardware is noisy. Qubits decohere. When a qubit interacts with its environment, the right description is no longer a single state vector: it is a density matrix.

Understanding density matrices lets you reason correctly about:

• What noise actually does to a quantum state
• Why measuring one qubit in an entangled pair leaves the other in a mixed state
• How error mitigation works mathematically
• What the output of DensityMatrixSimulator means in practice

Setup

pip install qiskit qiskit-aer

import numpy as np
from qiskit.quantum_info import DensityMatrix, partial_trace, state_fidelity, Statevector
from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel, depolarizing_error

Pure states vs mixed states

A pure state is one where you have complete knowledge of the quantum state. It is described by a state vector |ψ⟩. The density matrix of a pure state is:

ρ = |ψ⟩⟨ψ|

For the |0⟩ state:

# Pure state |0>
psi = Statevector([1, 0])  # |0>
rho_pure = DensityMatrix(psi)
print(rho_pure.data)
# [[1.+0.j, 0.+0.j],
#  [0.+0.j, 0.+0.j]]

For the |+⟩ = (|0⟩ + |1⟩)/√2 state:

psi_plus = Statevector([1/np.sqrt(2), 1/np.sqrt(2)])
rho_plus = DensityMatrix(psi_plus)
print(rho_plus.data)
# [[0.5+0.j, 0.5+0.j],
#  [0.5+0.j, 0.5+0.j]]

A mixed state represents a statistical mixture: you know the qubit is in one of several states but not which one. The density matrix is a probability-weighted average:

# Mixed state: 50% |0>, 50% |1>
# This is NOT the same as |+> = (|0> + |1>)/sqrt(2)
rho_0 = DensityMatrix([1, 0])  # |0>
rho_1 = DensityMatrix([0, 1])  # |1>

# Classical mixture: no coherence
rho_mixed = 0.5 * rho_0 + 0.5 * rho_1
print(rho_mixed.data)
# [[0.5+0.j, 0.+0.j],
#  [0.+0.j, 0.5+0.j]]

The key difference: the off-diagonal elements (coherences). The |+⟩ state has coherences of 0.5; the 50/50 classical mixture has coherences of 0. Decoherence is the process of these off-diagonal terms decaying to zero.

The purity check

A pure state satisfies Tr(ρ²) = 1. A mixed state satisfies Tr(ρ²) < 1. This is a quick diagnostic:

def purity(rho: DensityMatrix) -> float:
    return np.real(np.trace(rho.data @ rho.data))

print(f"Purity of |0>:      {purity(rho_pure):.4f}")    # 1.0
print(f"Purity of |+>:      {purity(rho_plus):.4f}")    # 1.0
print(f"Purity of mixture:  {purity(rho_mixed):.4f}")   # 0.5

# Maximum mixed state (fully decohered): identity / d
rho_max_mixed = DensityMatrix(np.eye(2) / 2)
print(f"Purity of maximally mixed: {purity(rho_max_mixed):.4f}")  # 0.5

Noise creates mixed states

Here is what depolarizing noise does in density matrix terms. Depolarizing noise with error rate p maps:

ρ → (1 - p) ρ + p * I/2

It mixes the pure state with the maximally mixed state, exactly what we saw above.

# Start with |+>
rho = DensityMatrix(Statevector([1/np.sqrt(2), 1/np.sqrt(2)]))
print("Before noise:")
print(np.round(rho.data, 4))

# Apply depolarizing noise manually
def depolarize(rho_data: np.ndarray, p: float) -> np.ndarray:
    """Apply single-qubit depolarizing noise with rate p."""
    I = np.eye(2, dtype=complex)
    return (1 - p) * rho_data + p * (I / 2)

p = 0.1
rho_noisy = depolarize(rho.data, p)
print(f"\nAfter depolarizing noise (p={p}):")
print(np.round(rho_noisy, 4))
print(f"Purity: {np.real(np.trace(rho_noisy @ rho_noisy)):.4f}")

The off-diagonal elements shrink by (1-p). With enough noise, ρ converges to I/2, the maximally mixed state where you know nothing about the qubit.

Simulating with density matrices in Qiskit

Qiskit’s AerSimulator can track density matrices through noisy circuits:

def run_density_matrix(circuit: QuantumCircuit, noise_model=None) -> np.ndarray:
    """Run a circuit and return the final density matrix."""
    qc = circuit.copy()
    qc.save_density_matrix()

    backend = AerSimulator(method='density_matrix', noise_model=noise_model)
    from qiskit import transpile
    result = backend.run(transpile(qc, backend), shots=1).result()
    return np.array(result.data()['density_matrix'])

# Bell state -- pure
qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)

rho_bell = run_density_matrix(qc)
print("Bell state density matrix (4x4):")
print(np.round(np.real(rho_bell), 4))
# Should have 0.5 in the [0,0], [0,3], [3,0], [3,3] positions

purity_bell = np.real(np.trace(rho_bell @ rho_bell))
print(f"Purity: {purity_bell:.4f}")  # Should be 1.0

# Same circuit with noise
noise = NoiseModel()
noise.add_all_qubit_quantum_error(depolarizing_error(0.05, 2), ['cx'])
rho_noisy = run_density_matrix(qc, noise_model=noise)
purity_noisy = np.real(np.trace(rho_noisy @ rho_noisy))
print(f"\nPurity with 5% CX error: {purity_noisy:.4f}")
# Should be < 1.0

The partial trace

When you have a multi-qubit system and want to describe only one subsystem, you use the partial trace. This traces out (averages over) the other qubits.

For the Bell state (|00⟩ + |11⟩)/√2, tracing out one qubit gives the maximally mixed state I/2 on the other, even though the joint state is pure. This is the quantum information signature of entanglement: subsystems of entangled pure states are mixed.

from qiskit.quantum_info import partial_trace

# Bell state
qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)
rho_bell = DensityMatrix(qc)

# Trace out qubit 1 to get the reduced state of qubit 0
rho_0_reduced = partial_trace(rho_bell, [1])  # trace out qubit index 1
print("Reduced state of qubit 0 (Bell state, trace out qubit 1):")
print(np.round(rho_0_reduced.data, 4))
# Should be [[0.5, 0], [0, 0.5]] -- maximally mixed

print(f"Purity of subsystem: {purity(rho_0_reduced):.4f}")  # 0.5

# Compare: product state |+>|0>
qc2 = QuantumCircuit(2)
qc2.h(0)
# qubit 1 stays in |0>
rho_product = DensityMatrix(qc2)
rho_0_reduced2 = partial_trace(rho_product, [1])
print("\nReduced state of qubit 0 (product state |+>|0>):")
print(np.round(rho_0_reduced2.data, 4))
print(f"Purity of subsystem: {purity(rho_0_reduced2):.4f}")  # 1.0 -- pure

If the global state is entangled, the subsystem is mixed. If the global state is a product state, each subsystem is pure. You can use this as a way to verify whether two qubits are entangled.

Expectation values from density matrices

The expectation value of an observable O in state ρ is:

⟨O⟩ = Tr(ρ O)

# Pauli Z matrix
Z = np.array([[1, 0], [0, -1]], dtype=complex)
I = np.eye(2, dtype=complex)

# Expectation of Z for |0>
rho_0 = np.array([[1, 0], [0, 0]], dtype=complex)
exp_z = np.real(np.trace(rho_0 @ Z))
print(f"<Z> for |0>: {exp_z:.4f}")   # 1.0

# Expectation of Z for |1>
rho_1 = np.array([[0, 0], [0, 1]], dtype=complex)
exp_z = np.real(np.trace(rho_1 @ Z))
print(f"<Z> for |1>: {exp_z:.4f}")   # -1.0

# Expectation of Z for maximally mixed state
rho_max = np.eye(2) / 2
exp_z = np.real(np.trace(rho_max @ Z))
print(f"<Z> for I/2: {exp_z:.4f}")   # 0.0 -- no information

# Depolarizing noise erases the signal proportionally
p = 0.3
rho_noisy = depolarize(rho_0, p)
exp_z = np.real(np.trace(rho_noisy @ Z))
print(f"<Z> for depolarized |0> (p=0.3): {exp_z:.4f}")  # (1-p) = 0.7

This is why depolarizing noise biases expectation values toward 0; it mixes in the I/2 state which contributes nothing to any Pauli expectation value. Zero-noise extrapolation exploits this linearity to recover the ideal value.

Fidelity: measuring how close two states are

State fidelity measures how similar two quantum states are. For pure states, F(ρ, σ) = |⟨ψ|φ⟩|². For mixed states it generalises:

from qiskit.quantum_info import state_fidelity

# Ideal Bell state
qc_ideal = QuantumCircuit(2)
qc_ideal.h(0)
qc_ideal.cx(0, 1)
rho_ideal = DensityMatrix(qc_ideal)

# Noisy Bell state (5% depolarizing on CX)
rho_noisy_dm = DensityMatrix(np.array(run_density_matrix(qc, noise_model=noise)))

fid = state_fidelity(rho_ideal, rho_noisy_dm)
print(f"Bell state fidelity under 5% CX depolarizing: {fid:.4f}")
# Typical value: 0.93 to 0.97

Fidelity of 1.0 means identical states. Fidelity of 0.5 for a single qubit means you are halfway to the maximally mixed state. Hardware benchmarks often report “process fidelity” or “state fidelity” of their prepared Bell states; now you know what that number means.

Putting it together: tracking decoherence

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Watch a qubit decohere as we add more noise
error_rates = np.linspace(0, 0.5, 20)
purities = []
coherences = []

rho_plus = np.array([[0.5, 0.5], [0.5, 0.5]], dtype=complex)  # |+>

for p in error_rates:
    rho_noisy = depolarize(rho_plus, p)
    purities.append(np.real(np.trace(rho_noisy @ rho_noisy)))
    coherences.append(np.real(rho_noisy[0, 1]))

# At p=0: purity=1, coherence=0.5
# At p=0.5: purity=0.625, coherence=0.25 (partially decohered)
print(f"Purity at p=0:   {purities[0]:.4f}")   # 1.0000
print(f"Purity at p=0.5: {purities[-1]:.4f}")  # 0.6250
print(f"Coherence at p=0:   {coherences[0]:.4f}")   # 0.5000
print(f"Coherence at p=0.5: {coherences[-1]:.4f}")  # 0.2500

When to use density matrices

Use density matrix formalism when:

• Modelling noise: noise channels (depolarizing, amplitude damping, phase damping) are naturally defined as density matrix transformations
• Subsystems of entangled states: use partial_trace to describe one qubit when it is part of a larger entangled system
• Verifying error mitigation: compare the density matrix before and after mitigation to quantify improvement
• Calculating fidelity: state fidelity is defined for density matrices and extends cleanly to mixed states

Stick to state vectors (Statevector) when:

• The circuit is noiseless (simulation only)
• You only need amplitudes and measurement probabilities for a small number of qubits
• Performance matters: density matrix simulation scales as 4ⁿ vs 2ⁿ for state vectors

The density matrix formulation is not more “correct” than state vectors for pure, noiseless systems; it is just more powerful, and more expensive. Use it when the physics demands it.