Nishimori phase transition

Usage estimate: 3 minutes on a Heron r2 processor (NOTE: This is an estimate only. Your runtime may vary.)

Learning outcomes

After going through this tutorial, users should expect the following outcomes:

Understand the Nishimori phase transition and how it manifests as the appearance of long-range entanglement in the random-bond Ising model.
Implement the generation of entanglement by measurement (GEM) protocol on quantum hardware using mid-circuit measurements and constant-depth circuits.
Characterize the transition by extracting the two-point correlation and the normalized variance of the magnetization from experimental data.

Prerequisites

We recommend familiarity with the following topics before going through this tutorial:

Measure qubits, in particular the section on mid-circuit measurement that the GEM protocol relies on.
Exact and noisy simulation with Qiskit Aer primitives, which is how the small-scale section is executed.
Long-range entanglement with dynamic circuits, a companion tutorial that uses the same measurement-based-entanglement paradigm.
Heavy hex lattice, the IBM hardware topology the plaquette lattice is built on.

Background

This tutorial demonstrates how to realize a Nishimori phase transition on a quantum processor. This experiment was originally described in Realizing the Nishimori transition across the error threshold for constant-depth quantum circuits.

The Nishimori phase transition refers to the transition between short- and long-range ordered phases in the random-bond Ising model. On a quantum computer, the long-range ordered phase manifests as a state in which qubits are entangled across the entire device. This highly entangled state is prepared using the generation of entanglement by measurement (GEM) protocol. By utilizing mid-circuit measurements, the GEM protocol is able to entangle qubits across the entire device using circuits of only constant depth. This tutorial uses the implementation of the GEM protocol from the GEM Suite software package.

Requirements

Before starting this tutorial, be sure you have the following installed:

Qiskit SDK v1.0 or later, with visualization support
Qiskit Runtime v0.22 or later (pip install qiskit-ibm-runtime)
Qiskit Aer v0.14 or later (pip install qiskit-aer)
GEM Suite (pip install gem-suite)

Setup

import matplotlib.pyplot as plt
import warnings

from collections import defaultdict

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_aer import AerSimulator

from qiskit.transpiler import generate_preset_pass_manager

from gem_suite import PlaquetteLattice
from gem_suite.experiments import GemExperiment

Small-scale simulator example

In this section, the full workflow is walked through on the noiseless AerSimulator. The plaquette lattice is restricted to a single plaquette (12 qubits) so the simulation stays small and fast, while still exercising every part of the GEM protocol: mid-circuit measurement, the $R_{ZZ}$ angle sweep, decoding, and the normalized-variance analysis. The same workflow is later scaled up to multiple plaquettes and the full lattice on real hardware.

Step 1: Map classical inputs to a quantum problem

The GEM protocol works on a quantum processor with qubit connectivity described by a lattice. Today's IBM quantum processors use the heavy hex lattice. The qubits of the processor are grouped into plaquettes based on which unit cell of the lattice they occupy. Because a qubit might occur in more than one unit cell, the plaquettes are not disjoint. On the heavy hex lattice, a plaquette contains 12 qubits. The plaquettes themselves also form lattice, where two plaquettes are connected if they share any qubits. On the heavy hex lattice, neighboring plaquettes share 3 qubits.

In the GEM Suite software package, the fundamental class for implementing the GEM protocol is PlaquetteLattice, which represents the lattice of plaquettes (which is distinct from the heavy hex lattice). A PlaquetteLattice can be initialized from a qubit coupling map. Currently, only heavy hex coupling maps are supported.

The following code cell initializes a plaquette lattice from the coupling map of an IBM quantum processor. The plaquette lattice does not always encompass the entire hardware. For example, ibm_torino has 133 total qubits but the largest plaquette lattice that fits on the device uses only 125 of them, comprising 18 plaquettes; ibm_pittsburgh (156 qubits) similarly fits 144 qubits into 21 plaquettes. The same pattern holds for other heavy-hex IBM Quantum® processors with different qubit counts.

