Quantum scrambling experiment example

View on QuantumAI Run in Google Colab View source on GitHub Download notebook

This tutorial shows an experiment example from the publication Information scrambling in computationally complex quantum circuits (henceforth "OTOC paper"), in particular how to recreate Figure 1.

Background on quantum scrambling

Quantum scrambling describes how interaction in a quantum system disperses local information into its many degrees of freedom. Analogous to classical chaos, scrambling manifests itself as a butterfly effect wherein a local perturbation is rapidly amplified over time. The perturbation is modeled by a local operator \(\hat{O}\) (the butterfly operator) acting on a qubit (the butterfly qubit). Under a dynamical process \(\hat{U} = \hat{U}(t)\), the butterfly operator evolves in the Heisenberg picture as

\[ \hat{O} (t) = \hat{U}^\dagger \hat{O} \hat{U} \]

Expanding the butterfly operator in a basis

\[ \hat{O} (t) = \sum_{i = 1}^{n_p} w_i B_i \]

where \(B_i = B^{(1)}_i \otimes \cdots \otimes B^{(n)}_i\) are basis elements and \(w_i\) are coefficients, we can describe the two different mechanisms of quantum scrambling as follows:

  • Operator spreading occurs when basis elements \(B_i\) require a higher weight (fewer non-identity terms) on average to describe \(\hat{O}(t)\).
  • Operator entanglement occurs when a larger \(n_p\) (summation index) is required to describe \(\hat{O}(t)\).

We experimentally evaluate the out-of-time order correlator (OTOC) \(C(t)\) between \(\hat{O}(t)\) and a measurement operator \(\hat{M}\) which is a Pauli operator acting on another qubit (the measurement qubit)

\[ C(t) := \langle \hat{O} ^\dagger (t) M^\dagger \hat{O} (t) M \rangle \]

over a collection of quantum circuits with microscopic differences (e.g., in the phases of gates). Operator scrambling is then reflected in the average OTOC value \(\bar{C}(t)\) over different circuits:

  • \(\bar{C}(t) = 1\) when \(\hat{O}(t)\) and \(\hat{M}\) do not overlap.
  • \(\bar{C}(t) < 1\) when \(\hat{O}(t)\) and \(\hat{M}\) overlap.
  • \(\bar{C}(t) \rightarrow 0\) in the fully scrambled limit where the commutation between \(\hat{O}(t)\) and \(\hat{M}\) is completely randomized.

Note that if operator entanglement is also present (\(n_p \gg 1\)), \(\bar{C}(t)\) approaches \(0\) for all circuits and their fluctuation vanishes as well.

In what follows, we show how to build OTOC circuits, run them and compute \(\bar{C}\), then visualize results and detect operator spreading.


    import recirq
except ImportError:
    print("Installing ReCirq...")
    !pip install git+https://github.com/quantumlib/recirq --quiet
    print("Installed ReCirq!")
import matplotlib.pyplot as plt
import numpy as np

import cirq
import cirq_google as cg
from recirq import otoc

No project_id provided and environment variable GOOGLE_CLOUD_PROJECT not set.
Using a noisy simulator.

Build the OTOC circuits

Picking qubits

First we specify a line of qubits to use. If running on the cirq_google.Engine, use the latest calibration report to identify a good line of qubits.

if not use_noisy_simulator:
    calibration = cg.get_engine_calibration(processor_id=processor_id)

Note that the ancilla qubit must be adjacent to the first qubit, and the first qubit is treated as the measurement qubit described above.

qubit_indices = [(4, 1), (4, 2), (5, 2), (5, 3), (6, 3), (7, 3), (7, 4), (7, 5)]
ancilla_index = (5, 1)

qubits = [cirq.GridQubit(*idx) for idx in qubit_indices]
ancilla = cirq.GridQubit(*ancilla_index)
num_qubits = len(qubits)

Qubit interactions and experiment parameters

Now we specify how qubits interact (in other words, the pattern of two-qubit gates). For the 1D chain in this example, only two layers of two-qubit gates (CZ is used here) are needed to enact all interactions.

