Quantum State Tomography with Qiskit

Quantum state tomography (QST) answers a specific question: given a device that prepares a quantum state, what state does it actually prepare? The answer requires many measurements in many bases, classical post-processing to reconstruct a density matrix, and some care about when the answer is meaningful.

This tutorial covers the full workflow using Qiskit’s StateTomography class, explains the math behind reconstruction, and shows how to interpret fidelity numbers in practice. It also covers process tomography, gate set tomography, and how tomography compares with randomized benchmarking for characterizing quantum hardware.

Setup

pip install qiskit qiskit-aer qiskit-experiments matplotlib

Note: All code in this tutorial requires qiskit-experiments. Install it before running any blocks.

import numpy as np
from qiskit import QuantumCircuit
from qiskit.quantum_info import (
    DensityMatrix, state_fidelity, Statevector,
    SparsePauliOp, Operator, partial_trace
)
from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel, depolarizing_error
from qiskit_experiments.library import StateTomography
from qiskit_experiments.framework import ExperimentData

The mathematics of density matrices

A density matrix (also called a density operator) provides the most general description of a quantum state. While a pure state can be written as a ket |ψ⟩, a mixed state requires a density matrix. The density matrix for a statistical ensemble of pure states is

$$\rho = \sum_i p_i \, |\psi_i\rangle\langle\psi_i|$$

where each $p_i \geq 0$ and $\sum_i p_i = 1$. The probabilities $p_i$ represent classical uncertainty about which pure state was prepared. A pure state has $\rho = |\psi\rangle\langle\psi|$ with a single term.

Three properties define a valid density matrix:

1. Hermitian: $\rho = \rho^\dagger$. This ensures all expectation values are real.
2. Positive semidefinite: all eigenvalues of $\rho$ are non-negative. This ensures all probabilities computed from $\rho$ are non-negative.
3. Unit trace: $\text{Tr}(\rho) = 1$. This ensures probabilities sum to one.

Any 2x2 matrix satisfying these three properties represents a valid single-qubit state.

Bloch vector representation

For a single qubit, the density matrix can be written in terms of the Pauli matrices:

$$\rho = \frac{1}{2}\left(I + r_x X + r_y Y + r_z Z\right)$$

where $\mathbf{r} = (r_x, r_y, r_z)$ is the Bloch vector. Each component is the expectation value of the corresponding Pauli operator:

• $r_x = \text{Tr}(\rho \, X) = \langle X \rangle$
• $r_y = \text{Tr}(\rho \, Y) = \langle Y \rangle$
• $r_z = \text{Tr}(\rho \, Z) = \langle Z \rangle$

The purity of a state is

$$\text{Tr}(\rho^2) = \frac{1 + |\mathbf{r}|^2}{2}$$

where $|\mathbf{r}|^2 = r_x^2 + r_y^2 + r_z^2$. For a pure state $|\mathbf{r}| = 1$ (on the surface of the Bloch sphere), and purity equals 1. For the maximally mixed state $\mathbf{r} = (0, 0, 0)$, purity equals 1/2. All physical states satisfy $|\mathbf{r}| \leq 1$.

This representation makes it clear what tomography must do: determine the three components $r_x$, $r_y$, $r_z$ from measurement data. Each component corresponds to measuring in one Pauli basis.

Measurement basis construction

A quantum computer typically measures in the computational (Z) basis. To measure in other bases, you apply a unitary rotation before measurement so that the target basis maps onto the computational basis.

Single-qubit rotations

For the Z basis, no rotation is needed. The eigenstates |0⟩ and |1⟩ are already the computational basis states.

For the X basis, apply a Hadamard gate before measurement. The Hadamard maps:

• $|+\rangle = (|0\rangle + |1\rangle)/\sqrt{2} \mapsto |0\rangle$
• $|-\rangle = (|0\rangle - |1\rangle)/\sqrt{2} \mapsto |1\rangle$

So measuring after Hadamard gives the X-basis outcome statistics. From these statistics you compute $r_x = P(+) - P(-) = P(0) - P(1)$.

For the Y basis, apply $S^\dagger$ followed by $H$ before measurement. The $S^\dagger H$ combination maps:

• $|{+i}\rangle = (|0\rangle + i|1\rangle)/\sqrt{2} \mapsto |0\rangle$
• $|{-i}\rangle = (|0\rangle - i|1\rangle)/\sqrt{2} \mapsto |1\rangle$

# Circuits that rotate each Pauli eigenbasis into the computational basis
qc_z = QuantumCircuit(1, 1)
qc_z.measure(0, 0)  # Z basis: no rotation needed