# QiskitRuntimeService.save_account(channel="ibm_quantum", token="<YOUR_API_KEY>", overwrite=True, set_as_default=True)
service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, simulator=False, min_num_qubits=127
)
aer_backend = AerSimulator.from_backend(backend)
plaquette_lattice = PlaquetteLattice.from_coupling_map(backend.coupling_map)

print(f"Number of qubits in backend: {backend.num_qubits}")
print(
    f"Number of qubits in plaquette lattice: {len(list(plaquette_lattice.qubits()))}"
)
print(f"Number of plaquettes: {len(list(plaquette_lattice.plaquettes()))}")

Output:

Number of qubits in backend: 156
Number of qubits in plaquette lattice: 144
Number of plaquettes: 21

You can visualize the plaquette lattice by generating a diagram of its graph representation. In the diagram, the plaquettes are represented by labeled hexagons, and two plaquettes are connected by an edge if they share qubits.

plaquette_lattice.draw_plaquettes()

Output:

You can retrieve information about individual plaquettes, such as the qubits they contain, using the plaquettes method.

# Get a list of the plaquettes
plaquettes = list(plaquette_lattice.plaquettes())
# Display information about plaquette 0
plaquettes[0]

Output:

PyPlaquette(index=0, qubits=[3, 4, 5, 6, 7, 16, 17, 23, 24, 25, 26, 27], neighbors=[4, 3, 1])

You can also produce a diagram of the underlying qubits that form the plaquette lattice.

plaquette_lattice.draw_qubits()

Output:

In addition to the qubit labels and the edges indicating which qubits are connected, the diagram contains three additional pieces of information that are relevant to the GEM protocol:

Each qubit is either shaded (gray) or unshaded. The shaded qubits are "site" qubits that represent the sites of the Ising model, and the unshaded qubits are "bond" qubits used to mediate interactions between the site qubits.
Each site qubit is labeled either (A) or (B), indicating one of two roles a site qubit can play in the GEM protocol (the roles are explained later).
Each edge is colored using one of six colors, thus partitioning the edges into six groups. This partitioning determines how two-qubit gates can be parallelized, as well as different scheduling patterns that are likely to incur different amounts of error on a noisy quantum processor. Because edges in a group are disjoint, a layer of two-qubit gates can be applied on those edges simultaneously. In fact, it is possible to partition the six colors into three groups of two colors such that the union of each group of two colors is still disjoint. Therefore, only three layers of two-qubit gates are needed to activate every edge. There are 12 ways to so partition the six colors, and each such partition yields a different 3-layer gate schedule.

Now that you have created a plaquette lattice, the next step is to initialize a GemExperiment object, passing both the plaquette lattice and the backend that you intend to run the experiment on. The GemExperiment class manages the actual implementation of the GEM protocol, including generating circuits, submitting jobs, and analyzing the data. The following code cell initializes the experiment class while restricting the plaquette lattice to a single plaquette (12 qubits), keeping the simulation small and fast. The full plaquette lattice is used later when scaling up to real hardware.

# Filter the plaquette lattice down to a single plaquette (12 qubits)
# so the AerSimulator run stays fast. The full lattice is used later
# in the large-scale hardware example.
gem_exp = GemExperiment(plaquette_lattice.filter([9]), backend=aer_backend)

# visualize the plaquette lattice after filtering
plaquette_lattice.filter([9]).draw_qubits()

Output:

A GEM protocol circuit is built using the following steps:

Prepare the all- $|+\rangle$ state by applying a Hadamard gate to every qubit.
Apply an $R_{ZZ}$ gate between every pair of connected qubits. This can be achieved using 3 layers of gates. Each $R_{ZZ}$ gate acts on a site qubit and a bond qubit. If the site qubit is labeled (B), then the angle is fixed to $\frac{\pi}{2}$ . If the site qubit is labeled (A), then the angle is allowed to vary, producing different circuits. By default, the range of angles is set to 21 equally spaced points between $0$ and $\frac{\pi}{2}$ , inclusive.
Measure each bond qubit in the Pauli $X$ basis. Since qubits are measured in the Pauli $Z$ basis, this can be accomplished by applying a Hadamard gate before measuring the qubit.

