Chapter 8: Programming Quantum Computers with Qiskit: Your First Quantum Program on a Real Quantum Processor

Learning Objectives

After completing this chapter, you will be able to:

• Understand the Qiskit ecosystem architecture (Terra, Aer, and the IBM Quantum platform)
• Set up an IBM Quantum account and configure your local environment
• Build, visualize, and simulate quantum circuits using Qiskit
• Choose between simulator backends and real quantum hardware
• Submit jobs to IBM Quantum processors and retrieve results
• Interpret measurement outcomes, error rates, and calibration data
• Visualize results using histograms and state tomography tools
• Navigate job queues and understand hardware execution models
• Use parameterized circuits for variational algorithms
• Apply error mitigation techniques to improve hardware results

8.1 The Qiskit Ecosystem

Qiskit (Quantum Information Science Kit) is IBM's open-source framework for quantum computing. It has evolved through several major versions; as of 2024, the architecture centers on Qiskit (the core SDK, formerly Terra) and Qiskit Runtime (the cloud execution layer). The legacy components — Aer (simulator), Ignis (error characterization), and Aqua (application algorithms) — have been integrated, replaced, or spun off into the broader ecosystem.

Historical Context. Qiskit was first released in 2017 by IBM Research, making it one of the earliest open-source quantum computing frameworks. It emerged from IBM's Quantum Experience, a cloud-based platform that allowed anyone to run programs on IBM's quantum processors. The name "Qiskit" reflects its goal: a "kit" of tools for "quantum information science." Over the years, Qiskit has undergone significant architectural changes, with the transition from Terra/Aer/Ignis/Aqua to the unified qiskit package being the most notable.

8.1.1 Core Components

What each component does:

• qiskit (Terra): The foundation. Provides the QuantumCircuit class, gate definitions, transpilation passes, and visualization tools. If you're building a circuit, you're using Terra.

• qiskit-aer: A high-performance simulator backend written in C++ with Python bindings. It supports multiple simulation methods: statevector (exact), density matrix (for noise), stabilizer (for Clifford circuits), and matrix product state (for low-entanglement circuits). Aer is essential for development and debugging.

• qiskit-ibm-runtime: The bridge to IBM Quantum hardware. Provides the Sampler and Estimator primitives for running circuits on real devices, and manages job submission, queuing, and result retrieval.

• qiskit-transpiler-service: An optional cloud-based transpilation service that uses AI-optimized routing and decomposition to produce shorter circuits than the local transpiler.

8.1.2 The Qiskit Programming Model

The workflow follows a consistent pattern:

1. Build — Construct quantum circuits using QuantumCircuit
2. Transpile — Optimize and map circuits to hardware constraints
3. Execute — Run on a simulator or real backend
4. Analyze — Process and visualize results

# The canonical Qiskit workflow
from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit.visualization import plot_histogram

# 1. Build
qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)
qc.measure_all()

# 2. Transpile (optional for simulators, essential for hardware)
# 3. Execute
simulator = AerSimulator()
job = simulator.run(qc, shots=1024)
result = job.result()

# 4. Analyze
counts = result.get_counts()
print(counts)

Recurring Theme. Quantum computing is linear algebra, not magic. The Qiskit workflow — Build, Transpile, Execute, Analyze — mirrors the mathematical process: define a unitary, decompose it into hardware-native operations, apply it to a state, and interpret the results. Each step has a mathematical counterpart.

8.1.3 The Qiskit Object Model

Understanding the key classes in Qiskit helps navigate the framework:

QuantumCircuit
├── QuantumRegister    (collection of qubits)
├── ClassicalRegister  (collection of clbits)
├── Instruction        (gate, measurement, barrier)
│   ├── Gate          (unitary operation)
│   ├── Measure       (projective measurement)
│   └── Barrier       (optimization fence)
└── DAGCircuit        (internal DAG representation for transpilation)

Backend
├── AerSimulator      (local simulator)
├── IBMBackend        (real hardware)
└── FakeBackend       (simulator mimicking real device properties)

Result
├── Counts            (measurement outcomes)
├── Statevector       (if saved)
└── Unitary           (if saved)

The QuantumCircuit is a container for instructions. When you add gates, you're appending Instruction objects to an ordered list. The transpiler converts this list into a DAGCircuit (directed acyclic graph) for optimization, then back to a QuantumCircuit for execution.

8.2 Setting Up Your Environment

8.2.1 Installation

pip install qiskit qiskit-aer qiskit-ibm-runtime matplotlib pylatexenc

The pylatexenc package enables LaTeX rendering in circuit diagrams. matplotlib is required for visualization.

Troubleshooting common installation issues:

• Version conflicts: Use pip install qiskit==1.0 qiskit-aer==0.14 qiskit-ibm-runtime==0.20 for a known-compatible set.
• Rust dependency: qiskit-aer requires a Rust compiler on some platforms. On macOS, install Xcode command line tools (xcode-select --install). On Linux, install build-essential and libopenblas-dev.
• Conda users: Use conda install -c conda-forge qiskit qiskit-aer for pre-compiled binaries.

8.2.2 IBM Quantum Account Setup

1. Create an account at https://quantum-computing.ibm.com
2. Navigate to your account dashboard and copy your API token
3. Save your credentials:

from qiskit_ibm_runtime import QiskitRuntimeService

# Save credentials to disk (one-time setup)
QiskitRuntimeService.save_account(
    channel="ibm_quantum",
    token="YOUR_API_TOKEN_HERE",
    overwrite=True
)

# Load saved credentials
service = QiskitRuntimeService()
print("Connected to IBM Quantum!")
print(service.active_account())

Common Misconception. "I need a quantum computer in my lab to run quantum programs." No! IBM Quantum provides free cloud access to real quantum processors. You only need a Python environment and an internet connection. The free tier gives you access to 5-127 qubit processors with monthly execution quotas.

8.2.3 Exploring Available Backends

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()

# List all backends
for backend in service.backends():
    status = backend.status()
    print(f"{backend.name:25s} | Qubits: {backend.num_qubits:3d} | "
          f"Pending: {status.pending_jobs:4d} | Operational: {status.operational}")

# Get details for a specific backend
backend = service.backend("ibm_brisbane")
print(f"\nBackend: {backend.name}")
print(f"Qubits: {backend.num_qubits}")
print(f"Coupling map: {backend.coupling_map}")
print(f"Basis gates: {backend.operation_names}")

Understanding backend properties:

Each backend provides detailed calibration data including:

• Gate error rates: The probability that a gate produces an incorrect result. Typical values: ~0.1% for single-qubit gates, ~1% for two-qubit gates.
• T1 times: The energy relaxation time. A qubit in $|1\rangle$ decays to $|0\rangle$ with time constant $T_1$.
• T2 times: The dephasing time. Superposition states lose phase coherence with time constant $T_2$.
• Readout error rates: The probability of measuring the wrong result.

# Detailed calibration data
props = backend.properties()

print("Gate error rates (selected):")
for gate in props.gates[:10]:  # Show first 10 for brevity
    for param in gate.parameters:
        if param.name == 'gate_error':
            qubits = ', '.join(str(q) for q in gate.qubits)
            print(f"  {gate.gate}({qubits}): {param.value:.4%}")

print(f"\nQubit T1 and T2 times:")
for i, qubit in enumerate(props.qubits[:5]):  # Show first 5
    t1 = next(p.value for p in qubit if p.name == 'T1')
    t2 = next(p.value for p in qubit if p.name == 'T2')
    print(f"  Qubit {i}: T1={t1:.1f} μs, T2={t2:.1f} μs")

8.3 Building Circuits with QuantumCircuit

8.3.1 The QuantumCircuit Class

QuantumCircuit is the central data structure. It stores a list of quantum and classical registers, and an ordered list of instructions (gates, measurements, barriers, and resets).

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister

# Method 1: Implicit registers
qc1 = QuantumCircuit(3, 3)  # 3 qubits, 3 classical bits

# Method 2: Explicit registers (preferred for complex circuits)
qr = QuantumRegister(3, 'q')
cr = ClassicalRegister(3, 'c')
qc2 = QuantumCircuit(qr, cr)

# Method 3: Empty circuit, add registers later
qc3 = QuantumCircuit()
qr3 = QuantumRegister(2, 'data')
cr3 = ClassicalRegister(2, 'meas')
qc3.add_register(qr3, cr3)

When to use each method:

• Method 1 (implicit): Quick prototyping, tutorials, simple circuits. The registers are named q and c by default.
• Method 2 (explicit): Production code, multi-register circuits, when you need named registers for clarity.
• Method 3 (dynamic): When you need to build circuits incrementally or compose sub-circuits.

8.3.2 Gate Methods

Qiskit provides methods for all standard gates:

qc = QuantumCircuit(3)

# Single-qubit gates
qc.h(0)           # Hadamard on qubit 0
qc.x(1)           # Pauli-X on qubit 1
qc.y(2)           # Pauli-Y on qubit 2
qc.z(0)           # Pauli-Z on qubit 0
qc.s(1)           # S gate
qc.t(2)           # T gate
qc.sdg(1)         # S-dagger
qc.tdg(2)         # T-dagger
qc.p(0.5, 0)      # Phase gate P(λ) with λ=0.5
qc.u(0.3, 0.2, 0.1, 1)  # U(θ, φ, λ)

# Two-qubit gates
qc.cx(0, 1)       # CNOT: control=0, target=1
qc.cz(1, 2)       # CZ: control=1, target=2
qc.cp(0.5, 0, 2)  # Controlled phase
qc.swap(0, 1)     # SWAP

# Three-qubit gates
qc.ccx(0, 1, 2)   # Toffoli (CCNOT)

# Special operations
qc.barrier()      # Prevent optimization across this point
qc.reset(0)       # Reset qubit 0 to |0⟩
qc.measure(0, 0)  # Measure qubit 0 → classical bit 0
qc.measure_all()  # Measure all qubits

