Simulating Superconducting Qubit Physics in Cirq

Google Sycamore and Willow Architecture

Google’s Sycamore and Willow processors arrange transmon qubits on a 2D square grid. Each qubit connects only to its four nearest neighbors (up, down, left, right). No qubit can interact directly with a non-adjacent qubit, so circuits must route two-qubit operations through a sequence of nearby gates.

The hardware-native two-qubit gate is the fSim (fermionic simulation) gate, which implements a specific combination of iSWAP-like and controlled-phase interaction. Cirq models all of this: the grid layout, the connectivity constraints, the fSim gate, and the noise characteristics of real superconducting hardware.

Transmon Physics Background

A transmon qubit is a Josephson junction shunted by a large capacitance, forming an anharmonic oscillator. The Josephson junction provides a nonlinear inductance that creates unequally spaced energy levels. The qubit is encoded in the two lowest energy levels: the ground state |0> and the first excited state |1>.

The critical property that makes the transmon usable as a qubit is anharmonicity. The |0> to |1> transition frequency differs from the |1> to |2> transition frequency by approximately 200 MHz (the anharmonicity). This separation allows microwave pulses tuned to the |0>-|1> frequency to drive transitions between those two levels without accidentally exciting the qubit into |2> or higher states. Without sufficient anharmonicity, you have a harmonic oscillator, and selective addressing of the computational subspace becomes impossible.

Single-qubit gates are implemented by applying calibrated microwave pulses at the qubit’s resonance frequency (typically 4 to 6 GHz). The pulse amplitude, duration, and phase determine the rotation axis and angle on the Bloch sphere. A pi-pulse around X produces an X gate; a pi/2-pulse produces a sqrt(X) gate. Phase shifts in the drive signal rotate the axis of rotation, enabling arbitrary single-qubit unitaries.

Two-qubit gates arise from capacitive coupling between adjacent transmons. When two qubits are brought into resonance (or near-resonance) through frequency tuning, the XX + YY exchange interaction produces the iSWAP-like component of fSim, while a conditional phase accumulation produces the CZ-like component. The total interaction is parameterized as fSim(theta, phi), where theta controls the swap angle and phi controls the conditional phase.

The 2D grid layout is not arbitrary. It minimizes parasitic cross-talk between non-adjacent qubits while maintaining sufficient connectivity for universal quantum computation. Each qubit couples only to its immediate neighbors, and the coupling strength to next-nearest neighbors is suppressed by several orders of magnitude. This nearest-neighbor constraint is a fundamental design choice that trades connectivity for noise isolation.

GridQubit: Addressing Qubits on the Grid

import cirq

# Create individual qubits by grid position
q00 = cirq.GridQubit(0, 0)
q01 = cirq.GridQubit(0, 1)
q10 = cirq.GridQubit(1, 0)

# Create a 3x3 patch of the grid
grid = cirq.GridQubit.rect(3, 3)
print(grid)
# [GridQubit(0,0), GridQubit(0,1), ..., GridQubit(2,2)]

# Check adjacency
print(cirq.GridQubit(0, 0).is_adjacent(cirq.GridQubit(0, 1)))  # True
print(cirq.GridQubit(0, 0).is_adjacent(cirq.GridQubit(1, 1)))  # False (diagonal)

Only adjacent qubits can have two-qubit operations. Attempting to apply a two-qubit gate between non-adjacent qubits on real hardware requires SWAP routing; Cirq’s compiler handles this automatically.

GridQubit Topology Operations

The GridQubit.rect(rows, cols) factory creates a rectangular patch of qubits. You can also specify an offset to place the patch at an arbitrary location on the grid:

import cirq

# Create a 4x5 grid patch starting at row=0, col=0
large_grid = cirq.GridQubit.rect(4, 5)
print(f"Total qubits: {len(large_grid)}")  # 20

# Create a patch offset to a different region
offset_grid = cirq.GridQubit.rect(2, 3, top=3, left=2)
print(offset_grid)
# [GridQubit(3,2), GridQubit(3,3), GridQubit(3,4),
#  GridQubit(4,2), GridQubit(4,3), GridQubit(4,4)]

# Sort qubits by row then column (default ordering)
qubits = cirq.GridQubit.rect(3, 3)
sorted_qubits = sorted(qubits, key=lambda q: (q.row, q.col))
for q in sorted_qubits:
    print(f"  ({q.row}, {q.col})", end="")
print()
# (0, 0) (0, 1) (0, 2) (1, 0) (1, 1) (1, 2) (2, 0) (2, 1) (2, 2)

For building entangling layers, you need to enumerate all adjacent pairs in a grid. This helper function returns every pair of adjacent GridQubits:

import cirq
from typing import List, Tuple

def get_adjacent_pairs(
    qubits: List[cirq.GridQubit],
) -> List[Tuple[cirq.GridQubit, cirq.GridQubit]]:
    """Return all pairs (q0, q1) where q0 and q1 are adjacent GridQubits.

    Each pair appears once with q0 < q1 in the default GridQubit ordering.
    """
    qubit_set = set(qubits)
    pairs = []
    for q in qubits:
        # Check right neighbor and down neighbor to avoid duplicates
        right = cirq.GridQubit(q.row, q.col + 1)
        down = cirq.GridQubit(q.row + 1, q.col)
        if right in qubit_set:
            pairs.append((q, right))
        if down in qubit_set:
            pairs.append((q, down))
    return pairs


# Example: find all adjacent pairs in a 3x3 grid
grid = cirq.GridQubit.rect(3, 3)
pairs = get_adjacent_pairs(grid)
print(f"Adjacent pairs in 3x3 grid: {len(pairs)}")  # 12
for a, b in pairs:
    print(f"  ({a.row},{a.col}) -- ({b.row},{b.col})")