Note that the paper cited in the introduction to this tutorial uses a different convention for the $R_{ZZ}$ angle, which differs from the convention used in this tutorial by a factor of 2.

In step 3, only the bond qubits are measured. To understand what state the site qubits remain in, it is instructive to consider the case that the $R_{ZZ}$ angle applied to site qubits (A) in step 2 is equal to $\frac{\pi}{2}$ . In this case, the site qubits are left in a highly entangled state similar to the GHZ state,

\lvert \text{GHZ} \rangle = \lvert 00 \cdots 00 \rangle + \lvert 11 \cdots 11 \rangle.

Due to the randomness in the measurement outcomes, the actual state of the site qubits might be a different state with long-range order, for example, $\lvert 00110 \rangle + \lvert 11001 \rangle$ . However, the GHZ state can be recovered by applying a decoding operation based on the measurement outcomes. When the $R_{ZZ}$ angle is tuned down from $\frac{\pi}{2}$ , the long-range order can still be recovered up until a critical angle, which in the absence of noise, is approximately $0.3 \pi$ . Below this angle, the resulting state no longer exhibits long-range entanglement. This transition between the presence and absence of long-range order is the Nishimori phase transition.

In the description above, the site qubits were left unmeasured, and the decoding operation can be performed by applying quantum gates. In the experiment as implemented in the GEM suite, which this tutorial follows, the site qubits are in fact measured, and the decoding operation is applied in a classical post-processing step.

In the description above, the decoding operation can be performed by applying quantum gates to the site qubits to recover the quantum state. However, if the goal is to immediately measure the state, for example, for characterization purposes, then the site qubits are measured together with the bond qubits, and the decoding operation can be applied in a classical post-processing step. This is how the experiment is implemented in the GEM suite, which this tutorial follows.

In addition to depending on the $R_{ZZ}$ angle in step 2, which by default sweeps across 21 values, the GEM protocol circuit also depends on the scheduling pattern used to implement the 3 layers of $R_{ZZ}$ gates. As discussed previously, there are 12 such scheduling patterns. Therefore, the total number of circuits in the experiment is $21 \times 12 = 252$ .

The circuits of the experiment can be generated using the circuits method of the GemExperiment class.

circuits = gem_exp.circuits()
print(f"Total number of circuits: {len(circuits)}")

Output:

Total number of circuits: 252

For the purposes of this tutorial, it is enough to consider just a single scheduling pattern. The following code cell restricts the experiment to the first scheduling pattern. As a result, the experiment only has 21 circuits, one for each $R_{ZZ}$ angle swept over.

# Restrict experiment to the first scheduling pattern
gem_exp.set_experiment_options(schedule_idx=0)

# There are less circuits now
circuits = gem_exp.circuits()
print(f"Total number of circuits: {len(circuits)}")

# Print the RZZ angles swept over
print(f"RZZ angles:\n{gem_exp.parameters()}")

Output:

Total number of circuits: 21
RZZ angles:
[0.         0.07853982 0.15707963 0.23561945 0.31415927 0.39269908
 0.4712389  0.54977871 0.62831853 0.70685835 0.78539816 0.86393798
 0.9424778  1.02101761 1.09955743 1.17809725 1.25663706 1.33517688
 1.41371669 1.49225651 1.57079633]

The following code cell draws a diagram of the circuit at index 5. To reduce the size of the diagram, the measurement gates at the end of the circuit are removed.

# Get the circuit at index 5
circuit = circuits[5]
# Remove the final measurements to ease visualization
circuit.remove_final_measurements()
# Draw the circuit
circuit.draw("mpl", fold=-1, scale=0.5)

