Noise Modeling in Cirq

Why Noise Modeling Matters

Real quantum hardware is noisy. Gates are imperfect, qubits lose coherence over time, and readout introduces errors. Before running a circuit on physical hardware, simulating its behavior under realistic noise helps you answer three questions:

• Will this algorithm produce useful output at this error rate?
• Which parts of the circuit accumulate the most error?
• Can error mitigation techniques recover the signal?

Cirq provides a clean set of noise channels and a density matrix simulator that lets you answer all three questions on a laptop. This tutorial covers the theory behind quantum noise, the full range of Cirq’s noise modeling tools, and practical techniques for benchmarking and mitigating errors.

Noise Channels in Cirq

A quantum noise channel is a completely positive, trace-preserving (CPTP) map. Unlike ideal unitary operations, channels can represent irreversible processes like energy loss and decoherence. Cirq implements the most common channels as first-class objects.

Channel	cirq class	Physical process
Depolarizing	cirq.depolarize(p)	Random Pauli error with probability p
Amplitude damping	cirq.amplitude_damp(gamma)	Energy loss (T1 decay)
Bit flip	cirq.bit_flip(p)	Random X error
Phase flip	cirq.phase_flip(p)	Random Z error
Phase damping	cirq.phase_damp(gamma)	Pure dephasing (T2 with no T1)

Each channel is specified by one or two parameters. The depolarizing channel depolarize(p) applies X, Y, or Z with probability p/3 each, leaving the state unchanged with probability 1-p. This is the most common channel used for generic noise benchmarking.

CPTP Maps: Why Quantum Noise Has Rules

Not every linear map on density matrices corresponds to a physically realizable process. Quantum mechanics imposes two strict constraints on any valid noise channel E:

Trace-preserving (TP): For any input state rho, Tr(E(rho)) = Tr(rho) = 1. This ensures that probabilities still sum to 1 after the noise acts. If a map is not trace-preserving, it creates or destroys probability, which is unphysical.

Completely positive (CP): The map E must be positive (it sends positive semidefinite matrices to positive semidefinite matrices), and it must remain positive even when extended to act on a larger system. Formally, for any ancilla system R of any dimension, the extended map (E tensor I_R) must also be positive. This is a stronger requirement than simple positivity, and it matters precisely because quantum systems can be entangled with other systems.

Why does complete positivity matter in practice? Consider a map that is positive but not completely positive. If your qubit is entangled with another qubit that the map does not touch, applying a merely positive map can produce a “density matrix” with negative eigenvalues. That would correspond to negative probabilities, which have no physical meaning.

The transposition map as a counterexample. The transpose map T, defined by T(rho) = rho^T, is positive: if rho has non-negative eigenvalues, so does rho^T (since transposition preserves eigenvalues). It is also trace-preserving, since Tr(rho^T) = Tr(rho). Yet T is not completely positive, and therefore it is not a valid quantum channel.

To see the failure, consider the maximally entangled Bell state on two qubits:

|Phi+> = (|00> + |11>) / sqrt(2)

Its density matrix is:

rho_AB = |Phi+><Phi+| = (1/2) * [[1,0,0,1],
                                   [0,0,0,0],
                                   [0,0,0,0],
                                   [1,0,0,1]]

Applying the transpose to qubit B only (the partial transpose) gives:

(I_A tensor T_B)(rho_AB) = (1/2) * [[1,0,0,0],
                                      [0,0,1,0],
                                      [0,1,0,0],
                                      [0,0,0,1]]

This matrix has eigenvalues {1/2, 1/2, 1/2, -1/2}. The negative eigenvalue proves that the output is not a valid density matrix. The transpose map is not physically realizable. This property is actually useful: the partial transpose test (the Peres-Horodecki criterion) uses exactly this fact to detect entanglement. If the partial transpose of a bipartite state has a negative eigenvalue, the state is entangled.

The following code verifies this in NumPy:

import numpy as np

# Bell state density matrix
bell = np.array([1, 0, 0, 1]) / np.sqrt(2)
rho_ab = np.outer(bell, bell.conj())

# Partial transpose on qubit B: swap (01) <-> (10) blocks
# For a 2x2 block structure, partial transpose on B transposes each 2x2 block
rho_pt = rho_ab.reshape(2, 2, 2, 2).transpose(0, 3, 2, 1).reshape(4, 4)

eigenvalues = np.linalg.eigvalsh(rho_pt)
print("Eigenvalues of partial transpose:", np.round(eigenvalues, 6))
# Output: [-0.5, 0.5, 0.5, 0.5]
# The negative eigenvalue proves the transpose map is not CP

Every noise channel in Cirq is CPTP by construction. When you build a custom channel from Kraus operators, Cirq verifies the completeness relation sum_i K_i^dagger K_i = I at construction time.

Kraus Operator Representation

Every CPTP map can be written in the Kraus form:

E(rho) = sum_i K_i rho K_i^dagger

where the operators {K_i} satisfy the completeness relation sum_i K_i^dagger K_i = I. This is the operator-sum representation. For amplitude damping with parameter gamma, the two Kraus operators are:

K_0 = [[1, 0], [0, sqrt(1-gamma)]]   (no decay)
K_1 = [[0, sqrt(gamma)], [0, 0]]      (decay to |0>)

You can verify the completeness relation: K_0^dagger K_0 + K_1^dagger K_1 = [[1,0],[0,1-gamma]] + [[0,0],[0,gamma]] = I.

Cirq’s built-in channels handle Kraus construction internally. Understanding the formalism matters when you define a custom channel, which we cover in later sections.

Stinespring Dilation: The Environment Model of Noise

The Kraus representation tells you how to compute the effect of noise, but the Stinespring dilation theorem tells you where noise comes from physically. Every CPTP map can be modeled as a unitary interaction between the system and an environment that starts in a known pure state, followed by tracing out the environment:

E(rho) = Tr_E[ U (rho tensor |0_E><0_E|) U^dagger ]

Here U acts on the joint system-environment Hilbert space, and |0_E> is the initial state of the environment. This is the Stinespring dilation. It says that all quantum noise is fundamentally unitary evolution on a larger system. The non-unitary, irreversible character of the noise appears only because we lose access to the environment degrees of freedom.

Amplitude damping as energy exchange with the environment. For a qubit coupled to a zero-temperature bath (the environment starts in |0>), the physical picture is spontaneous emission: the qubit in |1> can emit a photon into the environment, decaying to |0>. The joint unitary acts as:

U|0>|0_E> = |0>|0_E>                                         (ground state is stable)
U|1>|0_E> = sqrt(1-gamma)|1>|0_E> + sqrt(gamma)|0>|1_E>      (possible decay)

The environment state |1_E> represents “one photon emitted.” Writing this as a 4x4 matrix in the basis {|0,0_E>, |0,1_E>, |1,0_E>, |1,1_E>}:

U = [[1,         0,          0,            0        ],
     [0,         1,          0,            0        ],
     [0,         0,          sqrt(1-gamma), 0       ],
     [0,         0,          sqrt(gamma),   0       ]]

Wait, that is not unitary because the columns are not all normalized. The correct Stinespring unitary for amplitude damping is:

U = [[1,  0,            0,               0           ],
     [0,  0,            sqrt(gamma),     sqrt(1-gamma)],
     [0,  0,            sqrt(1-gamma),  -sqrt(gamma)  ],
     [0,  1,            0,               0            ]]

The following code constructs this unitary and verifies it reproduces the amplitude damping Kraus operators:

import numpy as np

gamma = 0.1

# Stinespring unitary for amplitude damping
# Basis order: |0,0_E>, |0,1_E>, |1,0_E>, |1,1_E>
U = np.array([
    [1, 0,                0,               0              ],
    [0, 0, np.sqrt(gamma),    np.sqrt(1-gamma)             ],
    [0, 0, np.sqrt(1-gamma), -np.sqrt(gamma)               ],
    [0, 1,                0,               0              ]
])

# Verify unitarity
print("U^dagger U = I?", np.allclose(U.conj().T @ U, np.eye(4)))

# Extract Kraus operators by projecting environment onto <0_E| and <1_E|
# K_k = <k_E| U |0_E>  (environment starts in |0_E>)
# System basis: {|0>, |1>}, Environment basis: {|0_E>, |1_E>}
# |0>|0_E> -> col 0, |0>|1_E> -> col 1, |1>|0_E> -> col 2, |1>|1_E> -> col 3

# K_0 = <0_E| U |0_E>: rows where env=0 (rows 0,2), cols where env=0 (cols 0,2)
K0_stinespring = U[np.ix_([0, 2], [0, 2])]
# K_1 = <1_E| U |0_E>: rows where env=1 (rows 1,3), cols where env=0 (cols 0,2)
K1_stinespring = U[np.ix_([1, 3], [0, 2])]

# Expected Kraus operators
K0_expected = np.array([[1, 0], [0, np.sqrt(1 - gamma)]])
K1_expected = np.array([[0, np.sqrt(gamma)], [0, 0]])

print("K0 matches?", np.allclose(K0_stinespring, K0_expected))
print("K1 matches?", np.allclose(K1_stinespring, K1_expected))

The Stinespring picture is valuable for building intuition. When you apply cirq.amplitude_damp(gamma) in your circuit, you can think of it as the qubit briefly interacting with its electromagnetic environment, with probability gamma of emitting a photon and falling to |0>.

DensityMatrixSimulator

Pure statevectors cannot represent mixed states. Noise creates statistical mixtures of quantum states that require the density matrix formalism: rho = sum_i p_i |psi_i><psi_i|. Cirq’s DensityMatrixSimulator tracks the full density matrix and applies channels exactly.

import cirq
import numpy as np

# Use DensityMatrixSimulator for any noisy circuit
sim = cirq.DensityMatrixSimulator()

The memory cost scales as 4^n complex numbers for n qubits, compared to 2^n for a pure statevector. In practice, density matrix simulation is tractable up to about 12 to 15 qubits on a standard workstation. At 15 qubits, the density matrix is 32768 x 32768 complex128, consuming about 8 GB of RAM. Plan your qubit counts accordingly.

Example 1: Depolarizing Noise on a Bell State

import cirq
import numpy as np

q0, q1 = cirq.LineQubit.range(2)
p_error = 0.02   # 2% depolarizing probability per gate

# Ideal Bell state circuit
ideal_circuit = cirq.Circuit([
    cirq.H(q0),
    cirq.CNOT(q0, q1),
])

# Noisy version: add depolarizing channel after each gate
noisy_circuit = cirq.Circuit([
    cirq.H(q0),
    cirq.depolarize(p_error)(q0),
    cirq.CNOT(q0, q1),
    cirq.depolarize(p_error)(q0),
    cirq.depolarize(p_error)(q1),
])

sim = cirq.DensityMatrixSimulator()

ideal_result = sim.simulate(ideal_circuit)
noisy_result = sim.simulate(noisy_circuit)

ideal_dm = ideal_result.final_density_matrix
noisy_dm = noisy_result.final_density_matrix

print("Ideal Bell state density matrix (diagonal):")
print(np.real(np.diag(ideal_dm)))

print("\nNoisy Bell state density matrix (diagonal):")
print(np.real(np.diag(noisy_dm)))

# Compute fidelity between ideal and noisy
fidelity = cirq.fidelity(ideal_dm, noisy_dm, qid_shape=(2, 2))
print(f"\nFidelity: {fidelity:.4f}")

The ideal Bell state has density matrix with 0.5 in the |00><00| and |11><11| positions. Depolarizing noise pushes probability mass into the off-diagonal elements and slightly populates |01> and |10>, reducing fidelity below 1.

Example 2: T1 Amplitude Damping

T1 is the relaxation time of a qubit: how long it takes to decay from |1> back to |0>. Given a gate time t_gate and a T1 time, the effective gamma per gate is:

gamma = 1 - exp(-t_gate / T1)

import cirq
import numpy as np

T1_ns = 50_000   # 50 microseconds, typical for superconducting qubit
t_gate_ns = 50   # 50 nanoseconds per gate
gamma = 1 - np.exp(-t_gate_ns / T1_ns)

q = cirq.LineQubit(0)

# Prepare |1> and apply several gates with T1 damping after each
circuit = cirq.Circuit()
circuit.append(cirq.X(q))   # Start in |1>

num_gates = 20
for _ in range(num_gates):
    circuit.append(cirq.I(q))                         # Identity gate (idle)
    circuit.append(cirq.amplitude_damp(gamma)(q))     # T1 decay per gate

sim = cirq.DensityMatrixSimulator()
result = sim.simulate(circuit)
dm = result.final_density_matrix