You can also check whether a set of qubits forms a connected subgraph. This matters when selecting qubit subsets for algorithms that require all qubits to be reachable from each other:

import cirq
from collections import deque
from typing import Set

def is_connected(qubits: Set[cirq.GridQubit]) -> bool:
    """Check if a set of GridQubits forms a connected subgraph."""
    if not qubits:
        return True
    visited = set()
    queue = deque([next(iter(qubits))])
    while queue:
        q = queue.popleft()
        if q in visited:
            continue
        visited.add(q)
        for neighbor in [
            cirq.GridQubit(q.row - 1, q.col),
            cirq.GridQubit(q.row + 1, q.col),
            cirq.GridQubit(q.row, q.col - 1),
            cirq.GridQubit(q.row, q.col + 1),
        ]:
            if neighbor in qubits and neighbor not in visited:
                queue.append(neighbor)
    return visited == qubits


# Connected subset
subset_a = {cirq.GridQubit(0, 0), cirq.GridQubit(0, 1), cirq.GridQubit(1, 1)}
print(f"Connected: {is_connected(subset_a)}")  # True

# Disconnected subset (gap in the middle)
subset_b = {cirq.GridQubit(0, 0), cirq.GridQubit(2, 2)}
print(f"Connected: {is_connected(subset_b)}")  # False

The fSim Gate

The fSim gate implements the unitary:

fSim(theta, phi) = [[1, 0,          0,          0],
                     [0, cos(theta), -i*sin(theta), 0],
                     [0, -i*sin(theta), cos(theta), 0],
                     [0, 0,          0,          exp(-i*phi)]]

Physical Origin

The fSim gate arises naturally from the XX + YY exchange interaction between capacitively coupled transmons. When two transmon qubits are coupled, the interaction Hamiltonian in the single-excitation subspace takes the form:

H_int = g * (|01><10| + |10><01|) + delta * |11><11|

The first term (proportional to coupling strength g) drives population exchange between |01> and |10>, producing the theta (iSWAP-like) component. The second term (conditional phase shift delta) accumulates a phase on |11>, producing the phi (CZ-like) component. By controlling the interaction time and detuning, the hardware tunes theta and phi independently.

Sycamore Gate Unitary

At theta=pi/2, phi=pi/6 the gate produces the Sycamore gate used in Google’s quantum supremacy experiments:

import numpy as np
import cirq

# The Sycamore gate: theta=pi/2, phi=pi/6
sycamore_gate = cirq.FSimGate(theta=np.pi / 2, phi=np.pi / 6)

q0, q1 = cirq.GridQubit(0, 0), cirq.GridQubit(0, 1)
circuit = cirq.Circuit(sycamore_gate(q0, q1))
print(circuit)

# Print the full 4x4 unitary matrix
U = cirq.unitary(sycamore_gate)
print("\nSycamore gate unitary:")
np.set_printoptions(precision=4, suppress=True)
print(U)
# Expected:
# [[1.+0.j     0.+0.j     0.+0.j     0.+0.j    ]
#  [0.+0.j     0.+0.j     0.-1.j     0.+0.j    ]
#  [0.+0.j     0.-1.j     0.+0.j     0.+0.j    ]
#  [0.+0.j     0.+0.j     0.+0.j     0.866-0.5j ]]

# Verify: cos(pi/2) = 0, -i*sin(pi/2) = -i, exp(-i*pi/6) = cos(pi/6) - i*sin(pi/6)
assert np.isclose(U[1, 1], 0.0)  # cos(pi/2)
assert np.isclose(U[1, 2], -1j)  # -i*sin(pi/2)
assert np.isclose(U[3, 3], np.exp(-1j * np.pi / 6))
print("\nAll matrix elements verified.")

Five Important Special Cases

The fSim gate family contains several gates of fundamental importance:

import numpy as np
import cirq

special_cases = {
    "iSWAP":     (np.pi / 2, 0),
    "CZ":        (0, np.pi),
    "SWAP-like": (np.pi / 2, np.pi),
    "Sycamore":  (np.pi / 2, np.pi / 6),
    "Identity":  (0, 0),
}

for name, (theta, phi) in special_cases.items():
    gate = cirq.FSimGate(theta=theta, phi=phi)
    U = cirq.unitary(gate)
    print(f"\n{name} gate: fSim(theta={theta:.4f}, phi={phi:.4f})")
    print(U)

iSWAP (theta=pi/2, phi=0): Swaps the amplitudes of |01> and |10> with a factor of -i. This gate arises from a half-period of the exchange interaction with no conditional phase. It is the natural “free evolution” gate for resonantly coupled transmons.

CZ (theta=0, phi=pi): Applies a phase of -1 to |11> and acts as identity on all other basis states. There is no population exchange. CZ is widely used as a universal entangling gate in many quantum computing architectures.

SWAP-like (theta=pi/2, phi=pi): Combines a full iSWAP swap with a CZ phase. This is equivalent to SWAP up to single-qubit phases. The population exchange swaps |01> and |10>, and the conditional phase flips the sign of |11>.

Sycamore (theta=pi/2, phi=pi/6): The specific gate calibrated on Google’s Sycamore processor. It provides maximal entangling power (theta=pi/2) with a moderate conditional phase. The pi/6 phase was chosen because it matches the natural interaction parameters of the hardware, minimizing calibration overhead.

Identity (theta=0, phi=0): No swap, no phase. The gate does nothing. This is useful as a baseline for verifying gate decompositions and as a limiting case in parameterized circuits.

Parameterized fSim Gates for Variational Algorithms

Variational algorithms like VQE and QAOA require circuits with tunable parameters. Cirq supports symbolic parameters in fSim gates through sympy:

import sympy
import numpy as np
import cirq

# Define symbolic parameters
theta = sympy.Symbol('theta')
phi = sympy.Symbol('phi')