Output:

Step 2: Optimize problem for quantum hardware execution

Transpiling quantum circuits for execution on hardware typically involves a number of stages. Typically, the stages that incur the most computational overhead are choosing the qubit layout, routing the two-qubit gates to conform to the qubit connectivity of the hardware, and optimizing the circuit to minimize its gate count and depth. In the GEM protocol, the layout and routing stages are unnecessary because the hardware connectivity is already incorporated into the design of the protocol. The circuits already have a qubit layout, and the two-qubit gates are already mapped onto native connections. Furthermore, in order to preserve the structure of the circuit as the $R_{ZZ}$ angle is varied, only very basic circuit optimization should be performed.

The GemExperiment class transparently transpiles circuits when executing the experiment. The layout and routing stages are already overridden by default to do nothing, and circuit optimization is performed at a level that only optimizes single-qubit gates. However, you can override or pass additional options using the set_transpile_options method. For the sake of visualization, the following code cell manually transpiles the circuit displayed previously, and draws the transpiled circuit.

# Demonstrate setting transpile options
gem_exp.set_transpile_options(
    optimization_level=1  # This is the default optimization level
)
pass_manager = generate_preset_pass_manager(
    backend=aer_backend,
    initial_layout=list(gem_exp.physical_qubits),
    **dict(gem_exp.transpile_options),
)
transpiled = pass_manager.run(circuit)
transpiled.draw("mpl", idle_wires=False, fold=-1, scale=0.5)

Output:

Step 3: Execute using Qiskit primitives

To execute the GEM protocol circuits on the hardware, call the run method of the GemExperiment object. You can specify the number of shots you want to sample from each circuit. The run method returns an ExperimentData object which you should save to a variable. Note that the run method only submits jobs without waiting for them to finish, so it is a non-blocking call.

exp_data = gem_exp.run(shots=10_000)

To wait for the results, call the block_for_results method of the ExperimentData object. This call will cause the interpreter to hang until the jobs are finished.

# The noiseless AerSimulator produces zero-variance UFloat objects in the
# analysis, which triggers a harmless warning from the `uncertainties`
# library. Suppress it so the output stays clean.
with warnings.catch_warnings():
    warnings.filterwarnings(
        "ignore", message="Using UFloat objects with std_dev==0"
    )
    exp_data.block_for_results()
exp_data

Output:

ExperimentData(GemExperiment, 90bf2a90-f729-4c4e-a6da-664aecb11039, job_ids=['04a7c405-47fd-46ca-aa4b-aaf7e339cfbe'], metadata=<5 items>, figure_names=['two_point_correlation.svg', 'normalized_variance.svg', 'plaquette_ops.svg', 'bond_ops.svg'])

Step 4: Post-process and return result in desired classical format

At an $R_{ZZ}$ angle of $\frac{\pi}{2}$ , the decoded state would be the GHZ state in the absence of noise. The long-range order of the GHZ state can be visualized by plotting the magnetization of the measured bitstrings. The magnetization $M$ is defined as the sum of the single-qubit Pauli $Z$ operators,

M = \sum_{j=1}^N Z_j,

where $N$ is the number of site qubits. Its value for a bitstring is equal to the difference between the number of zeros and the number of ones. Measuring the GHZ state yields the all zeros state or the all ones state with equal probability, so the magnetization would be $+N$ half of the time and $-N$ the other half of the time. In the presence of errors due to noise, other values would also appear, but if the noise is not too great, the distribution would still be peaked near $+N$ and $-N$ .

For the raw bitstrings before decoding, the distribution of the magnetization would be equivalent to that of uniformly random bitstrings, in the absence of noise.

The following code cell plots the magnetization of the raw bitstrings and the decoded bitstrings at the $R_{ZZ}$ angle of $\frac{\pi}{2}$ .