qc_x = QuantumCircuit(1, 1)
qc_x.h(0)           # X basis: Hadamard before measurement
qc_x.measure(0, 0)

qc_y = QuantumCircuit(1, 1)
qc_y.sdg(0)         # Y basis: S† then Hadamard before measurement
qc_y.h(0)
qc_y.measure(0, 0)

Multi-qubit measurement settings

For n qubits, you need to measure all tensor products of single-qubit Pauli operators. Each qubit independently gets one of the four Pauli operators {I, X, Y, Z}, but the identity operator does not require a separate measurement (it is always 1). The total number of independent measurement settings is $4^n - 1$.

For 2 qubits, this gives 15 independent settings. The possible Pauli products are: IX, IY, IZ, XI, XX, XY, XZ, YI, YX, YY, YZ, ZI, ZX, ZY, ZZ. The II component is always 1 (trace normalization) and does not need measurement.

For 3 qubits, this gives 63 independent settings. The Pauli products form a set of 64 three-fold tensor products (4 choices per qubit), minus the III identity. In practice, many of these can be measured simultaneously: for example, all 8 products involving only Z and I on each qubit (ZII, IZI, IIZ, ZZI, ZIZ, IZZ, ZZZ) are all diagonal in the computational basis and can be extracted from a single Z-basis measurement.

The StateTomography class in qiskit-experiments handles this grouping automatically, reducing the number of circuits below the raw $4^n - 1$ bound.

What tomography actually does

For a single qubit, any quantum state is fully described by a density matrix, a 2x2 Hermitian matrix with unit trace. It has 3 real degrees of freedom (the Bloch vector coordinates). To measure all three, you need expectation values in 3 bases: X, Y, and Z. For each basis you run many shots and compute the expectation value from the measurement statistics.

For n qubits the density matrix is 2ⁿ x 2ⁿ and requires 4ⁿ - 1 independent measurement settings. Two qubits: 15 settings. Three qubits: 63 settings. This exponential scaling is why tomography is only practical for small systems (up to about 4-5 qubits).

Single-qubit tomography

Start with a simple example: a qubit prepared in |+⟩ = (|0⟩ + |1⟩)/√2.

# Requires: qiskit_experiments
# The circuit we want to characterise
qc = QuantumCircuit(1)
qc.h(0)

# Perfect simulator -- should recover |+> exactly
backend_ideal = AerSimulator()

# Run tomography
exp = StateTomography(qc)
exp_data = exp.run(backend_ideal, shots=2000).block_for_results()

# Fit the density matrix
state_result = exp_data.analysis_results("state")
rho = state_result.value  # DensityMatrix object

print("Reconstructed density matrix:")
print(np.round(rho.data, 4))

# Ideal density matrix for |+>
ideal = DensityMatrix(Statevector([1/np.sqrt(2), 1/np.sqrt(2)]))
fidelity = state_fidelity(ideal, rho)
print(f"\nState fidelity: {fidelity:.4f}")
# Should be very close to 1.0 on an ideal simulator

Adding noise: what bad preparation looks like

# Requires: qiskit_experiments
# Noisy backend: 5% depolarizing on H gate
noise = NoiseModel()
noise.add_all_qubit_quantum_error(depolarizing_error(0.05, 1), ['h'])
backend_noisy = AerSimulator(noise_model=noise)

exp_noisy = StateTomography(qc)
exp_data_noisy = exp_noisy.run(backend_noisy, shots=4000).block_for_results()
rho_noisy = exp_data_noisy.analysis_results("state").value

print("Noisy reconstructed density matrix:")
print(np.round(rho_noisy.data, 4))

fidelity_noisy = state_fidelity(ideal, rho_noisy)
print(f"\nState fidelity (noisy): {fidelity_noisy:.4f}")
# Expect ~0.92-0.96 with 5% depolarizing
purity = np.real(np.trace(rho_noisy.data @ rho_noisy.data))
print(f"Purity: {purity:.4f}")
# Purity < 1 confirms the state is mixed (noise degraded it)

Two-qubit tomography: Bell state

Two-qubit tomography is the most common use case; verifying Bell state preparation fidelity is a standard hardware benchmark.

# Requires: qiskit_experiments
# Bell state: (|00> + |11>) / sqrt(2)
qc_bell = QuantumCircuit(2)
qc_bell.h(0)
qc_bell.cx(0, 1)

# Noise: 2% depolarizing on CX gate (realistic for superconducting hardware)
noise2q = NoiseModel()
noise2q.add_all_qubit_quantum_error(depolarizing_error(0.02, 2), ['cx'])
backend_2q = AerSimulator(noise_model=noise2q)

exp_bell = StateTomography(qc_bell)
exp_data_bell = exp_bell.run(backend_2q, shots=8000).block_for_results()
rho_bell = exp_data_bell.analysis_results("state").value