# Create a parameterized fSim gate
param_fsim = cirq.FSimGate(theta=theta, phi=phi)

# Build a variational ansatz using parameterized fSim gates
qubits = cirq.GridQubit.rect(1, 4)
circuit = cirq.Circuit()

# Layer of single-qubit rotations
for i, q in enumerate(qubits):
    circuit.append(cirq.ry(sympy.Symbol(f'ry_{i}')).on(q))

# Layer of parameterized fSim gates on adjacent pairs
for i in range(len(qubits) - 1):
    q0, q1 = qubits[i], qubits[i + 1]
    if q0.is_adjacent(q1):
        fsim_i = cirq.FSimGate(
            theta=sympy.Symbol(f'theta_{i}'),
            phi=sympy.Symbol(f'phi_{i}'),
        )
        circuit.append(fsim_i(q0, q1))

# Second layer of single-qubit rotations
for i, q in enumerate(qubits):
    circuit.append(cirq.ry(sympy.Symbol(f'ry2_{i}')).on(q))

circuit.append(cirq.measure(*qubits, key='result'))
print(circuit)

# Bind parameters to specific values using ParamResolver
resolver = cirq.ParamResolver({
    'ry_0': 0.3, 'ry_1': 0.5, 'ry_2': 0.7, 'ry_3': 0.1,
    'theta_0': np.pi / 4, 'phi_0': np.pi / 8,
    'theta_1': np.pi / 3, 'phi_1': np.pi / 5,
    'theta_2': np.pi / 6, 'phi_2': np.pi / 4,
    'ry2_0': 0.2, 'ry2_1': 0.4, 'ry2_2': 0.6, 'ry2_3': 0.8,
})

# Simulate with bound parameters
sim = cirq.Simulator()
resolved_circuit = cirq.resolve_parameters(circuit, resolver)
result = sim.run(resolved_circuit, repetitions=1000)
print(result.histogram(key='result'))

This pattern is directly relevant to running VQE on Google hardware: the classical optimizer updates the theta and phi values, and each evaluation resolves the symbolic parameters before execution. On real hardware, you submit the parameterized circuit once and sweep over parameter values efficiently using cirq.Sweepable.

Building a Random Circuit for XEB

Cross-entropy benchmarking (XEB) is the protocol Google used to demonstrate quantum supremacy. It generates random circuits, samples the output, and compares the distribution to the ideal (simulated) one. A high XEB fidelity means the hardware closely matches the ideal quantum computer.

XEB Protocol Step by Step

The XEB protocol follows a precise structure:

1. Generate a random circuit. Each cycle consists of a layer of random single-qubit gates (chosen from a fixed gate set, typically {sqrt(X), sqrt(Y), sqrt(W)}) followed by a layer of two-qubit fSim gates. The two-qubit gates are applied in a specific spatial pattern that alternates between four configurations (A, B, C, D layers) on the 2D grid. Each pattern covers a different subset of adjacent pairs, ensuring all qubits participate in entangling operations across cycles.

2. Simulate the ideal statevector. Using a classical simulator, compute the exact output state |psi> of the circuit. Extract the probability p(s) = |<s|psi>|^2 for every bitstring s. For n qubits, this requires storing 2^n amplitudes, which is the computational bottleneck that makes XEB hard classically for large n.

3. Run on hardware (or noisy simulator). Execute the circuit many times (typically thousands of shots) and collect the sampled bitstrings {s_1, s_2, …, s_K}.

4. Compute the linear cross-entropy benchmark. The XEB fidelity is:

   F_XEB = 2^n * mean(p_ideal(s_i)) - 1

   where the mean is taken over the K sampled bitstrings. For a perfect quantum computer, the sampled bitstrings are drawn from the ideal distribution, and Porter-Thomas statistics give F_XEB = 1. For a completely depolarized device that outputs uniformly random bitstrings, each bitstring samples p_ideal uniformly, giving mean(p_ideal) = 1/2^n, and therefore F_XEB = 0.

The statistical relationship between F_XEB and circuit fidelity is: F_XEB approximates the circuit fidelity (the probability that no error occurred during the entire circuit). This holds because depolarizing errors push the output distribution toward uniform, and F_XEB measures exactly this deviation.

import numpy as np
import cirq

def random_xeb_circuit(qubits, depth, seed=0):
    rng = np.random.default_rng(seed)
    single_qubit_gates = [cirq.X**0.5, cirq.Y**0.5, cirq.T, cirq.H]
    sycamore = cirq.FSimGate(np.pi / 2, np.pi / 6)

    moments = []
    for d in range(depth):
        # Single-qubit random gates on all qubits
        sq_ops = [rng.choice(single_qubit_gates)(q) for q in qubits]
        moments.append(cirq.Moment(sq_ops))

        # Two-qubit fSim on alternating adjacent pairs
        if d % 2 == 0:
            pairs = [(qubits[i], qubits[i+1])
                     for i in range(0, len(qubits)-1, 2)
                     if qubits[i].is_adjacent(qubits[i+1])]
        else:
            pairs = [(qubits[i], qubits[i+1])
                     for i in range(1, len(qubits)-1, 2)
                     if qubits[i].is_adjacent(qubits[i+1])]
        if pairs:
            tq_ops = [sycamore(a, b) for a, b in pairs]
            moments.append(cirq.Moment(tq_ops))

    moments.append(cirq.Moment([cirq.measure(*qubits, key="m")]))
    return cirq.Circuit(moments)

qubits = cirq.GridQubit.rect(2, 3)
circuit = random_xeb_circuit(qubits, depth=10)
print(circuit)

Computing XEB Fidelity

XEB fidelity compares the measured bitstring probabilities against the ideal probabilities:

def compute_xeb_fidelity(circuit, qubits, n_shots=10000):
    # Ideal simulation (statevector)
    ideal_sim = cirq.Simulator()
    ideal_result = ideal_sim.simulate(circuit[:-1])  # exclude measurement
    ideal_state = ideal_result.final_state_vector
    n = len(qubits)
    ideal_probs = np.abs(ideal_state) ** 2

    # Sampled result (noisy or ideal)
    sampled_result = ideal_sim.run(circuit, repetitions=n_shots)
    bitstrings = sampled_result.measurements["m"]

    # XEB fidelity: F_XEB = 2^n * <p_ideal(s)> - 1
    sampled_ideal_probs = []
    for bits in bitstrings:
        idx = int("".join(str(b) for b in bits), 2)
        sampled_ideal_probs.append(ideal_probs[idx])

    F_xeb = 2**n * np.mean(sampled_ideal_probs) - 1
    return F_xeb

fidelity = compute_xeb_fidelity(circuit, qubits)
print(f"XEB fidelity (ideal sim): {fidelity:.4f}")
# For ideal sampling, F_XEB = 2^n * sum(p^2) - 1, which is > 1 for a non-uniform
# distribution. F_XEB == 1 only on average over many random circuits (Porter-Thomas).
# For a noisy device, F_XEB approaches 0 as the output approaches uniform.

Adding a Noise Model

Real Sycamore-style hardware has T1 around 15 microseconds, T2 around 20 microseconds, and depolarizing error of roughly 0.1% per single-qubit gate and 0.5% per two-qubit gate. Cirq models this with DensityMatrixSimulator:

# Build a composite noise model
T1_ns = 15_000    # 15 microseconds in nanoseconds
T2_ns = 20_000
gate_time_ns = 25

class SycamoreNoise(cirq.NoiseModel):
    def noisy_moment(self, moment, system_qubits):
        noise_ops = []
        # Add depolarizing noise after each gate
        for op in moment.operations:
            n = len(op.qubits)
            if n == 1:
                noise_ops.append(cirq.depolarize(p=0.001).on(op.qubits[0]))
            elif n == 2:
                for q in op.qubits:
                    noise_ops.append(cirq.depolarize(p=0.005).on(q))
        if noise_ops:
            return [moment, cirq.Moment(noise_ops)]
        return [moment]

noisy_sim = cirq.DensityMatrixSimulator(noise=SycamoreNoise())
noisy_result = noisy_sim.run(circuit, repetitions=5000)
print(noisy_result.measurements["m"][:5])

T1 and T2 Noise Model

The depolarizing model above captures incoherent errors but does not faithfully represent the physics of energy relaxation (T1) and dephasing (T2). A more accurate model uses amplitude damping for T1 decay and phase damping for pure dephasing.

For a gate of duration t, the relevant error probabilities are:

• T1 amplitude damping probability: p_ad = 1 - exp(-t / T1)
• Pure dephasing probability: p_pd = 1 - exp(-t / T2) (after subtracting the T1 contribution, the pure dephasing rate is 1/T_phi = 1/T2 - 1/(2*T1))

Cirq provides cirq.AmplitudeDampingChannel(gamma) for T1 decay and cirq.PhaseDampingChannel(gamma) for pure dephasing:

import numpy as np
import cirq

class T1T2NoiseModel(cirq.NoiseModel):
    """Noise model combining amplitude damping, phase damping, and depolarizing error."""

    def __init__(
        self,
        t1_ns: float = 15_000,
        t2_ns: float = 20_000,
        single_qubit_gate_ns: float = 25,
        two_qubit_gate_ns: float = 40,
        depol_1q: float = 0.001,
        depol_2q: float = 0.005,
    ):
        self.t1_ns = t1_ns
        self.t2_ns = t2_ns
        self.single_qubit_gate_ns = single_qubit_gate_ns
        self.two_qubit_gate_ns = two_qubit_gate_ns
        self.depol_1q = depol_1q
        self.depol_2q = depol_2q

    def _noise_for_duration(self, duration_ns: float) -> list:
        """Compute amplitude damping and phase damping channels for a given gate duration."""
        # T1 amplitude damping
        p_ad = 1 - np.exp(-duration_ns / self.t1_ns)
        # Pure dephasing (T_phi contribution only)
        # 1/T_phi = 1/T2 - 1/(2*T1)
        t_phi = 1.0 / (1.0 / self.t2_ns - 1.0 / (2.0 * self.t1_ns))
        p_pd = 1 - np.exp(-duration_ns / t_phi)
        return p_ad, p_pd

    def noisy_moment(self, moment, system_qubits):
        ad_ops = []
        pd_ops = []
        dp_ops = []
        for op in moment.operations:
            n_qubits = len(op.qubits)
            if n_qubits == 1:
                duration = self.single_qubit_gate_ns
                p_ad, p_pd = self._noise_for_duration(duration)
                q = op.qubits[0]
                ad_ops.append(cirq.amplitude_damp(p_ad).on(q))
                pd_ops.append(cirq.phase_damp(p_pd).on(q))
                dp_ops.append(cirq.depolarize(self.depol_1q).on(q))
            elif n_qubits == 2:
                duration = self.two_qubit_gate_ns
                p_ad, p_pd = self._noise_for_duration(duration)
                for q in op.qubits:
                    ad_ops.append(cirq.amplitude_damp(p_ad).on(q))
                    pd_ops.append(cirq.phase_damp(p_pd).on(q))
                    dp_ops.append(cirq.depolarize(self.depol_2q).on(q))
        # Each channel type must be in its own Moment: a Moment cannot hold
        # multiple operations on the same qubit.
        result = [moment]
        if ad_ops:
            result.append(cirq.Moment(ad_ops))
        if pd_ops:
            result.append(cirq.Moment(pd_ops))
        if dp_ops:
            result.append(cirq.Moment(dp_ops))
        return result