excited_prob = np.real(dm[1, 1])
print(f"P(|1>) after {num_gates} gates with T1 damping: {excited_prob:.4f}")
print(f"Expected from T1: {np.exp(-num_gates * t_gate_ns / T1_ns):.4f}")

Thermal Relaxation: Combining T1 and T2

Real qubits experience both amplitude damping (T1) and phase damping (T2) simultaneously. The T2 time characterizes the total loss of phase coherence and satisfies the constraint T2 <= 2 * T1. When T2 < 2 * T1, there is “pure dephasing” on top of the T1-induced dephasing.

The combined thermal relaxation channel for a gate of duration t has the following structure. First, compute the amplitude damping parameter and the pure dephasing parameter:

gamma_1 = 1 - exp(-t / T1)              # amplitude damping rate
gamma_phi = 1 - exp(-t/T2 + t/(2*T1))   # pure dephasing rate (extra dephasing beyond T1)

The thermal relaxation channel applies amplitude damping with parameter gamma_1, followed by phase damping with the residual dephasing rate gamma_phi. In Kraus operator form, the combined channel has three operators that account for both decay and dephasing.

The following code builds the thermal relaxation channel from T1, T2, and gate time, then applies it in Cirq:

import cirq
import numpy as np

def thermal_relaxation_kraus(T1: float, T2: float, t: float):
    """Build Kraus operators for the thermal relaxation channel.

    Args:
        T1: Amplitude damping time (same units as t).
        T2: Total dephasing time, must satisfy T2 <= 2*T1.
        t: Gate duration (same units as T1, T2).

    Returns:
        List of 2x2 numpy arrays (Kraus operators).
    """
    assert T2 <= 2 * T1, f"T2={T2} must be <= 2*T1={2*T1}"

    # Amplitude damping parameter
    p_decay = 1 - np.exp(-t / T1)

    # Pure dephasing parameter (dephasing beyond what T1 causes)
    # T1-induced dephasing rate is 1/(2*T1), total dephasing rate is 1/T2
    # Pure dephasing rate is 1/T2 - 1/(2*T1)
    if T2 < 2 * T1:
        p_pure_dephase = 1 - np.exp(-t / T2 + t / (2 * T1))
    else:
        p_pure_dephase = 0.0

    # Amplitude damping Kraus operators
    K_ad_0 = np.array([[1, 0],
                        [0, np.sqrt(1 - p_decay)]])
    K_ad_1 = np.array([[0, np.sqrt(p_decay)],
                        [0, 0]])

    # Phase damping (pure dephasing) Kraus operators
    K_pd_0 = np.array([[1, 0],
                        [0, np.sqrt(1 - p_pure_dephase)]])
    K_pd_1 = np.array([[0, 0],
                        [0, np.sqrt(p_pure_dephase)]])

    # Combined channel: compose by taking all products K_pd_i @ K_ad_j
    kraus_ops = []
    for K_pd in [K_pd_0, K_pd_1]:
        for K_ad in [K_ad_0, K_ad_1]:
            K_combined = K_pd @ K_ad
            if np.linalg.norm(K_combined) > 1e-10:
                kraus_ops.append(K_combined)

    # Verify completeness: sum K_i^dagger K_i = I
    check = sum(K.conj().T @ K for K in kraus_ops)
    assert np.allclose(check, np.eye(2)), f"Completeness check failed:\n{check}"

    return kraus_ops

# Typical superconducting qubit parameters
T1 = 50e3    # 50 microseconds in nanoseconds
T2 = 70e3    # 70 microseconds (< 2*T1 = 100 microseconds)
t_gate = 50  # 50 ns gate time

kraus_ops = thermal_relaxation_kraus(T1, T2, t_gate)
print(f"Number of Kraus operators: {len(kraus_ops)}")

# Build a Cirq channel from these Kraus operators
thermal_channel = cirq.KrausChannel(kraus_ops, key='thermal_relax')

q = cirq.LineQubit(0)
circuit = cirq.Circuit()
circuit.append(cirq.X(q))  # Start in |1>

# Apply 100 idle gates with thermal relaxation
for _ in range(100):
    circuit.append(cirq.I(q))
    circuit.append(thermal_channel.on(q))

sim = cirq.DensityMatrixSimulator()
result = sim.simulate(circuit)
dm = result.final_density_matrix

print(f"P(|1>) after 100 idle gates: {np.real(dm[1, 1]):.4f}")
print(f"Off-diagonal |rho_01|: {np.abs(dm[0, 1]):.6f}")
print(f"Expected P(|1>) from T1 alone: {np.exp(-100 * t_gate / T1):.4f}")

The off-diagonal element decays faster than what T1 alone would predict, because the pure dephasing from T2 contributes additional decoherence. This distinction matters when simulating algorithms that rely on phase coherence, such as the Quantum Fourier Transform.

Readout Error Modeling

Gate noise is not the only source of error. The measurement process itself is imperfect. For superconducting qubits, typical readout error rates are:

• p(read 1 | state is 0): 0.5% to 1% (false positive)
• p(read 0 | state is 1): 1% to 3% (false negative, higher because |1> can decay to |0> during readout)

Cirq does not have a built-in readout error channel, but you can model it by inserting a bit-flip channel immediately before measurement. The trick is to use asymmetric bit-flip rates by constructing a custom channel.

For a simpler symmetric model, a single cirq.bit_flip(p) before measurement works:

import cirq
import numpy as np

q0, q1 = cirq.LineQubit.range(2)

# Readout error probabilities
p_readout = 0.02  # 2% symmetric readout error

# Bell state with readout noise
circuit = cirq.Circuit([
    cirq.H(q0),
    cirq.CNOT(q0, q1),
    # Insert bit-flip channels before measurement to model readout error
    cirq.bit_flip(p_readout)(q0),
    cirq.bit_flip(p_readout)(q1),
    cirq.measure(q0, q1, key='result'),
])

sim = cirq.DensityMatrixSimulator()
result = sim.run(circuit, repetitions=10000)
counts = result.histogram(key='result')

print("Measurement outcomes with readout error:")
for bitstring, count in sorted(counts.items()):
    label = format(bitstring, '02b')
    print(f"  |{label}>: {count} ({count/100:.1f}%)")

With no readout error, you expect roughly 50% |00> and 50% |11>. With 2% readout error per qubit, each qubit has a 2% chance of its bit being flipped, so you see a few percent of |01> and |10> outcomes contaminating the results.