# Ideal Bell state density matrix
bell_ideal = DensityMatrix(qc_bell)
fid_bell = state_fidelity(bell_ideal, rho_bell)
purity_bell = np.real(np.trace(rho_bell.data @ rho_bell.data))

print(f"Bell state fidelity: {fid_bell:.4f}")
print(f"Bell state purity:   {purity_bell:.4f}")

# Off-diagonal elements (coherences) indicate entanglement
# A high-fidelity Bell state has rho[0,3] = rho[3,0] = 0.5
print(f"\nCoherence rho[00,11]: {rho_bell.data[0,3]:.4f}")

Reading fidelity numbers

Fidelity ranges from 0 (completely wrong state) to 1 (perfect):

Fidelity	Interpretation
> 0.99	Excellent, top-tier hardware
0.95-0.99	Good, suitable for most NISQ algorithms
0.90-0.95	Moderate, usable but error mitigation needed
< 0.90	Poor for Bell states; hardware needs investigation

Published two-qubit gate fidelities from hardware providers (typically 99-99.9%) are process fidelities under ideal conditions. State tomography on circuits you prepare often gives lower numbers because it includes state preparation and measurement (SPAM) errors on top of gate errors.

Linear inversion vs. maximum likelihood estimation

Reconstructing a density matrix from measurement data is an inverse problem. There are two standard approaches, and qiskit-experiments uses both in sequence.

Linear inversion

Given the set of Pauli expectation values $\langle P_k \rangle$ estimated from measurement data, the density matrix can be written directly:

$$\rho_{\text{lin}} = \frac{1}{2^n} \sum_k \langle P_k \rangle \, P_k$$

where the sum runs over all $4^n$ Pauli operators (including the identity). This is a simple, closed-form solution. The problem is that with finite shots, the estimated expectation values carry statistical noise. This noise can push $\rho_{\text{lin}}$ outside the space of valid density matrices: some eigenvalues may become negative, violating the positive semidefiniteness requirement.

Linear inversion is fast and unbiased, but the result is not guaranteed to be physical.

Maximum likelihood estimation (MLE)

MLE finds the physical density matrix that best explains the observed data. The optimization problem is:

$$\text{maximize} \quad \sum_k n_k \log\left(\text{Tr}(\rho \, \Pi_k)\right)$$

subject to $\rho \geq 0$ (positive semidefinite) and $\text{Tr}(\rho) = 1$, where $n_k$ is the number of times outcome $k$ was observed and $\Pi_k$ is the corresponding measurement projector. This is equivalent to maximizing the log-likelihood of the observed data under the model parametrized by $\rho$.

In an equivalent formulation using Pauli expectation values, MLE solves a least-squares problem: find the valid $\rho$ that minimizes the sum of squared differences between the predicted and observed expectation values.

Why qiskit-experiments uses both

The StateTomography analysis pipeline in qiskit-experiments runs linear inversion first to get a fast initial estimate, then uses that estimate as the starting point for MLE optimization. This two-step approach combines the speed of linear inversion with the physicality guarantee of MLE. The linear inversion result provides a warm start that helps the MLE optimizer converge quickly.

# Requires: qiskit_experiments
# Demonstrate the difference: linear inversion can produce negative eigenvalues
# Run with very few shots to see the effect
exp_few = StateTomography(qc)
exp_data_few = exp_few.run(backend_noisy, shots=50).block_for_results()
rho_few = exp_data_few.analysis_results("state").value

eigenvalues = np.linalg.eigvalsh(rho_few.data)
print("Eigenvalues of MLE-reconstructed state (guaranteed non-negative):")
print(np.round(eigenvalues, 6))
# All eigenvalues should be >= 0 because MLE enforces physicality

Shot count requirements and statistical uncertainty

The precision of the reconstructed density matrix depends on the number of measurement shots. Each matrix element is estimated from counting statistics, so the standard deviation in each element scales as $O(1/\sqrt{N})$ where $N$ is the number of shots per measurement setting.

For state fidelity, the uncertainty also scales as $O(1/\sqrt{N})$. To achieve a fidelity uncertainty of $\delta$, you need roughly $N \sim 1/\delta^2$ shots per setting.

Practical shot budget

Qubits	Settings ($4^n - 1$)	Shots per setting for $\delta = 0.01$	Total shots
1	3	~10,000	~30,000
2	15	~10,000	~150,000
3	63	~10,000	~630,000
4	255	~10,000	~2,550,000

These numbers illustrate why tomography is impractical beyond a few qubits: both the number of settings and the shots per setting must be large.