# Instantiate the T1/T2 noise model with Sycamore-class parameters
sycamore_noise = T1T2NoiseModel(
    t1_ns=15_000,
    t2_ns=20_000,
    single_qubit_gate_ns=25,
    two_qubit_gate_ns=40,
    depol_1q=0.001,
    depol_2q=0.005,
)

noisy_sim_t1t2 = cirq.DensityMatrixSimulator(noise=sycamore_noise)

# Run the XEB circuit with the T1/T2 noise model
qubits = cirq.GridQubit.rect(2, 3)
circuit = random_xeb_circuit(qubits, depth=5)
result = noisy_sim_t1t2.run(circuit, repetitions=5000)
print("Sample measurements with T1/T2 noise:")
print(result.measurements["m"][:5])

Readout Error Modeling

Gate errors corrupt the quantum state during circuit execution, but readout errors corrupt the classical measurement outcomes. On Google hardware, typical readout fidelity is 95 to 98%, meaning that 2 to 5% of the time a qubit in |0> is read as 1 (or vice versa).

In Cirq, you model readout noise by inserting cirq.BitFlipChannel immediately before the measurement moment:

import cirq
import numpy as np

class ReadoutNoiseModel(cirq.NoiseModel):
    """Adds asymmetric readout noise before measurement operations."""

    def __init__(self, p_readout: float = 0.02):
        self.p_readout = p_readout

    def noisy_moment(self, moment, system_qubits):
        # Check if this moment contains measurements
        has_measurement = any(cirq.is_measurement(op) for op in moment.operations)
        if has_measurement:
            # Insert bit-flip noise before the measurement
            readout_ops = []
            for op in moment.operations:
                if cirq.is_measurement(op):
                    for q in op.qubits:
                        readout_ops.append(
                            cirq.bit_flip(self.p_readout).on(q)
                        )
            return [cirq.Moment(readout_ops), moment]
        return [moment]


# Combine gate noise and readout noise
class FullNoiseModel(cirq.NoiseModel):
    """Combines T1/T2 gate noise with readout error."""

    def __init__(self, gate_noise: cirq.NoiseModel, readout_noise: cirq.NoiseModel):
        self.gate_noise = gate_noise
        self.readout_noise = readout_noise

    def noisy_moment(self, moment, system_qubits):
        # Apply readout noise to measurement moments
        has_measurement = any(cirq.is_measurement(op) for op in moment.operations)
        if has_measurement:
            return self.readout_noise.noisy_moment(moment, system_qubits)
        # Apply gate noise to all other moments
        return self.gate_noise.noisy_moment(moment, system_qubits)


full_noise = FullNoiseModel(
    gate_noise=T1T2NoiseModel(),
    readout_noise=ReadoutNoiseModel(p_readout=0.03),
)

noisy_sim_full = cirq.DensityMatrixSimulator(noise=full_noise)

# Demonstrate the effect of readout noise on a simple circuit
q = cirq.GridQubit(0, 0)
simple_circuit = cirq.Circuit([
    cirq.X(q),  # Prepare |1>
    cirq.measure(q, key='m'),
])

result = noisy_sim_full.run(simple_circuit, repetitions=1000)
counts = result.histogram(key='m')
print(f"Prepared |1>, measured outcomes: {dict(counts)}")
print(f"Readout error rate: {counts.get(0, 0) / 1000:.3f}")
# Expect ~3% of outcomes to incorrectly read as 0

The distinction matters for error mitigation. Gate errors accumulate with circuit depth and affect the quantum state itself. Readout errors are independent of circuit depth and only corrupt the final classical outcomes. Readout error mitigation can be applied in post-processing using a confusion matrix, while gate errors require techniques like zero-noise extrapolation or probabilistic error cancellation.

Circuit Depth vs Fidelity Analysis

A critical property of NISQ circuits is that fidelity degrades exponentially with the number of noisy gates. For a circuit with n_2q two-qubit gates, each with error rate epsilon, the expected circuit fidelity is approximately:

F = (1 - epsilon)^n_2q

This exponential decay places a hard limit on useful circuit depth. The following code measures XEB fidelity at increasing circuit depths and fits the decay to extract the effective per-gate error rate:

import numpy as np
import cirq

def compute_xeb_fidelity_v2(circuit, qubits, simulator, n_shots=10000):
    ideal_sim = cirq.Simulator()
    ideal_result = ideal_sim.simulate(circuit[:-1])
    ideal_probs = np.abs(ideal_result.final_state_vector) ** 2

    sampled_result = simulator.run(circuit, repetitions=n_shots)
    bitstrings = sampled_result.measurements["m"]
    n = len(qubits)
    sampled_ideal_probs = [
        ideal_probs[int("".join(str(b) for b in bits), 2)]
        for bits in bitstrings
    ]
    return 2**n * np.mean(sampled_ideal_probs) - 1


def count_two_qubit_gates(circuit):
    """Count the number of two-qubit gates in a circuit (excluding measurements)."""
    count = 0
    for moment in circuit:
        for op in moment.operations:
            if len(op.qubits) == 2 and not cirq.is_measurement(op):
                count += 1
    return count


qubits = cirq.GridQubit.rect(2, 3)
noisy_sim = cirq.DensityMatrixSimulator(noise=T1T2NoiseModel())

depths = [1, 3, 5, 10, 20]
fidelities = []
gate_counts = []

for depth in depths:
    circ = random_xeb_circuit(qubits, depth=depth, seed=42)
    n_2q = count_two_qubit_gates(circ)
    fxeb = compute_xeb_fidelity_v2(circ, qubits, noisy_sim, n_shots=5000)
    fidelities.append(fxeb)
    gate_counts.append(n_2q)
    print(f"Depth {depth:2d}: {n_2q:3d} two-qubit gates, F_XEB = {fxeb:.4f}")

