Quantum State Tomography with Qiskit

Quantum state tomography (QST) is the procedure for reconstructing the full description of a quantum state from many repeated measurements. A qubit in an unknown state cannot be fully characterized by measuring it once; a single measurement collapses the state and returns only a classical bit (0 or 1). But if you prepare many identical copies of the state and measure each copy in a different basis, you can piece together the complete density matrix.

For a single qubit, you need to measure in the X, Y, and Z Pauli bases. Each basis measurement, repeated over many shots, gives you an expectation value. From these three expectation values you recover the three independent real parameters of the single-qubit density matrix (the fourth parameter is fixed by the normalization condition Tr(rho) = 1). For n qubits, you need measurements in all combinations of Pauli bases, which gives 3^n measurement settings and 4^n - 1 independent real parameters to estimate. This exponential scaling is the fundamental reason QST becomes impractical for more than a few qubits.

QST is an essential tool for debugging circuits, verifying entanglement, and characterizing quantum devices. It tells you everything about the state, not just one property of it.

When to Use State Tomography

QST is most useful when you need to:

• Verify that a circuit prepares the state you expect
• Compute state fidelity against an ideal target
• Diagnose SPAM (state preparation and measurement) errors
• Measure entanglement metrics like concurrence or purity
• Distinguish between coherent errors (wrong unitary) and incoherent errors (decoherence)

The cost scales as 3^n measurement settings for n qubits, so it is practical only for small systems (1-3 qubits in most real workflows). For larger systems, consider shadow tomography or randomized benchmarking instead.

What is a Density Matrix?

Most introductory quantum computing courses teach the statevector representation: a qubit is described by |psi> = alpha|0> + beta|1>, where alpha and beta are complex amplitudes. This works perfectly for pure states, where you have complete knowledge of the system. But quantum states can also be mixed, meaning the system is in a statistical ensemble of pure states rather than a single definite superposition.

Mixed states arise in three common situations:

1. Decoherence: When qubits interact with their environment, pure states evolve into mixed states. This is the dominant source of errors on real hardware.
2. Partial tracing: If two qubits are entangled and you only have access to one of them, that qubit is in a mixed state even though the joint system may be pure.
3. Noisy preparation: If your state preparation is imperfect and sometimes produces |psi_1> and sometimes |psi_2>, the resulting ensemble is a mixed state.

A statevector cannot represent mixed states. The density matrix can represent both pure and mixed states in a unified formalism. For a pure state |psi>, the density matrix is the outer product rho = |psi><psi|. For a mixed state that is |psi_1> with probability p_1 and |psi_2> with probability p_2, the density matrix is rho = p_1 |psi_1><psi_1| + p_2 |psi_2><psi_2|.

Key properties of a valid density matrix:

• Hermitian: rho = rho^dagger (ensures real expectation values)
• Positive semidefinite: all eigenvalues are >= 0 (ensures valid probabilities)
• Trace one: Tr(rho) = 1 (probabilities sum to 1)
• Pure states have rank 1: rho^2 = rho, so Tr(rho^2) = 1
• Mixed states have rank > 1: Tr(rho^2) < 1

For the Bell state |Phi+> = (|00> + |11>) / sqrt(2), the density matrix is a 4x4 matrix. Writing it out explicitly in the computational basis {|00>, |01>, |10>, |11>}:

rho = |Phi+><Phi+| = 1/2 * (|00> + |11>)(<00| + <11|)

     |00>  |01>  |10>  |11>
|00> [1/2    0     0    1/2]
|01> [ 0     0     0     0 ]
|10> [ 0     0     0     0 ]
|11> [1/2    0     0    1/2]

The off-diagonal elements (positions (0,3) and (3,0)) encode the quantum coherence between |00> and |11>. If those off-diagonal elements decay to zero (due to decoherence), the state becomes the classical mixture rho = 1/2 |00><00| + 1/2 |11><11|, which has the same measurement probabilities in the Z basis but no entanglement.

How State Tomography Works

Here is the step-by-step protocol for 2-qubit state tomography:

Step 1: Choose measurement bases. For 2 qubits, you measure in all combinations of the Pauli bases {X, Y, Z} for each qubit. This gives 3^2 = 9 measurement settings: XX, XY, XZ, YX, YY, YZ, ZX, ZY, ZZ. The notation “XY” means “measure qubit 0 in the X basis and qubit 1 in the Y basis.”

Measuring in the X basis means applying a Hadamard gate before measurement. Measuring in the Y basis means applying S-dagger then Hadamard before measurement. Measuring in the Z basis (the computational basis) requires no additional gates.

Step 2: Run shots for each basis. For each of the 9 measurement settings, prepare the state and measure many times (thousands of shots). Each run produces a histogram of bitstring counts. For example, measuring in the ZZ basis for the Bell state |Phi+> gives roughly 50% |00> and 50% |11>.

Step 3: Estimate expectation values. From the bitstring counts, compute the expectation value of each Pauli operator. A 2-qubit density matrix has 4^2 - 1 = 15 independent real parameters (the 16th is fixed by Tr(rho) = 1). The 9 measurement bases give you access to all 15 independent parameters because each basis yields information about multiple Pauli operators simultaneously.

Step 4: Reconstruct the density matrix. Use the expectation values to reconstruct rho. The standard approach is maximum likelihood estimation (MLE), which finds the density matrix that maximizes the probability of observing the measured data, subject to the constraint that rho is a valid density matrix (Hermitian, positive semidefinite, trace 1).

MLE is preferred over simple linear inversion because linear inversion can produce unphysical results (density matrices with negative eigenvalues). MLE always returns a valid quantum state.

How Qiskit Experiments handles this. The StateTomography class in qiskit-experiments automates all four steps. It generates the measurement circuits by appending the appropriate basis-change gates before measurement, submits them as a batch job, collects the results, and runs MLE reconstruction. You can inspect the measurement circuits it generates:

from qiskit import QuantumCircuit
from qiskit_experiments.library import StateTomography

# Prepare a Bell state
qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)

# Create the tomography experiment
qst = StateTomography(qc)

# See the measurement circuits
circuits = qst.circuits()
print(f"Number of measurement circuits: {len(circuits)}")
for i, circ in enumerate(circuits):
    print(f"\nCircuit {i}:")
    print(circ.draw())

This produces 9 circuits (for 2 qubits), each with the Bell state preparation followed by different basis-rotation gates before the final measurement.

Installation

pip install qiskit qiskit-aer qiskit-experiments

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

Running State Tomography

Qiskit Experiments provides a ready-made StateTomography class:

from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit_experiments.library import StateTomography

# Prepare a Bell state
qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)

# Create the tomography experiment
qst = StateTomography(qc)

# Run on the simulator
backend = AerSimulator()
job = qst.run(backend, shots=4096)
result = job.block_for_results()

# Extract the reconstructed density matrix
state_result = result.analysis_results("state")
dm = state_result.value
print(dm)

The returned value is a DensityMatrix object from qiskit.quantum_info. This object stores the full 4x4 complex matrix for a 2-qubit state and provides methods for computing fidelity, purity, and other properties.

Computing State Fidelity

Once you have the reconstructed density matrix, compare it to the ideal state. State fidelity F(rho, sigma) quantifies how similar two quantum states are. For two general density matrices, the fidelity is defined as:

F(rho, sigma) = (Tr(sqrt(sqrt(rho) * sigma * sqrt(rho))))^2

This formula reduces to a simpler expression when one of the states is pure. If sigma = |psi><psi| is a pure state, then:

F(rho, |psi>) = <psi| rho |psi>

which is just the expectation value of the density matrix in the target state. Fidelity ranges from 0 (orthogonal states) to 1 (identical states).

On an ideal noiseless simulator, the reconstructed Bell state fidelity is typically 0.97 to 0.99 rather than exactly 1.0. The deviation comes from statistical shot noise in the tomography: you are estimating expectation values from a finite number of measurements, and those estimates have sampling error. Increasing the number of shots reduces this error.

from qiskit.quantum_info import DensityMatrix, state_fidelity, Statevector