Fidelity vs. shots: a practical experiment

# Requires: qiskit_experiments
# Run Bell state tomography with increasing shot counts and observe fidelity convergence
shot_counts = [100, 1000, 4000, 16000]
bell_ideal_dm = DensityMatrix(qc_bell)

print(f"{'Shots':>8}  {'Fidelity':>10}  {'Purity':>8}  {'Fid. error':>10}")
print("-" * 44)

for shots in shot_counts:
    exp_s = StateTomography(qc_bell)
    data_s = exp_s.run(backend_2q, shots=shots).block_for_results()
    rho_s = data_s.analysis_results("state").value
    fid_s = state_fidelity(bell_ideal_dm, rho_s)
    pur_s = np.real(np.trace(rho_s.data @ rho_s.data))
    fid_err = abs(fid_s - 0.98)  # approximate expected fidelity with 2% depol
    print(f"{shots:>8}  {fid_s:>10.4f}  {pur_s:>8.4f}  {fid_err:>10.4f}")

# The fidelity error column should decrease roughly as 1/sqrt(shots)

Pauli expectation values from tomography data

The reconstructed density matrix contains all information about the state, including the expectation value of any observable. For Pauli operators this is straightforward.

Single-qubit example

For a single-qubit state $\rho$, the expectation value of Pauli $P$ is $\langle P \rangle = \text{Tr}(\rho \, P)$.

# Requires: qiskit_experiments
# Reconstruct |+> and extract Pauli expectations
exp_plus = StateTomography(qc)
data_plus = exp_plus.run(AerSimulator(), shots=8000).block_for_results()
rho_plus = data_plus.analysis_results("state").value

# Compute expectations manually
X = np.array([[0, 1], [1, 0]])
Y = np.array([[0, -1j], [1j, 0]])
Z = np.array([[1, 0], [0, -1]])

exp_x = np.real(np.trace(rho_plus.data @ X))
exp_y = np.real(np.trace(rho_plus.data @ Y))
exp_z = np.real(np.trace(rho_plus.data @ Z))

print(f"<X> = {exp_x:.4f}  (expect +1 for |+>)")
print(f"<Y> = {exp_y:.4f}  (expect  0 for |+>)")
print(f"<Z> = {exp_z:.4f}  (expect  0 for |+>)")

Two-qubit example using SparsePauliOp

For multi-qubit states, Qiskit’s SparsePauliOp makes this clean:

# Requires: qiskit_experiments
# Extract Pauli expectations from the Bell state density matrix
exp_bell2 = StateTomography(qc_bell)
data_bell2 = exp_bell2.run(AerSimulator(), shots=8000).block_for_results()
rho_bell2 = data_bell2.analysis_results("state").value

# For a Bell state (|00> + |11>)/sqrt(2), the expected values are:
# <XX> = +1, <YY> = -1, <ZZ> = +1, <XY> = 0, <ZI> = 0, <IZ> = 0
paulis = ["XX", "YY", "ZZ", "XY", "ZI", "IZ"]
print("Pauli expectations from reconstructed Bell state:")
for p in paulis:
    op = SparsePauliOp(p)
    val = np.real(rho_bell2.expectation_value(op))
    print(f"  <{p}> = {val:+.4f}")

The expectation values $\langle XX \rangle = +1$, $\langle YY \rangle = -1$, $\langle ZZ \rangle = +1$ are the signature of the Bell state $(|00\rangle + |11\rangle)/\sqrt{2}$. Any deviation from these ideal values quantifies how much noise has affected the state.

Concurrence and entanglement measure

For two-qubit states, the concurrence $C(\rho)$ provides a direct measure of entanglement. It ranges from 0 (separable) to 1 (maximally entangled). The formula due to Wootters is:

$$C(\rho) = \max\left(0,\; \sqrt{\lambda_1} - \sqrt{\lambda_2} - \sqrt{\lambda_3} - \sqrt{\lambda_4}\right)$$

where $\lambda_1 \geq \lambda_2 \geq \lambda_3 \geq \lambda_4$ are the eigenvalues of the matrix

$$R = \rho \, (Y \otimes Y) \, \rho^* \, (Y \otimes Y)$$

and $\rho^*$ denotes the complex conjugate (not the conjugate transpose) of $\rho$.