def magnetization_distribution(
    counts_dict: dict[str, int],
) -> dict[str, float]:
    """Compute magnetization distribution from counts dictionary."""
    # Construct dictionary from magnetization to count
    mag_dist = defaultdict(float)
    for bitstring, count in counts_dict.items():
        mag = bitstring.count("0") - bitstring.count("1")
        mag_dist[mag] += count
    # Normalize
    shots = sum(counts_dict.values())
    for mag in mag_dist:
        mag_dist[mag] /= shots
    return mag_dist


# Get counts dictionaries with and without decoding
data = exp_data.data()
# Get the last data point, which is at the angle for the GHZ state
raw_counts = data[-1]["counts"]
# Without decoding
site_indices = [
    i for i, q in enumerate(gem_exp.plaquettes.qubits()) if q.role == "Site"
]
site_raw_counts = defaultdict(int)
for key, val in raw_counts.items():
    site_str = "".join(key[-1 - i] for i in site_indices)
    site_raw_counts[site_str] += val
# With decoding
_, site_decoded_counts = gem_exp.plaquettes.decode_outcomes(
    raw_counts, return_counts=True
)

# Compute magnetization distribution
raw_magnetization = magnetization_distribution(site_raw_counts)
decoded_magnetization = magnetization_distribution(site_decoded_counts)

# Plot
plt.bar(*zip(*raw_magnetization.items()), label="raw")
plt.bar(*zip(*decoded_magnetization.items()), label="decoded", width=0.3)
plt.legend()
plt.xlabel("Magnetization")
plt.ylabel("Frequency")
plt.title("Magnetization distribution with and without decoding")

Output:

Text(0.5, 1.0, 'Magnetization distribution with and without decoding')

To more rigorously characterize the long-range order, you can examine the average two-point correlation $f$ , defined as

f = \frac{1}{N^2} \left(\langle M^2 \rangle - \langle M \rangle ^2\right).

A higher value indicates a greater degree of entanglement. The GemExperiment class automatically computes this value for the decoded bitstrings as part of processing the experimental data. It stores a figure that is accessible via the figure method of the experiment data class. In this case, the name of the figure is two_point_correlation.

exp_data.figure("two_point_correlation")

Output:

To determine the critical point of the Nishimori phase transition, you can look at the normalized variance of $M^2 / N$ , defined as

g = \frac{1}{N^3} \left(\langle M^4 \rangle - \langle M^2 \rangle^2\right),

which quantifies the amount of fluctuation in the squared magnetization. This value is maximized at the critical point of the Nishimori phase transition. In the absence of noise, the critical point occurs at approximately $0.3 \pi$ . In the presence of noise, the critical point is shifted higher, but the phase transition is still observed as long as the critical point occurs below $0.5 \pi$ .

exp_data.figure("normalized_variance")

Output:

Large-scale hardware example

Having validated the protocol on a simulator, the experiment is now scaled up and run on the real quantum hardware backend selected in the Setup section. Two larger problem sizes are used:

Six plaquettes (~49 qubits): a mid-size run that already shows the rightward shift of the critical point under hardware noise.
The full plaquette lattice: every plaquette the device's heavy-hex topology supports (for example, 18 plaquettes / 125 qubits on ibm_torino or 21 plaquettes / 144 qubits on ibm_pittsburgh), entangling qubits across the entire device with constant-depth circuits.

The single code cell below is self-contained: it builds the plaquette lattice from the backend's coupling map and runs both experiments, so this section can be executed after the Setup cells without first running the small-scale section.

# -------------------------Step 1-------------------------
# Initialize the runtime service, pick a real quantum hardware backend,
# and build the plaquette lattice from its coupling map. This is repeated
# from the small-scale example so this cell can run standalone after the
# Setup section. The full plaquette lattice is the "large-scale" target;
# a six-plaquette subset (range(3, 9)) is also used to show an intermediate
# scaling step.
service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, simulator=False, min_num_qubits=127
)
plaquette_lattice = PlaquetteLattice.from_coupling_map(backend.coupling_map)