For asymmetric readout errors (different rates for 0-to-1 and 1-to-0 flips), build a custom Kraus channel:

import cirq
import numpy as np

def asymmetric_readout_error(p_0_to_1: float, p_1_to_0: float):
    """Create an asymmetric readout error channel.

    p_0_to_1: probability of reading 1 when the state is 0.
    p_1_to_0: probability of reading 0 when the state is 1.
    """
    # K0: no flip, with state-dependent survival probabilities
    K0 = np.array([[np.sqrt(1 - p_0_to_1), 0],
                    [0, np.sqrt(1 - p_1_to_0)]])
    # K1: flip 0 -> 1
    K1 = np.array([[0, 0],
                    [np.sqrt(p_0_to_1), 0]])
    # K2: flip 1 -> 0
    K2 = np.array([[0, np.sqrt(p_1_to_0)],
                    [0, 0]])

    return cirq.KrausChannel([K0, K1, K2], key='readout_err')

q0, q1 = cirq.LineQubit.range(2)
readout_ch = asymmetric_readout_error(p_0_to_1=0.005, p_1_to_0=0.02)

circuit = cirq.Circuit([
    cirq.H(q0),
    cirq.CNOT(q0, q1),
    readout_ch.on(q0),
    readout_ch.on(q1),
    cirq.measure(q0, q1, key='result'),
])

result = cirq.DensityMatrixSimulator().run(circuit, repetitions=10000)
counts = result.histogram(key='result')
print("Asymmetric readout error results:")
for bitstring, count in sorted(counts.items()):
    label = format(bitstring, '02b')
    print(f"  |{label}>: {count} ({count/100:.1f}%)")

Example 3: Custom Noise Model

Subclass cirq.NoiseModel to apply gate-specific error rates. For example, a lower error rate on single-qubit gates than on two-qubit gates:

import cirq

class TwoRateNoiseModel(cirq.NoiseModel):
    def __init__(self, p1q: float, p2q: float):
        self.p1q = p1q
        self.p2q = p2q

    def noisy_operation(self, op: cirq.Operation) -> cirq.OP_TREE:
        yield op
        if len(op.qubits) == 1:
            yield cirq.depolarize(self.p1q)(op.qubits[0])
        elif len(op.qubits) == 2:
            for q in op.qubits:
                yield cirq.depolarize(self.p2q)(q)

q0, q1 = cirq.LineQubit.range(2)
base_circuit = cirq.Circuit([cirq.H(q0), cirq.CNOT(q0, q1)])
noise_model = TwoRateNoiseModel(p1q=0.001, p2q=0.01)
noisy_circuit = cirq.Circuit(noise_model.noisy_moments(base_circuit, base_circuit.all_qubits()))
result = cirq.DensityMatrixSimulator().simulate(noisy_circuit)
print(result.final_density_matrix.round(4))

ConstantQubitNoiseModel and Simplified Noise Injection

For uniform noise across all gates, Cirq provides cirq.ConstantQubitNoiseModel, which applies the same noise channel after every gate on every qubit. This is simpler than writing a full NoiseModel subclass when you want a uniform depolarizing error rate:

import cirq
import numpy as np

q0, q1 = cirq.LineQubit.range(2)

# Uniform depolarizing noise after every gate
uniform_noise = cirq.ConstantQubitNoiseModel(cirq.depolarize(0.02))

# Build a clean circuit, then generate the noisy version
clean_circuit = cirq.Circuit([
    cirq.H(q0),
    cirq.CNOT(q0, q1),
])

noisy_circuit = cirq.Circuit(
    uniform_noise.noisy_moments(clean_circuit, clean_circuit.all_qubits())
)

print("Noisy circuit with ConstantQubitNoiseModel:")
print(noisy_circuit)

sim = cirq.DensityMatrixSimulator()
result = sim.simulate(noisy_circuit)
fidelity_uniform = cirq.fidelity(
    sim.simulate(clean_circuit).final_density_matrix,
    result.final_density_matrix,
    qid_shape=(2, 2),
)
print(f"Fidelity (uniform noise): {fidelity_uniform:.4f}")

Compare this with the manual TwoRateNoiseModel from the previous section. ConstantQubitNoiseModel applies the same noise to every qubit involved in every operation, regardless of whether it is a single-qubit or two-qubit gate. The TwoRateNoiseModel lets you set different error rates for different gate types, which is more realistic: on real hardware, two-qubit gates typically have error rates 5x to 10x higher than single-qubit gates.

For more flexible noise injection based on gate type, you can use cirq.devices.InsertionNoiseModel, which lets you map specific gate types to specific noise channels:

import cirq
import numpy as np

# Map gate types to noise channels
noise_map = {
    cirq.ops.common_gates.HPowGate: cirq.depolarize(0.001),
    cirq.ops.common_gates.CNotPowGate: cirq.depolarize(0.01),
}

class GateTypeNoiseModel(cirq.NoiseModel):
    """Apply noise based on the gate type of each operation."""
    def __init__(self, gate_noise_map):
        self.gate_noise_map = gate_noise_map

    def noisy_operation(self, op: cirq.Operation) -> cirq.OP_TREE:
        yield op
        gate_type = type(op.gate)
        if gate_type in self.gate_noise_map:
            noise_channel = self.gate_noise_map[gate_type]
            for q in op.qubits:
                yield noise_channel(q)

q0, q1 = cirq.LineQubit.range(2)
clean = cirq.Circuit([cirq.H(q0), cirq.CNOT(q0, q1)])
model = GateTypeNoiseModel(noise_map)
noisy = cirq.Circuit(model.noisy_moments(clean, clean.all_qubits()))
print("Gate-type-specific noisy circuit:")
print(noisy)

Calibration-Based Noise Model from Device Data

Real quantum devices publish calibration data that reports per-qubit, per-gate error rates. You can use this data to build a noise model that reflects the actual performance of a specific chip. The following example simulates loading calibration data and constructing a per-qubit noise model:

import cirq
import numpy as np

# Simulated calibration data: (qubit, gate_type) -> error_rate
# In practice, this comes from a device calibration JSON or API
calibration_data = {
    (cirq.LineQubit(0), 'single'): 0.0008,
    (cirq.LineQubit(1), 'single'): 0.0012,
    (cirq.LineQubit(2), 'single'): 0.0010,
    (cirq.LineQubit(0), 'two'):    0.0085,
    (cirq.LineQubit(1), 'two'):    0.0120,
    (cirq.LineQubit(2), 'two'):    0.0095,
}