# Requires: qiskit_experiments
def concurrence(rho_matrix):
    """Compute the Wootters concurrence of a two-qubit density matrix."""
    # Pauli-Y matrix
    sy = np.array([[0, -1j], [1j, 0]])
    # Y tensor Y
    yy = np.kron(sy, sy)

    # R = rho * (Y⊗Y) * rho_conjugate * (Y⊗Y)
    rho_star = np.conjugate(rho_matrix)
    R = rho_matrix @ yy @ rho_star @ yy

    # Eigenvalues of R, sorted descending
    eigvals = np.sort(np.real(np.linalg.eigvals(R)))[::-1]
    # Clamp small negative values from numerical error
    eigvals = np.maximum(eigvals, 0)
    sqrt_eigvals = np.sqrt(eigvals)

    return max(0, sqrt_eigvals[0] - sqrt_eigvals[1] - sqrt_eigvals[2] - sqrt_eigvals[3])

# Bell state: expect concurrence near 1
exp_c_bell = StateTomography(qc_bell)
data_c_bell = exp_c_bell.run(AerSimulator(), shots=8000).block_for_results()
rho_c_bell = data_c_bell.analysis_results("state").value
C_bell = concurrence(rho_c_bell.data)
print(f"Concurrence of Bell state: {C_bell:.4f}  (expect ~1.0)")

# Product state |00>: expect concurrence = 0
qc_product = QuantumCircuit(2)  # |00> by default
exp_c_prod = StateTomography(qc_product)
data_c_prod = exp_c_prod.run(AerSimulator(), shots=8000).block_for_results()
rho_c_prod = data_c_prod.analysis_results("state").value
C_prod = concurrence(rho_c_prod.data)
print(f"Concurrence of |00>:       {C_prod:.4f}  (expect ~0.0)")

Concurrence is useful when you want a single number that quantifies “how entangled” a reconstructed state is. Note that it is defined only for two-qubit states; for larger systems, different entanglement measures (such as entanglement entropy or negativity) are needed.

Partial tomography: subsystem states

If you only care about a subsystem, you can prepare and measure only that part, then use partial_trace to get the reduced density matrix. This is more efficient than full tomography.

# Requires: qiskit_experiments
# Prepare a 3-qubit GHZ state but only characterise qubits 0 and 1
qc_ghz = QuantumCircuit(3)
qc_ghz.h(0)
qc_ghz.cx(0, 1)
qc_ghz.cx(0, 2)

# Full tomography on all 3 qubits
exp_ghz = StateTomography(qc_ghz)
exp_data_ghz = exp_ghz.run(AerSimulator(), shots=4000).block_for_results()
rho_ghz = exp_data_ghz.analysis_results("state").value

# Partial tomography: only on qubits [0, 1]
# qiskit-experiments supports targeting specific qubits
exp_partial = StateTomography(qc_ghz, measurement_indices=[0, 1])
exp_data_partial = exp_partial.run(AerSimulator(), shots=4000).block_for_results()
rho_01 = exp_data_partial.analysis_results("state").value

# Compare: manual partial trace from full state
rho_01_manual = partial_trace(rho_ghz, [2])  # trace out qubit 2

print("Subsystem [0,1] via partial tomography:")
print(np.round(rho_01.data, 3))
print("\nSubsystem [0,1] via partial trace of full state:")
print(np.round(rho_01_manual.data, 3))

Tomographic completeness and informationally complete measurements

A set of measurement operators is called informationally complete (IC) if the measurement statistics uniquely determine any density matrix. In other words, no two distinct density matrices can produce the same measurement outcome probabilities under an IC measurement.

Pauli measurements are IC

The standard approach in qiskit-experiments uses all $4^n - 1$ Pauli tensor products. This set is informationally complete because the Pauli operators (together with the identity) form a basis for the space of $2^n \times 2^n$ Hermitian matrices. Any density matrix can be uniquely expanded in this basis, and each coefficient is directly an expectation value.

However, Pauli measurements are overcomplete in a sense: you measure each qubit in three bases (X, Y, Z), producing $3^n$ distinct circuit configurations. Each configuration yields $2^n$ outcome probabilities. The total raw data is $3^n \times 2^n$ numbers, which is more than the $4^n - 1$ independent parameters. This redundancy is helpful because it provides statistical cross-checks on the data.

Symmetric IC (SIC) POVMs

The minimal IC measurement for a $d$-dimensional system requires exactly $d^2$ outcomes. For a single qubit ($d = 2$), this means 4 outcomes. A symmetric informationally complete POVM (SIC-POVM) achieves this minimum with measurement operators that are as “evenly spread” as possible in Hilbert space.

For a single qubit, the 4 SIC-POVM elements correspond to the vertices of a tetrahedron inscribed in the Bloch sphere. In contrast, Pauli measurements use 6 projectors (the vertices of an octahedron). SIC-POVMs are theoretically elegant and statistically optimal, but implementing them on hardware requires more complex gate sequences. For practical purposes, Pauli measurements are the standard choice in qiskit-experiments because they use only Clifford gates (H, S, etc.) that are natively available on most hardware platforms.

