Skip to main content
IBM Quantum Platform

Circuit cutting for periodic boundary conditions

Usage estimate: 2 minutes on ibm_boston (NOTE: This is an estimate only. Your runtime may vary.)

Learning outcomes

  • Use a virtual connection to transpile a periodic chain without routing its boundary gates.
  • Replace those gates with circuit-cutting experiments and reconstruct expectation values from Sampler results.
  • Compare cut and uncut circuits, accounting for circuit depth and sampling cost.

Prerequisites

You should be familiar with Qiskit patterns, transpilation, and the Sampler primitive.

Background

Periodic boundary conditions connect the last qubit in a chain to the first. If these qubits are far apart on a quantum processor, implementing their interaction requires routing, which can add two-qubit gates and increase circuit depth.

Gate cutting replaces a two-qubit gate with a weighted combination of single-qubit operations and measurements. Executing these experiments and combining their results reconstructs expectation values for the original circuit. This can reduce routing overhead, at the cost of additional sampling. See the circuit cutting explanation for the underlying quasiprobability decomposition.

You will first validate this workflow on a six-qubit noiseless simulator. Then you will apply it to a chain of at least 30 qubits on hardware and compare the cut and uncut results. The circuit uses parameters with analytically known expectation values, so the hardware comparison does not require a large statevector simulation.

Requirements

Before starting this tutorial, install:

  • Qiskit SDK v2.3 or later, with visualization support
  • Qiskit Runtime v0.46 or later (pip install 'qiskit-ibm-runtime~=0.46.1')
  • Circuit cutting Qiskit addon v0.10.0 (pip install 'qiskit-addon-cutting~=0.10.0')
  • Qiskit Aer v0.17.2 or later (pip install 'qiskit-aer~=0.17.2')

The hardware example also requires an IBM Quantum® Platform account with access to a QPU with at least 120 qubits and saved credentials. The simulator example runs without an account.

Setup

Import the tools used in both examples. Aer supports the intermediate measurements introduced by gate cutting.

import matplotlib.pyplot as plt
import numpy as np
import rustworkx as rx

from qiskit.circuit.equivalence_library import SessionEquivalenceLibrary
from qiskit.circuit.library import n_local
from qiskit.providers.fake_provider import GenericBackendV2
from qiskit.quantum_info import SparsePauliOp, Statevector
from qiskit.result import sampled_expectation_value
from qiskit.transpiler import CouplingMap, PassManager
from qiskit.transpiler.passes import (
    BasisTranslator,
    GateDirection,
    Optimize1qGatesDecomposition,
)
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit_aer.primitives import SamplerV2 as AerSampler
from qiskit_addon_cutting import (
    cut_gates,
    generate_cutting_experiments,
    reconstruct_expectation_values,
)
from qiskit_ibm_runtime import (
    Batch,
    QiskitRuntimeService,
    SamplerOptions,
    SamplerV2,
)

Small-scale simulator example

Use a six-qubit ring with two layers of entangling gates. The simulated device has linear connectivity, so the ring's closing edge requires routing unless it is cut.

Step 1: Map classical inputs to a quantum problem

The following function creates RX rotation layers interleaved with CX gates on alternating edges of the ring. Two repetitions give two boundary gates to cut. Set the first two rotation layers to zero and the final layer to π\pi. The ideal output is then ∣1⟩⊗n|1\rangle^{\otimes n}, up to a global phase.

Measure the average single-qubit Z expectation and the average adjacent ZZ expectation, including the periodic pair:

MZ=1n∑i=0n−1⟨Zi⟩,CZZ=1n∑i=0n−1⟨ZiZ(i+1) mod n⟩.M_Z = \frac{1}{n}\sum_{i=0}^{n-1}\langle Z_i\rangle, \qquad C_{ZZ} = \frac{1}{n}\sum_{i=0}^{n-1}\langle Z_i Z_{(i+1)\bmod n}\rangle.

Their ideal values are −1-1 and +1+1, respectively. These parameters provide a benchmark for the workflow; they do not demonstrate a classically difficult computation.

def build_periodic_problem(num_qubits, num_reps=2):
    """Build a periodic benchmark circuit and its Z and ZZ observables."""
    edges = [(i, (i + 1) % num_qubits) for i in range(num_qubits)]
    circuit = n_local(
        num_qubits,
        rotation_blocks="rx",
        entanglement_blocks="cx",
        entanglement=edges[::2] + edges[1::2],
        reps=num_reps,
    )
    parameters = np.zeros(circuit.num_parameters)
    parameters[-num_qubits:] = np.pi
    circuit = circuit.assign_parameters(parameters)

    terms = [("Z", [i], 1) for i in range(num_qubits)]
    terms += [("ZZ", [i, (i + 1) % num_qubits], 1) for i in range(num_qubits)]
    observables = SparsePauliOp.from_sparse_list(terms, num_qubits)
    return circuit, observables


