Simulating Google Sycamore Circuits with Cirq

In October 2019, Google announced that their Sycamore processor had achieved quantum computational advantage, sampling from the output distribution of specific random quantum circuits faster than any classical supercomputer could. The experiment used 53 qubits arranged in a 2D grid, with circuits composed of layers of single-qubit gates and the hardware-native two-qubit FSimGate. This tutorial walks through building and simulating those circuits in Cirq, implementing a noise model that approximates real device behavior, and computing the cross-entropy benchmarking (XEB) fidelity that Google used to certify the result.

The Sycamore Topology

Sycamore uses a 9x6 grid of qubits with some corners removed, giving 54 physical qubits (one defective, leaving 53 usable). Cirq represents grid qubits as GridQubit(row, col).

import cirq
import numpy as np
import matplotlib.pyplot as plt
from itertools import product

# Build the Sycamore qubit grid (simplified 3x3 for local simulation)
def sycamore_grid(rows=3, cols=3):
    """Return a rectangular grid of GridQubits."""
    return [cirq.GridQubit(r, c) for r in range(rows) for c in range(cols)]

# Define the coupling pattern: Sycamore uses nearest-neighbor connectivity
def get_neighbor_pairs(qubits, rows, cols):
    """Return all nearest-neighbor pairs for a rows x cols grid."""
    qubit_set = set(qubits)
    pairs = []
    for r, c in product(range(rows), range(cols)):
        q = cirq.GridQubit(r, c)
        if q not in qubit_set:
            continue
        for dr, dc in [(0, 1), (1, 0)]:
            neighbor = cirq.GridQubit(r + dr, c + dc)
            if neighbor in qubit_set:
                pairs.append((q, neighbor))
    return pairs

rows, cols = 3, 3
qubits = sycamore_grid(rows, cols)
pairs = get_neighbor_pairs(qubits, rows, cols)
print(f"Qubits: {len(qubits)}")
print(f"Coupling pairs: {len(pairs)}")

The FSimGate with Sycamore Parameters

The native two-qubit gate on Sycamore is the FSimGate(theta, phi):

$\text{FSim}(\theta, \phi) = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & \cos\theta & -i\sin\theta & 0 \\ 0 & -i\sin\theta & \cos\theta & 0 \\ 0 & 0 & 0 & e^{-i\phi} \end{pmatrix}$

The Sycamore gate uses $\theta = \pi/2$ (a full iSWAP-like interaction) and $\phi = \pi/6$. These parameters were chosen to maximize gate speed while maintaining high fidelity on the superconducting device.

# Sycamore FSimGate parameters
THETA_SYCAMORE = np.pi / 2
PHI_SYCAMORE = np.pi / 6

sycamore_gate = cirq.FSimGate(theta=THETA_SYCAMORE, phi=PHI_SYCAMORE)
print("Sycamore gate unitary:")
print(np.round(cirq.unitary(sycamore_gate), 3))

Generating Random Quantum Circuits

Google’s supremacy experiment used random circuits with a specific structure: alternating layers of random single-qubit gates (chosen from $\{W, X^{1/2}, Y^{1/2}, T\}$) and layers of two-qubit FSimGates applied to interleaved subsets of pairs.

SINGLE_QUBIT_GATES = [
    cirq.X**0.5,
    cirq.Y**0.5,
    cirq.PhasedXPowGate(phase_exponent=0.25, exponent=0.5),  # W gate
    cirq.T,
]

def random_single_qubit_layer(qubits, rng):
    """Apply a random single-qubit gate to each qubit."""
    return [rng.choice(SINGLE_QUBIT_GATES).on(q) for q in qubits]

def two_qubit_layer(pairs, subset_indices):
    """Apply FSimGate to a subset of pairs."""
    selected = [pairs[i] for i in subset_indices if i < len(pairs)]
    return [sycamore_gate.on(q0, q1) for q0, q1 in selected]