Visualizing the density matrix

Visualizing the real and imaginary parts of a density matrix helps build intuition about the state structure.

Heatmap visualization with matplotlib

# Requires: qiskit_experiments, matplotlib
import matplotlib.pyplot as plt

# Reconstruct a Bell state for visualization
exp_viz = StateTomography(qc_bell)
data_viz = exp_viz.run(AerSimulator(), shots=8000).block_for_results()
rho_viz = data_viz.analysis_results("state").value

labels = ["|00⟩", "|01⟩", "|10⟩", "|11⟩"]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# Real part
im1 = ax1.imshow(np.real(rho_viz.data), cmap="RdBu", vmin=-0.6, vmax=0.6)
ax1.set_xticks(range(4))
ax1.set_yticks(range(4))
ax1.set_xticklabels(labels)
ax1.set_yticklabels(labels)
ax1.set_title("Re(ρ)")
fig.colorbar(im1, ax=ax1, shrink=0.8)

# Imaginary part
im2 = ax2.imshow(np.imag(rho_viz.data), cmap="RdBu", vmin=-0.6, vmax=0.6)
ax2.set_xticks(range(4))
ax2.set_yticks(range(4))
ax2.set_xticklabels(labels)
ax2.set_yticklabels(labels)
ax2.set_title("Im(ρ)")
fig.colorbar(im2, ax=ax2, shrink=0.8)

plt.suptitle("Bell State Density Matrix")
plt.tight_layout()
plt.savefig("bell_density_matrix.png", dpi=150)
plt.show()

Text-based visualization

For environments without matplotlib, you can print the density matrix as a labeled array:

# Requires: qiskit_experiments
def print_density_matrix(rho_data, n_qubits):
    """Print a density matrix with row and column labels."""
    dim = 2 ** n_qubits
    labels = [format(i, f'0{n_qubits}b') for i in range(dim)]
    labels = [f"|{l}⟩" for l in labels]

    # Header
    header = "       " + "  ".join(f"{l:>10}" for l in labels)
    print(header)

    for i, row_label in enumerate(labels):
        row_str = f"{row_label:>5}  "
        for j in range(dim):
            val = rho_data[i, j]
            re, im = val.real, val.imag
            if abs(im) < 1e-4:
                row_str += f"  {re:>8.4f}  "
            else:
                sign = "+" if im >= 0 else "-"
                row_str += f"{re:>5.3f}{sign}{abs(im):.3f}j "
            row_str += ""
        print(row_str)

print("Bell state density matrix:")
print_density_matrix(rho_viz.data, 2)

For a perfect Bell state $(|00\rangle + |11\rangle)/\sqrt{2}$, you should see large values (around 0.5) at the four corners of the matrix: $\rho_{00,00}$, $\rho_{00,11}$, $\rho_{11,00}$, and $\rho_{11,11}$. The off-diagonal elements $\rho_{00,11}$ and $\rho_{11,00}$ are the coherences that signal entanglement. All other entries should be near zero.

Quantum process tomography (expanded)

State tomography tells you what state a preparation circuit makes. Process tomography tells you what a gate actually does: it characterizes the full quantum channel. For a single qubit, process tomography requires preparing 4 input states and running state tomography on each output, giving 4 x 3 = 12 measurement settings. For two qubits: 16 x 15 = 240 settings.

Full worked example: noisy CX gate

# Requires: qiskit_experiments
from qiskit_experiments.library import ProcessTomography

# Characterise a noisy CX gate
qc_cx = QuantumCircuit(2)
qc_cx.cx(0, 1)

# 2% depolarizing noise on CX
noise_cx = NoiseModel()
noise_cx.add_all_qubit_quantum_error(depolarizing_error(0.02, 2), ['cx'])
backend_cx = AerSimulator(noise_model=noise_cx)

exp_pt = ProcessTomography(qc_cx)
exp_data_pt = exp_pt.run(backend_cx, shots=4000).block_for_results()

# The reconstructed quantum channel as a Choi matrix
choi_result = exp_data_pt.analysis_results("state")
choi_matrix = choi_result.value  # This is the Choi representation of the channel

# Process fidelity
proc_fid = exp_data_pt.analysis_results("process_fidelity").value
print(f"CX process fidelity: {proc_fid:.4f}")

Average gate fidelity from process fidelity

Process fidelity $F_{\text{process}}$ and average gate fidelity $F_{\text{avg}}$ are related by:

$$F_{\text{avg}} = \frac{d \, F_{\text{process}} + 1}{d + 1}$$

where $d = 2^n$ is the Hilbert space dimension. For a two-qubit gate ($d = 4$):