The u gate in detail:

The $U(\theta, \phi, \lambda)$ gate is the most general single-qubit rotation:

$$U(\theta, \phi, \lambda) = \begin{pmatrix} \cos(\theta/2) & -e^{i\lambda}\sin(\theta/2) \\ e^{i\phi}\sin(\theta/2) & e^{i(\phi+\lambda)}\cos(\theta/2) \end{pmatrix}$$

All other single-qubit gates are special cases: - $H = U(\pi/2, 0, \pi)$ (up to global phase) - $X = U(\pi, 0, \pi)$ - $Y = U(\pi, 0, 0)$ (up to global phase) - $Z = U(0, 0, \pi)$ (up to global phase) - $S = U(0, 0, \pi/2)$ - $T = U(0, 0, \pi/4)$

8.3.3 Circuit Visualization

qc = QuantumCircuit(3)
qc.h(0)
qc.cx(0, 1)
qc.cx(1, 2)
qc.measure_all()

# Text-based drawing (works everywhere)
print(qc.draw('text'))

# Matplotlib drawing (requires matplotlib)
qc.draw('mpl')

# LaTeX drawing (requires pylatexenc)
qc.draw('latex')

Output (text):

     ┌───┐          ┌─┐
q_0: ┤ H ├──■───────┤M├
     └───┘┌─┴─┐     └╥┘
q_1: ─────┤ X ├──■───╫─
          └───┘┌─┴─┐ ║
q_2: ──────────┤ X ├─╫─
               └───┘ ║
c: 3/════════════════╩═
                      0

Customizing circuit diagrams:

# Reverse bit order (useful for matching physics convention)
print(qc.draw('text', reverse_bits=True))

# Fold long circuits
long_qc = QuantumCircuit(5)
for i in range(20):
    long_qc.h(i % 5)
print(long_qc.draw('text', fold=40))

# Style customization
style = {'displaycolor': {'h': '#FF0000', 'cx': '#0000FF'}}
qc.draw('mpl', style=style)

8.3.4 Circuit Properties and Inspection

qc = QuantumCircuit(3)
qc.h(0)
qc.cx(0, 1)
qc.cx(1, 2)

# Circuit metrics
print(f"Width (qubits): {qc.num_qubits}")
print(f"Depth: {qc.depth()}")
print(f"Gate count: {qc.size()}")
print(f"Gate count (non-local): {qc.num_nonlocal_gates()}")

# Detailed gate breakdown
print(f"Operations: {qc.count_ops()}")

# Check if circuit is unitary (no measurements or resets)
print(f"Is unitary: {qc.num_clbits == 0}")

# Get the unitary matrix
from qiskit.quantum_info import Operator
U = Operator(qc)
print(f"Unitary shape: {U.data.shape}")

Try it yourself: Create a 4-qubit GHZ state circuit and report its depth, gate count, and number of non-local (two-qubit) gates. Then transpile it for FakeBrisbane() and compare the metrics.

8.4 Simulators: The Safe Playground

8.4.1 Types of Simulators

Qiskit Aer provides several simulation methods:

Choosing the right simulator:

• Statevector: Use when you need the full quantum state. Limited to ~30 qubits (memory: $2^n \times 16$ bytes). Best for debugging and verifying circuits.

• QASM (shots-based): Use when you need measurement statistics. Can handle circuits with measurements and mid-circuit operations. The default mode.

• Density matrix: Use when modeling noise. Represents the mixed state $\rho$, requiring $2^{2n}$ memory. Limited to ~15 qubits.

• Matrix Product State (MPS): Use for circuits with low entanglement (e.g., shallow circuits, local Hamiltonians). Can handle up to ~100 qubits if entanglement is low.

• Stabilizer: Use for Clifford circuits (circuits with only H, S, CNOT, and measurements). Simulates in polynomial time. Limited to Clifford operations.

8.4.2 Statevector Simulation

from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit.quantum_info import Statevector
import numpy as np

qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)

# Method 1: Direct statevector computation
state = Statevector.from_instruction(qc)
print("Statevector:")
print(state.draw('latex'))

# Verify: should be (|00⟩ + |11⟩)/√2
print(f"\nProbabilities: {state.probabilities_dict()}")
print(f"Fidelity with |Φ+⟩: {state.fidelity(Statevector.from_label('phi+')):.6f}")

# Method 2: Via AerSimulator
sim = AerSimulator(method='statevector')
qc_sv = qc.copy()
qc_sv.save_statevector()
result = sim.run(qc_sv).result()
statevector = result.get_statevector()
print(f"\nStatevector: {statevector}")

# Method 3: Compute individual qubit states
from qiskit.quantum_info import partial_trace, DensityMatrix
rho = DensityMatrix(state)
rho_q0 = partial_trace(rho, [1])  # Trace out qubit 1
rho_q1 = partial_trace(rho, [0])  # Trace out qubit 0

print(f"\nQubit 0 reduced state:\n{rho_q0.data}")
print(f"Qubit 1 reduced state:\n{rho_q1.data}")

Worked Example 8.1: Computing the Bell state

Starting from $|00\rangle$:

$$|00\rangle \xrightarrow{H \otimes I} \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle) \otimes |0\rangle = \frac{1}{\sqrt{2}}(|00\rangle + |10\rangle)$$

$$\xrightarrow{\text{CNOT}_{0\to1}} \frac{1}{\sqrt{2}}(|00\rangle + |11\rangle) = |\Phi^+\rangle$$

The statevector is $\frac{1}{\sqrt{2}}(1, 0, 0, 1)^T$, with probabilities $\{00: 0.5, 11: 0.5\}$.

8.4.3 Sampling with Shots

from qiskit_aer import AerSimulator

qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)
qc.measure_all()

sim = AerSimulator(method='automatic')
job = sim.run(qc, shots=8192)
result = job.result()
counts = result.get_counts()

print("Measurement counts:")
for state, count in sorted(counts.items()):
    print(f"  |{state}⟩: {count:5d} ({count/8192*100:.1f}%)")

Expected output (approximately):

  |00⟩:  4096 (50.0%)
  |11⟩:  4096 (50.0%)

Understanding shot statistics: The results follow a binomial distribution. For a Bell state measured 8192 times:

• Expected count for $|00\rangle$: $\mu = 8192 \times 0.5 = 4096$
• Standard deviation: $\sigma = \sqrt{8192 \times 0.5 \times 0.5} \approx 64$
• 99.7% confidence interval: $4096 \pm 192$ (i.e., $3904 < \text{count} < 4288$)

# Statistical analysis of shots
import numpy as np

n_shots = 8192
p_00 = 0.5  # probability of |00⟩
std_dev = np.sqrt(n_shots * p_00 * (1 - p_00))
print(f"Expected: {n_shots * p_00:.0f} ± {3*std_dev:.0f} (3σ interval)")

# Run multiple times to see statistical variation
results = []
for _ in range(20):
    job = sim.run(qc, shots=8192)
    count_00 = job.result().get_counts().get('00', 0)
    results.append(count_00)

print(f"Mean count for |00⟩: {np.mean(results):.0f} ± {np.std(results):.0f}")
print(f"Expected: {n_shots * p_00:.0f} ± {std_dev:.0f}")

8.4.4 Noisy Simulation

from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel, depolarizing_error, thermal_relaxation_error

# Build a noise model with multiple error sources
noise_model = NoiseModel()
error_1q = depolarizing_error(0.001, 1)  # 0.1% single-qubit error
error_2q = depolarizing_error(0.01, 2)   # 1% two-qubit error

noise_model.add_all_qubit_quantum_error(error_1q, ['u', 'u2', 'u3', 'h', 'x', 's', 't'])
noise_model.add_all_qubit_quantum_error(error_2q, ['cx'])

# Run with noise
sim_noisy = AerSimulator(noise_model=noise_model)
job = sim_noisy.run(qc, shots=8192)
counts_noisy = job.result().get_counts()

print("Noisy counts:")
for state, count in sorted(counts_noisy.items()):
    print(f"  |{state}⟩: {count:5d} ({count/8192*100:.1f}%)")

Notice the appearance of $|01\rangle$ and $|10\rangle$ states — these are errors introduced by the noise model.

What is depolarizing error? The depolarizing channel with parameter $\lambda$ replaces the quantum state with the maximally mixed state with probability $\lambda$:

$$\mathcal{E}(\rho) = (1 - \lambda)\rho + \lambda \frac{I}{d}$$

For a single qubit ($d=2$), this is equivalent to:

$$\mathcal{E}(\rho) = (1 - \lambda)\rho + \frac{\lambda}{3}(X\rho X + Y\rho Y + Z\rho Z)$$

With probability $1 - \lambda$, the state is unchanged; with probability $\lambda/3$, a random Pauli error ($X$, $Y$, or $Z$) occurs.

Worked Example 8.2: Predicting error rates

For a Bell state circuit ($H$ + CNOT) with single-qubit error rate $\epsilon_1 = 0.001$ and two-qubit error rate $\epsilon_2 = 0.01$:

• Probability of no error: $(1 - \epsilon_1)^2 \times (1 - \epsilon_2) \approx (0.999)^2 \times 0.99 \approx 0.988$
• Probability of at least one error: $\approx 0.012$, or about 1.2%

Of the erroneous outcomes, the noise is roughly symmetric, so we'd expect $|01\rangle$ and $|10\rangle$ to each appear at about 0.6%.

8.4.5 Simulating from a Real Backend's Noise Model

from qiskit_aer import AerSimulator
from qiskit.providers.fake_provider import FakeBrisbane

# Load the noise model from a real backend
real_backend = FakeBrisbane()
noise_model = NoiseModel.from_backend(real_backend)
coupling_map = real_backend.configuration().coupling_map
basis_gates = noise_model.basis_gates