# Build a GemExperiment for the full plaquette lattice and one for the
# six-plaquette subset, each restricted to a single scheduling pattern so
# the experiment has one circuit per RZZ angle (21 circuits total).
gem_exp_full = GemExperiment(plaquette_lattice, backend=backend)
gem_exp_full.set_experiment_options(schedule_idx=0)
gem_exp_6 = GemExperiment(
    plaquette_lattice.filter(range(3, 9)), backend=backend
)
gem_exp_6.set_experiment_options(schedule_idx=0)

circuits = gem_exp_full.circuits()
print(f"Total number of circuits (full lattice): {len(circuits)}")

# -------------------------Step 2-------------------------
# GemExperiment transpiles internally for the target backend: the layout
# and routing stages are overridden because the plaquette lattice already
# matches the hardware connectivity, and optimization is restricted so the
# RZZ angle structure is preserved. The code below manually transpiles one
# circuit from the six-plaquette experiment with the same settings this
# experiment will use, and draws it for inspection. (The full-lattice
# transpiled circuit has too many qubits to visualize cleanly, so the
# six-plaquette circuit is used here as a representative example.)
gem_exp_6.set_transpile_options(optimization_level=1)
circuits_6 = gem_exp_6.circuits()
pass_manager = generate_preset_pass_manager(
    backend=backend,
    initial_layout=list(gem_exp_6.physical_qubits),
    **dict(gem_exp_6.transpile_options),
)
transpiled = pass_manager.run(circuits_6[5])
display(transpiled.draw("mpl", idle_wires=False, fold=-1, scale=0.5))

# -------------------------Step 3-------------------------
# Run both problem sizes on real hardware:
#   1. Six plaquettes (~49 qubits) — an intermediate scale-up.
#   2. The full plaquette lattice — every plaquette the device supports.
exp_data_6 = gem_exp_6.run(shots=10_000, job_tags=["TUT_NPT"])
exp_data_full = gem_exp_full.run(shots=10_000, job_tags=["TUT_NPT"])
exp_data_6.block_for_results()
exp_data_full.block_for_results()

# -------------------------Step 4-------------------------
# Plot the normalized variance at each scale. The peak marks the critical
# point of the Nishimori transition; as the system grows, hardware noise
# shifts the peak rightward.
display(exp_data_6.figure("normalized_variance"))
exp_data_full.figure("normalized_variance")

Output:

Total number of circuits (full lattice): 21

Note that, depending on the noise level of the backend used, the normalized-variance curves at the larger sizes may not show a clear peak within the swept angle range. In the runs above, the peak has been pushed all the way to $0.5 \pi$ , the right edge of the sweep (the analysis reports critical_angle = 0.5000 for both the six-plaquette and full-lattice runs). This means hardware noise has shifted the critical point to (or just past) the boundary of the protocol's physically meaningful angle range, so the transition is at the edge of what this sweep can resolve.

Conclusion

In this tutorial, you realized a Nishimori phase transition on a quantum processor using the GEM protocol. The metrics that you examined during post-processing, in particular the two-point correlation and the normalized variance, serve as benchmarks of the device's ability to generate long-range entangled states. These benchmarks extend the utility of the GEM protocol beyond probing interesting physics. As part of the protocol, you entangled qubits across the entire device using circuits of only constant depth. This feat is only possible due to the protocol's use of mid-circuit measurements. In this experiment, the entangled state was immediately measured, but an interesting avenue to explore would be to continue using this state in additional quantum processing!

Next steps

Recommendations

If you found this work interesting, you might be interested in the following material:

References

[1] E. H. Chen, G.-Y. Zhu, R. Verresen, A. Seif, E. Bäumer, D. Layden, N. Tantivasadakarn, G. Zhu, S. Sheldon, A. Vishwanath, S. Trebst, A. Kandala. Realizing the Nishimori transition across the error threshold for constant-depth quantum circuits. arXiv:2309.02863 (2023).

[2] GEM Suite software package.