class CalibrationNoiseModel(cirq.NoiseModel):
    """Noise model driven by per-qubit calibration data."""

    def __init__(self, cal_data: dict, default_1q: float = 0.001,
                 default_2q: float = 0.01):
        self.cal_data = cal_data
        self.default_1q = default_1q
        self.default_2q = default_2q

    def _get_error_rate(self, qubit, gate_type):
        return self.cal_data.get((qubit, gate_type), 
                                 self.default_1q if gate_type == 'single' 
                                 else self.default_2q)

    def noisy_operation(self, op: cirq.Operation) -> cirq.OP_TREE:
        yield op
        n_qubits = len(op.qubits)
        gate_type = 'single' if n_qubits == 1 else 'two'

        for q in op.qubits:
            p = self._get_error_rate(q, gate_type)
            if p > 0:
                yield cirq.depolarize(p)(q)

# Build and run a 3-qubit GHZ circuit with calibration-based noise
qubits = cirq.LineQubit.range(3)
ghz_circuit = cirq.Circuit([
    cirq.H(qubits[0]),
    cirq.CNOT(qubits[0], qubits[1]),
    cirq.CNOT(qubits[1], qubits[2]),
])

cal_model = CalibrationNoiseModel(calibration_data)
noisy_ghz = cirq.Circuit(
    cal_model.noisy_moments(ghz_circuit, ghz_circuit.all_qubits())
)

sim = cirq.DensityMatrixSimulator()
ideal_dm = sim.simulate(ghz_circuit).final_density_matrix
noisy_dm = sim.simulate(noisy_ghz).final_density_matrix

fidelity = cirq.fidelity(ideal_dm, noisy_dm, qid_shape=(2, 2, 2))
print(f"GHZ fidelity with calibration noise: {fidelity:.4f}")

# Show per-qubit error rates used
for (qubit, gate_type), rate in sorted(calibration_data.items(), 
                                         key=lambda x: (x[0][0], x[0][1])):
    print(f"  {qubit} {gate_type}: {rate*100:.2f}%")

This pattern scales to dozens of qubits. In production workflows with Google’s quantum hardware, cirq_google provides calibration loading utilities that populate similar data structures from the device’s most recent calibration report.

Crosstalk Noise Model

Crosstalk occurs when a gate applied to one set of qubits inadvertently affects neighboring qubits. This is a significant error source in multi-qubit devices, particularly for frequency-crowded superconducting architectures. The following noise model applies a small depolarizing error to neighboring qubits whenever a gate is executed:

import cirq
import numpy as np

class CrosstalkNoiseModel(cirq.NoiseModel):
    """Applies depolarizing crosstalk to neighboring qubits during gates."""

    def __init__(self, gate_error: float, crosstalk_error: float,
                 all_qubits: list):
        self.gate_error = gate_error
        self.crosstalk_error = crosstalk_error
        self.all_qubits = set(all_qubits)

    def _neighbors(self, qubit):
        """Return neighboring qubits on a linear chain."""
        if isinstance(qubit, cirq.LineQubit):
            candidates = [cirq.LineQubit(qubit.x - 1),
                          cirq.LineQubit(qubit.x + 1)]
            return [q for q in candidates if q in self.all_qubits]
        return []

    def noisy_operation(self, op: cirq.Operation) -> cirq.OP_TREE:
        yield op

        # Gate noise on operated qubits
        for q in op.qubits:
            yield cirq.depolarize(self.gate_error)(q)

        # Crosstalk noise on neighboring qubits not involved in the gate
        affected = set()
        for q in op.qubits:
            for neighbor in self._neighbors(q):
                if neighbor not in op.qubits:
                    affected.add(neighbor)
        for q in affected:
            yield cirq.depolarize(self.crosstalk_error)(q)

# Compare a 5-qubit circuit with and without crosstalk
qubits = cirq.LineQubit.range(5)

circuit = cirq.Circuit([
    cirq.H(qubits[0]),
    cirq.CNOT(qubits[0], qubits[1]),
    cirq.CNOT(qubits[1], qubits[2]),
    cirq.CNOT(qubits[2], qubits[3]),
    cirq.CNOT(qubits[3], qubits[4]),
])

# Model without crosstalk
no_crosstalk = TwoRateNoiseModel(p1q=0.001, p2q=0.01)
noisy_no_xt = cirq.Circuit(
    no_crosstalk.noisy_moments(circuit, circuit.all_qubits())
)

# Model with crosstalk
with_crosstalk = CrosstalkNoiseModel(
    gate_error=0.01,
    crosstalk_error=0.001,
    all_qubits=qubits,
)
noisy_xt = cirq.Circuit(
    with_crosstalk.noisy_moments(circuit, circuit.all_qubits())
)

sim = cirq.DensityMatrixSimulator()
ideal_dm = sim.simulate(circuit).final_density_matrix
dm_no_xt = sim.simulate(noisy_no_xt).final_density_matrix
dm_xt = sim.simulate(noisy_xt).final_density_matrix

shape = tuple(2 for _ in qubits)
fid_no_xt = cirq.fidelity(ideal_dm, dm_no_xt, qid_shape=shape)
fid_xt = cirq.fidelity(ideal_dm, dm_xt, qid_shape=shape)

print(f"Fidelity without crosstalk: {fid_no_xt:.4f}")
print(f"Fidelity with crosstalk:    {fid_xt:.4f}")
print(f"Fidelity reduction from crosstalk: {fid_no_xt - fid_xt:.4f}")

The crosstalk effect compounds with circuit width. For narrow circuits (2 qubits), crosstalk has little impact because there are few neighboring spectator qubits. For wide circuits (5+ qubits), the additional error on spectator qubits accumulates and noticeably reduces fidelity.

Fidelity Decay Analysis: Error Budgets for GHZ States

Understanding how fidelity scales with error rate and qubit count is essential for estimating whether a given algorithm will produce useful results on a target device. The GHZ state circuit (one Hadamard followed by n-1 CNOTs) is a good testbed because its depth grows linearly with n.

import cirq
import numpy as np