# Create a simulator with the real backend's noise model
sim_noisy = AerSimulator(noise_model=noise_model,
                         coupling_map=coupling_map,
                         basis_gates=basis_gates)

# Run the Bell state circuit
qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)
qc.measure_all()

result = sim_noisy.run(qc, shots=8192).result()
counts = result.get_counts()
print("Simulated hardware results:")
for state, count in sorted(counts.items()):
    print(f"  |{state}⟩: {count:5d} ({count/8192*100:.1f}%)")

This gives you a realistic preview of what you'd see on the actual device, without spending your quantum computing quota.

Recurring Theme. Noise is the enemy. Every simulation that includes noise shows degraded results compared to the ideal case. Understanding and quantifying this degradation is essential for designing algorithms that work on real hardware.

8.5 Running on Real Quantum Hardware

8.5.1 The Sampler Primitive

Qiskit Runtime introduces primitives — high-level interfaces for common quantum computing tasks. The Sampler primitive executes circuits and returns quasi-probability distributions.

from qiskit_ibm_runtime import QiskitRuntimeService, Sampler
from qiskit import QuantumCircuit, transpile

service = QiskitRuntimeService()
backend = service.backend("ibm_brisbane")

# Build a simple circuit
qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)
qc.measure_all()

# Transpile for the target backend
qc_transpiled = transpile(qc, backend=backend, optimization_level=3)

# Run using the Sampler primitive
sampler = Sampler(backend=backend)
job = sampler.run([qc_transpiled], shots=4096)

print(f"Job ID: {job.job_id()}")
print("Waiting for results...")

result = job.result()
quasi_dists = result.quasi_dists[0]

print("\nQuasi-probability distribution:")
for state, prob in sorted(quasi_dists.items()):
    print(f"  |{state:02b}⟩: {prob:.4f}")

What are quasi-probabilities? The Sampler returns a quasi-probability distribution, which is like a probability distribution but may contain negative or complex values (after error mitigation). To get standard probabilities, you can call .nearest_probability_distribution() on the result.

8.5.2 The Estimator Primitive

The Estimator primitive computes expectation values of observables:

from qiskit_ibm_runtime import QiskitRuntimeService, Estimator
from qiskit.quantum_info import SparsePauliOp
from qiskit import QuantumCircuit, transpile

service = QiskitRuntimeService()
backend = service.backend("ibm_brisbane")

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

# Observables: Z⊗Z, X⊗X, Y⊗Y
observable_zz = SparsePauliOp("ZZ")
observable_xx = SparsePauliOp("XX")
observable_yy = SparsePauliOp("YY")

qc_transpiled = transpile(qc, backend=backend)

estimator = Estimator(backend=backend)
job = estimator.run(
    [(qc_transpiled, observable_zz),
     (qc_transpiled, observable_xx),
     (qc_transpiled, observable_yy)],
    shots=4096
)

result = job.result()
print(f"⟨Z⊗Z⟩ = {result.values[0]:.4f}")  # Expected: 1.0
print(f"⟨X⊗X⟩ = {result.values[1]:.4f}")  # Expected: 1.0
print(f"⟨Y⊗Y⟩ = {result.values[2]:.4f}")  # Expected: -1.0

Why the Bell state has these expectation values:

The Bell state $|\Phi^+\rangle = \frac{1}{\sqrt{2}}(|00\rangle + |11\rangle)$ has:

$$\langle \Phi^+ | ZZ | \Phi^+ \rangle = \frac{1}{2}(1 + 1) = 1$$

$$\langle \Phi^+ | XX | \Phi^+ \rangle = \frac{1}{2}(1 + 1) = 1$$

$$\langle \Phi^+ | YY | \Phi^+ \rangle = \frac{1}{2}(-1 + (-1)) = -1$$

This is because $|\Phi^+\rangle$ is a +1 eigenstate of both $Z \otimes Z$ and $X \otimes X$, and a -1 eigenstate of $Y \otimes Y$.

8.5.3 Understanding Job Lifecycle

A job on IBM Quantum goes through several stages:

INITIALIZING → QUEUED → VALIDATING → RUNNING → COMPLETED
                                               → CANCELLED
                                               → ERROR

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()

# Retrieve a previous job
job = service.job("JOB_ID_HERE")
print(f"Status: {job.status()}")
print(f"Creation date: {job.creation_date()}")
print(f"Execution time: {job.metrics().get('execution_time', 'N/A')}")

# Monitor a running job
import time
while job.status() not in ['DONE', 'CANCELLED', 'ERROR']:
    print(f"Status: {job.status()}")
    time.sleep(10)
print(f"Final status: {job.status()}")

Understanding queue times: IBM Quantum uses a fair-share scheduling system. Your position in the queue depends on: 1. Your priority level (depends on your account tier) 2. The number of jobs you've submitted recently 3. The current load on the backend

Common Misconception. "I should always use the largest processor." No! Larger processors have longer queues and more qubits than you probably need. For a 2-qubit Bell state, use the smallest available processor. Only use large processors when your circuit requires many qubits.

8.5.4 Reading Error Rates and Calibration Data

backend = service.backend("ibm_brisbane")
props = backend.properties()

print("Gate error rates:")
for gate in props.gates:
    for param in gate.parameters:
        if param.name == 'gate_error':
            qubits = ', '.join(str(q) for q in gate.qubits)
            print(f"  {gate.gate}({qubits}): {param.value:.4%}")

print(f"\nQubit T1 times (μs):")
for qubit in props.qubits:
    t1 = next(p.value for p in qubit if p.name == 'T1')
    print(f"  Qubit {qubit.qubit}: {t1:.1f}")

print(f"\nQubit T2 times (μs):")
for qubit in props.qubits:
    t2 = next(p.value for p in qubit if p.name == 'T2')
    print(f"  Qubit {qubit.qubit}: {t2:.1f}")

How to interpret calibration data:

• T1 (energy relaxation time): How long a qubit in $|1\rangle$ takes to decay to $|0\rangle$. Longer is better. Typical values: 50-500 μs.
• T2 (dephasing time): How long a qubit maintains phase coherence. Longer is better. Typically T2 < T1.
• Gate error rates: Probability that a gate produces an incorrect result. Typical values: 0.01-0.1% for single-qubit gates, 0.5-5% for two-qubit gates.
• Readout error rates: Probability of measuring the wrong result. Typical values: 1-10%.

Try it yourself: Connect to IBM Quantum and identify the "best" qubit (lowest error rates) and the "best" two-qubit connection (lowest CNOT error). Use transpile with initial_layout to place your circuit on these qubits.

8.6 Visualizing Results

8.6.1 Histograms

from qiskit.visualization import plot_histogram
import matplotlib.pyplot as plt

# Simulated results
counts = {'00': 4123, '01': 45, '10': 38, '11': 3986}

plot_histogram(counts, title="Bell State Measurement Results",
               figsize=(8, 5), color='midnightblue')
plt.show()

# Comparing ideal vs noisy results
counts_ideal = {'00': 4096, '11': 4096}
counts_noisy = {'00': 3850, '01': 120, '10': 95, '11': 3931}

plot_histogram([counts_ideal, counts_noisy],
               legend=['Ideal', 'Noisy'],
               title="Bell State: Ideal vs Noisy",
               figsize=(10, 5))
plt.show()

8.6.2 State City (Density Matrix Visualization)

from qiskit.quantum_info import DensityMatrix
from qiskit.visualization import plot_state_city
import matplotlib.pyplot as plt

qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)

rho = DensityMatrix.from_instruction(qc)
plot_state_city(rho, title="Bell State Density Matrix")
plt.show()

# Compare with a mixed state
from qiskit.quantum_info import partial_trace
rho_mixed = 0.5 * DensityMatrix.from_label('00') + 0.5 * DensityMatrix.from_label('11')
plot_state_city(rho_mixed, title="Mixed State (Classical Mixture)")
plt.show()

Understanding the state city plot: The real and imaginary parts of the density matrix are displayed as 3D bar charts. For a pure Bell state, you'll see two bars (corresponding to the off-diagonal elements $\rho_{00,11}$ and $\rho_{11,00}$) in addition to the diagonal elements. For a mixed state, the off-diagonal elements are zero.

8.6.3 Bloch Sphere

from qiskit.visualization import plot_bloch_multivector
from qiskit.quantum_info import Statevector

# Single qubit states
qc = QuantumCircuit(1)
qc.h(0)
qc.s(0)

state = Statevector.from_instruction(qc)
plot_bloch_multivector(state)
plt.show()

# Compare multiple states on the Bloch sphere
from qiskit.visualization import plot_bloch_vector
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 4, figsize=(16, 4))
states = ['z', 'x', 'y', 'h']
labels = ['|0⟩', '|+⟩', '|+i⟩', '|-i⟩']
vectors = [(0, 0, 1), (1, 0, 0), (0, 1, 0), (0, -1, 0)]

for i, (vec, label) in enumerate(zip(vectors, labels)):
    plot_bloch_vector(vec, title=label, ax=axes[i])
plt.tight_layout()
plt.show()

8.6.4 Gate Map and Error Map

from qiskit.visualization import plot_gate_map, plot_error_map
from qiskit.providers.fake_provider import FakeBrisbane

backend = FakeBrisbane()

# Visualize the qubit connectivity
plot_gate_map(backend)
plt.show()

# Visualize error rates on the device
plot_error_map(backend)
plt.show()

The gate map shows which qubits are connected (where CNOT can be applied directly). The error map overlays error rates on each qubit and connection, helping you identify the best qubits for your circuit.

8.7 Complete Walkthrough: From Notebook to Quantum Processor

Let's put everything together in a complete, production-ready example:

"""
Complete Qiskit Workflow: Bell State on IBM Quantum Hardware
============================================================
This script demonstrates the full lifecycle:
Build → Transpile → Execute → Analyze
"""

from qiskit import QuantumCircuit, transpile
from qiskit_ibm_runtime import QiskitRuntimeService, Sampler
from qiskit.visualization import plot_histogram
from qiskit_aer import AerSimulator
from qiskit.quantum_info import Statevector, fidelity
import matplotlib.pyplot as plt
import numpy as np

# ─── Configuration ───────────────────────────────────────────
USE_HARDWARE = False  # Set to True to run on real hardware
BACKEND_NAME = "ibm_brisbane"
SHOTS = 4096

# ─── Step 1: Build the Circuit ───────────────────────────────
def create_bell_circuit() -> QuantumCircuit:
    """Create a circuit that prepares the Bell state |Φ⁺⟩."""
    qc = QuantumCircuit(2, 2)
    qc.h(0)
    qc.cx(0, 1)
    qc.barrier()
    qc.measure([0, 1], [0, 1])
    return qc

qc = create_bell_circuit()
print("Circuit:")
print(qc.draw('text'))
print(f"Depth: {qc.depth()}, Width: {qc.num_qubits}, Gates: {qc.size()}")

# ─── Step 2: Simulate Locally ────────────────────────────────
print("\n─── Local Simulation ───")
sim = AerSimulator()
job_sim = sim.run(qc, shots=SHOTS)
counts_sim = job_sim.result().get_counts()

print("Simulated counts:")
for state, count in sorted(counts_sim.items()):
    print(f"  |{state}⟩: {count:5d} ({count/SHOTS*100:.1f}%)")

# ─── Step 3: Noisy Simulation ────────────────────────────────
print("\n─── Noisy Simulation ───")
from qiskit_aer.noise import NoiseModel
from qiskit.providers.fake_provider import FakeBrisbane

real_backend = FakeBrisbane()
noise_model = NoiseModel.from_backend(real_backend)
sim_noisy = AerSimulator(noise_model=noise_model)

counts_noisy = sim_noisy.run(qc, shots=SHOTS).result().get_counts()
print("Noisy simulated counts:")
for state, count in sorted(counts_noisy.items()):
    print(f"  |{state}⟩: {count:5d} ({count/SHOTS*100:.1f}%)")

# ─── Step 4: Compute Fidelity ───────────────────────────────
ideal_probs = {'00': 0.5, '11': 0.5}

def compute_fidelity(counts, shots, ideal_probs):
    """Compute classical fidelity between measurement distribution and ideal."""
    total = sum(counts.values())
    fidelity = 0
    for state in ideal_probs:
        measured_prob = counts.get(state, 0) / total
        fidelity += np.sqrt(ideal_probs[state] * measured_prob)
    return fidelity ** 2

f_sim = compute_fidelity(counts_sim, SHOTS, ideal_probs)
f_noisy = compute_fidelity(counts_noisy, SHOTS, ideal_probs)

print(f"\nFidelity (ideal sim): {f_sim:.4f}")
print(f"Fidelity (noisy sim): {f_noisy:.4f}")

# ─── Step 5: Run on Hardware (if enabled) ────────────────────
if USE_HARDWARE:
    print(f"\n─── Running on {BACKEND_NAME} ───")

    service = QiskitRuntimeService()
    backend = service.backend(BACKEND_NAME)

    # Check backend status
    status = backend.status()
    print(f"Backend status: {'Operational' if status.operational else 'Down'}")
    print(f"Pending jobs: {status.pending_jobs}")

    # Transpile
    qc_transpiled = transpile(
        qc, backend=backend,
        optimization_level=3,
        seed_transpiler=42
    )
    print(f"\nTranspiled depth: {qc_transpiled.depth()}")
    print(f"Transpiled gates: {qc_transpiled.size()}")

    # Execute
    sampler = Sampler(backend=backend)
    job = sampler.run([qc_transpiled], shots=SHOTS)
    print(f"Job submitted: {job.job_id()}")

    # Wait for results
    result = job.result()
    quasi_dist = result.quasi_dists[0]

    # Convert quasi-distribution to counts-like format
    counts_hw = {f"{k:02b}": int(v * SHOTS) for k, v in quasi_dist.items()}
    print("\nHardware results:")
    for state, count in sorted(counts_hw.items()):
        print(f"  |{state}⟩: {count:5d}")

    # ─── Step 6: Compare and Analyze ─────────────────────────
    f_hw = compute_fidelity(counts_hw, SHOTS, ideal_probs)
    print(f"\nFidelity (hardware): {f_hw:.4f}")

    # Plot comparison
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    plot_histogram(counts_sim, ax=axes[0], title="Ideal Simulator")
    plot_histogram(counts_noisy, ax=axes[1], title="Noisy Simulator")
    plot_histogram(counts_hw, ax=axes[2], title=BACKEND_NAME)
    plt.suptitle("Bell State: Ideal vs Noisy vs Hardware")
    plt.tight_layout()
    plt.savefig("bell_state_comparison.png", dpi=150)
    print("\nPlot saved to bell_state_comparison.png")

8.8 Advanced Circuit Construction

8.8.1 Parameterized Circuits

Parameterized circuits are essential for variational algorithms (VQE, QAOA, etc.) where circuit parameters are optimized classically:

from qiskit.circuit import Parameter, ParameterVector
import numpy as np

# Single parameter
theta = Parameter('θ')
qc = QuantumCircuit(1)
qc.ry(theta, 0)

# Bind and execute
qc_bound = qc.assign_parameters({theta: np.pi/2})
print(qc_bound.draw('text'))

# Parameter vector
params = ParameterVector('θ', 4)
qc_variational = QuantumCircuit(2)
qc_variational.ry(params[0], 0)
qc_variational.ry(params[1], 1)
qc_variational.cx(0, 1)
qc_variational.ry(params[2], 0)
qc_variational.ry(params[3], 1)

# Bind all parameters
bound_circuit = qc_variational.assign_parameters(
    {p: np.random.random() * 2 * np.pi for p in params}
)
print(f"\nParameterized circuit depth: {qc_variational.depth()}")
print(f"Bound circuit depth: {bound_circuit.depth()}")

Worked Example 8.3: Variational quantum eigensolver setup

The parameterized circuit pattern is the core of VQE:

from qiskit.circuit import ParameterVector
from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp, Statevector
import numpy as np
from scipy.optimize import minimize

def create_ansatz(n_qubits, n_layers):
    """Create a hardware-efficient ansatz with parameterized rotations."""
    params = ParameterVector('θ', n_qubits * (2 * n_layers + 1))
    qc = QuantumCircuit(n_qubits)
    param_idx = 0

    # Initial rotations
    for i in range(n_qubits):
        qc.ry(params[param_idx], i)
        param_idx += 1

    for layer in range(n_layers):
        # Entanglement
        for i in range(n_qubits - 1):
            qc.cx(i, i + 1)

        # Rotations
        for i in range(n_qubits):
            qc.ry(params[param_idx], i)
            param_idx += 1
            qc.rz(params[param_idx], i)
            param_idx += 1

    return qc, params

# Example: minimize ⟨ψ(θ)|ZZ|ψ(θ)⟩ for 2 qubits
ansatz, params = create_ansatz(2, 1)

def objective(theta_values):
    """Compute ⟨ψ(θ)|ZZ|ψ(θ)⟩."""
    bound = ansatz.assign_parameters(dict(zip(params, theta_values)))
    state = Statevector.from_instruction(bound)
    obs = SparsePauliOp("ZZ")
    return state.expectation_value(obs).real

# Optimize
initial = np.random.random(len(params)) * 2 * np.pi
result = minimize(objective, initial, method='COBYLA')
print(f"Minimum eigenvalue: {result.fun:.6f}")
print(f"Expected (ZZ ground state): -1.0")

8.8.2 Composite Circuits

# Build sub-circuits
sub_qc1 = QuantumCircuit(2)
sub_qc1.h(0)
sub_qc1.cx(0, 1)

sub_qc2 = QuantumCircuit(2)
sub_qc2.x(0)
sub_qc2.z(1)

# Compose into a larger circuit
main_qc = QuantumCircuit(4)
main_qc.compose(sub_qc1, qubits=[0, 1], inplace=True)
main_qc.compose(sub_qc2, qubits=[2, 3], inplace=True)
main_qc.barrier()
main_qc.cx(1, 2)
main_qc.measure_all()

print(main_qc.draw('text'))

The compose method is powerful for building modular circuits:

• qc.compose(other, qubits=[...], clbits=[...], inplace=True) — appends other to qc
• qc.append(gate, qargs=[...]) — appends a single instruction
• qc += other — shorthand for compose with inplace=True

8.8.3 Custom Gates

from qiskit import QuantumCircuit
from qiskit.circuit.library import UnitaryGate
import numpy as np

# Define a custom unitary
theta = np.pi / 5
custom_matrix = np.array([
    [np.cos(theta), -np.sin(theta)],
    [np.sin(theta),  np.cos(theta)]
])

# Create a custom gate
custom_gate = UnitaryGate(custom_matrix, label='R(π/5)')

# Use in a circuit
qc = QuantumCircuit(2)
qc.append(custom_gate, [0])
qc.cx(0, 1)
qc.append(custom_gate, [1])

print(qc.draw('text'))

# Verify it's unitary
print(f"\nIs unitary? {np.allclose(custom_matrix @ custom_matrix.conj().T, np.eye(2))}")

8.8.4 Using the Circuit Library

Qiskit's circuit library provides pre-built building blocks:

from qiskit.circuit.library import (
    QFT, GroverOperator, PhaseEstimation,
    RealAmplitudes, EfficientSU2, TwoLocal,
    ZZFeatureMap, PauliFeatureMap
)