def build_random_circuit(qubits, pairs, n_cycles, rng=None):
    """
    Build a random quantum circuit in the style of Google's supremacy experiment.
    Alternates single-qubit layers with partitioned two-qubit layers.
    """
    if rng is None:
        rng = np.random.default_rng(42)

    # Partition pairs into 4 groups (A, B, C, D pattern)
    partitions = [[], [], [], []]
    for i, (q0, q1) in enumerate(pairs):
        partitions[i % 4].append(i)

    circuit = cirq.Circuit()
    cycle_order = [0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3]  # ABCDABCD...

    for cycle in range(n_cycles):
        # Single-qubit layer
        circuit.append(random_single_qubit_layer(qubits, rng))
        # Two-qubit layer
        partition = partitions[cycle_order[cycle % len(cycle_order)]]
        circuit.append(two_qubit_layer(pairs, partition))

    # Final single-qubit layer
    circuit.append(random_single_qubit_layer(qubits, rng))
    # Measure all
    circuit.append(cirq.measure(*qubits, key='result'))

    return circuit

rng = np.random.default_rng(1234)
circuit = build_random_circuit(qubits, pairs, n_cycles=6, rng=rng)
print(f"Circuit depth: {len(circuit)}")
print(f"Circuit:\n{circuit[:4]}")  # Print first 4 moments

Noiseless Statevector Simulation

# Noiseless simulation using Cirq's simulator
simulator = cirq.Simulator()
result = simulator.simulate(circuit[:-1])  # exclude measurement moment
state_vector = result.final_state_vector

print(f"State vector dimension: {len(state_vector)}")
print(f"Sum of probabilities: {np.sum(np.abs(state_vector)**2):.6f}")

# Sample from the ideal distribution
ideal_probs = np.abs(state_vector)**2
n_samples = 5000
ideal_samples = np.random.choice(len(ideal_probs), size=n_samples, p=ideal_probs)

Realistic Noise Model

Real quantum devices suffer from depolarizing noise (two-qubit gate errors dominate), readout errors, and decoherence. Cirq’s DensityMatrixSimulator handles mixed states.

# Build a realistic noise model
def build_sycamore_noise_model(
    single_qubit_error=0.001,
    two_qubit_error=0.006,
    readout_error=0.01,
):
    """
    Approximate Sycamore noise model.
    Error rates based on Google's reported values circa 2019.
    """
    return cirq.ConstantQubitNoiseModel(
        qubit_noise_gate=cirq.DepolarizingChannel(p=single_qubit_error)
    )

# For two-qubit gate noise, we insert noise manually
def build_noisy_circuit(circuit, two_qubit_error=0.006, readout_error=0.01):
    """Insert depolarizing noise after each two-qubit gate."""
    noisy_ops = []
    for moment in circuit[:-1]:  # skip measurement
        noisy_ops.append(moment)
        for op in moment.operations:
            if len(op.qubits) == 2:
                for q in op.qubits:
                    noisy_ops.append(
                        cirq.depolarize(two_qubit_error / 2).on(q)
                    )

    noisy_circuit = cirq.Circuit(noisy_ops)
    noisy_circuit.append(cirq.measure(*qubits, key='result'))
    return noisy_circuit

noisy_circuit = build_noisy_circuit(circuit)
noisy_sim = cirq.DensityMatrixSimulator()
noisy_result = noisy_sim.simulate(noisy_circuit[:-1])

# Get diagonal (probabilities) from density matrix
rho = noisy_result.final_density_matrix
noisy_probs = np.real(np.diag(rho))
noisy_probs = np.maximum(noisy_probs, 0)
noisy_probs /= noisy_probs.sum()

print(f"Max ideal probability: {ideal_probs.max():.6f}")
print(f"Max noisy probability: {noisy_probs.max():.6f}")
print(f"Uniform probability:   {1/len(ideal_probs):.6f}")

Cross-Entropy Benchmarking (XEB) Fidelity

XEB fidelity quantifies how closely a noisy device’s output distribution matches the ideal circuit output. Given bitstrings $\{x_i\}$ sampled from a device, and ideal probabilities $p(x)$ from classical simulation:

$F_{\text{XEB}} = 2^n \cdot \mathbb{E}_{\text{device}}[p(x)] - 1$