# Fit exponential decay: F = (1 - epsilon)^n_2q
# Take log: log(F) = n_2q * log(1 - epsilon)
# Linear fit in log space
fidelities_arr = np.array(fidelities)
gate_counts_arr = np.array(gate_counts)

# Filter out non-positive fidelities for log fitting
mask = fidelities_arr > 0
if np.sum(mask) >= 2:
    log_f = np.log(fidelities_arr[mask])
    n_2q_fit = gate_counts_arr[mask]
    # Linear regression: log(F) = slope * n_2q
    slope = np.polyfit(n_2q_fit, log_f, 1)[0]
    epsilon_fit = 1 - np.exp(slope)
    print(f"\nFitted per-gate error rate: {epsilon_fit:.4f}")
    print(f"This means each two-qubit gate reduces fidelity by ~{epsilon_fit*100:.2f}%")

The fitted per-gate error rate should be consistent with the depolarizing parameter in the noise model (approximately 0.5% for Sycamore-class noise). At depth 20 with 6 qubits, you see significant fidelity loss, matching the behavior observed on real hardware.

Qubit Routing with SWAP Networks

When your algorithm requires a two-qubit gate between non-adjacent qubits, SWAP gates must be inserted to move qubit states to adjacent positions. Each SWAP gate decomposes into three fSim gates, so routing adds significant overhead.

Manual SWAP Decomposition

A SWAP between adjacent qubits decomposes into three iSWAP gates (a special case of fSim):

import cirq
import numpy as np

q0 = cirq.GridQubit(0, 0)
q1 = cirq.GridQubit(0, 1)

# SWAP = three iSWAP gates with single-qubit corrections
# In practice, SWAP ~ fSim(pi/2, pi) which is a native-ish operation
swap_circuit = cirq.Circuit([
    cirq.SWAP(q0, q1),
])

# Decompose into native gates
print("SWAP decomposed:")
decomposed = cirq.optimize_for_target_gateset(
    swap_circuit,
    gateset=cirq.CZTargetGateset(),
)
print(decomposed)

Routing Non-Adjacent Operations

Consider a 3-qubit linear chain where you need a CNOT between qubits at positions (0,0) and (2,0), which are not adjacent. The routing compiler inserts SWAP gates through the intermediate qubit at (1,0):

import cirq
import networkx as nx

q0 = cirq.GridQubit(0, 0)
q1 = cirq.GridQubit(1, 0)
q2 = cirq.GridQubit(2, 0)

# Circuit that requires non-adjacent interaction
original_circuit = cirq.Circuit([
    cirq.H(q0),
    cirq.CNOT(q0, q2),  # q0 and q2 are NOT adjacent
    cirq.measure(q0, q1, q2, key='result'),
])
print("Original circuit:")
print(original_circuit)

# Build the device graph for the 3-qubit linear chain
device_qubits = [q0, q1, q2]
device_graph = nx.Graph()
device_graph.add_nodes_from(device_qubits)
device_graph.add_edge(q0, q1)
device_graph.add_edge(q1, q2)

# Use RouteCQC to insert SWAP gates for the grid topology
router = cirq.RouteCQC(device_graph)
routed_circuit, swap_map, initial_map = router.route_circuit(original_circuit)
print("\nRouted circuit (SWAP gates inserted):")
print(routed_circuit)

Routing is not free. Each inserted SWAP gate adds three two-qubit gates worth of noise. For algorithms with many long-range interactions, routing overhead can dominate the circuit depth. This is why hardware topology matters: a denser connectivity graph requires fewer SWAPs, but at the cost of increased cross-talk.

Willow Processor Improvements Over Sycamore

Google’s Willow processor (2024) represents a significant advance over Sycamore. Willow has 105 qubits (compared to Sycamore’s 53) and substantially improved noise characteristics:

Parameter	Sycamore (2019)	Willow (2024)
Qubit count	53	105
T1 coherence time	~15 microseconds	~100 microseconds
T2 coherence time	~20 microseconds	~80 microseconds
Two-qubit gate error	~0.5%	~0.3%
Single-qubit gate error	~0.1%	~0.05%
Readout error	~3%	~1%

The most significant achievement is demonstrating below-threshold error correction. “Below threshold” means that as you increase the code distance d of a surface code, the logical error rate actually decreases. This is the fundamental requirement for scalable quantum error correction. Above threshold, adding more physical qubits makes things worse (the additional noise from the extra qubits outweighs the error correction benefit). Below threshold, each increase in code distance exponentially suppresses logical errors.

To model Willow-class noise in Cirq, update the noise parameters:

import cirq
import numpy as np

# Willow-class noise model
willow_noise = T1T2NoiseModel(
    t1_ns=100_000,         # 100 microseconds
    t2_ns=80_000,          # 80 microseconds
    single_qubit_gate_ns=25,
    two_qubit_gate_ns=32,  # Slightly faster gates on Willow
    depol_1q=0.0005,       # 0.05% single-qubit error
    depol_2q=0.003,        # 0.3% two-qubit error
)

# Compare Sycamore vs Willow noise on the same circuit
sycamore_noise = T1T2NoiseModel(
    t1_ns=15_000,
    t2_ns=20_000,
    single_qubit_gate_ns=25,
    two_qubit_gate_ns=40,
    depol_1q=0.001,
    depol_2q=0.005,
)

qubits = cirq.GridQubit.rect(2, 3)
circuit = random_xeb_circuit(qubits, depth=10, seed=42)

sim_sycamore = cirq.DensityMatrixSimulator(noise=sycamore_noise)
sim_willow = cirq.DensityMatrixSimulator(noise=willow_noise)