# Quantum Fourier Transform
qft = QFT(num_qubits=4, do_swaps=True)
print("QFT circuit:")
print(qft.decompose().draw('text', fold=80))

# Variational ansatz
ansatz = EfficientSU2(num_qubits=4, su2_gates=['ry', 'rz'],
                      entanglement='circular', reps=2)
print(f"\nEfficientSU2 parameters: {ansatz.num_parameters}")

# Feature map for quantum machine learning
feature_map = ZZFeatureMap(feature_dimension=4, reps=2)
print(f"ZZFeatureMap depth: {feature_map.decompose().depth()}")

8.9 Error Mitigation Techniques

Real hardware introduces errors. Qiskit Runtime provides error mitigation:

8.9.1 Types of Errors

8.9.2 Error Mitigation Strategies

from qiskit_ibm_runtime import QiskitRuntimeService, Estimator, Options

service = QiskitRuntimeService()
backend = service.backend("ibm_brisbane")

# Configure error mitigation
options = Options()
options.resilience_level = 2  # 0=none, 1=readout, 2=gate+readout, 3=advanced
options.optimization_level = 3

estimator = Estimator(backend=backend, options=options)

# Run with error mitigation
qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)

from qiskit.quantum_info import SparsePauliOp
observable = SparsePauliOp("ZZ")

job = estimator.run([(qc, observable)], shots=4096)
result = job.result()
print(f"⟨Z⊗Z⟩ with error mitigation: {result.values[0]:.4f}")

Resilience levels explained:

How zero-noise extrapolation works:

1. Run the circuit at the original noise level
2. Run the circuit at artificially increased noise levels (by repeating gate sequences)
3. Fit a curve through the results as a function of noise level
4. Extrapolate to zero noise

# Manual zero-noise extrapolation example
from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel, depolarizing_error
import numpy as np

def run_with_noise(p_error, shots=4096):
    """Run Bell state circuit with specified depolarizing error."""
    noise_model = NoiseModel()
    noise_model.add_all_qubit_quantum_error(
        depolarizing_error(p_error, 1), ['h'])
    noise_model.add_all_qubit_quantum_error(
        depolarizing_error(p_error * 10, 2), ['cx'])

    qc = QuantumCircuit(2)
    qc.h(0)
    qc.cx(0, 1)
    qc.measure_all()

    sim = AerSimulator(noise_model=noise_model)
    counts = sim.run(qc, shots=shots).result().get_counts()
    return counts.get('00', 0) + counts.get('11', 0)

# Run at different noise levels
noise_levels = [0.001, 0.002, 0.004]
p_ideal = []
for noise in noise_levels:
    p = run_with_noise(noise) / 4096
    p_ideal.append(p)
    print(f"Noise {noise:.4f}: fidelity = {p:.4f}")

# Extrapolate to zero noise (Richardson extrapolation)
# For 3 noise levels λ₁ < λ₂ < λ₃:
# Extrapolated = (λ₂λ₃*p₁ - λ₁λ₃*p₂ + λ₁λ₂*p₃) / ((λ₂-λ₁)(λ₃-λ₂)) ... simplified
# Linear extrapolation:
from numpy.polynomial import polynomial as P
coeffs = np.polyfit(noise_levels, p_ideal, 2)
extrapolated = np.polyval(coeffs, 0)
print(f"\nExtrapolated (zero-noise) fidelity: {extrapolated:.4f}")
print(f"Expected: 1.0000")

8.9.3 Readout Error Mitigation

from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel, ReadoutError
import numpy as np

# Create a readout error model
# Confusion matrix: p(measure i | prepared j)
p0_given_0 = 0.98  # Correctly measure |0⟩ when in |0⟩
p1_given_0 = 0.02   # Incorrectly measure |1⟩ when in |0⟩
p0_given_1 = 0.03   # Incorrectly measure |0⟩ when in |1⟩
p1_given_1 = 0.97   # Correctly measure |1⟩ when in |1⟩

readout_error = ReadoutError([
    [p0_given_0, p1_given_0],
    [p0_given_1, p1_given_1]
])

noise_model = NoiseModel()
noise_model.add_all_qubit_readout_error(readout_error)

# Run Bell state with readout error
qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)
qc.measure_all()

sim = AerSimulator(noise_model=noise_model)
counts = sim.run(qc, shots=8192).result().get_counts()

print("With readout error:")
for state, count in sorted(counts.items()):
    print(f"  |{state}⟩: {count:5d} ({count/8192*100:.1f}%)")

# Simple readout mitigation: measure the confusion matrix and invert
confusion_matrix = np.array([
    [p0_given_0, p0_given_1],
    [p1_given_0, p1_given_1]
])

# Raw measurement vector (for 2 qubits)
raw = np.array([counts.get('00', 0), counts.get('01', 0),
                counts.get('10', 0), counts.get('11', 0)]) / 8192

# For 2 qubits, the full confusion matrix is the tensor product
full_confusion = np.kron(confusion_matrix, confusion_matrix)

# Mitigate by inverting the confusion matrix
mitigated = np.linalg.solve(full_confusion, raw)
mitigated = np.clip(mitigated, 0, None)
mitigated /= mitigated.sum()

print("\nMitigated probabilities:")
for i, label in enumerate(['00', '01', '10', '11']):
    print(f"  |{label}⟩: {mitigated[i]:.4f}")

8.10 Qiskit Transpiler Deep Dive

8.10.1 What the Transpiler Does

The transpiler takes a logical circuit and produces a physical circuit that respects hardware constraints:

1. Unrolling: Decompose custom gates into basis gates
2. Routing: Map logical qubits to physical qubits satisfying connectivity
3. SWAP insertion: Add SWAPs to enable non-adjacent two-qubit gates
4. Optimization: Cancel redundant gates, merge adjacent rotations

from qiskit import QuantumCircuit, transpile
from qiskit.providers.fake_provider import FakeBrisbane

backend = FakeBrisbane()

# A circuit that requires SWAPs
qc = QuantumCircuit(5)
qc.h(0)
qc.cx(0, 4)  # qubit 0 and 4 might not be adjacent!
qc.cx(1, 3)
qc.measure_all()

# Transpile at different optimization levels
for level in range(4):
    qc_t = transpile(qc, backend, optimization_level=level,
                      seed_transpiler=42)
    print(f"Level {level}: depth={qc_t.depth():3d}, "
          f"gates={qc_t.size():3d}, "
          f"swaps={qc_t.count_ops().get('swap', 0)}")

8.10.2 Custom Transpiler Passes

from qiskit.transpiler import PassManager, TransformationPass
from qiskit.dagcircuit import DAGCircuit
from qiskit.circuit.library import RZGate

class MergeConsecutiveRZ(TransformationPass):
    """Merge consecutive RZ gates on the same qubit."""

    def run(self, dag: DAGCircuit) -> DAGCircuit:
        for node in dag.op_nodes():
            if node.op.name == 'rz':
                # Find the next operation on the same qubit
                successors = list(dag.successors(node))
                for succ in successors:
                    if (hasattr(succ, 'op') and
                        succ.op.name == 'rz' and
                        succ.qargs == node.qargs):
                        # Merge: RZ(a) * RZ(b) = RZ(a+b)
                        new_angle = node.op.params[0] + succ.op.params[0]
                        new_gate = RZGate(new_angle)
                        dag.substitute_node(node, new_gate, propagate=False)
                        dag.remove_op_node(succ)
                        break
        return dag

# Use the custom pass
from qiskit.transpiler.passes import Optimize1qGatesDecomposition

pm = PassManager([MergeConsecutiveRZ(), Optimize1qGatesDecomposition()])

qc = QuantumCircuit(2)
qc.rz(0.3, 0)
qc.rz(0.5, 0)
qc.cx(0, 1)

qc_optimized = pm.run(qc)
print("Original:")
print(qc.draw('text'))
print("\nOptimized:")
print(qc_optimized.draw('text'))

8.11 Debugging Quantum Circuits

8.11.1 Common Bugs and How to Find Them

Bug 1: Bit order confusion

Qiskit uses little-endian ordering: qc.measure(0, 0) maps qubit 0 to classical bit 0, and the bit string is read right-to-left. The state |01⟩ means qubit 0 is in |1⟩ and qubit 1 is in |0⟩.

qc = QuantumCircuit(2, 2)
qc.x(0)  # Put qubit 0 in |1⟩
qc.measure([0, 1], [0, 1])

result = AerSimulator().run(qc, shots=100).result()
counts = result.get_counts()
print(counts)
# Output: {'01': 100}  ← qubit 0 is the rightmost bit!

Common Misconception. "The bit string 01 means qubit 0 is in state $|0\rangle$ and qubit 1 is in state $|1\rangle$." No! In Qiskit's little-endian convention, 01 means qubit 0 is $|1\rangle$ and qubit 1 is $|0\rangle$. The rightmost bit corresponds to qubit 0.

Bug 2: Forgetting measurement

A circuit without measurement won't produce counts:

qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)
# No measurement! This will fail:
# result = sim.run(qc, shots=100).result().get_counts()  # Error!

Fix: Add qc.measure_all() before running on a shot-based simulator.

Bug 3: Global phase matters for controlled operations

Two unitaries that differ only by global phase are equivalent as standalone operations, but they differ when used as controlled operations:

from qiskit.quantum_info import Operator
import numpy as np

# X gate
qc1 = QuantumCircuit(1)
qc1.x(0)

# X = -iY * i (same effect, different global phase)
qc2 = QuantumCircuit(1)
qc2.sdg(0)
qc2.y(0)
qc2.s(0)

# These are equivalent as standalone operations
print(f"Same unitary (up to phase)? "
      f"{np.allclose(np.abs(np.diag(Operator(qc1).data.conj().T @ Operator(qc2).data)), 1.0)}")