def ghz_fidelity(n_qubits: int, p_error: float) -> float:
    """Compute fidelity of a noisy GHZ state.

    Args:
        n_qubits: Number of qubits in the GHZ state.
        p_error: Depolarizing error probability per gate per qubit.

    Returns:
        Fidelity between ideal and noisy GHZ states.
    """
    qubits = cirq.LineQubit.range(n_qubits)

    # Build ideal GHZ circuit
    ideal = cirq.Circuit()
    ideal.append(cirq.H(qubits[0]))
    for i in range(n_qubits - 1):
        ideal.append(cirq.CNOT(qubits[i], qubits[i + 1]))

    # Build noisy version
    noise_model = TwoRateNoiseModel(p1q=p_error, p2q=p_error)
    noisy = cirq.Circuit(
        noise_model.noisy_moments(ideal, ideal.all_qubits())
    )

    sim = cirq.DensityMatrixSimulator()
    ideal_dm = sim.simulate(ideal).final_density_matrix
    noisy_dm = sim.simulate(noisy).final_density_matrix
    shape = tuple(2 for _ in qubits)
    return float(cirq.fidelity(ideal_dm, noisy_dm, qid_shape=shape))

# Sweep 1: Fixed n=5, varying error rate
print("Sweep 1: n=5 qubits, varying depolarizing error rate")
print(f"{'p_error':>10} {'Fidelity':>10}")
print("-" * 22)
for p in [0.001, 0.005, 0.01, 0.02, 0.03, 0.05]:
    fid = ghz_fidelity(5, p)
    print(f"{p:>10.3f} {fid:>10.4f}")

print()

# Sweep 2: Fixed p=0.01, varying number of qubits
print("Sweep 2: p=0.01 depolarizing, varying qubit count")
print(f"{'n_qubits':>10} {'Fidelity':>10}")
print("-" * 22)
for n in [2, 3, 4, 5, 6, 7, 8, 9, 10]:
    fid = ghz_fidelity(n, 0.01)
    print(f"{n:>10d} {fid:>10.4f}")

The results demonstrate two key scaling behaviors:

1. Error rate sweep (fixed n=5): Fidelity decreases roughly linearly for small p and more steeply as p increases. At p=0.001 (0.1% error), fidelity stays above 0.99. At p=0.05 (5% error), fidelity drops to around 0.6.

2. Qubit count sweep (fixed p=0.01): Fidelity decreases approximately exponentially with qubit count, because each additional qubit adds one more CNOT gate. This exponential scaling is the central challenge of NISQ computing: to maintain fidelity above a useful threshold (say 0.9), you need error rates that decrease as your circuit grows.

These scaling laws help you estimate the required error rate for a target fidelity. For example, if you need F > 0.95 for a 10-qubit GHZ state, you can use this sweep to identify the maximum tolerable error rate.

Noise Fingerprinting with Randomized Benchmarking

Randomized benchmarking (RB) is the standard experimental protocol for measuring the average error rate of a gate set. The core idea: apply a random sequence of m Clifford gates, then apply the inverse of the entire sequence, and measure whether the qubit returns to its initial state. With depolarizing noise of strength p per gate, the survival probability decays as (1 - p)^m.

The following code implements a simplified single-qubit RB experiment in Cirq:

import cirq
import numpy as np
from scipy.optimize import curve_fit

def random_clifford_circuit(qubit, num_cliffords, rng):
    """Generate a random single-qubit Clifford circuit and its inverse.

    The single-qubit Clifford group has 24 elements. We sample from
    a generating set {I, H, S, H.S, S.H, H.S.H, ...} by composing
    random generators.
    """
    # Generators of the single-qubit Clifford group
    generators = [cirq.I, cirq.H, cirq.S, cirq.X, cirq.Y, cirq.Z]

    ops = []
    composite_unitary = np.eye(2, dtype=complex)

    for _ in range(num_cliffords):
        # Pick a random sequence of 1-3 generators to form a Clifford
        n_gens = rng.integers(1, 4)
        for _ in range(n_gens):
            gate = generators[rng.integers(len(generators))]
            ops.append(gate(qubit))
            composite_unitary = cirq.unitary(gate) @ composite_unitary

    # Compute the inverse unitary
    inverse_unitary = composite_unitary.conj().T

    # Convert inverse to a gate
    inverse_gate = cirq.MatrixGate(inverse_unitary)
    ops.append(inverse_gate(qubit))

    return ops

def run_rb(p_depolarize: float, sequence_lengths: list, 
           num_samples: int = 50):
    """Run randomized benchmarking with depolarizing noise.

    Returns average survival probability for each sequence length.
    """
    qubit = cirq.LineQubit(0)
    sim = cirq.DensityMatrixSimulator()
    rng = np.random.default_rng(42)

    survival_probs = []

    for m in sequence_lengths:
        survivals = []
        for _ in range(num_samples):
            # Build random Clifford sequence with inverse
            ops = random_clifford_circuit(qubit, m, rng)

            # Create noisy circuit: depolarize after each operation
            circuit = cirq.Circuit()
            for op in ops:
                circuit.append(op)
                circuit.append(cirq.depolarize(p_depolarize)(qubit))

            result = sim.simulate(circuit)
            dm = result.final_density_matrix
            # Survival probability: probability of measuring |0>
            survivals.append(float(np.real(dm[0, 0])))

        survival_probs.append(np.mean(survivals))

    return survival_probs

# Run RB experiment
p_true = 0.005  # True depolarizing error rate
sequence_lengths = [1, 5, 10, 20, 50, 100]

print(f"True depolarizing error rate: {p_true}")
print(f"{'m (length)':>12} {'Survival prob':>15}")
print("-" * 29)

survival_data = run_rb(p_true, sequence_lengths, num_samples=50)
for m, surv in zip(sequence_lengths, survival_data):
    print(f"{m:>12d} {surv:>15.4f}")

# Fit exponential decay: P(m) = A * (1-p_fit)^m + B
# For single-qubit depolarizing: A=0.5, B=0.5, decay rate = 1-p
def rb_model(m, A, p_fit, B):
    return A * (1 - p_fit) ** m + B

try:
    popt, pcov = curve_fit(
        rb_model, sequence_lengths, survival_data,
        p0=[0.5, p_true, 0.5], bounds=([0, 0, 0], [1, 1, 1])
    )
    p_fitted = popt[1]
    print(f"\nFitted error rate per Clifford: {p_fitted:.5f}")
    print(f"True error rate:                {p_true:.5f}")
    print(f"Relative difference:            {abs(p_fitted - p_true)/p_true*100:.1f}%")