f_syc = compute_xeb_fidelity_v2(circuit, qubits, sim_sycamore, n_shots=5000)
f_wil = compute_xeb_fidelity_v2(circuit, qubits, sim_willow, n_shots=5000)

print(f"Sycamore-class XEB fidelity: {f_syc:.4f}")
print(f"Willow-class XEB fidelity:   {f_wil:.4f}")
print(f"Improvement factor:          {f_wil / max(f_syc, 1e-6):.2f}x")

The longer coherence times on Willow mean that circuits of the same depth suffer less from T1/T2 decay, and the lower gate error rates compound to give significantly higher circuit fidelity.

Comparing Ideal vs Noisy Circuit Output

To build concrete intuition about noise effects, consider a 3-qubit GHZ state circuit. The GHZ state is (|000> + |111>) / sqrt(2), which is straightforward to verify:

import numpy as np
import cirq

q0, q1, q2 = cirq.GridQubit(0, 0), cirq.GridQubit(0, 1), cirq.GridQubit(0, 2)

# Build GHZ circuit
ghz_circuit = cirq.Circuit([
    cirq.H(q0),
    cirq.CNOT(q0, q1),
    cirq.CNOT(q1, q2),
])

# --- Ideal statevector ---
ideal_sim = cirq.Simulator()
ideal_result = ideal_sim.simulate(ghz_circuit)
ideal_state = ideal_result.final_state_vector
print("Ideal statevector:")
for i, amp in enumerate(ideal_state):
    if abs(amp) > 1e-6:
        print(f"  |{i:03b}> : {amp:.4f}")

# Compute fidelity with the target GHZ state
ghz_target = np.zeros(8, dtype=complex)
ghz_target[0] = 1 / np.sqrt(2)  # |000>
ghz_target[7] = 1 / np.sqrt(2)  # |111>

ideal_fidelity = abs(np.vdot(ghz_target, ideal_state)) ** 2
print(f"\nIdeal fidelity with GHZ state: {ideal_fidelity:.6f}")
# Should be 1.0

# --- Noisy density matrix ---
noisy_sim = cirq.DensityMatrixSimulator(noise=T1T2NoiseModel())
noisy_dm_result = noisy_sim.simulate(ghz_circuit)
noisy_dm = noisy_dm_result.final_density_matrix

# Compute fidelity of noisy density matrix with ideal GHZ state
# F = <psi| rho |psi>
noisy_fidelity = np.real(ghz_target.conj() @ noisy_dm @ ghz_target)
print(f"Noisy density matrix fidelity: {noisy_fidelity:.6f}")

# --- Noisy sampling ---
ghz_with_meas = ghz_circuit + cirq.Circuit(
    cirq.measure(q0, q1, q2, key='m')
)
noisy_samples = noisy_sim.run(ghz_with_meas, repetitions=1000)
bitstrings = noisy_samples.measurements['m']

# Count outcomes
from collections import Counter
counts = Counter()
for bits in bitstrings:
    key = "".join(str(b) for b in bits)
    counts[key] += 1

print("\nNoisy sampling results (1000 shots):")
for bitstring in sorted(counts.keys()):
    pct = counts[bitstring] / 1000 * 100
    print(f"  |{bitstring}>: {counts[bitstring]:4d} ({pct:.1f}%)")

# The ideal GHZ state produces only |000> and |111>
# Noise introduces erroneous outcomes like |001>, |010>, |110>
correct_outcomes = counts.get("000", 0) + counts.get("111", 0)
print(f"\nFraction of correct outcomes: {correct_outcomes / 1000:.3f}")

The ideal simulation produces fidelity 1.0 with the GHZ state. The noisy density matrix typically shows fidelity between 0.85 and 0.95 depending on noise parameters. The noisy samples reveal incorrect bitstrings such as |001>, |010>, and |110> appearing at low frequency. These erroneous outcomes arise from bit-flip errors during the CNOT gates and from T1 decay during the circuit.

Validating Against Device Topology

Cirq can validate that a circuit respects the connectivity of a target device. For Google devices, cirq.google provides device objects:

try:
    import cirq_google
    device = cirq_google.Sycamore23
    device.validate_circuit(circuit)
    print("Circuit is valid for Sycamore23")
except Exception as e:
    print(f"Validation error: {e}")

If a gate is placed between non-adjacent qubits, validate_circuit raises an exception listing the offending operation.

Compiling for Sycamore

cirq.optimize_for_target_gateset() with SycamoreTargetGateset takes an arbitrary Cirq circuit and rewrites it using native Sycamore gates and qubit routing:

try:
    import cirq_google
    optimized = cirq.optimize_for_target_gateset(
        circuit, gateset=cirq_google.SycamoreTargetGateset()
    )
    print(f"Original depth: {len(circuit)}")
    print(f"Optimized depth: {len(optimized)}")
except ImportError:
    print("cirq_google not installed -- pip install cirq-google")

The compiler applies gate decompositions, merges adjacent single-qubit gates into a single PhasedXZGate, and routes two-qubit gates through the grid topology.

Comparing Ideal vs Noisy XEB

# Ideal fidelity: around 1.0 in expectation over many random circuits;
# individual circuits may differ (XEB can exceed 1.0 for a single circuit).
F_ideal = compute_xeb_fidelity(circuit, qubits, n_shots=10000)

# For noisy fidelity, replace ideal_sim.run with noisy_sim.run in the function
# (modify compute_xeb_fidelity to accept a simulator argument)
F_noisy = compute_xeb_fidelity_v2(circuit, qubits, noisy_sim)

print(f"Ideal XEB fidelity:  {F_ideal:.4f}")
print(f"Noisy XEB fidelity:  {F_noisy:.4f}")
print(f"Fidelity degradation: {(F_ideal - F_noisy) / F_ideal * 100:.1f}%")