# But controlled-X ≠ controlled-S†YS
qc1_ctrl = QuantumCircuit(2)
qc1_ctrl.append(qc1.to_gate().control(), [0, 1])

qc2_ctrl = QuantumCircuit(2)
qc2_ctrl.append(qc2.to_gate().control(), [0, 1])

print(f"Same controlled unitary? "
      f"{np.allclose(Operator(qc1_ctrl).data, Operator(qc2_ctrl).data)}")

Bug 4: CNOT direction

In Qiskit, qc.cx(control, target) means the control qubit is the first argument and the target is the second. Swapping them gives a different circuit.

8.12 Quantum Circuit Analysis and Metrics

8.12.1 Quantifying Circuit Quality

Beyond depth and gate count, several metrics characterize circuit quality:

Expressibility measures how uniformly a parameterized circuit covers the Hilbert space. A highly expressible circuit can represent a wide range of quantum states.

from qiskit import QuantumCircuit
from qiskit.circuit import ParameterVector
from qiskit.quantum_info import Statevector
import numpy as np

def compute_expressibility(n_qubits, n_params, circuit_fn, n_samples=10000):
    """Compute the expressibility of a parameterized circuit.

    Measures how close the circuit's output distribution is to
    the Haar (uniform) distribution over the state space.
    """
    # Sample random parameters
    states = []
    for _ in range(n_samples):
        params = np.random.uniform(0, 2*np.pi, n_params)
        qc = circuit_fn(params)
        sv = Statevector.from_instruction(qc)
        states.append(sv.data)

    # Compute pairwise fidelities
    fidelities = []
    for i in range(len(states)):
        for j in range(i+1, len(states)):
            fid = np.abs(np.conj(states[i]) @ states[j])**2
            fidelities.append(fid)

    # Compare with Haar distribution
    # For n qubits, Haar fidelity distribution is (2^n - 1) * (1-x)^(2^n - 2)
    d = 2**n_qubits
    haar_x = np.linspace(0, 1, 100)
    haar_y = (d - 1) * (1 - haar_x)**(d - 2)

    # Expressibility = distance between distributions (lower is better)
    hist, _ = np.histogram(fidelities, bins=50, range=(0, 1), density=True)

    return np.mean(fidelities), np.std(fidelities)

# Example: Compare two ansatzes
def hardware_efficient(params):
    qc = QuantumCircuit(2)
    p = 0
    for layer in range(2):
        for i in range(2):
            qc.ry(params[p], i); p += 1
            qc.rz(params[p], i); p += 1
        qc.cx(0, 1)
    return qc

def simple_ansatz(params):
    qc = QuantumCircuit(2)
    qc.h(0)
    qc.ry(params[0], 1)
    qc.cx(0, 1)
    return qc

mean_he, std_he = compute_expressibility(2, 8, hardware_efficient)
mean_simple, std_simple = compute_expressibility(2, 1, simple_ansatz)

print(f"Hardware-efficient: mean fidelity = {mean_he:.4f}")
print(f"Simple ansatz: mean fidelity = {mean_simple:.4f}")
print("(Lower mean fidelity = higher expressibility)")

8.12.2 Entangling Capability

Entangling capability measures how much entanglement a circuit can generate. A circuit with high entangling capability is more powerful for quantum advantage.

from qiskit import QuantumCircuit
from qiskit.quantum_info import DensityMatrix, partial_trace, entropy
import numpy as np

def compute_entangling_capability(n_qubits, circuit_fn, n_samples=1000):
    """Compute average entanglement entropy of a circuit's output.

    Higher entropy = more entanglement = potentially more powerful.
    """
    entropies = []

    for _ in range(n_samples):
        params = np.random.uniform(0, 2*np.pi, 4)
        qc = circuit_fn(params)

        # Get statevector
        sv = Statevector.from_instruction(qc)
        rho = DensityMatrix(sv)

        # Compute entanglement entropy (trace out qubit 0)
        rho_reduced = partial_trace(rho, [0])

        # Von Neumann entropy
        eigenvalues = np.linalg.eigvalsh(rho_reduced.data)
        eigenvalues = eigenvalues[eigenvalues > 1e-10]  # Remove zeros
        S = -np.sum(eigenvalues * np.log2(eigenvalues))
        entropies.append(S)

    return np.mean(entropies)

print("Entangling capability comparison:")
print(f"  Two-qubit entangling circuit: {compute_entangling_capability(2, lambda p: hardware_efficient(p)):.4f} bits")

8.12.3 Circuit Fidelity Estimation

Before running on hardware, you can estimate the expected fidelity of a circuit:

from qiskit import QuantumCircuit, transpile
from qiskit.providers.fake_provider import FakeBrisbane

def estimate_expected_fidelity(qc, backend):
    """Estimate the expected fidelity of a circuit on a given backend.

    Uses a simple model: F ≈ ∏(1 - ε_i) × e^{-t/T2}
    where ε_i are gate errors and t is the circuit duration.
    """
    props = backend.properties()

    # Transpile to get native gates
    qc_t = transpile(qc, backend, optimization_level=3)

    # Count gates by type
    ops = qc_t.count_ops()

    # Estimate fidelity
    fidelity = 1.0

    # Single-qubit gate errors
    for gate_name in ['rz', 'sx', 'x']:
        count = ops.get(gate_name, 0)
        if count > 0:
            avg_error = np.mean([props.gate_error(gate_name, q) 
                                  for q in range(backend.num_qubits)])
            fidelity *= (1 - avg_error) ** count

    # Two-qubit gate errors
    for gate_name in ['cx', 'ecr']:
        count = ops.get(gate_name, 0)
        if count > 0:
            # Get average error over all connected qubit pairs
            cx_errors = [next(p.value for p in g.parameters if p.name == 'gate_error')
                         for g in props.gates if g.gate == gate_name]
            avg_error = np.mean(cx_errors)
            fidelity *= (1 - avg_error) ** count

    # Readout errors
    readout_errors = [props.readout_error(q) for q in range(backend.num_qubits)]
    avg_readout_error = np.mean(readout_errors)
    measured_qubits = len([g for g in ops if g == 'measure'])
    fidelity *= (1 - avg_readout_error) ** measured_qubits

    return fidelity

# Test with a Bell state circuit
backend = FakeBrisbane()
qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)
qc.measure_all()

fidelity = estimate_expected_fidelity(qc, backend)
print(f"Expected fidelity for Bell state: {fidelity:.4f}")

# Test with a more complex circuit
qc2 = QuantumCircuit(5)
qc2.h(0)
for i in range(4):
    qc2.cx(i, i+1)
qc2.measure_all()

fidelity2 = estimate_expected_fidelity(qc2, backend)
print(f"Expected fidelity for 5-qubit GHZ: {fidelity2:.4f}")

8.11.6 Quantum Volume Benchmarking in Practice

Quantum Volume (QV) is the standard benchmark for quantum computers. A higher QV means the device can reliably execute larger, more complex circuits.

The QV test:

1. Generate random $n$-qubit circuits of depth $n$ using $SU(4)$ gates on random qubit pairs
2. Run each circuit and compare the output distribution to the ideal
3. Compute the "heavy output probability" — the fraction of outputs in the top half of the ideal probability distribution
4. The QV is $2^n$ where $n$ is the largest circuit size for which the heavy output probability exceeds 2/3 with 97.5% confidence

Historical QV milestones:

from qiskit.circuit.library import QuantumVolume
from qiskit.quantum_info import Statevector
import numpy as np

def compute_heavy_output_probability(qc, shots=8192):
    """Compute the heavy output probability for a Quantum Volume circuit."""
    # Get ideal distribution
    ideal_sv = Statevector.from_instruction(qc.remove_final_measurements(inplace=False))
    ideal_probs = ideal_sv.probabilities_dict()

    # Find the median probability
    median_prob = np.median(list(ideal_probs.values()))

    # Heavy outputs: outcomes with probability > median
    heavy_outputs = {k: v for k, v in ideal_probs.items() if v > median_prob}

    # Run on simulator
    from qiskit_aer import AerSimulator
    sim = AerSimulator()
    result = sim.run(qc, shots=shots).result()
    counts = result.get_counts()

    # Compute heavy output fraction
    heavy_count = sum(counts.get(k, 0) for k in heavy_outputs)
    heavy_fraction = heavy_count / shots

    return heavy_fraction

# Test QV for different sizes
print("Quantum Volume Heavy Output Test:")
print(f"{'Qubits':>8s}  {'Depth':>6s}  {'Heavy Output P':>15s}  {'Pass?':>6s}")
print("-" * 45)
for n in [2, 3, 4, 5]:
    qv_circuit = QuantumVolume(n, depth=n, seed=42)
    qv_circuit.measure_all()
    hop = compute_heavy_output_probability(qv_circuit)
    passed = hop > 2/3
    print(f"{n:>8d}  {n:>6d}  {hop:>15.4f}  {'Yes' if passed else 'No':>6s}")

Recurring Theme. We're at the beginning. Quantum Volume has been increasing roughly 10× per year, but we're still in the early stages. Achieving QV of $2^{50}$ (which would require ~50 high-fidelity qubits with low error rates) is a milestone that will enable useful quantum algorithms.

8.13 Working with Real Hardware: A Practical Guide

8.13.1 Choosing the Right Backend

IBM Quantum offers multiple backends with different characteristics. Choosing the right one is crucial:

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()

# List available backends with their key metrics
backends = service.backends()
print(f"{'Backend':<20s} {'Qubits':>6s} {'Queue':>6s} {'Status':>10s}")
print("-" * 50)
for b in sorted(backends, key=lambda x: x.num_qubits):
    status = b.status()
    print(f"{b.name:<20s} {b.num_qubits:>6d} {status.pending_jobs:>6d} "
          f"{'Operational' if status.operational else 'Down':>10s}")

Guidelines for choosing a backend:

1. Small circuits (≤5 qubits): Use the smallest available backend. Less qubits means fewer crosstalk errors.
2. Medium circuits (6-20 qubits): Use a mid-range backend. Check the error map for the best qubit subset.
3. Large circuits (>20 qubits): Use the largest available backend, but expect significant noise.
4. Circuit depth matters: A 2-qubit circuit with depth 100 may fail on the best backend. Always check that depth × error_rate < 1.

8.13.2 Qubit Selection and Layout Optimization

Not all qubits are created equal. On a real device, error rates vary significantly across qubits:

from qiskit import QuantumCircuit, transpile
from qiskit.providers.fake_provider import FakeBrisbane

backend = FakeBrisbane()

# Create a 3-qubit GHZ circuit
qc = QuantumCircuit(3)
qc.h(0)
qc.cx(0, 1)
qc.cx(1, 2)
qc.measure_all()

# Transpile with automatic qubit selection
qc_auto = transpile(qc, backend, optimization_level=3)
print("Auto-selected layout:")
print(f"  Depth: {qc_auto.depth()}")
print(f"  CNOT count: {qc_auto.count_ops().get('cx', 0)}")

# Transpile with manual qubit selection (best 3 connected qubits)
from qiskit.transpiler import CouplingMap

# Find the best 3-qubit chain in the coupling map
props = backend.properties()
coupling_map = backend.configuration().coupling_map

best_chain = None
best_score = float('inf')
for edge1 in coupling_map:
    for edge2 in coupling_map:
        if edge1[1] == edge2[0]:  # chain: a -> b -> c
            q0, q1, q2 = edge1[0], edge1[1], edge2[1]
            if q2 != q0:  # avoid self-loops
                try:
                    score = (props.gate_error('cx', [q0, q1]) +
                            props.gate_error('cx', [q1, q2]) +
                            props.gate_error('x', q0) +
                            props.gate_error('x', q1) +
                            props.gate_error('x', q2))
                    if score < best_score:
                        best_score = score
                        best_chain = [q0, q1, q2]
                except:
                    continue

print(f"\nBest 3-qubit chain: {best_chain} (score: {best_score:.4f})")

# Transpile with the optimal layout
qc_manual = transpile(qc, backend, optimization_level=3,
                       initial_layout=best_chain)
print(f"Manual layout:")
print(f"  Depth: {qc_manual.depth()}")
print(f"  CNOT count: {qc_manual.count_ops().get('cx', 0)}")

8.13.3 Understanding and Mitigating Readout Errors

Readout errors are among the largest error sources on current hardware. A readout error means measuring $|0\rangle$ and getting "1" (or vice versa).

from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel, ReadoutError
import numpy as np

# Model readout errors
p_0_given_0 = 0.97  # P(measure 1 | prepared 0) = 3%
p_1_given_1 = 0.95  # P(measure 0 | prepared 1) = 5%

# Confusion matrix
confusion = np.array([
    [p_0_given_0, 1 - p_0_given_0],  # [P(0|0), P(1|0)]
    [1 - p_1_given_1, p_1_given_1]     # [P(0|1), P(1|1)]
])

print("Readout confusion matrix:")
print(confusion)

# Build noise model with readout errors
noise_model = NoiseModel()
for qubit in range(2):
    noise_model.add_readout_error(ReadoutError(confusion), [qubit])

# Test: prepare |1⟩ and measure
qc = QuantumCircuit(1, 1)
qc.x(0)
qc.measure(0, 0)

sim = AerSimulator(noise_model=noise_model)
result = sim.run(qc, shots=10000).result()
counts = result.get_counts()
print(f"\nPrepared |1⟩, measured: {counts}")
print(f"Readout error rate: {counts.get('0', 0) / 10000:.4f} (expected: {1-p_1_given_1:.4f})")

8.13.4 Session Management for Multiple Circuits

When running many circuits on hardware, use sessions to reduce job overhead:

from qiskit_ibm_runtime import QiskitRuntimeService, Sampler, Session

service = QiskitRuntimeService()
backend = service.backend("ibm_brisbane")

# Create multiple circuits
circuits = []
for theta in np.linspace(0, np.pi, 10):
    qc = QuantumCircuit(1, 1)
    qc.ry(theta, 0)
    qc.measure(0, 0)
    circuits.append(qc)

# Run all circuits in a single session
with Session(backend=backend) as session:
    sampler = Sampler(session=session)
    job = sampler.run(circuits, shots=1024)
    result = job.result()

    for i, theta in enumerate(np.linspace(0, np.pi, 10)):
        quasi_dist = result.quasi_dists[i]
        prob_1 = quasi_dist.get(1, 0)
        print(f"θ={theta:.2f}: P(|1⟩) = {prob_1:.4f} (expected: {np.sin(theta/2)**2:.4f})")

Sessions keep the connection to the backend open between jobs, reducing the overhead from seconds to milliseconds. This is crucial when running parameter sweeps or variational algorithms that require many iterations.

8.13.5 Interpreting Calibration Data

IBM Quantum provides detailed calibration data for each backend. Understanding this data is essential for choosing qubits and estimating circuit fidelity:

from qiskit.providers.fake_provider import FakeBrisbane

backend = FakeBrisbane()
props = backend.properties()

print("=== Qubit Properties ===")
for i in range(min(5, backend.num_qubits)):
    qubit = props.qubits[i]
    t1 = next((p.value for p in qubit if p.name == 'T1'), None)
    t2 = next((p.value for p in qubit if p.name == 'T2'), None)
    freq = next((p.value for p in qubit if p.name == 'frequency'), None)
    readout_err = next((p.value for p in qubit if p.name == 'readout_error'), None)
    print(f"  Qubit {i}: T1={t1:.1f}μs, T2={t2:.1f}μs, "
          f"freq={freq:.3f}GHz, readout_err={readout_err:.4f}")

print("\n=== Two-Qubit Gate Errors (top 10) ===")
cx_gates = [g for g in props.gates if g.gate == 'cx']
cx_gates.sort(key=lambda g: next(p.value for p in g.parameters if p.name == 'gate_error'))
for gate in cx_gates[:10]:
    err = next(p.value for p in gate.parameters if p.name == 'gate_error')
    print(f"  CX({gate.qubits[0]},{gate.qubits[1]}): {err:.4%}")

How to estimate circuit fidelity: A rough estimate of circuit fidelity is:

$$F \approx \prod_{\text{gates}} (1 - \epsilon_{\text{gate}}) \times e^{-t_{\text{circuit}} / T_2}$$

where $\epsilon_{\text{gate}}$ is the error rate for each gate and $t_{\text{circuit}}$ is the circuit duration. For a circuit with $n_1$ single-qubit gates (error $\epsilon_1$), $n_2$ two-qubit gates (error $\epsilon_2$), and total time $t$:

$$F \approx (1-\epsilon_1)^{n_1} \times (1-\epsilon_2)^{n_2} \times e^{-t/T_2}$$

def estimate_circuit_fidelity(qc, backend, t2_avg=100):
    """Estimate the fidelity of a circuit on a given backend."""
    props = backend.properties()

    # Get average gate errors
    sx_errors = [next(p.value for p in g.parameters if p.name == 'gate_error')
                 for g in props.gates if g.gate == 'sx']
    cx_errors = [next(p.value for p in g.parameters if p.name == 'gate_error')
                 for g in props.gates if g.gate == 'cx']

    avg_1q_error = np.mean(sx_errors) if sx_errors else 0.001
    avg_2q_error = np.mean(cx_errors) if cx_errors else 0.01

    # Count gates
    ops = qc.count_ops()
    n_1q = sum(ops.get(g, 0) for g in ['h', 'x', 'y', 'z', 's', 't', 'sx', 'rz'])
    n_2q = ops.get('cx', 0) + ops.get('cz', 0) + ops.get('swap', 0) * 3

    # Estimate circuit time (rough: 20ns per 1q gate, 200ns per 2q gate)
    t_circuit = n_1q * 20 + n_2q * 200  # in nanoseconds

    # Fidelity estimate
    f_gates = (1 - avg_1q_error) ** n_1q * (1 - avg_2q_error) ** n_2q
    f_decoherence = np.exp(-t_circuit / (t2_avg * 1000))  # T2 in microseconds

    return f_gates * f_decoherence, {
        'n_1q': n_1q, 'n_2q': n_2q,
        'avg_1q_error': avg_1q_error, 'avg_2q_error': avg_2q_error,
        'f_gates': f_gates, 'f_decoherence': f_decoherence
    }

# Example
qc = QuantumCircuit(3)
qc.h(0)
qc.cx(0, 1)
qc.cx(1, 2)
qc.measure_all()

fidelity, details = estimate_circuit_fidelity(qc, FakeBrisbane())
print(f"Estimated fidelity: {fidelity:.4f}")
print(f"  1q gates: {details['n_1q']}, 2q gates: {details['n_2q']}")
print(f"  Gate fidelity: {details['f_gates']:.4f}")
print(f"  Decoherence fidelity: {details['f_decoherence']:.4f}")

8.14 Qiskit Pulse-Level Control

For advanced users, Qiskit provides access to pulse-level control of quantum hardware. This allows you to shape the microwave pulses that implement gates, enabling custom calibrations and optimal control experiments.

8.14.1 The Pulse Schedule

A pulse schedule is a sequence of waveforms applied to specific channels (drive, measure, control) at specific times:

from qiskit.pulse import (
    DriveChannel, MeasureChannel, ControlChannel,
    Waveform, SamplePulse, Schedule, Play, Delay
)
import numpy as np

# Create a simple Gaussian pulse
duration = 160  # in dt units (dt ≈ 0.222 ns for IBM devices)
sigma = 40
amp = 0.5

# Gaussian envelope
times = np.arange(duration)
gaussian = amp * np.exp(-((times - duration/2)**2) / (2*sigma**2))