# Construct the ideal Bell state directly from amplitudes.
# |Phi+> = (|00> + |11>) / sqrt(2)
# In Qiskit's little-endian ordering, the statevector is [1/sqrt(2), 0, 0, 1/sqrt(2)].
bell = Statevector([0.5**0.5, 0, 0, 0.5**0.5])
ideal_dm = DensityMatrix(bell)

# Compare the reconstructed density matrix (dm) to the ideal
fidelity = state_fidelity(dm, ideal_dm)
print(f"State fidelity: {fidelity:.4f}")

A fidelity of 1.0 means perfect agreement. On real hardware, Bell state fidelities typically range from 0.90 to 0.99 depending on the device, the qubit pair, and calibration quality.

You can also investigate how fidelity improves with more shots:

shot_counts = [1024, 2048, 4096, 8192, 16384]
fidelities = []

for shots in shot_counts:
    job = StateTomography(qc).run(backend, shots=shots)
    result = job.block_for_results()
    dm_test = result.analysis_results("state").value
    fidelities.append(state_fidelity(dm_test, ideal_dm))
    print(f"Shots: {shots:>6}, Fidelity: {fidelities[-1]:.4f}")

You should see fidelity converge toward 1.0 as the shot count increases, since the statistical sampling error decreases as 1/sqrt(shots).

Purity and Entanglement

The density matrix gives you access to several important physical quantities beyond fidelity.

Purity

Purity is defined as Tr(rho^2). It tells you how mixed a quantum state is:

• Pure states have purity = 1.0. The state is a single definite superposition with no classical uncertainty.
• Maximally mixed states have purity = 1/d, where d is the Hilbert space dimension. For a 2-qubit system (d = 4), the maximally mixed state has purity 0.25. This is the state of maximum ignorance, equivalent to a uniform classical distribution over all basis states.
• Intermediate values indicate partial decoherence. A purity below 0.9 for a state you expect to be pure is a clear sign of significant decoherence or noise.

Purity is useful as a quick diagnostic. If your circuit is supposed to produce a pure state and the reconstructed purity is 0.85, something is wrong with the preparation or the qubits are decohering during the circuit.

Concurrence

Concurrence C is an entanglement measure for 2-qubit states. It ranges from 0 to 1:

• C = 0: the state is a product state (separable), with no entanglement. Example: |00>.
• C = 1: the state is maximally entangled. Example: the Bell state |Phi+>.
• 0 < C < 1: the state has partial entanglement, which can arise from noisy preparation of an entangled state.

Concurrence is one of several entanglement measures. Others include entanglement entropy (the von Neumann entropy of the reduced density matrix) and negativity (based on the partial transpose). Concurrence is convenient for 2-qubit systems because it has a closed-form expression, but it does not generalize easily to more qubits.

from qiskit.quantum_info import DensityMatrix, concurrence
import numpy as np

# Purity: Tr(rho^2). Pure states have purity=1, mixed states <1
purity = np.real(np.trace(dm.data @ dm.data))
print(f"Purity: {purity:.4f}")

# For 2-qubit states, compute concurrence as an entanglement measure
c = concurrence(dm)
print(f"Concurrence: {c:.4f}")

# For comparison, run tomography on a product state
product_qc = QuantumCircuit(2)
product_qc.h(0)  # |+0> state, no entanglement

qst_product = StateTomography(product_qc)
job_product = qst_product.run(backend, shots=8192)
result_product = job_product.block_for_results()
dm_product = result_product.analysis_results("state").value

purity_product = np.real(np.trace(dm_product.data @ dm_product.data))
c_product = concurrence(dm_product)
print(f"\nProduct state |+0>:")
print(f"  Purity: {purity_product:.4f}")
print(f"  Concurrence: {c_product:.4f}")

For the Bell state, you should see purity close to 1.0 and concurrence close to 1.0. For the product state |+0>, you should see purity close to 1.0 but concurrence close to 0.0.

Including SPAM Error Analysis

SPAM errors are systematic biases from imperfect state preparation and measurement. Unlike shot noise (which averages out with more measurements), SPAM errors are present in every shot and do not decrease with more data. They come from two sources:

1. State preparation errors: the qubit does not start in exactly |0>. Residual thermal population or miscalibrated initialization leaves a small |1> component.
2. Measurement errors: the readout misidentifies |0> as |1> or vice versa. Typical readout error rates on superconducting hardware range from 0.5% to 5% per qubit.