With the noise parameters above, a 10-depth circuit over 6 qubits will show noticeable fidelity degradation. Increasing depth degrades fidelity further, matching the behavior observed on real Sycamore hardware.

Common Mistakes

Working with Cirq’s superconducting simulation involves several subtle pitfalls. This section covers the most frequent errors and how to avoid them.

Using Non-Adjacent GridQubits Without Routing

Applying a two-qubit gate between non-adjacent GridQubits compiles and simulates without error in Cirq’s software simulator. However, device.validate_circuit() raises a validation error, and submission to real hardware fails. Always validate your circuit against the target device before assuming it runs on hardware:

import cirq

q0 = cirq.GridQubit(0, 0)
q2 = cirq.GridQubit(0, 2)  # Not adjacent to q0

# This simulates fine but is invalid on hardware
bad_circuit = cirq.Circuit([cirq.CNOT(q0, q2)])
sim = cirq.Simulator()
result = sim.simulate(bad_circuit)  # Works in simulation
print("Simulation succeeded, but this circuit is invalid on real hardware.")

Wrong fSim Angle Convention

The theta parameter controls the iSWAP-like swap angle, and phi controls the CZ-like conditional phase. Confusing them produces a completely different gate. Remember: theta is the swap angle (iSWAP at pi/2), phi is the phase angle (CZ at pi).

import numpy as np
import cirq

# Correct: iSWAP has theta=pi/2, phi=0
iswap = cirq.FSimGate(theta=np.pi / 2, phi=0)

# WRONG: this gives CZ, not iSWAP
not_iswap = cirq.FSimGate(theta=0, phi=np.pi / 2)

print("Correct iSWAP unitary:")
print(np.round(cirq.unitary(iswap), 3))
print("\nIncorrect (swapped parameters):")
print(np.round(cirq.unitary(not_iswap), 3))

DensityMatrixSimulator Returns a Different Result Type

cirq.Simulator returns SimulationTrialResult with a final_state_vector attribute. cirq.DensityMatrixSimulator returns DensityMatrixTrialResult with a final_density_matrix attribute. Calling .final_state_vector on a density matrix result raises an AttributeError:

import cirq

q = cirq.GridQubit(0, 0)
circuit = cirq.Circuit([cirq.H(q)])

# Statevector simulator
sv_result = cirq.Simulator().simulate(circuit)
print(type(sv_result))  # SimulationTrialResult
print(sv_result.final_state_vector)  # Works

# Density matrix simulator
dm_result = cirq.DensityMatrixSimulator().simulate(circuit)
print(type(dm_result))  # DensityMatrixTrialResult
print(dm_result.final_density_matrix)  # Works
# dm_result.final_state_vector  # AttributeError!

Mixing GridQubits and LineQubits

GridQubit and LineQubit are different types. If you create some qubits as GridQubit and others as LineQubit, they do not interact as expected and cannot share a device topology:

import cirq

# These are different qubit types and should not be mixed
grid_q = cirq.GridQubit(0, 0)
line_q = cirq.LineQubit(0)

print(f"Same qubit? {grid_q == line_q}")  # False
# Mixing them in one circuit works in simulation but breaks device validation

Always use GridQubit consistently when targeting Google hardware.

Missing cirq-google Installation

The cirq_google module is not included in the base cirq package. Importing it without installing the separate package raises ImportError:

# The base package does not include Google-specific tools
pip install cirq         # Core cirq only
pip install cirq-google  # Required for Sycamore device, target gatesets, etc.

Always wrap import cirq_google in a try/except block or document the dependency.

Wrong Qubit Count in XEB Formula

The XEB fidelity formula uses a factor of 2^n, where n is the number of qubits. Using the wrong n (for example, using the total number of qubits in the device instead of the number of qubits in the circuit) produces completely wrong fidelity values:

import numpy as np

# Suppose you have 6 qubits and mean_p = 0.02
mean_p = 0.02
n_correct = 6
n_wrong = 53  # Accidentally used full Sycamore qubit count

f_correct = 2**n_correct * mean_p - 1  # 2^6 * 0.02 - 1 = 0.28
f_wrong = 2**n_wrong * mean_p - 1      # Astronomically wrong

print(f"Correct F_XEB (n=6):  {f_correct:.4f}")
print(f"Wrong F_XEB (n=53):   {f_wrong:.4e}")  # Nonsensical value

Always use n = len(qubits) where qubits is the list of qubits actually in the circuit.

Key Takeaways

• Transmon qubits use the anharmonicity of a Josephson junction to isolate a two-level computational subspace from higher energy levels.
• cirq.GridQubit models the 2D qubit layout of Google superconducting processors with nearest-neighbor connectivity only.
• cirq.FSimGate(theta, phi) is the hardware-native two-qubit gate; theta=pi/2, phi=pi/6 gives the Sycamore gate. The gate encompasses iSWAP, CZ, and SWAP as special cases.
• Parameterized fSim gates with sympy symbols enable variational algorithms on Google hardware.
• XEB fidelity measures how close a noisy device is to an ideal quantum computer, with F_XEB = 1 for ideal and F_XEB = 0 for fully depolarized output.
• A complete noise model combines amplitude damping (T1), phase damping (T2), depolarizing error, and readout noise.
• Circuit fidelity decays exponentially with the number of two-qubit gates: F = (1 - epsilon)^n_2q.
• Google’s Willow processor (2024) achieves ~100 microsecond T1 and 0.3% two-qubit error, enabling below-threshold error correction.
• cirq.DensityMatrixSimulator with a custom NoiseModel accurately simulates the effect of T1/T2 decay and depolarizing errors.
• cirq_google.optimized_for_sycamore() compiles arbitrary circuits for the hardware-native gate set and qubit topology.
• Always validate circuits against device topology before assuming hardware compatibility.