num_qubits = 6
circuit, observables = build_periodic_problem(num_qubits)
ideal_expvals = np.array(
    [
        Statevector(circuit).expectation_value(p).real
        for p in observables.paulis
    ]
)
print("Ideal averages:", ideal_expvals.reshape(2, num_qubits).mean(axis=1))
circuit.draw("mpl", style="iqp", fold=60)

Output:

Ideal averages: [-1.  1.]
Six-qubit periodic circuit with two boundary CX gates

Step 2: Optimize problem for quantum hardware execution

Create a device model with a linear coupling map. It supplies the compilation target; the execution in Step 3 will be noiseless.

simulation_backend = GenericBackendV2(
    num_qubits=num_qubits,
    basis_gates=["id", "rz", "sx", "x", "cx"],
    coupling_map=CouplingMap.from_line(num_qubits),
    seed=1234,
)
initial_layout = list(range(num_qubits))

Add a virtual connection between the physical qubits at the ends of the chain. Transpile once using this augmented coupling map and a CX basis. The virtual connection prevents the transpiler from routing the boundary gates, but the resulting circuit cannot run on the real device yet.

Find and cut the two boundary CX gates, apply the virtual layout to the observables, and generate all decomposition terms with num_samples=np.inf. This makes the decomposition coefficients exact; the measurements still use a finite number of shots.

Finally, translate the experiments into the device's native gates and correct their direction where needed. These passes preserve the physical qubit layout and avoid repeating layout selection and routing for every experiment. The helper returns the mapped observables needed for reconstruction.

def prepare_cut_experiments(circuit, observables, backend, initial_layout):
    """Compile with a virtual edge, cut it, and translate to native gates."""
    boundary = {initial_layout[0], initial_layout[-1]}
    virtual_map = CouplingMap(list(backend.coupling_map.get_edges()))
    virtual_map.make_symmetric()
    virtual_map.add_edge(initial_layout[-1], initial_layout[0])
    virtual_map.add_edge(initial_layout[0], initial_layout[-1])

    virtual_pm = generate_preset_pass_manager(
        optimization_level=1,
        coupling_map=virtual_map,
        basis_gates=["rz", "sx", "x", "cx"],
        initial_layout=initial_layout,
        seed_transpiler=1234,
    )
    virtual_circuit = virtual_pm.run(circuit)
    cut_indices = [
        i
        for i, instruction in enumerate(virtual_circuit.data)
        if instruction.operation.name == "cx"
        and {virtual_circuit.find_bit(q).index for q in instruction.qubits}
        == boundary
    ]
    mapped_observables = observables.apply_layout(virtual_circuit.layout)
    qpd_circuit, bases = cut_gates(virtual_circuit, cut_indices)
    experiments, coefficients = generate_cutting_experiments(
        circuits=qpd_circuit,
        observables=mapped_observables.paulis,
        num_samples=np.inf,
    )
    native_pm = PassManager(
        [
            BasisTranslator(
                SessionEquivalenceLibrary,
                target_basis=backend.operation_names,
                target=backend.target,
            ),
            GateDirection(backend.coupling_map, target=backend.target),
            BasisTranslator(
                SessionEquivalenceLibrary,
                target_basis=backend.operation_names,
                target=backend.target,
            ),
            Optimize1qGatesDecomposition(target=backend.target),
        ]
    )
    isa_experiments = native_pm.run(experiments)
    print(f"Boundary gates cut: {len(cut_indices)}")
    print(f"Number of experiments: {len(isa_experiments)}")
    print(
        f"Sampling overhead: {np.prod([basis.overhead for basis in bases])}"
    )
    return isa_experiments, coefficients, mapped_observables


isa_experiments, coefficients, mapped_observables = prepare_cut_experiments(
    circuit, observables, simulation_backend, initial_layout
)

Output:

Boundary gates cut: 2
Number of experiments: 36
Sampling overhead: 81.0

Each CX decomposition has six terms. With two cuts and one shared Z measurement basis, this produces 62=366^2=36 experiments. The sampling overhead is a different quantity: the squared sum of the absolute decomposition coefficients is 92=819^2=81. It describes the sampling cost of maintaining precision under quasiprobability reconstruction. See the sampling overhead explanation.