except RuntimeError:
    print("\nCurve fitting did not converge. Try more samples or sequence lengths.")

The fitted error rate should be close to the true depolarizing rate, though not exact because each “Clifford” in this simplified implementation consists of multiple generator gates, each of which accumulates noise. In a full RB implementation, you would compile each random Clifford into an optimal single gate, so the fitted rate reflects the error per compiled Clifford rather than per generator.

Process Tomography on a Noisy CNOT

Quantum process tomography (QPT) fully characterizes a noisy quantum channel by probing it with a complete set of input states and measuring the outputs. For a two-qubit gate, you need 16 input states (tensor products of {|0>, |1>, |+>, |+i>}) and full state tomography on each output.

The following code performs QPT on a CNOT gate with depolarizing noise and reconstructs the process in the Pauli transfer matrix representation:

import cirq
import numpy as np
from itertools import product

def prepare_state(qubit, state_label):
    """Return operations to prepare a qubit in the specified state."""
    if state_label == '0':
        return []  # Already |0>
    elif state_label == '1':
        return [cirq.X(qubit)]
    elif state_label == '+':
        return [cirq.H(qubit)]
    elif state_label == '+i':
        return [cirq.H(qubit), cirq.S(qubit)]
    raise ValueError(f"Unknown state: {state_label}")

def run_process_tomography(gate_ops, qubits, noise_ops=None):
    """Perform process tomography on a two-qubit gate.

    Args:
        gate_ops: List of operations implementing the gate.
        qubits: The two qubits.
        noise_ops: Optional noise operations to insert after the gate.

    Returns:
        Dictionary mapping (input_label, output_pauli) to expectation value.
    """
    sim = cirq.DensityMatrixSimulator()
    input_states = ['0', '1', '+', '+i']
    pauli_ops = {
        'II': np.eye(4),
        'IZ': np.kron(np.eye(2), np.array([[1,0],[0,-1]])),
        'ZI': np.kron(np.array([[1,0],[0,-1]]), np.eye(2)),
        'ZZ': np.kron(np.array([[1,0],[0,-1]]), np.array([[1,0],[0,-1]])),
        'XX': np.kron(np.array([[0,1],[1,0]]), np.array([[0,1],[1,0]])),
    }

    results = {}
    for s0, s1 in product(input_states, repeat=2):
        # Build circuit: prepare input, apply gate, optionally add noise
        circuit = cirq.Circuit()
        circuit.append(prepare_state(qubits[0], s0))
        circuit.append(prepare_state(qubits[1], s1))
        circuit.append(gate_ops)
        if noise_ops:
            circuit.append(noise_ops)

        dm = sim.simulate(circuit).final_density_matrix

        # Measure Pauli expectation values on output
        for pauli_label, pauli_matrix in pauli_ops.items():
            exp_val = float(np.real(np.trace(pauli_matrix @ dm)))
            results[(f"{s0}{s1}", pauli_label)] = exp_val

    return results

q0, q1 = cirq.LineQubit.range(2)

# Ideal CNOT
ideal_results = run_process_tomography(
    gate_ops=[cirq.CNOT(q0, q1)],
    qubits=[q0, q1],
)

# Noisy CNOT with 2% depolarizing
noisy_results = run_process_tomography(
    gate_ops=[cirq.CNOT(q0, q1)],
    qubits=[q0, q1],
    noise_ops=[cirq.depolarize(0.02)(q0), cirq.depolarize(0.02)(q1)],
)

# Compare a subset of Pauli expectation values
print(f"{'Input':>6} {'Pauli':>6} {'Ideal':>8} {'Noisy':>8} {'Diff':>8}")
print("-" * 40)
for (inp, pauli), ideal_val in sorted(ideal_results.items()):
    noisy_val = noisy_results[(inp, pauli)]
    if abs(ideal_val) > 0.01:  # Only show nonzero entries
        diff = abs(ideal_val - noisy_val)
        print(f"{inp:>6} {pauli:>6} {ideal_val:>8.3f} {noisy_val:>8.3f} {diff:>8.4f}")

For the ideal CNOT, all Pauli expectation values are exactly +1 or -1 (or 0). For the noisy CNOT, expectation values are attenuated toward zero. The degree of attenuation depends on the noise strength and the Pauli observable. This data forms the basis for reconstructing the full chi matrix (process matrix) of the noisy channel, though the full reconstruction requires solving a linear system with 256 unknowns for a two-qubit process.

Comparison to Qiskit and Mitiq Integration

Cirq makes noise explicit inside the circuit. Qiskit attaches a NoiseModel object at transpile time and applies it during AerSimulator execution. Both approaches produce equivalent density matrix trajectories; Cirq’s is more transparent for debugging, Qiskit’s scales further with GPU-accelerated Aer.

For error mitigation, Mitiq works directly with Cirq executors. Zero-noise extrapolation (ZNE) runs the circuit at scaled noise levels and extrapolates to the zero-noise limit:

import mitiq, cirq, numpy as np

def noisy_executor(circuit: cirq.Circuit) -> float:
    noisy = circuit + cirq.Circuit(cirq.depolarize(0.01)(q) for q in circuit.all_qubits())
    dm = cirq.DensityMatrixSimulator().simulate(noisy).final_density_matrix
    return float(np.real(dm[0, 0] + dm[1, 1] - dm[2, 2] - dm[3, 3]))

q0, q1 = cirq.LineQubit.range(2)
circuit = cirq.Circuit([cirq.H(q0), cirq.CNOT(q0, q1)])
mitigated = mitiq.zne.execute_with_zne(circuit, noisy_executor)
print(f"Mitigated expectation value: {mitigated:.4f}")

Mitiq PEC: Probabilistic Error Cancellation

While ZNE extrapolates to zero noise by running at multiple noise levels, Probabilistic Error Cancellation (PEC) takes a different approach: it represents the ideal (noiseless) gate as a linear combination of noisy implementable operations, then samples from that combination to produce an unbiased estimate of the noiseless result.

PEC requires knowing the noise model. You decompose each noisy gate into a quasi-probability distribution over implementable operations (the noisy gate itself plus Pauli corrections). The cost is that PEC requires more circuit executions than ZNE (the sampling overhead grows exponentially with circuit depth), but it produces an unbiased estimator.

import mitiq
import cirq
import numpy as np
from mitiq import pec