# Requires: qiskit_experiments
d = 4  # 2-qubit Hilbert space dimension
F_avg = (d * proc_fid + 1) / (d + 1)
print(f"Process fidelity:      {proc_fid:.4f}")
print(f"Average gate fidelity: {F_avg:.4f}")
# Average gate fidelity is always higher than process fidelity
# because it includes the trivial identity component

Extracting Kraus operators

The Choi matrix representation of a channel can be converted to Kraus operators, which give a physical picture of the noise. A channel $\mathcal{E}$ acts as $\mathcal{E}(\rho) = \sum_k K_k \rho K_k^\dagger$. For an ideal unitary gate, there is a single Kraus operator equal to the unitary itself.

# Requires: qiskit_experiments
from qiskit.quantum_info import Choi, Kraus

# Convert Choi matrix to Kraus representation
kraus_rep = Kraus(choi_matrix)
kraus_ops = kraus_rep.data

print(f"Number of Kraus operators: {len(kraus_ops)}")
print(f"\nDominant Kraus operator (should resemble CX):")
print(np.round(kraus_ops[0], 3))

# Check: for an ideal gate, the dominant Kraus operator has singular values all equal to 1
svd_vals = np.linalg.svd(kraus_ops[0], compute_uv=False)
print(f"\nSingular values of dominant Kraus op: {np.round(svd_vals, 4)}")

Process fidelity is what hardware providers report as “gate fidelity”; it is more informative than state fidelity for evaluating individual gates.

Gate set tomography: beyond QPT

Quantum process tomography has a fundamental limitation: it assumes that your state preparation and measurement (SPAM) operations are perfect. In reality, SPAM operations are themselves noisy. If you prepare a slightly wrong input state and then characterize a gate using QPT, the reconstructed channel is a convolution of the true gate errors with the SPAM errors. You cannot tell them apart.

Gate set tomography (GST) solves this problem by self-consistently characterizing everything at once: the state preparation, the measurement, and each gate. GST uses carefully designed gate sequences of varying lengths to separate SPAM errors from gate errors through a self-consistency condition.

The key insight is that GST uses “fiducial” circuits (short gate sequences that probe different aspects of the operations) and “germ” circuits (sequences designed to amplify specific error parameters). By running germs at increasing repetition counts, GST can isolate gate errors from SPAM errors.

# Requires: qiskit_experiments
# GST is available in qiskit-experiments for single-qubit gate sets
from qiskit_experiments.library import StandardRB  # For comparison

# Note: GateSetTomography setup requires specifying the target gate set,
# fiducials, and germs. For a single-qubit gate set {Idle, X_pi/2, Y_pi/2}:
#
# from qiskit_experiments.library.tomography import GateSetTomography
# gst_exp = GateSetTomography(physical_qubits=(0,), gate_set="default_1q")
# gst_data = gst_exp.run(backend, shots=4000).block_for_results()
#
# The result provides separately characterized SPAM parameters and gate
# process matrices, giving a more accurate picture than QPT alone.

GST is most valuable during hardware bring-up when you need to understand the detailed error structure of individual gates. It is significantly more expensive than QPT (requiring hundreds of circuits even for a single qubit), so it is not used for routine monitoring.

Randomized benchmarking vs. tomography

Tomography and randomized benchmarking (RB) serve different purposes. The table below provides a practical decision guide.

Criterion	State/Process Tomography	Randomized Benchmarking
Output	Full density matrix or process matrix	Single fidelity number (error per gate)
Information	Complete error characterization	Average error rate only
Scaling	Exponential in qubit count ($4^n$)	Polynomial in qubit count
SPAM sensitivity	Included in the result (a drawback for QPT)	SPAM errors are filtered out by design
Circuit depth	Fixed (short) circuits	Variable-depth random Clifford sequences
Practical qubit limit	4-5 qubits	10+ qubits
When to use	Diagnosing specific gate errors during bring-up	Routine quality checks in production
Error structure	Reveals coherent vs. incoherent, axis of rotation, etc.	Hidden; only the average rate is reported

Practical guidance

Use tomography when:

• You are characterizing a new gate calibration and need to know the specific error channel.
• You need to verify entanglement properties of a prepared state.
• You are debugging why a specific algorithm underperforms on hardware.
• You are working with 1-3 qubits and can afford the shot budget.

Use randomized benchmarking when:

• You need a quick health check of gate quality.
• You are comparing gate fidelity across many qubits on a device.
• You want SPAM-free error rates.
• You are monitoring hardware stability over time.

# Requires: qiskit_experiments
# Quick comparison: RB on the same noisy backend
from qiskit_experiments.library import StandardRB