# Create the waveform
waveform = Waveform(samples=gaussian, name='gaussian_pulse')

# Build a schedule
sched = Schedule()
sched = sched.append(Play(waveform, DriveChannel(0)), inplace=True)

print(f"Schedule duration: {sched.duration} dt")
print(f"Number of instructions: {len(sched.instructions)}")

8.14.2 Custom Gate Calibration

You can calibrate a custom gate by defining its pulse implementation:

from qiskit import QuantumCircuit, pulse
from qiskit.pulse import DriveChannel, Play, Waveform
import numpy as np

# Define a custom pulse for a "fast" X gate
duration = 120  # shorter than default
sigma = 30
amp = 0.6
times = np.arange(duration)
samples = amp * np.exp(-((times - duration/2)**2) / (2*sigma**2))

custom_x_pulse = Waveform(samples=samples, name='fast_x')

# Build a circuit with the custom gate
qc = QuantumCircuit(1, 1)
qc.x(0)  # This will use the custom calibration
qc.measure(0, 0)

# Add the calibration
with pulse.build() as custom_x_sched:
    pulse.play(custom_x_pulse, DriveChannel(0))

qc.add_calibration('x', [0], custom_x_sched)

print("Circuit with custom calibration:")
print(qc.draw('text'))

Pulse-level control is essential for: - Optimal control: Designing pulses that minimize gate time or error - Error suppression: Shaping pulses to reduce leakage and crosstalk - Custom gates: Implementing gates not in the native gate set - Characterization: Running Rabi, T1, T2, and other calibration experiments

Common Misconception. "Pulse-level control lets me do things the hardware can't normally do." Not exactly. Pulse control lets you optimize the implementation of existing operations, but you're still constrained by the hardware's frequency ranges, power limits, and bandwidth. You can't create fundamentally new operations — you can only implement existing unitaries more efficiently.

8.15 Dynamic Circuits and Mid-Circuit Measurements

Recent IBM Quantum devices support dynamic circuits — circuits where the result of a mid-circuit measurement determines subsequent operations. This enables feed-forward control, adaptive algorithms, and error correction.

8.15.1 If-Else Operations

from qiskit import QuantumCircuit

# Dynamic circuit: measure qubit 0, and if it's 1, flip qubit 1
qc = QuantumCircuit(2, 2)

qc.h(0)           # Put qubit 0 in superposition
qc.measure(0, 0)  # Mid-circuit measurement

# Conditional operation based on measurement result
with qc.if_test((0, 1)):  # If classical bit 0 == 1
    qc.x(1)

qc.measure(1, 1)  # Measure qubit 1

print(qc.draw('text'))

# Run on a dynamic circuit backend
from qiskit_aer import AerSimulator
sim = AerSimulator()
result = sim.run(qc, shots=1024).result()
counts = result.get_counts()
print(f"\nResults: {counts}")

8.15.2 Real-Time Classical Processing

Dynamic circuits also support real-time classical computation between measurements:

from qiskit import QuantumCircuit

qc = QuantumCircuit(3, 2)

# Prepare a GHZ-like state
qc.h(0)
qc.cx(0, 1)
qc.cx(0, 2)

# Measure qubits 0 and 1
qc.measure(0, 0)
qc.measure(1, 1)

# If the parity of the two measurements is odd, flip qubit 2
with qc.if_test((0, 1)):
    qc.x(2)
with qc.if_test((1, 1)):
    qc.x(2)

qc.measure(2, 0)  # Reuse classical bit 0 for final measurement
print(qc.draw('text'))

8.15.3 Applications of Dynamic Circuits

Dynamic circuits enable several important quantum computing techniques:

1. Quantum error correction: Measure syndromes mid-circuit and apply corrections based on the results
2. Adaptive algorithms: Adjust the circuit based on measurement outcomes (e.g., adaptive phase estimation)
3. Iterative amplitude estimation: Estimate amplitudes with fewer queries by adaptively narrowing the search range
4. Measurement-based quantum computing: Perform operations conditioned on measurement results

Try it yourself: Implement a simple quantum error correction code using dynamic circuits. Encode a logical qubit in 3 physical qubits, measure the syndromes, and apply corrections based on the syndrome results.

8.16 Quantum Volume and Benchmarking

8.16.1 What Is Quantum Volume?

Quantum Volume (QV) is a hardware-agnostic benchmark that measures the largest random circuit of equal width and depth that a quantum computer can successfully implement. It captures qubit count, gate fidelity, connectivity, and error rates in a single number.

The QV protocol: 1. Generate random $n$-qubit circuits of depth $n$ 2. Run each circuit on the hardware 3. Compute the heavy output probability (fraction of outputs in the top half of the ideal distribution) 4. The QV is $2^n$ where $n$ is the largest circuit size where the heavy output probability exceeds 2/3 with 97.5% confidence

from qiskit.circuit.library import QuantumVolume
from qiskit_aer import AerSimulator
import numpy as np

# Generate a Quantum Volume circuit
n_qubits = 4
qv_circuit = QuantumVolume(n_qubits, depth=n_qubits, seed=42)
qv_circuit.measure_all()

print(f"Quantum Volume circuit ({n_qubits} qubits, depth {n_qubits}):")
print(qv_circuit.draw('text', fold=80))

# Simulate
sim = AerSimulator()
result = sim.run(qv_circuit.decompose(), shots=8192).result()
counts = result.get_counts()

# Compute heavy output probability
from qiskit.quantum_info import Statevector
ideal_state = Statevector.from_instruction(qv_circuit.remove_final_measurements(inplace=False))
ideal_probs = ideal_state.probabilities_dict()

# Heavy outputs: outcomes with probability > median
median_prob = np.median(list(ideal_probs.values()))
heavy_outputs = {k: v for k, v in ideal_probs.items() if v > median_prob}
heavy_count = sum(counts.get(k, 0) for k in heavy_outputs)
heavy_prob = heavy_count / 8192

print(f"\nHeavy output probability: {heavy_prob:.4f}")
print(f"Threshold: 2/3 = {2/3:.4f}")
print(f"Quantum Volume: 2^{n_qubits} = {2**n_qubits}" if heavy_prob > 2/3 else "QV test FAILED")

8.16.2 Other Benchmarks

• CLOPS (Circuit Layer Operations Per Second): Measures how many QV layers a backend can execute per second
• Application-level benchmarks: Running specific algorithms (VQE, QAOA, QFT) and measuring accuracy
• Randomized benchmarking: Measuring average gate fidelity using random circuits
• Cross-entropy benchmarking: Comparing measured output distributions with ideal distributions

8.17 Best Practices for Quantum Programming

8.17.1 Code Organization

# GOOD: Modular, reusable circuit construction
def create_ghz_state(n_qubits):
    """Create an n-qubit GHZ state."""
    qc = QuantumCircuit(n_qubits)
    qc.h(0)
    for i in range(n_qubits - 1):
        qc.cx(i, i + 1)
    return qc

def create_bell_state():
    """Create a 2-qubit Bell state."""
    return create_ghz_state(2)

# BAD: Copy-pasting circuit construction everywhere
qc1 = QuantumCircuit(2)
qc1.h(0)
qc1.cx(0, 1)

qc2 = QuantumCircuit(3)
qc2.h(0)
qc2.cx(0, 1)
qc2.cx(1, 2)

8.17.2 Debugging Tips

1. Start with statevector simulation. Before running with shots, verify your circuit produces the correct statevector.

2. Use barriers to isolate errors. Insert qc.barrier() between sections of your circuit to prevent the transpiler from merging or reordering gates across sections.

3. Check unitaries. For small circuits (≤4 qubits), compute the unitary and verify it matches expectations.

4. Add assertions. Use qc.breakpoint() in dynamic circuits to pause execution and inspect intermediate states.

5. Simulate with noise models. Always test with noise before running on hardware. This sets realistic expectations.

6. Compare ideal vs. noisy results. The difference tells you how much noise affects your algorithm.

# Debugging pattern: verify circuit correctness
from qiskit.quantum_info import Statevector, Operator

qc = create_bell_state()
qc_no_measure = qc.remove_final_measurements(inplace=False)

# Check statevector
state = Statevector.from_instruction(qc_no_measure)
print(f"State: {state}")
print(f"Probabilities: {state.probabilities_dict()}")

# Check unitary
U = Operator(qc_no_measure)
print(f"Unitary shape: {U.data.shape}")
print(f"Is unitary: {np.allclose(U.data @ U.data.conj().T, np.eye(4))}")

8.17.3 Performance Tips

1. Minimize circuit depth. Shorter circuits have higher fidelity.
2. Use optimization_level=3 for production. It's slower but produces better circuits.
3. Choose good qubits. Use transpile with initial_layout to place your circuit on the best qubits.
4. Use sessions. Qiskit Runtime sessions reduce job overhead for multiple circuits.
5. Batch your circuits. Submit multiple circuits in a single job rather than separate jobs.
6. Read the calibration data. Choose qubits with the lowest error rates and CNOT connections with the lowest two-qubit errors.

# Choosing the best qubits
from qiskit.providers.fake_provider import FakeBrisbane

backend = FakeBrisbane()
props = backend.properties()

# Find the best qubit (lowest error rate)
best_qubit = min(range(backend.num_qubits),
                 key=lambda q: props.gate_error('x', q))
print(f"Best qubit: {best_qubit} (X error: {props.gate_error('x', best_qubit):.4%})")

# Find the best CNOT connection
best_cx = min(props.gates_of_type('cx'),
              key=lambda g: next(p.value for p in g.parameters if p.name == 'gate_error'))
print(f"Best CNOT: qubits {best_cx.qubits} (error: {next(p.value for p in best_cx.parameters if p.name == 'gate_error'):.4%})")