def noisy_executor(circuit: cirq.Circuit) -> float:
    """Execute a circuit with 1% depolarizing noise on all qubits."""
    noise_level = 0.01
    noisy_circuit = cirq.Circuit()
    for moment in circuit.moments:
        noisy_circuit.append(moment)
        for op in moment.operations:
            for q in op.qubits:
                noisy_circuit.append(cirq.depolarize(noise_level)(q))
    dm = cirq.DensityMatrixSimulator().simulate(noisy_circuit).final_density_matrix
    # Expectation value of ZZ observable
    zz = np.kron(np.array([[1,0],[0,-1]]), np.array([[1,0],[0,-1]]))
    return float(np.real(np.trace(zz @ dm)))

q0, q1 = cirq.LineQubit.range(2)
circuit = cirq.Circuit([cirq.H(q0), cirq.CNOT(q0, q1)])

# Build PEC representations for each gate at the known noise level
reps = pec.represent_operations_in_circuit_with_local_depolarizing_noise(
    ideal_circuit=circuit,
    noise_level=0.01,
)

# Run PEC
pec_result = pec.execute_with_pec(
    circuit,
    noisy_executor,
    representations=reps,
    num_samples=500,
    random_state=42,
)

# Compare all approaches
ideal_val = 1.0  # ZZ expectation for ideal Bell state |00>+|11> is +1
noisy_val = noisy_executor(circuit)
zne_val = mitiq.zne.execute_with_zne(circuit, noisy_executor)

print(f"Ideal ZZ expectation:     {ideal_val:.4f}")
print(f"Raw noisy result:         {noisy_val:.4f}")
print(f"ZNE mitigated result:     {zne_val:.4f}")
print(f"PEC mitigated result:     {pec_result:.4f}")
print(f"\nError (noisy):  {abs(ideal_val - noisy_val):.4f}")
print(f"Error (ZNE):    {abs(ideal_val - zne_val):.4f}")
print(f"Error (PEC):    {abs(ideal_val - pec_result):.4f}")

PEC generally produces a result closer to the ideal value than ZNE, but at the cost of more circuit executions. The sampling overhead for PEC scales as gamma^2 where gamma is the one-norm of the quasi-probability distribution, which grows exponentially with the number of noisy operations. For shallow circuits with moderate noise, PEC is practical. For deep circuits, ZNE may be the only feasible option.

Common Mistakes in Cirq Noise Modeling

Noise modeling has several subtle pitfalls that can silently produce wrong results. Here are the five most common mistakes.

1. Using cirq.Simulator instead of cirq.DensityMatrixSimulator for noisy circuits.

cirq.Simulator is a statevector simulator. It tracks pure states and cannot represent mixed states. If you add noise channels to a circuit and simulate with cirq.Simulator, the noise channels are silently ignored or produce incorrect results. Always use cirq.DensityMatrixSimulator when your circuit contains noise channels.

import cirq

q = cirq.LineQubit(0)
circuit = cirq.Circuit([
    cirq.X(q),
    cirq.depolarize(0.5)(q),  # Heavy noise: 50% depolarizing
])

# WRONG: statevector simulator does not handle noise correctly
sv_result = cirq.Simulator().simulate(circuit)
print("Statevector P(|1>):", abs(sv_result.final_state_vector[1])**2)

# CORRECT: density matrix simulator handles noise properly
dm_result = cirq.DensityMatrixSimulator().simulate(circuit)
print("DensityMatrix P(|1>):", float(dm_result.final_density_matrix[1, 1].real))

2. Placing noise channels before the gate instead of after.

Physically, noise occurs during or after a gate, not before it. Placing cirq.depolarize(p) before a Hadamard and then applying H models a different error channel than the intended one. The depolarized state gets rotated by H, which produces a different output than depolarizing the already-rotated state. Always insert noise after the gate it models.

import cirq
import numpy as np

q = cirq.LineQubit(0)
sim = cirq.DensityMatrixSimulator()

# WRONG: noise before gate
wrong = cirq.Circuit([cirq.depolarize(0.1)(q), cirq.H(q)])
# CORRECT: noise after gate
right = cirq.Circuit([cirq.H(q), cirq.depolarize(0.1)(q)])

dm_wrong = sim.simulate(wrong).final_density_matrix
dm_right = sim.simulate(right).final_density_matrix

print("Wrong (noise before H), diagonal:", np.real(np.diag(dm_wrong)).round(4))
print("Right (noise after H), diagonal: ", np.real(np.diag(dm_right)).round(4))
# These produce different density matrices

3. Confusing cirq.depolarize(p) parameter semantics.

In Cirq, cirq.depolarize(p) means the total probability of any Pauli error is p. Each of X, Y, Z occurs with probability p/3, and the identity (no error) occurs with probability 1-p. If your hardware reports a 1% gate error rate, use cirq.depolarize(0.01). Do not use p=0.03 thinking each Pauli should have probability 0.01; that would give a 3% total error rate.

4. Forgetting that DensityMatrixSimulator memory scales as 4^n.

The density matrix for n qubits is a 2^n x 2^n complex matrix, requiring 4^n complex numbers (16 bytes each for complex128). Here is the memory cost:

Qubits	Matrix size	Memory
10	1024 x 1024	16 MB
12	4096 x 4096	256 MB
14	16384 x 16384	4 GB
15	32768 x 32768	8 GB
16	65536 x 65536	32 GB

Plan your simulation qubit counts accordingly. For larger systems, consider using cirq.DensityMatrixSimulator with a reduced qubit count, or switch to Monte Carlo trajectory simulation (applying random Kraus operators to statevectors and averaging over many trajectories).

5. Not modeling measurement noise separately from gate noise.

Adding cirq.depolarize(p) after every gate handles gate errors, but measurement errors are a distinct physical process. Real devices have readout error rates that differ from gate error rates, and the errors are typically asymmetric (higher probability of reading 0 when the state is 1 than vice versa). If you only model gate noise and ignore readout error, your simulation underestimates the total error, particularly for algorithms where the output is a single bitstring read from many qubits. Always add a separate readout error model using bit-flip channels before measurement operations, as shown in the Readout Error Modeling section above.

Noise modeling in Cirq is explicit and composable. Channels are gates; they appear in the circuit and can be inspected, moved, and overridden. Start with depolarizing noise for a first pass, then switch to amplitude damping, thermal relaxation, or a calibration-based model when you need accuracy that matches a specific device’s characterization data.