# Staggered two-qubit interaction pattern.
int_sets = [
    {(qubits[i], qubits[i + 1]) for i in range(0, num_qubits - 1, 2)},
    {(qubits[i], qubits[i + 1]) for i in range(1, num_qubits - 1, 2)},

# Use CZ for the two-qubit interaction.
forward_ops = {
    (qubit_indices[i], qubit_indices[i + 1]): cirq.Circuit([cirq.CZ(qubits[i], qubits[i + 1])])
    for i in range(num_qubits - 1)

# CZ is self-inverse.
reverse_ops = forward_ops.copy()

Finally, we specify a number of trials to average over when computing \(\bar{C}(t)\), and a list of cycles which corresponds to discrete time values \(t\) at which to compute \(\bar{C}\).

num_trials = 4  # For averaging.
cycles = range(0, 20 + 1, 2)  # Number of circuit cycles.

Build the circuits

With the specified parameters, we use recirq.otoc.build_otoc_circuits to create the set of circuits. First we show an example for a given number of cycles and given butterfly qubit.

# Example number of cycles and butterfly qubit.
cycle = cycles[1]
q_b = qubits[2]

otoc_circuits = otoc.build_otoc_circuits(
    sq_gates=np.random.choice(8, (num_qubits, max(cycles))),

The output is an instance of recirq.otoc.otoc_circuits.OTOCCircuits which contains lists of circuits with \(X\), \(Y\), and \(Z\) butterfly operators on the butterfly qubit. Additionally, a list of circuits with no butterfly operator (equivalently, with the identity \(I\) as the butterfly operator) is included to renormalize results.

print("Ancilla qubit:", ancilla)
print("Measurement qubit:", qubits[0])
print("Butterfly qubit:", q_b)

print("\nExample OTOC circuit:\n")
Ancilla qubit: (5, 1)
Measurement qubit: (4, 1)
Butterfly qubit: (5, 2)

Example OTOC circuit:

Note that Pauli \(Y\) acts on the chosen butterfly qubit and the measurement is on the ancilla. Compare this circuit to Fig. 1(a) & (b) of the OTOC paper.

Now we create all OTOC circuits for the experiment, varying the number of cycles and the butterfly qubit.

np.random.seed(0)  # Random seed for OTOC circuits.

circuit_list = []
for i in range(num_trials):
    circuits_i = []

    for cycle in cycles:
        circuits_ic = []

        for k, q_b in enumerate(qubits[1:]):
            circs = otoc.build_otoc_circuits(
                sq_gates=np.random.choice(8, (num_qubits, max(cycles))),

            # If q_b is the measurement qubit, add both the normalization 
            # circuits and the circuits with X as the butterfly operator. 
            # Otherwise only add circuits with Y being the butterfly operator.
            if k == 0:


The circuit_list is a three-dimensional array where the first axis corresponds to the trial, the second axis corresponds to the number of cycles, and the third axis corresponds the butterfly operator / qubit (\(Q_b\)).

Run the circuits

Now we run the circuits and process the results.

# Compile to sqrt_iswap.
from tqdm import tqdm

progress = tqdm(total=np.array(circuit_list, dtype=object).size)
for i, circuits_trial in enumerate(circuit_list):
    for j, circuits_cycle in enumerate(circuits_trial):
        for k, circuit_otoc in enumerate(circuits_cycle):
            circuit_list[i][j][k] = cg.optimized_for_sycamore(circuit_otoc, optimizer_type="sqrt_iswap")
100%|██████████| 1408/1408 [14:50<00:00,  1.46s/it]
results = []
for i, circuits_i in enumerate(circuit_list):
    results_i = np.zeros((num_qubits, len(cycles)))

    for c, circuits_ic in enumerate(circuits_i):
        print("\r", f"On trial {i + 1} / {num_trials}, cycle {c + 1} / {len(cycles)}.", end="")

        job = sampler.run_batch(
            params_list=[{} for _ in range(len(circuits_ic))], 
            repetitions=int(2_000 + 10_000 * (c / max(cycles)) ** 3),

        # Compute ⟨σ_y⟩ (or normalization ⟨I⟩).
        for d in range(num_qubits):
            p =  np.mean(job[4 * d    ][0].measurements["z"])
            p -= np.mean(job[4 * d + 1][0].measurements["z"])
            p -= np.mean(job[4 * d + 2][0].measurements["z"])
            p += np.mean(job[4 * d + 3][0].measurements["z"])
            results_i[d, c] = 0.5 * p
On trial 4 / 4, cycle 11 / 11.

Plot the results

First we plot the normalization factor and \(\langle \sigma_y \rangle\) on each butterfly qubit.

# Average the data for the `num_trials` random circuits.
results_ave = np.mean(results, axis=0)

# Plot the result with no butterfly operator.
plt.plot(cycles, results_ave[0, :], "o-", color="black", mec="black", mfc="white", label="Normalization")

# Plot results for each butterfly qubit.
for i in range(1, num_qubits):
    plt.plot(cycles, results_ave[i, :], "o-", mec="black", label="$Q_b =$ {}".format(qubits[i]))

plt.title(f"$Q_a = {ancilla}$")
plt.ylabel(r"$\langle \sigma_y \rangle$")
plt.legend(bbox_to_anchor=(1, 1));


Compare this plot to Fig. 1(c) (left panel) of the OTOC paper. The values of \(\langle \sigma_y \rangle\) decay below 1 before \(\hat{O}\) and \(\hat{M}\) overlap due to noise. The curve labeled "normalization" contains no butterfly operator and its decay is solely due to noise. By renormalizing, we can analyze the effects due to quantum scrambling.

for i in range(1, num_qubits):
    plt.plot(cycles, results_ave[i, :] / results_ave[0, :], "o-", mec="black", label="$Q_b =$ {}".format(qubits[i]))

plt.title(f"$Q_a = {ancilla}$")
plt.ylabel(r"Average OTOC, $\bar{C}$")
plt.legend(bbox_to_anchor=(1, 1));


Note that for each location of \(Q_b\), \(\bar{C}\) retains values near 1 before sufficient circuit cycles have occured to allow an overlap between \(\hat{O}\) and \(\hat{M}\). Beyond these cycles, \(\bar{C}\) converges to \(0\), indicating that \(\hat{O}\) and \(\hat{M}\) have overlapped and no longer commute. These observations show behavior consistent with operator spreading as described at the start of this tutorial. Compare this plot to Fig. 1(c) (right panel) of the OTOC paper.