For a perfect device, $\mathbb{E}[p(x)] = \mathbb{E}_{p}[p(x)] = \sum_x p(x)^2 \approx 2/2^n$ (Porter-Thomas), giving $F_{\text{XEB}} = 1$. For a fully depolarized device sampling uniformly, $\mathbb{E}[p(x)] = 1/2^n$, giving $F_{\text{XEB}} = 0$.

def xeb_fidelity(device_samples, ideal_probs):
    """
    Compute XEB fidelity.
    device_samples: array of integer bitstring indices sampled from device
    ideal_probs: ideal probability distribution as numpy array
    """
    n = int(np.round(np.log2(len(ideal_probs))))
    N = 2**n
    avg_ideal_prob = np.mean(ideal_probs[device_samples])
    return N * avg_ideal_prob - 1

# Sample from noisy distribution
n_xeb_samples = 2000
noisy_samples = np.random.choice(len(noisy_probs), size=n_xeb_samples, p=noisy_probs)

F_xeb_noisy = xeb_fidelity(noisy_samples, ideal_probs)
print(f"XEB fidelity (noisy model): {F_xeb_noisy:.4f}")

# XEB with ideal samples (should be ~1)
ideal_samples_xeb = np.random.choice(len(ideal_probs), size=n_xeb_samples, p=ideal_probs)
F_xeb_ideal = xeb_fidelity(ideal_samples_xeb, ideal_probs)
print(f"XEB fidelity (ideal):       {F_xeb_ideal:.4f}")

# XEB with uniform samples (should be ~0)
uniform_samples = np.random.choice(len(ideal_probs), size=n_xeb_samples)
F_xeb_uniform = xeb_fidelity(uniform_samples, ideal_probs)
print(f"XEB fidelity (uniform):     {F_xeb_uniform:.4f}")

XEB vs Circuit Depth

As circuits grow deeper, noise accumulates and XEB fidelity drops exponentially. At each cycle, the fidelity is multiplied by approximately $(1 - \epsilon_{\text{2q}})^{n_{\text{gates}}}$, where $\epsilon_{\text{2q}}$ is the two-qubit gate error rate.

depths = range(2, 6, 2)
fidelities = []

for d in depths:
    c = build_random_circuit(qubits, pairs, n_cycles=d, rng=np.random.default_rng(42))
    noisy_c = build_noisy_circuit(c)
    sim_result = noisy_sim.simulate(noisy_c[:-1])
    rho_d = sim_result.final_density_matrix
    probs_d = np.real(np.diag(rho_d))
    probs_d = np.maximum(probs_d, 0) / probs_d.sum()

    ideal_result = simulator.simulate(c[:-1])
    ideal_p_d = np.abs(ideal_result.final_state_vector)**2

    samples_d = np.random.choice(len(probs_d), size=1000, p=probs_d)
    fidelities.append(xeb_fidelity(samples_d, ideal_p_d))

plt.figure(figsize=(7, 4))
plt.plot(list(depths), fidelities, 'o-')
plt.xlabel('Circuit Depth (cycles)')
plt.ylabel('XEB Fidelity')
plt.title('XEB Fidelity vs Depth (Sycamore-like noise)')
plt.axhline(0, color='gray', linestyle='--')
plt.tight_layout()
plt.savefig('xeb_vs_depth.png', dpi=150)

The Classical Simulation Wall

Google claimed their 53-qubit, 20-cycle circuits took 200 seconds to sample one million bitstrings on Sycamore, while the best classical algorithm at the time would take approximately 10,000 years. This calculation has since been contested, with classical algorithms improving rapidly.

The key difficulty for classical simulation is memory: an exact statevector for 53 qubits requires $2^{53} \times 16$ bytes (complex128) = 144 petabytes. Approximate methods such as tensor network contraction can reduce this but trade fidelity for speed.

For local simulation, 16-20 qubits is tractable with exact methods; 30+ qubits requires approximate techniques or distributed computing. Cirq’s qsim backend (via cirq_google) provides optimized statevector simulation for circuits up to ~30 qubits on modern hardware.

This tutorial demonstrates the core machinery behind the landmark experiment: random circuit construction on a 2D qubit grid, the hardware-native FSimGate, noise modeling with density matrices, and XEB fidelity as the certification metric.