# Run RB on 1 qubit with sequence lengths [1, 10, 20, 50, 100]
rb_exp = StandardRB(
    physical_qubits=(0,),
    lengths=[1, 10, 20, 50, 100],
    num_samples=5,  # number of random circuits per length
    seed=42
)
rb_data = rb_exp.run(backend_noisy, shots=1000).block_for_results()

# The error per Clifford (EPC) is the key metric
epc = rb_data.analysis_results("EPC").value
print(f"Error per Clifford (EPC): {epc:.5f}")
# For 5% depolarizing on H, each Clifford averages ~1.5 H gates,
# so expect EPC ~ 0.05 * 1.5 / 2 ~ 0.0375 (approximate)

When NOT to use tomography

State tomography is expensive and has specific failure modes:

Do not use it to verify algorithm correctness for large circuits. With n > 5 qubits the exponential shot cost makes it impractical. Use instead: spot-checking via known observable expectations, or algorithmic verification techniques.

Do not confuse state fidelity with algorithm performance. A prepared Bell state with 0.97 fidelity does not mean an algorithm using that state has 0.97 performance; errors compound through the circuit.

Tomography cannot detect all errors. If your state preparation always introduces the same coherent error (e.g., always over-rotates by the same angle), tomography will faithfully reconstruct the wrong state. For systematic gate errors you need process tomography, not state tomography.

The MLE fit can fail. Maximum likelihood estimation assumes the measurement data is consistent with a valid quantum state. With very few shots or pathological noise, the fit can produce an unphysical result. StateTomography uses LinearInversion followed by MLE to guarantee physicality, but the result can still have high uncertainty.

Common mistakes

These are the most frequent errors that practitioners encounter when using quantum state tomography. Avoiding them will save significant debugging time.

1. Forgetting that tomography requires qiskit-experiments

The StateTomography and ProcessTomography classes are not part of core Qiskit. They live in the qiskit-experiments package, which must be installed separately. If you see ModuleNotFoundError: No module named 'qiskit_experiments', run pip install qiskit-experiments.

2. Not running enough shots for the number of qubits

Each additional qubit quadruples the number of measurement settings. If you use the same total shot budget for 3 qubits as you used for 1 qubit, each individual setting gets roughly 20x fewer shots. The result is a noisy reconstruction with low confidence. A good rule of thumb: allocate at least 1000 shots per measurement setting, which means at least 4000 total shots for 1 qubit, 15,000 for 2 qubits, and 63,000 for 3 qubits.

3. Confusing state fidelity with process fidelity

State fidelity measures how close a prepared state is to a target state. Process fidelity measures how close a quantum channel is to a target channel. They are related by the formula $F_{\text{avg}} = (d \cdot F_{\text{process}} + 1)/(d+1)$, but they answer different questions. When comparing your tomography results with published gate fidelities from hardware providers, make sure you are comparing the same quantity.

4. Reporting purity < 1 as evidence of entanglement

Purity less than 1 means the state is mixed, not that it is entangled. A completely unentangled state can have low purity (e.g., a single qubit subjected to depolarizing noise). Conversely, a pure product state has purity 1 and zero entanglement. To quantify entanglement, use concurrence, negativity, or entanglement entropy, not purity.

5. Applying tomography to verify large circuit output

It can be tempting to use tomography to verify the output of a 10-qubit variational algorithm. With $4^{10} - 1 = 1,048,575$ independent settings, this is not feasible. For large circuits, measure specific observables that your algorithm cares about rather than reconstructing the full state. Qiskit’s Estimator primitive is designed for exactly this use case.

6. Forgetting to set measurement_indices for subsystem tomography

If you have a 4-qubit circuit but only care about the state of qubits 0 and 1, running StateTomography(circuit) will perform full 4-qubit tomography (255 settings). Instead, use StateTomography(circuit, measurement_indices=[0, 1]) to perform 2-qubit tomography on just the qubits you need (15 settings). This reduces both circuit count and shot budget by a large factor.

Practical summary

• Use state tomography to verify small (1-3 qubit) state preparation circuits; Bell state fidelity is the most common benchmark
• Use process tomography to characterize individual gates on real hardware
• Use gate set tomography when you need to separate SPAM errors from gate errors during device bring-up
• Use randomized benchmarking for routine, scalable quality checks
• Interpret fidelity > 0.95 for Bell states as “hardware is working well”
• Watch purity alongside fidelity: purity tells you how much incoherent noise (depolarizing, decoherence) is present vs. coherent errors
• Use concurrence to quantify entanglement in reconstructed two-qubit states
• For large circuits, prefer expectation-value benchmarks (quantum volume, mirror circuits) over tomography
• Always check that you have enough shots per measurement setting, not just enough total shots