For comparison, transpile a measured copy of the uncut circuit using the actual coupling map. measure_all() records the original logical qubits in the meas register; transpilation preserves that classical-bit mapping, so its counts can be evaluated with the original observables.

uncut_pm = generate_preset_pass_manager(
    optimization_level=1,
    backend=simulation_backend,
    initial_layout=initial_layout,
    seed_transpiler=1234,
)
isa_uncut = uncut_pm.run(circuit.measure_all(inplace=False))
print("Uncut circuit depth:", isa_uncut.depth())
print("Cut circuit depths:", sorted({qc.depth() for qc in isa_experiments}))

Output:

Uncut circuit depth: 15
Cut circuit depths: [6, 7, 8, 9, 10, 11, 12, 15]

Step 3: Execute using Qiskit primitives

Run both versions with the Aer Sampler V2 primitive. Each cut experiment receives 10,000 shots, as does the single uncut circuit. The cut workflow therefore uses 36 times as many shots in this comparison. Matching shots per circuit does not match the total execution budget.

shots = 10_000
sampler = AerSampler(seed=1234)
cut_result = sampler.run(isa_experiments, shots=shots).result()
uncut_result = sampler.run([isa_uncut], shots=shots).result()

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

Reconstruct each Pauli expectation value using the same mapped observables and experiment order used during generation. Then average the first nn values for MZM_Z and the remaining nn values for CZZC_{ZZ}. For the uncut circuit, calculate the expectation values directly from the measured counts.

cut_expvals = reconstruct_expectation_values(
    cut_result, coefficients, mapped_observables.paulis
)
uncut_counts = uncut_result[0].data.meas.get_counts()
uncut_expvals = [
    sampled_expectation_value(uncut_counts, p) for p in observables.paulis
]


def compare_results(cut_expvals, uncut_expvals, num_qubits):
    """Display average Z and periodic ZZ values against the ideal values."""
    cut_averages = np.asarray(cut_expvals).reshape(2, num_qubits).mean(axis=1)
    uncut_averages = (
        np.asarray(uncut_expvals).reshape(2, num_qubits).mean(axis=1)
    )
    print("Cut averages (Z, ZZ):", cut_averages)
    print("Uncut averages (Z, ZZ):", uncut_averages)
    print("Absolute cut errors:", np.abs(cut_averages - [-1, 1]))
    print("Absolute uncut errors:", np.abs(uncut_averages - [-1, 1]))

    fig, ax = plt.subplots(figsize=(7, 4))
    positions = np.arange(2)
    for offset, values, label, color in [
        (-0.25, [-1, 1], "Ideal", "#4589ff"),
        (0, cut_averages, "Cut", "#a56eff"),
        (0.25, uncut_averages, "Uncut", "#009d9a"),
    ]:
        ax.bar(
            positions + offset, values, width=0.25, label=label, color=color
        )
    ax.set_xticks(positions, ["Average Z", "Average periodic ZZ"])
    ax.set_ylabel("Expectation value")
    ax.axhline(0, color="black", linewidth=0.8)
    ax.legend()
    fig.tight_layout()
    plt.show()


compare_results(cut_expvals, uncut_expvals, num_qubits)

Output:

Cut averages (Z, ZZ): [-1.00100833  0.99778333]
Uncut averages (Z, ZZ): [-1.  1.]
Absolute cut errors: [0.00100833 0.00221667]
Absolute uncut errors: [0. 0.]
Ideal, cut, and uncut expectation values from noiseless simulation

The noiseless results should be close to the ideal values. Finite-shot reconstruction can fluctuate, and a reconstructed estimate can lie outside [−1,1][-1,1]. Cutting does not improve an already noiseless result; this example verifies the decomposition, layout mapping, and reconstruction before execution on hardware.

Large-scale hardware example

Reuse the workflow on the least-busy operational QPU with at least 120 qubits. Select a longest shortest path in its undirected coupling graph: each consecutive pair is connected, while the endpoints are far apart. This adapts the chain to the selected device instead of relying on a fixed Eagle layout. The code checks that the chain contains at least 30 qubits and that its endpoints are not neighbors.

The following block combines the four steps. Gate twirling, measurement twirling, and dynamical decoupling are disabled for both versions. Both jobs run in the same batch and receive the tutorial tag TUT_CCPBC. The cut experiments still use more total shots, so compare both accuracy and QPU usage.

