Concepts Intermediate Free 18/60 in series 25 min read By · Quantum Zeitgeist

Density Matrices and Mixed States in Qiskit

Understand density matrices, mixed states, and the partial trace. Essential for reasoning about noise, decoherence, and real hardware, with working Qiskit examples.

Prerequisites

  • Python proficiency
  • Beginner quantum computing concepts (superposition, entanglement)
  • Linear algebra basics
In this guide10 sections
  1. 01Setup
  2. 02Pure states vs mixed states
  3. 03The purity check
  4. 04Noise creates mixed states
  5. 05Simulating with density matrices in Qiskit
  6. 06The partial trace
  7. 07Expectation values from density matrices
  8. 08Fidelity: measuring how close two states are
  9. 09Putting it together: tracking decoherence
  10. 10When to use density matrices

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.

Was this tutorial helpful?

Continue learning with these guides