You can assess SPAM by running tomography on the simplest possible state (the ground state) and checking the reconstructed density matrix:

# Tomography on |0><0| - should give perfect result if no SPAM
ground_qc = QuantumCircuit(1)  # no gates, prepare |0>
qst_ground = StateTomography(ground_qc)
job = qst_ground.run(backend, shots=8192)
result = job.block_for_results()
ground_dm = result.analysis_results("state").value

ideal_ground = DensityMatrix.from_label("0")
spam_fidelity = state_fidelity(ground_dm, ideal_ground)
print(f"Ground state fidelity (SPAM proxy): {spam_fidelity:.4f}")

A low fidelity here indicates that the SPAM contribution is large, and you should interpret other tomography results with that in mind.

Using SPAM Fidelity to Interpret Other Results

The SPAM fidelity gives you a rough baseline for how much of the infidelity in other experiments comes from SPAM rather than from your circuit. Consider this example:

• Ground state |0> fidelity: 0.98 (2% loss from SPAM alone)
• Bell state fidelity: 0.95 (5% total loss)

Since the Bell state circuit contains two gates (H and CNOT) on top of the state preparation and measurement, the additional 3 percentage points of loss beyond the SPAM baseline come from gate errors and decoherence during the circuit. The 2 percentage points of SPAM loss appear in both experiments because SPAM affects every measurement equally.

This is a rough decomposition, not an exact one. SPAM errors can interact with gate errors in nonlinear ways, and multi-qubit SPAM is not simply the sum of single-qubit SPAM errors. For rigorous error budgeting, use interleaved randomized benchmarking or gate set tomography. But the SPAM baseline is a fast and useful sanity check.

# Run SPAM analysis on both qubits independently
for qubit_index in range(2):
    single_qc = QuantumCircuit(1)
    qst_single = StateTomography(single_qc)
    job = qst_single.run(backend, shots=8192)
    result = job.block_for_results()
    single_dm = result.analysis_results("state").value
    f = state_fidelity(single_dm, DensityMatrix.from_label("0"))
    print(f"Qubit {qubit_index} SPAM fidelity: {f:.4f}")

Reconstruction Methods

There are several algorithms for reconstructing a density matrix from measurement data. They differ in speed, physical validity of the result, and the additional information they provide.

Linear Inversion

The simplest approach. You directly solve for the density matrix elements from the measured expectation values using a linear system of equations. This is fast (a single matrix inversion), but the result is not guaranteed to be a valid density matrix. Sampling noise can cause the reconstructed matrix to have negative eigenvalues, which is unphysical. Linear inversion is useful for quick checks but should not be used for final results.

Maximum Likelihood Estimation (MLE)

MLE finds the density matrix that maximizes the likelihood of the observed measurement data, subject to the constraints that rho is Hermitian, positive semidefinite, and has trace 1. This always produces a physically valid density matrix. The optimization is convex (after a suitable parameterization) and converges reliably for small systems. Qiskit Experiments uses MLE by default, and this is the standard choice for most applications.

Bayesian Estimation

Bayesian methods treat the density matrix as a random variable with a prior distribution. Instead of a single point estimate, they produce a posterior distribution over density matrices, which gives you uncertainty estimates for every property you compute. The downside is computational cost: sampling the posterior requires Markov chain Monte Carlo or similar methods, which are much slower than MLE. Bayesian estimation is valuable when you need rigorous error bars.

Shadow Tomography (Classical Shadows)

A modern approach that uses random single-qubit rotations before measurement rather than the fixed Pauli bases of standard QST. From a modest number of random measurements, you can efficiently estimate many different properties of the state (expectation values, fidelities, entanglement witnesses) without ever reconstructing the full density matrix. Classical shadows require far fewer total measurements than full QST for predicting linear properties, making them practical for larger systems. The tradeoff is that you get estimates of specific properties rather than the complete density matrix.

You can select the reconstruction method in Qiskit Experiments by configuring the analysis:

from qiskit_experiments.library import StateTomography

qst = StateTomography(qc)

# Use linear inversion instead of the default MLE
qst.analysis.set_options(fitter="linear_inversion")