service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, simulator=False, min_num_qubits=120
)

# ------------------------- Step 1 -------------------------
graph = backend.coupling_map.graph.to_undirected(multigraph=False)
paths = rx.all_pairs_dijkstra_shortest_paths(graph, lambda _: 1)
initial_layout = list(
    max(
        (
            path
            for source_paths in paths.values()
            for path in source_paths.values()
        ),
        key=len,
    )
)
num_qubits = len(initial_layout)
if num_qubits < 30:
    raise ValueError(
        "Select a backend with a shortest path of at least 30 qubits."
    )
if graph.has_edge(initial_layout[0], initial_layout[-1]):
    raise ValueError("The chain endpoints must not be directly connected.")
print(f"Backend: {backend.name}; chain length: {num_qubits}")
print("Initial layout:", initial_layout)
circuit, observables = build_periodic_problem(num_qubits)

# ------------------------- Step 2 -------------------------
isa_experiments, coefficients, mapped_observables = prepare_cut_experiments(
    circuit, observables, backend, initial_layout
)
uncut_pm = generate_preset_pass_manager(
    optimization_level=1,
    backend=backend,
    initial_layout=initial_layout,
    seed_transpiler=1234,
)
isa_uncut = uncut_pm.run(circuit.measure_all(inplace=False))
print("Uncut circuit depth:", isa_uncut.depth())
print("Cut circuit depths:", sorted({qc.depth() for qc in isa_experiments}))

# ------------------------- Step 3 -------------------------
options = SamplerOptions()
options.dynamical_decoupling.enable = False
options.twirling.enable_gates = False
options.twirling.enable_measure = False
options.environment.job_tags = ["TUT_CCPBC"]
shots = 10_000
with Batch(backend=backend) as batch:
    sampler = SamplerV2(mode=batch, options=options)
    cut_job = sampler.run(isa_experiments, shots=shots)
    uncut_job = sampler.run([isa_uncut], shots=shots)
print("Cut job ID:", cut_job.job_id())
print("Uncut job ID:", uncut_job.job_id())
cut_result = cut_job.result()
uncut_result = uncut_job.result()

# ------------------------- Step 4 -------------------------
cut_expvals = reconstruct_expectation_values(
    cut_result, coefficients, mapped_observables.paulis
)
uncut_counts = uncut_result[0].data.meas.get_counts()
uncut_expvals = [
    sampled_expectation_value(uncut_counts, p) for p in observables.paulis
]
compare_results(cut_expvals, uncut_expvals, num_qubits)
qpu_seconds = sum(
    job.metrics()["usage"]["quantum_seconds"] for job in [cut_job, uncut_job]
)
print(
    f"Total QPU usage: {qpu_seconds:.1f} seconds ({qpu_seconds / 60:.2f} minutes)"
)

Output:

Backend: ibm_boston; chain length: 33
Initial layout: [0, 1, 2, 3, 16, 23, 24, 25, 37, 45, 44, 43, 56, 63, 64, 65, 77, 85, 86, 87, 97, 107, 106, 105, 117, 125, 124, 123, 136, 143, 142, 141, 140]
Boundary gates cut: 2
Number of experiments: 36
Sampling overhead: 81.0
Uncut circuit depth: 375
Cut circuit depths: [20]
Cut job ID: dagp7fj9k43c73admp40
Uncut job ID: dagp7fj9k43c73admp5g
Cut averages (Z, ZZ): [-0.92305152  0.89586667]
Uncut averages (Z, ZZ): [-0.61106061  0.51606061]
Absolute cut errors: [0.07694848 0.10413333]
Absolute uncut errors: [0.38893939 0.48393939]
Total QPU usage: 102.0 seconds (1.70 minutes)
Ideal, cut, and uncut expectation values measured on ibm_boston

In the recorded run on ibm_boston, the cut estimates were approximately MZ=−0.923M_Z=-0.923 and CZZ=0.896C_{ZZ}=0.896, compared with −0.611-0.611 and 0.5160.516 for the uncut circuit. Both cut estimates were closer to the ideal values of −1-1 and +1+1.

Compare each result with its ideal value and inspect the circuit depths and total QPU usage. Cutting removes the need to route the periodic gates, but introduces intermediate operations and additional statistical uncertainty. Hardware noise and calibration affect the comparison; an improvement is not guaranteed. Repeated runs are needed to assess whether an observed difference is statistically significant.

Next steps

Recommendations
Was this page helpful?
Report a bug, typo, or request content on GitHub.