# Run as usual
job = qst.run(backend, shots=8192)
result = job.block_for_results()
dm_linear = result.analysis_results("state").value
print(dm_linear)

For production use, stick with the default MLE unless you have a specific reason to change it.

Running on Real Hardware

Requires hardware credentials: This block connects to IBM Quantum. Skip if running locally without a registered IBM Quantum account.

On real hardware, use a fake backend first to verify your workflow, then switch to the IBM Runtime:

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False, min_num_qubits=2)

qst = StateTomography(qc)
job = qst.run(backend, shots=8192)
result = job.block_for_results()
dm = result.analysis_results("state").value
print(f"Fidelity: {state_fidelity(dm, ideal_dm):.4f}")

When running on real hardware, keep these points in mind:

• Use at least 8192 shots per basis. Real hardware has more noise than simulators, so you need more statistics to get reliable estimates.
• Check which qubit pairs have the best CNOT error rates before choosing your physical qubits. The backend.properties() method gives you gate error rates.
• Queue times vary. A tomography experiment submits multiple circuits as a single job, so it only occupies one slot in the queue.
• Transpilation matters. Let Qiskit transpile your circuit to the hardware’s native gate set and connectivity. The StateTomography class handles this automatically when you pass a real backend.

Common Mistakes

Using too few shots

With 1024 shots per basis, the statistical uncertainty in each expectation value is roughly 1/sqrt(1024) = 0.03, which propagates into a fidelity uncertainty of several percent. Use at least 4096 shots for 1-qubit tomography and 8192 or more for 2-qubit tomography. If you need fidelity precision below 0.01, you may need 32768 shots or more.

Ignoring statistical uncertainty

The fidelity, purity, and concurrence values you compute from a reconstructed density matrix are estimates, not exact values. They have statistical error bars from the finite number of shots. Qiskit Experiments provides confidence intervals in the analysis results:

state_result = result.analysis_results("state")
# The extra field contains fit metadata including standard error
print(f"State fidelity result: {result.analysis_results('state_fidelity')}")

Always report uncertainties when publishing or comparing tomography results.

Running QST on circuits with mid-circuit measurements

If your circuit includes a measurement gate before the end (a mid-circuit measurement), that measurement collapses the quantum state. The tomography protocol appends its own basis-change gates and measurements at the end of the circuit. If the state has already been partially collapsed by a mid-circuit measurement, the reconstruction will not reflect the pre-measurement quantum state. Remove all mid-circuit measurements from the circuit before passing it to StateTomography.

Using QST when a simpler protocol would suffice

Full state tomography reconstructs the entire density matrix, which is expensive. If you only need one number, there are cheaper alternatives:

• Direct fidelity estimation (DFE): Estimates fidelity to a target state using O(1) measurement settings (independent of system size), rather than the 3^n settings of QST. Use DFE when you only care about fidelity to a known target.
• Randomized benchmarking (RB): Estimates average gate fidelity by running random sequences of Clifford gates. Use RB when you want to characterize gate quality rather than a specific state.
• Pauli expectation values: If you only need to know a few specific expectation values (e.g., ), measure those directly instead of reconstructing the full state.

Choose the tool that matches your question. QST answers “what is the full quantum state?”, which is more information than most practical questions require.

Practical Limits

State tomography is exponentially expensive. For 1 qubit you need 3 Pauli bases; for 2 qubits, 9; for 3 qubits, 27; for n qubits, 3^n. The number of real parameters to estimate grows as 4^n - 1. Beyond 3-4 qubits, full state tomography becomes impractical in both measurement cost and classical post-processing time (the density matrix for 10 qubits is a 1024 x 1024 complex matrix with over a million parameters).

For larger systems, consider:

• Classical shadow tomography: Requires O(log(M) / epsilon^2) measurements to predict M linear properties to precision epsilon, regardless of system size.
• Matrix product state tomography: Exploits the structure of states with limited entanglement (low bond dimension) to reduce the number of measurements needed.
• Randomized benchmarking: Characterizes average gate quality without full state reconstruction.

QST is best used as a diagnostic tool during circuit development rather than as part of a production algorithm.